pkg/darwin/darwin_conserve_surface.F

8fbfd1f382 Oliv*0001 #include "DARWIN_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: DARWIN_CONSERVE_SURFACE
                0005 C !INTERFACE: ==========================================================
                0006       SUBROUTINE DARWIN_CONSERVE_SURFACE(myTime,myIter,myThid)
                0007
                0008 C !DESCRIPTION:
                0009 C     Fix conservation violation due to linear free surface
                0010
                0011 C !USES: ===============================================================
                0012       IMPLICIT NONE
                0013 #include "SIZE.h"
                0014 #include "EEPARAMS.h"
                0015 #include "PARAMS.h"
                0016 #include "GRID.h"
                0017 #include "PTRACERS_SIZE.h"
                0018 #include "PTRACERS_PARAMS.h"
                0019 #include "PTRACERS_FIELDS.h"
                0020 #include "DARWIN_SIZE.h"
                0021 #include "DARWIN_INDICES.h"
                0022 #include "DARWIN_PARAMS.h"
                0023 #include "DARWIN_TRAITS.h"
                0024 #include "DARWIN_FIELDS.h"
                0025
                0026 C !INPUT PARAMETERS: ===================================================
                0027       _RL myTime
                0028       INTEGER myIter
                0029       INTEGER myThid
                0030 CEOP
                0031
                0032 #ifdef ALLOW_DARWIN
                0033
                0034 C !LOCAL VARIABLES: ====================================================
                0035       INTEGER i,j,k,bi,bj
                0036       INTEGER np, iTr
                0037       _RL totflxC
                0038       _RL totflxP
                0039       _RL totflxN
                0040       _RL totflxFe
                0041       _RL totflxSi
                0042 #ifdef DARWIN_ALLOW_CARBON
                0043       _RL totflxA
                0044       _RL totflxO
                0045 #endif
                0046       _RL vol, voltot
8c49863325 Oliv*0047       _RL tilevol(nSx,nSy)
8fbfd1f382 Oliv*0048
                0049       IF ( DARWIN_linFSConserve ) THEN
                0050
                0051       totflxC = 0. _d 0
                0052       totflxP = 0. _d 0
                0053       totflxN = 0. _d 0
                0054       totflxFe = 0. _d 0
                0055       totflxSi = 0. _d 0
                0056 #ifdef DARWIN_ALLOW_CARBON
                0057       totflxA = 0. _d 0
                0058       totflxO = 0. _d 0
                0059 #endif
                0060
                0061 c carbon budget
                0062       totflxC=totflxC+totSurfCorPtr(iDIC)
                0063       totflxC=totflxC+totSurfCorPtr(iDOC)
                0064       totflxC=totflxC+totSurfCorPtr(iPOC)
                0065       totflxC=totflxC+totSurfCorPtr(iPIC)
                0066       DO np=1,nplank
                0067         iTr=ic+np-1
4cff8f255d Oliv*0068         totflxC=totflxC+totSurfCorPtr(iTr)*(1+R_PICPOC(np))
8fbfd1f382 Oliv*0069       ENDDO
c7b6c66d45 Oliv*0070 #ifdef DARWIN_ALLOW_CSTORE
                0071       DO np=1,nPhoto
                0072         iTr=ich+np-1
                0073         totflxC=totflxC+totSurfCorPtr(iTr)
                0074       ENDDO
                0075 #endif
8fbfd1f382 Oliv*0076 #ifdef DARWIN_ALLOW_CDOM
1c72adde5c Oliv*0077 # ifdef DARWIN_CDOM_UNITS_CARBON
                0078       totflxC=totflxC+totSurfCorPtr(iCDOM)
                0079 # else
8fbfd1f382 Oliv*0080       totflxC=totflxC+totSurfCorPtr(iCDOM)*R_CP_CDOM
1c72adde5c Oliv*0081 # endif
8fbfd1f382 Oliv*0082 #endif
                0083
                0084 c Phosphorus budget
                0085       totflxP=totflxP+totSurfCorPtr(iPO4)
                0086       totflxP=totflxP+totSurfCorPtr(iDOP)
                0087       totflxP=totflxP+totSurfCorPtr(iPOP)
                0088       DO np=1,nplank
                0089 #ifdef DARWIN_ALLOW_PQUOTA
                0090         iTr=ip+np-1
                0091         totflxP=totflxP+totSurfCorPtr(iTr)
                0092 #else
                0093         iTr=ic+np-1
                0094         totflxP=totflxP+totSurfCorPtr(iTr)*R_PC(np)
                0095 #endif
                0096       ENDDO
                0097 #ifdef DARWIN_ALLOW_CDOM
1c72adde5c Oliv*0098 # ifdef DARWIN_CDOM_UNITS_CARBON
                0099       totflxP=totflxP+totSurfCorPtr(iCDOM)*R_PC_CDOM
                0100 # else
8fbfd1f382 Oliv*0101       totflxP=totflxP+totSurfCorPtr(iCDOM)
1c72adde5c Oliv*0102 # endif
8fbfd1f382 Oliv*0103 #endif
                0104
                0105 c Nitrogen budget
                0106       totflxN=totflxN+totSurfCorPtr(iNO3)
                0107       totflxN=totflxN+totSurfCorPtr(iNH4)
                0108       totflxN=totflxN+totSurfCorPtr(iNO2)
                0109       totflxN=totflxN+totSurfCorPtr(iDON)
                0110       totflxN=totflxN+totSurfCorPtr(iPON)
                0111       DO np=1,nplank
                0112 #ifdef DARWIN_ALLOW_NQUOTA
                0113         iTr=in+np-1
                0114         totflxN=totflxN+totSurfCorPtr(iTr)
                0115 #else
                0116         iTr=ic+np-1
                0117         totflxN=totflxN+totSurfCorPtr(iTr)*R_NC(np)
                0118 #endif
                0119       ENDDO
                0120 #ifdef DARWIN_ALLOW_CDOM
1c72adde5c Oliv*0121 # ifdef DARWIN_CDOM_UNITS_CARBON
                0122       totflxN=totflxN+totSurfCorPtr(iCDOM)*R_NC_CDOM
                0123 # else
8fbfd1f382 Oliv*0124       totflxN=totflxN+totSurfCorPtr(iCDOM)*R_NP_CDOM
1c72adde5c Oliv*0125 # endif
8fbfd1f382 Oliv*0126 #endif
                0127
                0128 c Iron budget
                0129       totflxFe=totflxFe+totSurfCorPtr(iFeT)
                0130       totflxFe=totflxFe+totSurfCorPtr(iDOFe)
                0131       totflxFe=totflxFe+totSurfCorPtr(iPOFe)
                0132       DO np=1,nplank
                0133 #ifdef DARWIN_ALLOW_FEQUOTA
                0134         iTr=ife+np-1
                0135         totflxFe=totflxFe+totSurfCorPtr(iTr)
                0136 #else
                0137         iTr=ic+np-1
                0138         totflxFe=totflxFe+totSurfCorPtr(iTr)*R_FeC(np)
                0139 #endif
                0140       ENDDO
                0141 #ifdef DARWIN_ALLOW_CDOM
1c72adde5c Oliv*0142 # ifdef DARWIN_CDOM_UNITS_CARBON
                0143       totflxFe=totflxFe+totSurfCorPtr(iCDOM)*R_FeC_CDOM
                0144 # else
8fbfd1f382 Oliv*0145       totflxFe=totflxFe+totSurfCorPtr(iCDOM)*R_FeP_CDOM
1c72adde5c Oliv*0146 # endif
8fbfd1f382 Oliv*0147 #endif
                0148
                0149 c Silica budget
                0150       totflxSi=totflxSi+totSurfCorPtr(iSiO2)
                0151       totflxSi=totflxSi+totSurfCorPtr(iPOSi)
                0152       DO np=1,nplank
                0153 #ifdef DARWIN_ALLOW_SIQUOTA
                0154         iTr=isi+np-1
                0155         totflxSi=totflxSi+totSurfCorPtr(iTr)
                0156 #else
                0157         iTr=ic+np-1
                0158         totflxSi=totflxSi+totSurfCorPtr(iTr)*R_SiC(np)
                0159 #endif
                0160       ENDDO
                0161 c------------------------------------------------------------
                0162 #ifdef DARWIN_ALLOW_CARBON
                0163 c alkalinity budget !!!!! NOT FULL BUDGET!!!!!
                0164       totflxA=totflxA+totSurfCorPtr(iALK)
                0165 c oxygem budget !!!!!!!! NOT FULL BUDGET !!!!!!!!!!!
                0166       totflxO=totflxO+totSurfCorPtr(iO2)
                0167 #endif
                0168
                0169       DO bj=myByLo(myThid),myByHi(myThid)
                0170       DO bi=myBxLo(myThid),myBxHi(myThid)
8c49863325 Oliv*0171         tilevol(bi,bj) = 0. _d 0
8fbfd1f382 Oliv*0172         DO k=1,Nr
                0173         DO i=1,sNx
                0174         DO j=1,sNy
                0175           vol=rA(i,j,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj)
8c49863325 Oliv*0176           tilevol(bi,bj) = tilevol(bi,bj) + vol
8fbfd1f382 Oliv*0177         ENDDO
                0178         ENDDO
                0179         ENDDO
                0180       ENDDO
                0181       ENDDO
8c49863325 Oliv*0182       CALL GLOBAL_SUM_TILE_RL(tilevol,voltot,myThid)
8fbfd1f382 Oliv*0183
                0184 C--   Everyone else must wait
                0185       _BARRIER
                0186
                0187 C     mean tracer change
                0188       totflxC  = ptracers_dTLev(1)*totflxC /voltot
                0189       totflxN  = ptracers_dTLev(1)*totflxN /voltot
                0190       totflxP  = ptracers_dTLev(1)*totflxP /voltot
                0191       totflxFe = ptracers_dTLev(1)*totflxFe/voltot
                0192       totflxSi = ptracers_dTLev(1)*totflxSi/voltot
                0193 #ifdef DARWIN_ALLOW_CARBON
                0194       totflxA  = ptracers_dTLev(1)*totflxA /voltot
                0195       totflxO  = ptracers_dTLev(1)*totflxO /voltot
                0196 #endif
                0197
                0198 C     add neglected flux to tracer globally
                0199       DO bj=myByLo(myThid),myByHi(myThid)
                0200        DO bi=myBxLo(myThid),myBxHi(myThid)
                0201         DO k=1,Nr
                0202          DO i=1,sNx
                0203           DO j=1,sNy
                0204            Ptracer(i,j,k,bi,bj,iDIC)=Ptracer(i,j,k,bi,bj,iDIC)+
                0205      &                               totflxC*maskC(i,j,k,bi,bj)
                0206            Ptracer(i,j,k,bi,bj,iPO4)=Ptracer(i,j,k,bi,bj,iPO4)+
                0207      &                               totflxP *maskC(i,j,k,bi,bj)
                0208            Ptracer(i,j,k,bi,bj,iNO3)=Ptracer(i,j,k,bi,bj,iNO3)+
                0209      &                               totflxN *maskC(i,j,k,bi,bj)
                0210            Ptracer(i,j,k,bi,bj,iFeT)=Ptracer(i,j,k,bi,bj,iFeT)+
                0211      &                               totflxFe*maskC(i,j,k,bi,bj)
                0212            Ptracer(i,j,k,bi,bj,iSiO2)=Ptracer(i,j,k,bi,bj,iSiO2)+
                0213      &                               totflxSi*maskC(i,j,k,bi,bj)
                0214 #ifdef DARWIN_ALLOW_CARBON
                0215            Ptracer(i,j,k,bi,bj,iALK)=Ptracer(i,j,k,bi,bj,iALK)+
                0216      &                               totflxA*maskC(i,j,k,bi,bj)
                0217            Ptracer(i,j,k,bi,bj,iO2 )=Ptracer(i,j,k,bi,bj,iO2 )+
                0218      &                               totflxO*maskC(i,j,k,bi,bj)
                0219 #endif
                0220           ENDDO
                0221          ENDDO
                0222         ENDDO
                0223        ENDDO
                0224       ENDDO
                0225
                0226 C     DARWIN_linFSConserve
                0227       ENDIF
                0228
                0229 #endif /* ALLOW_DARWIN */
                0230
                0231       RETURN
                0232       END
                0233