/[MITgcm]/MITgcm/verification/global_with_CFC11/input/external_forcing_tr.F
ViewVC logotype

Annotation of /MITgcm/verification/global_with_CFC11/input/external_forcing_tr.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1.2.1 - (hide annotations) (download)
Thu May 1 18:18:29 2003 UTC (21 years ago) by dimitri
Branch: release1
Changes since 1.1: +297 -0 lines
release1_p15
o Added CFC-11 diagnostics, see
  verification/global_with_CFC11/README

1 dimitri 1.1.2.1 C $Header: $
2     C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: EXTERNAL_FORCING_TR
8     C !INTERFACE:
9     SUBROUTINE EXTERNAL_FORCING_TR(
10     I iMin, iMax, jMin, jMax,bi,bj,kLev,
11     I myTime,myThid)
12    
13     C !DESCRIPTION: \bv
14     C *==========================================================*
15     C | S/R EXTERNAL_FORCING_TR
16     C | o Contains problem specific forcing for tracer Tr1.
17     C *==========================================================*
18     C | Adds terms to gTr1 for forcing by external sources.
19     C | This routine is hardcoded for OCMIP CFC-11 experiment.
20     C | It assumes that myTime = on January 1, 1930.
21     C | See the OCMIP-2 CFC HOWTO for details.
22     C *==========================================================*
23     C \ev
24    
25     C !USES:
26     IMPLICIT NONE
27     C == Global data ==
28     #include "SIZE.h"
29     #include "EEPARAMS.h"
30     #include "PARAMS.h"
31     #include "GRID.h"
32     #include "DYNVARS.h"
33     #include "TR1.h"
34    
35     _RL sc_cfc, sol_cfc
36     EXTERNAL sc_cfc, sol_cfc
37    
38     C !INPUT/OUTPUT PARAMETERS:
39     C == Routine arguments ==
40     C iMin - Working range of tile for applying forcing.
41     C iMax
42     C jMin
43     C jMax
44     C kLev
45     INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
46     _RL myTime
47     INTEGER myThid
48    
49     C !LOCAL VARIABLES:
50     C == Local variables ==
51     C i,j :: Loop counters
52     C aWght, bWght :: Interpolation weights
53     INTEGER i,j,intime0,intime1,nForcingPeriods
54     _RL aWght,bWght,rdt
55     _RL fprd,fcyc,ftm,ftmold
56     _RL KW,CSAT,ALPHA,PCFC,PCFC1,PCFC2,TMP1,TMP2
57     CEOP
58    
59     IF ( kLev .EQ. 1 ) THEN
60    
61     C Find records and compute interpolation weights
62     fprd = 2629800
63     fcyc = 31557600
64     nForcingPeriods = fcyc / fprd
65     ftm = mod( myTime + fcyc - fprd / 2 , fcyc )
66     intime0 = int( ftm / fprd )
67     intime1 = mod( intime0 + 1, nForcingPeriods )
68     aWght = ( ftm - fprd * intime0 ) / ( fprd )
69     bWght = 1. - aWght
70     intime0 = intime0 + 1
71     intime1 = intime1 + 1
72    
73     C Now calculate whether it is time to update the forcing arrays
74     ftmold = mod( myTime - deltaTclock + fcyc - fprd / 2 , fcyc )
75     IF (
76     & (int(ftmold/fprd)+1) .NE. intime0
77     & .OR. myTime .EQ. startTime
78     & ) THEN
79     C If the above condition is met then we need to read in
80     C data for the period ahead and the period behind myTime.
81     _BEGIN_MASTER(myThid)
82     WRITE(*,*)
83     & 'S/R EXTERNAL_FORCING_TR: Reading new data',myTime
84     CALL MDSREADFIELD ( FiceFile, readBinaryPrec,
85     & 'RS', 1, fice0, intime0, myThid )
86     CALL MDSREADFIELD ( FiceFile, readBinaryPrec,
87     & 'RS', 1, fice1, intime1, myThid )
88     CALL MDSREADFIELD ( XkwFile, readBinaryPrec,
89     & 'RS', 1, xkw0, intime0, myThid )
90     CALL MDSREADFIELD ( XkwFile, readBinaryPrec,
91     & 'RS', 1, xkw1, intime1, myThid )
92     CALL MDSREADFIELD ( PatmFile, readBinaryPrec,
93     & 'RS', 1, patm0, intime0, myThid )
94     CALL MDSREADFIELD ( PatmFile, readBinaryPrec,
95     & 'RS', 1, patm1, intime1, myThid )
96     _END_MASTER(myThid)
97     _EXCH_XY_R4(fice0 , myThid )
98     _EXCH_XY_R4(fice1 , myThid )
99     _EXCH_XY_R4(xkw0 , myThid )
100     _EXCH_XY_R4(xkw1 , myThid )
101     _EXCH_XY_R4(patm0 , myThid )
102     _EXCH_XY_R4(patm1 , myThid )
103     ENDIF
104    
105     C-- Interpolate in time
106     DO j=1-Oly,sNy+Oly
107     DO i=1-Olx,sNx+Olx
108     fice(i,j,bi,bj) = bWght*fice0(i,j,bi,bj) +
109     & aWght*fice1(i,j,bi,bj)
110     if( fice(i,j,bi,bj) .LT. 0.2 ) fice(i,j,bi,bj) = 0.0
111     xkw(i,j,bi,bj) = bWght* xkw0(i,j,bi,bj) +
112     & aWght* xkw1(i,j,bi,bj)
113     patm(i,j,bi,bj) = bWght*patm0(i,j,bi,bj) +
114     & aWght*patm1(i,j,bi,bj)
115     ENDDO
116     ENDDO
117    
118     C Find records and interpolation weights for PCFC
119     fprd = 31557600
120     ftm = myTime + fprd / 2
121     intime0 = int( ftm / fprd )
122     intime1 = intime0 + 1
123     aWght = ( ftm - fprd * intime0 ) / ( fprd )
124     bWght = 1. - aWght
125     intime0 = intime0 + 1930
126     intime1 = intime1 + 1930
127     PCFC1 = bWght * CFCp11(intime0,1) + aWght * CFCp11(intime1,1)
128     PCFC2 = bWght * CFCp11(intime0,2) + aWght * CFCp11(intime1,2)
129    
130     C-- Forcing term: add tracer in top-layer
131     C TR1 is tracer concentration in mol/m^3
132     C OCMIP formula provides air-sea flux F in mol/m^2/s
133     C GTR1 = F / drF(1) in mol/m^3/s
134     C
135     DO j=jMin,jMax
136     DO i=iMin,iMax
137     TMP1 = theta(i,j,1,bi,bj)
138     TMP2 = salt(i,j,1,bi,bj)
139     KW = (1.-fice(i,j,bi,bj)) * xkw(i,j,bi,bj) *
140     & (sc_cfc(TMP1,11)/660)**(-1/2)
141     PCFC = pCFCw1(I,J,bi,bj) * PCFC1 +
142     & pCFCw2(I,J,bi,bj) * PCFC2
143     ALPHA = sol_cfc(TMP1,TMP2,11)
144     CSAT = ALPHA * PCFC * patm(i,j,bi,bj)
145     TMP1 = KW * ( CSAT - Tr1(i,j,1,bi,bj) )
146     gTr1(i,j,1,bi,bj) = gTr1(i,j,1,bi,bj) +
147     & TMP1 * recip_drF(1) * recip_hFacC(i,j,1,bi,bj)
148     ENDDO
149     ENDDO
150    
151     ENDIF
152    
153     RETURN
154     END
155    
156     C====================================================================
157    
158     c_ ---------------------------------------------------------------------
159     c_ RCS lines preceded by "c_ "
160     c_ ---------------------------------------------------------------------
161     c_
162     c_ $Source: /home/orr/ocmip/web/OCMIP/phase2/simulations/CFC/boundcond/RCS/sc_cfc.f,v $
163     c_ $Revision: 1.1 $
164     c_ $Date: 1998/07/07 15:22:00 $ ; $State: Exp $
165     c_ $Author: orr $ ; $Locker: $
166     c_
167     c_ ---------------------------------------------------------------------
168     c_ $Log: sc_cfc.f,v $
169     c_ Revision 1.1 1998/07/07 15:22:00 orr
170     c_ Initial revision
171     c_
172     c_ ---------------------------------------------------------------------
173     c_
174     _RL FUNCTION sc_cfc(t,kn)
175     c---------------------------------------------------
176     c CFC 11 and 12 schmidt number
177     c as a fonction of temperature.
178     c
179     c ref: Zheng et al (1998), JGR, vol 103,No C1
180     c
181     c t: temperature (degree Celcius)
182     c kn: = 11 for CFC-11, 12 for CFC-12
183     c
184     c J-C Dutay - LSCE
185     c---------------------------------------------------
186     IMPLICIT NONE
187     INTEGER kn
188     _RL a1 ( 11: 12), a2 ( 11: 12), a3 ( 11: 12), a4 ( 11: 12)
189     _RL t
190     c
191     c coefficients with t in degre Celcius
192     c ------------------------------------
193     a1(11) = 3501.8
194     a2(11) = -210.31
195     a3(11) = 6.1851
196     a4(11) = -0.07513
197     c
198     a1(12) = 3845.4
199     a2(12) = -228.95
200     a3(12) = 6.1908
201     a4(12) = -0.067430
202     c
203    
204     sc_cfc = a1(kn) + a2(kn) * t + a3(kn) *t*t
205     & + a4(kn) *t*t*t
206    
207     RETURN
208     END
209    
210     C====================================================================
211    
212     c_ ---------------------------------------------------------------------
213     c_ RCS lines preceded by "c_ "
214     c_ ---------------------------------------------------------------------
215     c_
216     c_ $Source: /www/html/ipsl/OCMIP/phase2/simulations/CFC/boundcond/RCS/sol_cfc.f,v $
217     c_ $Revision: 1.2 $
218     c_ $Date: 1998/07/17 07:37:02 $ ; $State: Exp $
219     c_ $Author: jomce $ ; $Locker: $
220     c_
221     c_ ---------------------------------------------------------------------
222     c_ $Log: sol_cfc.f,v $
223     c_ Revision 1.2 1998/07/17 07:37:02 jomce
224     c_ Fixed slight bug in units conversion: converted 1.0*e-12 to 1.0e-12
225     c_ following warning from Matthew Hecht at NCAR.
226     c_
227     c_ Revision 1.1 1998/07/07 15:22:00 orr
228     c_ Initial revision
229     c_
230     c_ ---------------------------------------------------------------------
231     c_
232     _RL FUNCTION sol_cfc(pt,ps,kn)
233     c-------------------------------------------------------------------
234     c
235     c CFC 11 and 12 Solubilities in seawater
236     c ref: Warner & Weiss (1985) , Deep Sea Research, vol32
237     c
238     c pt: temperature (degre Celcius)
239     c ps: salinity (o/oo)
240     c kn: 11 = CFC-11, 12 = CFC-12
241     c sol_cfc: in mol/m3/pptv
242     c 1 pptv = 1 part per trillion = 10^-12 atm = 1 picoatm
243    
244     c
245     c J-C Dutay - LSCE
246     c-------------------------------------------------------------------
247    
248     _RL pt, ps,ta,d
249     _RL a1 ( 11: 12), a2 ( 11: 12), a3 ( 11: 12), a4 ( 11: 12)
250     _RL b1 ( 11: 12), b2 ( 11: 12), b3 ( 11: 12)
251    
252     INTEGER kn
253    
254     cc
255     cc coefficient for solubility in mol/l/atm
256     cc ----------------------------------------
257     c
258     c for CFC 11
259     c ----------
260     a1 ( 11) = -229.9261
261     a2 ( 11) = 319.6552
262     a3 ( 11) = 119.4471
263     a4 ( 11) = -1.39165
264     b1 ( 11) = -0.142382
265     b2 ( 11) = 0.091459
266     b3 ( 11) = -0.0157274
267     c
268     c for CFC/12
269     c ----------
270     a1 ( 12) = -218.0971
271     a2 ( 12) = 298.9702
272     a3 ( 12) = 113.8049
273     a4 ( 12) = -1.39165
274     b1 ( 12) = -0.143566
275     b2 ( 12) = 0.091015
276     b3 ( 12) = -0.0153924
277     c
278    
279     ta = ( pt + 273.16)* 0.01
280     d = ( b3 ( kn)* ta + b2 ( kn))* ta + b1 ( kn)
281     c
282     c
283     sol_cfc
284     $ = exp ( a1 ( kn)
285     $ + a2 ( kn)/ ta
286     $ + a3 ( kn)* log ( ta )
287     $ + a4 ( kn)* ta * ta + ps* d )
288     c
289     c conversion from mol/(l * atm) to mol/(m^3 * atm)
290     c ------------------------------------------------
291     sol_cfc = 1000. * sol_cfc
292     c
293     c conversion from mol/(m^3 * atm) to mol/(m3 * pptv)
294     c --------------------------------------------------
295     sol_cfc = 1.0e-12 * sol_cfc
296    
297     END

  ViewVC Help
Powered by ViewVC 1.1.22