pkg/dic/calcite_saturation.F

8e447d7023 Step*0001 #include "DIC_OPTIONS.h"
                0002
                0003 CBOP
5c4b83b693 Jean*0004 C !ROUTINE: CALCITE_SATURATION
8e447d7023 Step*0005
                0006 C !INTERFACE: ==========================================================
5c4b83b693 Jean*0007       SUBROUTINE CALCITE_SATURATION( PTR_DIC, PTR_ALK, PTR_PO4,
8e447d7023 Step*0008      I           bi,bj,imin,imax,jmin,jmax,
2e3e8c330d Jona*0009      I           myTime, myIter, myThid )
8e447d7023 Step*0010
                0011 C !DESCRIPTION:
2e3e8c330d Jona*0012 C  Calculate carbonate fluxes:
                0013 C   - determine carbonate ion concentration through full domain
                0014 C   - calculate calcite saturation state
8e447d7023 Step*0015
                0016 C !USES: ===============================================================
                0017       IMPLICIT NONE
                0018 #include "SIZE.h"
                0019 #include "EEPARAMS.h"
                0020 #include "PARAMS.h"
                0021 #include "GRID.h"
5c4b83b693 Jean*0022 #include "DYNVARS.h"
2ef8966791 Davi*0023 #include "DIC_VARS.h"
8e447d7023 Step*0024
                0025 C !INPUT PARAMETERS: ===================================================
                0026 C  myTime               :: current time
2e3e8c330d Jona*0027 C  myIter               :: current timestep
5c4b83b693 Jean*0028 C  myThid               :: thread number
720e9330bd Step*0029        _RL  PTR_DIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0030        _RL  PTR_ALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0031        _RL  PTR_PO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
5c4b83b693 Jean*0032       INTEGER imin, imax, jmin, jmax, bi, bj
8e447d7023 Step*0033       _RL myTime
2e3e8c330d Jona*0034       INTEGER myIter
8e447d7023 Step*0035       INTEGER myThid
                0036
                0037 C !OUTPUT PARAMETERS: ===================================================
                0038
2e3e8c330d Jona*0039 #if ( defined DIC_BIOTIC && defined DIC_CALCITE_SAT )
8e447d7023 Step*0040
                0041 C !LOCAL VARIABLES: ====================================================
                0042 C  i,j,k                  :: loop indices
5c4b83b693 Jean*0043        INTEGER i,j,k
ae2be6150b Jona*0044        LOGICAL doIni_pH, debugPrt
a1d0e455fd Hann*0045        _RL carbonate(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8e447d7023 Step*0046        _RL calcium
a1d0e455fd Hann*0047        _RL sitlocal
8e447d7023 Step*0048        _RL po4local
                0049        _RL diclocal
                0050        _RL alklocal
a1d0e455fd Hann*0051        _RL pCO2local(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
5c4b83b693 Jean*0052        _RL locTemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0053        _RL locSalt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
ae2be6150b Jona*0054        INTEGER doIterCO3, CO3ITER
8e447d7023 Step*0055 CEOP
                0056
ae2be6150b Jona*0057        doIni_pH = ( myIter.EQ.(nIter0 + 1) ) .AND. .NOT.pH_isLoaded(2)
                0058 C-   set number of iterations for 3-D pH (for case selectPHsolver=0):
                0059        IF ( doIni_pH .OR.
                0060      &      ( calcOmegaCalciteFreq.GT.deltaTClock ) ) THEN
                0061 C    Since we start cold with an arbitrary pH, we must iterate pH solver
                0062 C    to ensure accurate estimate of CO3 at depth
                0063 C    OR we may call this for deep ocean infrequently and in this case
                0064 C    can afford to make several iterations
                0065          doIterCO3 = nIterCO3
                0066        ELSE
                0067          doIterCO3 = 1
                0068        ENDIF
                0069 c      IF ( bi+bj .EQ. 2 ) THEN
                0070 c       WRITE(standardMessageUnit,'(A,L5,I4,I10)')
                0071 c    &  '--- CALCITE_SATURATION: doIni_pH, doIterCO3, myIter=',
                0072 c    &                           doIni_pH, doIterCO3, myIter
                0073 c      ENDIF
                0074
5c4b83b693 Jean*0075        DO k=1,Nr
                0076         DO j=1-OLy,sNy+OLy
                0077          DO i=1-OLx,sNx+OLx
                0078             locTemp(i,j) = theta(i,j,k,bi,bj)
                0079             locSalt(i,j) = salt (i,j,k,bi,bj)
a1d0e455fd Hann*0080             carbonate(i,j) = 0. _d 0
                0081             pCO2local(i,j) = 0. _d 0
6acab690ae Jona*0082          ENDDO
                0083         ENDDO
7459d17a4c Jona*0084 #ifdef CARBONCHEM_SOLVESAPHE
6acab690ae Jona*0085 #ifdef ALLOW_DEBUG
a1d0e455fd Hann*0086         IF (debugMode) CALL DEBUG_CALL( 'DIC_COEFFS_DEEP',myThid )
6acab690ae Jona*0087 #endif
                0088 C Calculate carbon coefficients
                0089         CALL DIC_COEFFS_SURF(
5c4b83b693 Jean*0090      I                       locTemp, locSalt,
6acab690ae Jona*0091      I                       bi,bj,iMin,iMax,jMin,jMax,myThid)
                0092
                0093 C Now correct the coefficients for pressure dependence
                0094         CALL DIC_COEFFS_DEEP(
5c4b83b693 Jean*0095      I                       locTemp, locSalt,
6acab690ae Jona*0096      I                       bi,bj,iMin,iMax,jMin,jMax,
                0097      I                       k,myThid)
                0098 #else
                0099 #ifdef ALLOW_DEBUG
a1d0e455fd Hann*0100         IF (debugMode) CALL DEBUG_CALL(
                0101      &      'CARBON_COEFFS_PRESSURE_DEP',myThid)
6acab690ae Jona*0102 #endif
8e447d7023 Step*0103         CALL CARBON_COEFFS_PRESSURE_DEP(
5c4b83b693 Jean*0104      I                       locTemp, locSalt,
8e447d7023 Step*0105      I                       bi,bj,iMin,iMax,jMin,jMax,
5625485478 Jean*0106      I                       k,myThid)
6acab690ae Jona*0107 #endif
8e447d7023 Step*0108
6acab690ae Jona*0109         debugPrt = debugMode
51c3bf0077 Step*0110         DO j=jmin,jmax
                0111          DO i=imin,imax
8e447d7023 Step*0112
a1d0e455fd Hann*0113           IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
6acab690ae Jona*0114 #ifdef CARBONCHEM_SOLVESAPHE
2e3e8c330d Jona*0115 C calcium concentration already calculated in S/R DIC_COEFFS_SURF
a1d0e455fd Hann*0116            calcium = cat(i,j,bi,bj)
6acab690ae Jona*0117 #else
a1d0e455fd Hann*0118            calcium = 1.028 _d -2*salt(i,j,k,bi,bj)/35. _d 0
6acab690ae Jona*0119 #endif
8e447d7023 Step*0120
2e3e8c330d Jona*0121 C SilicaDeep filled with constant value of 0.03 mol/m3 unless
                0122 C     DIC_deepSilicaFile provided.
a1d0e455fd Hann*0123            sitlocal = silicaDeep(i,j,k,bi,bj)
                0124            po4local = PTR_PO4(i,j,k)
                0125            diclocal = PTR_DIC(i,j,k)
                0126            alklocal = PTR_ALK(i,j,k)
6acab690ae Jona*0127
                0128 #ifdef CARBONCHEM_SOLVESAPHE
a1d0e455fd Hann*0129            IF ( selectPHsolver.GT.0 ) THEN
6acab690ae Jona*0130 C Use Munhoven (2013) Solvesaphe routine to calculate pH and pCO2
ae2be6150b Jona*0131             IF ( doIni_pH ) THEN
6acab690ae Jona*0132 #ifdef ALLOW_DEBUG
ae2be6150b Jona*0133              IF (debugPrt) CALL DEBUG_CALL('AHINI_FOR_AT',myThid)
6acab690ae Jona*0134 #endif
                0135 C call AHINI_FOR_AT to get better initial guess of pH
ae2be6150b Jona*0136              CALL AHINI_FOR_AT(alklocal*permil,
                0137      I              diclocal*permil,
                0138      I              bt(i,j,bi,bj),
                0139      O              pH3d(i,j,k,bi,bj),
                0140      I              i,j,k,bi,bj,myIter,myThid )
                0141             ENDIF /* 3d pH is not loaded */
6acab690ae Jona*0142 #ifdef ALLOW_DEBUG
a1d0e455fd Hann*0143             IF (debugPrt) CALL DEBUG_CALL('CALC_PCO2_SOLVESAPHE',myThid)
6acab690ae Jona*0144 #endif
                0145             CALL CALC_PCO2_SOLVESAPHE(
5c4b83b693 Jean*0146      I          locTemp(i,j), locSalt(i,j),
6acab690ae Jona*0147      I          diclocal, po4local,
a1d0e455fd Hann*0148      I          sitlocal, alklocal,
ae2be6150b Jona*0149      U          pH3d(i,j,k,bi,bj), pCO2local(i,j), carbonate(i,j),
6acab690ae Jona*0150      I          i,j,k,bi,bj, debugPrt,myIter,myThid )
6b3f5dbce1 Hann*0151
                0152 C- convert carbonate to mol kg^-1-SW for calculation of saturation state
a1d0e455fd Hann*0153             carbonate(i,j) = carbonate(i,j)*permil
                0154            ELSE
6acab690ae Jona*0155 C Use the original Follows et al. (2006) solver
                0156 #endif /* CARBONCHEM_SOLVESAPHE */
                0157 #ifdef ALLOW_DEBUG
                0158             IF (debugPrt) CALL DEBUG_CALL('CALC_PCO2_APPROX',myThid)
                0159 #endif
ae2be6150b Jona*0160             DO CO3iter = 1, doIterCO3
                0161              CALL CALC_PCO2_APPROX(
5c4b83b693 Jean*0162      I          locTemp(i,j), locSalt(i,j),
8e447d7023 Step*0163      I          diclocal, po4local,
a1d0e455fd Hann*0164      I          sitlocal, alklocal,
8e447d7023 Step*0165      I          ak1(i,j,bi,bj),ak2(i,j,bi,bj),
                0166      I          ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
                0167      I          aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
e6a52874f7 Davi*0168      I          aksi(i,j,bi,bj),akf(i,j,bi,bj),
                0169      I          ak0(i,j,bi,bj), fugf(i,j,bi,bj), ff(i,j,bi,bj),
8e447d7023 Step*0170      I          bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
ae2be6150b Jona*0171      U          pH3d(i,j,k,bi,bj), pCO2local(i,j), carbonate(i,j),
e098791e51 Jean*0172      I          i,j,k,bi,bj,myIter,myThid )
5c4b83b693 Jean*0173             ENDDO
ae2be6150b Jona*0174
6acab690ae Jona*0175 #ifdef CARBONCHEM_SOLVESAPHE
a1d0e455fd Hann*0176            ENDIF
6acab690ae Jona*0177 #endif /* CARBONCHEM_SOLVESAPHE */
8e447d7023 Step*0178
a1d0e455fd Hann*0179            omegaC(i,j,k,bi,bj) = calcium * carbonate(i,j)
                0180      &                         / Ksp_TP_Calc(i,j,bi,bj)
2e3e8c330d Jona*0181            debugPrt = .FALSE.
8e447d7023 Step*0182
5c4b83b693 Jean*0183 Cmick...................................................
8e447d7023 Step*0184 c            if(omegaC(i,j,k,bi,bj) .eq. 0.) then
                0185 c             if(i .eq. 76 .and. j .eq. 36  .and. k .eq. 15) then
                0186 c               write(6,*)'i,j,k,KS,CO3,pHCa,T,S,hfacc,omega',
                0187 c     &                 i,j,k,
                0188 c     &                 Ksp_TP_Calc(i,j,bi,bj),
ae2be6150b Jona*0189 c     &                 carbonate(i,j), calcium, pH3d(i,j,k,bi,bj),
5c4b83b693 Jean*0190 c     &                 locTemp(i,j), locSalt(i,j),
8e447d7023 Step*0191 c     &                 hfacc(i,j,k,bi,bj),omegaC(i,j,k,bi,bj)
                0192 c              write(6,*)'Ksp_TP_Calc',
                0193 c     &                 Ksp_TP_Calc(i,j,bi,bj)
                0194 c               write(6,*)'dic, alk, po4 ',
                0195 c     &                 diclocal, alklocal,po4local
                0196 c               write(6,*)'k1, k2, k1p, k2p, k3p ',
                0197 c     &                 ak1(i,j,bi,bj),ak2(i,j,bi,bj),
                0198 c     &                ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj)
                0199 c               write(6,*)'ks, kb, kw, ksi ',
                0200 c     &                aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
                0201 c     &                aksi(i,j,bi,bj)
                0202 c               write(6,*)'akf, ff, bt, st, ft ',
                0203 c     &                akf(i,j,bi,bj),ff(i,j,bi,bj),
                0204 c     &                bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj)
                0205 c             end if
5c4b83b693 Jean*0206 Cmick....................................................
a1d0e455fd Hann*0207           ELSE
                0208 C-        else (case maskC = 0 , dry point)
                0209            omegaC(i,j,k,bi,bj) = 0. _d 0
                0210           ENDIF
8e447d7023 Step*0211
                0212          ENDDO
                0213         ENDDO
                0214
a1d0e455fd Hann*0215 #ifdef ALLOW_DIAGNOSTICS
ae2be6150b Jona*0216 C Fill up pCO2 and carbonate diagnostics from local arrays
a1d0e455fd Hann*0217         IF ( useDiagnostics ) THEN
                0218          CALL DIAGNOSTICS_FILL(pCO2local,'DIC3DPCO',k,1,2,bi,bj,myThid)
                0219          CALL DIAGNOSTICS_FILL(carbonate,'DIC3DCO3',k,1,2,bi,bj,myThid)
                0220         ENDIF
                0221 #endif /* ALLOW_DIAGNOSTICS */
                0222
5c4b83b693 Jean*0223 C-     end k loop
8e447d7023 Step*0224        ENDDO
5c4b83b693 Jean*0225
a1d0e455fd Hann*0226 #ifdef ALLOW_DIAGNOSTICS
ae2be6150b Jona*0227 C Fill up 3d pH, omegaC and silicate diagnostics
a1d0e455fd Hann*0228        IF ( useDiagnostics ) THEN
ae2be6150b Jona*0229         CALL DIAGNOSTICS_FILL(pH3d      ,'DIC3DPH ',0,Nr,1,bi,bj,myThid)
a1d0e455fd Hann*0230         CALL DIAGNOSTICS_FILL(omegaC    ,'OMEGAC  ',0,Nr,1,bi,bj,myThid)
                0231         CALL DIAGNOSTICS_FILL(silicaDeep,'DIC3DSIT',0,Nr,1,bi,bj,myThid)
                0232        ENDIF
                0233 #endif /* ALLOW_DIAGNOSTICS */
                0234
2e3e8c330d Jona*0235 #endif /* DIC_BIOTIC and DIC_CALCITE_SAT */
8e447d7023 Step*0236        RETURN
                0237        END