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

Contents 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 - (show annotations) (download)
Thu Aug 25 18:21:17 2005 UTC (18 years, 7 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 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
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 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 #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