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 |
|