!======================================================================== ! = = = = = = = = = = = = module DerivedFields = = = = = = = = = = = = !======================================================================== module DerivedFields private public derive_fields, mixing_ratio public PRECIPS ! Precipitation types type :: precipitations_t integer :: NONE = 0 integer :: RAIN = 1 integer :: SNOW = 2 integer :: OTHER = 3 integer :: CONVECTION = 4 end type precipitations_t type(precipitations_t), parameter :: PRECIPS = precipitations_t() contains !-----------------------------------------------------------------------+ subroutine derive_fields(imp_physics,t, rh, pres, hgt, totalWater, totalCond,& nz, topoK, hprcp, hcprcp, cin, cape, & ept, wbt, twp, pc, kx, lx, tott, prcpType) IMPLICIT NONE ! 3-D derived data: ! ept - equivalent potential temperature ! wbt - wet bulb temperature ! twp - total water path ! ! 2-D derived data: (indice for convective icing severity) ! kx - k index ! lx - lifted index ! tott - total totals ! ! 2-D derived data: ! pc - precipitation condensate ! prcpType - surface precipitation type integer, intent(in) :: imp_physics integer, intent(in) :: nz, topoK real, intent(in),dimension(nz) :: t, rh, pres, hgt real, intent(in),dimension(nz) :: totalWater,totalCond real, intent(in) :: hprcp, hcprcp, cin, cape real, intent(out) :: ept(nz), wbt(nz), twp(nz) real, intent(out) :: pc, kx, lx, tott integer, intent(out) :: prcpType real, allocatable :: td(:), tlc(:), wvm(:) integer :: k allocate(td(nz)) allocate(tlc(nz)) allocate(wvm(nz)) td(:) = getDewPointTemp(t(:), rh(:)) tlc(:) = get_tLCL(t(:), td(:)) wvm(:) = mixing_ratio(td(:), pres(:)) ept(:) = getThetae(pres(:), t(:), wvm(:), tlc(:)) wbt(:) = getThetaw(t(:), td(:), pres(:)) call calc_totalWaterPath(hgt, pres, t, totalWater, nz, twp) pc = getPrecipCond(totalCond, nz, topoK) ! indice for convective icing severity call calc_indice(t, td, pres, wvm, nz, topoK, kx, lx, tott) prcpType=getPrecipType(imp_physics,pres,hgt,t,rh,wbt,nz,& hprcp,hcprcp,pc,cin,cape,lx) deallocate(td) deallocate(tlc) deallocate(wvm) return end subroutine derive_fields !-----------------------------------------------------------------------+ ! dew point temperature in K ! t in K elemental real function getDewPointTemp(t, rh) IMPLICIT NONE real, intent(in) :: t, rh real vapr, rm ! actual water vapor presure in pascal ! 611.2 Pa is the saturated vapor presure at 0°C ! Pa = (rh/100) * Ps vapr = (max(1.0e-6,rh) / 100.0) * getVaporPres(t) rm = log(vapr/611.2) getDewPointTemp = 243.5 * rm / (17.67 - rm) + 273.15 return end function getDewPointTemp !-----------------------------------------------------------------------+ ! Equivalent potential temperature in K ! pres in Pa ! t in K ! wvm in kg/kg elemental real function getThetae(pres, t, wvm, tAtLCL) IMPLICIT NONE real, intent(in) :: pres, t, wvm, tatlcl real rmix, e, thtam rmix = max(0.0,wvm) ! in g/kg ! potential temperature e = (2.0/7.0) * ( 1.0 - 0.28 * rmix * 0.001) thtam = t * ((100000.0/pres)**e) getThetae = thtam * exp( (3.376/tAtLCL - 0.00254) * & (rmix * (1. + (.81 * rmix * 0.001)) )) return end function getThetae !-----------------------------------------------------------------------+ ! saturation vapor presure in Pa ! t in K elemental real function getVaporPres(t) IMPLICIT NONE real, intent(in) :: t real tc ! 611.2 Pa is the saturated vapor presure at 0°C tc = t-273.15 getVaporPres = 611.2 * exp( 17.67*tc/(tc+243.5)) return end function getVaporPres !-----------------------------------------------------------------------+ ! lifted condensation level temperature in K ! t and td in K elemental real function get_tLCL(t, td) IMPLICIT NONE real, intent(in) :: t, td get_tLCL = 1. / (1./(td - 56.) + log(t/td)/800.0) + 56.0 return end function get_tLCL !-----------------------------------------------------------------------+ ! mixing ratio in g/kg = water vapr/dry air ! td in K ! pres in Pa elemental real function mixing_ratio(td, pres) IMPLICIT NONE real, intent(in) :: td, pres real corr, e corr = 1.001 + ( (pres/100.0 - 100) /900.) * 0.0034 e = getVaporPres(td) * corr mixing_ratio = 0.62197 * (e/(pres-e)) * 1000.0 return end function mixing_ratio !-----------------------------------------------------------------------+ ! wet bulb temperature in K ! t and td in K ! pres in Pa elemental real function getThetaw(t, td, pres) IMPLICIT NONE real, intent(in) :: t, td, pres real :: Vl, Cp, dt, top, bottom, twet, twetNew real :: rmixr, smixr integer :: i Vl = 2.5e6 Cp = 1004.0 dt = 9.9e9 top = t bottom = td do i = 1, 100 twet = (top + bottom) / 2.0 rmixr = mixing_ratio(td, pres) / 1000.0 smixr = mixing_ratio(twet, pres) / 1000.0 twetNew = t - (Vl/Cp) * (smixr-rmixr) if(twetNew < twet) then top = twet else bottom = twet end if dt = abs(twet - twetNew) if (dt <= 0.1) exit end do ! assign a value anyway, don't need (dt < 1.) getThetaw = twet return end function getThetaw !-----------------------------------------------------------------------+ ! Precipitation Condensate in g/kg real function getPrecipCond(totalCond, nz, topoK) IMPLICIT NONE real, intent(in) :: totalCond(nz) integer, intent(in) :: nz, topoK integer :: k getPrecipCond = 0.0 do k = topoK, topoK-2, -1 ! three levels getPrecipCond = totalCond(k) + getPrecipCond end do return end function getPrecipCond !-----------------------------------------------------------------------+ ! These 2-D indice are used for convective icing severity subroutine calc_indice(t, td, pres, wvm, nz, topoK, & kIndex, liftedIndex, totalTotals) IMPLICIT NONE integer, intent(in) :: nz, topoK real, intent(in) :: t(nz), td(nz), pres(nz), wvm(nz) real, intent(out) :: kIndex, liftedIndex, totalTotals integer :: k real :: t500hPa, t700hPa, t850hPa, dpt700hPa, dpt850hPa real :: surfaceTemp,surfaceDewPtTemp,surfacePressure,surfaceMIXR real :: tempAtLCL, theta, pressAtLCL, thetaEAtLCL, tempFromThetaE, tem ! write(0,*)' nz=',nz,' pres=',pres(:) t500hPa = t(nz) t700hPa = t(nz) dpt700hPa = td(nz) t850hPa = t(nz) dpt850hPa = td(nz) ! do k = nz, 2, -1 ! use linear interpolation ! write(0,*)'k=',k,' pres=',pres(k) if ((pres(k)- 50000.0 >= 0.) .and. (pres(k-1)- 50000.0 < 0.) ) then if (abs(pres(k)- 50000.0) <= 0.1) then t500hPa = t(k) elseif (abs(pres(k-1)- 50000.0) <= 0.1) then t500hPa = t(k-1) else t500hPa = t(k) - ((pres(k)-50000.0)/(pres(k)-pres(k-1))) * (t(k)-t(k-1)) end if exit ! from surface to space, time to break the loop end if if ((pres(k)- 70000.0 >= 0.) .and. (pres(k-1)- 70000.0 < 0.) ) then if (abs(pres(k)- 70000.0) <= 0.1) then t700hPa = t(k) dpt700hPa = td(k) elseif (abs(pres(k-1)- 70000.0) <= 0.1) then t700hPa = t(k-1) dpt700hPa = td(k-1) else tem = (pres(k)-70000.0)/(pres(k)-pres(k-1)) t700hPa = t(k) - tem * (t(k)-t(k-1)) dpt700hPa = td(k) - tem * (td(k)-td(k-1)) end if endif ! write(0,*)'k=',k,' pres=',pres(k),pres(k-1) if ((pres(k)- 85000.0 >= 0.) .and. (pres(k-1)- 85000.0 < 0.) ) then if (abs(pres(k)- 85000.0) <= 0.1) then t850hPa = t(k) dpt850hPa = td(k) elseif (abs(pres(k-1)- 85000.0) <= 0.1) then t850hPa = t(k-1) dpt850hPa = td(k-1) else tem = (pres(k)-85000.0)/(pres(k)-pres(k-1)) t850hPa = t(k) - tem * (t(k)-t(k-1)) dpt850hPa = td(k) - tem * (td(k)-td(k-1)) end if endif end do ! kindex in C kIndex = (t850hPa-t500hPa) + (dpt850hPa-273.15) - (t700hPa-dpt700hPa) ! total total index totalTotals = t850hPa + dpt850hPa - t500hPa * 2.0 ! lifted index ! use topoK instead of 1 - Y Mao surfaceTemp = t(topoK) surfaceDewPtTemp = td(topoK) surfacePressure = pres(topoK) surfaceMIXR = wvm(topoK) tempAtLCL = get_tLCL(surfaceTemp, surfaceDewPtTemp) theta = surfaceTemp * (100000.0/surfacePressure)**0.286 pressAtLCL = 100000.0 * (tempAtLCL / theta)**3.4965 thetaEAtLCL = getThetae(pressAtLCL,tempAtLCL,surfaceMIXR,tempAtLCL) tempFromThetaE = getTemperature(thetaEAtLCL, 50000.0) liftedIndex = t500hPa - tempFromThetaE return end subroutine calc_indice !-----------------------------------------------------------------------+ real function getTemperature(thetaE, pres) IMPLICIT NONE real, intent(in) :: thetaE, pres real :: guess, w1, w2, thetaE1, thetaE2, cor integer :: i guess = (thetaE - 0.5 * (max(thetaE - 270.0, 0.0))**1.05) * & (pres/100000.0)**0.2 do i = 1, 100 w1 = calcRSubS(pres, guess) w2 = calcRSubS(pres, guess + 1.0) thetaE1 = getThetae(pres, guess, w1, guess) thetaE2 = getThetae(pres, guess + 1.0, w2, guess + 1.0) cor = (thetaE - thetaE1)/(thetaE2 - thetaE1) guess = guess + cor if (abs(cor) < 0.01) then exit end if end do getTemperature = guess return end function getTemperature !-----------------------------------------------------------------------+ elemental real function calcRSubS(pres, temp) IMPLICIT NONE real, intent(in) :: pres, temp real :: eSubs eSubS = calcESubS(temp) calcRSubS = (0.622*eSubS)/(pres - eSubS) return end function calcRSubS !-----------------------------------------------------------------------+ elemental real function calcESubS(temp) IMPLICIT NONE real, intent(in) :: temp real :: x real, parameter :: coeffs(9) = & (/ 610.5851, 44.40316, 1.430341, 0.02641412, 2.995057e-4,& 2.031998e-6, 6.936113e-9, 2.564861e-12, -3.704404e-14 /) x= max(-80.0, temp - 273.15) calcESubS = coeffs(1) + x*(coeffs(2) + x*(coeffs(3) + x*(coeffs(4) + & x*(coeffs(5) + x*(coeffs(6) + x*(coeffs(7) + & x*(coeffs(8) + x*coeffs(9)))))))) return end function calcESubS !-----------------------------------------------------------------------+ subroutine calc_totalWaterPath(hgt, pres, t, totalWater, nz, twp) IMPLICIT NONE integer, intent(in) :: nz real, intent(in) :: hgt(nz), pres(nz), t(nz), totalWater(nz) real, intent(out):: twp(nz) real :: condensateIntegrated integer :: topIdx, baseIdx integer :: k, m lp100: do k = 1, nz twp(k) = 0.0 topIdx = -1 baseIdx = -1 ! get the layer top and base for the total condensate layer lp200: do m = k, 2, -1 if (totalWater(m) > 0.001) then topIdx = m else exit end if end do lp200 lp300: do m = k, nz if (totalWater(m) > 0.001) then baseIdx = m else exit end if end do lp300 ! get the integrated condensate from the rep to the layer top if(topIdx == -1 .or. baseIdx == -1) cycle condensateIntegrated = 0.0 lp400: do m = topIdx, baseIdx if (t(m) <= 273.15) then condensateIntegrated = condensateIntegrated + & ( ((totalWater(m)*pres(m)) / (287.05 * t(m)) ) * & (hgt(m-1) - hgt(m) + 0.001)) end if end do lp400 twp(k) = condensateIntegrated end do lp100 return end subroutine calc_totalWaterPath !-----------------------------------------------------------------------+ integer function getPrecipType(imp_physics,pres, hgt, t, rh, wbt, nz, & hprcp, hcprcp, pc, cin, cape, lx) IMPLICIT NONE integer, intent(in) :: imp_physics integer, intent(in) :: nz real, intent(in) :: pres(nz), hgt(nz), t(nz), rh(nz), wbt(nz) real, intent(in) :: hprcp,hcprcp, pc, cin, cape, lx integer, allocatable :: wxType(:) integer :: idxMelt, idxRefz, iceidx real :: coldtemp, tColdArea, wetBuldArea real :: precip, precip_standard integer :: k, m allocate(wxType(nz)) idxMelt = nz ! melting idxRefz = nz ! refreezing if(imp_physics == 98 .or. imp_physics == 99) then precip = hprcp precip_standard =0.045 else if(imp_physics == 11 .or. imp_physics == 8) then precip = pc precip_standard = 0.01 end if ! Check for convection first if_conv: if ( hcprcp >= 1.0 .and. cape >= 1000.0 .and. cin > -100. & .and. lx < 0.0) then getPrecipType = PRECIPS% CONVECTION else ! find the height(index) where the snow melts, wetBulbTemp > 0C lp100: do k = 1, nz if( wbt(k) > 273.15) then idxMelt = k ! now look for a refreezing layer below warmnose do m = idxMelt, nz if (t(m) < 273.15) then idxRefz = m exit end if end do exit end if end do lp100 ! find the level that is below 150hPa or ! the min and max wet bulb temperature lp200: do k = nz, 1, -1 wxType(k) = 0 if_pres_wbt: if ( pres(k) >= 15000. .or. & (wbt(k) >= 200. .and. wbt(k) <= 1000.)) then iceIdx = k coldTemp = t(k) do m = k, 1, -1 ! look for level below cloud_top_pressure and ! cloud_top_altitude ! look for ctt if( (pres(m)>=30000. .or. (hgt(m)>hgt(nz)+700.))& .and. (rh(m) > 90.) .and. (t(m) < coldTemp)) then coldTemp = t(m) iceIdx = m end if end do ! found a cloud -- identify the precipitaiton type if_iceIdx: if (iceIdx /= k) then ! Sum the pos. area between wetBulbTemp and ! the 0 deg isotherm below coldTemp tColdArea = 0.0 do m = k, iceIdx, -1 if (wbt(m) >= 273.15) then tColdArea = tColdArea + (wbt(m)-273.15) * & (hgt(m-1) - hgt(m)) end if end do ! sum the total area between the wet bulb T and the 0 C isotherm ! in the lowest precip_type_altitude wetBuldArea = 0.0 do m = iceIdx, nz if ( (hgt(m) <= hgt(nz)+1500) .and. & (m > idxRefz)) then wetBuldArea = wetBuldArea + (wbt(m)-273.15) * & (hgt(m-1) - hgt(m)) end if end do ! get the probable Precip type if (coldTemp > 265.15) then wxType(k) = PRECIPS% OTHER else if (tColdArea < 350.) then if(t(k) <= 273.15) then wxType(k) = PRECIPS% SNOW else wxType(k) = PRECIPS% RAIN end if else if (wetBuldArea <= -250.) then wxType(k) = PRECIPS% OTHER else if(t(k) <= 273.15) then wxType(k) = PRECIPS% OTHER else wxType(k) = PRECIPS% RAIN end if else wxType(k) = PRECIPS% NONE end if if_iceIdx end if if_pres_wbt end do lp200 getPrecipType = PRECIPS% NONE ! if(hprcp >= 0.045) then ! if(pc >= 0.01) then ! g/kg if(precip >= precip_standard) then do k = nz, nz-1, -1 if (wxType(k) > getPrecipType) then getPrecipType = wxType(k) end if end do end if end if if_conv deallocate(wxType) return end function getPrecipType end module DerivedFields !======================================================================== ! = = = = = = = = = = = = = module CloudLayers = = = = = = = = = = = = = !======================================================================== module CloudLayers use DerivedFields, only : mixing_ratio private public calc_CloudLayers public clouds_t integer, parameter :: MaxLayers = 30 type :: clouds_t ! 2-D integer :: nLayers integer :: wmnIdx ! warm nose index real :: avv ! average vertical velocity ! 3-D, on model levels of nz real, allocatable :: layerQ(:) ! 3-D, of cloud layers integer :: topIdx(MaxLayers) integer :: baseIdx(MaxLayers) real :: ctt(MaxLayers) end type clouds_t contains !-----------------------------------------------------------------------+ subroutine calc_CloudLayers(rh,t,pres,ept,vv, nz, topoK, xlat, xlon,& region, clouds) IMPLICIT NONE integer, intent(in) :: nz, topoK real, intent(in) :: rh(nz),t(nz),pres(nz), ept(nz), vv(nz) real, intent(in) :: xlat,xlon integer, intent(out) :: region type(clouds_t), intent(inout) :: clouds ! layerQ has been allocated real :: t_rh integer :: in_cld,cur_base, num_lyr integer :: k, kk, n, kstart ! get the global region and set the RH thresholds if (abs(xlat)<23.5) then t_rh = 80.0 region = 1 elseif ( (abs(xlat)>=23.5).and.(abs(xlat)<66)) then t_rh = 75.0 region =2 else t_rh = 70.0 region =2 endif ! loop from the top (start at 2 so we can have n+1) ) ! bottom is nz and top is 1 num_lyr = 0 in_cld = 0 cur_base = 1 ! identify the very top layer if rh starts at high value if((rh(1) >= t_rh) .and. (rh(2) >= t_rh))then num_lyr = num_lyr + 1 in_cld = 1 ! find the cloud top and ctt clouds%topIdx(num_lyr) = 1 clouds%ctt(num_lyr) = t(1) ! find the cloud base do kk=2,topoK-1 if ((rh(kk-1)>=t_rh).and.(rh(kk)=t_rh).and.(rh(k+1)>=t_rh)) then num_lyr = num_lyr + 1 in_cld = 1 ! find the cloud top and ctt clouds%topIdx(num_lyr) = k clouds%ctt(num_lyr) = t(k) ! find the cloud base do kk=clouds%topIdx(num_lyr),topoK-1 if ((rh(kk-1)>=t_rh).and.(rh(kk) 0) then clouds%avv = getAverageVertVel(t, vv, nz, & clouds%topIdx(num_lyr), clouds%baseIdx(num_lyr)) else clouds%avv = 0 end if ! for severity call calc_layerQ(t, rh, pres, ept, nz, clouds) return end subroutine calc_CloudLayers !-----------------------------------------------------------------------+ ! From top to the surface, find the warmnoseIndex if existing ! | ! |T<=FREEZING ! ----|------------------warmnoseIndex ! | | ! | T > FREEZING| ! -------------- ! |T<=FREEZING ! ____|________ ! ///////////// surface !-----------------------------------------------------------------------+ integer function getWarmnoseIdx(t, nz) IMPLICIT NONE integer, intent(in) :: nz real, intent(in) :: t(nz) logical :: aboveFreezing integer :: tmpIndex integer :: k aboveFreezing = .false. tmpIndex = 0 !It doesn't matter of the initialized value getWarmnoseIdx = -1 ! search from top of model down to the surface do k = 1, nz if(t(k) <= 273.15) then if (aboveFreezing) then ! above 0 then below 0 getWarmnoseIdx = tmpIndex ! once warmnoseIndex is set, we can't get back ! here. let's break out of loop exit end if !aboveFreezing aboveFreezing = .false. else if(.not. aboveFreezing) tmpIndex = k aboveFreezing = .true. end if end do return end function getWarmnoseIdx !-----------------------------------------------------------------------+ ! Calculate the average VV in the dendritic layer (-13 to -16) if it's in ! the lowest cloud (clouds%nLayers). If the -13 to -16 layer falls ! between 2 levels then use the average VV from the levels on either side ! real function getAverageVertVel(t,vv,nz, topIdx_lowest,baseIdx_lowest) IMPLICIT NONE integer, intent(in) :: nz real, intent(in) :: t(nz), vv(nz) integer, intent(in) :: topIdx_lowest,baseIdx_lowest integer :: numVertVel, k, start_base real :: sumVertVel sumVertVel = 0.0 numVertVel = 0 ! The following loop starts at least nz-1 as the lowest model level if(baseIdx_lowest == nz) then start_base = nz - 1 else start_base = baseIdx_lowest end if do k = start_base, topIdx_lowest, -1 if ( t(k) <= 260.15 .and. t(k) >= 257.15) then sumVertVel = sumVertVel + vv(k) numVertVel = numVertVel + 1 else if (t(k) < 257.15) then sumVertVel = sumVertVel + vv(k) + vv(k+1) numVertVel = 2 exit end if end do if (numVertVel == 0) then getAverageVertVel = 0.0 else getAverageVertVel = sumVertVel / numVertVel end if return end function getAverageVertVel !-----------------------------------------------------------------------+ subroutine calc_layerQ(t, rh, pres, ept, nz, clouds) IMPLICIT NONE integer, intent(in) :: nz real, intent(in) :: t(nz), rh(nz), pres(nz), ept(nz) type(clouds_t), intent(inout) :: clouds real :: mean_rh, te_map, rh_map, adj_Q real :: totalQ, testThetaE, Q_top, Q_base, Q_layer, delta_TE, sum_rh integer :: base_k, test_k, num_layers integer :: k, m, n, kk clouds%layerQ(:) = 0.0 ! loop through each layer lp_nlayers: do n = 1, clouds%nLayers lp_k_inlayer: do k = clouds%topIdx(n), clouds%baseIdx(n)-1 totalQ = 0.0 testThetaE = ept(k) ! get base layer k value (first more than 4K colder than at k) base_k = clouds%baseIdx(n) do m = k, nz-1 if((testThetaE-ept(m)>4.) .and. (clouds%baseIdx(n)<=m))then base_k = m - 1 exit end if end do ! calculate delta thetaE and deltaQ and deltaRH for the layer lp_kk: do kk = k, base_k-1 Q_top = mixing_ratio(t(kk), pres(kk)) Q_base = mixing_ratio(t(kk+1), pres(kk+1)) Q_layer = Q_base - Q_top ! convert Q from g/kg to g/m**3 Q_layer = Q_layer * (pres(kk)/(287.05 * t(kk))) if (Q_layer < 0.0) Q_layer = 0.0 delta_TE = testThetaE - ept(kk+1) test_k = kk + 1 ! get the mean RH in the layer num_layers = 0 sum_rh = 0.0 do m = k, test_k num_layers = num_layers + 1 sum_rh = sum_rh + rh(m) end do if (num_layers == 0) then mean_rh = 0.0 else mean_rh = sum_rh / num_layers end if te_map = dq_delta_TE_map(delta_TE) rh_map = dq_rh_map(mean_rh) adj_Q = Q_layer * te_map * rh_map totalQ = totalQ + adj_Q end do lp_kk clouds%layerQ(k) = totalQ end do lp_k_inlayer end do lp_nlayers return end subroutine calc_layerQ !-----------------------------------------------------------------------+ ! 70.0 0.0, 100.0 1.0 real function dq_rh_map(rh) IMPLICIT NONE real, intent(in) :: rh if(rh <= 70.) then dq_rh_map = 0. else if( rh >= 100.) then dq_rh_map = 1.0 else dq_rh_map = (rh-70.)/30. end if return end function dq_rh_map !-----------------------------------------------------------------------+ ! 0.0 1.0, 4.0 0.0 real function dq_delta_TE_map(delTE) IMPLICIT NONE real, intent(in) :: delTE if (delTE <= 0.0) then dq_delta_TE_map = 1.0 else if (delTE <= 4.0) then dq_delta_TE_map = (4. - delTE) / 4.0 else dq_delta_TE_map = 0.0 end if return end function dq_delta_TE_map end module cloudlayers !======================================================================== ! = = = = = = = = = = = = module IcingPotential = = = = = = = = = = = = = !======================================================================== module IcingPotential use ctlblk_mod, only: me use DerivedFields, only : PRECIPS use CloudLayers, only : clouds_t private public icing_pot contains !-----------------------------------------------------------------------+ subroutine icing_pot(hgt, rh, t, liqCond, vv, nz, clouds, ice_pot) IMPLICIT NONE integer, intent(in) :: nz real, intent(in) :: hgt(nz), rh(nz), t(nz), liqCond(nz), vv(nz) type(clouds_t), intent(in) :: clouds real, intent(out) :: ice_pot(nz) real :: ctt, slw, slw_fac, vv_fac integer :: k, n ! variables for printout only real :: mapt,maprh,mapctt,mapvv,mapslw ice_pot(:) = 0.0 ! run the icing potential algorithm lp_n: do n = 1, clouds%nLayers ! apply the cloud top temperature to the cloud layer ctt = clouds%ctt(n) lp_k: do k = clouds%topIdx(n), clouds%baseIdx(n) ! convert the liquid water to slw if the CTT > -14C if ( ctt>=259.15 .and. ctt<=273.15 ) then slw = liqCond(k) else slw = 0.0 end if ice_pot(k) = tmap(t(k)) * rh_map(rh(k)) * ctt_map(ctt) ! add the VV map if (vv_map(vv(k))>=0.0) then vv_fac = (1.0-ice_pot(k)) * (0.25*vv_map(vv(k))) else vv_fac = ice_pot(k) * (0.25*vv_map(vv(k))) endif ! add the slw slw_fac = (1.0 - ice_pot(k))*(0.4*slw_map(slw)) ! calculate the final icing potential if (ice_pot(k)>0.001) then ice_pot(k) = ice_pot(k) + vv_fac + slw_fac endif ! make sure the values don't exceed 1.0 if (ice_pot(k)>1.0) then ice_pot(k) = 1.0 endif ! mapt=tmap(t(k)) ! maprh=rh_map(rh(k)) ! mapctt=ctt_map(ctt) ! mapvv=vv_map(vv(k)) ! mapslw=slw_map(slw) ! write(*,'(2x,I3,A,1x,A,I3,F7.3)')me,":","k,icip=",k,ice_pot(k) ! write(*,'(2x,I3,A,1x,F8.2,F7.3,F7.3,F7.3)') & ! me,":",t(k),mapt,rh(k),maprh ! write(*,'(2x,I3,A,1x,F8.2,F7.3,F10.6,F7.3)') & ! me,":",ctt,mapctt,vv(k),mapvv ! write(*,'(2x,I3,A,1x,F11.6,F11.6,F7.3)') & ! me,":",liqCond(k), slw, mapslw end do lp_k end do lp_n return end subroutine icing_pot !-----------------------------------------------------------------------+ real function tmap(temp) IMPLICIT NONE real, intent(in) :: temp if((temp>=248.15).and.(temp<=263.15)) then tmap=((temp-248.15)/15.0)**1.5 elseif((temp>263.15).and.(temp<268.15)) then tmap=1.0 elseif((temp>268.15).and.(temp<273.15)) then tmap=((273.15 - temp)/5.0)**0.75 else tmap=0.0 endif return end function tmap !-----------------------------------------------------------------------+ real function rh_map(rh) IMPLICIT NONE real, intent(in) :: rh if (rh>95.0) then rh_map=1.0 elseif ( (rh<=95.0).and.(rh>=50.0) ) then rh_map=((20./9.) * ((rh/100.0)-0.5))**3.0 else rh_map=0.0 endif return end function rh_map !-----------------------------------------------------------------------+ real function ctt_map(ctt) IMPLICIT NONE real, intent(in) :: ctt if((ctt>=261.0).and.(ctt<280.)) then ctt_map = 1.0 elseif((ctt>223.0).and.(ctt<261.0 )) then ctt_map = 0.2 + 0.8*(((ctt - 223.0)/38.0)**2.0) elseif(ctt<223.0) then ctt_map = 0.2 else ctt_map = 0.0 endif return end function ctt_map !-----------------------------------------------------------------------+ real function vv_map(vv) IMPLICIT NONE real, intent(in) :: vv if (vv>0.0) then vv_map = 0.0 elseif (vv<-0.5) then vv_map = 1.0 else vv_map = -1.0 * (vv/0.5) endif return end function vv_map !-----------------------------------------------------------------------+ real function slw_map(slw) IMPLICIT NONE real, intent(in) :: slw if(slw>0.2) then slw_map = 1.0 elseif (slw<=0.001) then slw_map = 0.0 else slw_map = slw/0.2 endif return end function slw_map end module IcingPotential !======================================================================== ! = = = = = = = = = = = = module SeverityMaps = = = = = = = = = = = = = !======================================================================== module SeverityMaps public type :: scenarios_t integer :: NO_PRECIPITAION = 0 integer :: PRECIPITAION_BELOW_WARMNOSE = 1 integer :: PRECIPITAION_ABOVE_WARMNOSE = 2 integer :: ALL_SNOW = 3 integer :: COLD_RAIN = 4 integer :: WARM_PRECIPITAION = 5 integer :: FREEZING_PRECIPITAION = 6 integer :: CONVECTION = 7 end type scenarios_t type(scenarios_t), parameter :: SCENARIOS = scenarios_t() contains !-----------------------------------------------------------------------+ ! scenario dependant !-----------------------------------------------------------------------+ real function twp_map(v, scenario) real, intent(in) :: v integer, intent(in) :: scenario twp_map = 0.0 select case (scenario) case(SCENARIOS%PRECIPITAION_ABOVE_WARMNOSE, SCENARIOS%ALL_SNOW, & SCENARIOS%COLD_RAIN, SCENARIOS%FREEZING_PRECIPITAION) ! 0.0 0.0, 1000.0 1.0 if(v <= 0.0) then twp_map = 0.0 else if( v < 1000.) then twp_map = v / 1000.0 else twp_map = 1. end if case( SCENARIOS%WARM_PRECIPITAION) ! 0.0 0.0, 500.0 1.0 if(v <= 0.0) then twp_map = 0.0 else if( v < 500.) then twp_map = v / 500.0 else twp_map = 1. end if end select return end function twp_map ! Only precip below warmnose has a different temperature map real function t_map(v, scenario) real, intent(in) :: v integer, intent(in) :: scenario t_map = 0. select case (scenario) case(SCENARIOS%PRECIPITAION_BELOW_WARMNOSE) ! 269.15 1.0, 273.15 0.0 if(v <= 269.15) then t_map = 1. elseif( v <= 273.15) then t_map = (273.15 - v ) / 4.0 else t_map = 0. end if case default ! 248.15 1.0, 248.65 0.8039, 249.15 0.7226, 251.15 0.5196, ! 253.15 0.3798, 255.15 0.2662, 257.15 0.1679, 259.15 0.0801, ! 261.15 0.0, 268.15 0.0, 273.15 1.0 if(v <= 248.15) then t_map = 1.0 elseif( v <= 248.65) then t_map = (49.162215 - 0.1961*v) * 2. elseif( v <= 249.15) then t_map = (20.617195 - 0.0813*v) * 2. elseif(v <= 251.15) then t_map = (52.02265 - 0.203*v) / 2.0 elseif(v <= 253.15) then t_map = (36.14997 - 0.1398*v) / 2.0 elseif(v <= 255.15) then t_map = (29.51744 - 0.1136*v) / 2.0 elseif(v <= 257.15) then t_map = (25.613645-0.0983*v) / 2.0 elseif(v <= 259.15) then t_map = (22.91357 - 0.0878*v) / 2.0 elseif( v <= 261.15) then t_map = (20.918115 - 0.0801*v) / 2. else t_map = 0.0 end if end select end function t_map ! Condensates near the surface take place of radar reflectivity in CIP real function prcpCondensate_map(v, scenario) real, intent(in) :: v integer, intent(in) :: scenario prcpCondensate_map = 0.0 select case (scenario) case(SCENARIOS%NO_PRECIPITAION) ! do nothing case(SCENARIOS%PRECIPITAION_BELOW_WARMNOSE, SCENARIOS%COLD_RAIN, & SCENARIOS%PRECIPITAION_ABOVE_WARMNOSE) ! 0.05 0.0, 0.2 1.0 if(v <= 0.05) then prcpCondensate_map = 0. elseif(v <= 0.2) then prcpCondensate_map = (v - 0.05) / 0.15 else prcpCondensate_map = 1. end if case(SCENARIOS%ALL_SNOW) ! 0.05 0.0, 0.25 1.0 if(v <= 0.05) then prcpCondensate_map = 0. elseif(v <= 0.25) then prcpCondensate_map = (v - 0.05) / 0.2 else prcpCondensate_map = 1.0 end if case(SCENARIOS%WARM_PRECIPITAION, SCENARIOS%FREEZING_PRECIPITAION) ! 0.05 0.0, 0.15 0.5 if(v <= 0.05) then prcpCondensate_map = 0. elseif(v <= 0.15) then prcpCondensate_map = (.5*v - 0.025) / 0.1 else prcpCondensate_map = 0.5 end if end select return end function prcpCondensate_map real function deltaZ_map(v, scenario) real, intent(in) :: v integer, intent(in) :: scenario deltaZ_map = 0.0 select case (scenario) case (SCENARIOS%NO_PRECIPITAION, SCENARIOS%WARM_PRECIPITAION) ! 0.0 0.0, 1828.8 1.0 if(v <= 0.0) then deltaZ_map = 0.0 else if(v <= 1828.8) then deltaZ_map = v /1828.8 else deltaZ_map = 1. end if case(SCENARIOS%PRECIPITAION_BELOW_WARMNOSE) ! 0.0 0.0, 30.49999 0.0, 30.5 1.0 if(v < 30.5) then deltaZ_map = 0.0 else deltaZ_map = 1. end if case(SCENARIOS%ALL_SNOW) ! 914.4 0.0, 2743.2 1.0 if(v <= 914.4) then deltaZ_map = 0.0 else if(v <= 2743.2) then deltaZ_map = (v - 914.4) / 1828.8 else deltaZ_map = 1. end if case(SCENARIOS%COLD_RAIN, SCENARIOS%FREEZING_PRECIPITAION) ! 1524.0 0.0, 4267.2 1.0 if(v <= 1524.) then deltaZ_map = 0.0 else if(v <= 4267.2) then deltaZ_map = (v - 1524.) / 2743.2 else deltaZ_map = 1. end if end select return end function deltaZ_map !-----------------------------------------------------------------------+ ! scenario independant. !-----------------------------------------------------------------------+ ! 223.15 0.8, 233.15 0.7446, 243.15 0.5784, 253.15 0.3014 ! 261.15 0.0, 280.15 0.0, 280.151 1.0 real function ctt_map(v) real, intent(in) :: v if(v <= 223.15) then ctt_map = 0.8 elseif(v <= 233.15) then ctt_map = (20.36251 - 0.0554*v) / 10. elseif(v <= 243.15) then ctt_map = (46.19553 - 0.1662*v) / 10. elseif(v <= 253.15) then ctt_map = (73.13655 - 0.277*v) / 10. elseif(v <= 261.15) then ctt_map = (78.71061 - 0.3014*v) / 8. else ctt_map = 0. end if end function ctt_map ! -0.5 1.0, 0.0 0.0 real function vv_map(v) real, intent(in) :: v if(v <= -0.5) then vv_map = 1. else if(v <= 0.) then vv_map = -2. * v else vv_map = 0. end if end function vv_map ! cloud top distance ! 609.6 1.0, 3048.0 0.0 real function cldTopDist_map(v) if( v <= 609.6) then cldTopDist_map = 1.0 elseif( v <= 3048.) then cldTopDist_map = (3048. - v) / 2438.4 else cldTopDist_map = 0. end if end function cldTopDist_map ! cloud base distance ! 304.8 1.0, 1524.0 0.0 real function cldBaseDist_map(v) if( v <= 304.8) then cldBaseDist_map = 1.0 elseif( v <= 1524.) then cldBaseDist_map = (1524. - v) / 1219.2 else cldBaseDist_map = 0. end if end function cldBaseDist_map ! 0.0 0.0, 1.0 1.0 real function deltaQ_map(v) if( v <= 0.) then deltaQ_map = 0 elseif( v <= 1.0) then deltaQ_map = v else deltaQ_map = 1. end if end function deltaQ_map real function moisture_map_cond(rh, liqCond, iceCond, pres, t) IMPLICIT NONE real, intent(in) :: rh, liqCond, iceCond, pres, t real :: liq_m3, ice_m3, rhFactor, liqFactor,iceFactor ! Convert liqCond, iceCond from g/kg to g/m^3 liq_m3 = (liqCond * pres) / (287.05 * t) ice_m3 = (iceCond * pres) / (287.05 * t) rhFactor = 0.5 * rh_map(rh) liqFactor = 0.6* condensate_map(liq_m3) iceFactor = 0.4 * condensate_map(ice_m3) moisture_map_cond = rhFactor + (1.0 - rhFactor) * (liqFactor + iceFactor) return end function moisture_map_cond ! If not identify liquid/ice condensate real function moisture_map_cwat(rh, cwat, pres, t) IMPLICIT NONE real, intent(in) :: rh, cwat, pres, t real :: cwat_m3, rhFactor, cwatFactor ! Convert cwat from g/kg to g/m^3 cwat_m3 = (cwat * pres) / (287.05 * t) rhFactor = 0.5 * rh_map(rh) cwatFactor = condensate_map(cwat_m3) moisture_map_cwat = rhFactor + (1.0 - rhFactor) * cwatFactor return end function moisture_map_cwat ! only called by moisture_map ! 70.0 0.0, 100.0 1.0 real function rh_map(v) real, intent(in) :: v if(v <= 70.) then rh_map = 0. else if(v <= 100) then rh_map = (v - 70.) / 30. else rh_map = 1. end if end function rh_map ! only called by moisture_map ! 0.00399 0.0, 0.004 0.0, 0.2 1.0 real function condensate_map(v) real, intent(in) :: v if(v <= 0.004) then condensate_map = 0. elseif( v <= 0.2) then condensate_map = (v - 0.004) / 0.196 else condensate_map = 1. end if end function condensate_map !-----------------------------------------------------------------------+ ! convective !-----------------------------------------------------------------------+ ! 243.150 0.0, 265.15 1.0, 269.15 1.0, 270.15 0.87 ! 271.15 0.71, 272.15 0.50, 273.15 0.0 real function convect_t_map(v) real, intent(in) :: v if(v <= 243.15) then convect_t_map = 0. else if(v <= 265.15) then convect_t_map = (v -243.15)/22. else if(v <= 269.15) then convect_t_map = 1. else if(v <= 270.15) then convect_t_map = -0.13*v + 35.9895 else if(v <= 271.15) then convect_t_map = -0.16*v + 44.094 else if(v <= 272.15) then convect_t_map = -0.21*v + 57.6515 else if(v <= 273.15) then convect_t_map = -0.50*v + 136.575 else convect_t_map = 0. end if end function convect_t_map ! 1.0 0.0, 3.0 1.0 real function convect_qpf_map(v) real, intent(in) :: v if(v <= 1.0) then convect_qpf_map = 0 else if(v <= 3.0) then convect_qpf_map = (v - 1.) / 2. else convect_qpf_map = 1. end if end function convect_qpf_map ! 1000.0 0.0, 2500.0 1.0 real function convect_cape_map(v) real, intent(in) :: v if (v <= 1000.0) then convect_cape_map = 0.0 else if(v <= 2500.) then convect_cape_map = (v - 1000.0) / 1400. else convect_cape_map = 1.0 end if end function convect_cape_map ! -10.0 1.0, 0.0 0.0 real function convect_liftedIdx_map(v) real, intent(in) :: v if(v <= -10.) then convect_liftedIdx_map = 1.0 else if(v <= 0.) then convect_liftedIdx_map = -0.1 * v else convect_liftedIdx_map = 0. end if end function convect_liftedIdx_map ! 20.0 0.0, 40.0 1.0 real function convect_kIdx_map(v) real, intent(in) :: v if(v <= 20.0) then convect_kIdx_map = 0 else if(v <= 40.0) then convect_kIdx_map = (v - 20.) / 20. else convect_kIdx_map = 1. end if end function convect_kIdx_map ! 20.0 0.0, 55.0 1.0 real function convect_totals_map(v) real, intent(in) :: v if(v <= 20.0) then convect_totals_map = 0 else if( v <= 55.0)then convect_totals_map = (v - 20) / 35. else convect_totals_map = 1.0 end if end function convect_totals_map end module SeverityMaps !======================================================================== ! = = = = = = = = = = = = module IcingSeverity = = = = = = = = = = = = = !======================================================================== module IcingSeverity use DerivedFields, only : PRECIPS use CloudLayers, only : clouds_t use SeverityMaps private public icing_sev real, parameter :: COLD_PRCP_T_THLD = 261.15 ! module members: variables of cloud information ! original real :: ctt, layerQ, avv integer :: nLayers, topIdx, baseIdx, wmnIdx ! derived real :: deltaZ, cloudTopDist ! whether it's the lowest cloud layer logical :: lowestCloud ! module member: derived variable real :: moistInt contains !-----------------------------------------------------------------------+ subroutine icing_sev(imp_physics,hgt, rh, t, pres, vv, liqCond, iceCond, twp, & ice_pot, nz, hcprcp, cape, lx, kx, tott, pc, prcpType, clouds, & iseverity) IMPLICIT NONE integer, intent(in) :: imp_physics integer, intent(in) :: nz real, intent(in) :: hgt(nz), rh(nz), t(nz), pres(nz), vv(nz) real, intent(in) :: liqCond(nz), iceCond(nz), twp(nz), ice_pot(nz) real, intent(in) :: hcprcp, cape, lx, kx, tott, pc integer, intent(in) :: prcpType type(clouds_t), intent(in) :: clouds real, intent(out) :: iseverity(nz) ! category severity real :: severity integer :: k, n iseverity(:) = 0.0 lowestCloud = .false. ! run the icing severity algorithm lp_n: do n = 1, clouds%nLayers ! extract cloud information wmnIdx = clouds%wmnIdx avv = clouds%avv ! only for scenario below warmnose if (n == clouds%nLayers) lowestCloud = .true. ! apply the cloud top temperature to the cloud layer ctt = clouds%ctt(n) topIdx = clouds%topIdx(n) baseIdx = clouds%baseIdx(n) lp_k: do k = topIdx, baseIdx if (ice_pot(k) < 0.001) cycle severity = 0.0 ! clouds information layerQ = clouds%layerQ(k) ! derived deltaZ = hgt(k) - hgt(baseIdx) cloudTopDist = hgt(topIdx) - hgt(k) ! derived variable if(imp_physics == 98 .or. imp_physics == 99) then moistInt = moisture_map_cwat(rh(k), liqCond(k)+iceCond(k), pres(k), t(k)) else moistInt = moisture_map_cond(rh(k), liqCond(k),iceCond(k), pres(k), t(k)) endif if(isConvection(prcpType)) then call convectionScenario & (t(k), hcprcp, cape, lx, kx, tott, severity) elseif(isClassicPrcpBlwWmn(prcpType, k)) then call classicPrcpBlwWmnScenario & (t(k), avv, pc, ice_pot(k), severity) elseif(isClassicPrcpAbvWmn(prcpType, k)) then call classicPrcpAbvWmnScenario & (twp(k), vv(k), t(k), pc, ice_pot(k),severity) elseif(isNoPrecip(prcpType)) then call noPrcpScenario & (vv(k), t(k), pc, ice_pot(k), severity) elseif(isSnow(prcpType)) then call snowScenario & (twp(k), vv(k), t(k), pc, ice_pot(k), severity) elseif(isColdRain(prcpType)) then call coldRainScenario & (twp(k), vv(k), t(k), pc, ice_pot(k), severity) elseif(isWarmPrecip(prcpType)) then call warmPrecipScenario & (twp(k), vv(k), t(k), pc, ice_pot(k), severity) elseif(isFreezingPrecip(prcpType)) then call freezingPrecipScenario & (twp(k), vv(k), t(k), pc, ice_pot(k), severity) end if ! severity category ! 0 = none (0, 0.08) ! 1 = trace [0.08, 0.21] ! 2 = light (0.21, 0.37] ! 3 = moderate (0.37, 0.67] ! 4 = heavy (0.67, 1] ! (0.0 0, 0.25 1, 0.425 2, 0.75 3, 1 4) ! (0.08 0, 0.21 1, 0.37 2, 0.67 3, 1 4) ! updated June 2015 ! however the official categories are: ! 0 none ! 1 light ! 2 moderate ! 3 severe (no value) ! 4 trace ! 5 heavy !http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-207.shtml ! move category defination out of GFIP3.f to MDL2P.f - July 2015 !if (severity < 0.08) then ! iseverity(k) = 0.0 !elseif (severity <= 0.21) then ! iseverity(k) = 4. !else if(severity <= 0.37) then ! iseverity(k) = 1.0 !else if(severity <= 0.67) then ! iseverity(k) = 2.0 !else ! iseverity(k) = 5.0 !endif ! make sure the values don't exceed 1.0 iseverity(k)=min(1., severity) iseverity(k)=max(0., severity) end do lp_k end do lp_n return end subroutine icing_sev !-----------------------------------------------------------------------+ logical function isConvection(prcpType) IMPLICIT NONE integer, intent(in) :: prcpType isConvection = prcpType == PRECIPS%CONVECTION return end function isConvection ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + subroutine convectionScenario(t, hcprcp, cape, lx, kx, tott, severity) IMPLICIT NONE real, intent(in) :: t, hcprcp, cape, lx, kx, tott real, intent(out) :: severity integer :: scenario real :: tInt, prcpInt, capeInt, lxInt, kxInt,tottInt real :: weights(6) = (/ 4.0, 3.0, 5.0, 5.0, 2.0, 2.0 /) real :: sumWeight sumWeight = sum(weights) tInt = convect_t_map(t) ** 0.5 !Quantitative Precipitation Forecast prcpInt = convect_qpf_map(hcprcp) capeInt = convect_cape_map(cape) lxInt = convect_liftedIdx_map(lx) kxInt = convect_kIdx_map(kx) tottInt = convect_totals_map(tott) scenario = SCENARIOS%CONVECTION severity = (weights(1) * tInt + weights(2) * prcpInt + & weights(3) * capeInt + weights(4) * lxInt + & weights(5) * kxInt + weights(6) * tottInt) / & sumWeight return end subroutine convectionScenario !-----------------------------------------------------------------------+ logical function isClassicPrcpBlwWmn(prcpType, k) IMPLICIT NONE integer, intent(in) :: prcpType, k ! wmnIdx, topIdx, ctt: module members (cloud info) if ((prcpType /= PRECIPS%NONE .and. prcpType /= PRECIPS%SNOW) .and.& (wmnIdx > 0 .and. k <= wmnIdx .and. topIdx > wmnIdx) .and. & (ctt < COLD_PRCP_T_THLD)) then isClassicPrcpBlwWmn = .true. else isClassicPrcpBlwWmn = .false. end if return end function isClassicPrcpBlwWmn ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + subroutine classicPrcpBlwWmnScenario(t, vv, pc, ice_pot, severity) IMPLICIT NONE real, intent(in) :: t, pc, vv, ice_pot real, intent(out) :: severity ! deltaZ, ctt: module members (cloud info) ! moistInt: module member real :: tInt, pcInt, dzInt real :: cttInt, vvInt integer :: scenario real :: weights(7) = (/ 3.0, 4.0, 3.0, 3.0, 3.5, 2.5, 2.0 /) real :: sumWeight sumWeight = sum(weights) scenario = SCENARIOS%PRECIPITAION_BELOW_WARMNOSE tInt = t_map(t, scenario) pcInt = prcpCondensate_map(pc, scenario) dzInt = deltaZ_map(deltaZ, scenario) cttInt = ctt_map(ctt) vvInt = vv_map(vv) severity = (weights(1) * tInt + weights(2) * pcInt + & weights(3) * dzInt + weights(4) * cttInt + & weights(5) * moistInt + weights(6) * vvInt + & weights(7) * ice_pot) / sumWeight return end subroutine classicPrcpBlwWmnScenario !-----------------------------------------------------------------------+ ! classical precipitation but not snow logical function isClassicPrcpAbvWmn(prcpType, k) IMPLICIT NONE integer, intent(in) :: prcpType, k ! wmnIdx,topIdx, ctt: module members (cloud info) if((prcpType /= PRECIPS%NONE .and. prcpType /= PRECIPS%SNOW) .and.& (wmnIdx > 0 .and. k > wmnIdx .and. topIdx > wmnIdx) .and. & (ctt < COLD_PRCP_T_THLD)) then isClassicPrcpAbvWmn = .true. else isClassicPrcpAbvWmn = .false. end if return end function isClassicPrcpAbvWmn ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + subroutine classicPrcpAbvWmnScenario(twp, vv, t, pc, ice_pot, severity) IMPLICIT NONE real, intent(in) :: twp, vv, t, pc, ice_pot real, intent(out) :: severity ! moistInt: module member real :: twpInt, vvInt real :: tempAdj, cttAdj, pcAdj integer :: scenario real :: weights(4) = (/ 3.0, 3.5, 4.0, 2.0 /) real :: sumWeight sumWeight = sum(weights) scenario = SCENARIOS%PRECIPITAION_ABOVE_WARMNOSE twpInt = twp_map(twp, scenario) vvInt = vv_map(vv) tempAdj = 0.0 cttAdj = 0.0 pcAdj = 0.0 ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info) call cal_DampingFactors(scenario, t, pc, tempAdj, cttAdj, pcAdj) severity = (weights(1) * twpInt + weights(2) * vvInt + & weights(3) * moistInt + weights(4) * ice_pot) / & (sumWeight + 6.5*tempAdj + 3.0*cttAdj + 3.0*pcAdj) return end subroutine classicPrcpAbvWmnScenario !-----------------------------------------------------------------------+ logical function isNoPrecip(prcpType) IMPLICIT NONE integer, intent(in) :: prcpType isNoPrecip = prcpType == PRECIPS%NONE return end function isNoPrecip ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + subroutine noPrcpScenario(vv, t, pc, ice_pot, severity) IMPLICIT NONE real, intent(in) :: vv, t, pc, ice_pot real, intent(out) :: severity ! layerQ, deltaZ: module members (cloud info) ! moistInt: module member real :: dqInt, dzInt, vvInt real :: tempAdj, cttAdj, pcAdj integer :: scenario real :: weights(5) = (/ 3.0, 3.5, 4.0, 4.0, 3.0 /) real :: sumWeight sumWeight = sum(weights) scenario = SCENARIOS%NO_PRECIPITAION dqInt = deltaQ_map(layerQ) dzInt = deltaZ_map(deltaZ, scenario) vvInt = vv_map(vv) tempAdj = 0.0 cttAdj = 0.0 pcAdj = 0.0 ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info) call cal_DampingFactors(scenario, t, pc, tempAdj, cttAdj, pcAdj) severity = (weights(1) * dqInt + weights(2) * dzInt + & weights(3) * vvInt + weights(4) * moistInt + & weights(5) * ice_pot) / & (sumWeight + 9.0*tempAdj + 4.5*cttAdj) return end subroutine noPrcpScenario !-----------------------------------------------------------------------+ logical function isSnow(prcpType) IMPLICIT NONE integer, intent(in) :: prcpType isSnow = prcpType == PRECIPS%SNOW return end function isSnow ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + subroutine snowScenario(twp, vv, t, pc, ice_pot, severity) IMPLICIT NONE real, intent(in) :: twp, vv, t, pc, ice_pot real, intent(out) :: severity ! deltaZ: module member (cloud info) ! moistInt: module member real :: twpInt, dzInt, vvInt real :: tempAdj, cttAdj, pcAdj integer :: scenario real :: weights(5) = (/ 3.0, 3.5, 3.5, 4.0, 3.0 /) real :: sumWeight sumWeight = sum(weights) scenario = SCENARIOS%ALL_SNOW twpInt = twp_map(twp, scenario) dzInt = deltaZ_map(deltaZ, scenario) vvInt = vv_map(vv) tempAdj = 0.0 cttAdj = 0.0 pcAdj = 0.0 ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info) call cal_DampingFactors(scenario, t, pc, tempAdj, cttAdj, pcAdj) severity = (weights(1) * twpInt + weights(2) * dzInt + & weights(3) * vvInt + weights(4) * moistInt + & weights(5) * ice_pot) / & (sumWeight + 9.0*tempAdj + 4.0*cttAdj + 4.0*pcAdj) return end subroutine snowScenario !-----------------------------------------------------------------------+ logical function isColdRain(prcpType) IMPLICIT NONE integer, intent(in) :: prcpType ! ctt: module members (cloud info) isColdRain = prcpType == PRECIPS%RAIN .and. ctt < COLD_PRCP_T_THLD return end function isColdRain ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + subroutine coldRainScenario(twp, vv, t, pc, ice_pot, severity) IMPLICIT NONE real, intent(in) :: twp, vv, t, pc, ice_pot real, intent(out) :: severity ! deltaZ: module member (cloud info) ! moistInt: module member real :: twpInt, dzInt, vvInt real :: tempAdj, cttAdj, pcAdj integer :: scenario real :: weights(5) = (/ 3.0, 3.5, 3.5, 4.0, 2.0 /) real :: sumWeight sumWeight = sum(weights) scenario = SCENARIOS%COLD_RAIN twpInt = twp_map(twp, scenario) dzInt = deltaZ_map(deltaZ, scenario) vvInt = vv_map(vv) tempAdj = 0.0 cttAdj = 0.0 pcAdj = 0.0 ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info) call cal_DampingFactors(scenario, t, pc, tempAdj, cttAdj, pcAdj) severity = (weights(1) * twpInt + weights(2) * dzInt + & weights(3) * vvInt + weights(4) * moistInt + & weights(5) * ice_pot) / & (sumWeight + 8.0*tempAdj + 4.0*cttAdj + 4.0*pcAdj) return end subroutine coldRainScenario !-----------------------------------------------------------------------+ ! The elseif clause is a result of sloppy thinking on the ! part of the scientitsts. They were trying to create/use ! two different tests for the same 'scenario.' There is ! implicit definition of a 'classic precipitation' scenario. logical function isWarmPrecip(prcpType) IMPLICIT NONE integer, intent(in) :: prcpType ! ctt, wmnIdx, topIdx: module members (cloud info) if((prcpType == PRECIPS%RAIN .or. prcpType == PRECIPS%OTHER) .and. & (ctt >= COLD_PRCP_T_THLD)) then isWarmPrecip = .true. elseif((prcpType /= PRECIPS%NONE .and. prcpType /= PRECIPS%SNOW).and.& (wmnIdx > 0 .and. topIdx <= wmnIdx) .and. & (ctt >= COLD_PRCP_T_THLD)) then isWarmPrecip = .true. else isWarmPrecip = .false. end if return end function isWarmPrecip ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + subroutine warmPrecipScenario(twp, vv, t, pc, ice_pot, severity) IMPLICIT NONE real, intent(in) :: twp, vv, t, pc, ice_pot real, intent(out) :: severity ! deltaZ, layerQ: module member (cloud info) ! moistInt: module member real :: twpInt, dzInt, vvInt, dqInt real :: tempAdj, cttAdj, pcAdj integer :: scenario real :: weights(6) = (/ 3.0, 3.0, 3.0, 3.5, 4.0, 2.0 /) real :: sumWeight sumWeight = sum(weights) scenario = SCENARIOS%WARM_PRECIPITAION twpInt = twp_map(twp, scenario) dzInt = deltaZ_map(deltaZ, scenario) vvInt = vv_map(vv) dqInt = deltaQ_map(layerQ) tempAdj = 0.0 cttAdj = 0.0 pcAdj = 0.0 ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info) call cal_DampingFactors(scenario, t, pc, tempAdj, cttAdj, pcAdj) severity = (weights(1) * twpInt + weights(2) * dzInt + & weights(3) * vvInt + weights(4) * moistInt + & weights(5) * ice_pot + weights(6) * dqInt) / & (sumWeight + 9.0*tempAdj + 4.5*cttAdj + 4.5*pcAdj) return end subroutine warmPrecipScenario !-----------------------------------------------------------------------+ logical function isFreezingPrecip(prcpType) IMPLICIT NONE integer, intent(in) :: prcpType ! ctt: module member (cloud info) if (prcpType == PRECIPS%OTHER .and. ctt < COLD_PRCP_T_THLD) then isFreezingPrecip = .true. else isFreezingPrecip = .false. end if return end function isFreezingPrecip ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + subroutine freezingPrecipScenario(twp, vv, t, pc, ice_pot, severity) IMPLICIT NONE real, intent(in) :: twp, vv, t, pc, ice_pot real, intent(out) :: severity ! deltaZ: module member (cloud info) ! moistInt: module member real :: twpInt, dzInt, vvInt real :: tempAdj, cttAdj, pcAdj integer :: scenario real :: weights(5) = (/ 3.0, 3.5, 3.5, 4.0, 2.0 /) real :: sumWeight sumWeight = sum(weights) scenario = SCENARIOS%FREEZING_PRECIPITAION twpInt = twp_map(twp, scenario) dzInt = deltaZ_map(deltaZ, scenario) vvInt = vv_map(vv) tempAdj = 0.0 cttAdj = 0.0 pcAdj = 0.0 ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info) call cal_DampingFactors(scenario, t, pc, tempAdj, cttAdj, pcAdj) severity = (weights(1) * twpInt + weights(2) * dzInt + & weights(3) * vvInt + weights(4) * moistInt + & weights(5) * ice_pot) / & (sumWeight + 8.0*tempAdj + 4.0*cttAdj + 4.0*pcAdj) return end subroutine freezingPrecipScenario !-----------------------------------------------------------------------+ ! The only scenario NOT applicable to: below the warm nose ! ! Temperature damping - can damp by a max of temp_adjust_factor ! ! CTT damping - can damp by a max of ctt_adjust_factor ! The CTT is multiplied by a map based on the height below cloud top ! For all clouds and it's interest increases the closer to cloud ! top we are. ! ! Condensate damping - can damp by a max of cond_adjust_factor ! Only in the lowest cloud layer and interest increases the closer to ! cloud base we are. Interest will be the max below cloud base as well. ! Maps differ for different scenarios ! subroutine cal_DampingFactors(scenario, t, pc, tempAdj,cttAdj,condAdj) IMPLICIT NONE integer, intent(in) :: scenario real, intent(in) :: t, pc real, intent(out) :: tempAdj, cttAdj, condAdj ! ctt, cloudTopDist, lowestCloud, deltaZ: module member (cloud info) tempAdj = t_map(t, scenario) cttAdj = ctt_map(ctt) * cldTopDist_map(cloudTopDist) if (lowestCloud) then condAdj = cldBaseDist_map(deltaZ) ! For no precipitation, condAdj is 0.0 condAdj = condAdj * prcpCondensate_map(pc, scenario) end if return end subroutine cal_DampingFactors end module IcingSeverity !======================================================================== ! = = = = = = = = = = = = = Icing Algorithm = = = = = = = = = = = = = !======================================================================== subroutine icing_algo(i,j,pres,temp,rh,hgt,omega,wh,& q,cwat,qqw,qqi,qqr,qqs,qqg,& nz,xlat,xlon,xalt,prate,cprate,& cape,cin, ice_pot, ice_sev) use ctlblk_mod, only: imp_physics, SPVAL, DTQ2,me use physcons,only: G => con_g, FV => con_fvirt, RD => con_rd use DerivedFields, only : derive_fields use CloudLayers, only : calc_CloudLayers, clouds_t use IcingPotential, only : icing_pot use IcingSeverity, only : icing_sev IMPLICIT NONE integer, external :: getTopoK !--------------------------------------------------------------------- ! The GFIP algorithm computes the probability of icing within a model ! column given the follow input data ! ! ! 2-D data: total precip rate (prate) ! convective precip rate (cprate) ! the topography height (xalt) ! the latitude and longitude (xlat and xlon) ! the number of vertical levels (nz) = 47 in my GFS file ! ===== for severity only ====== ! Convective Avail. Pot. Energy, pds7=65280 255-0 mb above ! gnd (cape) ! Convective inhibition, pds7=65280 255-0 mb above gnd (cin) ! 3-D data ! pressure(nz) in PA ! temperature(nz) in K ! rh(nz) in % ! hgt(nz) in GPM ! cwat(nz) in kg/x kg ! vv(nz) in Pa/sec !----------------------------------------------------------------------- ! !-----------------------------------------------------------------------+ ! This subroutine calculates the icing potential for a GFS column ! First the topography is mapped to the model's vertical coordinate ! in (topoK) ! Then derive related fields and cloud layers, ! up to 12 layers are possible ! The icing are computed in (icing_pot, ice_sev): ! The icing potential output should range from 0 - 1. ! The icing severity is in 4 categories: 1 2 3 4. ! !-----------------------------------------------------------------------+ integer, intent(in) :: i,j, nz real, intent(in),dimension(nz) :: pres,temp,rh,hgt,omega,wh ! q is used only for converting wh to omega real, intent(in),dimension(nz) :: q,cwat,qqw,qqi,qqr,qqs,qqg real, intent(in) :: xlat, xlon, xalt ! locations real, intent(in) :: prate, cprate ! precipitation rates real, intent(in) :: cape, cin real, intent(out) :: ice_pot(nz), ice_sev(nz) real :: hprcp, hcprcp integer :: topoK, region, prcpType real, allocatable :: vv(:), ept(:), wbt(:), twp(:) real, allocatable :: iceCond(:),liqCond(:), totWater(:),totCond(:) real :: pc, kx, lx, tott type(clouds_t) :: clouds integer :: k real :: tv allocate(vv(nz)) allocate(ept(nz)) allocate(wbt(nz)) allocate(twp(nz)) allocate(iceCond(nz),liqCond(nz)) allocate(totWater(nz),totCond(nz)) allocate(clouds%layerQ(nz)) ! 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 ! do k=1,nz ! 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) ! end do ! if(mod(i,300)==0 .and. mod(j,200)==0)then ! print*,'sample input to FIP ',i,j,nz,xlat,xlon,xalt,prate, cprate ! do k=1,nz ! print*,'k,P,T,RH,H,CWM,VV',k,pres(k),temp(k),rh(k),hgt(k),cwat(k),vv(k) ! end do ! end if topoK = getTopoK(hgt,xalt,nz) ! hourly accumulated precipitation hprcp = prate * 1000. / DTQ2 * 3600. ! large scale total precipitation in 1 hour hcprcp = cprate * 1000. / DTQ2 * 3600. ! large scale convective precipitation in 1 hour if(imp_physics == 98 .or. imp_physics == 99) then do k = 1, nz ! input is omega (pa/s) vv(k) = omega(k) if(cwat(k) == SPVAL) then liqCond(k) = SPVAL iceCond(k) = SPVAL totWater(k) = SPVAL totCond(k) = SPVAL cycle endif if(temp(k) >=259.15) then ! T>=-14C liqCond(k) = cwat(k) * 1000. iceCond(k) = 0. else liqCond(k) = 0. iceCond(k) = cwat(k) * 1000. end if totWater(k) = cwat(k) * 1000. totCond(k) = cwat(k) * 1000. end do else if(imp_physics == 11 .or. imp_physics == 8) then do k = 1, nz ! convert input w (m/s) to omega (pa/s) vv(k) = wh(k) if(vv(k) /= SPVAL) then tv = temp(k)*(1.+q(k)*FV) vv(k)=-vv(k)*G*(pres(k)/(RD*tv)) end if if(qqw(k) /= SPVAL .and. qqr(k) /= SPVAL) then liqCond(k) = (qqw(k) + qqr(k)) * 1000. else liqCond(k) = SPVAL end if if(qqi(k) /= SPVAL .and. qqs(k) /= SPVAL .and. qqg(k) /= SPVAL) then iceCond(k) = (qqi(k) + qqs(k) + qqg(k)) * 1000. else iceCond(k) = SPVAL end if if(liqCond(k) /= SPVAL .and. iceCond(k) /= SPVAL) then totWater(k) = liqCond(k) + iceCond(k) else totWater(k) = SPVAL end if if(qqr(k) /= SPVAL .and. qqs(k) /= SPVAL .and. qqg(k) /= SPVAL) then totCond(k) = (qqr(k) + qqs(k) + qqg(k)) * 1000. else totCond(k) = SPVAL end if end do else print *, "This schema is not supported. imp_physics=", imp_physics return end if call derive_fields(imp_physics,temp, rh, pres, hgt, totWater,totCond, & nz, topoK, hprcp, hcprcp, cin, cape,& ept, wbt, twp, pc, kx, lx, tott, prcpType) call calc_CloudLayers(rh, temp, pres, ept, vv, nz,topoK, xlat,xlon,& region, clouds) ! write(*,'(2x,I3,A,1x,A,2I6,2F8.2)')me,":","i,j,lat,lon=",i,j,xlat,xlon ! do k=1,8 ! write(*,'(2x,I3,A,1x,8F9.2)')me,":",pres((k-1)*8+1:(k-1)*8+8) ! end do call icing_pot(hgt, rh, temp, liqCond, vv, nz, clouds, ice_pot) call icing_sev(imp_physics, hgt, rh, temp, pres, vv, liqCond, iceCond, twp, & ice_pot, nz, hcprcp, cape, lx, kx, tott, pc, prcpType, clouds, & ice_sev) ! if(mod(i,300)==0 .and. mod(j,200)==0)then ! print*,'FIP: cin,cape,pc, kx, lx, tott,pcpTyp',cin,cape,pc, kx, lx, tott,prcpType ! print *, "FIP i,j,lat,lon=",i,j,xlat,xlon ! do k=nz/2,nz ! print*,'k,ept,wbt,twp,totalwater=',k,ept(k), wbt(k), twp(k),& ! totWater(k), ice_pot(k), ice_sev(k) ! print*,'clouds nLayers wmnIdx avv',clouds%nLayers,clouds%wmnIdx,clouds%avv ! if (clouds%nLayers> 0) then ! print*,'clouds topidx=', clouds%topIdx,'baseidx=',clouds%baseIdx,& ! 'ctt=', clouds%ctt ! end if ! end do ! end if deallocate(vv) deallocate(ept) deallocate(wbt) deallocate(twp) deallocate(iceCond,liqCond) deallocate(totWater,totCond) deallocate(clouds%layerQ) return end subroutine icing_algo !-----------------------------------------------------------------------+ integer function getTopoK(hgt, alt, nz) IMPLICIT NONE real, intent(in) :: hgt(nz) real, intent(in) :: alt integer, intent(in) :: nz integer :: k if(hgt(nz) >= alt) then getTopoK = nz else do k=nz,2,-1 if ((hgt(k-1) > alt) .and. (hgt(k) <= alt)) then getTopoK = k-1 exit endif end do end if return end function getTopoK