C $Header: /u/gcmpack/MITgcm/pkg/bling/bling_airseaflux.F,v 1.5 2017/05/11 13:32:24 mmazloff Exp $
C $Name: $
#include "BLING_OPTIONS.h"
#include "PTRACERS_OPTIONS.h"
CBOP
subroutine BLING_AIRSEAFLUX(
I PTR_DIC, PTR_ALK, PTR_O2, PTR_NO3, PTR_PO4,
O SGDIC, SGO2, FluxO2,
I bi, bj, imin, imax, jmin, jmax,
I myIter, myTime, myThid)
C =================================================================
C | subroutine bling_airseaflux
C | o Calculate the carbon and oxygen air-sea flux terms
C | Adapted from pkg/dic/dic_surfforcing.F
C | - Get atmospheric pCO2 value
C | Option 1: constant value, default 268.d-6, can be changed in
C | data.bling
C | Option 2: read 2D field using EXF pkg
C | - Update pCO2 and pH
C =================================================================
implicit none
C === Global variables ===
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "FFIELDS.h"
#ifdef ALLOW_EXF
# include "EXF_PARAM.h"
# include "EXF_FIELDS.h"
#endif
#include "BLING_VARS.h"
#ifdef ALLOW_AUTODIFF
# include "tamc.h"
#endif
C === Routine arguments ===
C myTime :: current time
C myIter :: current timestep
C myThid :: thread Id. number
_RL myTime
INTEGER myIter
INTEGER myThid
INTEGER iMin, iMax, jMin, jMax, bi, bj
C === Input ===
C PTR_DIC :: DIC tracer field
C PTR_ALK :: alkalinity tracer field
C PTR_NO3 :: nitrate tracer field
C PTR_PO4 :: phosphate tracer field
C PTR_O2 :: oxygen tracer field
_RL PTR_DIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL PTR_ALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL PTR_NO3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL PTR_PO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL PTR_O2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
C === Output ===
C SGDIC :: surface tendency of DIC due to air-sea exchange
C SGO2 :: surface tendency of O2 due to air-sea exchange
C FluxO2 :: air-sea flux of O2
C (FluxCO2 is a global variable)
_RL SGDIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL SGO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL FluxO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#ifdef ALLOW_PTRACERS
C === Local variables ===
C i,j :: Loop counters
INTEGER i,j,klev
C Number of iterations for pCO2 solvers
_RL co3dummy
_RL Kwexch_Pre (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C Solubility relation coefficients
_RL SchmidtNoDIC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL pCO2sat (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL Kwexch (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL pisvel (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C local variables for carbon chem
_RL surfalk (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL surfphos (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL surfsi (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL surftemp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL surfsalt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL surfdic (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C o2 solubility relation coefficients
_RL SchmidtNoO2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL O2sat (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL O2sat_percent(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL Kwexch_o2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL aTT
_RL aTK
_RL aTS
_RL aTS2
_RL aTS3
_RL aTS4
_RL aTS5
_RL o2s
_RL ttemp
_RL stemp
_RL oCnew
CEOP
C----------------------------------------------------------------------
C First, carbon
C----------------------------------------------------------------------
klev=1
C determine inorganic carbon chem coefficients
DO j=jmin,jmax
DO i=imin,imax
surfalk(i,j) = PTR_ALK(i,j,1)
& * maskC(i,j,1,bi,bj)
surfphos(i,j) = PTR_PO4(i,j,1)
& * maskC(i,j,1,bi,bj)
C FOR NON-INTERACTIVE Si
surfsi(i,j) = SILICA(i,j,bi,bj) * maskC(i,j,1,bi,bj)
surftemp(i,j) = theta(i,j,1,bi,bj)
surfsalt(i,j) = salt(i,j,1,bi,bj)
surfdic(i,j) = PTR_DIC(i,j,1)
ENDDO
ENDDO
CALL CARBON_COEFFS(
I surftemp,surfsalt,
I bi,bj,iMin,iMax,jMin,jMax,myThid)
DO j=jmin,jmax
DO i=imin,imax
C Compute Kwexch_Pre which is re-used for flux of O2
c Read EXF winds instead of value from file:
#ifdef ALLOW_EXF
wind(i,j,bi,bj) = wspeed(i,j,bi,bj)
#endif
C Pre-compute part of exchange coefficient: pisvel*(1-fice)
C Schmidt number is accounted for later
pisvel(i,j) = 0.337 _d 0 * wind(i,j,bi,bj)**2/3.6 _d 5
Kwexch_Pre(i,j) = pisvel(i,j)
& * (1. _d 0 - FIce(i,j,bi,bj))
ENDDO
ENDDO
c pCO2 solver...
CADJ STORE ph = comlev1, key = ikey_dynamics
C$TAF LOOP = parallel
DO j=jmin,jmax
C$TAF LOOP = parallel
DO i=imin,imax
IF ( maskC(i,j,klev,bi,bj).NE.0. _d 0 ) THEN
CALL CALC_PCO2_APPROX(
I surftemp(i,j),surfsalt(i,j),
I surfdic(i,j), surfphos(i,j),
I surfsi(i,j),surfalk(i,j),
I ak1(i,j,bi,bj),ak2(i,j,bi,bj),
I ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
I aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
I aksi(i,j,bi,bj),akf(i,j,bi,bj),
I ak0(i,j,bi,bj), fugf(i,j,bi,bj),
I ff(i,j,bi,bj),
I bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
U pH(i,j,klev,bi,bj),pCO2(i,j,bi,bj),co3dummy,
I i,j,klev,bi,bj,myIter,myThid )
ELSE
pCO2(i,j,bi,bj) = 0. _d 0
ENDIF
ENDDO
ENDDO
DO j=jmin,jmax
DO i=imin,imax
IF ( maskC(i,j,1,bi,bj).NE.0. _d 0 ) THEN
C calculate SCHMIDT NO. for CO2
SchmidtNoDIC(i,j) =
& sca1
& + sca2 * theta(i,j,1,bi,bj)
& + sca3 * theta(i,j,1,bi,bj)*theta(i,j,1,bi,bj)
& + sca4 * theta(i,j,1,bi,bj)*theta(i,j,1,bi,bj)
& *theta(i,j,1,bi,bj)
c make sure Schmidt number is not negative (will happen if temp>39C)
SchmidtNoDIC(i,j)=max(1.0 _d -2, SchmidtNoDIC(i,j))
C First determine local saturation pCO2
#ifdef USE_EXFCO2
pCO2sat(i,j) = apco2(i,j,bi,bj)
#else
pCO2sat(i,j) = bling_pCO2
#endif
c Correct for atmospheric pressure
#ifdef USE_EXF_ATMPRES
C Atm pressure in N/m2, convert to bars
pCO2sat(i,j) = pCO2sat(i,j)*(apressure(i,j,bi,bj)*0.00001)
#else
pCO2sat(i,j) = pCO2sat(i,j)*AtmosP(i,j,bi,bj)
#endif
C then account for Schmidt number
Kwexch(i,j) = Kwexch_Pre(i,j)
& / sqrt(SchmidtNoDIC(i,j)/660.0 _d 0)
C Calculate flux in terms of DIC units using K0, solubility
c Flux = kw*rho*(ff*pCO2atm-k0*FugFac*pCO2ocean)
FluxCO2(i,j,bi,bj) =
& Kwexch(i,j)*(
& ff(i,j,bi,bj)*pCO2sat(i,j) -
& pCO2(i,j,bi,bj)*fugf(i,j,bi,bj)
& *ak0(i,j,bi,bj) )
&
ELSE
FluxCO2(i,j,bi,bj) = 0. _d 0
ENDIF
C convert flux (mol kg-1 m s-1) to (mol m-2 s-1)
FluxCO2(i,j,bi,bj) = FluxCO2(i,j,bi,bj)/permil
ENDDO
ENDDO
C update tendency
DO j=jmin,jmax
DO i=imin,imax
SGDIC(i,j)= recip_drF(1)*recip_hFacC(i,j,1,bi,bj)
& *FluxCO2(i,j,bi,bj)
ENDDO
ENDDO
C----------------------------------------------------------------------
C Now oxygen
C----------------------------------------------------------------------
C calculate SCHMIDT NO. for O2
DO j=jmin,jmax
DO i=imin,imax
IF (maskC(i,j,1,bi,bj).NE.0.) THEN
ttemp = theta(i,j,1,bi,bj)
stemp = salt(i,j,1,bi,bj)
SchmidtNoO2(i,j) =
& sox1
& + sox2 * ttemp
& + sox3 * ttemp*ttemp
& + sox4 * ttemp*ttemp*ttemp
C Determine surface flux of O2
C exchange coeff accounting for ice cover and Schmidt no.
C Kwexch_Pre= pisvel*(1-fice): previously computed above
Kwexch_o2(i,j) = Kwexch_Pre(i,j)
& / sqrt(SchmidtNoO2(i,j)/660.0 _d 0)
C determine saturation O2
C using Garcia and Gordon (1992), L&O (mistake in original ?)
aTT = 298.15 _d 0 -ttemp
aTK = 273.15 _d 0 +ttemp
aTS = log(aTT/aTK)
aTS2 = aTS*aTS
aTS3 = aTS2*aTS
aTS4 = aTS3*aTS
aTS5 = aTS4*aTS
oCnew = oA0 + oA1*aTS + oA2*aTS2 + oA3*aTS3 +
& oA4*aTS4 + oA5*aTS5
& + stemp*(oB0 + oB1*aTS + oB2*aTS2 + oB3*aTS3)
& + oC0*(stemp*stemp)
o2s = EXP(oCnew)
c Convert from ml/l to mol/m^3
O2sat(i,j) = o2s/22391.6 _d 0 * 1. _d 3
cav
O2sat_percent(i,j) = PTR_O2(i,j,1)/O2sat(i,j)*100
C Determine flux, inc. correction for local atmos surface pressure
#ifdef USE_EXF_ATMPRES
C Atm pressure in N/m2, convert to bars
FluxO2(i,j) = Kwexch_o2(i,j)*(
& (apressure(i,j,bi,bj)*0.00001)
& *O2sat(i,j) - PTR_O2(i,j,1) )
#else
FluxO2(i,j) = Kwexch_o2(i,j)*
& (AtmosP(i,j,bi,bj)*O2sat(i,j)
& - PTR_O2(i,j,1))
#endif
ELSE
FluxO2(i,j) = 0. _d 0
ENDIF
ENDDO
ENDDO
C update surface tendencies
DO j=jmin,jmax
DO i=imin,imax
SGO2(i,j)= FluxO2(i,j)
& *recip_drF(1) * recip_hFacC(i,j,1,bi,bj)
ENDDO
ENDDO
_EXCH_XY_RL( pCO2, mythid)
_EXCH_XYZ_RL( pH, mythid)
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(O2sat_percent,'BLGO2SAT',0,1,2,bi,bj,
& myThid)
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
#endif /* ALLOW_PTRACER */
RETURN
END