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