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

Contents 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 - (show 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 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