pkg/darwin/darwin_surfforcing.F

8fbfd1f382 Oliv*0001 #include "DARWIN_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: DARWIN_SURFFORCING
                0005
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE DARWIN_SURFFORCING(
                0008      O           gDIC, gALK, gO2,
5e411acc9e Oliv*0009      I           dTsub,
8fbfd1f382 Oliv*0010      I           bi,bj,imin,imax,jmin,jmax,
                0011      I           myTime,myIter,myThid)
                0012
                0013 C !DESCRIPTION:
                0014 C Update tendency terms for alkalinity, oxygen and DIC from air-sea
                0015 C exchanges
                0016
                0017 C !USES: ===============================================================
                0018       IMPLICIT NONE
                0019 #include "SIZE.h"
                0020 #include "DYNVARS.h"
                0021 #include "EEPARAMS.h"
                0022 #include "PARAMS.h"
                0023 #include "GRID.h"
                0024 #include "FFIELDS.h"
                0025 #ifdef ALLOW_PTRACERS
                0026 #include "PTRACERS_SIZE.h"
                0027 #include "PTRACERS_FIELDS.h"
                0028 #endif
                0029 #ifdef ALLOW_DARWIN
                0030 #include "DARWIN_SIZE.h"
                0031 #include "DARWIN_INDICES.h"
                0032 #include "DARWIN_EXF_FIELDS.h"
                0033 #include "DARWIN_PARAMS.h"
                0034 #include "DARWIN_TRAITS.h"
                0035 #include "DARWIN_FIELDS.h"
                0036 #include "DARWIN_FLUX.h"
                0037 #endif
                0038
                0039 C !INPUT PARAMETERS: ===================================================
                0040 C  myThid               :: thread number
                0041 C  myIter               :: current timestep
                0042 C  myTime               :: current time
                0043 C  bi,bj                :: tile indices
5e411acc9e Oliv*0044       _RL dTsub(Nr)
8fbfd1f382 Oliv*0045       INTEGER iMin,iMax,jMin,jMax,bi,bj
                0046       INTEGER myIter, myThid
                0047       _RL myTime
                0048
                0049 C !OUTPUT PARAMETERS: ==================================================
                0050 C  gDIC :: DIC tendency due to air-sea exchange
                0051 C  gALK :: ALK tendency due to air-sea exchange
                0052 C  gO2  :: O2 tendency due to air-sea exchange
                0053       _RL  gDIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0054       _RL  gALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0055       _RL  gO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0056
                0057 #ifdef ALLOW_DARWIN
                0058 #ifdef DARWIN_ALLOW_CARBON
                0059
                0060 C !LOCAL VARIABLES: ====================================================
                0061       INTEGER i,j,k,ks
f0a72e2151 Oliv*0062       LOGICAL debugPrt
8fbfd1f382 Oliv*0063 C Number of iterations for pCO2 solvers...
                0064 C Solubility relation coefficients
                0065       _RL SchmidtNoDIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
9d53e28eac Oliv*0066       _RL apCO2sat(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8fbfd1f382 Oliv*0067       _RL Kwexch(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0068       _RL Kwexch0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0069       _RL pisvel0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0070 C local variables for carbon chem
                0071       _RL surfdic
                0072       _RL surfalk
                0073       _RL surfphos
                0074       _RL surfsi
                0075       _RL surftemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0076       _RL surfsalt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0077       _RL SchmidtNoO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0078       _RL O2sat(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
895c6145db Oliv*0079       _RL FluxO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0080       _RL pCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0081       _RL CO3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a97de47f46 Oliv*0082       _RL fugCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
895c6145db Oliv*0083       _RL FluxCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0084       _RL VFluxCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0085       _RL VFluxAlk(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8fbfd1f382 Oliv*0086       _RL aTT
                0087       _RL aTK
                0088       _RL aTS
                0089       _RL aTS2
                0090       _RL aTS3
                0091       _RL aTS4
                0092       _RL aTS5
                0093       _RL o2s
                0094       _RL ttemp
                0095       _RL stemp
                0096       _RL oCnew
1a75f0d7fd Oliv*0097       _RL calcium
8fbfd1f382 Oliv*0098 CEOP
                0099
                0100       ks = 1
                0101
                0102       DO j=jmin,jmax
                0103        DO i=imin,imax
                0104          gDIC(i,j) = 0.0 _d 0
                0105          gALK(i,j) = 0.0 _d 0
                0106          gO2(i,j) = 0.0 _d 0
                0107        ENDDO
                0108       ENDDO
                0109
                0110 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
                0111 C Compute AtmosP and Kwexch0 which are used for flux of CO2 and O2
                0112       DO j=jmin,jmax
                0113        DO i=imin,imax
                0114         IF (maskC(i,j,ks,bi,bj).NE.0. _d 0) THEN
                0115
bee566be5e Oliv*0116 #ifdef DARWIN_USE_PLOAD
8fbfd1f382 Oliv*0117 C Convert anomalous pressure pLoad (in Pa) from atmospheric model
                0118 C to total pressure (in Atm)
                0119 C Note: it is assumed the reference atmospheric pressure is 1Atm=1013mb
                0120 C       rather than the actual ref. pressure from Atm. model so that on
                0121 C       average AtmosP is about 1 Atm.
bee566be5e Oliv*0122          AtmosP(i,j,bi,bj)= (surf_pRef + pLoad(i,j,bi,bj))/Pa2Atm
8fbfd1f382 Oliv*0123 #endif
                0124
                0125 C Pre-compute part of exchange coefficient: pisvel*(1-fice)
                0126 C Schmidt number is accounted for later
                0127          pisvel0(i,j) = 0.337 _d 0 * windSpeed(i,j,bi,bj)**2/3.6 _d 5
                0128          Kwexch0(i,j) = pisvel0(i,j) *
                0129      &                   (1. _d 0 - iceFrac(i,j,bi,bj))
                0130
                0131         ENDIF
                0132        ENDDO
                0133       ENDDO
                0134
                0135 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
                0136 c flux of alkalinity
                0137
                0138 #ifdef ALLOW_OLD_VIRTUALFLUX
                0139       DO j=jmin,jmax
                0140        DO i=imin,imax
                0141         IF (maskC(i,j,ks,bi,bj).NE.0. _d 0) THEN
                0142 c calculate virtual flux
                0143 c EminusPforV = dS/dt*(1/Sglob)
                0144 C NOTE: Be very careful with signs here!
                0145 C Positive EminusPforV => loss of water to atmos and increase
                0146 C in salinity. Thus, also increase in other surface tracers
                0147 C (i.e. positive virtual flux into surface layer)
895c6145db Oliv*0148          VFluxAlk(i,j) = gsm_ALK*surfaceForcingS(i,j,bi,bj)/gsm_S
8fbfd1f382 Oliv*0149          gALK(i,j) =
895c6145db Oliv*0150      &     recip_drF(ks)*recip_hFacC(i,j,ks,bi,bj)*VFluxAlk(i,j)
d95627305e Oliv*0151 #ifdef DARWIN_ALLOW_CONS
5e411acc9e Oliv*0152          alkVirFlx(i,j,bi,bj) = alkVirFlx(i,j,bi,bj)
895c6145db Oliv*0153      &                        + VFluxAlk(i,j)*dTsub(ks)
d95627305e Oliv*0154 #endif
a59ebcf0a3 Oliv*0155         ELSE
895c6145db Oliv*0156          VFluxAlk(i,j) = 0 _d 0
8fbfd1f382 Oliv*0157         ENDIF
                0158        ENDDO
                0159       ENDDO
                0160 #endif /* ALLOW_OLD_VIRTUALFLUX */
                0161
                0162 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
                0163 c flux of O2
                0164
                0165 C calculate SCHMIDT NO. for O2
                0166       DO j=jmin,jmax
                0167        DO i=imin,imax
                0168         IF (maskC(i,j,ks,bi,bj).NE.0.) THEN
                0169          ttemp = theta(i,j,ks,bi,bj)
                0170          stemp = salt(i,j,ks,bi,bj)
                0171
                0172          SchmidtNoO2(i,j) =
                0173      &       sox1
                0174      &     + sox2 * ttemp
                0175      &     + sox3 * ttemp*ttemp
                0176      &     + sox4 * ttemp*ttemp*ttemp
                0177
                0178 C Determine surface flux of O2
                0179 C exchange coeff accounting for ice cover and Schmidt no.
                0180 C Kwexch0= pisvel*(1-fice): previously computed in dic_surfforcing.F
                0181
                0182          Kwexch(i,j) = Kwexch0(i,j)
                0183      &                / sqrt(SchmidtNoO2(i,j)/660.0 _d 0)
                0184
                0185 C determine saturation O2
                0186 C using Garcia and Gordon (1992), L&O (mistake in original ?)
                0187          aTT  = 298.15 _d 0 -ttemp
                0188          aTK  = 273.15 _d 0 +ttemp
                0189          aTS  = log(aTT/aTK)
                0190          aTS2 = aTS*aTS
                0191          aTS3 = aTS2*aTS
                0192          aTS4 = aTS3*aTS
                0193          aTS5 = aTS4*aTS
                0194
                0195          oCnew = oA0 + oA1*aTS + oA2*aTS2 + oA3*aTS3 +
                0196      &           oA4*aTS4 + oA5*aTS5 +
                0197      &           (oB0 + oB1*aTS + oB2*aTS2 + oB3*aTS3)*stemp +
                0198      &           oC0*stemp*stemp
                0199
                0200          o2s = EXP(oCnew)
                0201
                0202 c molar volume of O2: O2mol2L = 22.391 L mol-1 = 0.022391 L mmol-1
                0203 c o2s in ml/l = l/m3
                0204 c O2sat = o2s / O2mmol2L  (in mmol/m3)
                0205 c Convert from ml/l to mmol/m^3
                0206          O2sat(i,j) = o2s/22.3916 _d -3
                0207
                0208 C Determine flux, inc. correction for local atmos surface pressure
895c6145db Oliv*0209          FluxO2(i,j) = Kwexch(i,j)*
8fbfd1f382 Oliv*0210      &                     (AtmosP(i,j,bi,bj)*O2sat(i,j)
                0211      &                      - pTracer(i,j,ks,bi,bj,iO2))
d61a0bd101 Oliv*0212 #ifdef DARWIN_ALLOW_CONS
5e411acc9e Oliv*0213          oxySfcFlx(i,j,bi,bj) = oxySfcFlx(i,j,bi,bj)
895c6145db Oliv*0214      &                        + FluxO2(i,j)*dTsub(ks)
d61a0bd101 Oliv*0215 #endif
8fbfd1f382 Oliv*0216
                0217          gO2(i,j) =
895c6145db Oliv*0218      &      recip_drF(ks)*recip_hFacC(i,j,ks,bi,bj)*FluxO2(i,j)
3a1d328575 Oliv*0219         ELSE
                0220          FluxO2(i,j) = 0. _d 0
8fbfd1f382 Oliv*0221         ENDIF
                0222        ENDDO
                0223       ENDDO
                0224
                0225 C ======================================================================
                0226       DO k=1,Nr
                0227 C ======================================================================
                0228
                0229 C determine inorganic carbon chem coefficients
14f3186978 Oliv*0230        DO j=jmin,jmax
                0231         DO i=imin,imax
8fbfd1f382 Oliv*0232 c put bounds on tracers so pH solver doesn't blow up
14f3186978 Oliv*0233          surfsalt(i,j) = MAX(surfSaltMin,
                0234      &                   MIN(surfSaltMax, salt(i,j,k,bi,bj)))
                0235          surftemp(i,j) = MAX(surfTempMin,
                0236      &                   MIN(surfTempMax, theta(i,j,k,bi,bj)))
                0237         ENDDO
8fbfd1f382 Oliv*0238        ENDDO
                0239
f0a72e2151 Oliv*0240 #ifdef DARWIN_SOLVESAPHE
                0241 #ifdef ALLOW_DEBUG
14f3186978 Oliv*0242        IF (debugMode) CALL DEBUG_CALL(
                0243      &    'DARWIN_COEFFS_DEEP',myThid)
f0a72e2151 Oliv*0244 #endif
                0245 C Calculate carbon coefficients
14f3186978 Oliv*0246        CALL DARWIN_COEFFS_SURF(
                0247      I                          surftemp, surfsalt,
45ab1c3374 Oliv*0248      I                          bi,bj,iMin,iMax,jMin,jMax,
                0249      I                          k,myThid)
f0a72e2151 Oliv*0250
                0251 C Now correct the coefficients for pressure dependence
14f3186978 Oliv*0252        CALL DARWIN_COEFFS_DEEP(
                0253      I                          surftemp, surfsalt,
                0254      I                          bi,bj,iMin,iMax,jMin,jMax,
                0255      I                          k,myThid)
f0a72e2151 Oliv*0256 #else
                0257 #ifdef ALLOW_DEBUG
14f3186978 Oliv*0258        IF (debugMode) CALL DEBUG_CALL('DARWIN_CARBON_COEFFS',myThid)
f0a72e2151 Oliv*0259 #endif
14f3186978 Oliv*0260        CALL DARWIN_CARBON_COEFFS(
                0261      I                            surftemp,surfsalt,
                0262      I                            bi,bj,iMin,iMax,jMin,jMax,k,myThid)
f0a72e2151 Oliv*0263 #endif /* DARWIN_SOLVESAPHE */
8fbfd1f382 Oliv*0264
f0a72e2151 Oliv*0265 C====================================================================
                0266        debugPrt = debugMode
8fbfd1f382 Oliv*0267 c pCO2 solver...
14f3186978 Oliv*0268        DO j=jmin,jmax
                0269         DO i=imin,imax
                0270          IF ( maskC(i,j,k,bi,bj).NE.0. _d 0 ) THEN
f0a72e2151 Oliv*0271
8fbfd1f382 Oliv*0272 c put bounds on tracers so pH solver doesn't blow up
14f3186978 Oliv*0273           surfdic  = ptr2mol * MAX(surfDICMin,
                0274      &               MIN(surfDICMax, Ptracer(i,j,k,bi,bj,iDIC)))
                0275      &               * maskC(i,j,k,bi,bj)
                0276           surfalk  = ptr2mol * MAX(surfALKMin,
                0277      &               MIN(surfALKMax, Ptracer(i,j,k,bi,bj,iALK)))
                0278      &               * maskC(i,j,k,bi,bj)
                0279           surfphos = ptr2mol * MAX(surfPO4Min,
                0280      &               MIN(surfPO4Max, Ptracer(i,j,k,bi,bj,iPO4)))
                0281      &               * maskC(i,j,k,bi,bj)
20ff2335d4 Oliv*0282           surfsi   = ptr2mol * MAX(surfSiMin,
14f3186978 Oliv*0283      &               MIN(surfSiMax, Ptracer(i,j,k,bi,bj,iSiO2)))
                0284      &               * maskC(i,j,k,bi,bj)
f0a72e2151 Oliv*0285
                0286 #ifdef DARWIN_SOLVESAPHE
14f3186978 Oliv*0287           IF ( selectPHsolver.GT.0 ) THEN
f0a72e2151 Oliv*0288 C Use Munhoven (2013) Solvesaphe routine to calculate pH and pCO2
                0289 #ifdef ALLOW_DEBUG
                0290             IF (debugPrt) CALL DEBUG_CALL('CALC_PCO2_SOLVESAPHE',myThid)
                0291 #endif
                0292             CALL CALC_PCO2_SOLVESAPHE(
                0293      I          surftemp(i,j), surfsalt(i,j),
                0294      I          surfdic, surfphos, surfsi, surfalk,
895c6145db Oliv*0295      U          pH(i,j,k,bi,bj),
                0296      O          pCO2(i,j), CO3(i,j),
f0a72e2151 Oliv*0297      I          i,j,k,bi,bj, debugPrt,myIter,myThid )
                0298             debugPrt = .FALSE.
14f3186978 Oliv*0299           ELSE
f0a72e2151 Oliv*0300 C Use the original Follows et al. (2006) solver
                0301 #endif /* DARWIN_SOLVESAPHE */
                0302 #ifdef ALLOW_DEBUG
                0303             IF (debugPrt) CALL DEBUG_CALL('CALC_PCO2_APPROX',myThid)
                0304             debugPrt = .FALSE.
                0305 #endif
                0306             CALL DARWIN_CALC_PCO2_APPROX(
8fbfd1f382 Oliv*0307      I        surftemp(i,j), surfsalt(i,j),
                0308      I        surfdic, surfphos, surfsi, surfalk,
                0309      I        ak1(i,j,bi,bj), ak2(i,j,bi,bj),
                0310      I        ak1p(i,j,bi,bj), ak2p(i,j,bi,bj), ak3p(i,j,bi,bj),
                0311      I        aks(i,j,bi,bj), akb(i,j,bi,bj), akw(i,j,bi,bj),
                0312      I        aksi(i,j,bi,bj), akf(i,j,bi,bj),
                0313      I        ak0(i,j,bi,bj),  fugf(i,j,bi,bj),
                0314      I        ff(i,j,bi,bj),
                0315      I        bt(i,j,bi,bj), st(i,j,bi,bj), ft(i,j,bi,bj),
895c6145db Oliv*0316      U        pH(i,j,k,bi,bj),
                0317      O        pCO2(i,j), CO3(i,j),
8fbfd1f382 Oliv*0318      I        i,j,k,bi,bj,myIter,myThid )
f0a72e2151 Oliv*0319 #ifdef DARWIN_SOLVESAPHE
14f3186978 Oliv*0320           ENDIF
f0a72e2151 Oliv*0321 #endif /* DARWIN_SOLVESAPHE */
                0322
1a75f0d7fd Oliv*0323 C Compute calcite saturation ratio omegaC
                0324 #ifdef DARWIN_SOLVESAPHE
                0325           calcium = cat(i,j,bi,bj)
                0326 #else
                0327 C         calcium concentration in mol/kg-SW
                0328           calcium = 1.028 _d -2*salt(i,j,k,bi,bj)/35. _d 0
                0329 #endif
                0330 C         [calcium] = mol/kg
                0331 C         [CO3] = mol/kg
                0332 C         [Ksp_TP_Calc] = mol^2/kg^2
                0333           omegaC(i,j,k,bi,bj) = MAX(0 _d 0, calcium * CO3(i,j)) /
                0334      &                          Ksp_TP_Calc(i,j,bi,bj)
                0335
14f3186978 Oliv*0336          ELSE
895c6145db Oliv*0337           pCO2(i,j)=0. _d 0
                0338           CO3(i,j) = 0. _d 0
1a75f0d7fd Oliv*0339           omegaC(i,j,k,bi,bj) = 0 _d 0
14f3186978 Oliv*0340          ENDIF
1a75f0d7fd Oliv*0341
                0342 C      enddo i,j
14f3186978 Oliv*0343         ENDDO
8fbfd1f382 Oliv*0344        ENDDO
                0345
14f3186978 Oliv*0346        IF (k .EQ. ks) THEN
8fbfd1f382 Oliv*0347
14f3186978 Oliv*0348         DO j=jmin,jmax
                0349          DO i=imin,imax
                0350           IF ( maskC(i,j,ks,bi,bj).NE.0. _d 0 ) THEN
                0351            ttemp = theta(i,j,ks,bi,bj)
8fbfd1f382 Oliv*0352 C calculate SCHMIDT NO. for CO2
14f3186978 Oliv*0353            SchmidtNoDIC(i,j) =
                0354      &        sca1 +
                0355      &        sca2*ttemp +
                0356      &        sca3*ttemp*ttemp +
                0357      &        sca4*ttemp*ttemp*ttemp
8fbfd1f382 Oliv*0358 c make sure Schmidt number is not negative (will happen if temp>39C)
14f3186978 Oliv*0359            SchmidtNoDIC(i,j) = MAX(1.0 _d -2, SchmidtNoDIC(i,j))
8fbfd1f382 Oliv*0360
                0361 C Determine surface flux (FDIC)
9d53e28eac Oliv*0362 C first correct apCO2 for surface atmos pressure
14f3186978 Oliv*0363            apCO2sat(i,j) = AtmosP(i,j,bi,bj)*atmospCO2(i,j,bi,bj)
8fbfd1f382 Oliv*0364
                0365 C then account for Schmidt number
14f3186978 Oliv*0366            Kwexch(i,j) = Kwexch0(i,j)
                0367      &                   / sqrt(SchmidtNoDIC(i,j)/660.0 _d 0)
8fbfd1f382 Oliv*0368
a97de47f46 Oliv*0369 C compute fugacity of CO2
895c6145db Oliv*0370            fugCO2(i,j) = pCO2(i,j)*fugf(i,j,bi,bj)
9d53e28eac Oliv*0371
8fbfd1f382 Oliv*0372 c Flux = kw*rho*(ff*pCO2atm-k0*FugFac*pCO2ocean)
895c6145db Oliv*0373            FluxCO2(i,j) = Kwexch(i,j)*(
14f3186978 Oliv*0374      &      ff(i,j,bi,bj)*apCO2sat(i,j) - fugCO2(i,j)*ak0(i,j,bi,bj)
                0375      &     )
                0376           ELSE
3a1d328575 Oliv*0377            apCO2sat(i,j) = 0. _d 0
                0378            fugCO2(i,j) = 0. _d 0
895c6145db Oliv*0379            FluxCO2(i,j) = 0. _d 0
14f3186978 Oliv*0380           ENDIF
8fbfd1f382 Oliv*0381 C convert flux (mol kg-1 m s-1) to (mmol m-2 s-1)
895c6145db Oliv*0382           FluxCO2(i,j) = FluxCO2(i,j)/m3perkg/ptr2mol
d61a0bd101 Oliv*0383 #ifdef DARWIN_ALLOW_CONS
14f3186978 Oliv*0384           carbSfcFlx(i,j,bi,bj) = carbSfcFlx(i,j,bi,bj)
895c6145db Oliv*0385      &                          + FluxCO2(i,j)*dTsub(ks)
d61a0bd101 Oliv*0386 #endif
8fbfd1f382 Oliv*0387
                0388 #ifdef ALLOW_OLD_VIRTUALFLUX
14f3186978 Oliv*0389           IF (maskC(i,j,ks,bi,bj).NE.0. _d 0) THEN
8fbfd1f382 Oliv*0390 c calculate virtual flux
                0391 c EminusPforV = dS/dt*(1/Sglob)
                0392 C NOTE: Be very careful with signs here!
                0393 C Positive EminusPforV => loss of water to atmos and increase
                0394 C in salinity. Thus, also increase in other surface tracers
                0395 C (i.e. positive virtual flux into surface layer)
                0396 C ...so here, VirtualFLux = dC/dt!
895c6145db Oliv*0397            VFluxCO2(i,j)=gsm_DIC*surfaceForcingS(i,j,bi,bj)/gsm_S
d95627305e Oliv*0398 #ifdef DARWIN_ALLOW_CONS
14f3186978 Oliv*0399            carbVirFlx(i,j,bi,bj) = carbVirFlx(i,j,bi,bj)
895c6145db Oliv*0400      &                         + VFluxCO2(i,j)*dTsub(ks)
d95627305e Oliv*0401 #endif
8fbfd1f382 Oliv*0402 c OR
                0403 c let virtual flux be zero
                0404 c        VirtualFlux(i,j)=0.d0
                0405 c
14f3186978 Oliv*0406           ELSE
895c6145db Oliv*0407            VFluxCO2(i,j)=0. _d 0
14f3186978 Oliv*0408           ENDIF
8fbfd1f382 Oliv*0409 #endif /* ALLOW_OLD_VIRTUALFLUX */
14f3186978 Oliv*0410          ENDDO
                0411         ENDDO
8fbfd1f382 Oliv*0412
                0413 C update tendency
14f3186978 Oliv*0414         DO j=jmin,jmax
                0415          DO i=imin,imax
                0416           IF ( maskC(i,j,ks,bi,bj).NE.0. _d 0 ) THEN
                0417            gDIC(i,j) =
                0418      &       recip_drF(ks)*recip_hFacC(i,j,ks,bi,bj)*
895c6145db Oliv*0419      &                ( FluxCO2(i,j)
8fbfd1f382 Oliv*0420 #ifdef ALLOW_OLD_VIRTUALFLUX
895c6145db Oliv*0421      &                + VFluxCO2(i,j)
8fbfd1f382 Oliv*0422 #endif
14f3186978 Oliv*0423      &                )
                0424           ENDIF
                0425          ENDDO
                0426         ENDDO
8fbfd1f382 Oliv*0427
d8a8f9177f Oliv*0428 #ifdef ALLOW_DIAGNOSTICS
14f3186978 Oliv*0429         IF (useDiagnostics) THEN
                0430          CALL DIAGNOSTICS_FILL(atmospCO2,'apCO2   ',0,1,1,bi,bj,myThid)
895c6145db Oliv*0431          CALL DIAGNOSTICS_FILL(apCO2sat, 'apCO2sat',0,1,2,bi,bj,myThid)
                0432          CALL DIAGNOSTICS_FILL(fugf,     'fugfCO2 ',0,1,1,bi,bj,myThid)
                0433          CALL DIAGNOSTICS_FILL(fugCO2,   'fCO2    ',0,1,2,bi,bj,myThid)
                0434          CALL DIAGNOSTICS_FILL(FluxO2,   'fluxO2  ',0,1,2,bi,bj,myThid)
                0435          CALL DIAGNOSTICS_FILL(FluxCO2,  'fluxCO2 ',0,1,2,bi,bj,myThid)
ad941d6bd9 Oliv*0436 #ifdef ALLOW_OLD_VIRTUALFLUX
895c6145db Oliv*0437          CALL DIAGNOSTICS_FILL(VFluxCO2,'VfluxCO2',0,1,2,bi,bj,myThid)
                0438          CALL DIAGNOSTICS_FILL(VFluxAlk,'VfluxAlk',0,1,2,bi,bj,myThid)
d8a8f9177f Oliv*0439 #endif
14f3186978 Oliv*0440         ENDIF
ad941d6bd9 Oliv*0441 #endif
bdeda14b59 Oliv*0442
14f3186978 Oliv*0443 C      k is at surface
                0444        ENDIF
8fbfd1f382 Oliv*0445
895c6145db Oliv*0446 #ifdef ALLOW_DIAGNOSTICS
                0447        IF (useDiagnostics) THEN
2e7c9eb2a6 Oliv*0448         CALL DIAGNOSTICS_FILL(pCO2,       'pCO2    ',k,1,2,bi,bj,myThid)
                0449         CALL DIAGNOSTICS_FILL(CO3,        'CO3     ',k,1,2,bi,bj,myThid)
                0450         CALL DIAGNOSTICS_FILL(Ksp_TP_Calc,'KspTPClc',k,1,1,bi,bj,myThid)
895c6145db Oliv*0451        ENDIF
                0452 #endif
                0453
8fbfd1f382 Oliv*0454 C === k ================================================================
                0455       ENDDO
                0456 C ======================================================================
                0457
895c6145db Oliv*0458 #ifdef ALLOW_DIAGNOSTICS
                0459       IF (useDiagnostics) THEN
2e7c9eb2a6 Oliv*0460        CALL DIAGNOSTICS_FILL(pH,    'pH      ',0,Nr,1,bi,bj,myThid)
                0461        CALL DIAGNOSTICS_FILL(omegaC,'OmegaC  ',0,Nr,1,bi,bj,myThid)
895c6145db Oliv*0462       ENDIF
                0463 #endif
                0464
8fbfd1f382 Oliv*0465 #endif /* DARWIN_ALLOW_CARBON */
                0466 #endif /* ALLOW_DARWIN */
                0467
                0468       RETURN
                0469       END