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