Back to home page

darwin3
 
 

    

File indexing completed on 2024-12-17 18:33:12 UTC

view on githubraw file Latest commit a2844551 on 2020-07-25 02:56:46 UTC
c0d1c06c15 Matt*0001 #include "BLING_OPTIONS.h"
                0002 #include "PTRACERS_OPTIONS.h"
                0003 
                0004 CBOP
                0005       subroutine BLING_CARBONATE_INIT( myThid )
                0006 
                0007 C     ==========================================================
4ac06494d5 Matt*0008 C     | subroutine bling_carbonate_init
c0d1c06c15 Matt*0009 C     | o Calculate first guess of pH
4ac06494d5 Matt*0010 C     |   Adapted from pkg/dic/dic_surfforcing_init.F
c0d1c06c15 Matt*0011 C     ==========================================================
                0012 
e0f9a7ba0b Matt*0013       IMPLICIT NONE
                0014 
c0d1c06c15 Matt*0015 C     === Global variables ===
                0016 #include "SIZE.h"
                0017 #include "DYNVARS.h"
                0018 #include "EEPARAMS.h"
                0019 #include "PARAMS.h"
                0020 #include "GRID.h"
                0021 #include "FFIELDS.h"
                0022 #include "BLING_VARS.h"
                0023 #include "PTRACERS_SIZE.h"
                0024 #include "PTRACERS_PARAMS.h"
                0025 #include "PTRACERS_FIELDS.h"
                0026 #include "BLING_LOAD.h"
                0027 
                0028 C     === Routine arguments ===
                0029 C     myThid               :: thread Id. number
                0030       INTEGER  myThid
                0031 
                0032 #ifdef ALLOW_BLING
                0033 C     === Local variables ===
                0034        INTEGER i,j, k, it
                0035        INTEGER intimeP, intime0, intime1
                0036        _RL aWght, bWght, co3dummy
                0037 C Number of iterations for pCO2 solvers...
                0038 C Solubility relation coefficients
                0039 C local variables for carbon chem
                0040       INTEGER iMin,iMax,jMin,jMax, bi, bj
                0041       _RL alktmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0042       _RL phostmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0043       _RL sitmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0044       _RL thetatmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0045       _RL salttmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0046       _RL dictmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0047       LOGICAL pH_isLoaded
                0048 CEOP
                0049 
e0f9a7ba0b Matt*0050 #ifndef USE_SIBLING
c0d1c06c15 Matt*0051       IF ( periodicExternalForcing ) THEN
                0052 
                0053 c read in silica field
                0054          CALL LEF_ZERO( silica0,myThid )
                0055          CALL LEF_ZERO( silica1,myThid )
                0056 
                0057 C--   Now calculate whether it is time to update the forcing arrays
                0058        CALL GET_PERIODIC_INTERVAL(
                0059      O                   intimeP, intime0, intime1, bWght, aWght,
                0060      I                   externForcingCycle, externForcingPeriod,
e0f9a7ba0b Matt*0061      I                   deltaTClock, startTime, myThid )
c0d1c06c15 Matt*0062 
                0063        _BARRIER
                0064        _BEGIN_MASTER(myThid)
                0065         WRITE(standardMessageUnit,'(A,I10,A,2(2I5,A))')
                0066      &   ' BLING_SURFFORCING_INIT, it=', nIter0,
                0067      &   ' : Reading new data, i0,i1=', intime0, intime1
                0068        _END_MASTER(myThid)
                0069 
                0070        IF ( BLING_silicaFile .NE. ' '  ) THEN
                0071          CALL READ_REC_XY_RS( BLING_silicaFile,silica0,intime0,
                0072      &        nIter0,myThid )
                0073          CALL READ_REC_XY_RS( BLING_silicaFile,silica1,intime1,
                0074      &        nIter0,myThid )
                0075        ENDIF
                0076 
                0077 #ifdef ALLOW_OFFLINE
                0078        IF ( useOffLine ) THEN
                0079          CALL OFFLINE_FIELDS_LOAD( startTime, nIter0, myThid )
                0080        ENDIF
                0081 #endif
                0082 
                0083        _EXCH_XY_RS(silica0, myThid )
                0084        _EXCH_XY_RS(silica1, myThid )
                0085 
                0086        IF ( BLING_silicaFile .NE. ' '  ) THEN
                0087         DO bj = myByLo(myThid), myByHi(myThid)
                0088          DO bi = myBxLo(myThid), myBxHi(myThid)
e0f9a7ba0b Matt*0089           DO j=1-OLy,sNy+OLy
                0090            DO i=1-OLx,sNx+OLx
                0091              silica(i,j,bi,bj)= bWght*silica0(i,j,bi,bj)
c0d1c06c15 Matt*0092      &                        + aWght*silica1(i,j,bi,bj)
                0093            ENDDO
                0094           ENDDO
                0095          ENDDO
                0096         ENDDO
                0097        ENDIF
                0098 
                0099 c end periodicExternalForcing
                0100       ENDIF
e0f9a7ba0b Matt*0101 c end SiBLING
                0102 #endif
c0d1c06c15 Matt*0103 
                0104 C =================================================================
                0105 
                0106       jMin=1
                0107       jMax=sNy
                0108       iMin=1
                0109       iMax=sNx
                0110 
                0111       DO k=1,Nr
                0112       DO bj=myByLo(myThid),myByHi(myThid)
                0113        DO bi=myBxLo(myThid),myBxHi(myThid)
                0114          DO j=1-OLy,sNy+OLy
e0f9a7ba0b Matt*0115           DO i=1-OLx,sNx+OLx
c0d1c06c15 Matt*0116             pH(i,j,k,bi,bj) = 8. _d 0
                0117           ENDDO
                0118          ENDDO
                0119        ENDDO
                0120       ENDDO
                0121       ENDDO
                0122 
                0123       DO bj=myByLo(myThid),myByHi(myThid)
                0124        DO bi=myBxLo(myThid),myBxHi(myThid)
                0125         DO j=1-OLy,sNy+OLy
                0126          DO i=1-OLx,sNx+OLx
                0127           ak0(i,j,bi,bj)=0. _d 0
                0128           ak1(i,j,bi,bj)=0. _d 0
                0129           ak2(i,j,bi,bj)=0. _d 0
                0130           akw(i,j,bi,bj)=0. _d 0
                0131           akb(i,j,bi,bj)=0. _d 0
                0132           akf(i,j,bi,bj)=0. _d 0
                0133           ak1p(i,j,bi,bj)=0. _d 0
                0134           ak2p(i,j,bi,bj)=0. _d 0
                0135           ak3p(i,j,bi,bj)=0. _d 0
                0136           aksi(i,j,bi,bj)=0. _d 0
                0137           fugf(i,j,bi,bj)=0. _d 0
                0138           ff(i,j,bi,bj)=0. _d 0
                0139           ft(i,j,bi,bj)=0. _d 0
                0140           st(i,j,bi,bj)=0. _d 0
                0141           bt(i,j,bi,bj)=0. _d 0
a284455135 Matt*0142 #ifdef CARBONCHEM_SOLVESAPHE
                0143           cat(i,j,bi,bj)=0. _d 0
                0144           akn(i,j,bi,bj)=0. _d 0
                0145           akhs(i,j,bi,bj)=0. _d 0
                0146           aphscale(i,j,bi,bj)=0. _d 0
                0147 #endif
c0d1c06c15 Matt*0148          ENDDO
                0149         ENDDO
                0150        ENDDO
                0151       ENDDO
                0152 
                0153       pH_isLoaded = .FALSE.
                0154       IF ( nIter0.GT.PTRACERS_Iter0 .OR.
                0155      &    (nIter0.EQ.PTRACERS_Iter0 .AND. pickupSuff.NE.' ')
                0156      &   ) THEN
                0157 C       Read pH from a pickup file if needed
                0158         CALL BLING_READ_PICKUP(
                0159      O                        pH_isLoaded,
                0160      I                        nIter0, myThid )
                0161       ENDIF
                0162 
                0163       DO bj=myByLo(myThid),myByHi(myThid)
                0164        DO bi=myBxLo(myThid),myBxHi(myThid)
                0165 
                0166 C determine inorganic carbon chem coefficients
                0167 
                0168         IF ( .NOT.pH_isLoaded ) THEN
                0169 C set guess of pH for first step here
                0170 
                0171         DO k=1,Nr
e0f9a7ba0b Matt*0172          DO j=jMin,jMax
                0173           DO i=iMin,iMax
                0174              thetatmp(i,j) = theta(i,j,k,bi,bj)
                0175              salttmp(i,j) = salt(i,j,k,bi,bj)
                0176              dictmp(i,j) = PTRACER(i,j,k,bi,bj,1)
                0177      &                          * maskC(i,j,k,bi,bj)
c0d1c06c15 Matt*0178              alktmp(i,j) = PTRACER(i,j,k,bi,bj,2)
                0179      &                          * maskC(i,j,k,bi,bj)
a284455135 Matt*0180 #ifndef USE_BLING_V1
e0f9a7ba0b Matt*0181              phostmp(i,j) = PTRACER(i,j,k,bi,bj,5)
c0d1c06c15 Matt*0182      &                          * maskC(i,j,k,bi,bj)
a284455135 Matt*0183 #else
                0184              phostmp(i,j) = PTRACER(i,j,k,bi,bj,4)
                0185      &                          * maskC(i,j,k,bi,bj)
                0186 #endif
e0f9a7ba0b Matt*0187 #ifdef USE_SIBLING
                0188              sitmp(i,j) = PTRACER(i,j,k,bi,bj,9)
                0189      &                          * maskC(i,j,k,bi,bj)
                0190 #else
c0d1c06c15 Matt*0191 C FOR NON-INTERACTIVE Si
                0192              IF ( k.eq.1 ) THEN
e0f9a7ba0b Matt*0193               sitmp(i,j) = silica(i,j,bi,bj) * maskC(i,j,k,bi,bj)
c0d1c06c15 Matt*0194              ELSE
e0f9a7ba0b Matt*0195               sitmp(i,j) = 0.03 * maskC(i,j,k,bi,bj)
c0d1c06c15 Matt*0196              ENDIF
e0f9a7ba0b Matt*0197 #endif
                0198           ENDDO
c0d1c06c15 Matt*0199          ENDDO
                0200 
e0f9a7ba0b Matt*0201 #ifdef CARBONCHEM_SOLVESAPHE
                0202 #ifdef ALLOW_DEBUG
                0203       IF (debugMode) CALL DEBUG_CALL('DIC_COEFFS_DEEP',myThid)
                0204 #endif
                0205         CALL DIC_COEFFS_SURF(
                0206 cav changing surftheta to thetatmp
                0207      I                       thetatmp,salttmp,
                0208      I                       bi,bj,iMin,iMax,jMin,jMax,myThid)
                0209 C Now correct the coefficients for pressure dependence
                0210         CALL DIC_COEFFS_DEEP(
                0211      I                       thetatmp,salttmp,
                0212      I                       bi,bj,iMin,iMax,jMin,jMax,
                0213      I                       k,myThid)
                0214 #else
                0215 #ifdef ALLOW_DEBUG
                0216       IF (debugMode) CALL DEBUG_CALL('CARBON_COEFFS',myThid)
                0217 #endif
                0218          CALL CARBON_COEFFS_PRESSURE_DEP(
c0d1c06c15 Matt*0219      I                       thetatmp,salttmp,
                0220      I                       bi,bj,iMin,iMax,jMin,jMax,k,myThid)
e0f9a7ba0b Matt*0221 #endif
c0d1c06c15 Matt*0222 
                0223 C====================================================================
                0224 
                0225 c first approximation
                0226 
                0227           DO j=jMin,jMax
                0228            DO i=iMin,iMax
                0229             IF ( maskC(i,j,k,bi,bj) .NE. 0. _d 0) THEN
e0f9a7ba0b Matt*0230 
                0231 #ifdef CARBONCHEM_SOLVESAPHE
                0232              IF ( selectPHsolver.GT.0 ) THEN
                0233 C Use Munhoven (2013) Solvesaphe routine to initialize pH
                0234 #ifdef ALLOW_DEBUG
                0235               IF (debugMode) CALL DEBUG_CALL(
                0236      &            'AHINI_FOR_AT',myThid)
                0237 #endif
                0238 C call AHINI_FOR_AT to get better initial guess of pH
                0239               CALL AHINI_FOR_AT(
                0240      I           alktmp(i,j)*permil,
                0241      I           dictmp(i,j)*permil,
                0242      I           bt(i,j,bi,bj),
                0243      U           pH(i,j,k,bi,bj),
                0244      I           i,j,k,bi,bj,nIter0,myThid )
                0245 CAV C$TAF STORE pH(i,j,k,bi,bj)                            = dic_surf
                0246 CAV C$TAF STORE alktmp(i,j), phostmp(i,j), sitmp(i,j)      = dic_surf
                0247 #ifdef ALLOW_DEBUG
                0248               IF (debugMode) CALL DEBUG_CALL(
                0249      &           'CALC_PCO2_SOLVESAPHE',myThid)
                0250 #endif
                0251               CALL CALC_PCO2_SOLVESAPHE(
                0252      I        thetatmp(i,j),salttmp(i,j),
                0253      I        dictmp(i,j), phostmp(i,j),
                0254      I        sitmp(i,j),alktmp(i,j),
                0255      U        pH(i,j,k,bi,bj),pCO2(i,j,bi,bj),co3dummy,
                0256      I        i,j,k,bi,bj, it ,  myThid )
                0257              ELSE
                0258 C Use the original Follows et al. (2006) solver
                0259 #endif /* CARBONCHEM_SOLVESAPHE */
                0260 #ifdef ALLOW_DEBUG
                0261               IF (debugMode) CALL DEBUG_CALL(
                0262      &           'CALC_PCO2_APPROX',myThid)
                0263 #endif
                0264 
c0d1c06c15 Matt*0265              DO it=1,10
                0266               CALL CALC_PCO2_APPROX(
                0267      I        thetatmp(i,j),salttmp(i,j),
                0268      I        dictmp(i,j), phostmp(i,j),
                0269      I        sitmp(i,j),alktmp(i,j),
                0270      I        ak1(i,j,bi,bj),ak2(i,j,bi,bj),
                0271      I        ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
                0272      I        aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
                0273      I        aksi(i,j,bi,bj),akf(i,j,bi,bj),
                0274      I        ak0(i,j,bi,bj), fugf(i,j,bi,bj),
                0275      I        ff(i,j,bi,bj),
                0276      I        bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
                0277      U        pH(i,j,k,bi,bj),pCO2(i,j,bi,bj),co3dummy,
                0278      I        i,j,k,bi,bj, it ,  myThid )
                0279              ENDDO
                0280             ENDIF
e0f9a7ba0b Matt*0281 #ifdef CARBONCHEM_SOLVESAPHE
                0282              ENDIF
                0283 #endif /* CARBONCHEM_SOLVESAPHE */
c0d1c06c15 Matt*0284            ENDDO
                0285           ENDDO
                0286          ENDDO
                0287         ENDIF
                0288 
                0289 C     end bi,bj loops
                0290        ENDDO
                0291       ENDDO
                0292 
                0293 #endif /* ALLOW_BLING */
                0294 
                0295       RETURN
                0296       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/darwinproject/darwin3 by the 2.3.7-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!