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