pkg/bling/bling_carbonate_sys.F

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