pkg/dic/o2_surfforcing.F

C $Header: /u/gcmpack/MITgcm/pkg/dic/o2_surfforcing.F,v 1.10 2006/09/12 19:58:02 stephd Exp $
C $Name:  $

#include "DIC_OPTIONS.h"
#include "GCHEM_OPTIONS.h"

CBOP
C !ROUTINE: O2_SURFFORCING

C !INTERFACE: ==========================================================
      SUBROUTINE O2_SURFFORCING( PTR_O2, SGO2,
     I           bi,bj,iMin,iMax,jMin,jMax,
     I           myIter, myTime, myThid )

C !DESCRIPTION:
C Calculate the oxygen air-sea flux terms               

C !USES: ===============================================================
      IMPLICIT NONE
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "FFIELDS.h"
#ifdef DIC_BIOTIC
#include "DIC_ABIOTIC.h"
#endif


c  !INPUT PARAMETERS: ===================================================
C  myThid               :: thread number
C  myIter               :: current timestep
C  myTime               :: current time
C  PTR_O2               :: oxygen tracer field
      _RL myTime
      _RL  PTR_O2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      INTEGER iMin,iMax,jMin,jMax, bi, bj
      INTEGER myIter, myThid

c  !OUTPUT PARAMETERS: ===================================================
C  SGO2                  :: air-sea exchange of oxygen
      _RL  SGO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)

#ifdef ALLOW_PTRACERS


C !LOCAL VARIABLES: ===================================================
C I, J, K - Loop counters
      INTEGER I,J,K
C Solubility relation coefficients
      _RL SchmidtNoO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL O2sat(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL Kwexch(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL FluxO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL AtmosO2(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


      K=1

C calculate SCHMIDT NO. for O2
        DO j=jmin,jmax
          DO i=imin,imax
            IF (hFacC(i,j,k,bi,bj).NE.0.) THEN
              SchmidtNoO2(i,j) = 
     &            sox1 
     &          + sox2 * theta(i,j,k,bi,bj)
     &          + sox3 * theta(i,j,k,bi,bj)*theta(i,j,k,bi,bj)  
     &          + sox4 * theta(i,j,k,bi,bj)*theta(i,j,k,bi,bj) 
     &                *theta(i,j,k,bi,bj)

C Determine surface flux of O2 
C exchange coeff, accounting for ice cover and Schmidt no.
              pisvel(i,j,bi,bj)  =0.337*wind(i,j,bi,bj)**2/3.6d5
              Kwexch(i,j) = 
     &             pisvel(i,j,bi,bj)
     &             / sqrt(SchmidtNoO2(i,j)/660.0)
c ice influence
              Kwexch(i,j)  =(1.d0-Fice(i,j,bi,bj))*Kwexch(i,j)

              ttemp = theta(i,j,k,bi,bj)
              stemp = salt(i,j,k,bi,bj)
C determine saturation O2
C using Garcia and Gordon (1992), L&O (mistake in original???)
              aTT  = 298.15-ttemp
              aTK  = 273.15+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*1000.0

c Determine flux, inc. correction for local atmos surface pressure
cQQ PTR_O2?
              FluxO2(i,j) = maskC(i,j,k,bi,bj)*Kwexch(i,j)*
     &                     (atmosP(i,j,bi,bj)*O2sat(i,j) 
     &                      - PTR_O2(i,j,1)) 
            ELSE
              FluxO2(i,j) = 0.d0
            ENDIF


          END DO
        END DO

C update surface tendencies      
        DO j=jmin,jmax
          DO i=imin,imax
           SGO2(i,j)= maskC(i,j,1,bi,bj)*FluxO2(i,j)
     &         *recip_drF(1) * recip_hFacC(i,j,1,bi,bj)
          ENDDO
         ENDDO
#endif


        RETURN
        END

1	stephd	1.11	C $Header: /u/gcmpack/MITgcm/pkg/dic/o2_surfforcing.F,v 1.10 2006/09/12 19:58:02 stephd Exp $
2	jmc	1.6	C $Name: $
3
4	edhill	1.3	#include "DIC_OPTIONS.h"
5	stephd	1.1	#include "GCHEM_OPTIONS.h"
6
7	stephd	1.4	CBOP
8			C !ROUTINE: O2_SURFFORCING
9
10			C !INTERFACE: ==========================================================
11	stephd	1.5	SUBROUTINE O2_SURFFORCING( PTR_O2, SGO2,
12	stephd	1.1	I bi,bj,iMin,iMax,jMin,jMax,
13			I myIter, myTime, myThid )
14	stephd	1.4
15			C !DESCRIPTION:
16			C Calculate the oxygen air-sea flux terms
17
18			C !USES: ===============================================================
19	stephd	1.1	IMPLICIT NONE
20			#include "SIZE.h"
21			#include "DYNVARS.h"
22			#include "EEPARAMS.h"
23			#include "PARAMS.h"
24			#include "GRID.h"
25			#include "FFIELDS.h"
26			#ifdef DIC_BIOTIC
27			#include "DIC_ABIOTIC.h"
28			#endif
29
30
31	stephd	1.4	c !INPUT PARAMETERS: ===================================================
32			C myThid :: thread number
33			C myIter :: current timestep
34			C myTime :: current time
35			C PTR_O2 :: oxygen tracer field
36	stephd	1.1	_RL myTime
37			_RL PTR_O2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
38			INTEGER iMin,iMax,jMin,jMax, bi, bj
39	stephd	1.5	INTEGER myIter, myThid
40	stephd	1.1
41	stephd	1.4	c !OUTPUT PARAMETERS: ===================================================
42	stephd	1.5	C SGO2 :: air-sea exchange of oxygen
43			_RL SGO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
44	stephd	1.4
45	stephd	1.1	#ifdef ALLOW_PTRACERS
46
47
48	stephd	1.4	C !LOCAL VARIABLES: ===================================================
49	stephd	1.1	C I, J, K - Loop counters
50			INTEGER I,J,K
51			C Solubility relation coefficients
52			_RL SchmidtNoO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53			_RL O2sat(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54			_RL Kwexch(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
55			_RL FluxO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56			_RL AtmosO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57			_RL aTT
58			_RL aTK
59			_RL aTS
60			_RL aTS2
61			_RL aTS3
62			_RL aTS4
63			_RL aTS5
64			_RL o2s
65			_RL ttemp
66			_RL stemp
67			_RL oCnew
68	stephd	1.4	CEOP
69	stephd	1.1
70
71			K=1
72
73			C calculate SCHMIDT NO. for O2
74	stephd	1.9	DO j=jmin,jmax
75			DO i=imin,imax
76	stephd	1.1	IF (hFacC(i,j,k,bi,bj).NE.0.) THEN
77			SchmidtNoO2(i,j) =
78			& sox1
79			& + sox2 * theta(i,j,k,bi,bj)
80			& + sox3 * theta(i,j,k,bi,bj)*theta(i,j,k,bi,bj)
81			& + sox4 * theta(i,j,k,bi,bj)*theta(i,j,k,bi,bj)
82			& *theta(i,j,k,bi,bj)
83
84			C Determine surface flux of O2
85			C exchange coeff, accounting for ice cover and Schmidt no.
86	stephd	1.10	pisvel(i,j,bi,bj) =0.337wind(i,j,bi,bj)*2/3.6d5
87	stephd	1.1	Kwexch(i,j) =
88			& pisvel(i,j,bi,bj)
89			& / sqrt(SchmidtNoO2(i,j)/660.0)
90			c ice influence
91	stephd	1.7	Kwexch(i,j) =(1.d0-Fice(i,j,bi,bj))*Kwexch(i,j)
92	stephd	1.1
93			ttemp = theta(i,j,k,bi,bj)
94			stemp = salt(i,j,k,bi,bj)
95			C determine saturation O2
96			C using Garcia and Gordon (1992), L&O (mistake in original???)
97			aTT = 298.15-ttemp
98			aTK = 273.15+ttemp
99			aTS = log(aTT/aTK)
100			aTS2 = aTS*aTS
101			aTS3 = aTS2*aTS
102			aTS4 = aTS3*aTS
103			aTS5 = aTS4*aTS
104
105			oCnew = oA0 + oA1aTS + oA2aTS2 + oA3*aTS3 +
106			& oA4aTS4 + oA5aTS5
107			& + stemp(oB0 + oB1aTS + oB2aTS2 + oB3aTS3)
108			& + oC0(stempstemp)
109
110			o2s = EXP(oCnew)
111
112			c Convert from ml/l to mol/m^3
113			O2sat(i,j) = o2s/22391.6*1000.0
114
115			c Determine flux, inc. correction for local atmos surface pressure
116			cQQ PTR_O2?
117			FluxO2(i,j) = maskC(i,j,k,bi,bj)Kwexch(i,j)
118			& (atmosP(i,j,bi,bj)*O2sat(i,j)
119			& - PTR_O2(i,j,1))
120			ELSE
121			FluxO2(i,j) = 0.d0
122			ENDIF
123
124
125			END DO
126			END DO
127
128			C update surface tendencies
129	stephd	1.9	DO j=jmin,jmax
130			DO i=imin,imax
131	stephd	1.7	SGO2(i,j)= maskC(i,j,1,bi,bj)*FluxO2(i,j)
132			& recip_drF(1) recip_hFacC(i,j,1,bi,bj)
133	stephd	1.1	ENDDO
134			ENDDO
135			#endif
136
137
138			RETURN
139			END
140