pkg/bling/bling_carbonate_sys.F

c0d1c06c15 Matt*0001 #include "BLING_OPTIONS.h"
b26a461de7 Mart*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
c0d1c06c15 Matt*0005
                0006 CBOP
15ec8028f7 aver*0007       SUBROUTINE BLING_CARBONATE_SYS(
c0d1c06c15 Matt*0008      I           PTR_DIC, PTR_ALK, PTR_PO4,
e0f9a7ba0b Matt*0009 #ifdef USE_SIBLING
                0010      I           PTR_SI,
                0011 #endif
15ec8028f7 aver*0012      I           bi, bj, iMin, iMax, jMin, jMax,
                0013      I           myTime, myIter, myThid )
c0d1c06c15 Matt*0014
                0015 C     =================================================================
                0016 C     | subroutine bling_carbonate_sys
                0017 C     | o Calculate carbonate fluxes
                0018 C     |   Also update pH (3d field)
                0019 C     =================================================================
                0020
e0f9a7ba0b Matt*0021       IMPLICIT NONE
                0022
c0d1c06c15 Matt*0023 C     == GLobal variables ==
                0024 #include "SIZE.h"
                0025 #include "DYNVARS.h"
                0026 #include "EEPARAMS.h"
                0027 #include "PARAMS.h"
                0028 #include "GRID.h"
                0029 #include "BLING_VARS.h"
b26a461de7 Mart*0030 #ifdef ALLOW_AUTODIFF_TAMC
                0031 # include "tamc.h"
                0032 #endif
c0d1c06c15 Matt*0033
                0034 C     == Routine arguments ==
                0035 C     PTR_DIC              :: dissolved inorganic carbon
                0036 C     PTR_ALK              :: alkalinity
                0037 C     PTR_PO4              :: phosphate
                0038 C     myTime               :: current time
15ec8028f7 aver*0039 C     myIter               :: current timestep
                0040 C     myThid               :: thread Id. number
c0d1c06c15 Matt*0041       _RL  PTR_DIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0042       _RL  PTR_ALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0043       _RL  PTR_PO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
e0f9a7ba0b Matt*0044 #ifdef USE_SIBLING
                0045       _RL  PTR_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0046 #endif
15ec8028f7 aver*0047       INTEGER bi, bj, iMin, iMax, jMin, jMax
                0048       _RL  myTime
                0049       INTEGER myIter
                0050       INTEGER myThid
c0d1c06c15 Matt*0051
                0052 #ifdef ALLOW_PTRACERS
                0053
                0054 C     == Local variables ==
                0055 C     i,j,k             :: loop indices
                0056 C     carbonate         :: local value of calcium carbonate
                0057 C     calcium           :: local value of Ca
15ec8028f7 aver*0058 C     sitlocal          :: local value of Si
c0d1c06c15 Matt*0059 C     diclocal          :: local value of DIC
                0060 C     alklocal          :: local value of ALK
                0061 C     pCO2local         :: local value of pCO2
                0062 C     pHlocal           :: local value of pH
e0f9a7ba0b Matt*0063 C     CO3iter           :: iterations counter for CO3 ion calculation
                0064 C     CO3iterMax        :: total number of iterations
c0d1c06c15 Matt*0065        INTEGER i,j,k
                0066        _RL carbonate
                0067        _RL calcium
                0068        _RL po4local
15ec8028f7 aver*0069        _RL sitlocal
c0d1c06c15 Matt*0070        _RL diclocal
                0071        _RL alklocal
                0072        _RL pCO2local
                0073        _RL pHlocal
15ec8028f7 aver*0074        _RL locTemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0075        _RL locSalt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c0d1c06c15 Matt*0076
e0f9a7ba0b Matt*0077 c      INTEGER CO3iter
                0078 c      INTEGER CO3iterMax
b26a461de7 Mart*0079 #ifdef ALLOW_AUTODIFF_TAMC
                0080 C     tkey :: tape key (tile dependent)
                0081 C     kkey :: tape key (level and tile dependent)
                0082 C     ikey :: tape key (position, level, and tile dependent)
                0083       INTEGER tkey, kkey, ikey
                0084 #endif
c0d1c06c15 Matt*0085 CEOP
                0086
                0087 C  Since pH is now a 3D field and is solved for at every time step
                0088 C  few iterations are needed
e0f9a7ba0b Matt*0089 c      CO3iterMax = 1
c0d1c06c15 Matt*0090
                0091 C determine carbonate ion concentration through full domain
                0092 C determine calcite saturation state
                0093
b26a461de7 Mart*0094 #ifdef ALLOW_AUTODIFF_TAMC
                0095       tkey = bi + (bj - 1)*nSx + (ikey_dynamics - 1)*nSx*nSy
                0096 #endif
c0d1c06c15 Matt*0097 C$TAF LOOP = parallel
                0098        DO k=1,Nr
                0099
b26a461de7 Mart*0100 #ifdef ALLOW_AUTODIFF_TAMC
                0101         kkey = k + (tkey - 1)*Nr
                0102 #endif
15ec8028f7 aver*0103         DO j=jMin,jMax
                0104          DO i=iMin,iMax
                0105           locTemp(i,j) = theta(i,j,k,bi,bj)
                0106           locSalt(i,j) = salt (i,j,k,bi,bj)
                0107          ENDDO
                0108         ENDDO
c0d1c06c15 Matt*0109
e0f9a7ba0b Matt*0110 #ifdef CARBONCHEM_SOLVESAPHE
                0111 #ifdef ALLOW_DEBUG
                0112       IF (debugMode) CALL DEBUG_CALL(
                0113      &   'DIC_COEFFS_DEEP',myThid)
                0114 #endif
                0115 C Calculate carbon coefficients
                0116         CALL DIC_COEFFS_SURF(
15ec8028f7 aver*0117      I                       locTemp, locSalt,
e0f9a7ba0b Matt*0118      I                       bi,bj,iMin,iMax,jMin,jMax,myThid)
                0119
                0120 C Now correct the coefficients for pressure dependence
                0121         CALL DIC_COEFFS_DEEP(
15ec8028f7 aver*0122      I                       locTemp, locSalt,
e0f9a7ba0b Matt*0123      I                       bi,bj,iMin,iMax,jMin,jMax,
                0124      I                       k,myThid)
                0125
15ec8028f7 aver*0126 #else /* CARBONCHEM_SOLVESAPHE */
e0f9a7ba0b Matt*0127
                0128 #ifdef ALLOW_DEBUG
                0129       IF (debugMode) CALL DEBUG_CALL(
                0130      &   'CARBON_COEFFS_PRESSURE_DEP',myThid)
                0131 #endif
                0132
c0d1c06c15 Matt*0133 C  Get coefficients for carbonate calculations
15ec8028f7 aver*0134         CALL CARBON_COEFFS_PRESSURE_DEP(
                0135      I                       locTemp, locSalt,
                0136      I                       bi, bj, iMin, iMax, jMin, jMax,
                0137      I                       k, myThid )
                0138 #endif /* CARBONCHEM_SOLVESAPHE */
c0d1c06c15 Matt*0139
                0140 C--------------------------------------------------
                0141
                0142 C$TAF LOOP = parallel
15ec8028f7 aver*0143         DO j=jMin,jMax
c0d1c06c15 Matt*0144 C$TAF LOOP = parallel
15ec8028f7 aver*0145          DO i=iMin,iMax
c0d1c06c15 Matt*0146
15ec8028f7 aver*0147           IF ( hFacC(i,j,k,bi,bj) .GT. 0. _d 0) THEN
c0d1c06c15 Matt*0148
b26a461de7 Mart*0149 #ifdef ALLOW_AUTODIFF_TAMC
                0150            ikey = i + (j - 1)*sNx + (kkey - 1)*sNx*sNy
                0151 #endif
e0f9a7ba0b Matt*0152 #ifdef CARBONCHEM_SOLVESAPHE
15ec8028f7 aver*0153            calcium = cat(i,j,bi,bj)
e0f9a7ba0b Matt*0154 #else
c0d1c06c15 Matt*0155 C  Estimate calcium concentration from salinity
15ec8028f7 aver*0156            calcium = 1.028 _d -2*salt(i,j,k,bi,bj)/35. _d 0
e0f9a7ba0b Matt*0157 #endif
c0d1c06c15 Matt*0158
15ec8028f7 aver*0159            po4local = PTR_PO4(i,j,k)
                0160            diclocal = PTR_DIC(i,j,k)
                0161            alklocal = PTR_ALK(i,j,k)
                0162            pHlocal  = pH(i,j,k,bi,bj)
e0f9a7ba0b Matt*0163
                0164 C  Assume constant deep silica value
                0165 C  30 micromol = 0.03 mol m-3
                0166 C  unless SIBLING is used
                0167 #ifdef USE_SIBLING
15ec8028f7 aver*0168            sitlocal = PTR_SI(i,j,k)
e0f9a7ba0b Matt*0169 #else
15ec8028f7 aver*0170            sitlocal = 0.03 _d 0
e0f9a7ba0b Matt*0171 #endif
                0172
                0173 #ifdef CARBONCHEM_SOLVESAPHE
15ec8028f7 aver*0174            IF ( selectPHsolver.GT.0 ) THEN
e0f9a7ba0b Matt*0175 C Use Munhoven (2013) Solvesaphe routine to calculate pH and pCO2
                0176 #ifdef ALLOW_DEBUG
                0177          IF (debugMode) CALL DEBUG_CALL('AHINI_FOR_AT',myThid)
                0178 #endif
                0179 CAV since we carry pH, no need for an initial guess
                0180 C call AHINI_FOR_AT to get better initial guess of pH
                0181 c               CALL AHINI_FOR_AT(alklocal*permil,
                0182 c     I           diclocal*permil,
                0183 c     I           bt(i,j,bi,bj),
                0184 c     O           pHlocal,
                0185 c     I           i,j,k,bi,bj,myIter,myThid )
                0186 #ifdef ALLOW_DEBUG
                0187          IF (debugMode) CALL DEBUG_CALL('CALC_PCO2_SOLVESAPHE',myThid)
                0188 #endif
                0189             CALL CALC_PCO2_SOLVESAPHE(
15ec8028f7 aver*0190      I          locTemp(i,j), locSalt(i,j),
e0f9a7ba0b Matt*0191      I          diclocal, po4local,
15ec8028f7 aver*0192      I          sitlocal, alklocal,
                0193      U          pHlocal, pCO2local, carbonate,
e0f9a7ba0b Matt*0194      I          i,j,k,bi,bj,myIter,myThid )
15ec8028f7 aver*0195
                0196 C- convert carbonate to mol kg^-1-SW for calculation of saturation state
                0197             carbonate = carbonate*permil
                0198            ELSE
e0f9a7ba0b Matt*0199 C Use the original Follows et al. (2006) solver
                0200 #endif /* CARBONCHEM_SOLVESAPHE */
                0201 #ifdef ALLOW_DEBUG
                0202             IF (debugMode) CALL DEBUG_CALL('CALC_PCO2_APPROX',myThid)
                0203 #endif
c0d1c06c15 Matt*0204
                0205 C  Evaluate carbonate (CO3) ions concentration
                0206 C  iteratively
                0207
15ec8028f7 aver*0208 c           DO CO3iter = 1, CO3iterMax
c0d1c06c15 Matt*0209
                0210 C--------------------------------------------------
                0211
15ec8028f7 aver*0212              CALL CALC_PCO2_APPROX(
                0213      I          locTemp(i,j), locSalt(i,j),
c0d1c06c15 Matt*0214      I          diclocal, po4local,
15ec8028f7 aver*0215      I          sitlocal,alklocal,
c0d1c06c15 Matt*0216      I          ak1(i,j,bi,bj),ak2(i,j,bi,bj),
                0217      I          ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
                0218      I          aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
                0219      I          aksi(i,j,bi,bj),akf(i,j,bi,bj),
e0f9a7ba0b Matt*0220      I          ak0(i,j,bi,bj),fugf(i,j,bi,bj),ff(i,j,bi,bj),
c0d1c06c15 Matt*0221      I          bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
                0222      U          pHlocal,pCO2local,carbonate,
                0223      I          i,j,k,bi,bj,myIter,myThid )
15ec8028f7 aver*0224 c           ENDDO
c0d1c06c15 Matt*0225
e0f9a7ba0b Matt*0226 #ifdef CARBONCHEM_SOLVESAPHE
15ec8028f7 aver*0227            ENDIF
e0f9a7ba0b Matt*0228 #endif /* CARBONCHEM_SOLVESAPHE */
b26a461de7 Mart*0229 #ifdef ALLOW_AUTODIFF_TAMC
                0230 CADJ store carbonate = comlev1_bibj_ijk, key = ikey, kind = isbyte
                0231 #endif
15ec8028f7 aver*0232            pH(i,j,k,bi,bj) = pHlocal
c0d1c06c15 Matt*0233
e0f9a7ba0b Matt*0234 C  Calculate calcium carbonate (calcite and aragonite)
c0d1c06c15 Matt*0235 C  saturation state
15ec8028f7 aver*0236            omegaC(i,j,k,bi,bj) = calcium * carbonate /
c0d1c06c15 Matt*0237      &                          Ksp_TP_Calc(i,j,bi,bj)
15ec8028f7 aver*0238            omegaAr(i,j,k,bi,bj) = calcium * carbonate /
c0d1c06c15 Matt*0239      &                          Ksp_TP_Arag(i,j,bi,bj)
                0240
15ec8028f7 aver*0241           ELSE
c0d1c06c15 Matt*0242
15ec8028f7 aver*0243            pH(i,j,k,bi,bj) = 8. _d 0
                0244            omegaC(i,j,k,bi,bj)  = 0. _d 0
                0245            omegaAr(i,j,k,bi,bj) = 0. _d 0
c0d1c06c15 Matt*0246
15ec8028f7 aver*0247           ENDIF
c0d1c06c15 Matt*0248
                0249          ENDDO
                0250         ENDDO
                0251        ENDDO
                0252
                0253 #endif /* ALLOW_PTRACERS */
                0254        RETURN
                0255        END