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