pkg/dic/dic_surfforcing_init.F

e8625f0081 Step*0001 #include "DIC_OPTIONS.h"
29ad036528 Step*0002 #include "PTRACERS_OPTIONS.h"
                0003
08536d17ba Step*0004 CBOP
                0005 C !ROUTINE: DIC_SURFFORCING_INIT
                0006
                0007 C !INTERFACE: ==========================================================
29ad036528 Step*0008       SUBROUTINE DIC_SURFFORCING_INIT(
a1d0e455fd Hann*0009      I          myThid )
29ad036528 Step*0010
08536d17ba Step*0011 C !DESCRIPTION:
e28bbbf906 Jean*0012 C  Calculate first guess of pH
29ad036528 Step*0013
08536d17ba Step*0014 C !USES: ===============================================================
                0015       IMPLICIT NONE
29ad036528 Step*0016 #include "SIZE.h"
                0017 #include "DYNVARS.h"
                0018 #include "EEPARAMS.h"
                0019 #include "PARAMS.h"
                0020 #include "GRID.h"
                0021 #include "FFIELDS.h"
2ef8966791 Davi*0022 #include "DIC_VARS.h"
bcc34e2df6 Jean*0023 #include "PTRACERS_SIZE.h"
d800a455f8 Jean*0024 #include "PTRACERS_PARAMS.h"
e28bbbf906 Jean*0025 #include "PTRACERS_FIELDS.h"
29ad036528 Step*0026
08536d17ba Step*0027 C !INPUT PARAMETERS: ===================================================
                0028 C  myThid               :: thread number
29ad036528 Step*0029       INTEGER  myThid
                0030
946d3aa393 Jean*0031 #ifdef ALLOW_DIC
08536d17ba Step*0032
                0033 C !LOCAL VARIABLES: ====================================================
51e381e9c9 Jean*0034       INTEGER i,j, kLev, it
29ad036528 Step*0035       INTEGER iMin,iMax,jMin,jMax, bi, bj
51e381e9c9 Jean*0036        _RL co3dummy
                0037 C local variables for carbon chem
29ad036528 Step*0038       _RL surfalk(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0039       _RL surfphos(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0040       _RL surfsi(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8bf2c0e0ad Step*0041       _RL surftemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0042       _RL surfsalt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0043       _RL surfdic(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0719e4347a Jean*0044       INTEGER iprt,jprt
6acab690ae Jona*0045       LOGICAL debugPrt
51e381e9c9 Jean*0046 #ifdef ALLOW_DEBUG
6acab690ae Jona*0047 C     msgBuf     :: Informational/error message buffer
                0048       CHARACTER*(MAX_LEN_MBUF) msgBuf
51e381e9c9 Jean*0049 #endif
08536d17ba Step*0050 CEOP
29ad036528 Step*0051
                0052 C =================================================================
7f407c2fb7 Davi*0053
51e381e9c9 Jean*0054       kLev = 1
7f407c2fb7 Davi*0055       jMin=1
                0056       jMax=sNy
                0057       iMin=1
                0058       iMax=sNx
                0059
51e381e9c9 Jean*0060 C Solubility relation coefficients
3a677ccb28 Davi*0061       DO bj=myByLo(myThid),myByHi(myThid)
                0062        DO bi=myBxLo(myThid),myBxHi(myThid)
                0063         DO j=1-OLy,sNy+OLy
                0064          DO i=1-OLx,sNx+OLx
                0065           ak0(i,j,bi,bj)=0. _d 0
                0066           ak1(i,j,bi,bj)=0. _d 0
                0067           ak2(i,j,bi,bj)=0. _d 0
                0068           akw(i,j,bi,bj)=0. _d 0
                0069           akb(i,j,bi,bj)=0. _d 0
                0070           akf(i,j,bi,bj)=0. _d 0
                0071           ak1p(i,j,bi,bj)=0. _d 0
                0072           ak2p(i,j,bi,bj)=0. _d 0
                0073           ak3p(i,j,bi,bj)=0. _d 0
                0074           aksi(i,j,bi,bj)=0. _d 0
d0092a57ac Step*0075           fugf(i,j,bi,bj)=0. _d 0
3a677ccb28 Davi*0076           ff(i,j,bi,bj)=0. _d 0
                0077           ft(i,j,bi,bj)=0. _d 0
                0078           st(i,j,bi,bj)=0. _d 0
                0079           bt(i,j,bi,bj)=0. _d 0
                0080          ENDDO
                0081         ENDDO
                0082        ENDDO
                0083       ENDDO
                0084
d800a455f8 Jean*0085       DO bj=myByLo(myThid),myByHi(myThid)
                0086        DO bi=myBxLo(myThid),myBxHi(myThid)
29ad036528 Step*0087
                0088 C determine inorganic carbon chem coefficients
d800a455f8 Jean*0089         DO j=jMin,jMax
                0090          DO i=iMin,iMax
29ad036528 Step*0091 #ifdef DIC_BIOTIC
8bf2c0e0ad Step*0092 #ifdef DIC_BOUNDS
a1d0e455fd Hann*0093             surfalk(i,j)  = MAX( 0.4 _d 0,
                0094      &                      MIN( 10. _d 0, PTRACER(i,j,kLev,bi,bj,2) ) )
                0095      &                    * maskC(i,j,kLev,bi,bj)
                0096             surfphos(i,j) = MAX( 1.0 _d -11,
                0097      &                      MIN( 1. _d -1, PTRACER(i,j,kLev,bi,bj,3) ) )
                0098      &                    * maskC(i,j,kLev,bi,bj)
8bf2c0e0ad Step*0099 #else
a1d0e455fd Hann*0100             surfalk(i,j)  = PTRACER(i,j,kLev,bi,bj,2)
                0101      &                    * maskC(i,j,kLev,bi,bj)
                0102             surfphos(i,j) = PTRACER(i,j,kLev,bi,bj,3)
                0103      &                    * maskC(i,j,kLev,bi,bj)
8bf2c0e0ad Step*0104 #endif
6acab690ae Jona*0105 #else /* DIC_BIOTIC */
a1d0e455fd Hann*0106             surfalk(i,j)  = 2.366595 _d 0 *salt(i,j,kLev,bi,bj)/35. _d 0
                0107      &                    * maskC(i,j,kLev,bi,bj)
                0108             surfphos(i,j) = 5.1225 _d -4 * maskC(i,j,kLev,bi,bj)
6acab690ae Jona*0109 #endif /* DIC_BIOTIC */
a1d0e455fd Hann*0110 C for non-interactive Si
                0111             surfsi(i,j)   = silicaSurf(i,j,bi,bj)*maskC(i,j,kLev,bi,bj)
8bf2c0e0ad Step*0112 #ifdef DIC_BOUNDS
a1d0e455fd Hann*0113             surftemp(i,j) = MAX( -4. _d 0,
                0114      &                      MIN( 50. _d 0, theta(i,j,kLev,bi,bj) ) )
                0115             surfsalt(i,j) = MAX( 4. _d 0,
                0116      &                      MIN( 50. _d 0, salt(i,j,kLev,bi,bj) ) )
                0117             surfdic(i,j)  = MAX( 0.4 _d 0,
                0118      &                      MIN( 10. _d 0, PTRACER(i,j,kLev,bi,bj,1) ) )
8bf2c0e0ad Step*0119 #else
                0120             surftemp(i,j) = theta(i,j,kLev,bi,bj)
                0121             surfsalt(i,j) = salt(i,j,kLev,bi,bj)
dca99e3ae0 Jean*0122             surfdic(i,j)  = PTRACER(i,j,kLev,bi,bj,1)
                0123      &                    * maskC(i,j,kLev,bi,bj)
8bf2c0e0ad Step*0124 #endif
29ad036528 Step*0125          ENDDO
75e97f1e14 Davi*0126         ENDDO
29ad036528 Step*0127
6acab690ae Jona*0128 #ifdef CARBONCHEM_SOLVESAPHE
                0129 #ifdef ALLOW_DEBUG
                0130         IF (debugMode) CALL DEBUG_CALL('DIC_COEFFS_SURF',myThid)
                0131 #endif
                0132         CALL DIC_COEFFS_SURF(
                0133      I                       surftemp,surfsalt,
                0134      I                       bi,bj,iMin,iMax,jMin,jMax,myThid)
                0135 #else /* CARBONCHEM_SOLVESAPHE */
                0136 #ifdef ALLOW_DEBUG
                0137         IF (debugMode) CALL DEBUG_CALL('CARBON_COEFFS',myThid)
                0138 #endif
75e97f1e14 Davi*0139         CALL CARBON_COEFFS(
8bf2c0e0ad Step*0140      I                       surftemp,surfsalt,
5625485478 Jean*0141      I                       bi,bj,iMin,iMax,jMin,jMax,myThid)
a1d0e455fd Hann*0142
6acab690ae Jona*0143 #endif /* CARBONCHEM_SOLVESAPHE */
7f407c2fb7 Davi*0144
29ad036528 Step*0145 C====================================================================
                0146
51e381e9c9 Jean*0147         IF ( .NOT.pH_isLoaded(1) ) THEN
29ad036528 Step*0148 C set guess of pH for first step here
6acab690ae Jona*0149 #ifdef ALLOW_DEBUG
                0150           IF (debugMode) THEN
                0151             WRITE(msgBuf,'(A)') 'Initial pCO2 approximation method'
                0152             CALL DEBUG_MSG(msgBuf(1:33),myThid)
                0153           ENDIF
                0154 #endif
29ad036528 Step*0155
6acab690ae Jona*0156           debugPrt = debugMode
a1d0e455fd Hann*0157 C first approximation
29ad036528 Step*0158 C$TAF LOOP = parallel
d800a455f8 Jean*0159           DO j=jMin,jMax
29ad036528 Step*0160 C$TAF LOOP = parallel
d800a455f8 Jean*0161            DO i=iMin,iMax
0719e4347a Jean*0162             IF ( maskC(i,j,kLev,bi,bj) .NE. 0. _d 0) THEN
29ad036528 Step*0163 C$TAF init dic_surf = static, 10
6acab690ae Jona*0164 #ifdef CARBONCHEM_SOLVESAPHE
                0165              IF ( selectPHsolver.GT.0 ) THEN
                0166 C Use Munhoven (2013) Solvesaphe routine to initialize pH
                0167 #ifdef ALLOW_DEBUG
                0168               IF (debugPrt) CALL DEBUG_CALL('AHINI_FOR_AT',myThid)
                0169 #endif
                0170 C call AHINI_FOR_AT to get better initial guess of pH
                0171               CALL AHINI_FOR_AT(
                0172      I           surfalk(i,j)*permil,
                0173      I           surfdic(i,j)*permil,
                0174      I           bt(i,j,bi,bj),
                0175      U           pH(i,j,bi,bj),
                0176      I           i,j,kLev,bi,bj,nIter0,myThid )
                0177
dca99e3ae0 Jean*0178 C$TAF STORE pH(i,j,bi,bj)                              = dic_surf
29ad036528 Step*0179 C$TAF STORE surfalk(i,j), surfphos(i,j), surfsi(i,j)   = dic_surf
6acab690ae Jona*0180
                0181 #ifdef ALLOW_DEBUG
                0182               IF (debugPrt)
                0183      &          CALL DEBUG_CALL('CALC_PCO2_SOLVESAPHE',myThid)
                0184 #endif
                0185               CALL CALC_PCO2_SOLVESAPHE(
                0186      I         surftemp(i,j),surfsalt(i,j),
                0187      I         surfdic(i,j), surfphos(i,j),
                0188      I         surfsi(i,j),surfalk(i,j),
                0189      U         pH(i,j,bi,bj),pCO2(i,j,bi,bj),co3dummy,
                0190      I         i,j,kLev,bi,bj, debugPrt, nIter0, myThid )
                0191               debugPrt = .FALSE.
                0192              ELSE
                0193 C Use the original Follows et al. (2006) solver
                0194 #endif /* CARBONCHEM_SOLVESAPHE */
                0195 #ifdef ALLOW_DEBUG
                0196               IF (debugPrt) THEN
                0197                 CALL DEBUG_CALL('CALC_PCO2_APPROX',myThid)
                0198                 debugPrt = .FALSE.
                0199               ENDIF
                0200 #endif
                0201               DO it=1,10
7448700841 Mart*0202 cC$TAF STORE pH(i,j,bi,bj)                              = dic_surf
                0203 cC$TAF STORE surfalk(i,j), surfphos(i,j), surfsi(i,j)   = dic_surf
6acab690ae Jona*0204                CALL CALC_PCO2_APPROX(
                0205      I          surftemp(i,j),surfsalt(i,j),
                0206      I          surfdic(i,j), surfphos(i,j),
                0207      I          surfsi(i,j),surfalk(i,j),
                0208      I          ak1(i,j,bi,bj),ak2(i,j,bi,bj),
                0209      I          ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
                0210      I          aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
                0211      I          aksi(i,j,bi,bj),akf(i,j,bi,bj),
                0212      I          ak0(i,j,bi,bj), fugf(i,j,bi,bj),
                0213      I          ff(i,j,bi,bj),
                0214      I          bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
                0215      U          pH(i,j,bi,bj),pCO2(i,j,bi,bj),co3dummy,
                0216      I          i,j,kLev,bi,bj, it ,  myThid )
                0217               ENDDO
                0218 #ifdef CARBONCHEM_SOLVESAPHE
                0219              ENDIF
                0220 #endif /* CARBONCHEM_SOLVESAPHE */
0719e4347a Jean*0221             ENDIF
7f407c2fb7 Davi*0222            ENDDO
                0223           ENDDO
6acab690ae Jona*0224
                0225 #ifdef ALLOW_DEBUG
a1d0e455fd Hann*0226          IF (debugMode) THEN
0719e4347a Jean*0227           iprt = MIN(20,sNx)
                0228           jprt = MIN(20,sNy)
6acab690ae Jona*0229           WRITE(msgBuf,'(4(A,F9.6),2(A,F11.8),A,F9.6)')
                0230      &        ' first guess pH=', pH(iprt,jprt,bi,bj),
                0231      &        ', Temp=',theta(iprt,jprt,1,bi,bj),
                0232      &        ', Salt=',salt(iprt,jprt,1,bi,bj),
                0233      &        ', DIC=', surfdic(iprt,jprt),
                0234      &        ', PO4=', surfphos(iprt,jprt),
                0235      &        ', SiT=', surfsi(iprt,jprt),
                0236      &        ', ALK=', surfalk(iprt,jprt)
                0237           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0238      &                        SQUEEZE_RIGHT, myThid )
                0239 c         WRITE(standardMessageUnit,*) 'QQ first guess pH ',
                0240 c    &        pH(iprt,jprt,bi,bj),' Temp ',theta(iprt,jprt,1,bi,bj),
                0241 c    &        ' Salt ',salt(iprt,jprt,1,bi,bj),
                0242 c    &        ' DIC ' ,surfdic(iprt,jprt),
                0243 c    &        ' PO4 ' ,surfphos(iprt,jprt),
                0244 c    &        ' SiT ' ,surfsi(iprt,jprt),
                0245 c    &        ' ALK ' ,surfalk(iprt,jprt)
a1d0e455fd Hann*0246 c         CALL DEBUG_MSG(msgBuf,myThid)
                0247          ENDIF
6acab690ae Jona*0248 #endif
29ad036528 Step*0249
a1d0e455fd Hann*0250 C end if-block (.NOT.pH_isLoaded)
6acab690ae Jona*0251         ENDIF
0719e4347a Jean*0252 C     end bi,bj loops
                0253        ENDDO
                0254       ENDDO
7f407c2fb7 Davi*0255
946d3aa393 Jean*0256 #endif /* ALLOW_DIC */
0719e4347a Jean*0257       RETURN
                0258       END