UPP  11.0.0
 All Data Structures Files Functions Variables Pages
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
ctlblk_mod
Definition: CTLBLK.f:1
rqstfld_mod
Definition: RQSTFLD.f:1
fdlvl_uv
subroutine fdlvl_uv(ITYPE, NFD, HTFD, UFD, VFD)
Computes FD level for u,v.
Definition: FDLVL.f:453
grib2_module
Definition: grib2_module.f:1
mdl2std_p
subroutine mdl2std_p()
mdl2std_p() vertical interpolation of model levels to standard atmospheric pressure.
Definition: MDL2STD_P.f:22
relabel
real function relabel(p)
relabel() relabels the pressure level to reference (or standard atmospheric) pressure levels rather t...
Definition: MDL2STD_P.f:423
p2h
real function p2h(p)
P2H() converts pressure levels (hPa) to geopotential heights.
Definition: MDL2STD_P.f:399
fdlvl_mass
subroutine fdlvl_mass(ITYPE, NFD, PTFD, HTFD, NIN, QIN, QTYPE, QFD)
Computes FD level for mass variables.
Definition: FDLVL.f:784