| 3 |
michaesp |
1 |
program ptos
|
|
|
2 |
|
|
|
3 |
C ******************************************************************
|
|
|
4 |
C
|
|
|
5 |
C NAME:
|
|
|
6 |
C p2s
|
|
|
7 |
C
|
|
|
8 |
C Purpose
|
|
|
9 |
C -------
|
|
|
10 |
C
|
|
|
11 |
C Calculates secondary data files from primary data files
|
|
|
12 |
C (based upon IVE-routines).
|
|
|
13 |
C
|
|
|
14 |
C calling:
|
|
|
15 |
C p2s [-m] file variable-list [-s] [-o]
|
|
|
16 |
C
|
|
|
17 |
C example:
|
|
|
18 |
C p2s P911201_00 TH PV CH QX
|
|
|
19 |
C
|
|
|
20 |
C p2s -m #returns man-page
|
|
|
21 |
C
|
|
|
22 |
C Authors
|
|
|
23 |
C ------
|
|
|
24 |
C
|
|
|
25 |
C H. Wernli April 96
|
|
|
26 |
C D.N. Bresch 980311
|
|
|
27 |
C
|
|
|
28 |
C Modifications
|
|
|
29 |
C -------------
|
|
|
30 |
C completely rewritten by D.N. Bresch 980311 for F90
|
|
|
31 |
C and many more variables added (nearly all available in IVE)
|
|
|
32 |
C
|
|
|
33 |
C ADD YOUR OWN VARIABLES AT ALL PLACES WITH (++)
|
|
|
34 |
C for easy simple calculations, see VEL, M, B
|
|
|
35 |
C for complicated calculations, see VORT, QX and PVR
|
|
|
36 |
C
|
|
|
37 |
C Remarks
|
|
|
38 |
C -------
|
|
|
39 |
C ZB is only read once when needed, i.e. for first time...
|
|
|
40 |
C
|
|
|
41 |
C- ******************************************************************
|
|
|
42 |
|
|
|
43 |
integer,parameter :: ntmax=200,nzmax=200,nvarmax=100
|
|
|
44 |
real time(ntmax),time2(ntmax)
|
|
|
45 |
REAL,ALLOCATABLE, DIMENSION (:,:) :: sp,cl,tl,f,zb,t2m,td2m,vip,
|
|
|
46 |
> u10m,v10m,oro,gradpv
|
|
|
47 |
REAL,ALLOCATABLE, DIMENSION (:,:,:) :: var,th,pv,lpv,the,rh,dhr,
|
|
|
48 |
> tt,qq,uu,vv,ww,rho,alpha,zz,mm,zlay,ug,vg,fl,ipv
|
|
|
49 |
character*80 cdfnam,cstnam,outfnam
|
|
|
50 |
integer cdfid,cdfid1,cstid,ierr,ndim,vardim(4),stat
|
|
|
51 |
integer cdfid2,vardim2(4)
|
|
|
52 |
real dx,dy,mdv,varmin(4),varmax(4),stag(4)
|
|
|
53 |
real aklev(nzmax),bklev(nzmax),aklay(nzmax),bklay(nzmax),
|
|
|
54 |
> ak(nzmax),bk(nzmax)
|
|
|
55 |
integer nx,ny,nz,ntimes,ntimes2,i,j,k,n
|
|
|
56 |
integer stdate(5)
|
|
|
57 |
character*(80) qmode,arg,vnam(nvarmax)
|
|
|
58 |
integer mode,zdef
|
|
|
59 |
real rlat,rlon,lat
|
|
|
60 |
real pollon,pollat,yphys
|
|
|
61 |
real phstoph
|
|
|
62 |
logical prelev
|
|
|
63 |
|
|
|
64 |
real,parameter :: pi=3.141592654
|
|
|
65 |
|
|
|
66 |
integer PS_out,TH_out,TH_calc,RH_out,RH_calc
|
|
|
67 |
integer PV_out,PV_calc,THE_out,THE_calc,VIP_out,VIP_calc
|
|
|
68 |
integer GRADPV_out,GRADPV_calc
|
|
|
69 |
integer LPV_out,LPV_calc
|
|
|
70 |
integer CH_out,CH_calc,PVR_out,PVR_calc
|
|
|
71 |
integer THW_out,THW_calc,Z_outP,Z_out,Z_calc
|
|
|
72 |
integer DIVQU_out,DIVQU_calc
|
|
|
73 |
integer NSQ_out,NSQ_calc,RHO_out,RHO_calc,ALPHA_out,ALPHA_calc,VEL_out,VEL_calc
|
|
|
74 |
integer NSQM_out,NSQM_calc,W_out,W_calc,M_out,M_calc
|
|
|
75 |
integer VORT_out,VORT_calc,UG_out,UG_calc,VG_out,VG_calc
|
|
|
76 |
integer AVO_out,AVO_calc,CURVO_out,CURVO_calc
|
|
|
77 |
integer DTHDP_out,DTHDP_calc
|
|
|
78 |
integer COS_out
|
|
|
79 |
integer ZLAY_out,ZLAY_calc,UA_out,UA_calc,VA_out,VA_calc
|
|
|
80 |
integer P_out,P_calc,PLEV_out,PLEV_calc
|
|
|
81 |
integer QXF_out,QXF_calc,QYF_out,QYF_calc
|
|
|
82 |
integer QX_out,QX_calc,QY_out,QY_calc,PSRED_calc,PSRED_out
|
|
|
83 |
integer RI_calc,RI_out,BLH_calc,BLH_out
|
|
|
84 |
integer GRADTH_out,GRADTH_calc,B_out,B_calc !(++)
|
|
|
85 |
logical verbose
|
|
|
86 |
logical ZonP,TonP,PSonP,UonP,VonP,OMEGAonP,ZBonP,ZBneed
|
|
|
87 |
logical T2MonP,TD2MonP,U10MonP,V10MonP,PVonP,PV3onP
|
|
|
88 |
character*(80) zbfile
|
|
|
89 |
|
|
|
90 |
c set defaults:
|
|
|
91 |
verbose=.false.
|
|
|
92 |
ZonP=.false.
|
|
|
93 |
TonP=.false.
|
|
|
94 |
PSonP=.false.
|
|
|
95 |
ZBonP=.false.
|
|
|
96 |
UonP=.false.
|
|
|
97 |
VonP=.false.
|
|
|
98 |
OMEGAonP=.false.
|
|
|
99 |
ZBonP=.false.
|
|
|
100 |
T2MonP=.false.
|
|
|
101 |
TD2MonP=.false.
|
|
|
102 |
U10MonP=.false.
|
|
|
103 |
V10MonP=.false.
|
|
|
104 |
PVonP=.false.
|
|
|
105 |
PV3onP=.false. ! Lukas
|
|
|
106 |
ZBneed=.false.
|
|
|
107 |
zbfile=''
|
|
|
108 |
|
|
|
109 |
qmode='QNone'
|
|
|
110 |
PS_out=1
|
|
|
111 |
TH_out=0
|
|
|
112 |
TH_calc=0
|
|
|
113 |
RH_out=0
|
|
|
114 |
RH_calc=0
|
|
|
115 |
PV_out=0
|
|
|
116 |
PV_calc=0
|
|
|
117 |
LPV_out=0
|
|
|
118 |
LPV_calc=0
|
|
|
119 |
VIP_out=0
|
|
|
120 |
VIP_calc=0
|
|
|
121 |
THE_out=0
|
|
|
122 |
THE_calc=0
|
|
|
123 |
CH_out=0
|
|
|
124 |
CH_calc=0
|
|
|
125 |
PVR_out=0
|
|
|
126 |
PVR_calc=0
|
|
|
127 |
THW_out=0
|
|
|
128 |
THW_calc=0
|
|
|
129 |
DIVQU_out=0
|
|
|
130 |
DIVQU_calc=0
|
|
|
131 |
Z_calc=0
|
|
|
132 |
Z_out=0
|
|
|
133 |
Z_outP=0
|
|
|
134 |
NSQ_out=0
|
|
|
135 |
NSQ_calc=0
|
|
|
136 |
DTHDP_out=0
|
|
|
137 |
DTHDP_calc=0
|
|
|
138 |
NSQM_out=0
|
|
|
139 |
NSQM_calc=0
|
|
|
140 |
RHO_out=0
|
|
|
141 |
RHO_calc=0
|
|
|
142 |
ALPHA_out=0
|
|
|
143 |
ALPHA_calc=0
|
|
|
144 |
VEL_out=0
|
|
|
145 |
VEL_calc=0
|
|
|
146 |
M_out=0
|
|
|
147 |
M_calc=0
|
|
|
148 |
B_out=0
|
|
|
149 |
B_calc=0
|
|
|
150 |
W_out=0
|
|
|
151 |
W_calc=0
|
|
|
152 |
VORT_out=0
|
|
|
153 |
VORT_calc=0
|
|
|
154 |
AVO_out=0
|
|
|
155 |
AVO_calc=0
|
|
|
156 |
CURVO_out=0
|
|
|
157 |
CURVO_calc=0
|
|
|
158 |
COS_out=0
|
|
|
159 |
UG_out=0
|
|
|
160 |
UG_calc=0
|
|
|
161 |
VG_out=0
|
|
|
162 |
VG_calc=0
|
|
|
163 |
UA_out=0
|
|
|
164 |
UA_calc=0
|
|
|
165 |
VA_out=0
|
|
|
166 |
VA_calc=0
|
|
|
167 |
ZLAY_out=0
|
|
|
168 |
ZLAY_calc=0
|
|
|
169 |
P_out=0
|
|
|
170 |
P_calc=0
|
|
|
171 |
PLEV_out=0
|
|
|
172 |
PLEV_calc=0
|
|
|
173 |
FL_out=0
|
|
|
174 |
FL_calc=0
|
|
|
175 |
PSRED_out=0
|
|
|
176 |
PSRED_calc=0
|
|
|
177 |
RI_out=0
|
|
|
178 |
RI_calc=0
|
|
|
179 |
BLH_out=0
|
|
|
180 |
BLH_calc=0
|
|
|
181 |
QXF_out=0
|
|
|
182 |
QXF_calc=0
|
|
|
183 |
QYF_out=0
|
|
|
184 |
QYF_calc=0
|
|
|
185 |
QX_out=0
|
|
|
186 |
QX_calc=0
|
|
|
187 |
QY_out=0
|
|
|
188 |
QY_calc=0
|
|
|
189 |
GRADTH_out=0
|
|
|
190 |
GRADTH_calc=0
|
|
|
191 |
GRADPV_out=0
|
|
|
192 |
GRADPV_calc=0 !(++)
|
|
|
193 |
|
|
|
194 |
c get arguments:
|
|
|
195 |
c get parameters from command-line:
|
|
|
196 |
c COUNT THE ARGUMENTS:
|
|
|
197 |
if (iargc() .lt. 1) then
|
|
|
198 |
print*,'USAGE: p2s [-m] file variable-list [-s] ',
|
|
|
199 |
> '[-o] [-zb file]'
|
|
|
200 |
STOP
|
|
|
201 |
endif
|
|
|
202 |
|
|
|
203 |
c REQUESTD INPUT:
|
|
|
204 |
c ---------------
|
|
|
205 |
c GET WITH getarg DIRECTLY FROM SHELL:
|
|
|
206 |
call getarg(1,cdfnam)
|
|
|
207 |
if (trim(cdfnam).eq.'-m') then
|
|
|
208 |
print*,' '
|
|
|
209 |
print*,'computes derived variables from primary ones on'
|
|
|
210 |
print*,'the input-file'
|
|
|
211 |
print*,'if the output-file is already present, it will be'
|
|
|
212 |
print*,'updated'
|
|
|
213 |
print*,' '
|
|
|
214 |
print*,'check the source-code itself about details...'
|
|
|
215 |
print*,' '
|
|
|
216 |
print*,'file: netCDF file with basic variables on it'
|
|
|
217 |
print*,' requested are the variables needed to calculate'
|
|
|
218 |
print*,' the requested output (variable-list)'
|
|
|
219 |
print*,' If file is a P-file (starting with P), the output-'
|
|
|
220 |
print*,' file will be an S-file, otherwise the extension'
|
|
|
221 |
print*,' _out will be appended, unless -o is used'
|
|
|
222 |
print*,' if the S-file exists already, it is tried to '
|
|
|
223 |
print*,' append the new variable'
|
|
|
224 |
print*,'[P]date means that you can either give PYYMMDD_HH'
|
|
|
225 |
print*,' or YYMMDD_HH alone'
|
|
|
226 |
print*,'variable-list: a list of variables to be calculated'
|
|
|
227 |
print*,' and written to the S_file, available are:'
|
|
|
228 |
print*,' TH,PV,LPV,RH,THE,THW,CH,PVR,Z,ZonP,GRADTH,NSQ,NSQM'
|
|
|
229 |
print*,' M,B,W,RHO,VEL,VORT,AVO,CURVO,UG,VG,ZLAY,UA,VA,P'
|
|
|
230 |
print*,' PLEV,FL,QX,QY,QXF,QYF,PSRED,RI,BLH,GRADPV,DTHDP'
|
|
|
231 |
print*,' ALPHA'
|
|
|
232 |
print*,'-s: only small S-file, i.e. TH and PV'
|
|
|
233 |
print*,'-zb: file with ZB (for PSRED)'
|
|
|
234 |
print*,'-o output: filename of the output netCDF file'
|
|
|
235 |
STOP
|
|
|
236 |
endif
|
|
|
237 |
C check, if cdfnam is with or without P:
|
|
|
238 |
if (cdfnam(1:1).eq.'P') then
|
|
|
239 |
outfnam='S'//cdfnam(2:len_trim(cdfnam))
|
|
|
240 |
else
|
|
|
241 |
outfnam=trim(cdfnam)//'_out'
|
|
|
242 |
endif
|
|
|
243 |
i=2
|
|
|
244 |
do while (iargc().ge.i)
|
|
|
245 |
call getarg(i,arg)
|
|
|
246 |
i=i+1
|
|
|
247 |
if (arg.eq.'TH') TH_out=1
|
|
|
248 |
if (arg.eq.'THE') THE_out=1
|
|
|
249 |
if (arg.eq.'THW') THW_out=1
|
|
|
250 |
if (arg.eq.'DIVQU') DIVQU_out=1
|
|
|
251 |
if (arg.eq.'PV') PV_out=1
|
|
|
252 |
if (arg.eq.'LPV') LPV_out=1
|
|
|
253 |
if (arg.eq.'VIP') VIP_out=1
|
|
|
254 |
if (arg.eq.'RH') RH_out=1
|
|
|
255 |
if (arg.eq.'CH') CH_out=1
|
|
|
256 |
if (arg.eq.'PVR') PVR_out=1
|
|
|
257 |
if (arg.eq.'Z') Z_out=1
|
|
|
258 |
if (arg.eq.'ZonP') Z_outP=1
|
|
|
259 |
if (arg.eq.'GRADTH') GRADTH_out=1
|
|
|
260 |
if (arg.eq.'GRADPV') GRADPV_out=1
|
|
|
261 |
if (arg.eq.'VEL') VEL_out=1
|
|
|
262 |
if (arg.eq.'RHO') RHO_out=1
|
|
|
263 |
if (arg.eq.'ALPHA') ALPHA_out=1
|
|
|
264 |
if (arg.eq.'W') W_out=1
|
|
|
265 |
if (arg.eq.'M') M_out=1
|
|
|
266 |
if (arg.eq.'B') B_out=1
|
|
|
267 |
if (arg.eq.'VORT') VORT_out=1
|
|
|
268 |
if (arg.eq.'AVO') AVO_out=1
|
|
|
269 |
if (arg.eq.'CURVO') CURVO_out=1
|
|
|
270 |
if (arg.eq.'COS') COS_out=1
|
|
|
271 |
if (arg.eq.'UG') UG_out=1
|
|
|
272 |
if (arg.eq.'VG') VG_out=1
|
|
|
273 |
if (arg.eq.'UA') UA_out=1
|
|
|
274 |
if (arg.eq.'VA') VA_out=1
|
|
|
275 |
if (arg.eq.'QX') QX_out=1
|
|
|
276 |
if (arg.eq.'QY') QY_out=1
|
|
|
277 |
if (arg.eq.'QXF') QXF_out=1
|
|
|
278 |
if (arg.eq.'QYF') QYF_out=1
|
|
|
279 |
if (arg.eq.'ZLAY') ZLAY_out=1
|
|
|
280 |
if (arg.eq.'P') P_out=1
|
|
|
281 |
if (arg.eq.'PLEV') PLEV_out=1
|
|
|
282 |
if (arg.eq.'FL') FL_out=1
|
|
|
283 |
if (arg.eq.'PSRED') PSRED_out=1
|
|
|
284 |
if (arg.eq.'RI') RI_out=1
|
|
|
285 |
if (arg.eq.'BLH') BLH_out=1
|
|
|
286 |
if (arg.eq.'NSQ') NSQ_out=1
|
|
|
287 |
if (arg.eq.'DTHDP') DTHDP_out=1
|
|
|
288 |
if (arg.eq.'NSQM') NSQM_out=1 ! (++)
|
|
|
289 |
if (arg.eq.'-v') verbose=.true.
|
|
|
290 |
if (arg.eq.'-s') then
|
|
|
291 |
TH_out=1
|
|
|
292 |
PV_out=1
|
|
|
293 |
endif
|
|
|
294 |
if (arg.eq.'-o') then
|
|
|
295 |
if (iargc().ge.i) then
|
|
|
296 |
call getarg(i,outfnam)
|
|
|
297 |
i=i+1
|
|
|
298 |
else
|
|
|
299 |
print*,'option -o requires filename'
|
|
|
300 |
STOP
|
|
|
301 |
endif
|
|
|
302 |
endif
|
|
|
303 |
if (arg.eq.'-zb') then
|
|
|
304 |
if (iargc().ge.i) then
|
|
|
305 |
call getarg(i,zbfile)
|
|
|
306 |
ZBneed=.true.
|
|
|
307 |
i=i+1
|
|
|
308 |
else
|
|
|
309 |
print*,'option -zb requires filename'
|
|
|
310 |
STOP
|
|
|
311 |
endif
|
|
|
312 |
endif
|
|
|
313 |
if (arg.eq.'-nops') then
|
|
|
314 |
PS_out=0
|
|
|
315 |
endif
|
|
|
316 |
enddo
|
|
|
317 |
|
|
|
318 |
C force calculations of requested fields:
|
|
|
319 |
c ---------------------------------------
|
|
|
320 |
if (QX_out.eq.1) QX_calc=1
|
|
|
321 |
if (QY_out.eq.1) QY_calc=1
|
|
|
322 |
if (QXF_out.eq.1) QXF_calc=1
|
|
|
323 |
if (QYF_out.eq.1) QYF_calc=1
|
|
|
324 |
if (ZLAY_out.eq.1) ZLAY_calc=1
|
|
|
325 |
if (P_out.eq.1) P_calc=1
|
|
|
326 |
if (PLEV_out.eq.1) PLEV_calc=1
|
|
|
327 |
if (FL_out.eq.1) FL_calc=1
|
|
|
328 |
if (UG_out.eq.1) UG_calc=1
|
|
|
329 |
if (VG_out.eq.1) VG_calc=1
|
|
|
330 |
if (UA_out.eq.1) UA_calc=1
|
|
|
331 |
if (VA_out.eq.1) VA_calc=1
|
|
|
332 |
if (VORT_out.eq.1) VORT_calc=1
|
|
|
333 |
if (AVO_out.eq.1) AVO_calc=1
|
|
|
334 |
if (CURVO_out.eq.1) CURVO_calc=1
|
|
|
335 |
if (TH_out.eq.1) TH_calc=1
|
|
|
336 |
if (PV_out.eq.1) PV_calc=1
|
|
|
337 |
if (LPV_out.eq.1) LPV_calc=1
|
|
|
338 |
if (VIP_out.eq.1) VIP_calc=1
|
|
|
339 |
if (RH_out.eq.1) RH_calc=1
|
|
|
340 |
if (THE_out.eq.1) THE_calc=1
|
|
|
341 |
if (THW_out.eq.1) THW_calc=1
|
|
|
342 |
if (DIVQU_out.eq.1) DIVQU_calc=1
|
|
|
343 |
if (CH_out.eq.1) CH_calc=1
|
|
|
344 |
if (PVR_out.eq.1) PVR_calc=1
|
|
|
345 |
if (Z_out.eq.1) Z_calc=1
|
|
|
346 |
if (Z_outP.eq.1) Z_calc=1
|
|
|
347 |
if (GRADTH_out.eq.1) GRADTH_calc=1
|
|
|
348 |
if (GRADPV_out.eq.1) GRADPV_calc=1
|
|
|
349 |
if (RHO_out.eq.1) RHO_calc=1
|
|
|
350 |
if (ALPHA_out.eq.1) ALPHA_calc=1
|
|
|
351 |
if (RI_out.eq.1) RI_calc=1
|
|
|
352 |
if (BLH_out.eq.1) BLH_calc=1
|
|
|
353 |
if (VEL_out.eq.1) VEL_calc=1
|
|
|
354 |
if (M_out.eq.1) M_calc=1
|
|
|
355 |
if (B_out.eq.1) B_calc=1
|
|
|
356 |
if (NSQM_out.eq.1) NSQM_calc=1
|
|
|
357 |
if (W_out.eq.1) W_calc=1
|
|
|
358 |
if (DTHDP_out.eq.1) DTHDP_calc=1
|
|
|
359 |
if (NSQ_out.eq.1) NSQ_calc=1 !(++)
|
|
|
360 |
|
|
|
361 |
C make dependencies for variable calculations:
|
|
|
362 |
c --------------------------------------------
|
|
|
363 |
if (NSQ_calc.eq.1) TH_calc=1 !(++)
|
|
|
364 |
if (DTHDP_calc.eq.1) TH_calc=1
|
|
|
365 |
if (PSRED_out.eq.1) PSRED_calc=1
|
|
|
366 |
if (PSRED_calc.eq.1) ZBneed=.true.
|
|
|
367 |
if (QX_out.eq.1) UG_calc=1
|
|
|
368 |
if (QY_out.eq.1) VG_calc=1
|
|
|
369 |
if (UA_calc.eq.1) UG_calc=1
|
|
|
370 |
if (VA_calc.eq.1) VG_calc=1
|
|
|
371 |
if (UG_calc.eq.1) ZLAY_calc=1
|
|
|
372 |
if (VG_calc.eq.1) ZLAY_calc=1
|
|
|
373 |
if (ZLAY_calc.eq.1) Z_calc=1
|
|
|
374 |
if (B_calc.eq.1) M_calc=1
|
|
|
375 |
if (M_calc.eq.1) ZLAY_calc=1
|
|
|
376 |
if (NSQ_calc.eq.1) RHO_calc=1
|
|
|
377 |
if (NSQM_calc.eq.1) RHO_calc=1
|
|
|
378 |
if (NSQM_calc.eq.1) THE_calc=1
|
|
|
379 |
if (W_calc.eq.1) RHO_calc=1
|
|
|
380 |
if (GRADTH_calc.eq.1) TH_calc=1
|
|
|
381 |
if (THW_calc.eq.1) THE_calc=1
|
|
|
382 |
if (THE_calc.eq.1) TH_calc=1
|
|
|
383 |
if (VIP_calc.eq.1) PV_calc=1
|
|
|
384 |
if (PV_calc.eq.1) TH_calc=1
|
|
|
385 |
if (LPV_calc.eq.1) PV_calc=1
|
|
|
386 |
if (PVR_calc.eq.1) CH_calc=1
|
|
|
387 |
if (CH_calc.eq.1) TH_calc=1
|
|
|
388 |
if (CH_calc.eq.1) RH_calc=1
|
|
|
389 |
if (RI_calc.eq.1) TH_calc=1
|
|
|
390 |
if (ALPHA_calc.eq.1) RHO_calc=1
|
|
|
391 |
|
|
|
392 |
print*,'processing: ',trim(cdfnam)
|
|
|
393 |
|
|
|
394 |
C Open files and get infos about data domain
|
|
|
395 |
c ------------------------------------------
|
|
|
396 |
|
|
|
397 |
if (Z_outP.eq.1) then
|
|
|
398 |
call cdfwopn(trim(cdfnam),cdfid1,ierr)
|
|
|
399 |
else
|
|
|
400 |
call cdfopn(trim(cdfnam),cdfid1,ierr)
|
|
|
401 |
endif
|
|
|
402 |
if (ierr.ne.0) then
|
|
|
403 |
print*,'ERROR opening input file, stopped'
|
|
|
404 |
stop
|
|
|
405 |
endif
|
|
|
406 |
call getcfn(cdfid1,cstnam,ierr)
|
|
|
407 |
if (ierr /= 0) then
|
|
|
408 |
cstnam = cdfnam
|
|
|
409 |
cstid = cdfid1
|
|
|
410 |
else
|
|
|
411 |
call cdfopn(trim(cstnam),cstid,ierr)
|
|
|
412 |
endif
|
|
|
413 |
|
|
|
414 |
|
|
|
415 |
C Inquire the variables on the netCDF file:
|
|
|
416 |
c -----------------------------------------
|
|
|
417 |
C Inquire number of variables and variable names
|
|
|
418 |
call getvars(cdfid1,nvars,vnam,ierr)
|
|
|
419 |
C print*,nvars,vnam(1),vnam(2)
|
|
|
420 |
do i=1,nvars
|
|
|
421 |
if (trim(vnam(i)).eq.'Q') qmode='Q'
|
|
|
422 |
if (trim(vnam(i)).eq.'QD') qmode='QD'
|
|
|
423 |
if (trim(vnam(i)).eq.'Z') ZonP=.true.
|
|
|
424 |
if (trim(vnam(i)).eq.'T') TonP=.true.
|
|
|
425 |
if (trim(vnam(i)).eq.'PS') PSonP=.true.
|
|
|
426 |
if (trim(vnam(i)).eq.'U') UonP=.true.
|
|
|
427 |
if (trim(vnam(i)).eq.'V') VonP=.true.
|
|
|
428 |
if (trim(vnam(i)).eq.'ZB') ZBonP=.true.
|
|
|
429 |
if (trim(vnam(i)).eq.'T2M') T2MonP=.true.
|
|
|
430 |
if (trim(vnam(i)).eq.'TD2M') TD2MonP=.true.
|
|
|
431 |
if (trim(vnam(i)).eq.'U10M') U10MonP=.true.
|
|
|
432 |
if (trim(vnam(i)).eq.'V10M') V10MonP=.true.
|
|
|
433 |
if (trim(vnam(i)).eq.'OMEGA') OMEGAonP=.true.
|
|
|
434 |
if (trim(vnam(i)).eq.'PV') PVonP=.true.
|
|
|
435 |
if (trim(vnam(i)).eq.'PV3') PV3onP=.true.
|
|
|
436 |
enddo
|
|
|
437 |
|
|
|
438 |
|
|
|
439 |
if (TonP) then
|
|
|
440 |
call getdef(cdfid1,'T',ndim,mdv,vardim,varmin,varmax,stag,ierr)
|
|
|
441 |
else
|
|
|
442 |
call getdef(cdfid1,vnam(2),ndim,mdv,
|
|
|
443 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
444 |
endif
|
|
|
445 |
if (ierr.ne.0) goto 920
|
|
|
446 |
mdv=-999.98999
|
|
|
447 |
|
|
|
448 |
C Get the levels, pole, etc.
|
|
|
449 |
c --------------------------
|
|
|
450 |
|
|
|
451 |
nx=vardim(1)
|
|
|
452 |
ny=vardim(2)
|
|
|
453 |
nz=vardim(3)
|
|
|
454 |
|
|
|
455 |
C print*,nx,ny,nz
|
|
|
456 |
call getgrid(cstid,dx,dy,ierr)
|
|
|
457 |
call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
|
|
|
458 |
call getpole(cstid,pollon,pollat,ierr)
|
|
|
459 |
call getstart(cstid,stdate,ierr)
|
|
|
460 |
|
|
|
461 |
C Allocate all arrays:
|
|
|
462 |
C --------------------
|
|
|
463 |
allocate(sp(nx,ny),STAT=stat)
|
|
|
464 |
if (stat.ne.0) print*,'error allocating sp(nx,ny)'
|
|
|
465 |
allocate(oro(nx,ny),STAT=stat)
|
|
|
466 |
if (stat.ne.0) print*,'error allocating oro(nx,ny)'
|
|
|
467 |
allocate(cl(nx,ny),STAT=stat)
|
|
|
468 |
if (stat.ne.0) print*,'error allocating cl(nx,ny)'
|
|
|
469 |
allocate(tl(nx,ny),STAT=stat)
|
|
|
470 |
if (stat.ne.0) print*,'error allocating tl(nx,ny)'
|
|
|
471 |
allocate(f(nx,ny),STAT=stat)
|
|
|
472 |
if (stat.ne.0) print*,'error allocating f(nx,ny)'
|
|
|
473 |
|
|
|
474 |
c allocate memory for results
|
|
|
475 |
allocate(var(nx,ny,nz),STAT=stat)
|
|
|
476 |
if (stat.ne.0) print*,'error allocating var(nx,ny,nz)'
|
|
|
477 |
|
|
|
478 |
c allocate memory for variables on P-file (allocate all)
|
|
|
479 |
allocate(tt(nx,ny,nz),STAT=stat)
|
|
|
480 |
if (stat.ne.0) print*,'error allocating tt(nx,ny,nz)'
|
|
|
481 |
allocate(qq(nx,ny,nz),STAT=stat)
|
|
|
482 |
if (stat.ne.0) print*,'error allocating qq(nx,ny,nz)'
|
|
|
483 |
allocate(uu(nx,ny,nz),STAT=stat)
|
|
|
484 |
if (stat.ne.0) print*,'error allocating uu(nx,ny,nz)'
|
|
|
485 |
allocate(vv(nx,ny,nz),STAT=stat)
|
|
|
486 |
if (stat.ne.0) print*,'error allocating vv(nx,ny,nz)'
|
|
|
487 |
allocate(ww(nx,ny,nz),STAT=stat)
|
|
|
488 |
if (stat.ne.0) print*,'error allocating ww(nx,ny,nz)'
|
|
|
489 |
|
|
|
490 |
allocate(t2m(nx,ny),STAT=stat)
|
|
|
491 |
if (stat.ne.0) print*,'error allocating t2m(nx,ny)'
|
|
|
492 |
allocate(td2m(nx,ny),STAT=stat)
|
|
|
493 |
if (stat.ne.0) print*,'error allocating td2m(nx,ny)'
|
|
|
494 |
allocate(u10m(nx,ny),STAT=stat)
|
|
|
495 |
if (stat.ne.0) print*,'error allocating u10m(nx,ny)'
|
|
|
496 |
allocate(v10m(nx,ny),STAT=stat)
|
|
|
497 |
if (stat.ne.0) print*,'error allocating v10m(nx,ny)'
|
|
|
498 |
allocate(ipv(nx,ny,nz),STAT=stat)
|
|
|
499 |
if (stat.ne.0) print*,'error allocating ipv(nx,ny,nz)'
|
|
|
500 |
|
|
|
501 |
C Determine if data is on pressure or model levels
|
|
|
502 |
|
|
|
503 |
prelev=.true.
|
|
|
504 |
do k=1,nz
|
|
|
505 |
if (bklev(k).ne.0.) prelev=.false.
|
|
|
506 |
if (bklay(k).ne.0.) prelev=.false.
|
|
|
507 |
enddo
|
|
|
508 |
C print*,'prelev ',prelev
|
|
|
509 |
|
|
|
510 |
C Calculate cos(latitude) array and the coriolis parameter
|
|
|
511 |
|
|
|
512 |
if ((abs(pollon).gt.0.001).or.(abs(pollat-90.).gt.0.001)) then
|
|
|
513 |
do j=1,ny
|
|
|
514 |
rlat=varmin(2)+(j-1)*dy
|
|
|
515 |
do i=1,nx
|
|
|
516 |
rlon=varmin(1)+(i-1)*dx
|
|
|
517 |
yphys=phstoph(rlat,rlon,pollat,pollon)
|
|
|
518 |
C if I use sind(lat in deg): troubles at the N-pole
|
|
|
519 |
lat=2.*pi*yphys/360.
|
|
|
520 |
cl(i,j)=cosd(rlat)
|
|
|
521 |
tl(i,j)=tan(lat)
|
|
|
522 |
f(i,j)=0.000145444*sin(lat)
|
|
|
523 |
enddo
|
|
|
524 |
enddo
|
|
|
525 |
else
|
|
|
526 |
do j=1,ny
|
|
|
527 |
lat=varmin(2)+(j-1)*dy
|
|
|
528 |
lat=2.*pi*lat/360.
|
|
|
529 |
do i=1,nx
|
|
|
530 |
cl(i,j)=cos(lat)
|
|
|
531 |
f(i,j)=0.000145444*sin(lat)
|
|
|
532 |
enddo
|
|
|
533 |
enddo
|
|
|
534 |
endif
|
|
|
535 |
|
|
|
536 |
C Determine if data is on levels or layers
|
|
|
537 |
|
|
|
538 |
if (stag(3).eq.-0.5) then
|
|
|
539 |
ak=aklay
|
|
|
540 |
bk=bklay
|
|
|
541 |
else
|
|
|
542 |
ak=aklev
|
|
|
543 |
bk=bklev
|
|
|
544 |
endif
|
|
|
545 |
|
|
|
546 |
C Get all the fields
|
|
|
547 |
C ------------------
|
|
|
548 |
|
|
|
549 |
call gettimes(cdfid1,time,ntimes,ierr)
|
|
|
550 |
|
|
|
551 |
C Loop over all times
|
|
|
552 |
C -------------------
|
|
|
553 |
|
|
|
554 |
C print*,'loop over all times'
|
|
|
555 |
C print*,ntimes,vardim,nx,ny,nz
|
|
|
556 |
|
|
|
557 |
do n=1,ntimes
|
|
|
558 |
|
|
|
559 |
if (verbose) print*,'time=',time(n)
|
|
|
560 |
|
|
|
561 |
if (.not.prelev) then
|
|
|
562 |
if (PSonP) then
|
|
|
563 |
call getdat(cdfid1,'PS',time(n),0,sp,ierr)
|
|
|
564 |
if (ierr.ne.0) goto 921
|
|
|
565 |
else
|
|
|
566 |
if (verbose) print*,'PS not on P-file'
|
|
|
567 |
endif
|
|
|
568 |
endif
|
|
|
569 |
|
|
|
570 |
if (ZBonP) then
|
|
|
571 |
call getdat(cdfid1,'ZB',time(n),0,oro,ierr)
|
|
|
572 |
if (ierr.ne.0) goto 930
|
|
|
573 |
else
|
|
|
574 |
if (verbose) print*,'ZB not on P-file'
|
|
|
575 |
endif
|
|
|
576 |
|
|
|
577 |
if (TonP) then
|
|
|
578 |
call getdat(cdfid1,'T',time(n),0,tt,ierr)
|
|
|
579 |
if (ierr.ne.0) goto 920
|
|
|
580 |
else
|
|
|
581 |
if (verbose) print*,'T not on P-file'
|
|
|
582 |
endif
|
|
|
583 |
|
|
|
584 |
if (qmode.eq.'Q') then
|
|
|
585 |
call getdat(cdfid1,'Q',time(n),0,qq,ierr)
|
|
|
586 |
if (ierr.ne.0) goto 922
|
|
|
587 |
elseif (qmode.eq.'QD') then
|
|
|
588 |
call getdat(cdfid1,'QD',time(n),0,qq,ierr)
|
|
|
589 |
if (ierr.ne.0) goto 922
|
|
|
590 |
else
|
|
|
591 |
if (verbose) print*,'neither Q nor QD on P-file'
|
|
|
592 |
endif
|
|
|
593 |
|
|
|
594 |
if (ZonP) then
|
|
|
595 |
allocate(zz(nx,ny,nz),STAT=stat)
|
|
|
596 |
if (stat.ne.0) print*,'error allocating zz'
|
|
|
597 |
call getdat(cdfid1,'Z',time(n),0,zz,ierr)
|
|
|
598 |
if (ierr.ne.0) then
|
|
|
599 |
print*,'error reading Z'
|
|
|
600 |
STOP
|
|
|
601 |
else
|
|
|
602 |
Z_calc=0 !indicate Z read
|
|
|
603 |
endif
|
|
|
604 |
endif
|
|
|
605 |
|
|
|
606 |
if (ZBneed) then
|
|
|
607 |
ZBneed=.false.
|
|
|
608 |
allocate(zb(nx,ny),STAT=stat)
|
|
|
609 |
if (stat.ne.0) print*,'error allocating zb'
|
|
|
610 |
if (ZBfile.ne.'') then
|
|
|
611 |
call cdfopn(trim(ZBfile),cdfid2,ierr)
|
|
|
612 |
if (ierr.ne.0) then
|
|
|
613 |
print*,'error opening ',trim(ZBfile)
|
|
|
614 |
ZBneed=.true.
|
|
|
615 |
endif
|
|
|
616 |
call gettimes(cdfid2,time2,ntimes2,ierr)
|
|
|
617 |
call getdat(cdfid2,'ZB',time2(1),0,zb,ierr)
|
|
|
618 |
else if (ZBonP) then
|
|
|
619 |
call getdat(cdfid1,'ZB',time(1),0,zb,ierr)
|
|
|
620 |
else
|
|
|
621 |
print*,'ZB needed (use option -zb filename)'
|
|
|
622 |
ZBneed=.true.
|
|
|
623 |
endif
|
|
|
624 |
if (ierr.ne.0) then
|
|
|
625 |
print*,'error reading ZB'
|
|
|
626 |
ZBneed=.true.
|
|
|
627 |
endif
|
|
|
628 |
endif
|
|
|
629 |
|
|
|
630 |
if (UonP) then
|
|
|
631 |
call getdat(cdfid1,'U',time(n),0,uu,ierr)
|
|
|
632 |
if (ierr.ne.0) goto 923
|
|
|
633 |
else
|
|
|
634 |
if (verbose) print*,'U not on P-file'
|
|
|
635 |
endif
|
|
|
636 |
|
|
|
637 |
if (VonP) then
|
|
|
638 |
call getdat(cdfid1,'V',time(n),0,vv,ierr)
|
|
|
639 |
if (ierr.ne.0) goto 924
|
|
|
640 |
else
|
|
|
641 |
if (verbose) print*,'V not on P-file'
|
|
|
642 |
endif
|
|
|
643 |
|
|
|
644 |
if (OMEGAonP) then
|
|
|
645 |
call getdat(cdfid1,'OMEGA',time(n),0,ww,ierr)
|
|
|
646 |
if (ierr.ne.0) goto 925
|
|
|
647 |
else
|
|
|
648 |
if (verbose) print*,'OMEGA not on P-file'
|
|
|
649 |
endif
|
|
|
650 |
|
|
|
651 |
if (T2MonP) then
|
|
|
652 |
call getdat(cdfid1,'T2M',time(n),0,t2m,ierr)
|
|
|
653 |
if (ierr.ne.0) goto 926
|
|
|
654 |
else
|
|
|
655 |
if (verbose) print*,'T2M not on P-file'
|
|
|
656 |
endif
|
|
|
657 |
|
|
|
658 |
if (TD2MonP) then
|
|
|
659 |
call getdat(cdfid1,'TD2M',time(n),0,td2m,ierr)
|
|
|
660 |
if (ierr.ne.0) goto 926
|
|
|
661 |
else
|
|
|
662 |
if (verbose) print*,'TD2M not on P-file'
|
|
|
663 |
endif
|
|
|
664 |
|
|
|
665 |
if (U10MonP) then
|
|
|
666 |
call getdat(cdfid1,'U10M',time(n),0,u10m,ierr)
|
|
|
667 |
if (ierr.ne.0) goto 926
|
|
|
668 |
else
|
|
|
669 |
if (verbose) print*,'U10M not on P-file'
|
|
|
670 |
endif
|
|
|
671 |
|
|
|
672 |
if (V10MonP) then
|
|
|
673 |
call getdat(cdfid1,'V10M',time(n),0,v10m,ierr)
|
|
|
674 |
if (ierr.ne.0) goto 926
|
|
|
675 |
else
|
|
|
676 |
if (verbose) print*,'V10M not on P-file'
|
|
|
677 |
endif
|
|
|
678 |
|
|
|
679 |
if (PVonP) then
|
|
|
680 |
call getdat(cdfid1,'PV',time(n),0,ipv,ierr)
|
|
|
681 |
if (ierr.ne.0) goto 926
|
|
|
682 |
else
|
|
|
683 |
if (verbose) print*,'PV not on P-file'
|
|
|
684 |
endif
|
|
|
685 |
|
|
|
686 |
if (PV3onP) then
|
|
|
687 |
call getdat(cdfid1,'PV3',time(n),0,ipv,ierr)
|
|
|
688 |
if (ierr.ne.0) goto 926
|
|
|
689 |
else
|
|
|
690 |
if (verbose) print*,'PV3 not on P-file'
|
|
|
691 |
endif
|
|
|
692 |
|
|
|
693 |
C Calculation of the geopotential
|
|
|
694 |
|
|
|
695 |
if (Z_calc.eq.1) then
|
|
|
696 |
allocate(zz(nx,ny,nz),STAT=stat)
|
|
|
697 |
if (stat.ne.0) print*,'error allocating zz'
|
|
|
698 |
call zlev(zz,tt,qq,oro,sp,nx,ny,nz,aklev,bklev)
|
|
|
699 |
endif
|
|
|
700 |
|
|
|
701 |
if (Z_outP.eq.1) then
|
|
|
702 |
if (n.eq.1) then
|
|
|
703 |
stag(3)=0.0
|
|
|
704 |
call putdef(cdfid1,'Z',4,mdv,
|
|
|
705 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
706 |
write(*,*)'*** variable Z created on P-file'
|
|
|
707 |
stag(3)=-0.5
|
|
|
708 |
endif
|
|
|
709 |
call putdat(cdfid1,'Z',time(n),0,zz,ierr)
|
|
|
710 |
endif
|
|
|
711 |
|
|
|
712 |
C Create the secondary data file
|
|
|
713 |
|
|
|
714 |
if (n.eq.1) then
|
|
|
715 |
call cdfwopn(trim(outfnam),cdfid,ierr)
|
|
|
716 |
if (ierr.ne.0) then
|
|
|
717 |
call crecdf(trim(outfnam),cdfid,varmin,varmax,
|
|
|
718 |
> 3,trim(cstnam),ierr)
|
|
|
719 |
write(*,*)'*** NetCDF file ',trim(outfnam),' created'
|
|
|
720 |
else
|
|
|
721 |
write(*,*)'*** NetCDF file ',trim(outfnam),
|
|
|
722 |
> ' used for appending'
|
|
|
723 |
PS_out=0
|
|
|
724 |
endif
|
|
|
725 |
endif
|
|
|
726 |
|
|
|
727 |
|
|
|
728 |
C Put surface pressure on S-file
|
|
|
729 |
|
|
|
730 |
if ((.not.prelev).and.(PS_out.eq.1)) then
|
|
|
731 |
vardim(3)=1
|
|
|
732 |
if (n.eq.1) then
|
|
|
733 |
call putdef(cdfid,'PS',4,mdv,vardim,varmin,varmax,stag,ierr)
|
|
|
734 |
write(*,*)'*** variable PS created on ',trim(outfnam)
|
|
|
735 |
endif
|
|
|
736 |
call putdat(cdfid,'PS',time(n),0,sp,ierr)
|
|
|
737 |
vardim(3)=nz
|
|
|
738 |
endif
|
|
|
739 |
|
|
|
740 |
C Put geopotential on S-file
|
|
|
741 |
|
|
|
742 |
if (Z_out.eq.1) then
|
|
|
743 |
if (n.eq.1) then
|
|
|
744 |
stag(3)=0.0
|
|
|
745 |
call putdef(cdfid,'Z',4,mdv,
|
|
|
746 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
747 |
write(*,*)'*** variable Z created on ',trim(outfnam)
|
|
|
748 |
stag(3)=-0.5
|
|
|
749 |
endif
|
|
|
750 |
call putdat(cdfid,'Z',time(n),0,zz,ierr)
|
|
|
751 |
endif
|
|
|
752 |
|
|
|
753 |
C Calculate the secondary data variables
|
|
|
754 |
C --------------------------------------
|
|
|
755 |
|
|
|
756 |
C Calculation of potential temperature
|
|
|
757 |
|
|
|
758 |
if (TH_calc.eq.1) then
|
|
|
759 |
allocate(th(nx,ny,nz),STAT=stat)
|
|
|
760 |
if (stat.ne.0) print*,'error allocating th'
|
|
|
761 |
call pottemp(th,tt,sp,nx,ny,nz,ak,bk)
|
|
|
762 |
if (TH_out.eq.1) then
|
|
|
763 |
if (n.eq.1) then
|
|
|
764 |
call putdef(cdfid,'TH',4,mdv,
|
|
|
765 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
766 |
write(*,*)'*** variable TH created on ',trim(outfnam)
|
|
|
767 |
endif
|
|
|
768 |
call putdat(cdfid,'TH',time(n),0,th,ierr)
|
|
|
769 |
endif
|
|
|
770 |
endif
|
|
|
771 |
|
|
|
772 |
C Calculation of relative humidity
|
|
|
773 |
|
|
|
774 |
if (RH_calc.eq.1) then
|
|
|
775 |
allocate(rh(nx,ny,nz),STAT=stat)
|
|
|
776 |
if (stat.ne.0) print*,'error allocating rh'
|
|
|
777 |
call relhum(rh,qq,tt,sp,nx,ny,nz,ak,bk)
|
|
|
778 |
if (RH_out.eq.1) then
|
|
|
779 |
if (n.eq.1) then
|
|
|
780 |
call putdef(cdfid,'RH',4,mdv,
|
|
|
781 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
782 |
write(*,*)'*** variable RH created on ',trim(outfnam)
|
|
|
783 |
endif
|
|
|
784 |
call putdat(cdfid,'RH',time(n),0,rh,ierr)
|
|
|
785 |
endif
|
|
|
786 |
endif
|
|
|
787 |
|
|
|
788 |
C Calculation of potential vorticity (equals PV3 in IVE)
|
|
|
789 |
|
|
|
790 |
if (PV_calc.eq.1) then
|
|
|
791 |
allocate(pv(nx,ny,nz),STAT=stat)
|
|
|
792 |
if (stat.ne.0) print*,'error allocating pv'
|
|
|
793 |
call potvort(pv,uu,vv,th,sp,cl,f,nx,ny,nz,ak,bk,
|
|
|
794 |
> varmin,varmax)
|
|
|
795 |
if (PV_out.eq.1) then
|
|
|
796 |
if (n.eq.1) then
|
|
|
797 |
call putdef(cdfid,'PV',4,mdv,
|
|
|
798 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
799 |
write(*,*)'*** variable PV created on ',trim(outfnam)
|
|
|
800 |
endif
|
|
|
801 |
call putdat(cdfid,'PV',time(n),0,pv,ierr)
|
|
|
802 |
endif
|
|
|
803 |
endif
|
|
|
804 |
|
|
|
805 |
C Calculation of LAIT potential vorticity
|
|
|
806 |
|
|
|
807 |
if (LPV_calc.eq.1) then
|
|
|
808 |
allocate(lpv(nx,ny,nz),STAT=stat)
|
|
|
809 |
if (stat.ne.0) print*,'error allocating lpv'
|
|
|
810 |
call laitpv(lpv,pv,th,nx,ny,nz)
|
|
|
811 |
if (LPV_out.eq.1) then
|
|
|
812 |
if (n.eq.1) then
|
|
|
813 |
call putdef(cdfid,'LPV',4,mdv,
|
|
|
814 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
815 |
write(*,*)'*** variable LPV created on ',trim(outfnam)
|
|
|
816 |
endif
|
|
|
817 |
call putdat(cdfid,'LPV',time(n),0,lpv,ierr)
|
|
|
818 |
endif
|
|
|
819 |
endif
|
|
|
820 |
|
|
|
821 |
C Calculation of equivalent potential temperature
|
|
|
822 |
|
|
|
823 |
if (THE_calc.eq.1) then
|
|
|
824 |
allocate(the(nx,ny,nz),STAT=stat)
|
|
|
825 |
if (stat.ne.0) print*,'error allocating the'
|
|
|
826 |
call equpot(the,tt,qq,sp,nx,ny,nz,ak,bk)
|
|
|
827 |
if (THE_out.eq.1) then
|
|
|
828 |
if (n.eq.1) then
|
|
|
829 |
call putdef(cdfid,'THE',4,mdv,
|
|
|
830 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
831 |
write(*,*)'*** variable THE created on ',trim(outfnam)
|
|
|
832 |
endif
|
|
|
833 |
call putdat(cdfid,'THE',time(n),0,the,ierr)
|
|
|
834 |
endif
|
|
|
835 |
endif
|
|
|
836 |
|
|
|
837 |
C Calculation of wet-bulb potential temperature
|
|
|
838 |
|
|
|
839 |
if (THW_calc.eq.1) then
|
|
|
840 |
call wetbpt(var,the,sp,nx,ny,nz,ak,bk)
|
|
|
841 |
if (THW_out.eq.1) then
|
|
|
842 |
if (n.eq.1) then
|
|
|
843 |
call putdef(cdfid,'THW',4,mdv,
|
|
|
844 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
845 |
write(*,*)'*** variable THW created on ',trim(outfnam)
|
|
|
846 |
endif
|
|
|
847 |
call putdat(cdfid,'THW',time(n),0,var,ierr)
|
|
|
848 |
endif
|
|
|
849 |
endif
|
|
|
850 |
|
|
|
851 |
C Calculation of water flux divergence
|
|
|
852 |
|
|
|
853 |
if (DIVQU_calc.eq.1) then
|
|
|
854 |
call divqu(var,uu,vv,qq,cl,nx,ny,nz,varmin,varmax,mdv)
|
|
|
855 |
if (DIVQU_out.eq.1) then
|
|
|
856 |
if (n.eq.1) then
|
|
|
857 |
call putdef(cdfid,'DIVQU',4,mdv,
|
|
|
858 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
859 |
write(*,*)'*** variable DIVQU created on ',
|
|
|
860 |
> trim(outfnam)
|
|
|
861 |
endif
|
|
|
862 |
call putdat(cdfid,'DIVQU',time(n),0,var,ierr)
|
|
|
863 |
endif
|
|
|
864 |
endif
|
|
|
865 |
|
|
|
866 |
C Calculation of diabatic heating rate (also DHR in IVE)
|
|
|
867 |
|
|
|
868 |
if (CH_calc.eq.1) then
|
|
|
869 |
allocate(dhr(nx,ny,nz),STAT=stat)
|
|
|
870 |
if (stat.ne.0) print*,'error allocating dhr'
|
|
|
871 |
call diabheat(dhr,th,ww,rh,sp,nx,ny,nz,ak,bk)
|
|
|
872 |
if (CH_out.eq.1) then
|
|
|
873 |
if (n.eq.1) then
|
|
|
874 |
call putdef(cdfid,'CH',4,mdv,
|
|
|
875 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
876 |
write(*,*)'*** variable CH created on ',trim(outfnam)
|
|
|
877 |
endif
|
|
|
878 |
call putdat(cdfid,'CH',time(n),0,dhr,ierr)
|
|
|
879 |
endif
|
|
|
880 |
endif
|
|
|
881 |
|
|
|
882 |
C Calculation of diabatic PV rate (also DPVR in IVE)
|
|
|
883 |
|
|
|
884 |
if (PVR_calc.eq.1) then
|
|
|
885 |
call diabpvr(var,uu,vv,dhr,sp,cl,f,nx,ny,nz,ak,bk,
|
|
|
886 |
> varmin,varmax)
|
|
|
887 |
if (PVR_out.eq.1) then
|
|
|
888 |
if (n.eq.1) then
|
|
|
889 |
call putdef(cdfid,'PVR',4,mdv,
|
|
|
890 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
891 |
write(*,*)'*** variable PVR created on ',trim(outfnam)
|
|
|
892 |
endif
|
|
|
893 |
call putdat(cdfid,'PVR',time(n),0,var,ierr)
|
|
|
894 |
endif
|
|
|
895 |
endif
|
|
|
896 |
|
|
|
897 |
C Calculation of the theta gradient
|
|
|
898 |
|
|
|
899 |
if (GRADTH_calc.eq.1) then
|
|
|
900 |
call gradth(var,th,sp,prelev,cl,nx,ny,nz,ak,bk,varmin,varmax)
|
|
|
901 |
if (GRADTH_out.eq.1) then
|
|
|
902 |
if (n.eq.1) then
|
|
|
903 |
call putdef(cdfid,'GRADTH',4,mdv,
|
|
|
904 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
905 |
write(*,*)'*** variable GRADTH created on ',trim(outfnam)
|
|
|
906 |
endif
|
|
|
907 |
call putdat(cdfid,'GRADTH',time(n),0,var,ierr)
|
|
|
908 |
endif
|
|
|
909 |
endif
|
|
|
910 |
|
|
|
911 |
C Calculation of density RHO
|
|
|
912 |
|
|
|
913 |
if (RHO_calc.eq.1) then
|
|
|
914 |
allocate(rho(nx,ny,nz),STAT=stat)
|
|
|
915 |
if (stat.ne.0) print*,'error allocating rho'
|
|
|
916 |
call calc_rho(rho,tt,sp,nx,ny,nz,ak,bk)
|
|
|
917 |
if (RHO_out.eq.1) then
|
|
|
918 |
if (n.eq.1) then
|
|
|
919 |
call putdef(cdfid,'RHO',4,mdv,
|
|
|
920 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
921 |
write(*,*)'*** variable RHO created on ',trim(outfnam)
|
|
|
922 |
endif
|
|
|
923 |
call putdat(cdfid,'RHO',time(n),0,rho,ierr)
|
|
|
924 |
endif
|
|
|
925 |
endif
|
|
|
926 |
|
|
|
927 |
C Calculation of specific volume ALPHA
|
|
|
928 |
|
|
|
929 |
if (ALPHA_calc.eq.1) then
|
|
|
930 |
print*,'calculate alpha',n
|
|
|
931 |
allocate(alpha(nx,ny,nz),STAT=stat)
|
|
|
932 |
if (stat.ne.0) print*,'error allocating rho'
|
|
|
933 |
alpha = 1. / rho
|
|
|
934 |
if (ALPHA_out.eq.1) then
|
|
|
935 |
if (n.eq.1) then
|
|
|
936 |
call putdef(cdfid,'ALPHA',4,mdv,
|
|
|
937 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
938 |
write(*,*)'*** variable ALPHA created on ',
|
|
|
939 |
> trim(outfnam)
|
|
|
940 |
endif
|
|
|
941 |
call putdat(cdfid,'ALPHA',time(n),0,alpha,ierr)
|
|
|
942 |
endif
|
|
|
943 |
endif
|
|
|
944 |
|
|
|
945 |
C Calculation of Brunt-Vaisala frequ. N^2
|
|
|
946 |
|
|
|
947 |
if (NSQ_calc.eq.1) then
|
|
|
948 |
call nsq(var,rho,th,sp,nx,ny,nz,ak,bk)
|
|
|
949 |
if (NSQ_out.eq.1) then
|
|
|
950 |
if (n.eq.1) then
|
|
|
951 |
call putdef(cdfid,'NSQ',4,mdv,
|
|
|
952 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
953 |
write(*,*)'*** variable NSQ created on ',trim(outfnam)
|
|
|
954 |
endif
|
|
|
955 |
call putdat(cdfid,'NSQ',time(n),0,var,ierr)
|
|
|
956 |
endif
|
|
|
957 |
endif
|
|
|
958 |
|
|
|
959 |
C Calculation of Brunt-Vaisala frequ. N^2 with ThetaE
|
|
|
960 |
|
|
|
961 |
if (NSQM_calc.eq.1) then
|
|
|
962 |
call nsq(var,rho,the,sp,nx,ny,nz,ak,bk)
|
|
|
963 |
if (NSQM_out.eq.1) then
|
|
|
964 |
if (n.eq.1) then
|
|
|
965 |
call putdef(cdfid,'NSQM',4,mdv,
|
|
|
966 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
967 |
write(*,*)'*** variable NSQM created on ',trim(outfnam)
|
|
|
968 |
endif
|
|
|
969 |
call putdat(cdfid,'NSQM',time(n),0,var,ierr)
|
|
|
970 |
endif
|
|
|
971 |
endif
|
|
|
972 |
|
|
|
973 |
C Calculation of vertical gradient of potential temp. (-g d(th)/dp)
|
|
|
974 |
|
|
|
975 |
if (DTHDP_calc.eq.1) then
|
|
|
976 |
call dthetadp(var,th,sp,nx,ny,nz,ak,bk)
|
|
|
977 |
if (DTHDP_out.eq.1) then
|
|
|
978 |
if (n.eq.1) then
|
|
|
979 |
call putdef(cdfid,'DTHDP',4,mdv,
|
|
|
980 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
981 |
write(*,*)'*** variable DTHDP created on ',trim(outfnam)
|
|
|
982 |
endif
|
|
|
983 |
call putdat(cdfid,'DTHDP',time(n),0,var,ierr)
|
|
|
984 |
endif
|
|
|
985 |
endif
|
|
|
986 |
|
|
|
987 |
C Calculation of vertical velocity W
|
|
|
988 |
c get the vertical wind velocity in Cart. coordinates.
|
|
|
989 |
c Omega is given in hPa/s.
|
|
|
990 |
|
|
|
991 |
if (W_calc.eq.1) then
|
|
|
992 |
where (rho.ne.0.) var=-100.*ww/(9.80616*rho)
|
|
|
993 |
if (W_out.eq.1) then
|
|
|
994 |
if (n.eq.1) then
|
|
|
995 |
call putdef(cdfid,'W',4,mdv,
|
|
|
996 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
997 |
write(*,*)'*** variable W created on ',trim(outfnam)
|
|
|
998 |
endif
|
|
|
999 |
call putdat(cdfid,'W',time(n),0,var,ierr)
|
|
|
1000 |
endif
|
|
|
1001 |
endif
|
|
|
1002 |
|
|
|
1003 |
C Calculation of velocity VEL
|
|
|
1004 |
|
|
|
1005 |
if (VEL_calc.eq.1) then
|
|
|
1006 |
var=sqrt(uu**2+vv**2)
|
|
|
1007 |
if (VEL_out.eq.1) then
|
|
|
1008 |
if (n.eq.1) then
|
|
|
1009 |
call putdef(cdfid,'VEL',4,mdv,
|
|
|
1010 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1011 |
write(*,*)'*** variable VEL created on ',trim(outfnam)
|
|
|
1012 |
endif
|
|
|
1013 |
call putdat(cdfid,'VEL',time(n),0,var,ierr)
|
|
|
1014 |
endif
|
|
|
1015 |
endif
|
|
|
1016 |
|
|
|
1017 |
C Calculation of ZLAY
|
|
|
1018 |
|
|
|
1019 |
if (ZLAY_calc.eq.1) then
|
|
|
1020 |
allocate(zlay(nx,ny,nz),STAT=stat)
|
|
|
1021 |
if (stat.ne.0) print*,'error allocating zlay'
|
|
|
1022 |
call zlayer(zlay,tt,zz,sp,nx,ny,nz,aklev,bklev,aklay,bklay)
|
|
|
1023 |
if (ZLAY_out.eq.1) then
|
|
|
1024 |
if (n.eq.1) then
|
|
|
1025 |
call putdef(cdfid,'ZLAY',4,mdv,
|
|
|
1026 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1027 |
write(*,*)'*** variable ZLAY created on ',trim(outfnam)
|
|
|
1028 |
endif
|
|
|
1029 |
call putdat(cdfid,'ZLAY',time(n),0,zlay,ierr)
|
|
|
1030 |
endif
|
|
|
1031 |
endif
|
|
|
1032 |
|
|
|
1033 |
C Calculation of Montgomery-Potential M
|
|
|
1034 |
|
|
|
1035 |
if (M_calc.eq.1) then
|
|
|
1036 |
allocate(mm(nx,ny,nz),STAT=stat)
|
|
|
1037 |
if (stat.ne.0) print*,'error allocating mm'
|
|
|
1038 |
mm=1004.*(tt+273.15)+9.80616*1000.*zlay
|
|
|
1039 |
if (M_out.eq.1) then
|
|
|
1040 |
if (n.eq.1) then
|
|
|
1041 |
call putdef(cdfid,'M',4,mdv,
|
|
|
1042 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1043 |
write(*,*)'*** variable M created on ',trim(outfnam)
|
|
|
1044 |
endif
|
|
|
1045 |
call putdat(cdfid,'M',time(n),0,mm,ierr)
|
|
|
1046 |
endif
|
|
|
1047 |
endif
|
|
|
1048 |
|
|
|
1049 |
C Calculation of Bernoulli-Function B
|
|
|
1050 |
|
|
|
1051 |
if (B_calc.eq.1) then
|
|
|
1052 |
var=mm+0.5*(uu**2+vv**2)
|
|
|
1053 |
if (B_out.eq.1) then
|
|
|
1054 |
if (n.eq.1) then
|
|
|
1055 |
call putdef(cdfid,'B',4,mdv,
|
|
|
1056 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1057 |
write(*,*)'*** variable B created on ',trim(outfnam)
|
|
|
1058 |
endif
|
|
|
1059 |
call putdat(cdfid,'B',time(n),0,var,ierr)
|
|
|
1060 |
endif
|
|
|
1061 |
endif
|
|
|
1062 |
|
|
|
1063 |
C Calculation of relative vertical vorticity VORT
|
|
|
1064 |
|
|
|
1065 |
if (VORT_calc.eq.1) then
|
|
|
1066 |
call vort(var,uu,vv,sp,cl,nx,ny,nz,ak,bk,varmin,varmax)
|
|
|
1067 |
if (VORT_out.eq.1) then
|
|
|
1068 |
if (n.eq.1) then
|
|
|
1069 |
call putdef(cdfid,'VORT',4,mdv,
|
|
|
1070 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1071 |
write(*,*)'*** variable VORT created on ',trim(outfnam)
|
|
|
1072 |
endif
|
|
|
1073 |
call putdat(cdfid,'VORT',time(n),0,var,ierr)
|
|
|
1074 |
endif
|
|
|
1075 |
endif
|
|
|
1076 |
|
|
|
1077 |
C Calculation of absolute vertical vorticity AVO
|
|
|
1078 |
|
|
|
1079 |
if (AVO_calc.eq.1) then
|
|
|
1080 |
call avort(var,uu,vv,sp,cl,f,nx,ny,nz,ak,bk,varmin,varmax)
|
|
|
1081 |
if (AVO_out.eq.1) then
|
|
|
1082 |
if (n.eq.1) then
|
|
|
1083 |
call putdef(cdfid,'AVO',4,mdv,
|
|
|
1084 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1085 |
write(*,*)'*** variable AVO created on ',trim(outfnam)
|
|
|
1086 |
endif
|
|
|
1087 |
call putdat(cdfid,'AVO',time(n),0,var,ierr)
|
|
|
1088 |
endif
|
|
|
1089 |
endif
|
|
|
1090 |
|
|
|
1091 |
C Calculation of vertical curvature vorticity CURVO
|
|
|
1092 |
|
|
|
1093 |
if (CURVO_calc.eq.1) then
|
|
|
1094 |
call curvo(var,uu,vv,sp,cl,nx,ny,nz,ak,bk,varmin,varmax)
|
|
|
1095 |
if (CURVO_out.eq.1) then
|
|
|
1096 |
if (n.eq.1) then
|
|
|
1097 |
call putdef(cdfid,'CURVO',4,mdv,
|
|
|
1098 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1099 |
write(*,*)'*** variable CURVO created on ',trim(outfnam)
|
|
|
1100 |
endif
|
|
|
1101 |
call putdat(cdfid,'CURVO',time(n),0,var,ierr)
|
|
|
1102 |
endif
|
|
|
1103 |
endif
|
|
|
1104 |
|
|
|
1105 |
C Output of cos(latitude) array
|
|
|
1106 |
|
|
|
1107 |
if (COS_out.eq.1) then
|
|
|
1108 |
if (n.eq.1) then
|
|
|
1109 |
call putdef(cdfid,'COS',4,mdv,
|
|
|
1110 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1111 |
write(*,*)'*** variable COS created on ',trim(outfnam)
|
|
|
1112 |
endif
|
|
|
1113 |
call putdat(cdfid,'COS',time(n),0,cl,ierr)
|
|
|
1114 |
endif
|
|
|
1115 |
|
|
|
1116 |
C Output of vertically integrated PV (VIP)
|
|
|
1117 |
|
|
|
1118 |
if (VIP_out.eq.1) then
|
|
|
1119 |
allocate(vip(nx,ny),STAT=stat)
|
|
|
1120 |
if (stat.ne.0) print*,'error allocating vip'
|
|
|
1121 |
call vint(vip,pv,sp,150.,500.,nx,ny,nz,ak,bk)
|
|
|
1122 |
if (n.eq.1) then
|
|
|
1123 |
vardim(3)=1
|
|
|
1124 |
call putdef(cdfid,'VIP',4,mdv,
|
|
|
1125 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1126 |
vardim(3)=nz
|
|
|
1127 |
write(*,*)'*** variable VIP created on ',trim(outfnam)
|
|
|
1128 |
endif
|
|
|
1129 |
call putdat(cdfid,'VIP',time(n),0,vip,ierr)
|
|
|
1130 |
endif
|
|
|
1131 |
|
|
|
1132 |
C Calculation of PV gradient (GRADPV)
|
|
|
1133 |
|
|
|
1134 |
if (GRADPV_out.eq.1) then
|
|
|
1135 |
C write(*,*)'*** calling gradpv ***'
|
|
|
1136 |
allocate(gradpv(nx,ny),STAT=stat)
|
|
|
1137 |
if (stat.ne.0) print*,'error allocating gradpv'
|
|
|
1138 |
call calc_gradpv(gradpv,ipv,cl,nx,ny,nz,akd,bkd,
|
|
|
1139 |
> varmin,varmax)
|
|
|
1140 |
if (n.eq.1) then
|
|
|
1141 |
vardim(3)=1
|
|
|
1142 |
call putdef(cdfid,'GRADPV',4,mdv,
|
|
|
1143 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1144 |
vardim(3)=nz
|
|
|
1145 |
write(*,*)'*** variable GRADPV created on ',trim(outfnam)
|
|
|
1146 |
endif
|
|
|
1147 |
call putdat(cdfid,'GRADPV',time(n),0,gradpv,ierr)
|
|
|
1148 |
endif
|
|
|
1149 |
|
|
|
1150 |
C Calculation of P
|
|
|
1151 |
|
|
|
1152 |
if (P_calc.eq.1) then
|
|
|
1153 |
call pressure(var,sp,-0.5,nx,ny,nz,aklev,bklev,aklay,bklay)
|
|
|
1154 |
if (P_out.eq.1) then
|
|
|
1155 |
if (n.eq.1) then
|
|
|
1156 |
call putdef(cdfid,'P',4,mdv,
|
|
|
1157 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1158 |
write(*,*)'*** variable P created on ',trim(outfnam)
|
|
|
1159 |
endif
|
|
|
1160 |
call putdat(cdfid,'P',time(n),0,var,ierr)
|
|
|
1161 |
endif
|
|
|
1162 |
endif
|
|
|
1163 |
|
|
|
1164 |
C Calculation of PLEV
|
|
|
1165 |
|
|
|
1166 |
if (PLEV_calc.eq.1) then
|
|
|
1167 |
call pressure(var,sp,0.,nx,ny,nz,aklev,bklev,aklay,bklay)
|
|
|
1168 |
if (PLEV_out.eq.1) then
|
|
|
1169 |
if (n.eq.1) then
|
|
|
1170 |
call putdef(cdfid,'PLEV',4,mdv,
|
|
|
1171 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1172 |
write(*,*)'*** variable PLEV created on ',trim(outfnam)
|
|
|
1173 |
endif
|
|
|
1174 |
call putdat(cdfid,'PLEV',time(n),0,var,ierr)
|
|
|
1175 |
endif
|
|
|
1176 |
endif
|
|
|
1177 |
|
|
|
1178 |
C Calculation of FL
|
|
|
1179 |
|
|
|
1180 |
if (FL_calc.eq.1) then
|
|
|
1181 |
allocate(fl(nx,ny,nz),STAT=stat)
|
|
|
1182 |
if (stat.ne.0) print*,'error allocating fl'
|
|
|
1183 |
call pres(var,sp,nx,ny,nz,aklay,bklay)
|
|
|
1184 |
call flightlev(fl,var,nx,ny,nz)
|
|
|
1185 |
if (FL_out.eq.1) then
|
|
|
1186 |
if (n.eq.1) then
|
|
|
1187 |
call putdef(cdfid,'FL',4,mdv,
|
|
|
1188 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1189 |
write(*,*)'*** variable FL created on ',trim(outfnam)
|
|
|
1190 |
endif
|
|
|
1191 |
call putdat(cdfid,'FL',time(n),0,fl,ierr)
|
|
|
1192 |
endif
|
|
|
1193 |
endif
|
|
|
1194 |
|
|
|
1195 |
C Calculation of geostrophic wind UG
|
|
|
1196 |
|
|
|
1197 |
if (UG_calc.eq.1) then
|
|
|
1198 |
allocate(ug(nx,ny,nz),STAT=stat)
|
|
|
1199 |
if (stat.ne.0) print*,'error allocating ug'
|
|
|
1200 |
call calc_ug(ug,zlay,sp,cl,f,nx,ny,nz,ak,bk,varmin,varmax)
|
|
|
1201 |
if (UG_out.eq.1) then
|
|
|
1202 |
if (n.eq.1) then
|
|
|
1203 |
call putdef(cdfid,'UG',4,mdv,
|
|
|
1204 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1205 |
write(*,*)'*** variable UG created on ',trim(outfnam)
|
|
|
1206 |
endif
|
|
|
1207 |
call putdat(cdfid,'UG',time(n),0,ug,ierr)
|
|
|
1208 |
endif
|
|
|
1209 |
endif
|
|
|
1210 |
|
|
|
1211 |
C Calculation of geostrophic wind VG
|
|
|
1212 |
|
|
|
1213 |
if (VG_calc.eq.1) then
|
|
|
1214 |
allocate(vg(nx,ny,nz),STAT=stat)
|
|
|
1215 |
if (stat.ne.0) print*,'error allocating vg'
|
|
|
1216 |
call calc_vg(vg,zlay,sp,cl,f,nx,ny,nz,ak,bk,varmin,varmax)
|
|
|
1217 |
if (VG_out.eq.1) then
|
|
|
1218 |
if (n.eq.1) then
|
|
|
1219 |
call putdef(cdfid,'VG',4,mdv,
|
|
|
1220 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1221 |
write(*,*)'*** variable VG created on ',trim(outfnam)
|
|
|
1222 |
endif
|
|
|
1223 |
call putdat(cdfid,'VG',time(n),0,vg,ierr)
|
|
|
1224 |
endif
|
|
|
1225 |
endif
|
|
|
1226 |
|
|
|
1227 |
C Calculation of ageostrophic geostrophic wind UA
|
|
|
1228 |
|
|
|
1229 |
if (UA_calc.eq.1) then
|
|
|
1230 |
var=uu-ug
|
|
|
1231 |
if (UA_out.eq.1) then
|
|
|
1232 |
if (n.eq.1) then
|
|
|
1233 |
call putdef(cdfid,'UA',4,mdv,
|
|
|
1234 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1235 |
write(*,*)'*** variable UA created on ',trim(outfnam)
|
|
|
1236 |
endif
|
|
|
1237 |
call putdat(cdfid,'UA',time(n),0,var,ierr)
|
|
|
1238 |
endif
|
|
|
1239 |
endif
|
|
|
1240 |
|
|
|
1241 |
C Calculation of ageostrophic geostrophic wind VA
|
|
|
1242 |
|
|
|
1243 |
if (VA_calc.eq.1) then
|
|
|
1244 |
var=vv-vg
|
|
|
1245 |
if (VA_out.eq.1) then
|
|
|
1246 |
if (n.eq.1) then
|
|
|
1247 |
call putdef(cdfid,'VA',4,mdv,
|
|
|
1248 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1249 |
write(*,*)'*** variable VA created on ',trim(outfnam)
|
|
|
1250 |
endif
|
|
|
1251 |
call putdat(cdfid,'VA',time(n),0,var,ierr)
|
|
|
1252 |
endif
|
|
|
1253 |
endif
|
|
|
1254 |
|
|
|
1255 |
|
|
|
1256 |
C Calculation of x-component of the Q-vector (calc. with U)
|
|
|
1257 |
|
|
|
1258 |
if (QXF_calc.eq.1) then
|
|
|
1259 |
call calc_qx(var,uu,vv,tt,sp,cl,nx,ny,nz,ak,bk,varmin,varmax)
|
|
|
1260 |
if (QXF_out.eq.1) then
|
|
|
1261 |
if (n.eq.1) then
|
|
|
1262 |
call putdef(cdfid,'QXF',4,mdv,
|
|
|
1263 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1264 |
write(*,*)'*** variable QXF created on ',trim(outfnam)
|
|
|
1265 |
endif
|
|
|
1266 |
call putdat(cdfid,'QXF',time(n),0,var,ierr)
|
|
|
1267 |
endif
|
|
|
1268 |
endif
|
|
|
1269 |
|
|
|
1270 |
C Calculation of y-component of the Q-vector (calc. with V)
|
|
|
1271 |
|
|
|
1272 |
if (QYF_calc.eq.1) then
|
|
|
1273 |
call calc_qy(var,uu,vv,tt,sp,cl,tl,nx,ny,nz,ak,bk,varmin,varmax)
|
|
|
1274 |
if (QYF_out.eq.1) then
|
|
|
1275 |
if (n.eq.1) then
|
|
|
1276 |
call putdef(cdfid,'QYF',4,mdv,
|
|
|
1277 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1278 |
write(*,*)'*** variable QYF created on ',trim(outfnam)
|
|
|
1279 |
endif
|
|
|
1280 |
call putdat(cdfid,'QYF',time(n),0,var,ierr)
|
|
|
1281 |
endif
|
|
|
1282 |
endif
|
|
|
1283 |
|
|
|
1284 |
C Calculation of x-component of the Q-vector (calc. with UG)
|
|
|
1285 |
|
|
|
1286 |
if (QX_calc.eq.1) then
|
|
|
1287 |
call calc_qx(var,ug,vg,tt,sp,cl,nx,ny,nz,ak,bk,varmin,varmax)
|
|
|
1288 |
if (QX_out.eq.1) then
|
|
|
1289 |
if (n.eq.1) then
|
|
|
1290 |
call putdef(cdfid,'QX',4,mdv,
|
|
|
1291 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1292 |
write(*,*)'*** variable QX created on ',trim(outfnam)
|
|
|
1293 |
endif
|
|
|
1294 |
call putdat(cdfid,'QX',time(n),0,var,ierr)
|
|
|
1295 |
endif
|
|
|
1296 |
endif
|
|
|
1297 |
|
|
|
1298 |
C Calculation of y-component of the Q-vector (calc. with VG)
|
|
|
1299 |
|
|
|
1300 |
if (QY_calc.eq.1) then
|
|
|
1301 |
call calc_qy(var,ug,vg,tt,sp,cl,tl,nx,ny,nz,ak,bk,varmin,varmax)
|
|
|
1302 |
if (QY_out.eq.1) then
|
|
|
1303 |
if (n.eq.1) then
|
|
|
1304 |
call putdef(cdfid,'QY',4,mdv,
|
|
|
1305 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1306 |
write(*,*)'*** variable QY created on ',trim(outfnam)
|
|
|
1307 |
endif
|
|
|
1308 |
call putdat(cdfid,'QY',time(n),0,var,ierr)
|
|
|
1309 |
endif
|
|
|
1310 |
endif
|
|
|
1311 |
|
|
|
1312 |
C Calculation of reduced surface pressure PSRED
|
|
|
1313 |
|
|
|
1314 |
if (PSRED_calc.eq.1) then
|
|
|
1315 |
if (.not.ZBneed) then
|
|
|
1316 |
call calc_psred(var,sp,tt,zb,nx,ny,nz,aklay,bklay)
|
|
|
1317 |
if (PSRED_out.eq.1) then
|
|
|
1318 |
if (n.eq.1) then
|
|
|
1319 |
vardim2=vardim
|
|
|
1320 |
vardim2(3)=1 ! force one vertical-level
|
|
|
1321 |
call putdef(cdfid,'PSRED',4,mdv,
|
|
|
1322 |
> vardim2,varmin,varmax,stag,ierr)
|
|
|
1323 |
write(*,*)'*** variable PSRED created on ',trim(outfnam)
|
|
|
1324 |
endif
|
|
|
1325 |
call putdat(cdfid,'PSRED',time(n),0,var,ierr)
|
|
|
1326 |
endif
|
|
|
1327 |
else
|
|
|
1328 |
print*,'ZB needed to calculate PSRED'
|
|
|
1329 |
endif
|
|
|
1330 |
endif
|
|
|
1331 |
|
|
|
1332 |
C Calculation of Richardson number
|
|
|
1333 |
|
|
|
1334 |
if (RI_calc.eq.1) then
|
|
|
1335 |
call rich(var,sp,uu,vv,th,nx,ny,nz,ak,bk)
|
|
|
1336 |
if (RI_out.eq.1) then
|
|
|
1337 |
if (n.eq.1) then
|
|
|
1338 |
call putdef(cdfid,'RI',4,mdv,
|
|
|
1339 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1340 |
write(*,*)'*** variable RI created on ',trim(outfnam)
|
|
|
1341 |
endif
|
|
|
1342 |
call putdat(cdfid,'RI',time(n),0,var,ierr)
|
|
|
1343 |
endif
|
|
|
1344 |
endif
|
|
|
1345 |
|
|
|
1346 |
C Calculation of boundary layer height BLH
|
|
|
1347 |
|
|
|
1348 |
if (BLH_calc.eq.1) then
|
|
|
1349 |
call calc_blh(var,sp,tt,qq,uu,vv,t2m,td2m,u10m,v10m,
|
|
|
1350 |
> nx,ny,nz,aklay,bklay)
|
|
|
1351 |
if (BLH_out.eq.1) then
|
|
|
1352 |
if (n.eq.1) then
|
|
|
1353 |
vardim2=vardim
|
|
|
1354 |
vardim2(3)=1 ! force one vertical-level
|
|
|
1355 |
call putdef(cdfid,'BLH',4,mdv,
|
|
|
1356 |
> vardim2,varmin,varmax,stag,ierr)
|
|
|
1357 |
write(*,*)'*** variable BLH created on ',trim(outfnam)
|
|
|
1358 |
endif
|
|
|
1359 |
call putdat(cdfid,'BLH',time(n),0,var,ierr)
|
|
|
1360 |
endif
|
|
|
1361 |
endif
|
|
|
1362 |
|
|
|
1363 |
c add calculations of your own variables here: (++)
|
|
|
1364 |
|
|
|
1365 |
enddo
|
|
|
1366 |
|
|
|
1367 |
C Close the NetCDF files
|
|
|
1368 |
|
|
|
1369 |
call clscdf(cdfid,ierr)
|
|
|
1370 |
900 continue
|
|
|
1371 |
call clscdf(cdfid1,ierr)
|
|
|
1372 |
call clscdf(cstid,ierr)
|
|
|
1373 |
|
|
|
1374 |
goto 999
|
|
|
1375 |
|
|
|
1376 |
920 stop '*** error: variable T not found on P-file ***'
|
|
|
1377 |
921 stop '*** error: variable PS not found on P-file ***'
|
|
|
1378 |
922 stop '*** error: variable Q not found on P-file ***'
|
|
|
1379 |
923 stop '*** error: variable U not found on P-file ***'
|
|
|
1380 |
924 stop '*** error: variable V not found on P-file ***'
|
|
|
1381 |
925 stop '*** error: variable OMEGA not found on P-file ***'
|
|
|
1382 |
926 stop '*** error: variable T2M not found on P-file ***'
|
|
|
1383 |
927 stop '*** error: variable TD2M not found on P-file ***'
|
|
|
1384 |
928 stop '*** error: variable U10M not found on P-file ***'
|
|
|
1385 |
929 stop '*** error: variable V10M not found on P-file ***'
|
|
|
1386 |
930 stop '*** error: variable ZB not found on P-file ***'
|
|
|
1387 |
|
|
|
1388 |
996 stop '*** error: could not create S-file ***'
|
|
|
1389 |
999 continue
|
|
|
1390 |
end
|
|
|
1391 |
|
|
|
1392 |
c-----------------------------------------------------------------------
|
|
|
1393 |
c following subroutines
|
|
|
1394 |
c-----------------------------------------------------------------------
|
|
|
1395 |
|
|
|
1396 |
subroutine calc_rho(var,tt,sp,ie,je,ke,ak,bk)
|
|
|
1397 |
c ============================================
|
|
|
1398 |
c computation of density
|
|
|
1399 |
c
|
|
|
1400 |
c argument declaration
|
|
|
1401 |
integer ie,je,ke
|
|
|
1402 |
real,intent(OUT) :: var(ie,je,ke)
|
|
|
1403 |
real,intent(IN) :: tt(ie,je,ke),sp(ie,je)
|
|
|
1404 |
real,intent(IN) :: ak(ke),bk(ke)
|
|
|
1405 |
|
|
|
1406 |
integer i,j,k
|
|
|
1407 |
real,parameter :: tzero=273.15
|
|
|
1408 |
|
|
|
1409 |
c statement-functions for the computation of pressure
|
|
|
1410 |
real pp,psrf
|
|
|
1411 |
integer is
|
|
|
1412 |
pp(is)=ak(is)+bk(is)*psrf
|
|
|
1413 |
|
|
|
1414 |
c computation of potential temperature
|
|
|
1415 |
do i=1,ie
|
|
|
1416 |
do j=1,je
|
|
|
1417 |
psrf=sp(i,j)
|
|
|
1418 |
do k=1,ke
|
|
|
1419 |
c distinction of temperature in K and deg C
|
|
|
1420 |
if (tt(i,j,k).lt.100.) then
|
|
|
1421 |
var(i,j,k)=pp(k)/(2.87*(tt(i,j,k)+tzero))
|
|
|
1422 |
else
|
|
|
1423 |
var(i,j,k)=pp(k)/(2.87*tt(i,j,k))
|
|
|
1424 |
endif
|
|
|
1425 |
enddo
|
|
|
1426 |
enddo
|
|
|
1427 |
enddo
|
|
|
1428 |
end
|
|
|
1429 |
|
|
|
1430 |
|
|
|
1431 |
subroutine nsq(var,rho,th,sp,ie,je,ke,ak,bk)
|
|
|
1432 |
C ======================================================
|
|
|
1433 |
|
|
|
1434 |
c argument declaration
|
|
|
1435 |
integer,intent(IN) :: ie,je,ke
|
|
|
1436 |
real,intent(OUT) :: var(ie,je,ke)
|
|
|
1437 |
real,intent(IN) :: rho(ie,je,ke),th(ie,je,ke),sp(ie,je)
|
|
|
1438 |
real,intent(IN) :: ak(ke),bk(ke)
|
|
|
1439 |
|
|
|
1440 |
c variable declaration
|
|
|
1441 |
integer stat
|
|
|
1442 |
REAL,ALLOCATABLE, DIMENSION (:,:,:) :: dthdp !3D array
|
|
|
1443 |
|
|
|
1444 |
allocate(dthdp(ie,je,ke),STAT=stat)
|
|
|
1445 |
IF (stat.ne.0) PRINT*,'nsq: error allocating dthdp(ie,je,ke)'
|
|
|
1446 |
|
|
|
1447 |
call ddp(th,dthdp,sp,ie,je,ke,ak,bk)
|
|
|
1448 |
|
|
|
1449 |
where(th.ne.0.) var=-96.24*rho/th*dthdp
|
|
|
1450 |
|
|
|
1451 |
IF (ALLOCATED(dthdp)) DEALLOCATE(dthdp)
|
|
|
1452 |
|
|
|
1453 |
end
|
|
|
1454 |
|
|
|
1455 |
subroutine dthetadp(var,th,sp,ie,je,ke,ak,bk)
|
|
|
1456 |
C ======================================================
|
|
|
1457 |
|
|
|
1458 |
c argument declaration
|
|
|
1459 |
integer,intent(IN) :: ie,je,ke
|
|
|
1460 |
real,intent(OUT) :: var(ie,je,ke)
|
|
|
1461 |
real,intent(IN) :: th(ie,je,ke),sp(ie,je)
|
|
|
1462 |
real,intent(IN) :: ak(ke),bk(ke)
|
|
|
1463 |
|
|
|
1464 |
c variable declaration
|
|
|
1465 |
integer stat
|
|
|
1466 |
REAL,ALLOCATABLE, DIMENSION (:,:,:) :: dthdp !3D array
|
|
|
1467 |
|
|
|
1468 |
allocate(dthdp(ie,je,ke),STAT=stat)
|
|
|
1469 |
IF (stat.ne.0) PRINT*,'nsq: error allocating dthdp(ie,je,ke)'
|
|
|
1470 |
|
|
|
1471 |
call ddp(th,dthdp,sp,ie,je,ke,ak,bk)
|
|
|
1472 |
|
|
|
1473 |
c factor 100 is because VORT is in units of f and so VORT*DTHDP
|
|
|
1474 |
c is in pvu
|
|
|
1475 |
var=-9.80616*100.*dthdp
|
|
|
1476 |
|
|
|
1477 |
IF (ALLOCATED(dthdp)) DEALLOCATE(dthdp)
|
|
|
1478 |
|
|
|
1479 |
end
|
|
|
1480 |
|
|
|
1481 |
subroutine geopot(psi,q,t,oro,sp,ie,je,ke,ak,bk)
|
|
|
1482 |
c ================================================
|
|
|
1483 |
|
|
|
1484 |
c argument declaration
|
|
|
1485 |
integer ie,je,ke
|
|
|
1486 |
real psi(ie,je,ke),t(ie,je,ke),q(ie,je,ke),oro(ie,je),
|
|
|
1487 |
> sp(ie,je),ak(ke),bk(ke)
|
|
|
1488 |
|
|
|
1489 |
c variable declaration
|
|
|
1490 |
integer i,j,k
|
|
|
1491 |
real r,c,g
|
|
|
1492 |
data r,c,g /287.,0.608,9.80616/
|
|
|
1493 |
|
|
|
1494 |
c statement-functions for the computation of pressure
|
|
|
1495 |
real pp,psrf
|
|
|
1496 |
integer is
|
|
|
1497 |
pp(is)=ak(is)+bk(is)*psrf
|
|
|
1498 |
|
|
|
1499 |
c integration of geopotential height(special for first layer)
|
|
|
1500 |
do i=1,ie
|
|
|
1501 |
do j=1,je
|
|
|
1502 |
psrf=sp(i,j)
|
|
|
1503 |
psi(i,j,1)=1./g*(oro(i,j)
|
|
|
1504 |
> +r*(t(i,j,1)+273.15)*(1.+c*q(i,j,1))*
|
|
|
1505 |
> (psrf-pp(1))/(0.5*(psrf+pp(1))))
|
|
|
1506 |
enddo
|
|
|
1507 |
enddo
|
|
|
1508 |
do j=1,je
|
|
|
1509 |
do i=1,ie
|
|
|
1510 |
psrf=sp(i,j)
|
|
|
1511 |
do k=2,ke
|
|
|
1512 |
psi(i,j,k)=psi(i,j,k-1)+r/g*
|
|
|
1513 |
> ((t(i,j,k-1)+273.15)*(1.+c*q(i,j,k-1))+
|
|
|
1514 |
> (t(i,j,k)+273.15)*(1.+c*q(i,j,k)))*
|
|
|
1515 |
> (pp(k-1)-pp(k))/(pp(k-1)+pp(k))
|
|
|
1516 |
enddo
|
|
|
1517 |
enddo
|
|
|
1518 |
enddo
|
|
|
1519 |
end
|
|
|
1520 |
|
|
|
1521 |
|
|
|
1522 |
subroutine getoro(or,dx,dy,starty,lonmin,latmin,ie,je,ierr)
|
|
|
1523 |
c ===========================================================
|
|
|
1524 |
c reads the orography for the actual data domain from the file
|
|
|
1525 |
c "/home/henry/ecmwf/cdf/orogr", which contains the orography
|
|
|
1526 |
c for the whole northern hemisphere.
|
|
|
1527 |
|
|
|
1528 |
c argument declaration
|
|
|
1529 |
integer ie,je
|
|
|
1530 |
real or(ie,je)
|
|
|
1531 |
|
|
|
1532 |
c variable declaration
|
|
|
1533 |
integer i,j,cdfid,ierr,i0,j0,err,starty,nx
|
|
|
1534 |
real gloro(480*240)
|
|
|
1535 |
real dx,dy,dd,time,lonmin,latmin
|
|
|
1536 |
character*(80)filnam
|
|
|
1537 |
integer vardim(3),ndim
|
|
|
1538 |
real varmin(3),varmax(3),stag(3),mdv
|
|
|
1539 |
|
|
|
1540 |
time=0.
|
|
|
1541 |
|
|
|
1542 |
c open file with orography values and get them
|
|
|
1543 |
call ncpopt(NCVERBOS)
|
|
|
1544 |
|
|
|
1545 |
c q&d bug fix (dx, dy were changed below when calculating i0)
|
|
|
1546 |
dd=dx
|
|
|
1547 |
if ((dx.eq.0.75).and.(dy.eq.0.75)) then
|
|
|
1548 |
if (starty.lt.96) then
|
|
|
1549 |
if (starty.le.92) then
|
|
|
1550 |
filnam="/home/henry/cdf/oro/9209orogr0_75"
|
|
|
1551 |
else if (starty.le.93) then
|
|
|
1552 |
filnam="/home/henry/cdf/oro/9309orogr0_75"
|
|
|
1553 |
else if (starty.le.94) then
|
|
|
1554 |
filnam="/home/henry/cdf/oro/9411orogr0_75"
|
|
|
1555 |
else if (starty.le.95) then
|
|
|
1556 |
filnam="/home/henry/cdf/oro/9509orogr0_75"
|
|
|
1557 |
else
|
|
|
1558 |
filnam="/home/henry/cdf/oro/orogr0_75"
|
|
|
1559 |
endif
|
|
|
1560 |
else
|
|
|
1561 |
filnam="/home/henry/cdf/oro/95orogr0_75"
|
|
|
1562 |
endif
|
|
|
1563 |
else if ((dx.eq.1.0).and.(dy.eq.1.0)) then
|
|
|
1564 |
if (starty.lt.95) then
|
|
|
1565 |
filnam="/home/henry/cdf/oro/orogr1_00"
|
|
|
1566 |
else
|
|
|
1567 |
filnam="/home/henry/cdf/oro/gloro1_00"
|
|
|
1568 |
endif
|
|
|
1569 |
else if ((dx.eq.0.5).and.(dy.eq.0.5)) then
|
|
|
1570 |
filnam="/home/henry/cdf/oro/95orogr0_50"
|
|
|
1571 |
else
|
|
|
1572 |
print*,'unappropriate data for the orography file'
|
|
|
1573 |
print*,'grid interval should be 0.5, 0.75 or 1.00
|
|
|
1574 |
> and data only on the NH'
|
|
|
1575 |
ierr=2
|
|
|
1576 |
return
|
|
|
1577 |
endif
|
|
|
1578 |
print*,filnam
|
|
|
1579 |
call cdfopn(trim(filnam),cdfid,err)
|
|
|
1580 |
call getdat(cdfid,'ORO',time,0,gloro,err)
|
|
|
1581 |
call getdef(cdfid,'ORO',ndim,mdv,vardim,varmin,varmax,stag,err)
|
|
|
1582 |
* call clscdf(cdfid,err)
|
|
|
1583 |
|
|
|
1584 |
i0=nint((lonmin-varmin(1))/dd)
|
|
|
1585 |
j0=nint((latmin-varmin(2))/dd)
|
|
|
1586 |
if (j0.lt.0) j0=-j0
|
|
|
1587 |
|
|
|
1588 |
nx=nint(360./dd)
|
|
|
1589 |
do i=1,ie
|
|
|
1590 |
do j=1,je
|
|
|
1591 |
if (i0+i.le.nx) then
|
|
|
1592 |
or(i,j)=gloro(i0+i+(j0+j-1)*nx)
|
|
|
1593 |
else
|
|
|
1594 |
or(i,j)=gloro(i0+i-nx+(j0+j-1)*nx)
|
|
|
1595 |
endif
|
|
|
1596 |
if (or(i,j).gt.1.e20) print*,i,j,i0+i+(j0+j-1)*nx
|
|
|
1597 |
enddo
|
|
|
1598 |
enddo
|
|
|
1599 |
end
|
|
|
1600 |
|
|
|
1601 |
subroutine calc_qx(var,uu,vv,tt,sp,cl,ie,je,ke,ak,bk,vmin,vmax)
|
|
|
1602 |
C ===============================================================
|
|
|
1603 |
c calculate x-comp Q-vector: -9.8/273.*(D[U:X]*D[T:X]+D[V:X]*D[T:Y])
|
|
|
1604 |
|
|
|
1605 |
c argument declaration
|
|
|
1606 |
integer ie,je,ke
|
|
|
1607 |
real,intent(OUT) :: var(ie,je,ke)
|
|
|
1608 |
real,intent(IN) :: uu(ie,je,ke),vv(ie,je,ke),tt(ie,je,ke)
|
|
|
1609 |
real,intent(IN) :: sp(ie,je),cl(ie,je)
|
|
|
1610 |
real,intent(IN) :: ak(ke),bk(ke),vmin(4),vmax(4)
|
|
|
1611 |
|
|
|
1612 |
c variable declaration
|
|
|
1613 |
REAL,ALLOCATABLE, DIMENSION (:,:) :: dspdx,dspdy
|
|
|
1614 |
REAL,ALLOCATABLE, DIMENSION (:,:,:) :: dudx,dtdx,dvdx,dtdy
|
|
|
1615 |
integer stat
|
|
|
1616 |
integer k
|
|
|
1617 |
real mdv
|
|
|
1618 |
logical prelev
|
|
|
1619 |
|
|
|
1620 |
c test if pressure or model levels
|
|
|
1621 |
prelev=.true.
|
|
|
1622 |
do k=1,ke
|
|
|
1623 |
if (bk(k).ne.0.) prelev=.false.
|
|
|
1624 |
enddo
|
|
|
1625 |
mdv=-999.98999
|
|
|
1626 |
|
|
|
1627 |
if (.not.prelev) then
|
|
|
1628 |
allocate(dspdx(ie,je),STAT=stat)
|
|
|
1629 |
if (stat.ne.0) print*,'calc_qx: error allocating dspdx'
|
|
|
1630 |
allocate(dspdy(ie,je),STAT=stat)
|
|
|
1631 |
if (stat.ne.0) print*,'calc_qx: error allocating dspdy'
|
|
|
1632 |
endif
|
|
|
1633 |
|
|
|
1634 |
allocate(dudx(ie,je,ke),STAT=stat)
|
|
|
1635 |
if (stat.ne.0) print*,'calc_qx: error allocating dudx'
|
|
|
1636 |
allocate(dtdx(ie,je,ke),STAT=stat)
|
|
|
1637 |
if (stat.ne.0) print*,'calc_qx: error allocating dtdx'
|
|
|
1638 |
allocate(dvdx(ie,je,ke),STAT=stat)
|
|
|
1639 |
if (stat.ne.0) print*,'calc_qx: error allocating dvdx'
|
|
|
1640 |
allocate(dtdy(ie,je,ke),STAT=stat)
|
|
|
1641 |
if (stat.ne.0) print*,'calc_qx: error allocating dtdy'
|
|
|
1642 |
|
|
|
1643 |
if (prelev) then
|
|
|
1644 |
do k=1,ke
|
|
|
1645 |
call ddh2m(uu(1,1,k),dudx(1,1,k),cl,'X',ie,je,1,vmin,vmax,mdv)
|
|
|
1646 |
call ddh2m(tt(1,1,k),dtdx(1,1,k),cl,'X',ie,je,1,vmin,vmax,mdv)
|
|
|
1647 |
call ddh2m(vv(1,1,k),dvdx(1,1,k),cl,'X',ie,je,1,vmin,vmax,mdv)
|
|
|
1648 |
call ddh2m(tt(1,1,k),dtdy(1,1,k),cl,'Y',ie,je,1,vmin,vmax,mdv)
|
|
|
1649 |
enddo
|
|
|
1650 |
else
|
|
|
1651 |
call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)
|
|
|
1652 |
call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)
|
|
|
1653 |
call ddh3(uu,dudx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
|
|
|
1654 |
call ddh3(tt,dtdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
|
|
|
1655 |
call ddh3(vv,dvdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
|
|
|
1656 |
call ddh3(tt,dtdy,sp,dspdx,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
|
|
|
1657 |
endif
|
|
|
1658 |
|
|
|
1659 |
var=-1.e9*9.80616/273.*(dudx*dtdx+dvdx*dtdy)
|
|
|
1660 |
if (prelev) then
|
|
|
1661 |
do i=1,ie
|
|
|
1662 |
do j=1,je
|
|
|
1663 |
do k=1,ke
|
|
|
1664 |
if ((dudx(i,j,k).eq.mdv).or.
|
|
|
1665 |
> (dtdx(i,j,k).eq.mdv).or.
|
|
|
1666 |
> (dvdx(i,j,k).eq.mdv).or.
|
|
|
1667 |
> (dtdy(i,j,k).eq.mdv)) then
|
|
|
1668 |
var(i,j,k)=mdv
|
|
|
1669 |
endif
|
|
|
1670 |
enddo
|
|
|
1671 |
enddo
|
|
|
1672 |
enddo
|
|
|
1673 |
endif
|
|
|
1674 |
|
|
|
1675 |
IF (ALLOCATED(dspdx)) DEALLOCATE(dspdx)
|
|
|
1676 |
IF (ALLOCATED(dspdy)) DEALLOCATE(dspdy)
|
|
|
1677 |
IF (ALLOCATED(dudx)) DEALLOCATE(dudx)
|
|
|
1678 |
IF (ALLOCATED(dtdx)) DEALLOCATE(dtdx)
|
|
|
1679 |
IF (ALLOCATED(dvdx)) DEALLOCATE(dvdx)
|
|
|
1680 |
IF (ALLOCATED(dtdy)) DEALLOCATE(dtdy)
|
|
|
1681 |
|
|
|
1682 |
end
|
|
|
1683 |
|
|
|
1684 |
|
|
|
1685 |
subroutine calc_qy(var,uu,vv,tt,sp,cl,tl,ie,je,ke,ak,bk,vmin,vmax)
|
|
|
1686 |
C ===============================================================
|
|
|
1687 |
c calculate y-comp Q-vector: -9.8/273.*(D[U:Y]*D[T:X]+D[V:Y]*D[T:Y]
|
|
|
1688 |
c -U/6.37E6*TAN*D[T:X]-V/6.37E6*TAN*D[T:Y])
|
|
|
1689 |
|
|
|
1690 |
c argument declaration
|
|
|
1691 |
integer ie,je,ke
|
|
|
1692 |
real,intent(OUT) :: var(ie,je,ke)
|
|
|
1693 |
real,intent(IN) :: uu(ie,je,ke),vv(ie,je,ke),tt(ie,je,ke)
|
|
|
1694 |
real,intent(IN) :: sp(ie,je),cl(ie,je),tl(ie,je)
|
|
|
1695 |
real,intent(IN) :: ak(ke),bk(ke),vmin(4),vmax(4)
|
|
|
1696 |
|
|
|
1697 |
c variable declaration
|
|
|
1698 |
REAL,ALLOCATABLE, DIMENSION (:,:) :: dspdx,dspdy
|
|
|
1699 |
REAL,ALLOCATABLE, DIMENSION (:,:,:) :: dudy,dtdx,dvdy,dtdy,tl3d
|
|
|
1700 |
integer k,stat
|
|
|
1701 |
real mdv
|
|
|
1702 |
logical prelev
|
|
|
1703 |
|
|
|
1704 |
c test if pressure or model levels
|
|
|
1705 |
prelev=.true.
|
|
|
1706 |
do k=1,ke
|
|
|
1707 |
if (bk(k).ne.0.) prelev=.false.
|
|
|
1708 |
enddo
|
|
|
1709 |
mdv=-999.98999
|
|
|
1710 |
|
|
|
1711 |
if (.not.prelev) then
|
|
|
1712 |
allocate(dspdx(ie,je),STAT=stat)
|
|
|
1713 |
if (stat.ne.0) print*,'calc_qy: error allocating dspdx'
|
|
|
1714 |
allocate(dspdy(ie,je),STAT=stat)
|
|
|
1715 |
if (stat.ne.0) print*,'calc_qy: error allocating dspdy'
|
|
|
1716 |
endif
|
|
|
1717 |
|
|
|
1718 |
allocate(dudy(ie,je,ke),STAT=stat)
|
|
|
1719 |
if (stat.ne.0) print*,'calc_qy: error allocating dudy'
|
|
|
1720 |
allocate(dtdx(ie,je,ke),STAT=stat)
|
|
|
1721 |
if (stat.ne.0) print*,'calc_qy: error allocating dtdx'
|
|
|
1722 |
allocate(dvdy(ie,je,ke),STAT=stat)
|
|
|
1723 |
if (stat.ne.0) print*,'calc_qy: error allocating dvdy'
|
|
|
1724 |
allocate(dtdy(ie,je,ke),STAT=stat)
|
|
|
1725 |
if (stat.ne.0) print*,'calc_qy: error allocating dtdy'
|
|
|
1726 |
allocate(tl3d(ie,je,ke),STAT=stat)
|
|
|
1727 |
if (stat.ne.0) print*,'calc_qy: error allocating tl3d'
|
|
|
1728 |
|
|
|
1729 |
if (prelev) then
|
|
|
1730 |
do k=1,ke
|
|
|
1731 |
call ddh2m(uu(1,1,k),dudy(1,1,k),cl,'Y',ie,je,1,vmin,vmax,mdv)
|
|
|
1732 |
call ddh2m(tt(1,1,k),dtdx(1,1,k),cl,'X',ie,je,1,vmin,vmax,mdv)
|
|
|
1733 |
call ddh2m(vv(1,1,k),dvdy(1,1,k),cl,'Y',ie,je,1,vmin,vmax,mdv)
|
|
|
1734 |
call ddh2m(tt(1,1,k),dtdy(1,1,k),cl,'Y',ie,je,1,vmin,vmax,mdv)
|
|
|
1735 |
enddo
|
|
|
1736 |
else
|
|
|
1737 |
call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)
|
|
|
1738 |
call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)
|
|
|
1739 |
call ddh3(uu,dudy,sp,dspdx,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
|
|
|
1740 |
call ddh3(tt,dtdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
|
|
|
1741 |
call ddh3(vv,dvdy,sp,dspdx,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
|
|
|
1742 |
call ddh3(tt,dtdy,sp,dspdx,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
|
|
|
1743 |
endif
|
|
|
1744 |
|
|
|
1745 |
|
|
|
1746 |
do k=1,ke
|
|
|
1747 |
tl3d(1:ie,1:je,k)=tl(1:ie,1:je)
|
|
|
1748 |
enddo
|
|
|
1749 |
|
|
|
1750 |
var=-1.e9*9.80616/273.*(dudy*dtdx+dvdy*dtdy
|
|
|
1751 |
> -uu/6.37E6*tl3d*dtdx-vv/6.37E6*tl3d*dtdy)
|
|
|
1752 |
if (prelev) then
|
|
|
1753 |
do i=1,ie
|
|
|
1754 |
do j=1,je
|
|
|
1755 |
do k=1,ke
|
|
|
1756 |
if ((dudy(i,j,k).eq.mdv).or.
|
|
|
1757 |
> (dtdx(i,j,k).eq.mdv).or.
|
|
|
1758 |
> (dvdy(i,j,k).eq.mdv).or.
|
|
|
1759 |
> (dtdy(i,j,k).eq.mdv)) then
|
|
|
1760 |
var(i,j,k)=mdv
|
|
|
1761 |
endif
|
|
|
1762 |
enddo
|
|
|
1763 |
enddo
|
|
|
1764 |
enddo
|
|
|
1765 |
endif
|
|
|
1766 |
|
|
|
1767 |
IF (ALLOCATED(dspdx)) DEALLOCATE(dspdx)
|
|
|
1768 |
IF (ALLOCATED(dspdy)) DEALLOCATE(dspdy)
|
|
|
1769 |
IF (ALLOCATED(dudy)) DEALLOCATE(dudy)
|
|
|
1770 |
IF (ALLOCATED(dtdx)) DEALLOCATE(dtdx)
|
|
|
1771 |
IF (ALLOCATED(dvdy)) DEALLOCATE(dvdy)
|
|
|
1772 |
IF (ALLOCATED(dtdy)) DEALLOCATE(dtdy)
|
|
|
1773 |
IF (ALLOCATED(tl3d)) DEALLOCATE(tl3d)
|
|
|
1774 |
|
|
|
1775 |
end
|
|
|
1776 |
|
|
|
1777 |
|
|
|
1778 |
subroutine calc_ug(var,zz,sp,cl,f,ie,je,ke,ak,bk,vmin,vmax)
|
|
|
1779 |
C ===========================================================
|
|
|
1780 |
C calculate geostrophic wind
|
|
|
1781 |
|
|
|
1782 |
c argument declaration
|
|
|
1783 |
integer ie,je,ke
|
|
|
1784 |
real,intent(OUT) :: var(ie,je,ke)
|
|
|
1785 |
real,intent(IN) :: zz(ie,je,ke),
|
|
|
1786 |
> sp(ie,je),cl(ie,je),f(ie,je)
|
|
|
1787 |
real ak(ke),bk(ke),vmin(4),vmax(4)
|
|
|
1788 |
|
|
|
1789 |
c variable declaration
|
|
|
1790 |
REAL,ALLOCATABLE, DIMENSION (:,:) :: dspdy
|
|
|
1791 |
REAL,ALLOCATABLE, DIMENSION (:,:,:) :: dzzdy
|
|
|
1792 |
integer k,stat
|
|
|
1793 |
|
|
|
1794 |
allocate(dspdy(ie,je),STAT=stat)
|
|
|
1795 |
if (stat.ne.0) print*,'calc_ug: error allocating dspdy'
|
|
|
1796 |
allocate(dzzdy(ie,je,ke),STAT=stat)
|
|
|
1797 |
if (stat.ne.0) print*,'calc_ug: error allocating dzzdy'
|
|
|
1798 |
|
|
|
1799 |
call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)
|
|
|
1800 |
call ddh3(zz,dzzdy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
|
|
|
1801 |
|
|
|
1802 |
do k=1,ke
|
|
|
1803 |
var(1:ie,1:je,k)=-dzzdy(1:ie,1:je,k)
|
|
|
1804 |
> *9810.*(f(1:ie,1:je)/((ABS(f(1:ie,1:je))+1.e-12)**2))
|
|
|
1805 |
enddo
|
|
|
1806 |
|
|
|
1807 |
IF (ALLOCATED(dspdy)) DEALLOCATE(dspdy)
|
|
|
1808 |
IF (ALLOCATED(dzzdy)) DEALLOCATE(dzzdy)
|
|
|
1809 |
|
|
|
1810 |
end
|
|
|
1811 |
|
|
|
1812 |
|
|
|
1813 |
subroutine calc_vg(var,zz,sp,cl,f,ie,je,ke,ak,bk,vmin,vmax)
|
|
|
1814 |
C ===========================================================
|
|
|
1815 |
C calculate geostrophic wind
|
|
|
1816 |
|
|
|
1817 |
c argument declaration
|
|
|
1818 |
integer ie,je,ke
|
|
|
1819 |
real,intent(OUT) :: var(ie,je,ke)
|
|
|
1820 |
real,intent(IN) :: zz(ie,je,ke),
|
|
|
1821 |
> sp(ie,je),cl(ie,je),f(ie,je)
|
|
|
1822 |
real ak(ke),bk(ke),vmin(4),vmax(4)
|
|
|
1823 |
|
|
|
1824 |
c variable declaration
|
|
|
1825 |
REAL,ALLOCATABLE, DIMENSION (:,:) :: dspdx
|
|
|
1826 |
REAL,ALLOCATABLE, DIMENSION (:,:,:) :: dzzdx
|
|
|
1827 |
integer k,stat
|
|
|
1828 |
|
|
|
1829 |
allocate(dspdx(ie,je),STAT=stat)
|
|
|
1830 |
if (stat.ne.0) print*,'calc_vg: error allocating dspdx'
|
|
|
1831 |
allocate(dzzdx(ie,je,ke),STAT=stat)
|
|
|
1832 |
if (stat.ne.0) print*,'calc_vg: error allocating dzzdx'
|
|
|
1833 |
|
|
|
1834 |
call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)
|
|
|
1835 |
call ddh3(zz,dzzdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
|
|
|
1836 |
|
|
|
1837 |
do k=1,ke
|
|
|
1838 |
var(1:ie,1:je,k)=dzzdx(1:ie,1:je,k)
|
|
|
1839 |
> *9810.*(f(1:ie,1:je)/((ABS(f(1:ie,1:je))+1.e-12)**2))
|
|
|
1840 |
enddo
|
|
|
1841 |
|
|
|
1842 |
IF (ALLOCATED(dspdx)) DEALLOCATE(dspdx)
|
|
|
1843 |
IF (ALLOCATED(dzzdx)) DEALLOCATE(dzzdx)
|
|
|
1844 |
|
|
|
1845 |
end
|
|
|
1846 |
|
|
|
1847 |
subroutine zlayer(ap,t,z,sp,ie,je,ke,aklev,bklev,aklay,bklay)
|
|
|
1848 |
c =============================================================
|
|
|
1849 |
c argument declaration
|
|
|
1850 |
integer ie,je,ke
|
|
|
1851 |
real,intent(OUT) :: ap(ie,je,ke)
|
|
|
1852 |
real,intent(IN) :: t(ie,je,ke),z(ie,je,ke),sp(ie,je)
|
|
|
1853 |
real,intent(IN) :: aklev(ke),bklev(ke),aklay(ke),bklay(ke)
|
|
|
1854 |
|
|
|
1855 |
c variable declaration
|
|
|
1856 |
integer i,j,k
|
|
|
1857 |
real psrf
|
|
|
1858 |
real,parameter :: rdg=29.271,tzero=273.15
|
|
|
1859 |
|
|
|
1860 |
c statement-functions for the computation of pressure
|
|
|
1861 |
real prlev,prlay
|
|
|
1862 |
integer is
|
|
|
1863 |
prlev(is)=aklev(is)+bklev(is)*psrf
|
|
|
1864 |
prlay(is)=aklay(is)+bklay(is)*psrf
|
|
|
1865 |
|
|
|
1866 |
c computation of height of main levels
|
|
|
1867 |
do i=1,ie
|
|
|
1868 |
do j=1,je
|
|
|
1869 |
psrf=sp(i,j)
|
|
|
1870 |
do k=1,ke-1
|
|
|
1871 |
ap(i,j,k) = (z(i,j,k) - rdg*(t(i,j,k)+tzero)*
|
|
|
1872 |
+ alog(prlay(k)/prlev(k)))/1000.
|
|
|
1873 |
enddo
|
|
|
1874 |
ap(i,j,ke) = (z(i,j,ke-1) + rdg*(t(i,j,ke)+tzero)*
|
|
|
1875 |
+ alog(prlev(ke-1)/prlay(ke)))/1000.
|
|
|
1876 |
enddo
|
|
|
1877 |
enddo
|
|
|
1878 |
end
|
|
|
1879 |
|
|
|
1880 |
|
|
|
1881 |
subroutine calc_psred(psr,ps,t,zb,ie,je,ke,aklay,bklay)
|
|
|
1882 |
c =======================================================
|
|
|
1883 |
c calculate PSRED given T, ZB and ak,bk
|
|
|
1884 |
|
|
|
1885 |
c argument declaration
|
|
|
1886 |
integer ie,je,ke
|
|
|
1887 |
real,intent(OUT) :: psr(ie,je)
|
|
|
1888 |
real,intent(IN) :: ps(ie,je),t(ie,je,ke),zb(ie,je)
|
|
|
1889 |
real,intent(IN) :: aklay(ke),bklay(ke)
|
|
|
1890 |
|
|
|
1891 |
c variable declaration
|
|
|
1892 |
integer i,j
|
|
|
1893 |
real psrf
|
|
|
1894 |
real ztstar,zalpha,zt0
|
|
|
1895 |
real,parameter :: rdcp=0.286,tzero=273.15,r=287.05,g=9.80665
|
|
|
1896 |
|
|
|
1897 |
c statement-functions for the computation of pressure
|
|
|
1898 |
real prlay
|
|
|
1899 |
integer is
|
|
|
1900 |
prlay(is)=aklay(is)+bklay(is)*psrf
|
|
|
1901 |
|
|
|
1902 |
do i=1,ie
|
|
|
1903 |
do j=1,je
|
|
|
1904 |
psrf=ps(i,j)
|
|
|
1905 |
if (zb(i,j).lt.1.) then
|
|
|
1906 |
psr(i,j)=psrf
|
|
|
1907 |
else
|
|
|
1908 |
ztstar = (t(i,j,1)+tzero)*(1. + 0.0065*r/g*
|
|
|
1909 |
& (ps(i,j)/prlay(1) - 1.0))
|
|
|
1910 |
zalpha = 0.0065*r
|
|
|
1911 |
zt0 = ztstar + 0.0065*zb(i,j)
|
|
|
1912 |
if (zt0.gt.290.5) then
|
|
|
1913 |
if (ztstar.gt.290.5) then
|
|
|
1914 |
zalpha = 0.0
|
|
|
1915 |
ztstar = 0.5*(ztstar+290.5)
|
|
|
1916 |
else
|
|
|
1917 |
zalpha = r*(290.5-ztstar)/zb(i,j)
|
|
|
1918 |
endif
|
|
|
1919 |
else if (ztstar.lt.255.) then
|
|
|
1920 |
ztstar = 0.5*(255.0+ztstar)
|
|
|
1921 |
endif
|
|
|
1922 |
psr(i,j) = ps(i,j)* exp(g*zb(i,j)/(r*ztstar)*
|
|
|
1923 |
& (1.0 - 0.5*(zalpha*zb(i,j)/(r*ztstar)) +
|
|
|
1924 |
& 0.333*(zalpha*zb(i,j)/(r*ztstar))**2))
|
|
|
1925 |
endif
|
|
|
1926 |
enddo
|
|
|
1927 |
enddo
|
|
|
1928 |
end
|
|
|
1929 |
|
|
|
1930 |
subroutine calc_blh(blh,ps,t,q,u,v,t2m,td2m,u10m,v10m,
|
|
|
1931 |
> ie,je,ke,aklay,bklay)
|
|
|
1932 |
c ======================================================
|
|
|
1933 |
c calculate BLH with routines from Andi Stohl
|
|
|
1934 |
|
|
|
1935 |
c argument declaration
|
|
|
1936 |
integer ie,je,ke
|
|
|
1937 |
real,intent(OUT) :: blh(ie,je)
|
|
|
1938 |
real,intent(IN) :: ps(ie,je),t2m(ie,je),td2m(ie,je),
|
|
|
1939 |
> u10m(ie,je),v10m(ie,je)
|
|
|
1940 |
real,intent(IN) :: t(ie,je,ke),q(ie,je,ke),u(ie,je,ke),
|
|
|
1941 |
> v(ie,je,ke)
|
|
|
1942 |
real,intent(IN) :: aklay(ke),bklay(ke)
|
|
|
1943 |
|
|
|
1944 |
c variable declaration
|
|
|
1945 |
integer i,j
|
|
|
1946 |
real psrf
|
|
|
1947 |
real ztstar,zalpha,zt0
|
|
|
1948 |
real,parameter :: rdcp=0.286,tzero=273.15,r=287.05,g=9.80665
|
|
|
1949 |
|
|
|
1950 |
c statement-functions for the computation of pressure
|
|
|
1951 |
real prlay
|
|
|
1952 |
integer is
|
|
|
1953 |
prlay(is)=aklay(is)+bklay(is)*psrf
|
|
|
1954 |
|
|
|
1955 |
do i=1,ie
|
|
|
1956 |
do j=1,je
|
|
|
1957 |
call richardson(ps,ust,t(i,j,1),q(i,j,1),u(i,j,1),v(i,j,1),
|
|
|
1958 |
> ke,aklay,bklay,hf,t2m,td2m,blh(i,j),wst)
|
|
|
1959 |
enddo
|
|
|
1960 |
enddo
|
|
|
1961 |
end
|
|
|
1962 |
|
|
|
1963 |
subroutine richardson(psurf,ust,ttlev,qvlev,ulev,vlev,nuvz,
|
|
|
1964 |
+akz,bkz,hf,tt2,td2,h,wst)
|
|
|
1965 |
*****************************************************************************
|
|
|
1966 |
* *
|
|
|
1967 |
* Calculation of mixing height based on the critical Richardson number. *
|
|
|
1968 |
* Calculation of convective time scale. *
|
|
|
1969 |
* For unstable conditions, one iteration is performed. An excess *
|
|
|
1970 |
* temperature (dependent on hf and wst) is calculated, added to the *
|
|
|
1971 |
* temperature at the lowest model level. Then the procedure is repeated.*
|
|
|
1972 |
* *
|
|
|
1973 |
* Author: A. Stohl *
|
|
|
1974 |
* *
|
|
|
1975 |
* 22 August 1996 *
|
|
|
1976 |
* *
|
|
|
1977 |
* Literature: *
|
|
|
1978 |
* Vogelezang DHP and Holtslag AAM (1996): Evaluation and model impacts *
|
|
|
1979 |
* of alternative boundary-layer height formulations. Boundary-Layer *
|
|
|
1980 |
* Meteor. 81, 245-269. *
|
|
|
1981 |
* *
|
|
|
1982 |
* Update: 1999-02-01 by G. Wotawa *
|
|
|
1983 |
* *
|
|
|
1984 |
* Two meter level (temperature, humidity) is taken as reference level *
|
|
|
1985 |
* instead of first model level. *
|
|
|
1986 |
* New input variables tt2, td2 introduced. *
|
|
|
1987 |
* *
|
|
|
1988 |
*****************************************************************************
|
|
|
1989 |
* *
|
|
|
1990 |
* Variables: *
|
|
|
1991 |
* h mixing height [m] *
|
|
|
1992 |
* hf sensible heat flux *
|
|
|
1993 |
* psurf surface pressure at point (xt,yt) [Pa] *
|
|
|
1994 |
* tv virtual temperature *
|
|
|
1995 |
* wst convective velocity scale *
|
|
|
1996 |
* *
|
|
|
1997 |
* Constants: *
|
|
|
1998 |
* ric critical Richardson number *
|
|
|
1999 |
* *
|
|
|
2000 |
*****************************************************************************
|
|
|
2001 |
|
|
|
2002 |
* include 'includepar'
|
|
|
2003 |
|
|
|
2004 |
integer i,k,nuvz,itmax,iter
|
|
|
2005 |
real tv,tvold,zref,z,zold,pint,pold,theta,thetaref,wm,ri,riold
|
|
|
2006 |
real akz(nuvz),bkz(nuvz),ulev(nuvz),vlev(nuvz),hf,wst,tt2,td2,ew
|
|
|
2007 |
real psurf,ust,ttlev(nuvz),qvlev(nuvz),h,const,ric,b,excess,bs
|
|
|
2008 |
real thetaold,zl,ul,vl,thetal,ril
|
|
|
2009 |
parameter(r_air=287.05,ga=9.81)
|
|
|
2010 |
parameter(const=r_air/ga,ric=0.25,b=100.,bs=8.5,itmax=3)
|
|
|
2011 |
|
|
|
2012 |
|
|
|
2013 |
excess=0.0
|
|
|
2014 |
iter=0
|
|
|
2015 |
|
|
|
2016 |
C Compute virtual temperature and virtual potential temperature at
|
|
|
2017 |
C reference level (2 m)
|
|
|
2018 |
******************************************************************
|
|
|
2019 |
|
|
|
2020 |
30 iter=iter+1
|
|
|
2021 |
|
|
|
2022 |
pold=psurf
|
|
|
2023 |
tvold=tt2*(1.+0.378*esat(td2)/psurf)
|
|
|
2024 |
zold=2.0
|
|
|
2025 |
zref=zold
|
|
|
2026 |
|
|
|
2027 |
thetaref=tvold*(100000./pold)**(r_air/cpa)+excess
|
|
|
2028 |
riold=0.
|
|
|
2029 |
|
|
|
2030 |
|
|
|
2031 |
C Integrate z up to one level above zt
|
|
|
2032 |
**************************************
|
|
|
2033 |
|
|
|
2034 |
do 10 k=3,nuvz
|
|
|
2035 |
pint=akz(k)+bkz(k)*psurf ! pressure on model layers
|
|
|
2036 |
wm=qvlev(k)/(1.-qvlev(k))
|
|
|
2037 |
tv=ttlev(k)*(1.+0.378*wm/(wm+.622))
|
|
|
2038 |
|
|
|
2039 |
if (abs(tv-tvold).gt.0.2) then
|
|
|
2040 |
z=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
|
|
|
2041 |
else
|
|
|
2042 |
z=zold+const*log(pold/pint)*tv
|
|
|
2043 |
endif
|
|
|
2044 |
|
|
|
2045 |
theta=tv*(100000./pint)**(r_air/cpa)
|
|
|
2046 |
|
|
|
2047 |
|
|
|
2048 |
Calculate Richardson number at each level
|
|
|
2049 |
*****************************************
|
|
|
2050 |
|
|
|
2051 |
ri=ga/thetaref*(theta-thetaref)*(z-zref)/
|
|
|
2052 |
+ max(((ulev(k)-ulev(2))**2+(vlev(k)-vlev(2))**2+b*ust**2),0.1)
|
|
|
2053 |
|
|
|
2054 |
if (ri.gt.ric) goto 20
|
|
|
2055 |
|
|
|
2056 |
riold=ri
|
|
|
2057 |
tvold=tv
|
|
|
2058 |
pold=pint
|
|
|
2059 |
thetaold=theta
|
|
|
2060 |
10 zold=z
|
|
|
2061 |
|
|
|
2062 |
20 continue
|
|
|
2063 |
|
|
|
2064 |
C Determine Richardson number between the critical levels
|
|
|
2065 |
*********************************************************
|
|
|
2066 |
|
|
|
2067 |
do 15 i=1,20
|
|
|
2068 |
zl=zold+float(i)/20.*(z-zold)
|
|
|
2069 |
ul=ulev(k-1)+float(i)/20.*(ulev(k)-ulev(k-1))
|
|
|
2070 |
vl=vlev(k-1)+float(i)/20.*(vlev(k)-vlev(k-1))
|
|
|
2071 |
vl=vlev(k-1)+float(i)/20.*(vlev(k)-vlev(k-1))
|
|
|
2072 |
thetal=thetaold+float(i)/20.*(theta-thetaold)
|
|
|
2073 |
ril=ga/thetaref*(thetal-thetaref)*(zl-zref)/
|
|
|
2074 |
+ max(((ul-ulev(2))**2+(vl-vlev(2))**2+b*ust**2),0.1)
|
|
|
2075 |
if (ril.gt.ric) goto 25
|
|
|
2076 |
15 continue
|
|
|
2077 |
|
|
|
2078 |
25 continue
|
|
|
2079 |
h=zl
|
|
|
2080 |
|
|
|
2081 |
|
|
|
2082 |
C Calculate convective velocity scale
|
|
|
2083 |
*************************************
|
|
|
2084 |
|
|
|
2085 |
if (hf.lt.0.) then
|
|
|
2086 |
wst=(-h*ga/thetaref*hf/cpa)**0.333
|
|
|
2087 |
excess=-bs*hf/cpa/wst
|
|
|
2088 |
if (iter.lt.itmax) goto 30
|
|
|
2089 |
else
|
|
|
2090 |
wst=0.
|
|
|
2091 |
endif
|
|
|
2092 |
|
|
|
2093 |
return
|
|
|
2094 |
end
|