pkg/darwin/darwin_surfforcing_init.F

8fbfd1f382 Oliv*0001 #include "DARWIN_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: DARWIN_SURFFORCING_INIT
                0005
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE DARWIN_SURFFORCING_INIT(myThid)
                0008
                0009 C !DESCRIPTION:
                0010 C  Calculate first guess of pH
                0011
                0012 C !USES: ===============================================================
                0013       IMPLICIT NONE
                0014 #include "SIZE.h"
                0015 #include "DYNVARS.h"
                0016 #include "EEPARAMS.h"
                0017 #include "PARAMS.h"
                0018 #include "GRID.h"
                0019 #include "FFIELDS.h"
                0020 #include "PTRACERS_SIZE.h"
                0021 #include "PTRACERS_PARAMS.h"
                0022 #include "PTRACERS_FIELDS.h"
                0023 #ifdef ALLOW_DARWIN
                0024 #include "DARWIN_SIZE.h"
                0025 #include "DARWIN_INDICES.h"
                0026 #include "DARWIN_EXF_FIELDS.h"
                0027 #include "DARWIN_PARAMS.h"
                0028 #include "DARWIN_TRAITS.h"
                0029 #include "DARWIN_FIELDS.h"
                0030 #endif
                0031
                0032 C !INPUT PARAMETERS: ===================================================
                0033 C  myThid               :: thread number
                0034       INTEGER  myThid
                0035
                0036 #ifdef ALLOW_DARWIN
                0037 #ifdef DARWIN_ALLOW_CARBON
                0038
                0039 C !LOCAL VARIABLES: ====================================================
f0a72e2151 Oliv*0040        LOGICAL debugPrt
8fbfd1f382 Oliv*0041        INTEGER i,j,k,it
895c6145db Oliv*0042        _RL pco2dummy
8fbfd1f382 Oliv*0043        _RL co3dummy
                0044 C Number of iterations for pCO2 solvers...
                0045 C Solubility relation coefficients
                0046 C local variables for carbon chem
                0047       INTEGER iMin,iMax,jMin,jMax, bi, bj
                0048       _RL surfdic
                0049       _RL surfalk
                0050       _RL surfphos
                0051       _RL surfsi
                0052       _RL surftemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0053       _RL surfsalt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0054 CEOP
                0055
                0056 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
                0057
                0058 #ifdef ALLOW_OFFLINE
                0059       IF ( useOffLine ) THEN
                0060         CALL OFFLINE_FIELDS_LOAD( startTime, nIter0, myThid )
                0061       ENDIF
                0062 #endif
                0063
                0064 C =================================================================
                0065
                0066       jMin=1
                0067       jMax=sNy
                0068       iMin=1
                0069       iMax=sNx
                0070
f17484eeb1 Oliv*0071       IF (nIter0.EQ.PTRACERS_Iter0 .AND. darwin_pickupSuff.EQ.' ') THEN
8fbfd1f382 Oliv*0072 C set guess of pH for first step here
                0073
                0074        DO bj=myByLo(myThid),myByHi(myThid)
                0075         DO bi=myBxLo(myThid),myBxHi(myThid)
                0076
                0077          DO k=1,Nr
                0078
                0079 C determine inorganic carbon chem coefficients
                0080           DO j=jMin,jMax
                0081            DO i=iMin,iMax
                0082             surfsalt(i,j) = MAX(surfSaltMinInit,
                0083      &                      MIN(surfSaltMaxInit, salt(i,j,k,bi,bj)))
                0084             surftemp(i,j) = MAX(surfTempMinInit,
                0085      &                      MIN(surfTempMaxInit, theta(i,j,k,bi,bj)))
                0086            ENDDO
                0087           ENDDO
f0a72e2151 Oliv*0088
                0089           DO j=1-OLy,sNy+OLy
                0090            DO i=1-OLx,sNx+OLx
                0091             ak0(i,j,bi,bj)=0. _d 0
                0092             ak1(i,j,bi,bj)=0. _d 0
                0093             ak2(i,j,bi,bj)=0. _d 0
                0094             akw(i,j,bi,bj)=0. _d 0
                0095             akb(i,j,bi,bj)=0. _d 0
                0096             akf(i,j,bi,bj)=0. _d 0
                0097             ak1p(i,j,bi,bj)=0. _d 0
                0098             ak2p(i,j,bi,bj)=0. _d 0
                0099             ak3p(i,j,bi,bj)=0. _d 0
                0100             aksi(i,j,bi,bj)=0. _d 0
                0101             fugf(i,j,bi,bj)=0. _d 0
                0102             ff(i,j,bi,bj)=0. _d 0
                0103             ft(i,j,bi,bj)=0. _d 0
                0104             st(i,j,bi,bj)=0. _d 0
                0105             bt(i,j,bi,bj)=0. _d 0
                0106            ENDDO
                0107           ENDDO
                0108
                0109 #ifdef DARWIN_SOLVESAPHE
                0110 #ifdef ALLOW_DEBUG
                0111           IF (debugMode) CALL DEBUG_CALL('DARWIN_COEFFS_SURF',myThid)
                0112 #endif
                0113           CALL DARWIN_COEFFS_SURF(
                0114      I                       surftemp,surfsalt,
45ab1c3374 Oliv*0115      I                       bi,bj,iMin,iMax,jMin,jMax,
                0116      I                       k,myThid)
f0a72e2151 Oliv*0117
                0118 C Now correct the coefficients for pressure dependence
                0119           CALL DARWIN_COEFFS_DEEP(
                0120      I                       surftemp, surfsalt,
                0121      I                       bi,bj,iMin,iMax,jMin,jMax,
                0122      I                       k,myThid)
                0123 #else /* DARWIN_SOLVESAPHE */
                0124 #ifdef ALLOW_DEBUG
                0125           IF (debugMode) CALL DEBUG_CALL('CARBON_COEFFS',myThid)
                0126 #endif
8fbfd1f382 Oliv*0127           CALL DARWIN_CARBON_COEFFS(
                0128      I                       surftemp,surfsalt,
                0129      I                       bi,bj,iMin,iMax,jMin,jMax,k,myThid)
f0a72e2151 Oliv*0130 #endif /* DARWIN_SOLVESAPHE */
8fbfd1f382 Oliv*0131
                0132           DO j=1-OLy,sNy+OLy
                0133            DO i=1-OLx,sNx+OLx
                0134             pH(i,j,k,bi,bj) = 8. _d 0
                0135            ENDDO
                0136           ENDDO
                0137
f0a72e2151 Oliv*0138           debugPrt = debugMode
                0139
8fbfd1f382 Oliv*0140 c first approximation
                0141 C$TAF LOOP = parallel
                0142           DO j=jMin,jMax
                0143 C$TAF LOOP = parallel
                0144            DO i=iMin,iMax
                0145             IF ( maskC(i,j,k,bi,bj) .NE. 0. _d 0) THEN
                0146 c put bounds on tracers so pH solver doesn't blow up
                0147              surfdic  = ptr2mol * MAX(surfDICMinInit,
                0148      &                  MIN(surfDICMaxInit, Ptracer(i,j,k,bi,bj,iDIC)))
                0149      &                  * maskC(i,j,k,bi,bj)
                0150              surfalk  = ptr2mol * MAX(surfALKMinInit,
                0151      &                  MIN(surfALKMaxInit, Ptracer(i,j,k,bi,bj,iALK)))
                0152      &                  * maskC(i,j,k,bi,bj)
                0153              surfphos = ptr2mol * MAX(surfPO4MinInit,
                0154      &                  MIN(surfPO4MaxInit, Ptracer(i,j,k,bi,bj,iPO4)))
                0155      &                  * maskC(i,j,k,bi,bj)
20ff2335d4 Oliv*0156              surfsi   = ptr2mol * MAX(surfSiMinInit,
8fbfd1f382 Oliv*0157      &                  MIN(surfSiMaxInit, Ptracer(i,j,k,bi,bj,iSiO2)))
                0158      &                  * maskC(i,j,k,bi,bj)
                0159 C$TAF init dic_surf = static, 10
f0a72e2151 Oliv*0160 #ifdef DARWIN_SOLVESAPHE
                0161              IF ( selectPHsolver.GT.0 ) THEN
                0162 C Use Munhoven (2013) Solvesaphe routine to initialize pH
                0163 #ifdef ALLOW_DEBUG
                0164               IF (debugPrt) CALL DEBUG_CALL('AHINI_FOR_AT',myThid)
                0165 #endif
                0166 C call AHINI_FOR_AT to get better initial guess of pH
                0167               CALL AHINI_FOR_AT(
9270923c0b Oliv*0168      I           surfalk*m3perkg,
                0169      I           surfdic*m3perkg,
f0a72e2151 Oliv*0170      I           bt(i,j,bi,bj),
                0171      U           pH(i,j,k,bi,bj),
                0172      I           i,j,k,bi,bj,nIter0,myThid )
                0173
                0174 C$TAF STORE pH(i,j,k,bi,bj)                            = dic_surf
                0175 C$TAF STORE surfalk, surfphos, surfsi                  = dic_surf
                0176
                0177 #ifdef ALLOW_DEBUG
                0178               IF (debugPrt)
                0179      &          CALL DEBUG_CALL('CALC_PCO2_SOLVESAPHE',myThid)
                0180 #endif
                0181               CALL CALC_PCO2_SOLVESAPHE(
                0182      I         surftemp(i,j),surfsalt(i,j),
                0183      I         surfdic, surfphos,
                0184      I         surfsi,surfalk,
895c6145db Oliv*0185      U         pH(i,j,k,bi,bj),
                0186      O         pco2dummy,co3dummy,
f0a72e2151 Oliv*0187      I         i,j,k,bi,bj, debugPrt, nIter0, myThid )
                0188               debugPrt = .FALSE.
                0189              ELSE
                0190 C Use the original Follows et al. (2006) solver
                0191 #endif /* DARWIN_SOLVESAPHE */
                0192 #ifdef ALLOW_DEBUG
                0193               IF (debugPrt) THEN
                0194                 CALL DEBUG_CALL('CALC_PCO2_APPROX',myThid)
                0195                 debugPrt = .FALSE.
                0196               ENDIF
                0197 #endif
8fbfd1f382 Oliv*0198              DO it=1,10
f0a72e2151 Oliv*0199 C$TAF STORE pH(i,j,k,bi,bj)                            = dic_surf
                0200 C$TAF STORE surfalk, surfphos, surfs                   = dic_surf
8fbfd1f382 Oliv*0201               CALL DARWIN_CALC_PCO2_APPROX(
                0202      I        surftemp(i,j), surfsalt(i,j),
                0203      I        surfdic, surfphos, surfsi, surfalk,
                0204      I        ak1(i,j,bi,bj), ak2(i,j,bi,bj),
                0205      I        ak1p(i,j,bi,bj), ak2p(i,j,bi,bj), ak3p(i,j,bi,bj),
                0206      I        aks(i,j,bi,bj), akb(i,j,bi,bj), akw(i,j,bi,bj),
                0207      I        aksi(i,j,bi,bj), akf(i,j,bi,bj),
                0208      I        ak0(i,j,bi,bj), fugf(i,j,bi,bj),
                0209      I        ff(i,j,bi,bj),
                0210      I        bt(i,j,bi,bj), st(i,j,bi,bj), ft(i,j,bi,bj),
895c6145db Oliv*0211      U        pH(i,j,k,bi,bj),
                0212      O        pco2dummy, co3dummy,
8fbfd1f382 Oliv*0213      I        i,j,k,bi,bj, it ,  myThid )
                0214              ENDDO
f0a72e2151 Oliv*0215 #ifdef DARWIN_SOLVESAPHE
                0216              ENDIF
                0217 #endif /* DARWIN_SOLVESAPHE */
8fbfd1f382 Oliv*0218             ENDIF
                0219            ENDDO
                0220           ENDDO
                0221
                0222 C        end k loop
                0223          ENDDO
                0224
                0225 C      end bi,bj loops
                0226         ENDDO
                0227        ENDDO
                0228
                0229 C     time to initialize
                0230       ENDIF
                0231
                0232 #endif /* DARWIN_ALLOW_CARBON */
                0233 #endif /* ALLOW_DARWIN */
                0234
                0235       RETURN
                0236       END