pkg/dic/o2_surfforcing.F

C $Header: /u/gcmpack/MITgcm/pkg/dic/o2_surfforcing.F,v 1.20 2008/04/04 21:37:06 dfer Exp $
C $Name:  $

#include "DIC_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"
#include "DIC_VARS.h"

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

#ifdef ALLOW_O2

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  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 (maskC(i,j,k,bi,bj).NE.0.) THEN
              ttemp = theta(i,j,k,bi,bj)
              stemp = salt(i,j,k,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 in dic_surfforcing.F

              Kwexch(i,j) = Kwexch_Pre(i,j,bi,bj)
     &                    / 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

C Determine flux, inc. correction for local atmos surface pressure
              FluxO2(i,j) = Kwexch(i,j)*
     &                     (AtmosP(i,j,bi,bj)*O2sat(i,j) 
     &                      - PTR_O2(i,j,K)) 
            ELSE
              FluxO2(i,j) = 0. _d 0
            ENDIF


          END DO
        END DO

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


        RETURN
        END

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