UPP v11.0.0
Loading...
Searching...
No Matches
GFIP3.f
1!========================================================================
2! = = = = = = = = = = = = module DerivedFields = = = = = = = = = = = =
3!========================================================================
4module 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
22contains
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
577end module derivedfields
578
579!========================================================================
580! = = = = = = = = = = = = = module CloudLayers = = = = = = = = = = = = =
581!========================================================================
582module 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
606contains
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
915end module cloudlayers
916
917
918!========================================================================
919! = = = = = = = = = = = = module IcingPotential = = = = = = = = = = = = =
920!========================================================================
921module 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
933contains
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
1097end module icingpotential
1098
1099
1100!========================================================================
1101! = = = = = = = = = = = = module SeverityMaps = = = = = = = = = = = = =
1102!========================================================================
1103module 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
1120contains
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
1519end module severitymaps
1520
1521!========================================================================
1522! = = = = = = = = = = = = module IcingSeverity = = = = = = = = = = = = =
1523!========================================================================
1524module 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
1548contains
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
2124end module icingseverity
2125
2126!========================================================================
2127! = = = = = = = = = = = = = Icing Algorithm = = = = = = = = = = = = =
2128!========================================================================
2129subroutine 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
2330end subroutine icing_algo
2331
2332!-----------------------------------------------------------------------+
2333integer 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
2353end function gettopok
cloudlayers::clouds_t
Definition GFIP3.f:593
severitymaps::scenarios_t
Definition GFIP3.f:1108