UPP  V11.0.0
 All Data Structures Files Functions Pages
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
grib2_module
Definition: grib2_module.f:1
upp_physics
calcape() computes CAPE/CINS and other storm related variables.
Definition: UPP_PHYSICS.f:27