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