UPP (upp-srw-2.2.0)
Loading...
Searching...
No Matches
MDL2STD_P.f
Go to the documentation of this file.
1
19!--------------------------------------------------------------------------------------
22 SUBROUTINE mdl2std_p()
23
24!
25 use vrbls3d, only: pint, pmid, zmid
26 use vrbls3d, only: t, q, uh, vh, omga, cwm, qqw, qqi, qqr, qqs, qqg
27
28 use vrbls3d, only: icing_gfip, icing_gfis, catedr, mwt, gtg
29 use ctlblk_mod, only: grib, cfld, fld_info, datapd, im, jsta, jend, jm, &
30 lm, htfd, spval, nfd, me,&
31 jsta_2l, jend_2u, modelname,&
32 ista, iend, ista_2l, iend_2u
33 use rqstfld_mod, only: iget, lvls, iavblfld, lvlsxml
34 use grib2_module, only: pset
35
36!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
37!
38 implicit none
39
40 real, external :: P2H, relabel
41
42 real,dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: grid1
43 real,dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: EGRID1,EGRID2,EGRID3,EGRID4
44
45!
46 integer I,J,ii,jj,L,ITYPE,IFD,ITYPEFDLVL(NFD)
47
48! Variables introduced to allow FD levels from control file - Y Mao
49 integer :: N,NFDCTL
50 REAL, allocatable :: HTFDCTL(:)
51 integer, allocatable :: ITYPEFDLVLCTL(:)
52 real, allocatable :: QIN(:,:,:,:), QFD(:,:,:,:)
53 character, allocatable :: QTYPE(:)
54 real, allocatable :: VAR3D1(:,:,:), VAR3D2(:,:,:)
55
56 integer, parameter :: NFDMAX=50 ! Max number of fields with the same HTFDCTL
57 integer :: IDS(NFDMAX) ! All field IDs with the same HTFDCTL
58 integer :: nFDS ! How many fields with the same HTFDCTL in the control file
59 integer :: iID ! which field with HTFDCTL
60 integer :: N1, N2
61!
62!******************************************************************************
63!
64! START MDL2STD_P.
65!
66
67! --------------WAFS block----------------------
68! 479 ICSEV
69! 481 ICIP
70! 476 EDPARM
71! 477 CAT
72! 478 MWTURB
73! 518 HGT
74! 519 TMP
75! 520 UGRD
76! 521 VGRD
77 IF(iget(479)>0 .or. iget(481)>0 .or. &
78 iget(476)>0 .or. iget(477)>0 .or. iget(478)>0 .or. &
79 iget(518)>0 .or. iget(519)>0 .or. iget(520)>0 .or. &
80 iget(521)>0) then
81
82! STEP 1 -- U V (POSSIBLE FOR ABSV) INTERPLOCATION
83 IF(iget(520)>0 .or. iget(521)>0 ) THEN
84! U/V are always paired, use any for HTFDCTL
85 iid=520
86 n = iavblfld(iget(iid))
87 nfdctl=size(pset%param(n)%level)
88 if(allocated(itypefdlvlctl)) deallocate(itypefdlvlctl)
89 allocate(itypefdlvlctl(nfdctl))
90 DO ifd = 1,nfdctl
91 itypefdlvlctl(ifd)=lvls(ifd,iget(iid))
92 ENDDO
93 if(allocated(htfdctl)) deallocate(htfdctl)
94 allocate(htfdctl(nfdctl))
95 htfdctl=pset%param(n)%level
96 DO i = 1, nfdctl
97 htfdctl(i)=p2h(htfdctl(i)/100.)
98 ENDDO
99 if(allocated(var3d1)) deallocate(var3d1)
100 if(allocated(var3d2)) deallocate(var3d2)
101 allocate(var3d1(ista_2l:iend_2u,jsta_2l:jend_2u,nfdctl))
102 allocate(var3d2(ista_2l:iend_2u,jsta_2l:jend_2u,nfdctl))
103 var3d1=spval
104 var3d2=spval
105
106 call fdlvl_uv(itypefdlvlctl,nfdctl,htfdctl,var3d1,var3d2)
107
108 DO ifd = 1,nfdctl
109 ! U
110 IF (lvls(ifd,iget(520)) > 0) THEN
111!$omp parallel do private(i,j)
112 DO j=jsta,jend
113 DO i=ista,iend
114 grid1(i,j)=var3d1(i,j,ifd)
115 ENDDO
116 ENDDO
117 if(grib=='grib2') then
118 cfld=cfld+1
119 fld_info(cfld)%ifld=iavblfld(iget(520))
120 fld_info(cfld)%lvl=lvlsxml(ifd,iget(520))
121!$omp parallel do private(i,j,ii,jj)
122 do j=1,jend-jsta+1
123 jj = jsta+j-1
124 do i=1,iend-ista+1
125 ii = ista+i-1
126 datapd(i,j,cfld) = grid1(ii,jj)
127 enddo
128 enddo
129 endif
130 ENDIF
131 ! V
132 IF (lvls(ifd,iget(521)) > 0) THEN
133!$omp parallel do private(i,j)
134 DO j=jsta,jend
135 DO i=ista,iend
136 grid1(i,j)=var3d2(i,j,ifd)
137 ENDDO
138 ENDDO
139 if(grib=='grib2') then
140 cfld=cfld+1
141 fld_info(cfld)%ifld=iavblfld(iget(521))
142 fld_info(cfld)%lvl=lvlsxml(ifd,iget(521))
143!$omp parallel do private(i,j,ii,jj)
144 do j=1,jend-jsta+1
145 jj = jsta+j-1
146 do i=1,iend-ista+1
147 ii = ista+i-1
148 datapd(i,j,cfld) = grid1(ii,jj)
149 enddo
150 enddo
151 endif
152 ENDIF
153 ENDDO
154
155 deallocate(var3d1)
156 deallocate(var3d2)
157
158 ENDIF
159
160! STEP 2 -- MASS FIELDS INTERPOLATION EXCEPT:
161! HGT(TO BE FIXED VALUES)
162! RH ABSV (TO BE CACULATED)
163
164 if(allocated(qin)) deallocate(qin)
165 if(allocated(qtype)) deallocate(qtype)
166 ALLOCATE(qin(ista:iend,jsta:jend,lm,nfdmax))
167 ALLOCATE(qtype(nfdmax))
168
169! INITIALIZE INPUTS
170 nfds = 0
171 IF(iget(479) > 0) THEN
172 nfds = nfds + 1
173 ids(nfds) = 479
174 qin(ista:iend,jsta:jend,1:lm,nfds)=icing_gfip(ista:iend,jsta:jend,1:lm)
175 qtype(nfds)="O"
176 end if
177 IF(iget(481) > 0) THEN
178 nfds = nfds + 1
179 ids(nfds) = 481
180 qin(ista:iend,jsta:jend,1:lm,nfds)=icing_gfis(ista:iend,jsta:jend,1:lm)
181 qtype(nfds)="O"
182 end if
183 IF(iget(476) > 0) THEN
184 nfds = nfds + 1
185 ids(nfds) = 476
186 qin(ista:iend,jsta:jend,1:lm,nfds)=gtg(ista:iend,jsta:jend,1:lm)
187 qtype(nfds)="O"
188 end if
189 IF(iget(477) > 0) THEN
190 nfds = nfds + 1
191 ids(nfds) = 477
192 qin(ista:iend,jsta:jend,1:lm,nfds)=catedr(ista:iend,jsta:jend,1:lm)
193 qtype(nfds)="O"
194 end if
195 IF(iget(478) > 0) THEN
196 nfds = nfds + 1
197 ids(nfds) = 478
198 qin(ista:iend,jsta:jend,1:lm,nfds)=mwt(ista:iend,jsta:jend,1:lm)
199 qtype(nfds)="O"
200 end if
201 IF(iget(519) > 0) THEN
202 nfds = nfds + 1
203 ids(nfds) = 519
204 qin(ista:iend,jsta:jend,1:lm,nfds)=t(ista:iend,jsta:jend,1:lm)
205 qtype(nfds)="T"
206 end if
207
208! FOR WAFS, ALL LEVLES OF DIFFERENT VARIABLES ARE THE SAME, USE ANY
209 iid=ids(1)
210 n = iavblfld(iget(iid))
211 nfdctl=size(pset%param(n)%level)
212 if(allocated(itypefdlvlctl)) deallocate(itypefdlvlctl)
213 allocate(itypefdlvlctl(nfdctl))
214 DO ifd = 1,nfdctl
215 itypefdlvlctl(ifd)=lvls(ifd,iget(iid))
216 ENDDO
217 if(allocated(htfdctl)) deallocate(htfdctl)
218 allocate(htfdctl(nfdctl))
219 htfdctl=pset%param(n)%level
220 DO i = 1, nfdctl
221 htfdctl(i)=p2h(htfdctl(i)/100.)
222 ENDDO
223
224 if(allocated(qfd)) deallocate(qfd)
225 ALLOCATE(qfd(ista:iend,jsta:jend,nfdctl,nfds))
226 qfd=spval
227
228 call fdlvl_mass(itypefdlvlctl,nfdctl,pset%param(n)%level,htfdctl,nfds,qin,qtype,qfd)
229
230! Adjust values before output
231 n1 = -1
232 DO n=1,nfds
233 iid=ids(n)
234
235! Icing Potential
236 if(iid==481) then
237 n1=n
238 DO ifd = 1,nfdctl
239 DO j=jsta,jend
240 DO i=ista,iend
241 if(qfd(i,j,ifd,n) < spval) then
242 qfd(i,j,ifd,n)=max(0.0,qfd(i,j,ifd,n))
243 qfd(i,j,ifd,n)=min(1.0,qfd(i,j,ifd,n))
244 endif
245 ENDDO
246 ENDDO
247 ENDDO
248 endif
249
250! Icing severity categories
251! 0 = none (0, 0.08)
252! 1 = trace [0.08, 0.21]
253! 2 = light (0.21, 0.37]
254! 3 = moderate (0.37, 0.67]
255! 4 = severe (0.67, 1]
256! http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-207.shtml
257 if(iid==479) then
258 DO ifd = 1,nfdctl
259 DO j=jsta,jend
260 DO i=ista,iend
261 if(n1 > 0) then
262 ! Icing severity is 0 when icing potential is too small
263 if(qfd(i,j,ifd,n1) < 0.001) qfd(i,j,ifd,n)=0.
264 endif
265 if(qfd(i,j,ifd,n) == spval) cycle
266 if (qfd(i,j,ifd,n) < 0.08) then
267 qfd(i,j,ifd,n) = 0.0
268 elseif (qfd(i,j,ifd,n) <= 0.21) then
269 qfd(i,j,ifd,n) = 1.
270 else if(qfd(i,j,ifd,n) <= 0.37) then
271 qfd(i,j,ifd,n) = 2.0
272 else if(qfd(i,j,ifd,n) <= 0.67) then
273 qfd(i,j,ifd,n) = 3.0
274 else
275 qfd(i,j,ifd,n) = 4.0
276 endif
277 ENDDO
278 ENDDO
279 ENDDO
280 endif
281
282! GTG turbulence: EDRPARM, CAT, MWTURB
283 if(iid==476 .or. iid==477 .or. iid==478) then
284 DO ifd = 1,nfdctl
285 DO j=jsta,jend
286 DO i=ista,iend
287 if(qfd(i,j,ifd,n) < spval) then
288 qfd(i,j,ifd,n)=max(0.0,qfd(i,j,ifd,n))
289 qfd(i,j,ifd,n)=min(1.0,qfd(i,j,ifd,n))
290 endif
291 ENDDO
292 ENDDO
293 ENDDO
294 endif
295
296 ENDDO
297
298! Output
299 DO n=1,nfds
300 iid=ids(n)
301 DO ifd = 1,nfdctl
302 IF (lvls(ifd,iget(iid)) > 0) THEN
303!$omp parallel do private(i,j)
304 DO j=jsta,jend
305 DO i=ista,iend
306 grid1(i,j)=qfd(i,j,ifd,n)
307 ENDDO
308 ENDDO
309 if(grib=='grib2') then
310 cfld=cfld+1
311 fld_info(cfld)%ifld=iavblfld(iget(iid))
312 fld_info(cfld)%lvl=lvlsxml(ifd,iget(iid))
313!$omp parallel do private(i,j,ii,jj)
314 do j=1,jend-jsta+1
315 jj = jsta+j-1
316 do i=1,iend-ista+1
317 ii = ista+i-1
318 datapd(i,j,cfld) = grid1(ii,jj)
319 enddo
320 enddo
321 endif
322 ENDIF
323 ENDDO
324 ENDDO
325
326 DEALLOCATE(qin,qfd)
327 DEALLOCATE(qtype)
328
329! STEP 3 - MASS FIELDS CALCULATION
330! HGT(TO BE FIXED VALUES)
331! RH ABSV (TO BE CACULATED)
332 ! HGT
333 IF(iget(518) > 0) THEN
334 iid=518
335 n = iavblfld(iget(iid))
336 nfdctl=size(pset%param(n)%level)
337 if(allocated(htfdctl)) deallocate(htfdctl)
338 allocate(htfdctl(nfdctl))
339 htfdctl=pset%param(n)%level
340 DO i = 1, nfdctl
341 htfdctl(i)=p2h(htfdctl(i)/100.)
342 ENDDO
343
344 DO ifd = 1,nfdctl
345 IF (lvls(ifd,iget(iid)) > 0) THEN
346!$omp parallel do private(i,j)
347 DO j=jsta,jend
348 DO i=ista,iend
349 grid1(i,j)=htfdctl(ifd)
350 ENDDO
351 ENDDO
352 if(grib=='grib2') then
353 cfld=cfld+1
354 fld_info(cfld)%ifld=iavblfld(iget(iid))
355 fld_info(cfld)%lvl=lvlsxml(ifd,iget(iid))
356!$omp parallel do private(i,j,ii,jj)
357 do j=1,jend-jsta+1
358 jj = jsta+j-1
359 do i=1,iend-ista+1
360 ii = ista+i-1
361 datapd(i,j,cfld) = grid1(ii,jj)
362 enddo
363 enddo
364 endif
365 ENDIF
366 ENDDO
367 ENDIF
368
369 ! Relabel the pressure level to reference levels
370! IDS = 0
371 ids = (/ 481,479,476,477,478,518,519,520,521,(0,i=10,50) /)
372 do i = 1, nfdmax
373 iid=ids(i)
374 if(iid == 0) exit
375 n = iavblfld(iget(iid))
376 nfdctl=size(pset%param(n)%level)
377 do j = 1, nfdctl
378 pset%param(n)%level(j) = relabel(pset%param(n)%level(j))
379 end do
380 end do
381
382 ENDIF
383
384!
385! END OF ROUTINE.
386!
387 RETURN
388 END
389
390!--------------------------------------------------------------------------------------
398
399 FUNCTION p2h(p)
400 implicit none
401 real, intent(in) :: p
402 real :: p2h
403 real, parameter :: lapse = 0.0065
404 real, parameter :: surf_temp = 288.15
405 real, parameter :: gravity = 9.80665
406 real, parameter :: moles_dry_air = 0.02896442
407 real, parameter :: gas_const = 8.31432
408 real, parameter :: surf_pres = 1013.25
409 real, parameter :: power_const = (gravity * moles_dry_air) &
410 / (gas_const * lapse)
411
412 p2h = (surf_temp/lapse)*(1-(p/surf_pres)**(1/power_const))
413 END
414
415!--------------------------------------------------------------------------------------
422
423 function relabel(p)
424 implicit none
425 real, intent(in) :: p
426 real :: relabel
427 relabel=p
428 if(p == 10040.) relabel=10000
429 if(p == 12770.) relabel=12500
430 if(p == 14750.) relabel=15000
431 if(p == 17870.) relabel=17500
432 if(p == 19680.) relabel=20000
433 if(p == 22730.) relabel=22500
434 if(p == 27450.) relabel=27500
435 if(p == 30090.) relabel=30000
436 if(p == 34430.) relabel=35000
437 if(p == 39270.) relabel=40000
438 if(p == 44650.) relabel=45000
439 if(p == 50600.) relabel=50000
440 if(p == 59520.) relabel=60000
441 if(p == 69680.) relabel=70000
442 if(p == 75260.) relabel=75000
443 if(p == 81200.) relabel=80000
444 if(p == 84310.) relabel=85000
445 END
fdlvl_uv
subroutine fdlvl_uv(itype, nfd, htfd, ufd, vfd)
Computes FD level for u,v.
Definition FDLVL.f:454
fdlvl_mass
subroutine fdlvl_mass(itype, nfd, ptfd, htfd, nin, qin, qtype, qfd)
Computes FD level for mass variables.
Definition FDLVL.f:785
relabel
real function relabel(p)
relabel() relabels the pressure level to reference (or standard atmospheric) pressure levels rather t...
Definition MDL2STD_P.f:424
mdl2std_p
subroutine mdl2std_p()
mdl2std_p() vertical interpolation of model levels to standard atmospheric pressure.
Definition MDL2STD_P.f:23
p2h
real function p2h(p)
P2H() converts pressure levels (hPa) to geopotential heights.
Definition MDL2STD_P.f:400