/[MITgcm]/MITgcm_contrib/AITCZ/code/mitphys_do_inphys.F
ViewVC logotype

Annotation of /MITgcm_contrib/AITCZ/code/mitphys_do_inphys.F

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


Revision 1.1 - (hide annotations) (download)
Wed Aug 20 15:24:59 2003 UTC (21 years, 11 months ago) by czaja
Branch: MAIN
CVS Tags: HEAD
Initial creation of Arnaud's simple coupled simulation.

1 czaja 1.1 C $Header: /u/gcmpack/models/MITgcmUV/pkg/aim/aim_do_inphys.F,v 1.3.2.1 2001/04/05 03:58:19 jamous Exp $
2     C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5    
6     SUBROUTINE MITPHYS_INIT( myThid )
7     C /==================================================================\
8     C | S/R MITPHYS_INIT |
9     C |==================================================================|
10     C | Interface between MITPHYS atmospheric physics package and the |
11     C | dynamical model for initialisation. |
12     C \==================================================================/
13     IMPLICIT NONE
14    
15     C -------------- Global variables ------------------------------------
16     #include "atparam.h"
17     #include "atparam1.h"
18    
19    
20     #include "EEPARAMS.h"
21     #include "PARAMS.h"
22     #include "DYNVARS.h"
23     #include "GRID.h"
24    
25     #include "MITPHYS_DIAGS.h"
26     #include "MITPHYS_PARAMS.h"
27     INTEGER NGP
28     INTEGER NLON
29     INTEGER NLAT
30     INTEGER NLEV
31     PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT )
32    
33     #include "com_mitphysvar.h"
34     #include "com_forcing1.h"
35     C == Routine arguments ==
36     C myThid - Number of this instance
37     INTEGER myThid
38    
39     #ifdef ALLOW_MITPHYS
40     C == Local variables ==
41     C FSG - Cell mid-point in vertical
42     C HSG - Cell face in vertical
43     C RLAT - Call mid-point latitude
44     C pGround - Lower boundary pressure
45     C J, K, bi,bj - Loop counters
46     _RL FSG( Nr)
47     _RL HSG( Nr+1)
48     _RL RLAT(sNy)
49     _RL pGround
50     INTEGER I,I2,J, K
51     INTEGER Katm
52     INTEGER bi,bj
53     _RL LAT_CYCLE, DIST
54     _RS FMASK2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
55    
56    
57     pground = 1012.5D2
58    
59     c omp: Single thread per processor:
60     bi = 1
61     bj = 1
62    
63     DO K=1,Nr
64     Katm = K
65     FSG(Katm ) = rC(K)/pGround
66     HSG(Katm+1) = rF(K)/pGround
67     ENDDO
68     HSG(1) = rF(Nr+1)/pGround
69    
70    
71     DO J = 1, sNy
72     DO I = 1, sNx
73     I2 = (sNx)*(J-1)+I
74     if (UsingSphericalPolarGrid) THEN
75     C latitude and logitude, assuming a spherical coordinate system
76     LAT_G(I2) = yC(I,J,bi,bj)
77     LON_G(I2) = xC(I,J,bi,bj)
78     ELSE
79     C For Cartesian Grid, need to specify - equatorial value for now...
80     LAT_G(I2) = 0.0
81     LON_G(I2) = 0.0
82     END IF
83     CBMFG1 (I2)= 0.0
84     END DO
85     END DO
86    
87    
88     C DO J=1,sNy
89     C RLAT(J) = yC(1,J,1,1)*deg2rad
90     C ENDDO
91    
92     CALL MITPHYS_READPARMS(myThid)
93    
94    
95     C 2. SST / STL
96     C modif omp: the background SST is defined thourgh the MITPHYS namelist
97     C modif acz: introduce analytical STL
98    
99     C omp: for spherical geometry
100     IF (UsingSphericalPolarGrid) THEN
101    
102     LAT_CYCLE = 1.0
103     DO J=1,sNy
104     DO I=1,sNx
105     I2 = (sNx)*(J-1)+I
106    
107     SST_REF(I2) = SST_BACK
108     STL_REF(I2) = SST_BACK
109    
110     C (acz): NTA positive SST anomaly
111    
112     dist = SQRT ( ( LAT_G(I2) - LAT_PERT_NTA) ** 2)
113    
114     if (dist .LE. DEL_LAT_NTA) then
115     SST_REF(I2) = SST_REF(I2) + DELT_PERT_NTA
116     & * cos ( pi * dist / 2./DEL_LAT_NTA ) **2
117     end if
118    
119    
120     C Equator to pole SST gradient
121    
122     SST_REF(I2) = SST_REF(I2) - DELT_EQ_PO *
123     & sin( LAT_G(I2) * pi / 180. ) **2
124    
125    
126     C Equator to northern boundary SST gradient
127    
128     SST_REF(I2) = SST_REF(I2) - DELT_EQ_BND *
129     & sin( LAT_G(I2) / abs(Phimin)) **2
130    
131     C Lindzen and Hou type gradient:
132    
133     SST_REF(I2) = SST_REF(I2) - DELT_LH * (
134     & 3 * ( sin ( PI * LAT_G(I2) / 180.d0 )
135     & - sin(PI * LAT_CYCLE * LAT_LH / 180.d0)
136     & ) ** 2)
137    
138    
139     CC (acz) The SST used for TOSA1:
140     CC Zonally symmetric SST perturbation
141    
142    
143     dist = SQRT ( ( LAT_G(I2) - LAT_PERT*LAT_CYCLE) ** 2)
144    
145     if (dist .LE. DEL_LAT) then
146     SST_REF(I2) = SST_REF(I2) + DELT_PERT
147     & * cos ( pi * dist / 2./DEL_LAT ) **2
148     end if
149    
150     CC (acz) STL used for TOSA1:
151     CC
152     CC 1 - Warm ring centered at the equator with same
153     CC meridional extent (DEL_LAT) than the above (TOSA1) SST
154     CC The profile is further time-dependent:
155     C strongest at equinox, zero at solstice
156     C so it is set at zero for this initialization
157    
158     dist = ABS ( LAT_G(I2) )
159    
160     if (dist .LE. DEL_LAT) then
161     STL_REF(I2) = STL_REF(I2) + DELT_STL_EQU
162     & * cos ( pi * dist / 2./DEL_LAT ) **2
163     & * cos ( 2 * pi * (0. / 24./ 3600.)
164     & / 360. + 0.5 * pi ) **2
165     end if
166    
167     C 2- Opposite sign variations north and south of equator
168     C North is warm at t=0,
169    
170     dist = ABS ( LAT_G(I2) )
171    
172     if (dist .LE. DEL_LAT_STL) then
173     STL_REF(I2) = STL_REF(I2) + DELT_STL_HEM
174     & * sin ( pi * LAT_G(I2) / DEL_LAT_STL )
175     end if
176    
177     ENDDO
178     ENDDO
179    
180     ELSE
181    
182     C Cartesian Grid
183     CC (acz) nothing for the moment ...
184    
185     END IF
186    
187     C 3. SOIL MOISTURE
188     C (acz): start with full bucket
189    
190     C for spherical geometry
191     IF (UsingSphericalPolarGrid) THEN
192    
193     DO J=1,sNy
194     DO I=1,sNx
195     I2 = (sNx)*(J-1)+I
196     BUCKET(I2) = BKT0
197     END DO
198     END DO
199    
200     CC(acz)
201     CC CALL WRITE_FLD_XY_RS('BKTini','0000000000',BUCKET,0,myThid)
202     CC write(6,*) 'end of do_inphys: BKT0=', BKT0
203    
204     ELSE
205    
206     C Cartesian Grid: nothing is done for now...
207    
208     END IF
209    
210    
211    
212     C 4. LAND-SEA MASK
213     C (acz): FMASK1 is defined in com_forcing1.h It is determined analytically here
214     C FMASK1 = 0 over ocean
215     C FMASK1 = 1 over land
216    
217     C for spherical geometry
218     IF (UsingSphericalPolarGrid) THEN
219    
220     DO J=1,sNy
221     DO I=1,sNx
222     I2 = (sNx)*(J-1)+I
223    
224     FMASK1(I2) = 0
225    
226     cc ATL1 experiment: AMERICA - ATLANTIC - AFRICA
227     cc each block is defined within two meridians
228     cc bounded by 30S and 60N (maximum northward extent
229     cc of the model)
230    
231     FMASK1(I2) = 0.0
232     if ( (LON_G(I2) .GE. 105.) .AND. (LON_G(I2) .LE. 150.) .AND. (LAT_G(I2) .GE. -30) ) then
233     FMASK1(I2) = 1.0
234     end if
235     if ( (LON_G(I2) .GE. 210.) .AND. (LON_G(I2) .LE. 255.) .AND. (LAT_G(I2) .GE. -30) ) then
236     FMASK1(I2) = 1.0
237     end if
238    
239    
240     END DO
241     END DO
242    
243    
244     ELSE
245    
246     C Cartesian Grid: nothing is done for now...
247    
248     END IF
249    
250    
251     C 5. QFLUX
252     C
253     C for spherical geometry
254     IF (UsingSphericalPolarGrid) THEN
255    
256     DO J=1,sNy
257     DO I=1,sNx
258     I2 = (sNx)*(J-1)+I
259    
260     cc QSTAR: `equinox', centered on the equator
261     cc with meridional scale set by DELL_STAR and strength
262     cc by DELQ_STAR. This is added to a background cst value QSTAR_BACK
263    
264     QSTAR(I2) = QSTAR_BACK
265    
266     dist = ABS ( LAT_G(I2) )
267    
268     if (dist .LE. DELL_STAR) then
269     QSTAR(I2) = QSTAR(I2) + DELQ_STAR * cos ( pi * dist / 2./DELL_STAR ) **2
270     end if
271    
272     cc
273     cc QOCEAN: MOC and / or ENSO forcing
274     cc
275    
276     if (FMASK1(I2) .EQ. 1) then
277    
278     QOCEAN(I2) = 0.
279    
280     else
281    
282     QOCEAN(I2) = 0.
283    
284     cc MOC: meridional dipole across the equator only over the Atlantic
285    
286     if ( (LON_G(I2) .GE. 150.) .AND. (LON_G(I2) .LE. 210.) ) then
287     dist = ABS ( LAT_G(I2) )
288     if (dist .LE. DELL_OCEAN) then
289     QOCEAN(I2) = DELQ_OCEAN * sin ( pi * LAT_G(I2) / DELL_OCEAN )
290     else
291     QOCEAN(I2) = 0.
292     end if
293     end if
294    
295     cc ENSO: Nino3 monopole in the `tropical Pacific'
296     cc with meridional scale DELL_OCEAN, non zero between longitudes
297     cc LOE_OCEAN and LOW_OCEAN, and amplitude DELQ_OCEAN
298     cc
299     cc !!!!!! check that QOCEAN is not reset when entering here...
300     cc + define a new DELQ_OCEAN for ENSO...
301     cc
302    
303     cc if ( (LON_G(I2) .LE. LOW_OCEAN) .AND. (LON_G(I2) .GE. LOE_OCEAN) ) then
304     cc QOCEAN(I2) = DELQ_OCEAN * sin( pi*(LON_G(I2)-LOE_OCEAN)
305     cc & / (LOW_OCEAN-LOE_OCEAN) )
306     cc end if
307     cc dist = ABS ( LAT_G(I2) )
308     cc if (dist .LE. DELL_OCEAN) then
309     cc QOCEAN(I2) = QOCEAN(I2) * cos ( pi * dist / 2./DELL_OCEAN ) **2
310     cc else
311     cc QOCEAN(I2) = 0.
312     cc end if
313    
314     end if
315    
316     QFLUX(I2) = QSTAR(I2) + QOCEAN(I2)
317    
318     END DO
319     END DO
320    
321     CC(acz)
322     CC CALL WRITE_FLD_XY_RS('QFLUX','0000000000',QFLUX,0,myThid)
323    
324     ELSE
325    
326     C Cartesian Grid: nothing is done for now...
327    
328     END IF
329    
330    
331     C DO J=1,sNy
332     C DO I=1,sNx
333     C I2 = (sNx)*(J-1)+I
334     C SST1 (I2)= SST_REF(I2)
335     C CBMF_DYN(I,J,bi,bj) = CBMFG1(I2)
336     C SST_DYN(I,J,bi,bj) = SST_REF(I2)
337     C END DO
338     C END DO
339     C
340     C _EXCH_XY_R8(CBMF_DYN,myThid)
341     C _EXCH_XY_R8(SST_DYN,myThid)
342    
343    
344     c -- for the MIT physics, this is not necessary anymore:
345     ccc CALL INPHYS( FSG, HSG, RLAT )
346     c --
347    
348    
349     #ifdef ALLOW_TIMEAVE
350     C Initialise diagnostic counters
351     C (these are cleared on model starti i.e. not
352     C loaded from history file for now ).
353     DO bj = myByLo(myThid), myByHi(myThid)
354     DO bi = myBxLo(myThid), myBxHi(myThid)
355     CALL TIMEAVE_RESET(TSWtave, 1, bi, bj, myThid)
356     CALL TIMEAVE_RESET(TLWtave, 1, bi, bj, myThid)
357     CALL TIMEAVE_RESET(SSWtave, 1, bi, bj, myThid)
358     CALL TIMEAVE_RESET(SLWtave, 1, bi, bj, myThid)
359     CALL TIMEAVE_RESET(SINStave, 1, bi, bj, myThid)
360     CALL TIMEAVE_RESET(BKTtave, 1, bi, bj, myThid)
361     CALL TIMEAVE_RESET(SHFtave, 1, bi, bj, myThid)
362     CALL TIMEAVE_RESET(LHFtave, 1, bi, bj, myThid)
363     CALL TIMEAVE_RESET(PRECNVtave, 1, bi, bj, myThid)
364     CALL TIMEAVE_RESET(PRECLStave, 1, bi, bj, myThid)
365     CALL TIMEAVE_RESET(CLOUDCtave, 1, bi, bj, myThid)
366     CALL TIMEAVE_RESET(SSTtave, 1, bi, bj, myThid)
367     CALL TIMEAVE_RESET(PRWtave, 1, bi, bj, myThid)
368     CALL TIMEAVE_RESET(TSW0tave, 1, bi, bj, myThid)
369     CALL TIMEAVE_RESET(TLW0tave, 1, bi, bj, myThid)
370     CALL TIMEAVE_RESET(CBMFtave, 1, bi, bj, myThid)
371    
372     CALL TIMEAVE_RESET(RHtave, Nr, bi, bj, myThid)
373     CALL TIMEAVE_RESET(CLDFtave, Nr, bi, bj, myThid)
374     CALL TIMEAVE_RESET(CLDQtave, Nr, bi, bj, myThid)
375     CALL TIMEAVE_RESET(CLDQCtave, Nr, bi, bj, myThid)
376     CALL TIMEAVE_RESET(RCtave, Nr, bi, bj, myThid)
377     CALL TIMEAVE_RESET(CRLWtave, Nr, bi, bj, myThid)
378     CALL TIMEAVE_RESET(CRSWtave, Nr, bi, bj, myThid)
379     CALL TIMEAVE_RESET(MUPtave, Nr, bi, bj, myThid)
380     CALL TIMEAVE_RESET(MDNtave, Nr, bi, bj, myThid)
381     CALL TIMEAVE_RESET(MDN0tave, Nr, bi, bj, myThid)
382     CALL TIMEAVE_RESET(ENTtave, Nr, bi, bj, myThid)
383     CALL TIMEAVE_RESET(DETtave, Nr, bi, bj, myThid)
384     CALL TIMEAVE_RESET(Utave, Nr, bi, bj, myThid)
385     CALL TIMEAVE_RESET(Vtave, Nr, bi, bj, myThid)
386     CALL TIMEAVE_RESET(Wtave, Nr, bi, bj, myThid)
387     CALL TIMEAVE_RESET(Ttave, Nr, bi, bj, myThid)
388     CALL TIMEAVE_RESET(Qtave, Nr, bi, bj, myThid)
389     CALL TIMEAVE_RESET(PHItave, Nr, bi, bj, myThid)
390     CALL TIMEAVE_RESET(SXtave, Nr, bi, bj, myThid)
391     CALL TIMEAVE_RESET(HXtave, Nr, bi, bj, myThid)
392     CALL TIMEAVE_RESET(DTCNVtave, Nr, bi, bj, myThid)
393     CALL TIMEAVE_RESET(DTLSCtave, Nr, bi, bj, myThid)
394     CALL TIMEAVE_RESET(DTPBLtave, Nr, bi, bj, myThid)
395    
396     DO K = 1, Nr
397     MITPHYS_TimeAve(K,bi,bj) = 0.
398     ENDDO
399    
400     ENDDO
401     ENDDO
402     #endif /* ALLOW_TIMEAVE */
403    
404     #endif /* ALLOW_MITPHYS */
405    
406     RETURN
407     END
408    
409    
410    
411    
412    
413    

  ViewVC Help
Powered by ViewVC 1.1.22