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

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

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

revision 1.1 by dimitri, Thu Aug 25 16:22:17 2005 UTC revision 1.1.2.1 by dimitri, Thu Aug 25 16:22:17 2005 UTC
# Line 0  Line 1 
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=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 NEW_INTERP( FiceFile,readBinaryPrec,fice0,intime0,
94         &           xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,mythid )
95                CALL NEW_INTERP( FiceFile,readBinaryPrec,fice1,intime1,
96         &           xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,mythid )
97                CALL NEW_INTERP( XkwFile,readBinaryPrec,xkw0,intime0,
98         &           xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,mythid )
99                CALL NEW_INTERP( XkwFile,readBinaryPrec,xkw1,intime1,
100         &           xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,mythid )
101                CALL NEW_INTERP( PatmFile,readBinaryPrec,patm0,intime0,
102         &           xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,mythid )
103                CALL NEW_INTERP( PatmFile,readBinaryPrec,patm1,intime1,
104         &           xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,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
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.1.2.1
  ViewVC Help
Powered by ViewVC 1.1.22