UPP  001
 All Data Structures Files Functions Pages
GFIP3.f
1 !========================================================================
2 ! = = = = = = = = = = = = module DerivedFields = = = = = = = = = = = =
3 !========================================================================
4 module derivedfields
5 
6  implicit none
7 
8  private
9  public derive_fields, mixing_ratio
10  public precips
11 
12  ! Precipitation types
13  type :: precipitations_t
14  integer :: none = 0
15  integer :: rain = 1
16  integer :: snow = 2
17  integer :: other = 3
18  integer :: convection = 4
19  end type precipitations_t
20  type(precipitations_t), parameter :: precips = precipitations_t()
21 
22 contains
23 
24 !-----------------------------------------------------------------------+
25  subroutine derive_fields(imp_physics,t, rh, pres, hgt, totalWater, totalCond,&
26  nz, topok, hprcp, hcprcp, cin, cape, &
27  ept, wbt, twp, pc, kx, lx, tott, prcptype)
28  IMPLICIT NONE
29 ! 3-D derived data:
30 ! ept - equivalent potential temperature
31 ! wbt - wet bulb temperature
32 ! twp - total water path
33 !
34 ! 2-D derived data: (indice for convective icing severity)
35 ! kx - k index
36 ! lx - lifted index
37 ! tott - total totals
38 !
39 ! 2-D derived data:
40 ! pc - precipitation condensate
41 ! prcpType - surface precipitation type
42 
43  integer, intent(in) :: imp_physics
44  integer, intent(in) :: nz, topok
45  real, intent(in),dimension(nz) :: t, rh, pres, hgt
46  real, intent(in),dimension(nz) :: totalwater,totalcond
47  real, intent(in) :: hprcp, hcprcp, cin, cape
48  real, intent(out) :: ept(nz), wbt(nz), twp(nz)
49  real, intent(out) :: pc, kx, lx, tott
50  integer, intent(out) :: prcptype
51 
52  real, allocatable :: td(:), tlc(:), wvm(:)
53  integer :: k
54 
55  allocate(td(nz))
56  allocate(tlc(nz))
57  allocate(wvm(nz))
58 
59  td(:) = getdewpointtemp(t(:), rh(:))
60  tlc(:) = get_tlcl(t(:), td(:))
61  wvm(:) = mixing_ratio(td(:), pres(:))
62 
63  ept(:) = getthetae(pres(:), t(:), wvm(:), tlc(:))
64  wbt(:) = getthetaw(t(:), td(:), pres(:))
65 
66  call calc_totalwaterpath(hgt, pres, t, totalwater, nz, twp)
67 
68 
69  pc = getprecipcond(totalcond, nz, topok)
70 
71  ! indice for convective icing severity
72  call calc_indice(t, td, pres, wvm, nz, topok, kx, lx, tott)
73 
74  prcptype=getpreciptype(imp_physics,pres,hgt,t,rh,wbt,nz,&
75  hprcp,hcprcp,pc,cin,cape,lx)
76 
77  deallocate(td)
78  deallocate(tlc)
79  deallocate(wvm)
80 
81  return
82  end subroutine derive_fields
83 
84 !-----------------------------------------------------------------------+
85 ! dew point temperature in K
86 ! t in K
87  elemental real function getdewpointtemp(t, rh)
88  IMPLICIT NONE
89  real, intent(in) :: t, rh
90 
91  real vapr, rm
92 
93  ! actual water vapor presure in pascal
94  ! 611.2 Pa is the saturated vapor presure at 0°C
95  ! Pa = (rh/100) * Ps
96  vapr = (max(1.0e-6,rh) / 100.0) * getvaporpres(t)
97  rm = log(vapr/611.2)
98 
99  getdewpointtemp = 243.5 * rm / (17.67 - rm) + 273.15
100 
101  return
102  end function getdewpointtemp
103 
104 !-----------------------------------------------------------------------+
105 ! Equivalent potential temperature in K
106 ! pres in Pa
107 ! t in K
108 ! wvm in kg/kg
109  elemental real function getthetae(pres, t, wvm, tAtLCL)
110  IMPLICIT NONE
111  real, intent(in) :: pres, t, wvm, tatlcl
112 
113  real rmix, e, thtam
114 
115  rmix = max(0.0,wvm) ! in g/kg
116  ! potential temperature
117  e = (2.0/7.0) * ( 1.0 - 0.28 * rmix * 0.001)
118  thtam = t * ((100000.0/pres)**e)
119  getthetae = thtam * exp( (3.376/tatlcl - 0.00254) * &
120  (rmix * (1. + (.81 * rmix * 0.001)) ))
121 
122  return
123  end function getthetae
124 
125 !-----------------------------------------------------------------------+
126 ! saturation vapor presure in Pa
127 ! t in K
128  elemental real function getvaporpres(t)
129  IMPLICIT NONE
130  real, intent(in) :: t
131 
132  real tc
133 
134  ! 611.2 Pa is the saturated vapor presure at 0°C
135  tc = t-273.15
136  getvaporpres = 611.2 * exp( 17.67*tc/(tc+243.5))
137 
138  return
139  end function getvaporpres
140 
141 !-----------------------------------------------------------------------+
142 ! lifted condensation level temperature in K
143 ! t and td in K
144  elemental real function get_tlcl(t, td)
145  IMPLICIT NONE
146  real, intent(in) :: t, td
147 
148  get_tlcl = 1. / (1./(td - 56.) + log(t/td)/800.0) + 56.0
149 
150  return
151  end function get_tlcl
152 
153 !-----------------------------------------------------------------------+
154 ! mixing ratio in g/kg = water vapr/dry air
155 ! td in K
156 ! pres in Pa
157  elemental real function mixing_ratio(td, pres)
158  IMPLICIT NONE
159  real, intent(in) :: td, pres
160 
161  real corr, e
162 
163  corr = 1.001 + ( (pres/100.0 - 100) /900.) * 0.0034
164  e = getvaporpres(td) * corr
165  mixing_ratio = 0.62197 * (e/(pres-e)) * 1000.0
166 
167  return
168  end function mixing_ratio
169 
170 !-----------------------------------------------------------------------+
171 ! wet bulb temperature in K
172 ! t and td in K
173 ! pres in Pa
174  elemental real function getthetaw(t, td, pres)
175  IMPLICIT NONE
176  real, intent(in) :: t, td, pres
177 
178  real :: vl, cp, dt, top, bottom, twet, twetnew
179  real :: rmixr, smixr
180  integer :: i
181 
182  vl = 2.5e6
183  cp = 1004.0
184 
185  dt = 9.9e9
186  top = t
187  bottom = td
188 
189  do i = 1, 100
190 
191  twet = (top + bottom) / 2.0
192  rmixr = mixing_ratio(td, pres) / 1000.0
193  smixr = mixing_ratio(twet, pres) / 1000.0
194 
195  twetnew = t - (vl/cp) * (smixr-rmixr)
196  if(twetnew < twet) then
197  top = twet
198  else
199  bottom = twet
200  end if
201 
202  dt = abs(twet - twetnew)
203  if (dt <= 0.1) exit
204 
205  end do
206 
207  ! assign a value anyway, don't need (dt < 1.)
208  getthetaw = twet
209 
210  return
211  end function getthetaw
212 
213 !-----------------------------------------------------------------------+
214 ! Precipitation Condensate in g/kg
215  real function getprecipcond(totalCond, nz, topoK)
216  IMPLICIT NONE
217  real, intent(in) :: totalcond(nz)
218  integer, intent(in) :: nz, topok
219 
220  integer :: k
221 
222  getprecipcond = 0.0
223 
224  do k = topok, topok-2, -1 ! three levels
225  getprecipcond = totalcond(k) + getprecipcond
226  end do
227 
228  return
229  end function getprecipcond
230 
231 !-----------------------------------------------------------------------+
232 ! These 2-D indice are used for convective icing severity
233  subroutine calc_indice(t, td, pres, wvm, nz, topoK, &
234  kindex, liftedindex, totaltotals)
235  IMPLICIT NONE
236  integer, intent(in) :: nz, topok
237  real, intent(in) :: t(nz), td(nz), pres(nz), wvm(nz)
238  real, intent(out) :: kindex, liftedindex, totaltotals
239 
240  integer :: k
241  real :: t500hpa, t700hpa, t850hpa, dpt700hpa, dpt850hpa
242 
243  real :: surfacetemp,surfacedewpttemp,surfacepressure,surfacemixr
244  real :: tempatlcl, theta, pressatlcl, thetaeatlcl, tempfromthetae, tem
245 
246 ! write(0,*)' nz=',nz,' pres=',pres(:)
247  t500hpa = t(nz)
248  t700hpa = t(nz)
249  dpt700hpa = td(nz)
250  t850hpa = t(nz)
251  dpt850hpa = td(nz)
252 !
253  do k = nz, 2, -1
254 
255  ! use linear interpolation
256 
257 ! write(0,*)'k=',k,' pres=',pres(k)
258  if ((pres(k)- 50000.0 >= 0.) .and. (pres(k-1)- 50000.0 < 0.) ) then
259  if (abs(pres(k)- 50000.0) <= 0.1) then
260  t500hpa = t(k)
261  elseif (abs(pres(k-1)- 50000.0) <= 0.1) then
262  t500hpa = t(k-1)
263  else
264  t500hpa = t(k) - ((pres(k)-50000.0)/(pres(k)-pres(k-1))) * (t(k)-t(k-1))
265  end if
266  exit ! from surface to space, time to break the loop
267  end if
268 
269  if ((pres(k)- 70000.0 >= 0.) .and. (pres(k-1)- 70000.0 < 0.) ) then
270  if (abs(pres(k)- 70000.0) <= 0.1) then
271  t700hpa = t(k)
272  dpt700hpa = td(k)
273  elseif (abs(pres(k-1)- 70000.0) <= 0.1) then
274  t700hpa = t(k-1)
275  dpt700hpa = td(k-1)
276  else
277  tem = (pres(k)-70000.0)/(pres(k)-pres(k-1))
278  t700hpa = t(k) - tem * (t(k)-t(k-1))
279  dpt700hpa = td(k) - tem * (td(k)-td(k-1))
280  end if
281  endif
282 
283 ! write(0,*)'k=',k,' pres=',pres(k),pres(k-1)
284  if ((pres(k)- 85000.0 >= 0.) .and. (pres(k-1)- 85000.0 < 0.) ) then
285  if (abs(pres(k)- 85000.0) <= 0.1) then
286  t850hpa = t(k)
287  dpt850hpa = td(k)
288  elseif (abs(pres(k-1)- 85000.0) <= 0.1) then
289  t850hpa = t(k-1)
290  dpt850hpa = td(k-1)
291  else
292  tem = (pres(k)-85000.0)/(pres(k)-pres(k-1))
293  t850hpa = t(k) - tem * (t(k)-t(k-1))
294  dpt850hpa = td(k) - tem * (td(k)-td(k-1))
295  end if
296  endif
297 
298  end do
299 
300  ! kindex in C
301  kindex = (t850hpa-t500hpa) + (dpt850hpa-273.15) - (t700hpa-dpt700hpa)
302 
303  ! total total index
304  totaltotals = t850hpa + dpt850hpa - t500hpa * 2.0
305 
306  ! lifted index
307  ! use topoK instead of 1 - Y Mao
308  surfacetemp = t(topok)
309  surfacedewpttemp = td(topok)
310  surfacepressure = pres(topok)
311  surfacemixr = wvm(topok)
312 
313  tempatlcl = get_tlcl(surfacetemp, surfacedewpttemp)
314  theta = surfacetemp * (100000.0/surfacepressure)**0.286
315  pressatlcl = 100000.0 * (tempatlcl / theta)**3.4965
316  thetaeatlcl = getthetae(pressatlcl,tempatlcl,surfacemixr,tempatlcl)
317 
318  tempfromthetae = gettemperature(thetaeatlcl, 50000.0)
319 
320  liftedindex = t500hpa - tempfromthetae
321 
322  return
323  end subroutine calc_indice
324 
325 !-----------------------------------------------------------------------+
326  real function gettemperature(thetaE, pres)
327  IMPLICIT NONE
328  real, intent(in) :: thetae, pres
329 
330  real :: guess, w1, w2, thetae1, thetae2, cor
331  integer :: i
332 
333  guess = (thetae - 0.5 * (max(thetae - 270.0, 0.0))**1.05) * &
334  (pres/100000.0)**0.2
335 
336  do i = 1, 100
337  w1 = calcrsubs(pres, guess)
338  w2 = calcrsubs(pres, guess + 1.0)
339  thetae1 = getthetae(pres, guess, w1, guess)
340  thetae2 = getthetae(pres, guess + 1.0, w2, guess + 1.0)
341  cor = (thetae - thetae1)/(thetae2 - thetae1)
342  guess = guess + cor
343 
344  if (abs(cor) < 0.01) then
345  exit
346  end if
347  end do
348 
349  gettemperature = guess
350 
351  return
352  end function gettemperature
353 
354 !-----------------------------------------------------------------------+
355  elemental real function calcrsubs(pres, temp)
356  IMPLICIT NONE
357  real, intent(in) :: pres, temp
358 
359  real :: esubs
360 
361  esubs = calcesubs(temp)
362  calcrsubs = (0.622*esubs)/(pres - esubs)
363 
364  return
365  end function calcrsubs
366 
367 !-----------------------------------------------------------------------+
368  elemental real function calcesubs(temp)
369  IMPLICIT NONE
370  real, intent(in) :: temp
371 
372  real :: x
373  real, parameter :: coeffs(9) = &
374  (/ 610.5851, 44.40316, 1.430341, 0.02641412, 2.995057e-4,&
375  2.031998e-6, 6.936113e-9, 2.564861e-12, -3.704404e-14 /)
376 
377 
378  x= max(-80.0, temp - 273.15)
379 
380  calcesubs = coeffs(1) + x*(coeffs(2) + x*(coeffs(3) + x*(coeffs(4) + &
381  x*(coeffs(5) + x*(coeffs(6) + x*(coeffs(7) + &
382  x*(coeffs(8) + x*coeffs(9))))))))
383  return
384  end function calcesubs
385 
386 !-----------------------------------------------------------------------+
387  subroutine calc_totalwaterpath(hgt, pres, t, totalWater, nz, twp)
388  IMPLICIT NONE
389  integer, intent(in) :: nz
390  real, intent(in) :: hgt(nz), pres(nz), t(nz), totalwater(nz)
391  real, intent(out):: twp(nz)
392 
393  real :: condensateintegrated
394  integer :: topidx, baseidx
395  integer :: k, m
396 
397  lp100: do k = 1, nz
398 
399  twp(k) = 0.0
400  topidx = -1
401  baseidx = -1
402 
403  ! get the layer top and base for the total condensate layer
404  lp200: do m = k, 2, -1
405  if (totalwater(m) > 0.001) then
406  topidx = m
407  else
408  exit
409  end if
410  end do lp200
411 
412  lp300: do m = k, nz
413  if (totalwater(m) > 0.001) then
414  baseidx = m
415  else
416  exit
417  end if
418  end do lp300
419 
420  ! get the integrated condensate from the rep to the layer top
421  if(topidx == -1 .or. baseidx == -1) cycle
422  condensateintegrated = 0.0
423  lp400: do m = topidx, baseidx
424  if (t(m) <= 273.15) then
425  condensateintegrated = condensateintegrated + &
426  ( ((totalwater(m)*pres(m)) / (287.05 * t(m)) ) * &
427  (hgt(m-1) - hgt(m) + 0.001))
428  end if
429  end do lp400
430 
431  twp(k) = condensateintegrated
432 
433  end do lp100
434 
435  return
436  end subroutine calc_totalwaterpath
437 
438 !-----------------------------------------------------------------------+
439  integer function getpreciptype(imp_physics,pres, hgt, t, rh, wbt, nz, &
440  hprcp, hcprcp, pc, cin, cape, lx)
441  IMPLICIT NONE
442  integer, intent(in) :: imp_physics
443  integer, intent(in) :: nz
444  real, intent(in) :: pres(nz), hgt(nz), t(nz), rh(nz), wbt(nz)
445  real, intent(in) :: hprcp,hcprcp, pc, cin, cape, lx
446 
447  integer, allocatable :: wxtype(:)
448 
449  integer :: idxmelt, idxrefz, iceidx
450  real :: coldtemp, tcoldarea, wetbuldarea
451 
452  real :: precip, precip_standard
453 
454  integer :: k, m
455 
456  allocate(wxtype(nz))
457 
458  idxmelt = nz ! melting
459  idxrefz = nz ! refreezing
460 
461  if(imp_physics == 98 .or. imp_physics == 99) then
462  precip = hprcp
463  precip_standard =0.045
464  else if(imp_physics == 11 .or. imp_physics == 8) then
465  precip = pc
466  precip_standard = 0.01
467  end if
468 
469  ! Check for convection first
470  if_conv: if ( hcprcp >= 1.0 .and. cape >= 1000.0 .and. cin > -100. &
471  .and. lx < 0.0) then
472  getpreciptype = precips% CONVECTION
473  else
474 
475  ! find the height(index) where the snow melts, wetBulbTemp > 0C
476  lp100: do k = 1, nz
477  if( wbt(k) > 273.15) then
478  idxmelt = k
479  ! now look for a refreezing layer below warmnose
480  do m = idxmelt, nz
481  if (t(m) < 273.15) then
482  idxrefz = m
483  exit
484  end if
485  end do
486  exit
487  end if
488  end do lp100
489 
490  ! find the level that is below 150hPa or
491  ! the min and max wet bulb temperature
492  lp200: do k = nz, 1, -1
493 
494  wxtype(k) = 0
495 
496  if_pres_wbt: if ( pres(k) >= 15000. .or. &
497  (wbt(k) >= 200. .and. wbt(k) <= 1000.)) then
498  iceidx = k
499  coldtemp = t(k)
500 
501  do m = k, 1, -1
502  ! look for level below cloud_top_pressure and
503  ! cloud_top_altitude
504  ! look for ctt
505  if( (pres(m)>=30000. .or. (hgt(m)>hgt(nz)+700.))&
506  .and. (rh(m) > 90.) .and. (t(m) < coldtemp)) then
507  coldtemp = t(m)
508  iceidx = m
509  end if
510  end do
511 
512  ! found a cloud -- identify the precipitaiton type
513  if_iceidx: if (iceidx /= k) then
514 
515  ! Sum the pos. area between wetBulbTemp and
516  ! the 0 deg isotherm below coldTemp
517  tcoldarea = 0.0
518  do m = k, iceidx, -1
519  if (wbt(m) >= 273.15) then
520  tcoldarea = tcoldarea + (wbt(m)-273.15) * &
521  (hgt(m-1) - hgt(m))
522  end if
523  end do
524 
525  ! sum the total area between the wet bulb T and the 0 C isotherm
526  ! in the lowest precip_type_altitude
527  wetbuldarea = 0.0
528  do m = iceidx, nz
529  if ( (hgt(m) <= hgt(nz)+1500) .and. &
530  (m > idxrefz)) then
531  wetbuldarea = wetbuldarea + (wbt(m)-273.15) * &
532  (hgt(m-1) - hgt(m))
533  end if
534  end do
535 
536  ! get the probable Precip type
537  if (coldtemp > 265.15) then
538  wxtype(k) = precips% OTHER
539  else if (tcoldarea < 350.) then
540  if(t(k) <= 273.15) then
541  wxtype(k) = precips% SNOW
542  else
543  wxtype(k) = precips% RAIN
544  end if
545  else if (wetbuldarea <= -250.) then
546  wxtype(k) = precips% OTHER
547  else if(t(k) <= 273.15) then
548  wxtype(k) = precips% OTHER
549  else
550  wxtype(k) = precips% RAIN
551  end if
552 
553  else
554  wxtype(k) = precips% NONE
555  end if if_iceidx
556  end if if_pres_wbt
557 
558  end do lp200
559 
560  getpreciptype = precips% NONE
561 ! if(hprcp >= 0.045) then
562 ! if(pc >= 0.01) then ! g/kg
563  if(precip >= precip_standard) then
564  do k = nz, nz-1, -1
565  if (wxtype(k) > getpreciptype) then
566  getpreciptype = wxtype(k)
567  end if
568  end do
569  end if
570  end if if_conv
571 
572  deallocate(wxtype)
573 
574  return
575  end function getpreciptype
576 
577 end module derivedfields
578 
579 !========================================================================
580 ! = = = = = = = = = = = = = module CloudLayers = = = = = = = = = = = = =
581 !========================================================================
582 module cloudlayers
583 
584  use derivedfields, only : mixing_ratio
585 
586  implicit none
587 
588  private
589  public calc_cloudlayers
590  public clouds_t
591 
592  integer, parameter :: maxlayers = 30
593  type :: clouds_t
594  ! 2-D
595  integer :: nlayers
596  integer :: wmnidx ! warm nose index
597  real :: avv ! average vertical velocity
598  ! 3-D, on model levels of nz
599  real, allocatable :: layerq(:)
600  ! 3-D, of cloud layers
601  integer :: topidx(maxlayers)
602  integer :: baseidx(maxlayers)
603  real :: ctt(maxlayers)
604  end type clouds_t
605 
606 contains
607 
608 !-----------------------------------------------------------------------+
609  subroutine calc_cloudlayers(rh,t,pres,ept,vv, nz, topoK, xlat, xlon,&
610  region, clouds)
611  IMPLICIT NONE
612  integer, intent(in) :: nz, topok
613  real, intent(in) :: rh(nz),t(nz),pres(nz), ept(nz), vv(nz)
614  real, intent(in) :: xlat,xlon
615  integer, intent(out) :: region
616  type(clouds_t), intent(inout) :: clouds ! layerQ has been allocated
617 
618  real :: t_rh
619  integer :: in_cld,cur_base, num_lyr
620 
621  integer :: k, kk, n, kstart
622 
623  ! get the global region and set the RH thresholds
624  if (abs(xlat)<23.5) then
625  t_rh = 80.0
626  region = 1
627  elseif ( (abs(xlat)>=23.5).and.(abs(xlat)<66)) then
628  t_rh = 75.0
629  region =2
630  else
631  t_rh = 70.0
632  region =2
633  endif
634 
635  ! loop from the top (start at 2 so we can have n+1) )
636  ! bottom is nz and top is 1
637 
638  num_lyr = 0
639  in_cld = 0
640  cur_base = 1
641 
642  ! identify the very top layer if rh starts at high value
643  if((rh(1) >= t_rh) .and. (rh(2) >= t_rh))then
644  num_lyr = num_lyr + 1
645  in_cld = 1
646  ! find the cloud top and ctt
647  clouds%topIdx(num_lyr) = 1
648  clouds%ctt(num_lyr) = t(1)
649 
650  ! find the cloud base
651  do kk=2,topok-1
652  if ((rh(kk-1)>=t_rh).and.(rh(kk)<t_rh)) then
653  if ((rh(kk+1)<t_rh).and.(in_cld==1)) then
654  clouds%baseIdx(num_lyr) = kk
655  cur_base = kk
656  in_cld = 0
657  endif
658  endif
659  end do
660  end if
661  kstart = cur_base + 1
662 
663  do k = kstart, topok-1
664  if (cur_base<=k) then
665  if ((rh(k-1)<t_rh).and.(in_cld==0)) then
666  if ((rh(k)>=t_rh).and.(rh(k+1)>=t_rh)) then
667  num_lyr = num_lyr + 1
668  in_cld = 1
669  ! find the cloud top and ctt
670  clouds%topIdx(num_lyr) = k
671  clouds%ctt(num_lyr) = t(k)
672 
673  ! find the cloud base
674  do kk=clouds%topIdx(num_lyr),topok-1
675  if ((rh(kk-1)>=t_rh).and.(rh(kk)<t_rh)) then
676  if ((rh(kk+1)<t_rh).and.(in_cld==1)) then
677  clouds%baseIdx(num_lyr) = kk
678  cur_base = kk
679  in_cld = 0
680  endif
681  endif
682  end do
683 
684  endif
685  endif
686  endif
687  end do
688 
689  ! if the loop exits still in cloud make the cloud base the surface
690  if (in_cld==1) then
691  clouds%baseIdx(num_lyr) = topok
692  endif
693 
694  clouds%nLayers = num_lyr
695 
696  clouds%wmnIdx = getwarmnoseidx(t, nz)
697 
698  ! only for the lowest cloud layer
699  ! only used for severity scenario where precip is below warm nose
700  if(num_lyr > 0) then
701  clouds%avv = getaveragevertvel(t, vv, nz, &
702  clouds%topIdx(num_lyr), clouds%baseIdx(num_lyr))
703  else
704  clouds%avv = 0
705  end if
706 
707  ! for severity
708  call calc_layerq(t, rh, pres, ept, nz, clouds)
709 
710  return
711  end subroutine calc_cloudlayers
712 
713 
714 !-----------------------------------------------------------------------+
715 ! From top to the surface, find the warmnoseIndex if existing
716 ! |
717 ! |T<=FREEZING
718 ! ----|------------------warmnoseIndex
719 ! | |
720 ! | T > FREEZING|
721 ! --------------
722 ! |T<=FREEZING
723 ! ____|________
724 ! ///////////// surface
725 !-----------------------------------------------------------------------+
726  integer function getwarmnoseidx(t, nz)
727  IMPLICIT NONE
728  integer, intent(in) :: nz
729  real, intent(in) :: t(nz)
730 
731  logical :: abovefreezing
732  integer :: tmpindex
733 
734  integer :: k
735 
736  abovefreezing = .false.
737  tmpindex = 0 !It doesn't matter of the initialized value
738 
739  getwarmnoseidx = -1
740 
741  ! search from top of model down to the surface
742  do k = 1, nz
743 
744  if(t(k) <= 273.15) then
745  if (abovefreezing) then ! above 0 then below 0
746  getwarmnoseidx = tmpindex
747  ! once warmnoseIndex is set, we can't get back
748  ! here. let's break out of loop
749  exit
750  end if !aboveFreezing
751  abovefreezing = .false.
752  else
753  if(.not. abovefreezing) tmpindex = k
754  abovefreezing = .true.
755  end if
756 
757  end do
758 
759  return
760  end function getwarmnoseidx
761 
762 !-----------------------------------------------------------------------+
763 ! Calculate the average VV in the dendritic layer (-13 to -16) if it's in
764 ! the lowest cloud (clouds%nLayers). If the -13 to -16 layer falls
765 ! between 2 levels then use the average VV from the levels on either side
766 !
767  real function getaveragevertvel(t,vv,nz, topIdx_lowest,baseIdx_lowest)
768  IMPLICIT NONE
769  integer, intent(in) :: nz
770  real, intent(in) :: t(nz), vv(nz)
771  integer, intent(in) :: topidx_lowest,baseidx_lowest
772 
773  integer :: numvertvel, k, start_base
774 
775  real :: sumvertvel
776 
777  sumvertvel = 0.0
778  numvertvel = 0
779 
780  ! The following loop starts at least nz-1 as the lowest model level
781  if(baseidx_lowest == nz) then
782  start_base = nz - 1
783  else
784  start_base = baseidx_lowest
785  end if
786 
787  do k = start_base, topidx_lowest, -1
788  if ( t(k) <= 260.15 .and. t(k) >= 257.15) then
789  sumvertvel = sumvertvel + vv(k)
790  numvertvel = numvertvel + 1
791  else if (t(k) < 257.15) then
792  sumvertvel = sumvertvel + vv(k) + vv(k+1)
793  numvertvel = 2
794  exit
795  end if
796  end do
797 
798  if (numvertvel == 0) then
799  getaveragevertvel = 0.0
800  else
801  getaveragevertvel = sumvertvel / numvertvel
802  end if
803 
804  return
805  end function getaveragevertvel
806 
807 !-----------------------------------------------------------------------+
808  subroutine calc_layerq(t, rh, pres, ept, nz, clouds)
809  IMPLICIT NONE
810  integer, intent(in) :: nz
811  real, intent(in) :: t(nz), rh(nz), pres(nz), ept(nz)
812  type(clouds_t), intent(inout) :: clouds
813 
814 
815  real :: mean_rh, te_map, rh_map, adj_q
816  real :: totalq, testthetae, q_top, q_base, q_layer, delta_te, sum_rh
817  integer :: base_k, test_k, num_layers
818 
819  integer :: k, m, n, kk
820 
821  clouds%layerQ(:) = 0.0
822 
823  ! loop through each layer
824  lp_nlayers: do n = 1, clouds%nLayers
825  lp_k_inlayer: do k = clouds%topIdx(n), clouds%baseIdx(n)-1
826 
827  totalq = 0.0
828  testthetae = ept(k)
829 
830  ! get base layer k value (first more than 4K colder than at k)
831  base_k = clouds%baseIdx(n)
832  do m = k, nz-1
833  if((testthetae-ept(m)>4.) .and. (clouds%baseIdx(n)<=m))then
834  base_k = m - 1
835  exit
836  end if
837  end do
838 
839  ! calculate delta thetaE and deltaQ and deltaRH for the layer
840  lp_kk: do kk = k, base_k-1
841  q_top = mixing_ratio(t(kk), pres(kk))
842  q_base = mixing_ratio(t(kk+1), pres(kk+1))
843  q_layer = q_base - q_top
844 
845  ! convert Q from g/kg to g/m**3
846  q_layer = q_layer * (pres(kk)/(287.05 * t(kk)))
847 
848  if (q_layer < 0.0) q_layer = 0.0
849 
850  delta_te = testthetae - ept(kk+1)
851  test_k = kk + 1
852 
853  ! get the mean RH in the layer
854  num_layers = 0
855  sum_rh = 0.0
856  do m = k, test_k
857  num_layers = num_layers + 1
858  sum_rh = sum_rh + rh(m)
859  end do
860 
861  if (num_layers == 0) then
862  mean_rh = 0.0
863  else
864  mean_rh = sum_rh / num_layers
865  end if
866 
867  te_map = dq_delta_te_map(delta_te)
868  rh_map = dq_rh_map(mean_rh)
869  adj_q = q_layer * te_map * rh_map
870 
871  totalq = totalq + adj_q
872  end do lp_kk
873 
874  clouds%layerQ(k) = totalq
875  end do lp_k_inlayer
876  end do lp_nlayers
877 
878  return
879  end subroutine calc_layerq
880 
881 !-----------------------------------------------------------------------+
882  ! 70.0 0.0, 100.0 1.0
883  real function dq_rh_map(rh)
884  IMPLICIT NONE
885  real, intent(in) :: rh
886 
887  if(rh <= 70.) then
888  dq_rh_map = 0.
889  else if( rh >= 100.) then
890  dq_rh_map = 1.0
891  else
892  dq_rh_map = (rh-70.)/30.
893  end if
894 
895  return
896  end function dq_rh_map
897 
898 !-----------------------------------------------------------------------+
899  ! 0.0 1.0, 4.0 0.0
900  real function dq_delta_te_map(delTE)
901  IMPLICIT NONE
902  real, intent(in) :: delte
903 
904  if (delte <= 0.0) then
905  dq_delta_te_map = 1.0
906  else if (delte <= 4.0) then
907  dq_delta_te_map = (4. - delte) / 4.0
908  else
909  dq_delta_te_map = 0.0
910  end if
911 
912  return
913  end function dq_delta_te_map
914 
915 end module cloudlayers
916 
917 
918 !========================================================================
919 ! = = = = = = = = = = = = module IcingPotential = = = = = = = = = = = = =
920 !========================================================================
921 module icingpotential
922 
923  use ctlblk_mod, only: me
924 
925  use derivedfields, only : precips
926  use cloudlayers, only : clouds_t
927 
928  implicit none
929 
930  private
931  public icing_pot
932 
933 contains
934 
935 !-----------------------------------------------------------------------+
936  subroutine icing_pot(hgt, rh, t, liqCond, vv, nz, clouds, ice_pot)
937  IMPLICIT NONE
938  integer, intent(in) :: nz
939  real, intent(in) :: hgt(nz), rh(nz), t(nz), liqcond(nz), vv(nz)
940  type(clouds_t), intent(in) :: clouds
941  real, intent(out) :: ice_pot(nz)
942 
943  real :: ctt, slw, slw_fac, vv_fac
944  integer :: k, n
945 
946  ! variables for printout only
947  real :: mapt,maprh,mapctt,mapvv,mapslw
948 
949  ice_pot(:) = 0.0
950 
951  ! run the icing potential algorithm
952  lp_n: do n = 1, clouds%nLayers
953 
954  ! apply the cloud top temperature to the cloud layer
955  ctt = clouds%ctt(n)
956 
957  lp_k: do k = clouds%topIdx(n), clouds%baseIdx(n)
958 
959  ! convert the liquid water to slw if the CTT > -14C
960  if ( ctt>=259.15 .and. ctt<=273.15 ) then
961  slw = liqcond(k)
962  else
963  slw = 0.0
964  end if
965 
966  ice_pot(k) = tmap(t(k)) * rh_map(rh(k)) * ctt_map(ctt)
967 
968  ! add the VV map
969  if (vv_map(vv(k))>=0.0) then
970  vv_fac = (1.0-ice_pot(k)) * (0.25*vv_map(vv(k)))
971  else
972  vv_fac = ice_pot(k) * (0.25*vv_map(vv(k)))
973  endif
974 
975  ! add the slw
976  slw_fac = (1.0 - ice_pot(k))*(0.4*slw_map(slw))
977 
978  ! calculate the final icing potential
979  if (ice_pot(k)>0.001) then
980  ice_pot(k) = ice_pot(k) + vv_fac + slw_fac
981  endif
982 
983  ! make sure the values don't exceed 1.0
984  if (ice_pot(k)>1.0) then
985  ice_pot(k) = 1.0
986  endif
987 
988 ! mapt=tmap(t(k))
989 ! maprh=rh_map(rh(k))
990 ! mapctt=ctt_map(ctt)
991 ! mapvv=vv_map(vv(k))
992 ! mapslw=slw_map(slw)
993 ! write(*,'(2x,I3,A,1x,A,I3,F7.3)')me,":","k,icip=",k,ice_pot(k)
994 ! write(*,'(2x,I3,A,1x,F8.2,F7.3,F7.3,F7.3)') &
995 ! me,":",t(k),mapt,rh(k),maprh
996 ! write(*,'(2x,I3,A,1x,F8.2,F7.3,F10.6,F7.3)') &
997 ! me,":",ctt,mapctt,vv(k),mapvv
998 ! write(*,'(2x,I3,A,1x,F11.6,F11.6,F7.3)') &
999 ! me,":",liqCond(k), slw, mapslw
1000 
1001  end do lp_k
1002  end do lp_n
1003 
1004  return
1005  end subroutine icing_pot
1006 
1007 !-----------------------------------------------------------------------+
1008  real function tmap(temp)
1009  IMPLICIT NONE
1010  real, intent(in) :: temp
1011 
1012  if((temp>=248.15).and.(temp<=263.15)) then
1013  tmap=((temp-248.15)/15.0)**1.5
1014  elseif((temp>263.15).and.(temp<268.15)) then
1015  tmap=1.0
1016  elseif((temp>268.15).and.(temp<273.15)) then
1017  tmap=((273.15 - temp)/5.0)**0.75
1018  else
1019  tmap=0.0
1020  endif
1021 
1022  return
1023  end function tmap
1024 
1025 
1026 !-----------------------------------------------------------------------+
1027  real function rh_map(rh)
1028  IMPLICIT NONE
1029  real, intent(in) :: rh
1030 
1031  if (rh>95.0) then
1032  rh_map=1.0
1033  elseif ( (rh<=95.0).and.(rh>=50.0) ) then
1034  rh_map=((20./9.) * ((rh/100.0)-0.5))**3.0
1035  else
1036  rh_map=0.0
1037  endif
1038 
1039  return
1040  end function rh_map
1041 
1042 
1043 !-----------------------------------------------------------------------+
1044  real function ctt_map(ctt)
1045  IMPLICIT NONE
1046  real, intent(in) :: ctt
1047 
1048  if((ctt>=261.0).and.(ctt<280.)) then
1049  ctt_map = 1.0
1050  elseif((ctt>223.0).and.(ctt<261.0 )) then
1051  ctt_map = 0.2 + 0.8*(((ctt - 223.0)/38.0)**2.0)
1052  elseif(ctt<223.0) then
1053  ctt_map = 0.2
1054  else
1055  ctt_map = 0.0
1056  endif
1057 
1058  return
1059  end function ctt_map
1060 
1061 
1062 !-----------------------------------------------------------------------+
1063 
1064  real function vv_map(vv)
1065  IMPLICIT NONE
1066  real, intent(in) :: vv
1067 
1068  if (vv>0.0) then
1069  vv_map = 0.0
1070  elseif (vv<-0.5) then
1071  vv_map = 1.0
1072  else
1073  vv_map = -1.0 * (vv/0.5)
1074  endif
1075 
1076  return
1077  end function vv_map
1078 
1079 
1080 !-----------------------------------------------------------------------+
1081 
1082  real function slw_map(slw)
1083  IMPLICIT NONE
1084  real, intent(in) :: slw
1085 
1086  if(slw>0.2) then
1087  slw_map = 1.0
1088  elseif (slw<=0.001) then
1089  slw_map = 0.0
1090  else
1091  slw_map = slw/0.2
1092  endif
1093 
1094  return
1095  end function slw_map
1096 
1097 end module icingpotential
1098 
1099 
1100 !========================================================================
1101 ! = = = = = = = = = = = = module SeverityMaps = = = = = = = = = = = = =
1102 !========================================================================
1103 module severitymaps
1104  implicit none
1105  public
1106  public scenarios
1107 
1108  type :: scenarios_t
1109  integer :: no_precipitaion = 0
1110  integer :: precipitaion_below_warmnose = 1
1111  integer :: precipitaion_above_warmnose = 2
1112  integer :: all_snow = 3
1113  integer :: cold_rain = 4
1114  integer :: warm_precipitaion = 5
1115  integer :: freezing_precipitaion = 6
1116  integer :: convection = 7
1117  end type scenarios_t
1118  type(scenarios_t), parameter :: scenarios = scenarios_t()
1119 
1120 contains
1121 
1122 !-----------------------------------------------------------------------+
1123 ! scenario dependant
1124 !-----------------------------------------------------------------------+
1125 
1126  real function twp_map(v, scenario)
1127  implicit none
1128  real, intent(in) :: v
1129  integer, intent(in) :: scenario
1130  twp_map = 0.0
1131  select case (scenario)
1132  case(scenarios%PRECIPITAION_ABOVE_WARMNOSE, scenarios%ALL_SNOW, &
1133  scenarios%COLD_RAIN, scenarios%FREEZING_PRECIPITAION)
1134  ! 0.0 0.0, 1000.0 1.0
1135  if(v <= 0.0) then
1136  twp_map = 0.0
1137  else if( v < 1000.) then
1138  twp_map = v / 1000.0
1139  else
1140  twp_map = 1.
1141  end if
1142  case( scenarios%WARM_PRECIPITAION)
1143  ! 0.0 0.0, 500.0 1.0
1144  if(v <= 0.0) then
1145  twp_map = 0.0
1146  else if( v < 500.) then
1147  twp_map = v / 500.0
1148  else
1149  twp_map = 1.
1150  end if
1151  end select
1152  return
1153  end function twp_map
1154 
1155  ! Only precip below warmnose has a different temperature map
1156  real function t_map(v, scenario)
1157  implicit none
1158  real, intent(in) :: v
1159  integer, intent(in) :: scenario
1160  t_map = 0.
1161  select case (scenario)
1162  case(scenarios%PRECIPITAION_BELOW_WARMNOSE)
1163  ! 269.15 1.0, 273.15 0.0
1164  if(v <= 269.15) then
1165  t_map = 1.
1166  elseif( v <= 273.15) then
1167  t_map = (273.15 - v ) / 4.0
1168  else
1169  t_map = 0.
1170  end if
1171  case default
1172  ! 248.15 1.0, 248.65 0.8039, 249.15 0.7226, 251.15 0.5196,
1173  ! 253.15 0.3798, 255.15 0.2662, 257.15 0.1679, 259.15 0.0801,
1174  ! 261.15 0.0, 268.15 0.0, 273.15 1.0
1175  if(v <= 248.15) then
1176  t_map = 1.0
1177  elseif( v <= 248.65) then
1178  t_map = (49.162215 - 0.1961*v) * 2.
1179  elseif( v <= 249.15) then
1180  t_map = (20.617195 - 0.0813*v) * 2.
1181  elseif(v <= 251.15) then
1182  t_map = (52.02265 - 0.203*v) / 2.0
1183  elseif(v <= 253.15) then
1184  t_map = (36.14997 - 0.1398*v) / 2.0
1185  elseif(v <= 255.15) then
1186  t_map = (29.51744 - 0.1136*v) / 2.0
1187  elseif(v <= 257.15) then
1188  t_map = (25.613645-0.0983*v) / 2.0
1189  elseif(v <= 259.15) then
1190  t_map = (22.91357 - 0.0878*v) / 2.0
1191  elseif( v <= 261.15) then
1192  t_map = (20.918115 - 0.0801*v) / 2.
1193  else
1194  t_map = 0.0
1195  end if
1196  end select
1197  end function t_map
1198 
1199 
1200  ! Condensates near the surface take place of radar reflectivity in CIP
1201  real function prcpcondensate_map(v, scenario)
1202  implicit none
1203  real, intent(in) :: v
1204  integer, intent(in) :: scenario
1205  prcpcondensate_map = 0.0
1206  select case (scenario)
1207  case(scenarios%NO_PRECIPITAION)
1208  ! do nothing
1209  case(scenarios%PRECIPITAION_BELOW_WARMNOSE, scenarios%COLD_RAIN, &
1210  scenarios%PRECIPITAION_ABOVE_WARMNOSE)
1211  ! 0.05 0.0, 0.2 1.0
1212  if(v <= 0.05) then
1213  prcpcondensate_map = 0.
1214  elseif(v <= 0.2) then
1215  prcpcondensate_map = (v - 0.05) / 0.15
1216  else
1217  prcpcondensate_map = 1.
1218  end if
1219  case(scenarios%ALL_SNOW)
1220  ! 0.05 0.0, 0.25 1.0
1221  if(v <= 0.05) then
1222  prcpcondensate_map = 0.
1223  elseif(v <= 0.25) then
1224  prcpcondensate_map = (v - 0.05) / 0.2
1225  else
1226  prcpcondensate_map = 1.0
1227  end if
1228  case(scenarios%WARM_PRECIPITAION, scenarios%FREEZING_PRECIPITAION)
1229  ! 0.05 0.0, 0.15 0.5
1230  if(v <= 0.05) then
1231  prcpcondensate_map = 0.
1232  elseif(v <= 0.15) then
1233  prcpcondensate_map = (.5*v - 0.025) / 0.1
1234  else
1235  prcpcondensate_map = 0.5
1236  end if
1237  end select
1238  return
1239  end function prcpcondensate_map
1240 
1241 
1242  real function deltaz_map(v, scenario)
1243  implicit none
1244  real, intent(in) :: v
1245  integer, intent(in) :: scenario
1246  deltaz_map = 0.0
1247  select case (scenario)
1248  case (scenarios%NO_PRECIPITAION, scenarios%WARM_PRECIPITAION)
1249  ! 0.0 0.0, 1828.8 1.0
1250  if(v <= 0.0) then
1251  deltaz_map = 0.0
1252  else if(v <= 1828.8) then
1253  deltaz_map = v /1828.8
1254  else
1255  deltaz_map = 1.
1256  end if
1257  case(scenarios%PRECIPITAION_BELOW_WARMNOSE)
1258  ! 0.0 0.0, 30.49999 0.0, 30.5 1.0
1259  if(v < 30.5) then
1260  deltaz_map = 0.0
1261  else
1262  deltaz_map = 1.
1263  end if
1264  case(scenarios%ALL_SNOW)
1265  ! 914.4 0.0, 2743.2 1.0
1266  if(v <= 914.4) then
1267  deltaz_map = 0.0
1268  else if(v <= 2743.2) then
1269  deltaz_map = (v - 914.4) / 1828.8
1270  else
1271  deltaz_map = 1.
1272  end if
1273  case(scenarios%COLD_RAIN, scenarios%FREEZING_PRECIPITAION)
1274  ! 1524.0 0.0, 4267.2 1.0
1275  if(v <= 1524.) then
1276  deltaz_map = 0.0
1277  else if(v <= 4267.2) then
1278  deltaz_map = (v - 1524.) / 2743.2
1279  else
1280  deltaz_map = 1.
1281  end if
1282  end select
1283  return
1284  end function deltaz_map
1285 
1286 !-----------------------------------------------------------------------+
1287 ! scenario independant.
1288 !-----------------------------------------------------------------------+
1289 
1290  ! 223.15 0.8, 233.15 0.7446, 243.15 0.5784, 253.15 0.3014
1291  ! 261.15 0.0, 280.15 0.0, 280.151 1.0
1292  real function ctt_map(v)
1293  implicit none
1294  real, intent(in) :: v
1295  if(v <= 223.15) then
1296  ctt_map = 0.8
1297  elseif(v <= 233.15) then
1298  ctt_map = (20.36251 - 0.0554*v) / 10.
1299  elseif(v <= 243.15) then
1300  ctt_map = (46.19553 - 0.1662*v) / 10.
1301  elseif(v <= 253.15) then
1302  ctt_map = (73.13655 - 0.277*v) / 10.
1303  elseif(v <= 261.15) then
1304  ctt_map = (78.71061 - 0.3014*v) / 8.
1305  else
1306  ctt_map = 0.
1307  end if
1308  end function ctt_map
1309 
1310  ! -0.5 1.0, 0.0 0.0
1311  real function vv_map(v)
1312  implicit none
1313  real, intent(in) :: v
1314  if(v <= -0.5) then
1315  vv_map = 1.
1316  else if(v <= 0.) then
1317  vv_map = -2. * v
1318  else
1319  vv_map = 0.
1320  end if
1321  end function vv_map
1322 
1323 
1324  ! cloud top distance
1325  ! 609.6 1.0, 3048.0 0.0
1326  real function cldtopdist_map(v)
1327  implicit none
1328  real, intent(in) :: v
1329  if( v <= 609.6) then
1330  cldtopdist_map = 1.0
1331  elseif( v <= 3048.) then
1332  cldtopdist_map = (3048. - v) / 2438.4
1333  else
1334  cldtopdist_map = 0.
1335  end if
1336  end function cldtopdist_map
1337 
1338 
1339  ! cloud base distance
1340  ! 304.8 1.0, 1524.0 0.0
1341  real function cldbasedist_map(v)
1342  implicit none
1343  real, intent(in) :: v
1344  if( v <= 304.8) then
1345  cldbasedist_map = 1.0
1346  elseif( v <= 1524.) then
1347  cldbasedist_map = (1524. - v) / 1219.2
1348  else
1349  cldbasedist_map = 0.
1350  end if
1351  end function cldbasedist_map
1352 
1353  ! 0.0 0.0, 1.0 1.0
1354  real function deltaq_map(v)
1355  implicit none
1356  real, intent(in) :: v
1357  if( v <= 0.) then
1358  deltaq_map = 0
1359  elseif( v <= 1.0) then
1360  deltaq_map = v
1361  else
1362  deltaq_map = 1.
1363  end if
1364  end function deltaq_map
1365 
1366  real function moisture_map_cond(rh, liqCond, iceCond, pres, t)
1367  IMPLICIT NONE
1368  real, intent(in) :: rh, liqcond, icecond, pres, t
1369  real :: liq_m3, ice_m3, rhfactor, liqfactor,icefactor
1370  ! Convert liqCond, iceCond from g/kg to g/m^3
1371  liq_m3 = (liqcond * pres) / (287.05 * t)
1372  ice_m3 = (icecond * pres) / (287.05 * t)
1373  rhfactor = 0.5 * rh_map(rh)
1374  liqfactor = 0.6* condensate_map(liq_m3)
1375  icefactor = 0.4 * condensate_map(ice_m3)
1376  moisture_map_cond = rhfactor + (1.0 - rhfactor) * (liqfactor + icefactor)
1377  return
1378  end function moisture_map_cond
1379 
1380  ! If not identify liquid/ice condensate
1381  real function moisture_map_cwat(rh, cwat, pres, t)
1382  IMPLICIT NONE
1383  real, intent(in) :: rh, cwat, pres, t
1384  real :: cwat_m3, rhfactor, cwatfactor
1385  ! Convert cwat from g/kg to g/m^3
1386  cwat_m3 = (cwat * pres) / (287.05 * t)
1387  rhfactor = 0.5 * rh_map(rh)
1388  cwatfactor = condensate_map(cwat_m3)
1389  moisture_map_cwat = rhfactor + (1.0 - rhfactor) * cwatfactor
1390  return
1391  end function moisture_map_cwat
1392 
1393  ! only called by moisture_map
1394  ! 70.0 0.0, 100.0 1.0
1395  real function rh_map(v)
1396  implicit none
1397  real, intent(in) :: v
1398  if(v <= 70.) then
1399  rh_map = 0.
1400  else if(v <= 100) then
1401  rh_map = (v - 70.) / 30.
1402  else
1403  rh_map = 1.
1404  end if
1405  end function rh_map
1406 
1407  ! only called by moisture_map
1408  ! 0.00399 0.0, 0.004 0.0, 0.2 1.0
1409  real function condensate_map(v)
1410  implicit none
1411  real, intent(in) :: v
1412  if(v <= 0.004) then
1413  condensate_map = 0.
1414  elseif( v <= 0.2) then
1415  condensate_map = (v - 0.004) / 0.196
1416  else
1417  condensate_map = 1.
1418  end if
1419  end function condensate_map
1420 
1421 
1422 !-----------------------------------------------------------------------+
1423 ! convective
1424 !-----------------------------------------------------------------------+
1425 
1426  ! 243.150 0.0, 265.15 1.0, 269.15 1.0, 270.15 0.87
1427  ! 271.15 0.71, 272.15 0.50, 273.15 0.0
1428  real function convect_t_map(v)
1429  implicit none
1430  real, intent(in) :: v
1431  if(v <= 243.15) then
1432  convect_t_map = 0.
1433  else if(v <= 265.15) then
1434  convect_t_map = (v -243.15)/22.
1435  else if(v <= 269.15) then
1436  convect_t_map = 1.
1437  else if(v <= 270.15) then
1438  convect_t_map = -0.13*v + 35.9895
1439  else if(v <= 271.15) then
1440  convect_t_map = -0.16*v + 44.094
1441  else if(v <= 272.15) then
1442  convect_t_map = -0.21*v + 57.6515
1443  else if(v <= 273.15) then
1444  convect_t_map = -0.50*v + 136.575
1445  else
1446  convect_t_map = 0.
1447  end if
1448  end function convect_t_map
1449 
1450 
1451  ! 1.0 0.0, 3.0 1.0
1452  real function convect_qpf_map(v)
1453  implicit none
1454  real, intent(in) :: v
1455  if(v <= 1.0) then
1456  convect_qpf_map = 0
1457  else if(v <= 3.0) then
1458  convect_qpf_map = (v - 1.) / 2.
1459  else
1460  convect_qpf_map = 1.
1461  end if
1462  end function convect_qpf_map
1463 
1464  ! 1000.0 0.0, 2500.0 1.0
1465  real function convect_cape_map(v)
1466  implicit none
1467  real, intent(in) :: v
1468 
1469  if (v <= 1000.0) then
1470  convect_cape_map = 0.0
1471  else if(v <= 2500.) then
1472  convect_cape_map = (v - 1000.0) / 1400.
1473  else
1474  convect_cape_map = 1.0
1475  end if
1476  end function convect_cape_map
1477 
1478 
1479  ! -10.0 1.0, 0.0 0.0
1480  real function convect_liftedidx_map(v)
1481  implicit none
1482  real, intent(in) :: v
1483  if(v <= -10.) then
1484  convect_liftedidx_map = 1.0
1485  else if(v <= 0.) then
1486  convect_liftedidx_map = -0.1 * v
1487  else
1488  convect_liftedidx_map = 0.
1489  end if
1490  end function convect_liftedidx_map
1491 
1492 
1493  ! 20.0 0.0, 40.0 1.0
1494  real function convect_kidx_map(v)
1495  implicit none
1496  real, intent(in) :: v
1497  if(v <= 20.0) then
1498  convect_kidx_map = 0
1499  else if(v <= 40.0) then
1500  convect_kidx_map = (v - 20.) / 20.
1501  else
1502  convect_kidx_map = 1.
1503  end if
1504  end function convect_kidx_map
1505 
1506  ! 20.0 0.0, 55.0 1.0
1507  real function convect_totals_map(v)
1508  implicit none
1509  real, intent(in) :: v
1510  if(v <= 20.0) then
1511  convect_totals_map = 0
1512  else if( v <= 55.0)then
1513  convect_totals_map = (v - 20) / 35.
1514  else
1515  convect_totals_map = 1.0
1516  end if
1517  end function convect_totals_map
1518 
1519 end module severitymaps
1520 
1521 !========================================================================
1522 ! = = = = = = = = = = = = module IcingSeverity = = = = = = = = = = = = =
1523 !========================================================================
1524 module icingseverity
1525  use derivedfields, only : precips
1526  use cloudlayers, only : clouds_t
1527  use severitymaps
1528 
1529  implicit none
1530 
1531  private
1532  public icing_sev
1533 
1534  real, parameter :: cold_prcp_t_thld = 261.15
1535 
1536  ! module members: variables of cloud information
1537  ! original
1538  real :: ctt, layerq, avv
1539  integer :: nlayers, topidx, baseidx, wmnidx
1540  ! derived
1541  real :: deltaz, cloudtopdist
1542  ! whether it's the lowest cloud layer
1543  logical :: lowestcloud
1544 
1545  ! module member: derived variable
1546  real :: moistint
1547 
1548 contains
1549 
1550 !-----------------------------------------------------------------------+
1551  subroutine icing_sev(imp_physics,hgt, rh, t, pres, vv, liqCond, iceCond, twp, &
1552  ice_pot, nz, hcprcp, cape, lx, kx, tott, pc, prcptype, clouds, &
1553  iseverity)
1554  IMPLICIT NONE
1555  integer, intent(in) :: imp_physics
1556  integer, intent(in) :: nz
1557  real, intent(in) :: hgt(nz), rh(nz), t(nz), pres(nz), vv(nz)
1558  real, intent(in) :: liqcond(nz), icecond(nz), twp(nz), ice_pot(nz)
1559  real, intent(in) :: hcprcp, cape, lx, kx, tott, pc
1560  integer, intent(in) :: prcptype
1561  type(clouds_t), intent(in) :: clouds
1562  real, intent(out) :: iseverity(nz) ! category severity
1563 
1564  real :: severity
1565  integer :: k, n
1566 
1567  real :: moistint
1568 
1569  iseverity(:) = 0.0
1570 
1571  lowestcloud = .false.
1572 
1573  ! run the icing severity algorithm
1574  lp_n: do n = 1, clouds%nLayers
1575 
1576  ! extract cloud information
1577  wmnidx = clouds%wmnIdx
1578  avv = clouds%avv ! only for scenario below warmnose
1579  if (n == clouds%nLayers) lowestcloud = .true.
1580  ! apply the cloud top temperature to the cloud layer
1581  ctt = clouds%ctt(n)
1582  topidx = clouds%topIdx(n)
1583  baseidx = clouds%baseIdx(n)
1584 
1585  lp_k: do k = topidx, baseidx
1586  if (ice_pot(k) < 0.001) cycle
1587 
1588  severity = 0.0
1589 
1590  ! clouds information
1591  layerq = clouds%layerQ(k)
1592  ! derived
1593  deltaz = hgt(k) - hgt(baseidx)
1594  cloudtopdist = hgt(topidx) - hgt(k)
1595 
1596  ! derived variable
1597  if(imp_physics == 98 .or. imp_physics == 99) then
1598  moistint = moisture_map_cwat(rh(k), liqcond(k)+icecond(k), pres(k), t(k))
1599  else
1600  moistint = moisture_map_cond(rh(k), liqcond(k),icecond(k), pres(k), t(k))
1601  endif
1602 
1603  if(isconvection(prcptype)) then
1604  call convectionscenario &
1605  (t(k), hcprcp, cape, lx, kx, tott, severity)
1606  elseif(isclassicprcpblwwmn(prcptype, k)) then
1607  call classicprcpblwwmnscenario &
1608  (t(k), avv, pc, ice_pot(k), severity)
1609  elseif(isclassicprcpabvwmn(prcptype, k)) then
1610  call classicprcpabvwmnscenario &
1611  (twp(k), vv(k), t(k), pc, ice_pot(k),severity)
1612  elseif(isnoprecip(prcptype)) then
1613  call noprcpscenario &
1614  (vv(k), t(k), pc, ice_pot(k), severity)
1615  elseif(issnow(prcptype)) then
1616  call snowscenario &
1617  (twp(k), vv(k), t(k), pc, ice_pot(k), severity)
1618  elseif(iscoldrain(prcptype)) then
1619  call coldrainscenario &
1620  (twp(k), vv(k), t(k), pc, ice_pot(k), severity)
1621  elseif(iswarmprecip(prcptype)) then
1622  call warmprecipscenario &
1623  (twp(k), vv(k), t(k), pc, ice_pot(k), severity)
1624  elseif(isfreezingprecip(prcptype)) then
1625  call freezingprecipscenario &
1626  (twp(k), vv(k), t(k), pc, ice_pot(k), severity)
1627  end if
1628 
1629  ! severity category
1630  ! 0 = none (0, 0.08)
1631  ! 1 = trace [0.08, 0.21]
1632  ! 2 = light (0.21, 0.37]
1633  ! 3 = moderate (0.37, 0.67]
1634  ! 4 = heavy (0.67, 1]
1635  ! (0.0 0, 0.25 1, 0.425 2, 0.75 3, 1 4)
1636  ! (0.08 0, 0.21 1, 0.37 2, 0.67 3, 1 4) ! updated June 2015
1637 
1638  ! however the official categories are:
1639  ! 0 none
1640  ! 1 light
1641  ! 2 moderate
1642  ! 3 severe (no value)
1643  ! 4 trace
1644  ! 5 heavy
1645  !http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-207.shtml
1646 
1647  ! move category defination out of GFIP3.f to MDL2P.f - July 2015
1648 
1649  !if (severity < 0.08) then
1650  ! iseverity(k) = 0.0
1651  !elseif (severity <= 0.21) then
1652  ! iseverity(k) = 4.
1653  !else if(severity <= 0.37) then
1654  ! iseverity(k) = 1.0
1655  !else if(severity <= 0.67) then
1656  ! iseverity(k) = 2.0
1657  !else
1658  ! iseverity(k) = 5.0
1659  !endif
1660 
1661  ! make sure the values don't exceed 1.0
1662  iseverity(k)=min(1., severity)
1663  iseverity(k)=max(0., severity)
1664 
1665  end do lp_k
1666  end do lp_n
1667 
1668  return
1669  end subroutine icing_sev
1670 
1671 !-----------------------------------------------------------------------+
1672  logical function isconvection(prcpType)
1673  IMPLICIT NONE
1674  integer, intent(in) :: prcptype
1675 
1676  isconvection = prcptype == precips%CONVECTION
1677 
1678  return
1679  end function isconvection
1680 
1681 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +
1682  subroutine convectionscenario(t, hcprcp, cape, lx, kx, tott, severity)
1683  IMPLICIT NONE
1684  real, intent(in) :: t, hcprcp, cape, lx, kx, tott
1685  real, intent(out) :: severity
1686 
1687  integer :: scenario
1688 
1689  real :: tint, prcpint, capeint, lxint, kxint,tottint
1690 
1691  real :: weights(6) = (/ 4.0, 3.0, 5.0, 5.0, 2.0, 2.0 /)
1692  real :: sumweight
1693  sumweight = sum(weights)
1694 
1695  tint = convect_t_map(t) ** 0.5
1696  !Quantitative Precipitation Forecast
1697  prcpint = convect_qpf_map(hcprcp)
1698  capeint = convect_cape_map(cape)
1699  lxint = convect_liftedidx_map(lx)
1700  kxint = convect_kidx_map(kx)
1701  tottint = convect_totals_map(tott)
1702 
1703  scenario = scenarios%CONVECTION
1704 
1705  severity = (weights(1) * tint + weights(2) * prcpint + &
1706  weights(3) * capeint + weights(4) * lxint + &
1707  weights(5) * kxint + weights(6) * tottint) / &
1708  sumweight
1709  return
1710  end subroutine convectionscenario
1711 
1712 
1713 !-----------------------------------------------------------------------+
1714  logical function isclassicprcpblwwmn(prcpType, k)
1715  IMPLICIT NONE
1716  integer, intent(in) :: prcptype, k
1717  ! wmnIdx, topIdx, ctt: module members (cloud info)
1718 
1719  if ((prcptype /= precips%NONE .and. prcptype /= precips%SNOW) .and.&
1720  (wmnidx > 0 .and. k <= wmnidx .and. topidx > wmnidx) .and. &
1721  (ctt < cold_prcp_t_thld)) then
1722  isclassicprcpblwwmn = .true.
1723  else
1724  isclassicprcpblwwmn = .false.
1725  end if
1726 
1727  return
1728  end function isclassicprcpblwwmn
1729 
1730 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +
1731  subroutine classicprcpblwwmnscenario(t, vv, pc, ice_pot, severity)
1732  IMPLICIT NONE
1733  real, intent(in) :: t, pc, vv, ice_pot
1734  real, intent(out) :: severity
1735  ! deltaZ, ctt: module members (cloud info)
1736  ! moistInt: module member
1737 
1738  real :: tint, pcint, dzint
1739  real :: cttint, vvint
1740 
1741  integer :: scenario
1742 
1743  real :: weights(7) = (/ 3.0, 4.0, 3.0, 3.0, 3.5, 2.5, 2.0 /)
1744  real :: sumweight
1745  sumweight = sum(weights)
1746 
1747  scenario = scenarios%PRECIPITAION_BELOW_WARMNOSE
1748 
1749  tint = t_map(t, scenario)
1750  pcint = prcpcondensate_map(pc, scenario)
1751  dzint = deltaz_map(deltaz, scenario)
1752  cttint = ctt_map(ctt)
1753  vvint = vv_map(vv)
1754 
1755  severity = (weights(1) * tint + weights(2) * pcint + &
1756  weights(3) * dzint + weights(4) * cttint + &
1757  weights(5) * moistint + weights(6) * vvint + &
1758  weights(7) * ice_pot) / sumweight
1759 
1760  return
1761  end subroutine classicprcpblwwmnscenario
1762 
1763 
1764 !-----------------------------------------------------------------------+
1765 ! classical precipitation but not snow
1766  logical function isclassicprcpabvwmn(prcpType, k)
1767  IMPLICIT NONE
1768  integer, intent(in) :: prcptype, k
1769  ! wmnIdx,topIdx, ctt: module members (cloud info)
1770 
1771  if((prcptype /= precips%NONE .and. prcptype /= precips%SNOW) .and.&
1772  (wmnidx > 0 .and. k > wmnidx .and. topidx > wmnidx) .and. &
1773  (ctt < cold_prcp_t_thld)) then
1774  isclassicprcpabvwmn = .true.
1775  else
1776  isclassicprcpabvwmn = .false.
1777  end if
1778 
1779  return
1780  end function isclassicprcpabvwmn
1781 
1782 
1783 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +
1784  subroutine classicprcpabvwmnscenario(twp, vv, t, pc, ice_pot, severity)
1785  IMPLICIT NONE
1786  real, intent(in) :: twp, vv, t, pc, ice_pot
1787  real, intent(out) :: severity
1788  ! moistInt: module member
1789 
1790  real :: twpint, vvint
1791  real :: tempadj, cttadj, pcadj
1792 
1793  integer :: scenario
1794 
1795  real :: weights(4) = (/ 3.0, 3.5, 4.0, 2.0 /)
1796  real :: sumweight
1797  sumweight = sum(weights)
1798 
1799  scenario = scenarios%PRECIPITAION_ABOVE_WARMNOSE
1800 
1801  twpint = twp_map(twp, scenario)
1802  vvint = vv_map(vv)
1803 
1804  tempadj = 0.0
1805  cttadj = 0.0
1806  pcadj = 0.0
1807  ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info)
1808  call cal_dampingfactors(scenario, t, pc, tempadj, cttadj, pcadj)
1809 
1810  severity = (weights(1) * twpint + weights(2) * vvint + &
1811  weights(3) * moistint + weights(4) * ice_pot) / &
1812  (sumweight + 6.5*tempadj + 3.0*cttadj + 3.0*pcadj)
1813 
1814  return
1815  end subroutine classicprcpabvwmnscenario
1816 
1817 
1818 
1819 !-----------------------------------------------------------------------+
1820  logical function isnoprecip(prcpType)
1821  IMPLICIT NONE
1822  integer, intent(in) :: prcptype
1823 
1824  isnoprecip = prcptype == precips%NONE
1825 
1826  return
1827  end function isnoprecip
1828 
1829 
1830 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +
1831 
1832  subroutine noprcpscenario(vv, t, pc, ice_pot, severity)
1833  IMPLICIT NONE
1834  real, intent(in) :: vv, t, pc, ice_pot
1835  real, intent(out) :: severity
1836  ! layerQ, deltaZ: module members (cloud info)
1837  ! moistInt: module member
1838 
1839  real :: dqint, dzint, vvint
1840  real :: tempadj, cttadj, pcadj
1841 
1842  integer :: scenario
1843 
1844  real :: weights(5) = (/ 3.0, 3.5, 4.0, 4.0, 3.0 /)
1845  real :: sumweight
1846  sumweight = sum(weights)
1847 
1848  scenario = scenarios%NO_PRECIPITAION
1849 
1850  dqint = deltaq_map(layerq)
1851  dzint = deltaz_map(deltaz, scenario)
1852  vvint = vv_map(vv)
1853 
1854  tempadj = 0.0
1855  cttadj = 0.0
1856  pcadj = 0.0
1857  ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info)
1858  call cal_dampingfactors(scenario, t, pc, tempadj, cttadj, pcadj)
1859 
1860  severity = (weights(1) * dqint + weights(2) * dzint + &
1861  weights(3) * vvint + weights(4) * moistint + &
1862  weights(5) * ice_pot) / &
1863  (sumweight + 9.0*tempadj + 4.5*cttadj)
1864 
1865  return
1866  end subroutine noprcpscenario
1867 
1868 
1869 !-----------------------------------------------------------------------+
1870 
1871  logical function issnow(prcpType)
1872  IMPLICIT NONE
1873  integer, intent(in) :: prcptype
1874 
1875  issnow = prcptype == precips%SNOW
1876 
1877  return
1878  end function issnow
1879 
1880 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +
1881 
1882  subroutine snowscenario(twp, vv, t, pc, ice_pot, severity)
1883  IMPLICIT NONE
1884  real, intent(in) :: twp, vv, t, pc, ice_pot
1885  real, intent(out) :: severity
1886  ! deltaZ: module member (cloud info)
1887  ! moistInt: module member
1888 
1889  real :: twpint, dzint, vvint
1890  real :: tempadj, cttadj, pcadj
1891 
1892  integer :: scenario
1893 
1894  real :: weights(5) = (/ 3.0, 3.5, 3.5, 4.0, 3.0 /)
1895  real :: sumweight
1896  sumweight = sum(weights)
1897 
1898  scenario = scenarios%ALL_SNOW
1899 
1900  twpint = twp_map(twp, scenario)
1901  dzint = deltaz_map(deltaz, scenario)
1902  vvint = vv_map(vv)
1903 
1904  tempadj = 0.0
1905  cttadj = 0.0
1906  pcadj = 0.0
1907  ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info)
1908  call cal_dampingfactors(scenario, t, pc, tempadj, cttadj, pcadj)
1909 
1910  severity = (weights(1) * twpint + weights(2) * dzint + &
1911  weights(3) * vvint + weights(4) * moistint + &
1912  weights(5) * ice_pot) / &
1913  (sumweight + 9.0*tempadj + 4.0*cttadj + 4.0*pcadj)
1914 
1915  return
1916  end subroutine snowscenario
1917 
1918 
1919 !-----------------------------------------------------------------------+
1920  logical function iscoldrain(prcpType)
1921  IMPLICIT NONE
1922  integer, intent(in) :: prcptype
1923  ! ctt: module members (cloud info)
1924 
1925  iscoldrain = prcptype == precips%RAIN .and. ctt < cold_prcp_t_thld
1926 
1927  return
1928  end function iscoldrain
1929 
1930 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +
1931 
1932  subroutine coldrainscenario(twp, vv, t, pc, ice_pot, severity)
1933  IMPLICIT NONE
1934  real, intent(in) :: twp, vv, t, pc, ice_pot
1935  real, intent(out) :: severity
1936  ! deltaZ: module member (cloud info)
1937  ! moistInt: module member
1938 
1939  real :: twpint, dzint, vvint
1940  real :: tempadj, cttadj, pcadj
1941 
1942  integer :: scenario
1943 
1944  real :: weights(5) = (/ 3.0, 3.5, 3.5, 4.0, 2.0 /)
1945  real :: sumweight
1946  sumweight = sum(weights)
1947 
1948  scenario = scenarios%COLD_RAIN
1949 
1950  twpint = twp_map(twp, scenario)
1951  dzint = deltaz_map(deltaz, scenario)
1952  vvint = vv_map(vv)
1953 
1954  tempadj = 0.0
1955  cttadj = 0.0
1956  pcadj = 0.0
1957  ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info)
1958  call cal_dampingfactors(scenario, t, pc, tempadj, cttadj, pcadj)
1959 
1960  severity = (weights(1) * twpint + weights(2) * dzint + &
1961  weights(3) * vvint + weights(4) * moistint + &
1962  weights(5) * ice_pot) / &
1963  (sumweight + 8.0*tempadj + 4.0*cttadj + 4.0*pcadj)
1964 
1965  return
1966  end subroutine coldrainscenario
1967 
1968 
1969 !-----------------------------------------------------------------------+
1970 ! The elseif clause is a result of sloppy thinking on the
1971 ! part of the scientitsts. They were trying to create/use
1972 ! two different tests for the same 'scenario.' There is
1973 ! implicit definition of a 'classic precipitation' scenario.
1974  logical function iswarmprecip(prcpType)
1975  IMPLICIT NONE
1976  integer, intent(in) :: prcptype
1977  ! ctt, wmnIdx, topIdx: module members (cloud info)
1978 
1979  if((prcptype == precips%RAIN .or. prcptype == precips%OTHER) .and. &
1980  (ctt >= cold_prcp_t_thld)) then
1981  iswarmprecip = .true.
1982  elseif((prcptype /= precips%NONE .and. prcptype /= precips%SNOW).and.&
1983  (wmnidx > 0 .and. topidx <= wmnidx) .and. &
1984  (ctt >= cold_prcp_t_thld)) then
1985  iswarmprecip = .true.
1986  else
1987  iswarmprecip = .false.
1988  end if
1989 
1990  return
1991  end function iswarmprecip
1992 
1993 
1994 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +
1995  subroutine warmprecipscenario(twp, vv, t, pc, ice_pot, severity)
1996  IMPLICIT NONE
1997  real, intent(in) :: twp, vv, t, pc, ice_pot
1998  real, intent(out) :: severity
1999  ! deltaZ, layerQ: module member (cloud info)
2000  ! moistInt: module member
2001 
2002  real :: twpint, dzint, vvint, dqint
2003  real :: tempadj, cttadj, pcadj
2004 
2005  integer :: scenario
2006 
2007  real :: weights(6) = (/ 3.0, 3.0, 3.0, 3.5, 4.0, 2.0 /)
2008  real :: sumweight
2009  sumweight = sum(weights)
2010 
2011  scenario = scenarios%WARM_PRECIPITAION
2012 
2013  twpint = twp_map(twp, scenario)
2014  dzint = deltaz_map(deltaz, scenario)
2015  vvint = vv_map(vv)
2016  dqint = deltaq_map(layerq)
2017 
2018  tempadj = 0.0
2019  cttadj = 0.0
2020  pcadj = 0.0
2021  ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info)
2022  call cal_dampingfactors(scenario, t, pc, tempadj, cttadj, pcadj)
2023 
2024  severity = (weights(1) * twpint + weights(2) * dzint + &
2025  weights(3) * vvint + weights(4) * moistint + &
2026  weights(5) * ice_pot + weights(6) * dqint) / &
2027  (sumweight + 9.0*tempadj + 4.5*cttadj + 4.5*pcadj)
2028 
2029  return
2030  end subroutine warmprecipscenario
2031 
2032 
2033 !-----------------------------------------------------------------------+
2034  logical function isfreezingprecip(prcpType)
2035  IMPLICIT NONE
2036  integer, intent(in) :: prcptype
2037  ! ctt: module member (cloud info)
2038 
2039  if (prcptype == precips%OTHER .and. ctt < cold_prcp_t_thld) then
2040  isfreezingprecip = .true.
2041  else
2042  isfreezingprecip = .false.
2043  end if
2044 
2045  return
2046  end function isfreezingprecip
2047 
2048 
2049 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +
2050 
2051  subroutine freezingprecipscenario(twp, vv, t, pc, ice_pot, severity)
2052  IMPLICIT NONE
2053  real, intent(in) :: twp, vv, t, pc, ice_pot
2054  real, intent(out) :: severity
2055  ! deltaZ: module member (cloud info)
2056  ! moistInt: module member
2057 
2058  real :: twpint, dzint, vvint
2059  real :: tempadj, cttadj, pcadj
2060 
2061  integer :: scenario
2062 
2063  real :: weights(5) = (/ 3.0, 3.5, 3.5, 4.0, 2.0 /)
2064  real :: sumweight
2065  sumweight = sum(weights)
2066 
2067  scenario = scenarios%FREEZING_PRECIPITAION
2068 
2069  twpint = twp_map(twp, scenario)
2070  dzint = deltaz_map(deltaz, scenario)
2071  vvint = vv_map(vv)
2072 
2073  tempadj = 0.0
2074  cttadj = 0.0
2075  pcadj = 0.0
2076  ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info)
2077  call cal_dampingfactors(scenario, t, pc, tempadj, cttadj, pcadj)
2078 
2079  severity = (weights(1) * twpint + weights(2) * dzint + &
2080  weights(3) * vvint + weights(4) * moistint + &
2081  weights(5) * ice_pot) / &
2082  (sumweight + 8.0*tempadj + 4.0*cttadj + 4.0*pcadj)
2083 
2084  return
2085  end subroutine freezingprecipscenario
2086 
2087 
2088 !-----------------------------------------------------------------------+
2089 ! The only scenario NOT applicable to: below the warm nose
2090 !
2091 ! Temperature damping - can damp by a max of temp_adjust_factor
2092 !
2093 ! CTT damping - can damp by a max of ctt_adjust_factor
2094 ! The CTT is multiplied by a map based on the height below cloud top
2095 ! For all clouds and it's interest increases the closer to cloud
2096 ! top we are.
2097 !
2098 ! Condensate damping - can damp by a max of cond_adjust_factor
2099 ! Only in the lowest cloud layer and interest increases the closer to
2100 ! cloud base we are. Interest will be the max below cloud base as well.
2101 ! Maps differ for different scenarios
2102 !
2103  subroutine cal_dampingfactors(scenario, t, pc, tempAdj,cttAdj,condAdj)
2104  IMPLICIT NONE
2105  integer, intent(in) :: scenario
2106  real, intent(in) :: t, pc
2107  real, intent(out) :: tempadj, cttadj, condadj
2108  ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info)
2109 
2110  tempadj = t_map(t, scenario)
2111 
2112  cttadj = ctt_map(ctt) * cldtopdist_map(cloudtopdist)
2113 
2114  if (lowestcloud) then
2115  condadj = cldbasedist_map(deltaz)
2116  ! For no precipitation, condAdj is 0.0
2117  condadj = condadj * prcpcondensate_map(pc, scenario)
2118  end if
2119 
2120  return
2121  end subroutine cal_dampingfactors
2122 
2123 
2124 end module icingseverity
2125 
2126 !========================================================================
2127 ! = = = = = = = = = = = = = Icing Algorithm = = = = = = = = = = = = =
2128 !========================================================================
2129 subroutine icing_algo(i,j,pres,temp,rh,hgt,omega,wh,&
2130  q,cwat,qqw,qqi,qqr,qqs,qqg,&
2131  nz,xlat,xlon,xalt,prate,cprate,&
2132  cape,cin, ice_pot, ice_sev)
2133  use ctlblk_mod, only: imp_physics, spval, dtq2,me
2134  use physcons_post,only: g => con_g, fv => con_fvirt, rd => con_rd
2135 
2136  use derivedfields, only : derive_fields
2137  use cloudlayers, only : calc_cloudlayers, clouds_t
2138  use icingpotential, only : icing_pot
2139  use icingseverity, only : icing_sev
2140 
2141  IMPLICIT NONE
2142 
2143  integer, external :: gettopok
2144 !---------------------------------------------------------------------
2145 ! The GFIP algorithm computes the probability of icing within a model
2146 ! column given the follow input data
2147 !
2148 !
2149 ! 2-D data: total precip rate (prate)
2150 ! convective precip rate (cprate)
2151 ! the topography height (xalt)
2152 ! the latitude and longitude (xlat and xlon)
2153 ! the number of vertical levels (nz) = 47 in my GFS file
2154 ! ===== for severity only ======
2155 ! Convective Avail. Pot. Energy, pds7=65280 255-0 mb above
2156 ! gnd (cape)
2157 ! Convective inhibition, pds7=65280 255-0 mb above gnd (cin)
2158 ! 3-D data
2159 ! pressure(nz) in PA
2160 ! temperature(nz) in K
2161 ! rh(nz) in %
2162 ! hgt(nz) in GPM
2163 ! cwat(nz) in kg/x kg
2164 ! vv(nz) in Pa/sec
2165 !-----------------------------------------------------------------------
2166 !
2167 !-----------------------------------------------------------------------+
2168 ! This subroutine calculates the icing potential for a GFS column
2169 ! First the topography is mapped to the model's vertical coordinate
2170 ! in (topoK)
2171 ! Then derive related fields and cloud layers,
2172 ! up to 12 layers are possible
2173 ! The icing are computed in (icing_pot, ice_sev):
2174 ! The icing potential output should range from 0 - 1.
2175 ! The icing severity is in 4 categories: 1 2 3 4.
2176 !
2177 !-----------------------------------------------------------------------+
2178  integer, intent(in) :: i,j, nz
2179  real, intent(in),dimension(nz) :: pres,temp,rh,hgt,omega,wh
2180  ! q is used only for converting wh to omega
2181  real, intent(in),dimension(nz) :: q,cwat,qqw,qqi,qqr,qqs,qqg
2182  real, intent(in) :: xlat, xlon, xalt ! locations
2183  real, intent(in) :: prate, cprate ! precipitation rates
2184  real, intent(in) :: cape, cin
2185  real, intent(out) :: ice_pot(nz), ice_sev(nz)
2186 
2187  real :: hprcp, hcprcp
2188  integer :: topok, region, prcptype
2189  real, allocatable :: vv(:), ept(:), wbt(:), twp(:)
2190  real, allocatable :: icecond(:),liqcond(:), totwater(:),totcond(:)
2191  real :: pc, kx, lx, tott
2192  type(clouds_t) :: clouds
2193 
2194  integer :: k
2195 
2196  real :: tv
2197 
2198  allocate(vv(nz))
2199  allocate(ept(nz))
2200  allocate(wbt(nz))
2201  allocate(twp(nz))
2202  allocate(icecond(nz),liqcond(nz))
2203  allocate(totwater(nz),totcond(nz))
2204 
2205  allocate(clouds%layerQ(nz))
2206 
2207 ! write(*,'(2x,I3,A,1x,A,3I6,6F8.2)')me,":","iphy,i,j,lat,lon,prec,cprate,cape,cin=",imp_physics,i,j,xlat,xlon,prate,cprate,cape,cin
2208 ! do k=1,nz
2209 ! write(*,'(2x,I3,A,1x,I2,12F9.2)')me,":",k,pres(k),temp(k),rh(k),hgt(k),wh(k),q(k),cwat(k),qqw(k),qqi(k),qqr(k),qqs(k),qqg(k)
2210 ! end do
2211 
2212 ! if(mod(i,300)==0 .and. mod(j,200)==0)then
2213 ! print*,'sample input to FIP ',i,j,nz,xlat,xlon,xalt,prate, cprate
2214 ! do k=1,nz
2215 ! print*,'k,P,T,RH,H,CWM,VV',k,pres(k),temp(k),rh(k),hgt(k),cwat(k),vv(k)
2216 ! end do
2217 ! end if
2218 
2219  topok = gettopok(hgt,xalt,nz)
2220 
2221  ! hourly accumulated precipitation
2222  hprcp = prate * 1000. / dtq2 * 3600. ! large scale total precipitation in 1 hour
2223  hcprcp = cprate * 1000. / dtq2 * 3600. ! large scale convective precipitation in 1 hour
2224 
2225  if(imp_physics == 98 .or. imp_physics == 99) then
2226  do k = 1, nz
2227 
2228  ! input is omega (pa/s)
2229  vv(k) = omega(k)
2230 
2231  if(cwat(k) == spval) then
2232  liqcond(k) = spval
2233  icecond(k) = spval
2234  totwater(k) = spval
2235  totcond(k) = spval
2236  cycle
2237  endif
2238  if(temp(k) >=259.15) then ! T>=-14C
2239  liqcond(k) = cwat(k) * 1000.
2240  icecond(k) = 0.
2241  else
2242  liqcond(k) = 0.
2243  icecond(k) = cwat(k) * 1000.
2244  end if
2245 
2246  totwater(k) = cwat(k) * 1000.
2247 
2248  totcond(k) = cwat(k) * 1000.
2249  end do
2250  else if(imp_physics == 11 .or. imp_physics == 8) then
2251  do k = 1, nz
2252 
2253  ! convert input w (m/s) to omega (pa/s)
2254  vv(k) = wh(k)
2255  if(vv(k) /= spval) then
2256  tv = temp(k)*(1.+q(k)*fv)
2257  vv(k)=-vv(k)*g*(pres(k)/(rd*tv))
2258  end if
2259 
2260  if(qqw(k) /= spval .and. qqr(k) /= spval) then
2261  liqcond(k) = (qqw(k) + qqr(k)) * 1000.
2262  else
2263  liqcond(k) = spval
2264  end if
2265 
2266  if(qqi(k) /= spval .and. qqs(k) /= spval .and. qqg(k) /= spval) then
2267  icecond(k) = (qqi(k) + qqs(k) + qqg(k)) * 1000.
2268  else
2269  icecond(k) = spval
2270  end if
2271 
2272  if(liqcond(k) /= spval .and. icecond(k) /= spval) then
2273  totwater(k) = liqcond(k) + icecond(k)
2274  else
2275  totwater(k) = spval
2276  end if
2277 
2278  if(qqr(k) /= spval .and. qqs(k) /= spval .and. qqg(k) /= spval) then
2279  totcond(k) = (qqr(k) + qqs(k) + qqg(k)) * 1000.
2280  else
2281  totcond(k) = spval
2282  end if
2283  end do
2284  else
2285  print *, "This schema is not supported. imp_physics=", imp_physics
2286  return
2287  end if
2288 
2289  call derive_fields(imp_physics,temp, rh, pres, hgt, totwater,totcond, &
2290  nz, topok, hprcp, hcprcp, cin, cape,&
2291  ept, wbt, twp, pc, kx, lx, tott, prcptype)
2292 
2293  call calc_cloudlayers(rh, temp, pres, ept, vv, nz,topok, xlat,xlon,&
2294  region, clouds)
2295 
2296 ! write(*,'(2x,I3,A,1x,A,2I6,2F8.2)')me,":","i,j,lat,lon=",i,j,xlat,xlon
2297 ! do k=1,8
2298 ! write(*,'(2x,I3,A,1x,8F9.2)')me,":",pres((k-1)*8+1:(k-1)*8+8)
2299 ! end do
2300 
2301  call icing_pot(hgt, rh, temp, liqcond, vv, nz, clouds, ice_pot)
2302 
2303 
2304  call icing_sev(imp_physics, hgt, rh, temp, pres, vv, liqcond, icecond, twp, &
2305  ice_pot, nz, hcprcp, cape, lx, kx, tott, pc, prcptype, clouds, &
2306  ice_sev)
2307 ! if(mod(i,300)==0 .and. mod(j,200)==0)then
2308 ! print*,'FIP: cin,cape,pc, kx, lx, tott,pcpTyp',cin,cape,pc, kx, lx, tott,prcpType
2309 ! print *, "FIP i,j,lat,lon=",i,j,xlat,xlon
2310 ! do k=nz/2,nz
2311 ! print*,'k,ept,wbt,twp,totalwater=',k,ept(k), wbt(k), twp(k),&
2312 ! totWater(k), ice_pot(k), ice_sev(k)
2313 ! print*,'clouds nLayers wmnIdx avv',clouds%nLayers,clouds%wmnIdx,clouds%avv
2314 ! if (clouds%nLayers> 0) then
2315 ! print*,'clouds topidx=', clouds%topIdx,'baseidx=',clouds%baseIdx,&
2316 ! 'ctt=', clouds%ctt
2317 ! end if
2318 ! end do
2319 ! end if
2320 
2321  deallocate(vv)
2322  deallocate(ept)
2323  deallocate(wbt)
2324  deallocate(twp)
2325  deallocate(icecond,liqcond)
2326  deallocate(totwater,totcond)
2327  deallocate(clouds%layerQ)
2328 
2329  return
2330 end subroutine icing_algo
2331 
2332 !-----------------------------------------------------------------------+
2333 integer function gettopok(hgt, alt, nz)
2334  IMPLICIT NONE
2335  real, intent(in) :: hgt(nz)
2336  real, intent(in) :: alt
2337  integer, intent(in) :: nz
2338 
2339  integer :: k
2340 
2341  if(hgt(nz) >= alt) then
2342  gettopok = nz
2343  else
2344  do k=nz,2,-1
2345  if ((hgt(k-1) > alt) .and. (hgt(k) <= alt)) then
2346  gettopok = k-1
2347  exit
2348  endif
2349  end do
2350  end if
2351 
2352  return
2353 end function gettopok
physcons_post
Definition: physcons.f:1