GFIP3.f

Go to the documentation of this file.

!========================================================================
! = = = = = = = = = = = = module DerivedFields  = = = = = = = = = = = =
!========================================================================
module derivedfields

  implicit none

  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


    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

!-----------------------------------------------------------------------+
  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(*,*)' 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(*,*)'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(*,*)'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

  implicit none

  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)) then  
             if ((rh(kk+1)<t_rh).and.(in_cld==1)) then 
                clouds%baseIdx(num_lyr) = kk
                cur_base = kk
                in_cld = 0
             endif
          endif
       end do
    end if
    kstart = cur_base + 1

    do k = kstart, topok-1 
       if (cur_base<=k) then
          if ((rh(k-1)<t_rh).and.(in_cld==0)) then
             if ((rh(k)>=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)<t_rh)) then  
                      if ((rh(kk+1)<t_rh).and.(in_cld==1)) then 
                         clouds%baseIdx(num_lyr) = kk
                         cur_base = kk
                         in_cld = 0
                      endif
                   endif
                end do

             endif
          endif
       endif
    end do

    ! if the loop exits still in cloud make the cloud base the surface
    if (in_cld==1) then 
       clouds%baseIdx(num_lyr) = topok
    endif

    clouds%nLayers = num_lyr

    clouds%wmnIdx = getwarmnoseidx(t, nz)

    ! only for the lowest cloud layer
    ! only used for severity scenario where precip is below warm nose
    if(num_lyr > 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

  implicit none

  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
  implicit none
  public
  public scenarios

  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)
    implicit none
    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)
    implicit none
    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)
    implicit none
    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)
    implicit none
    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)
    implicit none
    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)
    implicit none
    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)
    implicit none
    real, intent(in) :: 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)
    implicit none
    real, intent(in) :: 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)
    implicit none
    real, intent(in) :: 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)
    implicit none
    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)
    implicit none
    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)
    implicit none
    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)
     implicit none
     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)
     implicit none
     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)
     implicit none
     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)
     implicit none
     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)
     implicit none
     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

  implicit none

  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

    real :: moistint

    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_post,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
!---------------------------------------------------------------------
! icing_algo() 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