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

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

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


Revision 1.1.2.2 - (hide annotations) (download)
Thu Aug 25 18:21:17 2005 UTC (15 years, 10 months ago) by dimitri
Branch: release1_50yr
Changes since 1.1.2.1: +14 -14 lines
updating verification/global_with_CFC11/code1x1/external_forcing_tr.F
and verification/global_with_CFC11/code1x1/ini_tr1.F
for release1_50yr

1 dimitri 1.1.2.2 C $Header: /u/gcmpack/MITgcm/verification/global_with_CFC11/code1x1/Attic/external_forcing_tr.F,v 1.1.2.1 2005/08/25 16:22:17 dimitri Exp $
2     C $Name: $
3 dimitri 1.1.2.1
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=0 on January 1 for FICE, XKW,
21     C | and PATM, and data.cal provides the date for cfc1112.atm.
22     C | See the OCMIP-2 CFC HOWTO for details.
23     C *==========================================================*
24     C \ev
25    
26     C !USES:
27     IMPLICIT NONE
28     C == Global data ==
29     #include "SIZE.h"
30     #include "EEPARAMS.h"
31     #include "PARAMS.h"
32     #include "GRID.h"
33     #include "DYNVARS.h"
34     #include "TR1.h"
35    
36     _RL sc_cfc, sol_cfc
37     EXTERNAL sc_cfc, sol_cfc
38    
39     C !INPUT/OUTPUT PARAMETERS:
40     C == Routine arguments ==
41     C iMin - Working range of tile for applying forcing.
42     C iMax
43     C jMin
44     C jMax
45     C kLev
46     INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
47     _RL myTime
48     INTEGER myThid
49    
50     C !LOCAL VARIABLES:
51     C == Local variables ==
52     C i,j :: Loop counters
53     C aWght, bWght :: Interpolation weights
54     INTEGER i,j,intime0,intime1,nForcingPeriods
55     INTEGER currentdate(4), tempdate(4), timediff(4)
56     _RL aWght,bWght,rdt,timeint
57     _RL fprd,fcyc,ftm,ftmold
58     _RL KW,CSAT,ALPHA,PCFC,PCFC1,PCFC2,TMP1,TMP2
59     CEOP
60    
61     IF ( kLev .EQ. 1 ) THEN
62    
63     C Find records and compute interpolation weights
64     fprd = 2629800
65     fcyc = 31557600
66     nForcingPeriods = fcyc / fprd
67     ftm = mod( myTime + fcyc - fprd / 2 , fcyc )
68     intime0 = int( ftm / fprd )
69     intime1 = mod( intime0 + 1, nForcingPeriods )
70     aWght = ( ftm - fprd * intime0 ) / ( fprd )
71     bWght = 1. - aWght
72     intime0 = intime0 + 1
73     intime1 = intime1 + 1
74    
75     C Now calculate whether it is time to update the forcing arrays
76     ftmold = mod( myTime - deltaTclock + fcyc - fprd / 2 , fcyc )
77     IF (
78     & (int(ftmold/fprd)+1) .NE. intime0
79     & .OR. myTime .EQ. startTime
80     & ) THEN
81     C If the above condition is met then we need to read in
82     C data for the period ahead and the period behind myTime.
83     _BEGIN_MASTER(myThid)
84     WRITE(*,*)
85     & 'S/R EXTERNAL_FORCING_TR: Reading new data',myTime
86    
87     #define CFC_USE_INTERP
88     C-- CFC_USE_INTERP option is useful for high-resolution integrations.
89     C For low-resolution, less that 1-deg, it's best to generate files
90     C separately because of sea-ice fraction.
91     C Caution: CFC_USE_INTERP as used is not thread-safe.
92     #ifdef CFC_USE_INTERP
93 dimitri 1.1.2.2 CALL EXF_INTERP( FiceFile,readBinaryPrec,fice0,intime0,
94     & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid )
95     CALL EXF_INTERP( FiceFile,readBinaryPrec,fice1,intime1,
96     & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid )
97     CALL EXF_INTERP( XkwFile,readBinaryPrec,xkw0,intime0,
98     & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid )
99     CALL EXF_INTERP( XkwFile,readBinaryPrec,xkw1,intime1,
100     & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid )
101     CALL EXF_INTERP( PatmFile,readBinaryPrec,patm0,intime0,
102     & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid )
103     CALL EXF_INTERP( PatmFile,readBinaryPrec,patm1,intime1,
104     & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid )
105 dimitri 1.1.2.1 #else
106     CALL MDSREADFIELD ( FiceFile, readBinaryPrec,
107     & 'RS', 1, fice0, intime0, myThid )
108     CALL MDSREADFIELD ( FiceFile, readBinaryPrec,
109     & 'RS', 1, fice1, intime1, myThid )
110     CALL MDSREADFIELD ( XkwFile, readBinaryPrec,
111     & 'RS', 1, xkw0, intime0, myThid )
112     CALL MDSREADFIELD ( XkwFile, readBinaryPrec,
113     & 'RS', 1, xkw1, intime1, myThid )
114     CALL MDSREADFIELD ( PatmFile, readBinaryPrec,
115     & 'RS', 1, patm0, intime0, myThid )
116     CALL MDSREADFIELD ( PatmFile, readBinaryPrec,
117     & 'RS', 1, patm1, intime1, myThid )
118     #endif
119    
120     _END_MASTER(myThid)
121     _EXCH_XY_R4(fice0 , myThid )
122     _EXCH_XY_R4(fice1 , myThid )
123     _EXCH_XY_R4(xkw0 , myThid )
124     _EXCH_XY_R4(xkw1 , myThid )
125     _EXCH_XY_R4(patm0 , myThid )
126     _EXCH_XY_R4(patm1 , myThid )
127     ENDIF
128    
129     C-- Interpolate in time
130     DO j=1-Oly,sNy+Oly
131     DO i=1-Olx,sNx+Olx
132     fice(i,j,bi,bj) = bWght*fice0(i,j,bi,bj) +
133     & aWght*fice1(i,j,bi,bj)
134     if( fice(i,j,bi,bj) .LT. 0.2 ) fice(i,j,bi,bj) = 0.0
135     xkw(i,j,bi,bj) = bWght* xkw0(i,j,bi,bj) +
136     & aWght* xkw1(i,j,bi,bj)
137     patm(i,j,bi,bj) = bWght*patm0(i,j,bi,bj) +
138     & aWght*patm1(i,j,bi,bj)
139     ENDDO
140     ENDDO
141    
142     C Find records and interpolation weights for PCFC
143     call cal_FullDate( 19310101, 0, tempdate, mythid )
144     call cal_GetDate( 0, mytime, currentdate, mythid )
145     call cal_SubDates( currentdate, tempdate, timediff, mythid )
146     if( timediff(1) .GT. 0 ) then
147     call cal_ToSeconds( timediff, timeint, mythid )
148     else
149     timeint=0
150     endif
151     fprd = 31557600
152     ftm = timeint + fprd / 2
153     intime0 = int( ftm / fprd )
154     intime1 = intime0 + 1
155     aWght = ( ftm - fprd * intime0 ) / ( fprd )
156     bWght = 1. - aWght
157     intime0 = intime0 + 1930
158     intime1 = intime1 + 1930
159     PCFC1 = bWght * CFCp11(intime0,1) + aWght * CFCp11(intime1,1)
160     PCFC2 = bWght * CFCp11(intime0,2) + aWght * CFCp11(intime1,2)
161    
162     C-- Forcing term: add tracer in top-layer
163     C TR1 is tracer concentration in mol/m^3
164     C OCMIP formula provides air-sea flux F in mol/m^2/s
165     C GTR1 = F / drF(1) in mol/m^3/s
166     C
167     DO j=jMin,jMax
168     DO i=iMin,iMax
169     TMP1 = theta(i,j,1,bi,bj)
170     TMP2 = salt(i,j,1,bi,bj)
171     KW = (1.-fice(i,j,bi,bj)) * xkw(i,j,bi,bj) *
172     & (sc_cfc(TMP1,11)/660)**(-1/2)
173     PCFC = pCFCw1(I,J,bi,bj) * PCFC1 +
174     & pCFCw2(I,J,bi,bj) * PCFC2
175     ALPHA = sol_cfc(TMP1,TMP2,11)
176     CSAT = ALPHA * PCFC * patm(i,j,bi,bj)
177     TMP1 = KW * ( CSAT - Tr1(i,j,1,bi,bj) )
178     surfaceTendencyTr1(i,j,bi,bj) =
179     & TMP1 * recip_drF(1) * recip_hFacC(i,j,1,bi,bj)
180     ENDDO
181     ENDDO
182    
183     ENDIF
184    
185     RETURN
186     END
187    
188     C====================================================================
189    
190     _RL FUNCTION sc_cfc(t,kn)
191    
192     c---------------------------------------------------
193     c CFC 11 and 12 schmidt number
194     c as a fonction of temperature.
195     c
196     c ref: Zheng et al (1998), JGR, vol 103,No C1
197     c
198     c t: temperature (degree Celcius)
199     c kn: = 11 for CFC-11, 12 for CFC-12
200     c
201     c J-C Dutay - LSCE
202     c---------------------------------------------------
203    
204     IMPLICIT NONE
205     INTEGER kn
206     _RL a1 ( 11: 12), a2 ( 11: 12), a3 ( 11: 12), a4 ( 11: 12)
207     _RL t
208     c
209     c coefficients with t in degre Celcius
210     c ------------------------------------
211     a1(11) = 3501.8
212     a2(11) = -210.31
213     a3(11) = 6.1851
214     a4(11) = -0.07513
215     c
216     a1(12) = 3845.4
217     a2(12) = -228.95
218     a3(12) = 6.1908
219     a4(12) = -0.067430
220     c
221    
222     sc_cfc = a1(kn) + a2(kn) * t + a3(kn) *t*t
223     & + a4(kn) *t*t*t
224    
225     RETURN
226     END
227    
228     C====================================================================
229    
230     _RL FUNCTION sol_cfc(pt,ps,kn)
231    
232     c-------------------------------------------------------------------
233     c
234     c CFC 11 and 12 Solubilities in seawater
235     c ref: Warner & Weiss (1985) , Deep Sea Research, vol32
236     c
237     c pt: temperature (degre Celcius)
238     c ps: salinity (o/oo)
239     c kn: 11 = CFC-11, 12 = CFC-12
240     c sol_cfc: in mol/m3/pptv
241     c 1 pptv = 1 part per trillion = 10^-12 atm = 1 picoatm
242     c
243     c J-C Dutay - LSCE
244     c-------------------------------------------------------------------
245    
246     _RL pt, ps,ta,d
247     _RL a1 ( 11: 12), a2 ( 11: 12), a3 ( 11: 12), a4 ( 11: 12)
248     _RL b1 ( 11: 12), b2 ( 11: 12), b3 ( 11: 12)
249    
250     INTEGER kn
251    
252     cc
253     cc coefficient for solubility in mol/l/atm
254     cc ----------------------------------------
255     c
256     c for CFC 11
257     c ----------
258     a1 ( 11) = -229.9261
259     a2 ( 11) = 319.6552
260     a3 ( 11) = 119.4471
261     a4 ( 11) = -1.39165
262     b1 ( 11) = -0.142382
263     b2 ( 11) = 0.091459
264     b3 ( 11) = -0.0157274
265     c
266     c for CFC/12
267     c ----------
268     a1 ( 12) = -218.0971
269     a2 ( 12) = 298.9702
270     a3 ( 12) = 113.8049
271     a4 ( 12) = -1.39165
272     b1 ( 12) = -0.143566
273     b2 ( 12) = 0.091015
274     b3 ( 12) = -0.0153924
275     c
276    
277     ta = ( pt + 273.16)* 0.01
278     d = ( b3 ( kn)* ta + b2 ( kn))* ta + b1 ( kn)
279     c
280     c
281     sol_cfc
282     $ = exp ( a1 ( kn)
283     $ + a2 ( kn)/ ta
284     $ + a3 ( kn)* log ( ta )
285     $ + a4 ( kn)* ta * ta + ps* d )
286     c
287     c conversion from mol/(l * atm) to mol/(m^3 * atm)
288     c ------------------------------------------------
289     sol_cfc = 1000. * sol_cfc
290     c
291     c conversion from mol/(m^3 * atm) to mol/(m3 * pptv)
292     c --------------------------------------------------
293     sol_cfc = 1.0e-12 * sol_cfc
294    
295     END

  ViewVC Help
Powered by ViewVC 1.1.22