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
6acab690ae Jona*0044        LOGICAL 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)
                0052        _RL pHlocal(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
5c4b83b693 Jean*0053        _RL locTemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0054        _RL locSalt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
7459d17a4c Jona*0055
8e447d7023 Step*0056        INTEGER CO3ITER
                0057 CEOP
                0058
5c4b83b693 Jean*0059        DO k=1,Nr
                0060         DO j=1-OLy,sNy+OLy
                0061          DO i=1-OLx,sNx+OLx
                0062             locTemp(i,j) = theta(i,j,k,bi,bj)
                0063             locSalt(i,j) = salt (i,j,k,bi,bj)
a1d0e455fd Hann*0064             carbonate(i,j) = 0. _d 0
                0065             pCO2local(i,j) = 0. _d 0
                0066             pHlocal(i,j)   = 0. _d 0
6acab690ae Jona*0067          ENDDO
                0068         ENDDO
7459d17a4c Jona*0069 #ifdef CARBONCHEM_SOLVESAPHE
6acab690ae Jona*0070 #ifdef ALLOW_DEBUG
a1d0e455fd Hann*0071         IF (debugMode) CALL DEBUG_CALL( 'DIC_COEFFS_DEEP',myThid )
6acab690ae Jona*0072 #endif
                0073 C Calculate carbon coefficients
                0074         CALL DIC_COEFFS_SURF(
5c4b83b693 Jean*0075      I                       locTemp, locSalt,
6acab690ae Jona*0076      I                       bi,bj,iMin,iMax,jMin,jMax,myThid)
                0077
                0078 C Now correct the coefficients for pressure dependence
                0079         CALL DIC_COEFFS_DEEP(
5c4b83b693 Jean*0080      I                       locTemp, locSalt,
6acab690ae Jona*0081      I                       bi,bj,iMin,iMax,jMin,jMax,
                0082      I                       k,myThid)
                0083 #else
                0084 #ifdef ALLOW_DEBUG
a1d0e455fd Hann*0085         IF (debugMode) CALL DEBUG_CALL(
                0086      &      'CARBON_COEFFS_PRESSURE_DEP',myThid)
6acab690ae Jona*0087 #endif
8e447d7023 Step*0088         CALL CARBON_COEFFS_PRESSURE_DEP(
5c4b83b693 Jean*0089      I                       locTemp, locSalt,
8e447d7023 Step*0090      I                       bi,bj,iMin,iMax,jMin,jMax,
5625485478 Jean*0091      I                       k,myThid)
6acab690ae Jona*0092 #endif
8e447d7023 Step*0093
6acab690ae Jona*0094         debugPrt = debugMode
51c3bf0077 Step*0095         DO j=jmin,jmax
                0096          DO i=imin,imax
8e447d7023 Step*0097
a1d0e455fd Hann*0098           IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
6acab690ae Jona*0099 #ifdef CARBONCHEM_SOLVESAPHE
2e3e8c330d Jona*0100 C calcium concentration already calculated in S/R DIC_COEFFS_SURF
a1d0e455fd Hann*0101            calcium = cat(i,j,bi,bj)
6acab690ae Jona*0102 #else
a1d0e455fd Hann*0103            calcium = 1.028 _d -2*salt(i,j,k,bi,bj)/35. _d 0
6acab690ae Jona*0104 #endif
8e447d7023 Step*0105
2e3e8c330d Jona*0106 C SilicaDeep filled with constant value of 0.03 mol/m3 unless
                0107 C     DIC_deepSilicaFile provided.
a1d0e455fd Hann*0108            sitlocal = silicaDeep(i,j,k,bi,bj)
                0109            po4local = PTR_PO4(i,j,k)
                0110            diclocal = PTR_DIC(i,j,k)
                0111            alklocal = PTR_ALK(i,j,k)
6acab690ae Jona*0112
                0113 #ifdef CARBONCHEM_SOLVESAPHE
a1d0e455fd Hann*0114            IF ( selectPHsolver.GT.0 ) THEN
6acab690ae Jona*0115 C Use Munhoven (2013) Solvesaphe routine to calculate pH and pCO2
                0116 #ifdef ALLOW_DEBUG
a1d0e455fd Hann*0117             IF (debugPrt) CALL DEBUG_CALL('AHINI_FOR_AT',myThid)
6acab690ae Jona*0118 #endif
                0119 C call AHINI_FOR_AT to get better initial guess of pH
                0120             CALL AHINI_FOR_AT(alklocal*permil,
                0121      I           diclocal*permil,
                0122      I           bt(i,j,bi,bj),
a1d0e455fd Hann*0123      O           pHlocal(i,j),
6acab690ae Jona*0124      I           i,j,k,bi,bj,myIter,myThid )
                0125
                0126 #ifdef ALLOW_DEBUG
a1d0e455fd Hann*0127             IF (debugPrt) CALL DEBUG_CALL('CALC_PCO2_SOLVESAPHE',myThid)
6acab690ae Jona*0128 #endif
                0129             CALL CALC_PCO2_SOLVESAPHE(
5c4b83b693 Jean*0130      I          locTemp(i,j), locSalt(i,j),
6acab690ae Jona*0131      I          diclocal, po4local,
a1d0e455fd Hann*0132      I          sitlocal, alklocal,
                0133      U          pHlocal(i,j), pCO2local(i,j), carbonate(i,j),
6acab690ae Jona*0134      I          i,j,k,bi,bj, debugPrt,myIter,myThid )
6b3f5dbce1 Hann*0135
                0136 C- convert carbonate to mol kg^-1-SW for calculation of saturation state
a1d0e455fd Hann*0137             carbonate(i,j) = carbonate(i,j)*permil
                0138            ELSE
6acab690ae Jona*0139 C Use the original Follows et al. (2006) solver
                0140 #endif /* CARBONCHEM_SOLVESAPHE */
                0141 #ifdef ALLOW_DEBUG
                0142             IF (debugPrt) CALL DEBUG_CALL('CALC_PCO2_APPROX',myThid)
                0143 #endif
a1d0e455fd Hann*0144             pHlocal(i,j) = 7.9 _d 0
                0145
2e3e8c330d Jona*0146 C Since we start cold with an arbitrary pH, we must iterate pH solver
                0147 C    to ensure accurate estimate of co3 at depth
                0148 C    Because we may call this for deep ocean infrequently can afford to make
                0149 C    several iterations
a1d0e455fd Hann*0150
2e3e8c330d Jona*0151             DO CO3iter = 1, nIterCO3
6acab690ae Jona*0152               CALL CALC_PCO2_APPROX(
5c4b83b693 Jean*0153      I          locTemp(i,j), locSalt(i,j),
8e447d7023 Step*0154      I          diclocal, po4local,
a1d0e455fd Hann*0155      I          sitlocal, alklocal,
8e447d7023 Step*0156      I          ak1(i,j,bi,bj),ak2(i,j,bi,bj),
                0157      I          ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
                0158      I          aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
e6a52874f7 Davi*0159      I          aksi(i,j,bi,bj),akf(i,j,bi,bj),
                0160      I          ak0(i,j,bi,bj), fugf(i,j,bi,bj), ff(i,j,bi,bj),
8e447d7023 Step*0161      I          bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
a1d0e455fd Hann*0162      U          pHlocal(i,j), pCO2local(i,j), carbonate(i,j),
e098791e51 Jean*0163      I          i,j,k,bi,bj,myIter,myThid )
5c4b83b693 Jean*0164             ENDDO
6acab690ae Jona*0165 #ifdef CARBONCHEM_SOLVESAPHE
a1d0e455fd Hann*0166            ENDIF
6acab690ae Jona*0167 #endif /* CARBONCHEM_SOLVESAPHE */
8e447d7023 Step*0168
a1d0e455fd Hann*0169            omegaC(i,j,k,bi,bj) = calcium * carbonate(i,j)
                0170      &                         / Ksp_TP_Calc(i,j,bi,bj)
2e3e8c330d Jona*0171            debugPrt = .FALSE.
8e447d7023 Step*0172
5c4b83b693 Jean*0173 Cmick...................................................
8e447d7023 Step*0174 c            if(omegaC(i,j,k,bi,bj) .eq. 0.) then
                0175 c             if(i .eq. 76 .and. j .eq. 36  .and. k .eq. 15) then
                0176 c               write(6,*)'i,j,k,KS,CO3,pHCa,T,S,hfacc,omega',
                0177 c     &                 i,j,k,
                0178 c     &                 Ksp_TP_Calc(i,j,bi,bj),
a1d0e455fd Hann*0179 c     &                 carbonate(i,j), calcium, pHlocal(i,j),
5c4b83b693 Jean*0180 c     &                 locTemp(i,j), locSalt(i,j),
8e447d7023 Step*0181 c     &                 hfacc(i,j,k,bi,bj),omegaC(i,j,k,bi,bj)
                0182 c              write(6,*)'Ksp_TP_Calc',
                0183 c     &                 Ksp_TP_Calc(i,j,bi,bj)
                0184 c               write(6,*)'dic, alk, po4 ',
                0185 c     &                 diclocal, alklocal,po4local
                0186 c               write(6,*)'k1, k2, k1p, k2p, k3p ',
                0187 c     &                 ak1(i,j,bi,bj),ak2(i,j,bi,bj),
                0188 c     &                ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj)
                0189 c               write(6,*)'ks, kb, kw, ksi ',
                0190 c     &                aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
                0191 c     &                aksi(i,j,bi,bj)
                0192 c               write(6,*)'akf, ff, bt, st, ft ',
                0193 c     &                akf(i,j,bi,bj),ff(i,j,bi,bj),
                0194 c     &                bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj)
                0195 c             end if
5c4b83b693 Jean*0196 Cmick....................................................
a1d0e455fd Hann*0197           ELSE
                0198 C-        else (case maskC = 0 , dry point)
                0199            omegaC(i,j,k,bi,bj) = 0. _d 0
                0200           ENDIF
8e447d7023 Step*0201
                0202          ENDDO
                0203         ENDDO
                0204
a1d0e455fd Hann*0205 C Fill up some diagnostics of deep pH, pCO2, and carbonate local arrays
                0206 #ifdef ALLOW_DIAGNOSTICS
                0207         IF ( useDiagnostics ) THEN
                0208          CALL DIAGNOSTICS_FILL(pHlocal  ,'DIC3DPH ',k,1,2,bi,bj,myThid)
                0209          CALL DIAGNOSTICS_FILL(pCO2local,'DIC3DPCO',k,1,2,bi,bj,myThid)
                0210          CALL DIAGNOSTICS_FILL(carbonate,'DIC3DCO3',k,1,2,bi,bj,myThid)
                0211         ENDIF
                0212 #endif /* ALLOW_DIAGNOSTICS */
                0213
5c4b83b693 Jean*0214 C-     end k loop
8e447d7023 Step*0215        ENDDO
5c4b83b693 Jean*0216
a1d0e455fd Hann*0217 #ifdef ALLOW_DIAGNOSTICS
                0218        IF ( useDiagnostics ) THEN
                0219         CALL DIAGNOSTICS_FILL(omegaC    ,'OMEGAC  ',0,Nr,1,bi,bj,myThid)
                0220         CALL DIAGNOSTICS_FILL(silicaDeep,'DIC3DSIT',0,Nr,1,bi,bj,myThid)
                0221        ENDIF
                0222 #endif /* ALLOW_DIAGNOSTICS */
                0223
2e3e8c330d Jona*0224 #endif /* DIC_BIOTIC and DIC_CALCITE_SAT */
8e447d7023 Step*0225        RETURN
                0226        END