pkg/bling/bling_airseaflux.F

c0d1c06c15 Matt*0001 #include "BLING_OPTIONS.h"
c8c54387f0 Matt*0002 #ifdef ALLOW_EXF
                0003 # include "EXF_OPTIONS.h"
                0004 #endif
a284455135 Matt*0005 #ifdef ALLOW_AUTODIFF
                0006 # include "AUTODIFF_OPTIONS.h"
                0007 #endif
c0d1c06c15 Matt*0008
                0009 CBOP
a284455135 Matt*0010       SUBROUTINE BLING_AIRSEAFLUX(
e0f9a7ba0b Matt*0011      I           PTR_DIC, PTR_ALK, PTR_O2, PTR_PO4,
                0012 #ifdef USE_SIBLING
                0013      I           PTR_SI,
                0014 #endif
                0015      O           SGDIC, SGO2, FluxO2,
c0d1c06c15 Matt*0016      I           bi, bj, imin, imax, jmin, jmax,
e0f9a7ba0b Matt*0017      I           myTime, myIter, myThid)
c0d1c06c15 Matt*0018
                0019 C     =================================================================
                0020 C     | subroutine bling_airseaflux
                0021 C     | o Calculate the carbon and oxygen air-sea flux terms
4ac06494d5 Matt*0022 C     |   Adapted from pkg/dic/dic_surfforcing.F
c0d1c06c15 Matt*0023 C     | - Get atmospheric pCO2 value
e0f9a7ba0b Matt*0024 C     |   Option 1: constant value, default 268.d-6, can be changed
                0025 C     |             in data.bling
c0d1c06c15 Matt*0026 C     |   Option 2: read 2D field using EXF pkg
4ac06494d5 Matt*0027 C     | - Update pCO2 and pH
c0d1c06c15 Matt*0028 C     =================================================================
                0029
e0f9a7ba0b Matt*0030       IMPLICIT NONE
                0031
c0d1c06c15 Matt*0032 C     === Global variables ===
                0033 #include "SIZE.h"
                0034 #include "DYNVARS.h"
                0035 #include "EEPARAMS.h"
                0036 #include "PARAMS.h"
                0037 #include "GRID.h"
                0038 #include "FFIELDS.h"
                0039 #ifdef ALLOW_EXF
9f0da36f91 Jean*0040 # include "EXF_INTERP_SIZE.h"
c0d1c06c15 Matt*0041 #endif
079948e6a6 Matt*0042 #include "BLING_VARS.h"
a284455135 Matt*0043 #ifdef ALLOW_AUTODIFF_TAMC
c0d1c06c15 Matt*0044 # include "tamc.h"
                0045 #endif
                0046
                0047 C     === Routine arguments ===
                0048 C     myTime           :: current time
                0049 C     myIter           :: current timestep
                0050 C     myThid           :: thread Id. number
                0051       _RL myTime
                0052       INTEGER myIter
                0053       INTEGER myThid
                0054       INTEGER iMin, iMax, jMin, jMax, bi, bj
                0055 C     === Input ===
                0056 C     PTR_DIC          :: DIC tracer field
                0057 C     PTR_ALK          :: alkalinity tracer field
                0058 C     PTR_PO4          :: phosphate tracer field
                0059 C     PTR_O2           :: oxygen tracer field
                0060       _RL  PTR_DIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0061       _RL  PTR_ALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0062       _RL  PTR_PO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0063       _RL  PTR_O2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
e0f9a7ba0b Matt*0064 #ifdef USE_SIBLING
                0065       _RL  PTR_SI (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0066 #endif
c0d1c06c15 Matt*0067 C     === Output ===
e0f9a7ba0b Matt*0068 C     SGDIC            :: surface tendency of DIC due to air-sea exchange
                0069 C     SGO2             :: surface tendency of O2 due to air-sea exchange
                0070 C     FluxO2           :: air-sea flux of O2
                0071 C     (FluxCO2 is a global variable)
c0d1c06c15 Matt*0072       _RL  SGDIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0073       _RL  SGO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
b00a067069 Matt*0074       _RL  FluxO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c0d1c06c15 Matt*0075
                0076 #ifdef ALLOW_PTRACERS
                0077
                0078 C     === Local variables ===
                0079 C     i,j              :: Loop counters
                0080       INTEGER i,j,klev
                0081 C Number of iterations for pCO2 solvers
                0082       _RL co3dummy
b00a067069 Matt*0083       _RL Kwexch_Pre   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c0d1c06c15 Matt*0084 C Solubility relation coefficients
b00a067069 Matt*0085       _RL SchmidtNoDIC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0086       _RL pCO2sat      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0087       _RL Kwexch       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0088       _RL pisvel       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c0d1c06c15 Matt*0089 C local variables for carbon chem
b00a067069 Matt*0090       _RL surfalk      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0091       _RL surfphos     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0092       _RL surfsi       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0093       _RL surftemp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0094       _RL surfsalt     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0095       _RL surfdic      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c0d1c06c15 Matt*0096 C o2 solubility relation coefficients
b00a067069 Matt*0097       _RL SchmidtNoO2  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0098       _RL O2sat        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0099       _RL O2sat_percent(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0100       _RL Kwexch_o2    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c0d1c06c15 Matt*0101       _RL aTT
                0102       _RL aTK
                0103       _RL aTS
                0104       _RL aTS2
                0105       _RL aTS3
                0106       _RL aTS4
                0107       _RL aTS5
                0108       _RL o2s
                0109       _RL ttemp
                0110       _RL stemp
                0111       _RL oCnew
7c50f07931 Mart*0112 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0113 C     tkey :: tape key (tile dependent)
                0114       INTEGER tkey
7c50f07931 Mart*0115 #endif
c0d1c06c15 Matt*0116 CEOP
                0117
a284455135 Matt*0118 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0119       tkey = bi + (bj - 1)*nSx + (ikey_dynamics - 1)*nSx*nSy
a284455135 Matt*0120 #endif
                0121
c0d1c06c15 Matt*0122 C----------------------------------------------------------------------
                0123 C First, carbon
                0124 C----------------------------------------------------------------------
                0125         klev=1
                0126 C determine inorganic carbon chem coefficients
                0127         DO j=jmin,jmax
                0128          DO i=imin,imax
                0129
                0130              surfalk(i,j)  = PTR_ALK(i,j,1)
                0131      &                          * maskC(i,j,1,bi,bj)
                0132              surfphos(i,j) = PTR_PO4(i,j,1)
                0133      &                          * maskC(i,j,1,bi,bj)
                0134
                0135 C FOR NON-INTERACTIVE Si
                0136              surftemp(i,j) = theta(i,j,1,bi,bj)
                0137              surfsalt(i,j) = salt(i,j,1,bi,bj)
                0138              surfdic(i,j)  = PTR_DIC(i,j,1)
e0f9a7ba0b Matt*0139 #ifdef USE_SIBLING
                0140              surfsi(i,j)   = PTR_SI(i,j,1)
                0141 #else
                0142              surfsi(i,j)   = silica(i,j,bi,bj) * maskC(i,j,1,bi,bj)
                0143 #endif
c0d1c06c15 Matt*0144
                0145           ENDDO
                0146          ENDDO
                0147
e0f9a7ba0b Matt*0148 #ifdef CARBONCHEM_SOLVESAPHE
                0149 #ifdef ALLOW_DEBUG
                0150       IF (debugMode) CALL DEBUG_CALL('CARBON_COEFFS_SOLVESAPHE',myThid)
                0151 #endif
                0152         CALL DIC_COEFFS_SURF(
                0153      I                       surftemp,surfsalt,
                0154      I                       bi,bj,iMin,iMax,jMin,jMax,myThid)
                0155 #else
                0156 #ifdef ALLOW_DEBUG
                0157       IF (debugMode) CALL DEBUG_CALL('CARBON_COEFFS',myThid)
                0158 #endif
c0d1c06c15 Matt*0159          CALL CARBON_COEFFS(
                0160      I                       surftemp,surfsalt,
                0161      I                       bi,bj,iMin,iMax,jMin,jMax,myThid)
e0f9a7ba0b Matt*0162 #endif
c0d1c06c15 Matt*0163
                0164        DO j=jmin,jmax
                0165         DO i=imin,imax
                0166
                0167 C Pre-compute part of exchange coefficient: pisvel*(1-fice)
e0f9a7ba0b Matt*0168 C which is re-used for flux of O2
c0d1c06c15 Matt*0169 C Schmidt number is accounted for later
e0f9a7ba0b Matt*0170
c0d1c06c15 Matt*0171               pisvel(i,j) = 0.337 _d 0 * wind(i,j,bi,bj)**2/3.6 _d 5
e0f9a7ba0b Matt*0172
                0173               Kwexch_Pre(i,j) = pisvel(i,j)* (1. _d 0 - FIce(i,j,bi,bj))
c0d1c06c15 Matt*0174
                0175         ENDDO
                0176        ENDDO
                0177
                0178 c pCO2 solver...
                0179
a284455135 Matt*0180 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0181 CADJ STORE ph(:,:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
a284455135 Matt*0182 #endif
c0d1c06c15 Matt*0183
                0184 C$TAF LOOP = parallel
                0185        DO j=jmin,jmax
                0186 C$TAF LOOP = parallel
                0187         DO i=imin,imax
                0188
                0189           IF ( maskC(i,j,klev,bi,bj).NE.0. _d 0 ) THEN
e0f9a7ba0b Matt*0190 #ifdef CARBONCHEM_SOLVESAPHE
                0191              IF ( selectPHsolver.GT.0 ) THEN
                0192 C Use Munhoven (2013) Solvesaphe routine to calculate pH and pCO2
                0193 #ifdef ALLOW_DEBUG
                0194                 IF (debugMode) CALL DEBUG_CALL(
                0195      &     'CALC_PCO2_SOLVESAPHE from DIC_SURFFORCING',myThid)
                0196 #endif
                0197                 CALL CALC_PCO2_SOLVESAPHE(
                0198      I        surftemp(i,j),surfsalt(i,j),
                0199      I        surfdic(i,j), surfphos(i,j),
                0200      I        surfsi(i,j),surfalk(i,j),
                0201      U        pH(i,j,klev,bi,bj),pCO2(i,j,bi,bj),co3dummy,
                0202      I        i,j,klev,bi,bj,myIter, myThid )
                0203              ELSE
                0204 C Use the original Follows et al. (2006) solver
                0205 #endif /* CARBONCHEM_SOLVESAPHE */
                0206 #ifdef ALLOW_DEBUG
                0207                 IF (debugMode) CALL DEBUG_CALL(
                0208      &        'CALC_PCO2_APPROX',myThid)
                0209 #endif
                0210
c0d1c06c15 Matt*0211             CALL CALC_PCO2_APPROX(
                0212      I        surftemp(i,j),surfsalt(i,j),
                0213      I        surfdic(i,j), surfphos(i,j),
                0214      I        surfsi(i,j),surfalk(i,j),
                0215      I        ak1(i,j,bi,bj),ak2(i,j,bi,bj),
                0216      I        ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
                0217      I        aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
                0218      I        aksi(i,j,bi,bj),akf(i,j,bi,bj),
                0219      I        ak0(i,j,bi,bj), fugf(i,j,bi,bj),
                0220      I        ff(i,j,bi,bj),
                0221      I        bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
                0222      U        pH(i,j,klev,bi,bj),pCO2(i,j,bi,bj),co3dummy,
                0223      I        i,j,klev,bi,bj,myIter,myThid )
e0f9a7ba0b Matt*0224
                0225 #ifdef CARBONCHEM_SOLVESAPHE
                0226              ENDIF
                0227 #endif /* CARBONCHEM_SOLVESAPHE */
                0228
c0d1c06c15 Matt*0229           ELSE
                0230             pCO2(i,j,bi,bj) = 0. _d 0
                0231           ENDIF
                0232
                0233         ENDDO
                0234        ENDDO
b26a461de7 Mart*0235 #ifdef ALLOW_AUTODIFF_TAMC
                0236 CADJ STORE pCO2(:,:,bi,bj) = comlev1_bibj, key = tkey, kind = isbyte
                0237 #endif
c0d1c06c15 Matt*0238
                0239        DO j=jmin,jmax
                0240         DO i=imin,imax
                0241
                0242           IF ( maskC(i,j,1,bi,bj).NE.0. _d 0 ) THEN
                0243 C calculate SCHMIDT NO. for CO2
                0244               SchmidtNoDIC(i,j) =
                0245      &            sca1
                0246      &          + sca2 * theta(i,j,1,bi,bj)
                0247      &          + sca3 * theta(i,j,1,bi,bj)*theta(i,j,1,bi,bj)
                0248      &          + sca4 * theta(i,j,1,bi,bj)*theta(i,j,1,bi,bj)
                0249      &                *theta(i,j,1,bi,bj)
                0250 c make sure Schmidt number is not negative (will happen if temp>39C)
e0f9a7ba0b Matt*0251              SchmidtNoDIC(i,j) = max(1.0 _d -2, SchmidtNoDIC(i,j))
c0d1c06c15 Matt*0252
                0253 C First determine local saturation pCO2
                0254 c Correct for atmospheric pressure
a284455135 Matt*0255               pCO2sat(i,j) = apco2(i,j,bi,bj) * AtmosP(i,j,bi,bj)
c0d1c06c15 Matt*0256
                0257 C then account for Schmidt number
                0258               Kwexch(i,j) = Kwexch_Pre(i,j)
                0259      &                    / sqrt(SchmidtNoDIC(i,j)/660.0 _d 0)
                0260
                0261 C Calculate flux in terms of DIC units using K0, solubility
                0262 c Flux = kw*rho*(ff*pCO2atm-k0*FugFac*pCO2ocean)
                0263                FluxCO2(i,j,bi,bj) =
                0264      &          Kwexch(i,j)*(
                0265      &            ff(i,j,bi,bj)*pCO2sat(i,j) -
                0266      &            pCO2(i,j,bi,bj)*fugf(i,j,bi,bj)
                0267      &            *ak0(i,j,bi,bj) )
                0268      &
                0269           ELSE
                0270               FluxCO2(i,j,bi,bj) = 0. _d 0
                0271           ENDIF
e0f9a7ba0b Matt*0272
c0d1c06c15 Matt*0273 C convert flux (mol kg-1 m s-1) to (mol m-2 s-1)
                0274           FluxCO2(i,j,bi,bj) = FluxCO2(i,j,bi,bj)/permil
                0275
                0276           ENDDO
                0277          ENDDO
                0278
                0279 C update tendency
                0280          DO j=jmin,jmax
                0281           DO i=imin,imax
                0282            SGDIC(i,j)= recip_drF(1)*recip_hFacC(i,j,1,bi,bj)
                0283      &              *FluxCO2(i,j,bi,bj)
                0284           ENDDO
                0285          ENDDO
                0286
                0287 C----------------------------------------------------------------------
                0288 C Now oxygen
                0289 C----------------------------------------------------------------------
                0290
                0291 C calculate SCHMIDT NO. for O2
                0292         DO j=jmin,jmax
                0293           DO i=imin,imax
                0294             IF (maskC(i,j,1,bi,bj).NE.0.) THEN
                0295               ttemp = theta(i,j,1,bi,bj)
                0296               stemp = salt(i,j,1,bi,bj)
                0297
                0298               SchmidtNoO2(i,j) =
                0299      &            sox1
                0300      &          + sox2 * ttemp
                0301      &          + sox3 * ttemp*ttemp
                0302      &          + sox4 * ttemp*ttemp*ttemp
                0303
                0304 C Determine surface flux of O2
                0305 C exchange coeff accounting for ice cover and Schmidt no.
                0306 C Kwexch_Pre= pisvel*(1-fice): previously computed above
                0307
                0308               Kwexch_o2(i,j) = Kwexch_Pre(i,j)
                0309      &                    / sqrt(SchmidtNoO2(i,j)/660.0 _d 0)
                0310
                0311 C determine saturation O2
                0312 C using Garcia and Gordon (1992), L&O (mistake in original ?)
                0313               aTT  = 298.15 _d 0 -ttemp
                0314               aTK  = 273.15 _d 0 +ttemp
                0315               aTS  = log(aTT/aTK)
                0316               aTS2 = aTS*aTS
                0317               aTS3 = aTS2*aTS
                0318               aTS4 = aTS3*aTS
                0319               aTS5 = aTS4*aTS
                0320
                0321               oCnew  = oA0 + oA1*aTS + oA2*aTS2 + oA3*aTS3 +
                0322      &            oA4*aTS4 + oA5*aTS5
                0323      &          + stemp*(oB0 + oB1*aTS + oB2*aTS2 + oB3*aTS3)
                0324      &          + oC0*(stemp*stemp)
                0325
                0326               o2s = EXP(oCnew)
                0327
                0328 c Convert from ml/l to mol/m^3
                0329               O2sat(i,j) = o2s/22391.6 _d 0 * 1. _d 3
e0f9a7ba0b Matt*0330 c Calculate percent saturation for diagnostic
b00a067069 Matt*0331               O2sat_percent(i,j) = PTR_O2(i,j,1)/O2sat(i,j)*100
                0332
e0f9a7ba0b Matt*0333 C Determine flux, incl. correction for local atmos surface pressure
c0d1c06c15 Matt*0334               FluxO2(i,j) = Kwexch_o2(i,j)*
                0335      &                     (AtmosP(i,j,bi,bj)*O2sat(i,j)
                0336      &                      - PTR_O2(i,j,1))
                0337             ELSE
                0338               FluxO2(i,j) = 0. _d 0
                0339             ENDIF
                0340
                0341           ENDDO
                0342         ENDDO
                0343
                0344 C update surface tendencies
                0345         DO j=jmin,jmax
                0346           DO i=imin,imax
                0347            SGO2(i,j)= FluxO2(i,j)
                0348      &         *recip_drF(1) * recip_hFacC(i,j,1,bi,bj)
                0349           ENDDO
                0350         ENDDO
                0351
e0f9a7ba0b Matt*0352         _EXCH_XY_RL( pCO2, myThid )
                0353         _EXCH_XYZ_RL( pH, myThid )
c0d1c06c15 Matt*0354
b00a067069 Matt*0355 #ifdef ALLOW_DIAGNOSTICS
                0356       IF ( useDiagnostics ) THEN
                0357         CALL DIAGNOSTICS_FILL(O2sat_percent,'BLGO2SAT',0,1,2,bi,bj,
                0358      &       myThid)
                0359       ENDIF
                0360 #endif /* ALLOW_DIAGNOSTICS */
                0361
c0d1c06c15 Matt*0362 #endif /* ALLOW_PTRACER */
                0363
                0364         RETURN
                0365         END