gfdl_cloud_microphys_mod.F90

 
!***********************************************************************
!*                   GNU Lesser General Public License
!*
!* This file is part of the GFDL Cloud Microphysics.
!*
!* The GFDL Cloud Microphysics is free software: you can
!* redistribute it and/or modify it under the terms of the
!* GNU Lesser General Public License as published by the
!* Free Software Foundation, either version 3 of the License, or
!* (at your option) any later version.
!*
!* The GFDL Cloud Microphysics is distributed in the hope it will be
!* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
!* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
!* See the GNU General Public License for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with the GFDL Cloud Microphysics.
!* If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************
! =======================================================================
 
module gfdl_cloud_microphys_mod
 
    ! use mpp_mod, only: stdlog, mpp_pe, mpp_root_pe, mpp_clock_id, &
    ! mpp_clock_begin, mpp_clock_end, clock_routine, &
    ! input_nml_file
    ! use diag_manager_mod, only: register_diag_field, send_data
    ! use time_manager_mod, only: time_type, get_time
    ! use constants_mod, only: grav, rdgas, rvgas, cp_air, hlv, hlf, pi => pi_8
    ! use fms_mod, only: write_version_number, open_namelist_file, &
    ! check_nml_error, file_exist, close_file
 
  ! -----------------------------------------------------------------------
   use module_mp_radar
   use module_gfdlmp_param, only: read_gfdlmp_nml, mp_time, t_min, t_sub, &
        tau_r2g, tau_smlt, tau_g2r, dw_land, dw_ocean, vi_fac, vr_fac,    &
        vs_fac, vg_fac, ql_mlt, do_qa, fix_negative, vi_max, vs_max,      &
        vg_max, vr_max, qs_mlt, qs0_crt, qi_gen, ql0_max, qi0_max,        &
        qi0_crt, qr0_crt, fast_sat_adj, rh_inc, rh_ins, rh_inr, const_vi, &
        const_vs, const_vg, const_vr, use_ccn, rthresh, ccn_l, ccn_o,     &
        qc_crt, tau_g2v, tau_v2g, sat_adj0, c_piacr, tau_imlt, tau_v2l,   &
        tau_l2v, tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci,       &
        c_pgacs, z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, &
        tice, rad_snow, rad_graupel, rad_rain, cld_min, use_ppm,          &
        mono_prof, do_sedi_heat, sedi_transport, do_sedi_w, de_ice,       &
        icloud_f, irain_f, mp_print, reiflag, rewmin, rewmax, reimin,     &
        reimax, rermin, rermax, resmin, resmax, regmin, regmax, tintqs,   &
        do_hail
 
   implicit none
 
   private
 
   public gfdl_cloud_microphys_mod_driver, gfdl_cloud_microphys_mod_init, &
          gfdl_cloud_microphys_mod_end, cloud_diagnosis
!   public wqs1, wqs2, qs_blend, wqsat_moist, wqsat2_moist
!   public qsmith_init, qsmith, es2_table1d, es3_table1d, esw_table1d
!   public setup_con, wet_bulb
 
   real :: missing_value = - 1.e10
 
   logical :: module_is_initialized = .false.
   logical :: qsmith_tables_initialized = .false.
 
   character (len = 20) :: mod_name = 'gfdl_cloud_microphys'
 
   real, parameter :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6
   real, parameter :: rhos = 0.1e3, rhog = 0.4e3
   real,                 parameter :: grav = 9.80665
   real,                 parameter :: rdgas = 287.05
   real,                 parameter :: rvgas = 461.50
   real,                 parameter :: cp_air = 1004.6
   real,                 parameter :: hlv = 2.5e6          
   real,                 parameter :: hlf = 3.3358e5       
   real,                 parameter :: pi = 3.1415926535897931
 
   ! real, parameter :: rdgas = 287.04                     !< gfdl: gas constant for dry air
 
   ! real, parameter :: cp_air = rdgas * 7. / 2.           ! 1004.675, heat capacity of dry air at constant pressure
   real, parameter :: cp_vap = 4.0 * rvgas                 
   ! real, parameter :: cv_air = 717.56                    ! satoh value
   real, parameter :: cv_air = cp_air - rdgas              
   ! real, parameter :: cv_vap = 1410.0                    ! emanuel value
   real, parameter :: cv_vap = 3.0 * rvgas                 
 
   ! the following two are from emanuel's book "atmospheric convection"
   ! real, parameter :: c_ice = 2106.0                     ! heat capacity of ice at 0 deg c: c = c_ice + 7.3 * (t - tice)
   ! real, parameter :: c_liq = 4190.0                     ! heat capacity of water at 0 deg c
 
   real, parameter :: c_ice = 1972.0
   real, parameter :: c_liq = 4185.5
   ! real, parameter :: c_liq = 4218.0                     ! ifs: heat capacity of liquid at 0 deg c
 
   real, parameter :: eps = rdgas / rvgas                  
   real, parameter :: zvir = rvgas / rdgas - 1.            
 
   real, parameter :: t_ice = 273.16
   real, parameter :: table_ice = 273.16
 
   ! real, parameter :: e00 = 610.71                       ! gfdl: saturation vapor pressure at 0 deg c
   real, parameter :: e00 = 611.21
 
   real, parameter :: dc_vap = cp_vap - c_liq              
   real, parameter :: dc_ice = c_liq - c_ice               
 
   real, parameter :: hlv0 = hlv                           
   ! real, parameter :: hlv0 = 2.501e6                     ! emanuel appendix - 2
   real, parameter :: hlf0 = hlf                           
   ! real, parameter :: hlf0 = 3.337e5                     ! emanuel
 
   real, parameter :: lv0 = hlv0 - dc_vap * t_ice          
   real, parameter :: li00 = hlf0 - dc_ice * t_ice         
 
   real, parameter :: d2ice = dc_vap + dc_ice              
   real, parameter :: li2 = lv0 + li00                     
 
   real, parameter :: qrmin = 1.e-8                        
   real, parameter :: qvmin = 1.e-20                       
   real, parameter :: qcmin = 1.e-12                       
 
   real, parameter :: vr_min = 1.e-3                       
   real, parameter :: vf_min = 1.e-5                       
 
   real, parameter :: dz_min = 1.e-2                       
 
   real, parameter :: sfcrho = 1.2
   real, parameter :: rhor = 1.e3                          
    ! intercept parameters
 
   real, parameter :: rnzr = 8.0e6 ! lin83
   real, parameter :: rnzs = 3.0e6 ! lin83
   real, parameter :: rnzg = 4.0e6 ! rh84
   real, parameter :: rnzh = 4.0e4 ! lin83 --- lmh 29 sep 17
 
    ! density parameters
 
   real, parameter :: rhoh = 0.917e3 ! lin83 --- lmh 29 sep 17
 
   public rhor, rhos, rhog, rhoh, rnzr, rnzs, rnzg, rnzh
   real :: cracs, csacr, cgacr, cgacs, csacw, craci, csaci, cgacw, cgaci, cracw
   real :: acco (3, 4)
   real :: cssub (5), cgsub (5), crevp (5), cgfr (2), csmlt (5), cgmlt (5)
 
   real :: es0, ces0
   real :: pie, rgrav, fac_rc
   real :: c_air, c_vap
 
   real :: lati, latv, lats, lat2, lcp, icp, tcp
 
   real :: d0_vap
   real :: lv00
 
   ! cloud microphysics switchers
   logical :: do_setup = .true.                            
   logical :: p_nonhydro = .false.                         
 
   real, allocatable :: table (:), table2 (:), table3 (:), tablew (:)
   real, allocatable :: des (:), des2 (:), des3 (:), desw (:)
 
    logical :: tables_are_initialized = .false.
 
   ! logical :: master
   ! integer :: id_rh, id_vtr, id_vts, id_vtg, id_vti, id_rain, id_snow, id_graupel, &
   ! id_ice, id_prec, id_cond, id_var, id_droplets
   real, parameter :: dt_fr = 8.
   ! minimum temperature water can exist (moore & molinero nov. 2011, nature)
   ! dt_fr can be considered as the error bar
 
   real :: p_min = 100.
 
   ! slj, the following parameters are for cloud - resolving resolution: 1 - 5 km
 
   ! qi0_crt = 0.8e-4
   ! qs0_crt = 0.6e-3
   ! c_psaci = 0.1
   ! c_pgacs = 0.1
   
   real :: log_10, tice0, t_wfr
 
contains
 
! -----------------------------------------------------------------------
! the driver of the gfdl cloud microphysics
! -----------------------------------------------------------------------
 
subroutine gfdl_cloud_microphys_mod_driver (                                    &
            iis, iie, jjs, jje, kks, kke, ktop, kbot,                           &
            qv, ql, qr, qi, qs, qg, qa, qn,                                     &
            qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, pt_dt, pt, w,      &
            uin, vin, udt, vdt, dz, delp, area, dt_in, land,                    &
            rain, snow, ice, graupel, hydrostatic, phys_hydrostatic,            &
            p, lradar, refl_10cm, reset, pfils, pflls)
 
   implicit none
 
   ! Interface variables
   integer, intent (in) :: iis, iie, jjs, jje ! physics window
   integer, intent (in) :: kks, kke ! vertical dimension
   integer, intent (in) :: ktop, kbot ! vertical compute domain
 
   real, intent (in) :: dt_in ! physics time step
 
   real, intent (in), dimension (iis:iie, jjs:jje) :: area ! cell area
   real, intent (in), dimension (iis:iie, jjs:jje) :: land ! land fraction
 
   real, intent (in), dimension (iis:iie, jjs:jje, kks:kke) :: delp, dz, uin, vin
   real, intent (in), dimension (iis:iie, jjs:jje, kks:kke) :: pt, qv, ql, qr, qg, qa, qn
 
   real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: qi, qs
   real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: pt_dt, qa_dt, udt, vdt, w
   real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: qv_dt, ql_dt, qr_dt
   real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: qi_dt, qs_dt, qg_dt
 
   real, intent (out), dimension (iis:iie, jjs:jje) :: rain, snow, ice, graupel
 
   logical, intent (in) :: hydrostatic, phys_hydrostatic
 
   !integer, intent (in) :: seconds
   real, intent (in), dimension (iis:iie, jjs:jje, kks:kke) :: p
   logical, intent (in) :: lradar
   real, intent (out), dimension (iis:iie, jjs:jje, kks:kke) :: refl_10cm
   logical, intent (in) :: reset
   real, intent (out), dimension (iis:iie, jjs:jje, kks:kke) :: pfils, pflls
 
   ! Local variables
   logical :: melti = .false.
 
   real :: mpdt, rdt, dts, convt, tot_prec
 
   integer :: i, j, k
   integer :: is, ie, js, je ! physics window
   integer :: ks, ke ! vertical dimension
   integer :: days, ntimes, kflip
 
   real, dimension (iie-iis+1, jje-jjs+1) :: prec_mp, prec1, cond, w_var, rh0
 
   real, dimension (iie-iis+1, jje-jjs+1,kke-kks+1) :: vt_r, vt_s, vt_g, vt_i, qn2
 
   real, dimension (size(pt,1), size(pt,3)) :: m2_rain, m2_sol
 
   real :: allmax
!+---+-----------------------------------------------------------------+
!For 3D reflectivity calculations
   real, dimension(ktop:kbot):: qv1d, t1d, p1d, qr1d, qs1d, qg1d, dbz
!+---+-----------------------------------------------------------------+
 
   is = 1
   js = 1
   ks = 1
   ie = iie - iis + 1
   je = jje - jjs + 1
   ke = kke - kks + 1
   ! call mpp_clock_begin (gfdl_mp_clock)
 
   ! -----------------------------------------------------------------------
   ! define heat capacity of dry air and water vapor based on hydrostatical property
   ! -----------------------------------------------------------------------
 
   if (phys_hydrostatic .or. hydrostatic) then
       c_air = cp_air
       c_vap = cp_vap
       p_nonhydro = .false.
   else
       c_air = cv_air
       c_vap = cv_vap
       p_nonhydro = .true.
   endif
   d0_vap = c_vap - c_liq
   lv00 = hlv0 - d0_vap * t_ice
 
   if (hydrostatic) do_sedi_w = .false.
 
   ! -----------------------------------------------------------------------
   ! define latent heat coefficient used in wet bulb and Bigg mechanism
   ! -----------------------------------------------------------------------
 
   latv = hlv
   lati = hlf
   lats = latv + lati
   lat2 = lats * lats
 
   lcp = latv / cp_air
   icp = lati / cp_air
   tcp = (latv + lati) / cp_air
 
   ! tendency zero out for am moist processes should be done outside the driver
 
   ! -----------------------------------------------------------------------
   ! define cloud microphysics sub time step
   ! -----------------------------------------------------------------------
 
   mpdt   = min(dt_in, mp_time)
   rdt    = 1. / dt_in
   ntimes = nint(dt_in / mpdt)
 
   ! small time step:
   dts = dt_in / real(ntimes)
 
   ! call get_time (time, seconds, days)
 
   ! -----------------------------------------------------------------------
   ! initialize precipitation
   ! -----------------------------------------------------------------------
 
   do j = js, je
       do i = is, ie
           graupel(i, j) = 0.
           rain(i, j) = 0.
           snow(i, j) = 0.
           ice(i, j) = 0.
           cond(i, j) = 0.
       enddo
   enddo
 
   pfils = 0.
   pflls = 0.
 
   ! -----------------------------------------------------------------------
   ! major cloud microphysics
   ! -----------------------------------------------------------------------
 
   do j = js, je
       call mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, qg,&
           qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes,  &
           rain(:, j), snow(:, j), graupel(:, j), ice(:, j), m2_rain,     &
           m2_sol, cond(:, j), area(:, j), land(:, j), udt, vdt, pt_dt,    &
           qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, w_var, vt_r,      &
           vt_s, vt_g, vt_i, qn2)
       do k = ktop, kbot
           do i = is, ie
               pfils(i, j, k) = m2_sol(i, k)
               pflls(i, j, k) = m2_rain(i, k)
           enddo
       enddo
   enddo
 
   ! -----------------------------------------------------------------------
   ! no clouds allowed above ktop
   ! -----------------------------------------------------------------------
 
   if (ks < ktop) then
       do k = ks, ktop
           if (do_qa) then
               do j = js, je
                   do i = is, ie
                       qa_dt(i, j, k) = 0.
                   enddo
               enddo
           else
               do j = js, je
                   do i = is, ie
                       ! qa_dt (i, j, k) = - qa (i, j, k) * rdt
                       qa_dt(i, j, k) = 0. ! gfs
                   enddo
               enddo
           endif
       enddo
   endif
 
   ! -----------------------------------------------------------------------
   ! diagnostic output
   ! -----------------------------------------------------------------------
 
   ! if (id_vtr > 0) then
   ! used = send_data (id_vtr, vt_r, time, is_in = iis, js_in = jjs)
   ! endif
 
   ! if (id_vts > 0) then
   ! used = send_data (id_vts, vt_s, time, is_in = iis, js_in = jjs)
   ! endif
 
   ! if (id_vtg > 0) then
   ! used = send_data (id_vtg, vt_g, time, is_in = iis, js_in = jjs)
   ! endif
 
   ! if (id_vti > 0) then
   ! used = send_data (id_vti, vt_i, time, is_in = iis, js_in = jjs)
   ! endif
 
   ! if (id_droplets > 0) then
   ! used = send_data (id_droplets, qn2, time, is_in = iis, js_in = jjs)
   ! endif
 
   ! if (id_var > 0) then
   ! used = send_data (id_var, w_var, time, is_in = iis, js_in = jjs)
   ! endif
 
   ! convert to mm / day
 
   convt = 86400. * rdt * rgrav
   do j = js, je
       do i = is, ie
           rain(i, j) = rain(i, j) * convt
           snow(i, j) = snow(i, j) * convt
           ice(i, j) = ice(i, j) * convt
           graupel(i, j) = graupel(i, j) * convt
           prec_mp(i, j) = rain(i, j) + snow(i, j) + ice(i, j) + graupel(i, j)
       enddo
   enddo
 
   ! if (id_cond > 0) then
   ! do j = js, je
   ! do i = is, ie
   ! cond (i, j) = cond (i, j) * rgrav
   ! enddo
   ! enddo
   ! used = send_data (id_cond, cond, time, is_in = iis, js_in = jjs)
   ! endif
 
   ! if (id_snow > 0) then
   ! used = send_data (id_snow, snow, time, iis, jjs)
   ! used = send_data (id_snow, snow, time, is_in = iis, js_in = jjs)
   ! if (mp_print .and. seconds == 0) then
   ! tot_prec = g_sum (snow, is, ie, js, je, area, 1)
   ! if (master) write (*, *) 'mean snow = ', tot_prec
   ! endif
   ! endif
   !
   ! if (id_graupel > 0) then
   ! used = send_data (id_graupel, graupel, time, iis, jjs)
   ! used = send_data (id_graupel, graupel, time, is_in = iis, js_in = jjs)
   ! if (mp_print .and. seconds == 0) then
   ! tot_prec = g_sum (graupel, is, ie, js, je, area, 1)
   ! if (master) write (*, *) 'mean graupel = ', tot_prec
   ! endif
   ! endif
   !
   ! if (id_ice > 0) then
   ! used = send_data (id_ice, ice, time, iis, jjs)
   ! used = send_data (id_ice, ice, time, is_in = iis, js_in = jjs)
   ! if (mp_print .and. seconds == 0) then
   ! tot_prec = g_sum (ice, is, ie, js, je, area, 1)
   ! if (master) write (*, *) 'mean ice_mp = ', tot_prec
   ! endif
   ! endif
   !
   ! if (id_rain > 0) then
   ! used = send_data (id_rain, rain, time, iis, jjs)
   ! used = send_data (id_rain, rain, time, is_in = iis, js_in = jjs)
   ! if (mp_print .and. seconds == 0) then
   ! tot_prec = g_sum (rain, is, ie, js, je, area, 1)
   ! if (master) write (*, *) 'mean rain = ', tot_prec
   ! endif
   ! endif
   !
   ! if (id_rh > 0) then !not used?
   ! used = send_data (id_rh, rh0, time, iis, jjs)
   ! used = send_data (id_rh, rh0, time, is_in = iis, js_in = jjs)
   ! endif
   !
   !
   ! if (id_prec > 0) then
   ! used = send_data (id_prec, prec_mp, time, iis, jjs)
   ! used = send_data (id_prec, prec_mp, time, is_in = iis, js_in = jjs)
   ! endif
 
   ! if (mp_print) then
   ! prec1 (:, :) = prec1 (:, :) + prec_mp (:, :)
   ! if (seconds == 0) then
   ! prec1 (:, :) = prec1 (:, :) * dt_in / 86400.
   ! tot_prec = g_sum (prec1, is, ie, js, je, area, 1)
   ! if (master) write (*, *) 'daily prec_mp = ', tot_prec
   ! prec1 (:, :) = 0.
   ! endif
   ! endif
 
   ! call mpp_clock_end (gfdl_mp_clock)
    if(lradar) then
       ! Only set melti to true at the output times
       if (reset) then
         melti=.true.
       else
         melti=.false.
       endif
       do j = js, je
          do i = is, ie
             do k = ktop,kbot
               kflip = kbot-ktop+1-k+1
               t1d(k)  = pt(i,j,kflip)
               p1d(k)  = p(i,j,kflip)
               qv1d(k) = qv(i,j,kflip)/(1-qv(i,j,kflip))
               qr1d(k) = qr(i,j,kflip)
               qs1d(k) = qs(i,j,kflip)
               qg1d(k) = qg(i,j,kflip)
             enddo
             call refl10cm_gfdl (qv1d, qr1d, qs1d, qg1d,                 &
                       t1d, p1d, dbz, ktop, kbot, i, j, melti)
             do k = ktop,kbot
                kflip = kbot-ktop+1-k+1
                refl_10cm(i,j,kflip) = max(-35., dbz(k))
             enddo
          enddo
       enddo
    endif
 
end subroutine gfdl_cloud_microphys_mod_driver
 
! -----------------------------------------------------------------------
subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs,     &
        qg, qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes, &
        rain, snow, graupel, ice, m2_rain, m2_sol, cond, area1, land,         &
        u_dt, v_dt, pt_dt, qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt,   &
        w_var, vt_r, vt_s, vt_g, vt_i, qn2)
 
    implicit none
 
    logical, intent (in) :: hydrostatic
 
    integer, intent (in) :: j, is, ie, js, je, ks, ke
    integer, intent (in) :: ntimes, ktop, kbot
 
    real, intent (in) :: dt_in
 
    real, intent (in), dimension (is:) :: area1, land
 
    real, intent (in), dimension (is:, js:, ks:) :: uin, vin, delp, pt, dz
    real, intent (in), dimension (is:, js:, ks:) :: qv, ql, qr, qg, qa, qn
 
    real, intent (inout), dimension (is:, js:, ks:) :: qi, qs
    real, intent (inout), dimension (is:, js:, ks:) :: u_dt, v_dt, w, pt_dt, qa_dt
    real, intent (inout), dimension (is:, js:, ks:) :: qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt
 
    real, intent (inout), dimension (is:) :: rain, snow, ice, graupel, cond
 
    real, intent (out), dimension (is:, js:) :: w_var
 
    real, intent (out), dimension (is:, js:, ks:) :: vt_r, vt_s, vt_g, vt_i, qn2
 
    real, intent (out), dimension (is:, ks:) :: m2_rain, m2_sol
 
    real, dimension (ktop:kbot) :: qvz, qlz, qrz, qiz, qsz, qgz, qaz
    real, dimension (ktop:kbot) :: vtiz, vtsz, vtgz, vtrz
    real, dimension (ktop:kbot) :: dp0, dp1, dz0, dz1
    real, dimension (ktop:kbot) :: qv0, ql0, qr0, qi0, qs0, qg0, qa0
    real, dimension (ktop:kbot) :: t0, den, den0, tz, p1, denfac
    real, dimension (ktop:kbot) :: ccn, c_praut, m1_rain, m1_sol, m1
    real, dimension (ktop:kbot) :: u0, v0, u1, v1, w1
 
    real :: cpaut, rh_adj, rh_rain
    real :: r1, s1, i1, g1, rdt, ccn0
    real :: dt_rain, dts
    real :: s_leng, t_land, t_ocean, h_var
    real :: cvm, tmp, omq
    real :: dqi, qio, qin
 
    integer :: i, k, n
 
    dts = dt_in / real(ntimes)
    dt_rain = dts * 0.5
    rdt = 1. / dt_in
 
    ! -----------------------------------------------------------------------
    ! use local variables
    ! -----------------------------------------------------------------------
    
    !GJF: assign values to intent(out) variables that are commented out
    w_var = 0.0
    vt_r = 0.0
    vt_s = 0.0
    vt_g = 0.0
    vt_i = 0.0
    qn2 = 0.0
    !GJF
    
    do i = is, ie
 
        do k = ktop, kbot
            qiz(k) = qi(i, j, k)
            qsz(k) = qs(i, j, k)
        enddo
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        if (de_ice) then
            do k = ktop, kbot
                qio = qiz(k) - dt_in * qi_dt(i, j, k) ! original qi before phys
                qin = max(qio, qi0_max) ! adjusted value
                if (qiz(k) > qin) then
                    qsz(k) = qsz(k) + qiz(k) - qin
                    qiz(k) = qin
                    dqi = (qin - qio) * rdt ! modified qi tendency
                    qs_dt(i, j, k) = qs_dt(i, j, k) + qi_dt(i, j, k) - dqi
                    qi_dt(i, j, k) = dqi
                    qi(i, j, k) = qiz(k)
                    qs(i, j, k) = qsz(k)
                endif
            enddo
        endif
 
        do k = ktop, kbot
 
            t0(k) = pt(i, j, k)
            tz(k) = t0(k)
            dp1(k) = delp(i, j, k)
            dp0(k) = dp1(k) ! moist air mass * grav
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            qvz(k) = qv(i, j, k)
            qlz(k) = ql(i, j, k)
            qrz(k) = qr(i, j, k)
            qgz(k) = qg(i, j, k)
 
            ! dp1: dry air_mass
            ! dp1 (k) = dp1 (k) * (1. - (qvz (k) + qlz (k) + qrz (k) + qiz (k) + qsz (k) + qgz (k)))
            dp1(k) = dp1(k) * (1. - qvz(k)) ! gfs
            omq = dp0(k) / dp1(k)
 
            qvz(k) = qvz(k) * omq
            qlz(k) = qlz(k) * omq
            qrz(k) = qrz(k) * omq
            qiz(k) = qiz(k) * omq
            qsz(k) = qsz(k) * omq
            qgz(k) = qgz(k) * omq
 
            qa0(k) = qa(i, j, k)
            qaz(k) = 0.
            dz0(k) = dz(i, j, k)
 
            den0(k) = - dp1(k) / (grav * dz0(k)) ! density of dry air
            p1(k) = den0(k) * rdgas * t0(k) ! dry air pressure
 
            ! -----------------------------------------------------------------------
            ! save a copy of old value for computing tendencies
            ! -----------------------------------------------------------------------
 
            qv0(k) = qvz(k)
            ql0(k) = qlz(k)
            qr0(k) = qrz(k)
            qi0(k) = qiz(k)
            qs0(k) = qsz(k)
            qg0(k) = qgz(k)
 
            ! -----------------------------------------------------------------------
            ! for sedi_momentum
            ! -----------------------------------------------------------------------
 
            m1(k) = 0.
            u0(k) = uin(i, j, k)
            v0(k) = vin(i, j, k)
            u1(k) = u0(k)
            v1(k) = v0(k)
 
        enddo
 
        if (do_sedi_w) then
            do k = ktop, kbot
                w1(k) = w(i, j, k)
            enddo
        ! Set to -999 if not used
        else
            w1(:) = -999.0
        endif
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        cpaut = c_paut * 0.104 * grav / 1.717e-5
 
        if (prog_ccn) then
            do k = ktop, kbot
                ! convert # / cc to # / m^3
                ccn(k) = qn(i, j, k) * 1.e6
                c_praut(k) = cpaut * (ccn(k) * rhor) ** (- 1. / 3.)
            enddo
            use_ccn = .false.
        else
            ccn0 = (ccn_l * land(i) + ccn_o * (1. - land(i))) * 1.e6
            if (use_ccn) then
                ! -----------------------------------------------------------------------
                ! ccn is formulted as ccn = ccn_surface * (den / den_surface)
                ! -----------------------------------------------------------------------
                ccn0 = ccn0 * rdgas * tz(kbot) / p1(kbot)
            endif
            tmp = cpaut * (ccn0 * rhor) ** (- 1. / 3.)
            do k = ktop, kbot
                c_praut(k) = tmp
                ccn(k) = ccn0
            enddo
        endif
 
        ! -----------------------------------------------------------------------
        ! total water subgrid deviation in horizontal direction
        ! default area dependent form: use dx ~ 100 km as the base
        ! -----------------------------------------------------------------------
 
        s_leng = sqrt(sqrt(area1(i) / 1.e10))
        t_land = dw_land * s_leng
        t_ocean = dw_ocean * s_leng
        h_var = t_land * land(i) + t_ocean * (1. - land(i))
        h_var = min(0.20, max(0.01, h_var))
        ! if (id_var > 0) w_var (i, j) = h_var
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        rh_adj = 1. - h_var - rh_inc
        rh_rain = max(0.35, rh_adj - rh_inr) ! rh_inr = 0.25
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        if (fix_negative) &
            call neg_adj (ktop, kbot, tz, dp1, qvz, qlz, qrz, qiz, qsz, qgz)
 
        m2_rain(i, :) = 0.
        m2_sol(i, :) = 0.
 
        do n = 1, ntimes
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            if (p_nonhydro) then
                do k = ktop, kbot
                    dz1(k) = dz0(k)
                    den(k) = den0(k) ! dry air density remains the same
                    denfac(k) = sqrt(sfcrho / den(k))
                enddo
            else
                do k = ktop, kbot
                    dz1(k) = dz0(k) * tz(k) / t0(k) ! hydrostatic balance
                    den(k) = den0(k) * dz0(k) / dz1(k)
                    denfac(k) = sqrt(sfcrho / den(k))
                enddo
            endif
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            call warm_rain (dt_rain, ktop, kbot, dp1, dz1, tz, qvz, qlz, qrz, qiz, qsz, &
                qgz, den, denfac, ccn, c_praut, rh_rain, vtrz, r1, m1_rain, w1, h_var)
 
            rain(i) = rain(i) + r1
 
            do k = ktop, kbot
                m2_rain(i, k) = m2_rain(i, k) + m1_rain(k)
                m1(k) = m1(k) + m1_rain(k)
            enddo
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            call fall_speed (ktop, kbot, den, qsz, qiz, qgz, qlz, tz, vtsz, vtiz, vtgz)
 
            call terminal_fall (dts, ktop, kbot, tz, qvz, qlz, qrz, qgz, qsz, qiz, &
                dz1, dp1, den, vtgz, vtsz, vtiz, r1, g1, s1, i1, m1_sol, w1)
 
            rain(i) = rain(i) + r1 ! from melted snow & ice that reached the ground
            snow(i) = snow(i) + s1
            graupel(i) = graupel(i) + g1
            ice(i) = ice(i) + i1
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            if (do_sedi_heat) &
                call sedi_heat (ktop, kbot, dp1, m1_sol, dz1, tz, qvz, qlz, qrz, qiz, &
                    qsz, qgz, c_ice)
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            call warm_rain (dt_rain, ktop, kbot, dp1, dz1, tz, qvz, qlz, qrz, qiz, qsz, &
                qgz, den, denfac, ccn, c_praut, rh_rain, vtrz, r1, m1_rain, w1, h_var)
 
            rain(i) = rain(i) + r1
 
            do k = ktop, kbot
                m2_rain(i, k) = m2_rain(i, k) + m1_rain(k)
                m2_sol(i, k) = m2_sol(i, k) + m1_sol(k)
                m1(k) = m1(k) + m1_rain(k) + m1_sol(k)
            enddo
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            call icloud (ktop, kbot, tz, p1, qvz, qlz, qrz, qiz, qsz, qgz, dp1, den, &
                denfac, vtsz, vtgz, vtrz, qaz, rh_adj, rh_rain, dts, h_var)
 
        enddo
 
        ! convert units from Pa*kg/kg to kg/m^2/s
        m2_rain(i, :) = m2_rain(i, :) * rdt * rgrav
        m2_sol(i, :) = m2_sol(i, :) * rdt * rgrav
 
        ! -----------------------------------------------------------------------
        ! note: dp1 is dry mass; dp0 is the old moist (total) mass
        ! -----------------------------------------------------------------------
 
        if (sedi_transport) then
            do k = ktop + 1, kbot
                u1(k) = (dp0(k) * u1(k) + m1(k - 1) * u1(k - 1)) / (dp0(k) + m1(k - 1))
                v1(k) = (dp0(k) * v1(k) + m1(k - 1) * v1(k - 1)) / (dp0(k) + m1(k - 1))
                u_dt(i, j, k) = u_dt(i, j, k) + (u1(k) - u0(k)) * rdt
                v_dt(i, j, k) = v_dt(i, j, k) + (v1(k) - v0(k)) * rdt
            enddo
        endif
 
        if (do_sedi_w) then
            do k = ktop, kbot
                w(i, j, k) = w1(k)
            enddo
        endif
 
        ! -----------------------------------------------------------------------
        ! convert to dry mixing ratios
        ! -----------------------------------------------------------------------
 
        do k = ktop, kbot
            omq = dp1(k) / dp0(k)
            qv_dt(i, j, k) = qv_dt(i, j, k) + rdt * (qvz(k) - qv0(k)) * omq
            ql_dt(i, j, k) = ql_dt(i, j, k) + rdt * (qlz(k) - ql0(k)) * omq
            qr_dt(i, j, k) = qr_dt(i, j, k) + rdt * (qrz(k) - qr0(k)) * omq
            qi_dt(i, j, k) = qi_dt(i, j, k) + rdt * (qiz(k) - qi0(k)) * omq
            qs_dt(i, j, k) = qs_dt(i, j, k) + rdt * (qsz(k) - qs0(k)) * omq
            qg_dt(i, j, k) = qg_dt(i, j, k) + rdt * (qgz(k) - qg0(k)) * omq
            cvm = c_air + qvz(k) * c_vap + (qrz(k) + qlz(k)) * c_liq + (qiz(k) + qsz(k) + qgz(k)) * c_ice
            pt_dt(i, j, k) = pt_dt(i, j, k) + rdt * (tz(k) - t0(k)) * cvm / cp_air
        enddo
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        do k = ktop, kbot
            if (do_qa) then
                qa_dt(i, j, k) = 0.
            else
                qa_dt(i, j, k) = qa_dt(i, j, k) + rdt * (qaz(k) / real(ntimes) - qa0(k))
            endif
        enddo
 
        ! -----------------------------------------------------------------------
        ! fms diagnostics:
        ! -----------------------------------------------------------------------
 
        ! if (id_cond > 0) then
        ! do k = ktop, kbot ! total condensate
        ! cond (i) = cond (i) + dp1 (k) * (qlz (k) + qrz (k) + qsz (k) + qiz (k) + qgz (k))
        ! enddo
        ! endif
        !
        ! if (id_vtr > 0) then
        ! do k = ktop, kbot
        ! vt_r (i, j, k) = vtrz (k)
        ! enddo
        ! endif
        !
        ! if (id_vts > 0) then
        ! do k = ktop, kbot
        ! vt_s (i, j, k) = vtsz (k)
        ! enddo
        ! endif
        !
        ! if (id_vtg > 0) then
        ! do k = ktop, kbot
        ! vt_g (i, j, k) = vtgz (k)
        ! enddo
        ! endif
        !
        ! if (id_vts > 0) then
        ! do k = ktop, kbot
        ! vt_i (i, j, k) = vtiz (k)
        ! enddo
        ! endif
        !
        ! if (id_droplets > 0) then
        ! do k = ktop, kbot
        ! qn2 (i, j, k) = ccn (k)
        ! enddo
        ! endif
 
    enddo
 
end subroutine mpdrv
 
! -----------------------------------------------------------------------
subroutine sedi_heat (ktop, kbot, dm, m1, dz, tz, qv, ql, qr, qi, qs, qg, cw)
 
    implicit none
 
    ! input q fields are dry mixing ratios, and dm is dry air mass
 
    integer, intent (in) :: ktop, kbot
 
    real, intent (in), dimension (ktop:kbot) :: dm, m1, dz, qv, ql, qr, qi, qs, qg
 
    real, intent (inout), dimension (ktop:kbot) :: tz
 
    real, intent (in) :: cw ! heat capacity
 
    real, dimension (ktop:kbot) :: dgz, cvn
 
    real :: tmp
 
    integer :: k
 
    do k = ktop, kbot
        dgz(k) = - 0.5 * grav * dz(k) ! > 0
        cvn(k) = dm(k) * (cv_air + qv(k) * cv_vap + (qr(k) + ql(k)) * &
            c_liq + (qi(k) + qs(k) + qg(k)) * c_ice)
    enddo
 
    ! -----------------------------------------------------------------------
    ! sjl, july 2014
    ! assumption: the ke in the falling condensates is negligible compared to the potential energy
    ! that was unaccounted for. local thermal equilibrium is assumed, and the loss in pe is transformed
    ! into internal energy (to heat the whole grid box)
    ! backward time - implicit upwind transport scheme:
    ! dm here is dry air mass
    ! -----------------------------------------------------------------------
 
    k = ktop
    tmp = cvn(k) + m1(k) * cw
    tz(k) = (tmp * tz(k) + m1(k) * dgz(k)) / tmp
 
    ! -----------------------------------------------------------------------
    ! implicit algorithm: can't be vectorized
    ! needs an inner i - loop for vectorization
    ! -----------------------------------------------------------------------
 
    do k = ktop + 1, kbot
        tz(k) = ((cvn(k) + cw * (m1(k) - m1(k - 1))) * tz(k) + m1(k - 1) * &
            cw * tz(k - 1) + dgz(k) * (m1(k - 1) + m1(k))) / (cvn(k) + cw * m1(k))
    enddo
 
end subroutine sedi_heat
 
! -----------------------------------------------------------------------
subroutine warm_rain (dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, &
        den, denfac, ccn, c_praut, rh_rain, vtr, r1, m1_rain, w1, h_var)
 
    implicit none
 
    integer, intent (in) :: ktop, kbot
 
    real, intent (in) :: dt ! time step (s)
    real, intent (in) :: rh_rain, h_var
 
    real, intent (in), dimension (ktop:kbot) :: dp, dz, den
    real, intent (in), dimension (ktop:kbot) :: denfac, ccn, c_praut
 
    real, intent (inout), dimension (ktop:kbot) :: tz, vtr
    real, intent (inout), dimension (ktop:kbot) :: qv, ql, qr, qi, qs, qg
    real, intent (inout), dimension (ktop:kbot) :: m1_rain, w1
 
    real, intent (out) :: r1
 
    real, parameter :: so3 = 7. / 3.
 
    real, dimension (ktop:kbot) :: dl, dm
    real, dimension (ktop:kbot + 1) :: ze, zt
 
    real :: sink, dq, qc0, qc
    real :: qden
    real :: zs = 0.
    real :: dt5
 
    integer :: k
 
    ! fall velocity constants:
 
    real, parameter :: vconr = 2503.23638966667
    real, parameter :: normr = 25132741228.7183
    real, parameter :: thr = 1.e-8
 
    logical :: no_fall
 
    dt5 = 0.5 * dt
 
    ! -----------------------------------------------------------------------
    ! -----------------------------------------------------------------------
 
    m1_rain(:) = 0.
 
    call check_column (ktop, kbot, qr, no_fall)
 
    if (no_fall) then
        vtr(:) = vf_min
        r1 = 0.
    else
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        if (const_vr) then
            vtr(:) = vr_fac ! ifs_2016: 4.0
        else
            do k = ktop, kbot
                qden = qr(k) * den(k)
                if (qr(k) < thr) then
                    vtr(k) = vr_min
                else
                    vtr(k) = vr_fac * vconr * sqrt(min(10., sfcrho / den(k))) * &
                        exp(0.2 * log(qden / normr))
                    vtr(k) = min(vr_max, max(vr_min, vtr(k)))
                endif
            enddo
        endif
 
        ze(kbot + 1) = zs
        do k = kbot, ktop, - 1
            ze(k) = ze(k + 1) - dz(k) ! dz < 0
        enddo
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        ! if (.not. fast_sat_adj) &
        call revap_racc (ktop, kbot, dt5, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
 
        if (do_sedi_w) then
            do k = ktop, kbot
                dm(k) = dp(k) * (1. + qv(k) + ql(k) + qr(k) + qi(k) + qs(k) + qg(k))
            enddo
        endif
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        if (use_ppm) then
            zt(ktop) = ze(ktop)
            do k = ktop + 1, kbot
                zt(k) = ze(k) - dt5 * (vtr(k - 1) + vtr(k))
            enddo
            zt(kbot + 1) = zs - dt * vtr(kbot)
 
            do k = ktop, kbot
                if (zt(k + 1) >= zt(k)) zt(k + 1) = zt(k) - dz_min
            enddo
            call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qr, r1, m1_rain, mono_prof)
        else
            call implicit_fall (dt, ktop, kbot, ze, vtr, dp, qr, r1, m1_rain)
        endif
 
        ! -----------------------------------------------------------------------
        !  - Calculate vertical velocity transportation during sedimentation.
        ! (do_sedi_w =.true. to turn on vertical motion tranport during sedimentation
        ! .false. by default)
        ! -----------------------------------------------------------------------
 
        if (do_sedi_w) then
            w1(ktop) = (dm(ktop) * w1(ktop) + m1_rain(ktop) * vtr(ktop)) / (dm(ktop) - m1_rain(ktop))
            do k = ktop + 1, kbot
                w1(k) = (dm(k) * w1(k) - m1_rain(k - 1) * vtr(k - 1) + m1_rain(k) * vtr(k)) &
                     / (dm(k) + m1_rain(k - 1) - m1_rain(k))
            enddo
        endif
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        if (do_sedi_heat) &
            call sedi_heat (ktop, kbot, dp, m1_rain, dz, tz, qv, ql, qr, qi, qs, qg, c_liq)
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        call revap_racc (ktop, kbot, dt5, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
 
    endif
 
    ! -----------------------------------------------------------------------
    ! -----------------------------------------------------------------------
 
    if (irain_f /= 0) then
 
        ! -----------------------------------------------------------------------
        ! no subgrid varaibility
        ! -----------------------------------------------------------------------
 
        do k = ktop, kbot
            qc0 = fac_rc * ccn(k)
            if (tz(k) > t_wfr) then
                if (use_ccn) then
                    ! -----------------------------------------------------------------------
                    ! ccn is formulted as ccn = ccn_surface * (den / den_surface)
                    ! -----------------------------------------------------------------------
                    qc = qc0
                else
                    qc = qc0 / den(k)
                endif
                dq = ql(k) - qc
                if (dq > 0.) then
                    sink = min(dq, dt * c_praut(k) * den(k) * exp(so3 * log(ql(k))))
                    ql(k) = ql(k) - sink
                    qr(k) = qr(k) + sink
                endif
            endif
        enddo
 
    else
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        call linear_prof (kbot - ktop + 1, ql(ktop), dl(ktop), z_slope_liq, h_var)
 
        do k = ktop, kbot
            qc0 = fac_rc * ccn(k)
            if (tz(k) > t_wfr + dt_fr) then
                dl(k) = min(max(1.e-6, dl(k)), 0.5 * ql(k))
                ! --------------------------------------------------------------------
                ! as in klein's gfdl am2 stratiform scheme (with subgrid variations)
                ! --------------------------------------------------------------------
                if (use_ccn) then
                    ! --------------------------------------------------------------------
                    ! ccn is formulted as ccn = ccn_surface * (den / den_surface)
                    ! --------------------------------------------------------------------
                    qc = qc0
                else
                    qc = qc0 / den(k)
                endif
                dq = 0.5 * (ql(k) + dl(k) - qc)
                ! --------------------------------------------------------------------
                ! dq = dl if qc == q_minus = ql - dl
                ! dq = 0 if qc == q_plus = ql + dl
                ! --------------------------------------------------------------------
                if (dq > 0.) then ! q_plus > qc
                    ! --------------------------------------------------------------------
                    ! --------------------------------------------------------------------
                    sink = min(1., dq / dl(k)) * dt * c_praut(k) * den(k) * exp(so3 * log(ql(k)))
                    ql(k) = ql(k) - sink
                    qr(k) = qr(k) + sink
                endif
            endif
        enddo
    endif
 
end subroutine warm_rain
 
! -----------------------------------------------------------------------
subroutine revap_racc (ktop, kbot, dt, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
 
    implicit none
 
    integer, intent (in) :: ktop, kbot
 
    real, intent (in) :: dt ! time step (s)
    real, intent (in) :: rh_rain, h_var
 
    real, intent (in), dimension (ktop:kbot) :: den, denfac
 
    real, intent (inout), dimension (ktop:kbot) :: tz, qv, qr, ql, qi, qs, qg
 
    real, dimension (ktop:kbot) :: lhl, cvm, q_liq, q_sol, lcpk
 
    real :: dqv, qsat, dqsdt, evap, t2, qden, q_plus, q_minus, sink
    real :: qpz, dq, dqh, tin
 
    integer :: k
 
    do k = ktop, kbot
 
        if (tz(k) > t_wfr .and. qr(k) > qrmin) then
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            lhl(k) = lv00 + d0_vap * tz(k)
            q_liq(k) = ql(k) + qr(k)
            q_sol(k) = qi(k) + qs(k) + qg(k)
            cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
            lcpk(k) = lhl(k) / cvm(k)
 
            tin = tz(k) - lcpk(k) * ql(k) ! presence of clouds suppresses the rain evap
            qpz = qv(k) + ql(k)
            qsat = wqs2(tin, den(k), dqsdt)
            dqh = max(ql(k), h_var * max(qpz, qcmin))
            dqh = min(dqh, 0.2 * qpz) ! new limiter
            dqv = qsat - qv(k) ! use this to prevent super - sat the gird box
            q_minus = qpz - dqh
            q_plus = qpz + dqh
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            if (dqv > qvmin .and. qsat > q_minus) then
                if (qsat > q_plus) then
                    dq = qsat - qpz
                else
                    ! -----------------------------------------------------------------------
                    ! q_minus < qsat < q_plus
                    ! dq == dqh if qsat == q_minus
                    ! -----------------------------------------------------------------------
                    dq = 0.25 * (q_minus - qsat) ** 2 / dqh
                endif
                qden = qr(k) * den(k)
                t2 = tin * tin
                evap = crevp(1) * t2 * dq * (crevp(2) * sqrt(qden) + crevp(3) * &
                    exp(0.725 * log(qden))) / (crevp(4) * t2 + crevp(5) * qsat * den(k))
                evap = min(qr(k), dt * evap, dqv / (1. + lcpk(k) * dqsdt))
                ! -----------------------------------------------------------------------
                ! alternative minimum evap in dry environmental air
                ! sink = min (qr (k), dim (rh_rain * qsat, qv (k)) / (1. + lcpk (k) * dqsdt))
                ! evap = max (evap, sink)
                ! -----------------------------------------------------------------------
                qr(k) = qr(k) - evap
                qv(k) = qv(k) + evap
                q_liq(k) = q_liq(k) - evap
                cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
                tz(k) = tz(k) - evap * lhl(k) / cvm(k)
            endif
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            ! if (qr (k) > qrmin .and. ql (k) > 1.e-7 .and. qsat < q_plus) then
            if (qr(k) > qrmin .and. ql(k) > 1.e-6 .and. qsat < q_minus) then
                sink = dt * denfac(k) * cracw * exp(0.95 * log(qr(k) * den(k)))
                sink = sink / (1. + sink) * ql(k)
                ql(k) = ql(k) - sink
                qr(k) = qr(k) + sink
            endif
 
        endif ! warm - rain
    enddo
 
end subroutine revap_racc
 
! -----------------------------------------------------------------------
! qi -- > ql & ql -- > qr
! edges: qe == qbar + / - dm
subroutine linear_prof (km, q, dm, z_var, h_var)
 
    implicit none
 
    integer, intent (in) :: km
 
    real, intent (in) :: q (km), h_var
 
    real, intent (out) :: dm (km)
 
    logical, intent (in) :: z_var
 
    real :: dq (km)
 
    integer :: k
 
    if (z_var) then
        do k = 2, km
            dq(k) = 0.5 * (q(k) - q(k - 1))
        enddo
        dm(1) = 0.
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        do k = 2, km - 1
            dm(k) = 0.5 * min(abs(dq(k) + dq(k + 1)), 0.5 * q(k))
            if (dq(k) * dq(k + 1) <= 0.) then
                if (dq(k) > 0.) then ! local max
                    dm(k) = min(dm(k), dq(k), - dq(k + 1))
                else
                    dm(k) = 0.
                endif
            endif
        enddo
        dm(km) = 0.
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        do k = 1, km
            dm(k) = max(dm(k), qvmin, h_var * q(k))
        enddo
    else
        do k = 1, km
            dm(k) = max(qvmin, h_var * q(k))
        enddo
    endif
 
end subroutine linear_prof
 
! =======================================================================
subroutine icloud (ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, &
        den, denfac, vts, vtg, vtr, qak, rh_adj, rh_rain, dts, h_var)
 
    implicit none
 
    integer, intent (in) :: ktop, kbot
 
    real, intent (in), dimension (ktop:kbot) :: p1, dp1, den, denfac, vts, vtg, vtr
 
    real, intent (inout), dimension (ktop:kbot) :: tzk, qvk, qlk, qrk, qik, qsk, qgk, qak
 
    real, intent (in) :: rh_adj, rh_rain, dts, h_var
 
    real, dimension (ktop:kbot) :: lcpk, icpk, tcpk, di, lhl, lhi
    real, dimension (ktop:kbot) :: cvm, q_liq, q_sol
 
    real :: rdts, fac_g2v, fac_v2g, fac_i2s, fac_imlt
    real :: tz, qv, ql, qr, qi, qs, qg, melt
    real :: pracs, psacw, pgacw, psacr, pgacr, pgaci, praci, psaci
    real :: pgmlt, psmlt, pgfr, pgaut, psaut, pgsub
    real :: tc, tsq, dqs0, qden, qim, qsm
    real :: dt5, factor, sink, qi_crt
    real :: tmp, qsw, qsi, dqsdt, dq
    real :: dtmp, qc, q_plus, q_minus
 
    integer :: k
 
    dt5 = 0.5 * dts
 
    rdts = 1. / dts
 
    ! -----------------------------------------------------------------------
    ! -----------------------------------------------------------------------
 
    fac_i2s = 1. - exp(- dts / tau_i2s)
    fac_g2v = 1. - exp(- dts / tau_g2v)
    fac_v2g = 1. - exp(- dts / tau_v2g)
 
    fac_imlt = 1. - exp(- dt5 / tau_imlt)
 
    ! -----------------------------------------------------------------------
    ! -----------------------------------------------------------------------
 
    do k = ktop, kbot
        lhi(k) = li00 + dc_ice * tzk(k)
        q_liq(k) = qlk(k) + qrk(k)
        q_sol(k) = qik(k) + qsk(k) + qgk(k)
        cvm(k) = c_air + qvk(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
        icpk(k) = lhi(k) / cvm(k)
    enddo
 
    ! -----------------------------------------------------------------------
    ! sources of cloud ice: pihom, cold rain, and the sat_adj
    ! (initiation plus deposition)
    ! sources of snow: cold rain, auto conversion + accretion (from cloud ice)
    ! sat_adj (deposition; requires pre - existing snow) ; initial snow comes from auto conversion
    ! -----------------------------------------------------------------------
 
    do k = ktop, kbot
        if (tzk(k) > tice .and. qik(k) > qcmin) then
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            melt = min(qik(k), fac_imlt * (tzk(k) - tice) / icpk(k))
            tmp = min(melt, dim(ql_mlt, qlk(k))) ! max ql amount
            qlk(k) = qlk(k) + tmp
            qrk(k) = qrk(k) + melt - tmp
            qik(k) = qik(k) - melt
            q_liq(k) = q_liq(k) + melt
            q_sol(k) = q_sol(k) - melt
            cvm(k) = c_air + qvk(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
            tzk(k) = tzk(k) - melt * lhi(k) / cvm(k)
 
        elseif (tzk(k) < t_wfr .and. qlk(k) > qcmin) then
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            dtmp = t_wfr - tzk(k)
            factor = min(1., dtmp / dt_fr)
            sink = min(qlk(k) * factor, dtmp / icpk(k))
            qi_crt = qi_gen * min(qi_lim, 0.1 * (tice - tzk(k))) / den(k)
            tmp = min(sink, dim(qi_crt, qik(k)))
            qlk(k) = qlk(k) - sink
            qsk(k) = qsk(k) + sink - tmp
            qik(k) = qik(k) + tmp
            q_liq(k) = q_liq(k) - sink
            q_sol(k) = q_sol(k) + sink
            cvm(k) = c_air + qvk(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
            tzk(k) = tzk(k) + sink * lhi(k) / cvm(k)
 
        endif
    enddo
 
    ! -----------------------------------------------------------------------
    ! -----------------------------------------------------------------------
 
    call linear_prof (kbot - ktop + 1, qik(ktop), di(ktop), z_slope_ice, h_var)
 
    ! -----------------------------------------------------------------------
    ! -----------------------------------------------------------------------
 
    do k = ktop, kbot
        lhl(k) = lv00 + d0_vap * tzk(k)
        lhi(k) = li00 + dc_ice * tzk(k)
        lcpk(k) = lhl(k) / cvm(k)
        icpk(k) = lhi(k) / cvm(k)
        tcpk(k) = lcpk(k) + icpk(k)
    enddo
 
    do k = ktop, kbot
 
        ! -----------------------------------------------------------------------
        ! do nothing above p_min
        ! -----------------------------------------------------------------------
 
        if (p1(k) < p_min) cycle
 
        tz = tzk(k)
        qv = qvk(k)
        ql = qlk(k)
        qi = qik(k)
        qr = qrk(k)
        qs = qsk(k)
        qg = qgk(k)
 
        pgacr = 0.
        pgacw = 0.
        tc = tz - tice
 
        if (tc .ge. 0.) then
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            dqs0 = ces0 / p1(k) - qv
 
            if (qs > qcmin) then
 
                ! -----------------------------------------------------------------------
                ! only rate is used (for snow melt) since tc > 0.
                ! -----------------------------------------------------------------------
 
                if (ql > qrmin) then
                    factor = denfac(k) * csacw * exp(0.8125 * log(qs * den(k)))
                    psacw = factor / (1. + dts * factor) * ql ! rate
                else
                    psacw = 0.
                endif
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
 
                if (qr > qrmin) then
                    psacr = min(acr3d(vts(k), vtr(k), qr, qs, csacr, acco(1, 2), &
                        den(k)), qr * rdts)
                    pracs = acr3d(vtr(k), vts(k), qs, qr, cracs, acco(1, 1), den(k))
                else
                    psacr = 0.
                    pracs = 0.
                endif
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
 
                psmlt = max(0., smlt(tc, dqs0, qs * den(k), psacw, psacr, csmlt, &
                    den(k), denfac(k)))
                sink = min(qs, dts * (psmlt + pracs), tc / icpk(k))
                qs = qs - sink
                ! sjl, 20170321:
                tmp = min(sink, dim(qs_mlt, ql)) ! max ql due to snow melt
                ql = ql + tmp
                qr = qr + sink - tmp
                ! qr = qr + sink
                ! sjl, 20170321:
                q_liq(k) = q_liq(k) + sink
                q_sol(k) = q_sol(k) - sink
                cvm(k) = c_air + qv * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
                tz = tz - sink * lhi(k) / cvm(k)
                tc = tz - tice
 
            endif
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            lhi(k) = li00 + dc_ice * tz
            icpk(k) = lhi(k) / cvm(k)
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            if (qg > qcmin .and. tc > 0.) then
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
 
                if (qr > qrmin) &
                    pgacr = min(acr3d(vtg(k), vtr(k), qr, qg, cgacr, acco(1, 3), &
                        den(k)), rdts * qr)
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
 
                qden = qg * den(k)
                if (ql > qrmin) then
                    factor = cgacw * qden / sqrt(den(k) * sqrt(sqrt(qden)))
                    pgacw = factor / (1. + dts * factor) * ql ! rate
                endif
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
 
                pgmlt = dts * gmlt(tc, dqs0, qden, pgacw, pgacr, cgmlt, den(k))
                pgmlt = min(max(0., pgmlt), qg, tc / icpk(k))
                qg = qg - pgmlt
                qr = qr + pgmlt
                q_liq(k) = q_liq(k) + pgmlt
                q_sol(k) = q_sol(k) - pgmlt
                cvm(k) = c_air + qv * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
                tz = tz - pgmlt * lhi(k) / cvm(k)
 
            endif
 
        else
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            if (qi > 3.e-7) then ! cloud ice sink terms
 
                if (qs > 1.e-7) then
                    ! -----------------------------------------------------------------------
                    ! sjl added (following lin eq. 23) the temperature dependency
                    ! to reduce accretion, use esi = exp (0.05 * tc) as in hong et al 2004
                    ! -----------------------------------------------------------------------
                    factor = dts * denfac(k) * csaci * exp(0.05 * tc + 0.8125 * log(qs * den(k)))
                    psaci = factor / (1. + factor) * qi
                else
                    psaci = 0.
                endif
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
                if (qi0_crt < 0.) then
                    qim = - qi0_crt
                else
                    qim = qi0_crt / den(k)
                endif
                ! -----------------------------------------------------------------------
                ! assuming linear subgrid vertical distribution of cloud ice
                ! the mismatch computation following lin et al. 1994, mwr
                ! -----------------------------------------------------------------------
 
                if (const_vi) then
                    tmp = fac_i2s
                else
                    tmp = fac_i2s * exp(0.025 * tc)
                endif
 
                di(k) = max(di(k), qrmin)
                q_plus = qi + di(k)
                if (q_plus > (qim + qrmin)) then
                    if (qim > (qi - di(k))) then
                        dq = (0.25 * (q_plus - qim) ** 2) / di(k)
                    else
                        dq = qi - qim
                    endif
                    psaut = tmp * dq
                else
                    psaut = 0.
                endif
                ! -----------------------------------------------------------------------
                ! sink is no greater than 75% of qi
                ! -----------------------------------------------------------------------
                sink = min(0.75 * qi, psaci + psaut)
                qi = qi - sink
                qs = qs + sink
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
 
                if (qg > 1.e-6) then
                    ! -----------------------------------------------------------------------
                    ! factor = dts * cgaci / sqrt (den (k)) * exp (0.05 * tc + 0.875 * log (qg * den (k)))
                    ! simplified form: remove temp dependency & set the exponent "0.875" -- > 1
                    ! -----------------------------------------------------------------------
                    factor = dts * cgaci * sqrt(den(k)) * qg
                    pgaci = factor / (1. + factor) * qi
                    qi = qi - pgaci
                    qg = qg + pgaci
                endif
 
            endif
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            ! -----------------------------------------------------------------------
            ! rain to ice, snow, graupel processes:
            ! -----------------------------------------------------------------------
 
            tc = tz - tice
 
            if (qr > 1.e-7 .and. tc < 0.) then
 
                ! -----------------------------------------------------------------------
                ! * sink * terms to qr: psacr + pgfr
                ! source terms to qs: psacr
                ! source terms to qg: pgfr
                ! -----------------------------------------------------------------------
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
 
                if (qs > 1.e-7) then ! if snow exists
                    psacr = dts * acr3d(vts(k), vtr(k), qr, qs, csacr, acco(1, 2), den(k))
                else
                    psacr = 0.
                endif
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
 
                pgfr = dts * cgfr(1) / den(k) * (exp(- cgfr(2) * tc) - 1.) * &
                    exp(1.75 * log(qr * den(k)))
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
 
                sink = psacr + pgfr
                factor = min(sink, qr, - tc / icpk(k)) / max(sink, qrmin)
 
                psacr = factor * psacr
                pgfr = factor * pgfr
 
                sink = psacr + pgfr
                qr = qr - sink
                qs = qs + psacr
                qg = qg + pgfr
                q_liq(k) = q_liq(k) - sink
                q_sol(k) = q_sol(k) + sink
                cvm(k) = c_air + qv * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
                tz = tz + sink * lhi(k) / cvm(k)
 
            endif
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            lhi(k) = li00 + dc_ice * tz
            icpk(k) = lhi(k) / cvm(k)
 
            ! -----------------------------------------------------------------------
            ! -----------------------------------------------------------------------
 
            if (qs > 1.e-7) then
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
 
                if (qg > qrmin) then
                    sink = dts * acr3d(vtg(k), vts(k), qs, qg, cgacs, acco(1, 4), den(k))
                else
                    sink = 0.
                endif
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
 
                qsm = qs0_crt / den(k)
                if (qs > qsm) then
                    factor = dts * 1.e-3 * exp(0.09 * (tz - tice))
                    sink = sink + factor / (1. + factor) * (qs - qsm)
                endif
                sink = min(qs, sink)
                qs = qs - sink
                qg = qg + sink
 
            endif ! snow existed
 
            if (qg > 1.e-7 .and. tz < tice0) then
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
 
                if (ql > 1.e-6) then
                    qden = qg * den(k)
                    factor = dts * cgacw * qden / sqrt(den(k) * sqrt(sqrt(qden)))
                    pgacw = factor / (1. + factor) * ql
                else
                    pgacw = 0.
                endif
 
                ! -----------------------------------------------------------------------
                ! -----------------------------------------------------------------------
 
                if (qr > 1.e-6) then
                    pgacr = min(dts * acr3d(vtg(k), vtr(k), qr, qg, cgacr, acco(1, 3), &
                        den(k)), qr)
                else
                    pgacr = 0.
                endif
 
                sink = pgacr + pgacw
                factor = min(sink, dim(tice, tz) / icpk(k)) / max(sink, qrmin)
                pgacr = factor * pgacr
                pgacw = factor * pgacw
 
                sink = pgacr + pgacw
                qg = qg + sink
                qr = qr - pgacr
                ql = ql - pgacw
                q_liq(k) = q_liq(k) - sink
                q_sol(k) = q_sol(k) + sink
                cvm(k) = c_air + qv * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
                tz = tz + sink * lhi(k) / cvm(k)
 
            endif
 
        endif
 
        tzk(k) = tz
        qvk(k) = qv
        qlk(k) = ql
        qik(k) = qi
        qrk(k) = qr
        qsk(k) = qs
        qgk(k) = qg
 
    enddo
 
    ! -----------------------------------------------------------------------
    ! -----------------------------------------------------------------------
 
    call subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, rh_adj, tzk, qvk, &
        qlk, qrk, qik, qsk, qgk, qak, h_var, rh_rain)
 
end subroutine icloud
 
! =======================================================================
subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, rh_adj, tz, qv, &
    ql, qr, qi, qs, qg, qa, h_var, rh_rain)
 
    implicit none
 
    integer, intent (in) :: ktop, kbot
 
    real, intent (in), dimension (ktop:kbot) :: p1, den, denfac
 
    real, intent (in) :: dts, rh_adj, h_var, rh_rain
 
    real, intent (inout), dimension (ktop:kbot) :: tz, qv, ql, qr, qi, qs, qg, qa
 
    real, dimension (ktop:kbot) :: lcpk, icpk, tcpk, tcp3, lhl, lhi
    real, dimension (ktop:kbot) :: cvm, q_liq, q_sol, q_cond
 
    real :: fac_v2l, fac_l2v
 
    real :: pidep, qi_crt
 
    ! -----------------------------------------------------------------------
    ! qstar over water may be accurate only down to - 80 deg c with ~10% uncertainty
    ! must not be too large to allow psc
    ! -----------------------------------------------------------------------
 
    real :: rh, rqi, tin, qsw, qsi, qpz, qstar
    real :: dqsdt, dwsdt, dq, dq0, factor, tmp
    real :: q_plus, q_minus, dt_evap, dt_pisub
    real :: evap, sink, tc, pisub, q_adj, dtmp
    real :: pssub, pgsub, tsq, qden, fac_g2v, fac_v2g
 
    integer :: k
 
    if (fast_sat_adj) then
        dt_evap = 0.5 * dts
    else
        dt_evap = dts
    endif
 
    ! -----------------------------------------------------------------------
    ! -----------------------------------------------------------------------
 
    fac_v2l = 1. - exp(- dt_evap / tau_v2l)
    fac_l2v = 1. - exp(- dt_evap / tau_l2v)
 
    fac_g2v = 1. - exp(- dts / tau_g2v)
    fac_v2g = 1. - exp(- dts / tau_v2g)
 
    ! -----------------------------------------------------------------------
    ! -----------------------------------------------------------------------
 
    do k = ktop, kbot
        lhl(k) = lv00 + d0_vap * tz(k)
        lhi(k) = li00 + dc_ice * tz(k)
        q_liq(k) = ql(k) + qr(k)
        q_sol(k) = qi(k) + qs(k) + qg(k)
        cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
        lcpk(k) = lhl(k) / cvm(k)
        icpk(k) = lhi(k) / cvm(k)
        tcpk(k) = lcpk(k) + icpk(k)
        tcp3(k) = lcpk(k) + icpk(k) * min(1., dim(tice, tz(k)) / (tice - t_wfr))
    enddo
 
    do k = ktop, kbot
 
        if (p1(k) < p_min) cycle
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        if (tz(k) < t_min) then
            sink = dim(qv(k), 1.e-7)
            qv(k) = qv(k) - sink
            qi(k) = qi(k) + sink
            q_sol(k) = q_sol(k) + sink
            cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
            tz(k) = tz(k) + sink * (lhl(k) + lhi(k)) / cvm(k)
            if (.not. do_qa) qa(k) = qa(k) + 1. ! air fully saturated; 100 % cloud cover
            cycle
        endif
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        lhl(k) = lv00 + d0_vap * tz(k)
        lhi(k) = li00 + dc_ice * tz(k)
        lcpk(k) = lhl(k) / cvm(k)
        icpk(k) = lhi(k) / cvm(k)
        tcpk(k) = lcpk(k) + icpk(k)
        tcp3(k) = lcpk(k) + icpk(k) * min(1., dim(tice, tz(k)) / (tice - t_wfr))
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        qpz = qv(k) + ql(k) + qi(k)
        tin = tz(k) - (lhl(k) * (ql(k) + qi(k)) + lhi(k) * qi(k)) / (c_air + &
            qpz * c_vap + qr(k) * c_liq + (qs(k) + qg(k)) * c_ice)
        if (tin > t_sub + 6.) then
            rh = qpz / iqs1(tin, den(k))
            if (rh < rh_adj) then ! qpz / rh_adj < qs
                tz(k) = tin
                qv(k) = qpz
                ql(k) = 0.
                qi(k) = 0.
                cycle ! cloud free
            endif
        endif
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        qsw = wqs2(tz(k), den(k), dwsdt)
        dq0 = qsw - qv(k)
        if (dq0 > 0.) then
            ! SJL 20170703 added ql factor to prevent the situation of high ql and low RH
            ! factor = min (1., fac_l2v * sqrt (max (0., ql (k)) / 1.e-5) * 10. * dq0 / qsw)
            ! factor = fac_l2v
            ! factor = 1
            factor = min(1., fac_l2v * (10. * dq0 / qsw)) ! the rh dependent factor = 1 at 90%
            evap = min(ql(k), factor * dq0 / (1. + tcp3(k) * dwsdt))
        else ! condensate all excess vapor into cloud water
            ! -----------------------------------------------------------------------
            ! evap = fac_v2l * dq0 / (1. + tcp3 (k) * dwsdt)
            ! sjl, 20161108
            ! -----------------------------------------------------------------------
            evap = dq0 / (1. + tcp3(k) * dwsdt)
        endif
        qv(k) = qv(k) + evap
        ql(k) = ql(k) - evap
        q_liq(k) = q_liq(k) - evap
        cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
        tz(k) = tz(k) - evap * lhl(k) / cvm(k)
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        lhi(k) = li00 + dc_ice * tz(k)
        icpk(k) = lhi(k) / cvm(k)
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        dtmp = t_wfr - tz(k) ! [ - 40, - 48]
        if (dtmp > 0. .and. ql(k) > qcmin) then
            sink = min(ql(k), ql(k) * dtmp * 0.125, dtmp / icpk(k))
            ql(k) = ql(k) - sink
            qi(k) = qi(k) + sink
            q_liq(k) = q_liq(k) - sink
            q_sol(k) = q_sol(k) + sink
            cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
            tz(k) = tz(k) + sink * lhi(k) / cvm(k)
        endif
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        lhi(k) = li00 + dc_ice * tz(k)
        icpk(k) = lhi(k) / cvm(k)
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        if (fast_sat_adj) then
            dt_pisub = 0.5 * dts
        else
            dt_pisub = dts
            tc = tice - tz(k)
            if (ql(k) > qrmin .and. tc > 0.) then
                sink = 3.3333e-10 * dts * (exp(0.66 * tc) - 1.) * den(k) * ql(k) * ql(k)
                sink = min(ql(k), tc / icpk(k), sink)
                ql(k) = ql(k) - sink
                qi(k) = qi(k) + sink
                q_liq(k) = q_liq(k) - sink
                q_sol(k) = q_sol(k) + sink
                cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
                tz(k) = tz(k) + sink * lhi(k) / cvm(k)
            endif ! significant ql existed
        endif
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        lhl(k) = lv00 + d0_vap * tz(k)
        lhi(k) = li00 + dc_ice * tz(k)
        lcpk(k) = lhl(k) / cvm(k)
        icpk(k) = lhi(k) / cvm(k)
        tcpk(k) = lcpk(k) + icpk(k)
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        if (tz(k) < tice) then
            qsi = iqs2(tz(k), den(k), dqsdt)
            dq = qv(k) - qsi
            sink = dq / (1. + tcpk(k) * dqsdt)
            if (qi(k) > qrmin) then
                ! eq 9, hong et al. 2004, mwr
                ! for a and b, see dudhia 1989: page 3103 eq (b7) and (b8)
                pidep = dt_pisub * dq * 349138.78 * exp(0.875 * log(qi(k) * den(k))) &
                     / (qsi * den(k) * lat2 / (0.0243 * rvgas * tz(k) ** 2) + 4.42478e4)
            else
                pidep = 0.
            endif
            if (dq > 0.) then ! vapor - > ice
                tmp = tice - tz(k)
                ! 20160912: the following should produce more ice at higher altitude
                ! qi_crt = 4.92e-11 * exp (1.33 * log (1.e3 * exp (0.1 * tmp))) / den (k)
                qi_crt = qi_gen * min(qi_lim, 0.1 * tmp) / den(k)
                sink = min(sink, max(qi_crt - qi(k), pidep), tmp / tcpk(k))
            else ! ice -- > vapor
                pidep = pidep * min(1., dim(tz(k), t_sub) * 0.2)
                sink = max(pidep, sink, - qi(k))
            endif
            qv(k) = qv(k) - sink
            qi(k) = qi(k) + sink
            q_sol(k) = q_sol(k) + sink
            cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
            tz(k) = tz(k) + sink * (lhl(k) + lhi(k)) / cvm(k)
        endif
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        lhl(k) = lv00 + d0_vap * tz(k)
        lhi(k) = li00 + dc_ice * tz(k)
        lcpk(k) = lhl(k) / cvm(k)
        icpk(k) = lhi(k) / cvm(k)
        tcpk(k) = lcpk(k) + icpk(k)
 
        ! -----------------------------------------------------------------------
        ! this process happens for all temp rage
        ! -----------------------------------------------------------------------
 
        if (qs(k) > qrmin) then
            qsi = iqs2(tz(k), den(k), dqsdt)
            qden = qs(k) * den(k)
            tmp = exp(0.65625 * log(qden))
            tsq = tz(k) * tz(k)
            dq = (qsi - qv(k)) / (1. + tcpk(k) * dqsdt)
            pssub = cssub(1) * tsq * (cssub(2) * sqrt(qden) + cssub(3) * tmp * &
                sqrt(denfac(k))) / (cssub(4) * tsq + cssub(5) * qsi * den(k))
            pssub = (qsi - qv(k)) * dts * pssub
            if (pssub > 0.) then ! qs -- > qv, sublimation
                pssub = min(pssub * min(1., dim(tz(k), t_sub) * 0.2), qs(k))
            else
                if (tz(k) > tice) then
                    pssub = 0. ! no deposition
                else
                    pssub = max(pssub, dq, (tz(k) - tice) / tcpk(k))
                endif
            endif
            qs(k) = qs(k) - pssub
            qv(k) = qv(k) + pssub
            q_sol(k) = q_sol(k) - pssub
            cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
            tz(k) = tz(k) - pssub * (lhl(k) + lhi(k)) / cvm(k)
        endif
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        lhl(k) = lv00 + d0_vap * tz(k)
        lhi(k) = li00 + dc_ice * tz(k)
        lcpk(k) = lhl(k) / cvm(k)
        icpk(k) = lhi(k) / cvm(k)
        tcpk(k) = lcpk(k) + icpk(k)
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        if (qg(k) > qrmin) then
            qsi = iqs2(tz(k), den(k), dqsdt)
            dq = (qv(k) - qsi) / (1. + tcpk(k) * dqsdt)
            pgsub = (qv(k) / qsi - 1.) * qg(k)
            if (pgsub > 0.) then ! deposition
                if (tz(k) > tice) then
                    pgsub = 0. ! no deposition
                else
                    pgsub = min(fac_v2g * pgsub, 0.2 * dq, ql(k) + qr(k), &
                        (tice - tz(k)) / tcpk(k))
                endif
            else ! submilation
                pgsub = max(fac_g2v * pgsub, dq) * min(1., dim(tz(k), t_sub) * 0.1)
            endif
            qg(k) = qg(k) + pgsub
            qv(k) = qv(k) - pgsub
            q_sol(k) = q_sol(k) + pgsub
            cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
            tz(k) = tz(k) + pgsub * (lhl(k) + lhi(k)) / cvm(k)
        endif
 
#ifdef USE_MIN_EVAP
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        lhl(k) = lv00 + d0_vap * tz(k)
        lcpk(k) = lhl(k) / cvm(k)
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        if (qr(k) > qcmin) then
            qsw = wqs2(tz(k), den(k), dqsdt)
            sink = min(qr(k), dim(rh_rain * qsw, qv(k)) / (1. + lcpk(k) * dqsdt))
            qv(k) = qv(k) + sink
            qr(k) = qr(k) - sink
            q_liq(k) = q_liq(k) - sink
            cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
            tz(k) = tz(k) - sink * lhl(k) / cvm(k)
        endif
#endif
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        lhl(k) = lv00 + d0_vap * tz(k)
        cvm(k) = c_air + (qv(k) + q_liq(k) + q_sol(k)) * c_vap
        lcpk(k) = lhl(k) / cvm(k)
 
        ! -----------------------------------------------------------------------
        ! compute cloud fraction
        ! -----------------------------------------------------------------------
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        if (do_qa) cycle
 
        if (rad_snow) then
            q_sol(k) = qi(k) + qs(k)
        else
            q_sol(k) = qi(k)
        endif
        if (rad_rain) then
            q_liq(k) = ql(k) + qr(k)
        else
            q_liq(k) = ql(k)
        endif
        q_cond(k) = q_liq(k) + q_sol(k)
 
        qpz = qv(k) + q_cond(k) ! qpz is conserved
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        tin = tz(k) - (lcpk(k) * q_cond(k) + icpk(k) * q_sol(k)) ! minimum temperature
        ! tin = tz (k) - ((lv00 + d0_vap * tz (k)) * q_cond (k) + &
        ! (li00 + dc_ice * tz (k)) * q_sol (k)) / (c_air + qpz * c_vap)
 
        ! -----------------------------------------------------------------------
        ! determine saturated specific humidity
        ! -----------------------------------------------------------------------
 
        if (tin <= t_wfr) then
            ! ice phase:
            qstar = iqs1(tin, den(k))
        elseif (tin >= tice) then
            ! liquid phase:
            qstar = wqs1(tin, den(k))
        else
            ! mixed phase:
            qsi = iqs1(tin, den(k))
            qsw = wqs1(tin, den(k))
            if (q_cond(k) > 3.e-6) then
                rqi = q_sol(k) / q_cond(k)
            else
                ! -----------------------------------------------------------------------
                ! mostly liquid water q_cond (k) at initial cloud development stage
                ! -----------------------------------------------------------------------
                rqi = (tice - tin) / (tice - t_wfr)
            endif
            qstar = rqi * qsi + (1. - rqi) * qsw
        endif
 
        ! -----------------------------------------------------------------------
        ! -----------------------------------------------------------------------
 
        if (qpz > qrmin) then
            ! partial cloudiness by pdf:
            dq = max(qcmin, h_var * qpz)
            q_plus = qpz + dq ! cloud free if qstar > q_plus
            q_minus = qpz - dq
            if (qstar < q_minus) then
                qa(k) = qa(k) + 1. ! air fully saturated; 100 % cloud cover
            elseif (qstar < q_plus .and. q_cond(k) > qc_crt) then
                qa(k) = qa(k) + (q_plus - qstar) / (dq + dq) ! partial cloud cover
                ! qa (k) = sqrt (qa (k) + (q_plus - qstar) / (dq + dq))
            endif
        endif
 
    enddo
 
end subroutine subgrid_z_proc
 
! =======================================================================
subroutine revap_rac1 (hydrostatic, is, ie, dt, tz, qv, ql, qr, qi, qs, qg, den, hvar)
 
    implicit none
 
    logical, intent (in) :: hydrostatic
 
    integer, intent (in) :: is, ie
 
    real, intent (in) :: dt ! time step (s)
 
    real, intent (in), dimension (is:ie) :: den, hvar, qi, qs, qg
 
    real, intent (inout), dimension (is:ie) :: tz, qv, qr, ql
 
    real, dimension (is:ie) :: lcp2, denfac, q_liq, q_sol, cvm, lhl
 
    real :: dqv, qsat, dqsdt, evap, qden, q_plus, q_minus, sink
    real :: tin, t2, qpz, dq, dqh
 
    integer :: i
 
    ! -----------------------------------------------------------------------
    ! define latend heat coefficient
    ! -----------------------------------------------------------------------
 
    do i = is, ie
        lhl(i) = lv00 + d0_vap * tz(i)
        q_liq(i) = ql(i) + qr(i)
        q_sol(i) = qi(i) + qs(i) + qg(i)
        cvm(i) = c_air + qv(i) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
        lcp2(i) = lhl(i) / cvm(i)
        ! denfac (i) = sqrt (sfcrho / den (i))
    enddo
 
    do i = is, ie
        if (qr(i) > qrmin .and. tz(i) > t_wfr) then
            qpz = qv(i) + ql(i)
            tin = tz(i) - lcp2(i) * ql(i) ! presence of clouds suppresses the rain evap
            qsat = wqs2(tin, den(i), dqsdt)
            dqh = max(ql(i), hvar(i) * max(qpz, qcmin))
            dqv = qsat - qv(i)
            q_minus = qpz - dqh
            q_plus = qpz + dqh
 
            ! -----------------------------------------------------------------------
            ! qsat must be > q_minus to activate evaporation
            ! qsat must be < q_plus to activate accretion
            ! -----------------------------------------------------------------------
 
            ! -----------------------------------------------------------------------
            ! rain evaporation
            ! -----------------------------------------------------------------------
 
            if (dqv > qvmin .and. qsat > q_minus) then
                if (qsat > q_plus) then
                    dq = qsat - qpz
                else
                    ! q_minus < qsat < q_plus
                    ! dq == dqh if qsat == q_minus
                    dq = 0.25 * (q_minus - qsat) ** 2 / dqh
                endif
                qden = qr(i) * den(i)
                t2 = tin * tin
                evap = crevp(1) * t2 * dq * (crevp(2) * sqrt(qden) + crevp(3) * exp(0.725 * log(qden))) &
                     / (crevp(4) * t2 + crevp(5) * qsat * den(i))
                evap = min(qr(i), dt * evap, dqv / (1. + lcp2(i) * dqsdt))
                qr(i) = qr(i) - evap
                qv(i) = qv(i) + evap
                q_liq(i) = q_liq(i) - evap
                cvm(i) = c_air + qv(i) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
                tz(i) = tz(i) - evap * lhl(i) / cvm(i)
            endif
 
            ! -----------------------------------------------------------------------
            ! accretion: pracc
            ! -----------------------------------------------------------------------
 
            if (qr(i) > qrmin .and. ql(i) > 1.e-8 .and. qsat < q_plus) then
                denfac(i) = sqrt(sfcrho / den(i))
                sink = dt * denfac(i) * cracw * exp(0.95 * log(qr(i) * den(i)))
                sink = sink / (1. + sink) * ql(i)
                ql(i) = ql(i) - sink
                qr(i) = qr(i) + sink
            endif
        endif
    enddo
 
end subroutine revap_rac1
 
! =======================================================================
subroutine terminal_fall (dtm, ktop, kbot, tz, qv, ql, qr, qg, qs, qi, dz, dp, &
        den, vtg, vts, vti, r1, g1, s1, i1, m1_sol, w1)
 
    implicit none
 
    integer, intent (in) :: ktop, kbot
 
    real, intent (in) :: dtm ! time step (s)
 
    real, intent (in), dimension (ktop:kbot) :: vtg, vts, vti, den, dp, dz
 
    real, intent (inout), dimension (ktop:kbot) :: qv, ql, qr, qg, qs, qi, tz, m1_sol, w1
 
    real, intent (out) :: r1, g1, s1, i1
 
    real, dimension (ktop:kbot + 1) :: ze, zt
 
    real :: qsat, dqsdt, dt5, evap, dtime
    real :: factor, frac
    real :: tmp, precip, tc, sink
 
    real, dimension (ktop:kbot) :: lcpk, icpk, cvm, q_liq, q_sol, lhl, lhi
    real, dimension (ktop:kbot) :: m1, dm
 
    real :: zs = 0.
    real :: fac_imlt
 
    integer :: k, k0, m
 
    logical :: no_fall
 
    dt5 = 0.5 * dtm
    fac_imlt = 1. - exp(- dt5 / tau_imlt)
 
    ! -----------------------------------------------------------------------
    ! define heat capacity and latend heat coefficient
    ! -----------------------------------------------------------------------
 
    do k = ktop, kbot
        m1_sol(k) = 0.
        lhl(k) = lv00 + d0_vap * tz(k)
        lhi(k) = li00 + dc_ice * tz(k)
        q_liq(k) = ql(k) + qr(k)
        q_sol(k) = qi(k) + qs(k) + qg(k)
        cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
        lcpk(k) = lhl(k) / cvm(k)
        icpk(k) = lhi(k) / cvm(k)
    enddo
 
    ! -----------------------------------------------------------------------
    ! find significant melting level
    ! -----------------------------------------------------------------------
 
    k0 = kbot
    do k = ktop, kbot - 1
        if (tz(k) > tice) then
            k0 = k
            exit
        endif
    enddo
 
    ! -----------------------------------------------------------------------
    ! melting of cloud_ice (before fall) :
    ! -----------------------------------------------------------------------
 
    do k = k0, kbot
        tc = tz(k) - tice
        if (qi(k) > qcmin .and. tc > 0.) then
            sink = min(qi(k), fac_imlt * tc / icpk(k))
            tmp = min(sink, dim(ql_mlt, ql(k)))
            ql(k) = ql(k) + tmp
            qr(k) = qr(k) + sink - tmp
            qi(k) = qi(k) - sink
            q_liq(k) = q_liq(k) + sink
            q_sol(k) = q_sol(k) - sink
            cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
            tz(k) = tz(k) - sink * lhi(k) / cvm(k)
            tc = tz(k) - tice
        endif
    enddo
 
    ! -----------------------------------------------------------------------
    ! turn off melting when cloud microphysics time step is small
    ! -----------------------------------------------------------------------
 
    if (dtm < 60.) k0 = kbot
 
    ! sjl, turn off melting of falling cloud ice, snow and graupel
    k0 = kbot
    ! sjl, turn off melting of falling cloud ice, snow and graupel
 
    ze(kbot + 1) = zs
    do k = kbot, ktop, - 1
        ze(k) = ze(k + 1) - dz(k) ! dz < 0
    enddo
 
    zt(ktop) = ze(ktop)
 
    ! -----------------------------------------------------------------------
    ! update capacity heat and latend heat coefficient
    ! -----------------------------------------------------------------------
 
    do k = k0, kbot
        lhi(k) = li00 + dc_ice * tz(k)
        icpk(k) = lhi(k) / cvm(k)
    enddo
 
    ! -----------------------------------------------------------------------
    ! melting of falling cloud ice into rain
    ! -----------------------------------------------------------------------
 
    call check_column (ktop, kbot, qi, no_fall)
 
    if (vi_fac < 1.e-5 .or. no_fall) then
        i1 = 0.
    else
 
        do k = ktop + 1, kbot
            zt(k) = ze(k) - dt5 * (vti(k - 1) + vti(k))
        enddo
        zt(kbot + 1) = zs - dtm * vti(kbot)
 
        do k = ktop, kbot
            if (zt(k + 1) >= zt(k)) zt(k + 1) = zt(k) - dz_min
        enddo
 
        if (k0 < kbot) then
            do k = kbot - 1, k0, - 1
                if (qi(k) > qrmin) then
                    do m = k + 1, kbot
                        if (zt(k + 1) >= ze(m)) exit
                        if (zt(k) < ze(m + 1) .and. tz(m) > tice) then
                            dtime = min(1.0, (ze(m) - ze(m + 1)) / (max(vr_min, vti(k)) * tau_imlt))
                            sink = min(qi(k) * dp(k) / dp(m), dtime * (tz(m) - tice) / icpk(m))
                            tmp = min(sink, dim(ql_mlt, ql(m)))
                            ql(m) = ql(m) + tmp
                            qr(m) = qr(m) - tmp + sink
                            tz(m) = tz(m) - sink * icpk(m)
                            qi(k) = qi(k) - sink * dp(m) / dp(k)
                        endif
                    enddo
                endif
            enddo
        endif
 
        if (do_sedi_w) then
            do k = ktop, kbot
                dm(k) = dp(k) * (1. + qv(k) + ql(k) + qr(k) + qi(k) + qs(k) + qg(k))
            enddo
        endif
 
        if (use_ppm) then
            call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qi, i1, m1_sol, mono_prof)
        else
            call implicit_fall (dtm, ktop, kbot, ze, vti, dp, qi, i1, m1_sol)
        endif
 
        if (do_sedi_w) then
            w1(ktop) = (dm(ktop) * w1(ktop) + m1_sol(ktop) * vti(ktop)) / (dm(ktop) - m1_sol(ktop))
            do k = ktop + 1, kbot
                w1(k) = (dm(k) * w1(k) - m1_sol(k - 1) * vti(k - 1) + m1_sol(k) * vti(k)) &
                     / (dm(k) + m1_sol(k - 1) - m1_sol(k))
            enddo
        endif
 
    endif
 
    ! -----------------------------------------------------------------------
    ! melting of falling snow into rain
    ! -----------------------------------------------------------------------
 
    r1 = 0.
 
    call check_column (ktop, kbot, qs, no_fall)
 
    if (no_fall) then
        s1 = 0.
    else
 
        do k = ktop + 1, kbot
            zt(k) = ze(k) - dt5 * (vts(k - 1) + vts(k))
        enddo
        zt(kbot + 1) = zs - dtm * vts(kbot)
 
        do k = ktop, kbot
            if (zt(k + 1) >= zt(k)) zt(k + 1) = zt(k) - dz_min
        enddo
 
        if (k0 < kbot) then
            do k = kbot - 1, k0, - 1
                if (qs(k) > qrmin) then
                    do m = k + 1, kbot
                        if (zt(k + 1) >= ze(m)) exit
                        dtime = min(dtm, (ze(m) - ze(m + 1)) / (vr_min + vts(k)))
                        if (zt(k) < ze(m + 1) .and. tz(m) > tice) then
                            dtime = min(1.0, dtime / tau_smlt)
                            sink = min(qs(k) * dp(k) / dp(m), dtime * (tz(m) - tice) / icpk(m))
                            tz(m) = tz(m) - sink * icpk(m)
                            qs(k) = qs(k) - sink * dp(m) / dp(k)
                            if (zt(k) < zs) then
                                r1 = r1 + sink * dp(m) ! precip as rain
                            else
                                ! qr source here will fall next time step (therefore, can evap)
                                qr(m) = qr(m) + sink
                            endif
                        endif
                        if (qs(k) < qrmin) exit
                    enddo
                endif
            enddo
        endif
 
        if (do_sedi_w) then
            do k = ktop, kbot
                dm(k) = dp(k) * (1. + qv(k) + ql(k) + qr(k) + qi(k) + qs(k) + qg(k))
            enddo
        endif
 
        if (use_ppm) then
            call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qs, s1, m1, mono_prof)
        else
            call implicit_fall (dtm, ktop, kbot, ze, vts, dp, qs, s1, m1)
        endif
 
        do k = ktop, kbot
            m1_sol(k) = m1_sol(k) + m1(k)
        enddo
 
        if (do_sedi_w) then
            w1(ktop) = (dm(ktop) * w1(ktop) + m1(ktop) * vts(ktop)) / (dm(ktop) - m1(ktop))
            do k = ktop + 1, kbot
                w1(k) = (dm(k) * w1(k) - m1(k - 1) * vts(k - 1) + m1(k) * vts(k)) &
                     / (dm(k) + m1(k - 1) - m1(k))
            enddo
        endif
 
    endif
 
    ! ----------------------------------------------
    ! melting of falling graupel into rain
    ! ----------------------------------------------
 
    call check_column (ktop, kbot, qg, no_fall)
 
    if (no_fall) then
        g1 = 0.
    else
 
        do k = ktop + 1, kbot
            zt(k) = ze(k) - dt5 * (vtg(k - 1) + vtg(k))
        enddo
        zt(kbot + 1) = zs - dtm * vtg(kbot)
 
        do k = ktop, kbot
            if (zt(k + 1) >= zt(k)) zt(k + 1) = zt(k) - dz_min
        enddo
 
        if (k0 < kbot) then
            do k = kbot - 1, k0, - 1
                if (qg(k) > qrmin) then
                    do m = k + 1, kbot
                        if (zt(k + 1) >= ze(m)) exit
                        dtime = min(dtm, (ze(m) - ze(m + 1)) / vtg(k))
                        if (zt(k) < ze(m + 1) .and. tz(m) > tice) then
                            dtime = min(1., dtime / tau_g2r)
                            sink = min(qg(k) * dp(k) / dp(m), dtime * (tz(m) - tice) / icpk(m))
                            tz(m) = tz(m) - sink * icpk(m)
                            qg(k) = qg(k) - sink * dp(m) / dp(k)
                            if (zt(k) < zs) then
                                r1 = r1 + sink * dp(m)
                            else
                                qr(m) = qr(m) + sink
                            endif
                        endif
                        if (qg(k) < qrmin) exit
                    enddo
                endif
            enddo
        endif
 
        if (do_sedi_w) then
            do k = ktop, kbot
                dm(k) = dp(k) * (1. + qv(k) + ql(k) + qr(k) + qi(k) + qs(k) + qg(k))
            enddo
        endif
 
        if (use_ppm) then
            call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qg, g1, m1, mono_prof)
        else
            call implicit_fall (dtm, ktop, kbot, ze, vtg, dp, qg, g1, m1)
        endif
 
        do k = ktop, kbot
            m1_sol(k) = m1_sol(k) + m1(k)
        enddo
 
        if (do_sedi_w) then
            w1(ktop) = (dm(ktop) * w1(ktop) + m1(ktop) * vtg(ktop)) / (dm(ktop) - m1(ktop))
            do k = ktop + 1, kbot
                w1(k) = (dm(k) * w1(k) - m1(k - 1) * vtg(k - 1) + m1(k) * vtg(k)) &
                     / (dm(k) + m1(k - 1) - m1(k))
            enddo
        endif
 
    endif
 
end subroutine terminal_fall
 
! =======================================================================
subroutine check_column (ktop, kbot, q, no_fall)
 
    implicit none
 
    integer, intent (in) :: ktop, kbot
 
    real, intent (in) :: q (ktop:kbot)
 
    logical, intent (out) :: no_fall
 
    integer :: k
 
    no_fall = .true.
 
    do k = ktop, kbot
        if (q(k) > qrmin) then
            no_fall = .false.
            exit
        endif
    enddo
 
end subroutine check_column
 
! =======================================================================
subroutine implicit_fall (dt, ktop, kbot, ze, vt, dp, q, precip, m1)
 
    implicit none
 
    integer, intent (in) :: ktop, kbot
 
    real, intent (in) :: dt
 
    real, intent (in), dimension (ktop:kbot + 1) :: ze
 
    real, intent (in), dimension (ktop:kbot) :: vt, dp
 
    real, intent (inout), dimension (ktop:kbot) :: q
 
    real, intent (out), dimension (ktop:kbot) :: m1
 
    real, intent (out) :: precip
 
    real, dimension (ktop:kbot) :: dz, qm, dd
 
    integer :: k
 
    do k = ktop, kbot
        dz(k) = ze(k) - ze(k + 1)
        dd(k) = dt * vt(k)
        q(k) = q(k) * dp(k)
    enddo
 
    ! -----------------------------------------------------------------------
    ! sedimentation: non - vectorizable loop
    ! -----------------------------------------------------------------------
 
    qm(ktop) = q(ktop) / (dz(ktop) + dd(ktop))
    do k = ktop + 1, kbot
        qm(k) = (q(k) + dd(k - 1) * qm(k - 1)) / (dz(k) + dd(k))
    enddo
 
    ! -----------------------------------------------------------------------
    ! qm is density at this stage
    ! -----------------------------------------------------------------------
 
    do k = ktop, kbot
        qm(k) = qm(k) * dz(k)
    enddo
 
    ! -----------------------------------------------------------------------
    ! output mass fluxes: non - vectorizable loop
    ! -----------------------------------------------------------------------
 
    m1(ktop) = q(ktop) - qm(ktop)
    do k = ktop + 1, kbot
        m1(k) = m1(k - 1) + q(k) - qm(k)
    enddo
    precip = m1(kbot)
 
    ! -----------------------------------------------------------------------
    ! update:
    ! -----------------------------------------------------------------------
 
    do k = ktop, kbot
        q(k) = qm(k) / dp(k)
    enddo
 
end subroutine implicit_fall
 
! =======================================================================
subroutine lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, q, precip, m1, mono)
 
    implicit none
 
    integer, intent (in) :: ktop, kbot
 
    real, intent (in) :: zs
 
    logical, intent (in) :: mono
 
    real, intent (in), dimension (ktop:kbot + 1) :: ze, zt
 
    real, intent (in), dimension (ktop:kbot) :: dp
 
    ! m1: flux
    real, intent (inout), dimension (ktop:kbot) :: q, m1
 
    real, intent (out) :: precip
 
    real, dimension (ktop:kbot) :: qm, dz
 
    real :: a4 (4, ktop:kbot)
 
    real :: pl, pr, delz, esl
 
    integer :: k, k0, n, m
 
    real, parameter :: r3 = 1. / 3., r23 = 2. / 3.
 
    ! -----------------------------------------------------------------------
    ! density:
    ! -----------------------------------------------------------------------
 
    do k = ktop, kbot
        dz(k) = zt(k) - zt(k + 1) ! note: dz is positive
        q(k) = q(k) * dp(k)
        a4(1, k) = q(k) / dz(k)
        qm(k) = 0.
    enddo
 
    ! -----------------------------------------------------------------------
    ! construct vertical profile with zt as coordinate
    ! -----------------------------------------------------------------------
 
    call cs_profile (a4(1, ktop), dz(ktop), kbot - ktop + 1, mono)
 
    k0 = ktop
    do k = ktop, kbot
        do n = k0, kbot
            if (ze(k) <= zt(n) .and. ze(k) >= zt(n + 1)) then
                pl = (zt(n) - ze(k)) / dz(n)
                if (zt(n + 1) <= ze(k + 1)) then
                    ! entire new grid is within the original grid
                    pr = (zt(n) - ze(k + 1)) / dz(n)
                    qm(k) = a4(2, n) + 0.5 * (a4(4, n) + a4(3, n) - a4(2, n)) * (pr + pl) - &
                        a4(4, n) * r3 * (pr * (pr + pl) + pl ** 2)
                    qm(k) = qm(k) * (ze(k) - ze(k + 1))
                    k0 = n
                    goto 555
                else
                    qm(k) = (ze(k) - zt(n + 1)) * (a4(2, n) + 0.5 * (a4(4, n) + &
                        a4(3, n) - a4(2, n)) * (1. + pl) - a4(4, n) * (r3 * (1. + pl * (1. + pl))))
                    if (n < kbot) then
                        do m = n + 1, kbot
                            ! locate the bottom edge: ze (k + 1)
                            if (ze(k + 1) < zt(m + 1)) then
                                qm(k) = qm(k) + q(m)
                            else
                                delz = zt(m) - ze(k + 1)
                                esl = delz / dz(m)
                                qm(k) = qm(k) + delz * (a4(2, m) + 0.5 * esl * &
                                     (a4(3, m) - a4(2, m) + a4(4, m) * (1. - r23 * esl)))
                                k0 = m
                                goto 555
                            endif
                        enddo
                    endif
                    goto 555
                endif
            endif
        enddo
        555 continue
    enddo
 
    m1(ktop) = q(ktop) - qm(ktop)
    do k = ktop + 1, kbot
        m1(k) = m1(k - 1) + q(k) - qm(k)
    enddo
    precip = m1(kbot)
 
    ! convert back to * dry * mixing ratio:
    ! dp must be dry air_mass (because moist air mass will be changed due to terminal fall) .
 
    do k = ktop, kbot
        q(k) = qm(k) / dp(k)
    enddo
 
end subroutine lagrangian_fall_ppm
 
subroutine cs_profile (a4, del, km, do_mono)
 
    implicit none
 
    integer, intent (in) :: km ! vertical dimension
 
    real, intent (in) :: del (km)
 
    logical, intent (in) :: do_mono
 
    real, intent (inout) :: a4 (4, km)
 
    real, parameter :: qp_min = 1.e-6
 
    real :: gam (km)
    real :: q (km + 1)
    real :: d4, bet, a_bot, grat, pmp, lac
    real :: pmp_1, lac_1, pmp_2, lac_2
    real :: da1, da2, a6da
 
    integer :: k
 
    logical extm (km)
 
    grat = del(2) / del(1) ! grid ratio
    bet = grat * (grat + 0.5)
    q(1) = (2. * grat * (grat + 1.) * a4(1, 1) + a4(1, 2)) / bet
    gam(1) = (1. + grat * (grat + 1.5)) / bet
 
    do k = 2, km
        d4 = del(k - 1) / del(k)
        bet = 2. + 2. * d4 - gam(k - 1)
        q(k) = (3. * (a4(1, k - 1) + d4 * a4(1, k)) - q(k - 1)) / bet
        gam(k) = d4 / bet
    enddo
 
    a_bot = 1. + d4 * (d4 + 1.5)
    q(km + 1) = (2. * d4 * (d4 + 1.) * a4(1, km) + a4(1, km - 1) - a_bot * q(km)) &
         / (d4 * (d4 + 0.5) - a_bot * gam(km))
 
    do k = km, 1, - 1
        q(k) = q(k) - gam(k) * q(k + 1)
    enddo
 
    ! -----------------------------------------------------------------------
    ! apply constraints
    ! -----------------------------------------------------------------------
 
    do k = 2, km
        gam(k) = a4(1, k) - a4(1, k - 1)
    enddo
 
    ! -----------------------------------------------------------------------
    ! apply large - scale constraints to all fields if not local max / min
    ! -----------------------------------------------------------------------
 
    ! -----------------------------------------------------------------------
    ! top:
    ! -----------------------------------------------------------------------
 
    q(1) = max(q(1), 0.)
    q(2) = min(q(2), max(a4(1, 1), a4(1, 2)))
    q(2) = max(q(2), min(a4(1, 1), a4(1, 2)), 0.)
 
    ! -----------------------------------------------------------------------
    ! interior:
    ! -----------------------------------------------------------------------
 
    do k = 3, km - 1
        if (gam(k - 1) * gam(k + 1) > 0.) then
            q(k) = min(q(k), max(a4(1, k - 1), a4(1, k)))
            q(k) = max(q(k), min(a4(1, k - 1), a4(1, k)))
        else
            if (gam(k - 1) > 0.) then
                ! there exists a local max
                q(k) = max(q(k), min(a4(1, k - 1), a4(1, k)))
            else
                ! there exists a local min
                q(k) = min(q(k), max(a4(1, k - 1), a4(1, k)))
                q(k) = max(q(k), 0.0)
            endif
        endif
    enddo
 
    ! -----------------------------------------------------------------------
    ! bottom :
    ! -----------------------------------------------------------------------
 
    q(km) = min(q(km), max(a4(1, km - 1), a4(1, km)))
    q(km) = max(q(km), min(a4(1, km - 1), a4(1, km)), 0.)
    ! q (km + 1) = max (q (km + 1), 0.)
 
    ! -----------------------------------------------------------------------
    ! f (s) = al + s * [ (ar - al) + a6 * (1 - s) ] (0 <= s <= 1)
    ! -----------------------------------------------------------------------
 
    do k = 1, km - 1
        a4(2, k) = q(k)
        a4(3, k) = q(k + 1)
    enddo
 
    do k = 2, km - 1
        if (gam(k) * gam(k + 1) > 0.0) then
            extm(k) = .false.
        else
            extm(k) = .true.
        endif
    enddo
 
    if (do_mono) then
        do k = 3, km - 2
            if (extm(k)) then
                ! positive definite constraint only if true local extrema
                if (a4(1, k) < qp_min .or. extm(k - 1) .or. extm(k + 1)) then
                    a4(2, k) = a4(1, k)
                    a4(3, k) = a4(1, k)
                endif
            else
                a4(4, k) = 6. * a4(1, k) - 3. * (a4(2, k) + a4(3, k))
                if (abs(a4(4, k)) > abs(a4(2, k) - a4(3, k))) then
                    ! check within the smooth region if subgrid profile is non - monotonic
                    pmp_1 = a4(1, k) - 2.0 * gam(k + 1)
                    lac_1 = pmp_1 + 1.5 * gam(k + 2)
                    a4(2, k) = min(max(a4(2, k), min(a4(1, k), pmp_1, lac_1)), &
                        max(a4(1, k), pmp_1, lac_1))
                    pmp_2 = a4(1, k) + 2.0 * gam(k)
                    lac_2 = pmp_2 - 1.5 * gam(k - 1)
                    a4(3, k) = min(max(a4(3, k), min(a4(1, k), pmp_2, lac_2)), &
                        max(a4(1, k), pmp_2, lac_2))
                endif
            endif
        enddo
    else
        do k = 3, km - 2
            if (extm(k)) then
                if (a4(1, k) < qp_min .or. extm(k - 1) .or. extm(k + 1)) then
                    a4(2, k) = a4(1, k)
                    a4(3, k) = a4(1, k)
                endif
            endif
        enddo
    endif
 
    do k = 1, km - 1
        a4(4, k) = 6. * a4(1, k) - 3. * (a4(2, k) + a4(3, k))
    enddo
 
    k = km - 1
    if (extm(k)) then
        a4(2, k) = a4(1, k)
        a4(3, k) = a4(1, k)
        a4(4, k) = 0.
    else
        da1 = a4(3, k) - a4(2, k)
        da2 = da1 ** 2
        a6da = a4(4, k) * da1
        if (a6da < - da2) then
            a4(4, k) = 3. * (a4(2, k) - a4(1, k))
            a4(3, k) = a4(2, k) - a4(4, k)
        elseif (a6da > da2) then
            a4(4, k) = 3. * (a4(3, k) - a4(1, k))
            a4(2, k) = a4(3, k) - a4(4, k)
        endif
    endif
 
    call cs_limiters (km - 1, a4)
 
    ! -----------------------------------------------------------------------
    ! bottom layer:
    ! -----------------------------------------------------------------------
 
    a4(2, km) = a4(1, km)
    a4(3, km) = a4(1, km)
    a4(4, km) = 0.
 
end subroutine cs_profile
 
subroutine cs_limiters (km, a4)
 
    implicit none
 
    integer, intent (in) :: km
 
    real, intent (inout) :: a4 (4, km) ! ppm array
 
    real, parameter :: r12 = 1. / 12.
 
    integer :: k
 
    ! -----------------------------------------------------------------------
    ! positive definite constraint
    ! -----------------------------------------------------------------------
 
    do k = 1, km
        if (abs(a4(3, k) - a4(2, k)) < - a4(4, k)) then
            if ((a4(1, k) + 0.25 * (a4(3, k) - a4(2, k)) ** 2 / a4(4, k) + a4(4, k) * r12) < 0.) then
                if (a4(1, k) < a4(3, k) .and. a4(1, k) < a4(2, k)) then
                    a4(3, k) = a4(1, k)
                    a4(2, k) = a4(1, k)
                    a4(4, k) = 0.
                elseif (a4(3, k) > a4(2, k)) then
                    a4(4, k) = 3. * (a4(2, k) - a4(1, k))
                    a4(3, k) = a4(2, k) - a4(4, k)
                else
                    a4(4, k) = 3. * (a4(3, k) - a4(1, k))
                    a4(2, k) = a4(3, k) - a4(4, k)
                endif
            endif
        endif
    enddo
 
end subroutine cs_limiters
 
! =======================================================================
subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg)
 
    implicit none
 
    integer, intent (in) :: ktop, kbot
 
    real, intent (in), dimension (ktop:kbot) :: den, qs, qi, qg, ql, tk
    real, intent (out), dimension (ktop:kbot) :: vts, vti, vtg
 
    ! fall velocity constants:
 
    real, parameter :: thi = 1.0e-8 
    real, parameter :: thg = 1.0e-8
    real, parameter :: ths = 1.0e-8
 
    ! coefficient for the parameterization of mass weighted fall velocity
    ! as a function of temperature and IWC.
    ! Table 1 in Deng and Mace (2008) \cite deng_and_mace_2008.
    real, parameter :: aa = - 4.14122e-5
    real, parameter :: bb = - 0.00538922
    real, parameter :: cc = - 0.0516344
    real, parameter :: dd = 0.00216078
    real, parameter :: ee = 1.9714
 
    ! marshall - palmer constants
 
    real, parameter :: vcons = 6.6280504
    real, parameter :: vcong = 87.2382675
    real, parameter :: vconh = vcong * sqrt(rhoh / rhog)
    real, parameter :: norms = 942477796.076938
    real, parameter :: normg = 5026548245.74367
    real, parameter :: normh = pi * rhoh * rnzh
 
    real, dimension (ktop:kbot) :: qden, tc, rhof
 
    real :: vi0
 
    integer :: k
 
    ! -----------------------------------------------------------------------
    ! marshall - palmer formula
    ! -----------------------------------------------------------------------
 
    ! -----------------------------------------------------------------------
    ! try the local air density -- for global model; the true value could be
    ! much smaller than sfcrho over high mountains
    ! -----------------------------------------------------------------------
 
    do k = ktop, kbot
        rhof(k) = sqrt(min(10., sfcrho / den(k)))
    enddo
 
    ! -----------------------------------------------------------------------
    ! -----------------------------------------------------------------------
 
    if (const_vi) then
        vti(:) = vi_fac
    else
        ! -----------------------------------------------------------------------
        ! use deng and mace (2008, grl), which gives smaller fall speed than hd90 formula
        ! -----------------------------------------------------------------------
        vi0 = 0.01 * vi_fac
        do k = ktop, kbot
            if (qi(k) < thi) then ! this is needed as the fall - speed maybe problematic for small qi
                vti(k) = vf_min
            else
                tc(k) = tk(k) - tice
                vti(k) = (3. + log10(qi(k) * den(k))) * (tc(k) * (aa * tc(k) + bb) + cc) + dd * tc(k) + ee
                vti(k) = vi0 * exp(log_10 * vti(k)) * 0.9
                vti(k) = min(vi_max, max(vf_min, vti(k)))
            endif
        enddo
    endif
 
    ! -----------------------------------------------------------------------
    ! -----------------------------------------------------------------------
 
    if (const_vs) then
        vts(:) = vs_fac ! 1. ifs_2016
    else
        do k = ktop, kbot
            if (qs(k) < ths) then
                vts(k) = vf_min
            else
                vts(k) = vs_fac * vcons * rhof(k) * exp(0.0625 * log(qs(k) * den(k) / norms))
                vts(k) = min(vs_max, max(vf_min, vts(k)))
            endif
        enddo
    endif
 
    ! -----------------------------------------------------------------------
    ! -----------------------------------------------------------------------
    if (const_vg) then
        vtg(:) = vg_fac ! 2.
    else
        if (do_hail) then
            do k = ktop, kbot
                if (qg(k) < thg) then
                    vtg(k) = vf_min
                else
                    vtg(k) = vg_fac * vconh * rhof(k) * sqrt(sqrt(sqrt(qg(k) * den(k) / normh)))
                    vtg(k) = min(vg_max, max(vf_min, vtg(k)))
                endif
            enddo
        else
        do k = ktop, kbot
            if (qg(k) < thg) then
                vtg(k) = vf_min
            else
                vtg(k) = vg_fac * vcong * rhof(k) * sqrt(sqrt(sqrt(qg(k) * den(k) / normg)))
                vtg(k) = min(vg_max, max(vf_min, vtg(k)))
            endif
        enddo
    endif
    endif
 
end subroutine fall_speed
 
! =======================================================================
subroutine setupm
 
    implicit none
 
    real :: gcon, cd, scm3, pisq, act (8)
    real :: vdifu, tcond
    real :: visk
    real :: ch2o, hltf
    real :: hlts, hltc, ri50
 
    real, parameter :: gam263 = 1.456943, gam275 = 1.608355, gam290 = 1.827363, &
        gam325 = 2.54925, gam350 = 3.323363, gam380 = 4.694155, &
        gam425 = 8.285063, gam450 = 11.631769, gam480 = 17.837789, &
        gam625 = 184.860962, gam680 = 496.604067
 
    ! intercept parameters
 
!    real, parameter :: rnzr = 8.0e6 ! lin83
!    real, parameter :: rnzs = 3.0e6 ! lin83
!    real, parameter :: rnzg = 4.0e6 ! rh84
 
    ! density parameters
 
!    real, parameter :: rhos = 0.1e3 !< lin83 (snow density; 1 / 10 of water)
!    real, parameter :: rhog = 0.4e3 !< rh84 (graupel density)
    real, parameter :: acc (3) = (/ 5.0, 2.0, 0.5 /)
 
    real den_rc
 
    integer :: i, k
 
    pie = 4. * atan(1.0)
 
    ! s. klein's formular (eq 16) from am2
 
    fac_rc = (4. / 3.) * pie * rhor * rthresh ** 3
 
    if (prog_ccn) then
        ! if (master) write (*, *) 'prog_ccn option is .t.'
    else
        den_rc = fac_rc * ccn_o * 1.e6
        ! if (master) write (*, *) 'mp: for ccn_o = ', ccn_o, 'ql_rc = ', den_rc
        den_rc = fac_rc * ccn_l * 1.e6
        ! if (master) write (*, *) 'mp: for ccn_l = ', ccn_l, 'ql_rc = ', den_rc
    endif
 
    vdifu = 2.11e-5
    tcond = 2.36e-2
 
    visk = 1.259e-5
    hlts = 2.8336e6
    hltc = 2.5e6
    hltf = 3.336e5
 
    ch2o = 4.1855e3
    ri50 = 1.e-4
 
    pisq = pie * pie
    scm3 = (visk / vdifu) ** (1. / 3.)
 
    cracs = pisq * rnzr * rnzs * rhos
    csacr = pisq * rnzr * rnzs * rhor
    if (do_hail) then
        cgacr = pisq * rnzr * rnzh * rhor
        cgacs = pisq * rnzh * rnzs * rhos
    else
        cgacr = pisq * rnzr * rnzg * rhor
        cgacs = pisq * rnzg * rnzs * rhos
    endif
    cgacs = cgacs * c_pgacs
 
    ! act: 1 - 2:racs (s - r) ; 3 - 4:sacr (r - s) ;
    ! 5 - 6:gacr (r - g) ; 7 - 8:gacs (s - g)
 
    act(1) = pie * rnzs * rhos
    act(2) = pie * rnzr * rhor
    if (do_hail) then
        act(6) = pie * rnzh * rhoh
    else
        act(6) = pie * rnzg * rhog
    endif
    act(3) = act(2)
    act(4) = act(1)
    act(5) = act(2)
    act(7) = act(1)
    act(8) = act(6)
 
    do i = 1, 3
        do k = 1, 4
            acco(i, k) = acc(i) / (act(2 * k - 1) ** ((7 - i) * 0.25) * act(2 * k) ** (i * 0.25))
        enddo
    enddo
 
    gcon = 40.74 * sqrt(sfcrho) ! 44.628
 
    csacw = pie * rnzs * clin * gam325 / (4. * act(1) ** 0.8125)
    ! decreasing csacw to reduce cloud water --- > snow
 
    craci = pie * rnzr * alin * gam380 / (4. * act(2) ** 0.95)
    csaci = csacw * c_psaci
 
    if (do_hail) then
        cgacw = pie * rnzh * gam350 * gcon / (4. * act(6) ** 0.875)
    else
        cgacw = pie * rnzg * gam350 * gcon / (4. * act(6) ** 0.875)
    endif
    ! cgaci = cgacw * 0.1
 
    ! sjl, may 28, 2012
    cgaci = cgacw * 0.05
    ! sjl, may 28, 2012
 
    cracw = craci ! cracw = 3.27206196043822
    cracw = c_cracw * cracw
 
    ! subl and revp: five constants for three separate processes
 
    cssub(1) = 2. * pie * vdifu * tcond * rvgas * rnzs
    if (do_hail) then
        cgsub(1) = 2. * pie * vdifu * tcond * rvgas * rnzh
    else
        cgsub(1) = 2. * pie * vdifu * tcond * rvgas * rnzg
    endif
    crevp(1) = 2. * pie * vdifu * tcond * rvgas * rnzr
    cssub(2) = 0.78 / sqrt(act(1))
    cgsub(2) = 0.78 / sqrt(act(6))
    crevp(2) = 0.78 / sqrt(act(2))
    cssub(3) = 0.31 * scm3 * gam263 * sqrt(clin / visk) / act(1) ** 0.65625
    cgsub(3) = 0.31 * scm3 * gam275 * sqrt(gcon / visk) / act(6) ** 0.6875
    crevp(3) = 0.31 * scm3 * gam290 * sqrt(alin / visk) / act(2) ** 0.725
    cssub(4) = tcond * rvgas
    cssub(5) = hlts ** 2 * vdifu
    cgsub(4) = cssub(4)
    crevp(4) = cssub(4)
    cgsub(5) = cssub(5)
    crevp(5) = hltc ** 2 * vdifu
 
    cgfr(1) = 20.e2 * pisq * rnzr * rhor / act(2) ** 1.75
    cgfr(2) = 0.66
 
    ! smlt: five constants (lin et al. 1983)
 
    csmlt(1) = 2. * pie * tcond * rnzs / hltf
    csmlt(2) = 2. * pie * vdifu * rnzs * hltc / hltf
    csmlt(3) = cssub(2)
    csmlt(4) = cssub(3)
    csmlt(5) = ch2o / hltf
 
    ! gmlt: five constants
 
    if (do_hail) then
        cgmlt(1) = 2. * pie * tcond * rnzh / hltf
        cgmlt(2) = 2. * pie * vdifu * rnzh * hltc / hltf
    else
        cgmlt(1) = 2. * pie * tcond * rnzg / hltf
        cgmlt(2) = 2. * pie * vdifu * rnzg * hltc / hltf
    endif
    cgmlt(3) = cgsub(2)
    cgmlt(4) = cgsub(3)
    cgmlt(5) = ch2o / hltf
 
    es0 = 6.107799961e2 ! ~6.1 mb
    ces0 = eps * es0
 
end subroutine setupm
 
! =======================================================================
! initialization of gfdl cloud microphysics
subroutine gfdl_cloud_microphys_mod_init (me, master, nlunit, input_nml_file, logunit, &
     fn_nml, errmsg, errflg)
 
    implicit none
 
    integer, intent (in) :: me
    integer, intent (in) :: master
    integer, intent (in) :: nlunit
    integer, intent (in) :: logunit
 
    character (len = 64), intent (in) :: fn_nml
    character (len = *),  intent (in) :: input_nml_file(:)
    character(len=*),     intent(out) :: errmsg
    integer,              intent(out) :: errflg
 
    integer :: ios
    logical :: exists
 
    ! integer, intent (in) :: id, jd, kd
    ! integer, intent (in) :: axes (4)
    ! type (time_type), intent (in) :: time
 
    ! integer :: unit, io, ierr, k, logunit
    ! logical :: flag
    ! real :: tmp, q1, q2
 
    ! master = (mpp_pe () .eq.mpp_root_pe ())
 
    ! Initialize CCPP error-handling
    errflg = 0
    errmsg = ''
 
    ! -----------------------------------------------------------------------
    ! Read namelist
    ! -----------------------------------------------------------------------
    call read_gfdlmp_nml(errmsg = errmsg, errflg = errflg, unit = nlunit,   &
         input_nml_file = input_nml_file, fn_nml = fn_nml, version=1,       &
         iostat = ios)
    if(errflg/=0) return
 
    ! write version number and namelist to log file
    if (me == master) then
        write (logunit, *) " ================================================================== "
        write (logunit, *) "gfdl_cloud_microphysics_nml"
    endif
 
    !
    if (do_setup) then
        call setup_con
        call setupm
        do_setup = .false.
    endif
 
    log_10 = log(10.)
 
    tice0 = tice - 0.01
    t_wfr = tice - 40.0 ! supercooled water can exist down to - 48 c, which is the "absolute"
 
    ! if (master) write (logunit, nml = gfdl_cloud_microphys_nml)
    !
    ! id_vtr = register_diag_field (mod_name, 'vt_r', axes (1:3), time, &
    ! 'rain fall speed', 'm / s', missing_value = missing_value)
    ! id_vts = register_diag_field (mod_name, 'vt_s', axes (1:3), time, &
    ! 'snow fall speed', 'm / s', missing_value = missing_value)
    ! id_vtg = register_diag_field (mod_name, 'vt_g', axes (1:3), time, &
    ! 'graupel fall speed', 'm / s', missing_value = missing_value)
    ! id_vti = register_diag_field (mod_name, 'vt_i', axes (1:3), time, &
    ! 'ice fall speed', 'm / s', missing_value = missing_value)
 
    ! id_droplets = register_diag_field (mod_name, 'droplets', axes (1:3), time, &
    ! 'droplet number concentration', '# / m3', missing_value = missing_value)
    ! id_rh = register_diag_field (mod_name, 'rh_lin', axes (1:2), time, &
    ! 'relative humidity', 'n / a', missing_value = missing_value)
 
    ! id_rain = register_diag_field (mod_name, 'rain_lin', axes (1:2), time, &
    ! 'rain_lin', 'mm / day', missing_value = missing_value)
    ! id_snow = register_diag_field (mod_name, 'snow_lin', axes (1:2), time, &
    ! 'snow_lin', 'mm / day', missing_value = missing_value)
    ! id_graupel = register_diag_field (mod_name, 'graupel_lin', axes (1:2), time, &
    ! 'graupel_lin', 'mm / day', missing_value = missing_value)
    ! id_ice = register_diag_field (mod_name, 'ice_lin', axes (1:2), time, &
    ! 'ice_lin', 'mm / day', missing_value = missing_value)
    ! id_prec = register_diag_field (mod_name, 'prec_lin', axes (1:2), time, &
    ! 'prec_lin', 'mm / day', missing_value = missing_value)
 
    ! if (master) write (*, *) 'prec_lin diagnostics initialized.', id_prec
 
    ! id_cond = register_diag_field (mod_name, 'cond_lin', axes (1:2), time, &
    ! 'total condensate', 'kg / m ** 2', missing_value = missing_value)
    ! id_var = register_diag_field (mod_name, 'var_lin', axes (1:2), time, &
    ! 'subgrid variance', 'n / a', missing_value = missing_value)
 
    ! call qsmith_init
 
    ! testing the water vapor tables
 
    ! if (mp_debug .and. master) then
    ! write (*, *) 'testing water vapor tables in gfdl_cloud_microphys'
    ! tmp = tice - 90.
    ! do k = 1, 25
    ! q1 = wqsat_moist (tmp, 0., 1.e5)
    ! q2 = qs1d_m (tmp, 0., 1.e5)
    ! write (*, *) nint (tmp - tice), q1, q2, 'dq = ', q1 - q2
    ! tmp = tmp + 5.
    ! enddo
    ! endif
 
    ! if (master) write (*, *) 'gfdl_cloud_micrphys diagnostics initialized.'
 
    module_is_initialized = .true.
 
!+---+-----------------------------------------------------------------+
!..Set these variables needed for computing radar reflectivity.  These
!.. get used within radar_init to create other variables used in the
!.. radar module.
 
   xam_r = pi*rhor/6.
   xbm_r = 3.
   xmu_r = 0.
   xam_s = pi*rhos/6.
   xbm_s = 3.
   xmu_s = 0.
   xam_g = pi*rhog/6.
   xbm_g = 3.
   xmu_g = 0.
 
   call radar_init
 
end subroutine gfdl_cloud_microphys_mod_init
 
! =======================================================================
! end of gfdl cloud microphysics
subroutine gfdl_cloud_microphys_mod_end()
 
    implicit none
 
    deallocate (table)
    deallocate (table2)
    deallocate (table3)
    deallocate (tablew)
    deallocate (des)
    deallocate (des2)
    deallocate (des3)
    deallocate (desw)
 
    tables_are_initialized = .false.
 
end subroutine gfdl_cloud_microphys_mod_end
 
! =======================================================================
! qsmith table initialization
subroutine setup_con
 
    implicit none
 
    ! master = (mpp_pe () .eq.mpp_root_pe ())
 
    rgrav = 1. / grav
 
    if (.not. qsmith_tables_initialized) call qsmith_init
 
    qsmith_tables_initialized = .true.
 
end subroutine setup_con
 
! =======================================================================
! =======================================================================
 
real function acr3d (v1, v2, q1, q2, c, cac, rho)
 
    implicit none
 
    real, intent (in) :: v1, v2, c, rho
    real, intent (in) :: q1, q2 ! mixing ratio!!!
    real, intent (in) :: cac (3)
 
    real :: t1, s1, s2
 
    ! integer :: k
    !
    ! real :: a
    !
    ! a = 0.0
    ! do k = 1, 3
    ! a = a + cac (k) * ((q1 * rho) ** ((7 - k) * 0.25) * (q2 * rho) ** (k * 0.25))
    ! enddo
    ! acr3d = c * abs (v1 - v2) * a / rho
 
    ! optimized
 
    t1 = sqrt(q1 * rho)
    s1 = sqrt(q2 * rho)
    s2 = sqrt(s1) ! s1 = s2 ** 2
    acr3d = c * abs(v1 - v2) * q1 * s2 * (cac(1) * t1 + cac(2) * sqrt(t1) * s2 + cac(3) * s1)
 
end function acr3d
 
! =======================================================================
! =======================================================================
 
real function smlt (tc, dqs, qsrho, psacw, psacr, c, rho, rhofac)
 
    implicit none
 
    real, intent (in) :: tc, dqs, qsrho, psacw, psacr, c (5), rho, rhofac
 
    smlt = (c(1) * tc / rho - c(2) * dqs) * (c(3) * sqrt(qsrho) + &
        c(4) * qsrho ** 0.65625 * sqrt(rhofac)) + c(5) * tc * (psacw + psacr)
 
end function smlt
 
! =======================================================================
! =======================================================================
 
real function gmlt (tc, dqs, qgrho, pgacw, pgacr, c, rho)
 
    implicit none
 
    real, intent (in) :: tc, dqs, qgrho, pgacw, pgacr, c (5), rho
 
    gmlt = (c(1) * tc / rho - c(2) * dqs) * (c(3) * sqrt(qgrho) + &
        c(4) * qgrho ** 0.6875 / rho ** 0.25) + c(5) * tc * (pgacw + pgacr)
 
end function gmlt
 
! =======================================================================
! initialization
! prepare saturation water vapor pressure tables
! =======================================================================
! =======================================================================
subroutine qsmith_init
 
    implicit none
 
    integer, parameter :: length = 2621
 
    integer :: i
 
    if (.not. tables_are_initialized) then
 
       ! master = (mpp_pe () .eq. mpp_root_pe ())
       ! if (master) print *, ' gfdl mp: initializing qs tables'
 
       ! debug code
       ! print *, mpp_pe (), allocated (table), allocated (table2), &
       ! allocated (table3), allocated (tablew), allocated (des), &
       ! allocated (des2), allocated (des3), allocated (desw)
       ! end debug code
 
       ! generate es table (dt = 0.1 deg. c)
 
       allocate (table(length))
       allocate (table2(length))
       allocate (table3(length))
       allocate (tablew(length))
       allocate (des(length))
       allocate (des2(length))
       allocate (des3(length))
       allocate (desw(length))
 
       call qs_table (length)
       call qs_table2 (length)
       call qs_table3 (length)
       call qs_tablew (length)
 
       do i = 1, length - 1
           des(i) = max(0., table(i + 1) - table(i))
           des2(i) = max(0., table2(i + 1) - table2(i))
           des3(i) = max(0., table3(i + 1) - table3(i))
           desw(i) = max(0., tablew(i + 1) - tablew(i))
       enddo
       des(length) = des(length - 1)
       des2(length) = des2(length - 1)
       des3(length) = des3(length - 1)
       desw(length) = desw(length - 1)
 
        tables_are_initialized = .true.
 
    endif
 
end subroutine qsmith_init
 
! =======================================================================
! compute the saturated specific humidity for table ii
real function wqs1 (ta, den)
 
    implicit none
 
 
    real, intent (in) :: ta, den
 
    real :: es, ap1, tmin
 
    integer :: it
 
    tmin = table_ice - 160.
    ap1 = 10. * dim(ta, tmin) + 1.
    ap1 = min(2621., ap1)
    it = ap1
    es = tablew(it) + (ap1 - it) * desw(it)
    wqs1 = es / (rvgas * ta * den)
 
end function wqs1
 
! =======================================================================
! compute the gradient of saturated specific humidity for table ii
! =======================================================================
 
real function wqs2 (ta, den, dqdt)
 
    implicit none
 
    ! pure water phase; universal dry / moist formular using air density
    ! input "den" can be either dry or moist air density
 
    real, intent (in) :: ta, den
 
    real, intent (out) :: dqdt
 
    real :: es, ap1, tmin
 
    integer :: it
 
    tmin = table_ice - 160.
 
    if (.not. tables_are_initialized) call qsmith_init
 
    ap1 = 10. * dim(ta, tmin) + 1.
    ap1 = min(2621., ap1)
    it = ap1
    es = tablew(it) + (ap1 - it) * desw(it)
    wqs2 = es / (rvgas * ta * den)
    it = ap1 - 0.5
    ! finite diff, del_t = 0.1:
    dqdt = 10. * (desw(it) + (ap1 - it) * (desw(it + 1) - desw(it))) / (rvgas * ta * den)
 
end function wqs2
 
! =======================================================================
! compute wet buld temperature
! =======================================================================
 
real function wet_bulb (q, t, den)
 
    implicit none
 
    real, intent (in) :: t, q, den
 
    real :: qs, tp, dqdt
 
    wet_bulb = t
    qs = wqs2(wet_bulb, den, dqdt)
    tp = 0.5 * (qs - q) / (1. + lcp * dqdt) * lcp
    wet_bulb = wet_bulb - tp
 
    ! tp is negative if super - saturated
    if (tp > 0.01) then
        qs = wqs2(wet_bulb, den, dqdt)
        tp = (qs - q) / (1. + lcp * dqdt) * lcp
        wet_bulb = wet_bulb - tp
    endif
 
end function wet_bulb
 
! =======================================================================
! =======================================================================
 
real function iqs1 (ta, den)
 
    implicit none
 
 
    real, intent (in) :: ta, den
 
    real :: es, ap1, tmin
 
    integer :: it
 
    tmin = table_ice - 160.
    ap1 = 10. * dim(ta, tmin) + 1.
    ap1 = min(2621., ap1)
    it = ap1
    es = table2(it) + (ap1 - it) * des2(it)
    iqs1 = es / (rvgas * ta * den)
 
end function iqs1
 
! =======================================================================
! =======================================================================
 
real function iqs2 (ta, den, dqdt)
 
    implicit none
 
    ! water - ice phase; universal dry / moist formular using air density
    ! input "den" can be either dry or moist air density
 
    real, intent (in) :: ta, den
 
    real, intent (out) :: dqdt
 
    real :: es, ap1, tmin
 
    integer :: it
 
    tmin = table_ice - 160.
    ap1 = 10. * dim(ta, tmin) + 1.
    ap1 = min(2621., ap1)
    it = ap1
    es = table2(it) + (ap1 - it) * des2(it)
    iqs2 = es / (rvgas * ta * den)
    it = ap1 - 0.5
    dqdt = 10. * (des2(it) + (ap1 - it) * (des2(it + 1) - des2(it))) / (rvgas * ta * den)
 
end function iqs2
 
! =======================================================================
! =======================================================================
 
real function qs1d_moist (ta, qv, pa, dqdt)
 
    implicit none
 
    real, intent (in) :: ta, pa, qv
 
    real, intent (out) :: dqdt
 
    real :: es, ap1, tmin, eps10
 
    integer :: it
 
    tmin = table_ice - 160.
    eps10 = 10. * eps
    ap1 = 10. * dim(ta, tmin) + 1.
    ap1 = min(2621., ap1)
    it = ap1
    es = table2(it) + (ap1 - it) * des2(it)
    qs1d_moist = eps * es * (1. + zvir * qv) / pa
    it = ap1 - 0.5
    dqdt = eps10 * (des2(it) + (ap1 - it) * (des2(it + 1) - des2(it))) * (1. + zvir * qv) / pa
 
end function qs1d_moist
 
! =======================================================================
! compute the gradient of saturated specific humidity for table ii
! =======================================================================
 
real function wqsat2_moist (ta, qv, pa, dqdt)
 
    implicit none
 
    real, intent (in) :: ta, pa, qv
 
    real, intent (out) :: dqdt
 
    real :: es, ap1, tmin, eps10
 
    integer :: it
 
    tmin = table_ice - 160.
    eps10 = 10. * eps
    ap1 = 10. * dim(ta, tmin) + 1.
    ap1 = min(2621., ap1)
    it = ap1
    es = tablew(it) + (ap1 - it) * desw(it)
    wqsat2_moist = eps * es * (1. + zvir * qv) / pa
    it = ap1 - 0.5
    dqdt = eps10 * (desw(it) + (ap1 - it) * (desw(it + 1) - desw(it))) * (1. + zvir * qv) / pa
 
end function wqsat2_moist
 
! =======================================================================
! compute the saturated specific humidity for table ii
! =======================================================================
 
real function wqsat_moist (ta, qv, pa)
 
    implicit none
 
    real, intent (in) :: ta, pa, qv
 
    real :: es, ap1, tmin
 
    integer :: it
 
    tmin = table_ice - 160.
    ap1 = 10. * dim(ta, tmin) + 1.
    ap1 = min(2621., ap1)
    it = ap1
    es = tablew(it) + (ap1 - it) * desw(it)
    wqsat_moist = eps * es * (1. + zvir * qv) / pa
 
end function wqsat_moist
 
! =======================================================================
! =======================================================================
 
real function qs1d_m (ta, qv, pa)
 
    implicit none
 
    real, intent (in) :: ta, pa, qv
 
    real :: es, ap1, tmin
 
    integer :: it
 
    tmin = table_ice - 160.
    ap1 = 10. * dim(ta, tmin) + 1.
    ap1 = min(2621., ap1)
    it = ap1
    es = table2(it) + (ap1 - it) * des2(it)
    qs1d_m = eps * es * (1. + zvir * qv) / pa
 
end function qs1d_m
 
! =======================================================================
! =======================================================================
 
real function d_sat (ta, den)
 
    implicit none
 
    real, intent (in) :: ta, den
 
    real :: es_w, es_i, ap1, tmin
 
    integer :: it
 
    tmin = table_ice - 160.
    ap1 = 10. * dim(ta, tmin) + 1.
    ap1 = min(2621., ap1)
    it = ap1
    es_w = tablew(it) + (ap1 - it) * desw(it)
    es_i = table2(it) + (ap1 - it) * des2(it)
    d_sat = dim(es_w, es_i) / (rvgas * ta * den) ! take positive difference
 
end function d_sat
 
! =======================================================================
! =======================================================================
 
real function esw_table (ta)
 
    implicit none
 
    real, intent (in) :: ta
 
    real :: ap1, tmin
 
    integer :: it
 
    tmin = table_ice - 160.
    ap1 = 10. * dim(ta, tmin) + 1.
    ap1 = min(2621., ap1)
    it = ap1
    esw_table = tablew(it) + (ap1 - it) * desw(it)
 
end function esw_table
 
! =======================================================================
! =======================================================================
 
real function es2_table (ta)
 
    implicit none
 
    real, intent (in) :: ta
 
    real :: ap1, tmin
 
    integer :: it
 
    tmin = table_ice - 160.
    ap1 = 10. * dim(ta, tmin) + 1.
    ap1 = min(2621., ap1)
    it = ap1
    es2_table = table2(it) + (ap1 - it) * des2(it)
 
end function es2_table
 
! =======================================================================
subroutine esw_table1d (ta, es, n)
 
    implicit none
 
    integer, intent (in) :: n
 
    real, intent (in) :: ta (n)
 
    real, intent (out) :: es (n)
 
    real :: ap1, tmin
 
    integer :: i, it
 
    tmin = table_ice - 160.
 
    do i = 1, n
        ap1 = 10. * dim(ta(i), tmin) + 1.
        ap1 = min(2621., ap1)
        it = ap1
        es(i) = tablew(it) + (ap1 - it) * desw(it)
    enddo
 
end subroutine esw_table1d
 
! =======================================================================
subroutine es2_table1d (ta, es, n)
 
    implicit none
 
    integer, intent (in) :: n
 
    real, intent (in) :: ta (n)
 
    real, intent (out) :: es (n)
 
    real :: ap1, tmin
 
    integer :: i, it
 
    tmin = table_ice - 160.
 
    do i = 1, n
        ap1 = 10. * dim(ta(i), tmin) + 1.
        ap1 = min(2621., ap1)
        it = ap1
        es(i) = table2(it) + (ap1 - it) * des2(it)
    enddo
 
end subroutine es2_table1d
 
! =======================================================================
subroutine es3_table1d (ta, es, n)
 
    implicit none
 
    integer, intent (in) :: n
 
    real, intent (in) :: ta (n)
 
    real, intent (out) :: es (n)
 
    real :: ap1, tmin
 
    integer :: i, it
 
    tmin = table_ice - 160.
 
    do i = 1, n
        ap1 = 10. * dim(ta(i), tmin) + 1.
        ap1 = min(2621., ap1)
        it = ap1
        es(i) = table3(it) + (ap1 - it) * des3(it)
    enddo
 
end subroutine es3_table1d
 
! =======================================================================
! 1 - phase table
subroutine qs_tablew (n)
 
    implicit none
 
    integer, intent (in) :: n
 
    real :: delt = 0.1
    real :: tmin, tem, fac0, fac1, fac2
 
    integer :: i
 
    tmin = table_ice - 160.
 
    ! -----------------------------------------------------------------------
    ! compute es over water
    ! -----------------------------------------------------------------------
 
    do i = 1, n
        tem = tmin + delt * real(i - 1)
        fac0 = (tem - t_ice) / (tem * t_ice)
        fac1 = fac0 * lv0
        fac2 = (dc_vap * log(tem / t_ice) + fac1) / rvgas
        tablew(i) = e00 * exp(fac2)
    enddo
 
end subroutine qs_tablew
 
! =======================================================================
! 2 - phase table
subroutine qs_table2 (n)
 
    implicit none
 
    integer, intent (in) :: n
 
    real :: delt = 0.1
    real :: tmin, tem0, tem1, fac0, fac1, fac2
 
    integer :: i, i0, i1
 
    tmin = table_ice - 160.
 
    do i = 1, n
        tem0 = tmin + delt * real(i - 1)
        fac0 = (tem0 - t_ice) / (tem0 * t_ice)
        if (i <= 1600) then
            ! -----------------------------------------------------------------------
            ! compute es over ice between - 160 deg c and 0 deg c.
            ! -----------------------------------------------------------------------
            fac1 = fac0 * li2
            fac2 = (d2ice * log(tem0 / t_ice) + fac1) / rvgas
        else
            ! -----------------------------------------------------------------------
            ! compute es over water between 0 deg c and 102 deg c.
            ! -----------------------------------------------------------------------
            fac1 = fac0 * lv0
            fac2 = (dc_vap * log(tem0 / t_ice) + fac1) / rvgas
        endif
        table2(i) = e00 * exp(fac2)
    enddo
 
    ! -----------------------------------------------------------------------
    ! smoother around 0 deg c
    ! -----------------------------------------------------------------------
 
    i0 = 1600
    i1 = 1601
    tem0 = 0.25 * (table2(i0 - 1) + 2. * table(i0) + table2(i0 + 1))
    tem1 = 0.25 * (table2(i1 - 1) + 2. * table(i1) + table2(i1 + 1))
    table2(i0) = tem0
    table2(i1) = tem1
 
end subroutine qs_table2
 
! =======================================================================
! 2 - phase table with " - 2 c" as the transition point
subroutine qs_table3 (n)
 
    implicit none
 
    integer, intent (in) :: n
 
    real :: delt = 0.1
    real :: esbasw, tbasw, esbasi, tmin, tem, aa, b, c, d, e
    real :: tem0, tem1
 
    integer :: i, i0, i1
 
    esbasw = 1013246.0
    tbasw = table_ice + 100.
    esbasi = 6107.1
    tmin = table_ice - 160.
 
    do i = 1, n
        tem = tmin + delt * real(i - 1)
        ! if (i <= 1600) then
        if (i <= 1580) then ! change to - 2 c
            ! -----------------------------------------------------------------------
            ! compute es over ice between - 160 deg c and 0 deg c.
            ! see smithsonian meteorological tables page 350.
            ! -----------------------------------------------------------------------
            aa = - 9.09718 * (table_ice / tem - 1.)
            b = - 3.56654 * alog10(table_ice / tem)
            c = 0.876793 * (1. - tem / table_ice)
            e = alog10(esbasi)
            table3(i) = 0.1 * 10 ** (aa + b + c + e)
        else
            ! -----------------------------------------------------------------------
            ! compute es over water between - 2 deg c and 102 deg c.
            ! see smithsonian meteorological tables page 350.
            ! -----------------------------------------------------------------------
            aa = - 7.90298 * (tbasw / tem - 1.)
            b = 5.02808 * alog10(tbasw / tem)
            c = - 1.3816e-7 * (10 ** ((1. - tem / tbasw) * 11.344) - 1.)
            d = 8.1328e-3 * (10 ** ((tbasw / tem - 1.) * (- 3.49149)) - 1.)
            e = alog10(esbasw)
            table3(i) = 0.1 * 10 ** (aa + b + c + d + e)
        endif
    enddo
 
    ! -----------------------------------------------------------------------
    ! smoother around - 2 deg c
    ! -----------------------------------------------------------------------
 
    i0 = 1580
    i1 = 1581
    tem0 = 0.25 * (table3(i0 - 1) + 2. * table(i0) + table3(i0 + 1))
    tem1 = 0.25 * (table3(i1 - 1) + 2. * table(i1) + table3(i1 + 1))
    table3(i0) = tem0
    table3(i1) = tem1
 
end subroutine qs_table3
 
! =======================================================================
! compute the saturated specific humidity for table
! note: this routine is based on "moist" mixing ratio
real function qs_blend (t, p, q)
 
    implicit none
 
    real, intent (in) :: t, p, q
 
    real :: es, ap1, tmin
 
    integer :: it
 
    tmin = table_ice - 160.
    ap1 = 10. * dim(t, tmin) + 1.
    ap1 = min(2621., ap1)
    it = ap1
    es = table(it) + (ap1 - it) * des(it)
    qs_blend = eps * es * (1. + zvir * q) / p
 
end function qs_blend
 
! =======================================================================
! 3 - phase table
subroutine qs_table (n)
 
    implicit none
 
    integer, intent (in) :: n
 
    real :: delt = 0.1
    real :: tmin, tem, esh20
    real :: wice, wh2o, fac0, fac1, fac2
    real :: esupc (200)
 
    integer :: i
 
    tmin = table_ice - 160.
 
    ! -----------------------------------------------------------------------
    ! compute es over ice between - 160 deg c and 0 deg c.
    ! -----------------------------------------------------------------------
 
    do i = 1, 1600
        tem = tmin + delt * real(i - 1)
        fac0 = (tem - t_ice) / (tem * t_ice)
        fac1 = fac0 * li2
        fac2 = (d2ice * log(tem / t_ice) + fac1) / rvgas
        table(i) = e00 * exp(fac2)
    enddo
 
    ! -----------------------------------------------------------------------
    ! compute es over water between - 20 deg c and 102 deg c.
    ! -----------------------------------------------------------------------
 
    do i = 1, 1221
        tem = 253.16 + delt * real(i - 1)
        fac0 = (tem - t_ice) / (tem * t_ice)
        fac1 = fac0 * lv0
        fac2 = (dc_vap * log(tem / t_ice) + fac1) / rvgas
        esh20 = e00 * exp(fac2)
        if (i <= 200) then
            esupc(i) = esh20
        else
            table(i + 1400) = esh20
        endif
    enddo
 
    ! -----------------------------------------------------------------------
    ! derive blended es over ice and supercooled water between - 20 deg c and 0 deg c
    ! -----------------------------------------------------------------------
 
    do i = 1, 200
        tem = 253.16 + delt * real(i - 1)
        wice = 0.05 * (table_ice - tem)
        wh2o = 0.05 * (tem - 253.16)
        table(i + 1400) = wice * table(i + 1400) + wh2o * esupc(i)
    enddo
 
end subroutine qs_table
 
! =======================================================================
! compute the saturated specific humidity and the gradient of saturated specific humidity
! input t in deg k, p in pa; p = rho rdry tv, moist pressure
!@details It als oincludes the option for computing des/dT.
subroutine qsmith (im, km, ks, t, p, q, qs, dqdt)
 
    implicit none
 
    integer, intent (in) :: im, km, ks
 
    real, intent (in), dimension (im, km) :: t, p, q
 
    real, intent (out), dimension (im, km) :: qs
 
    real, intent (out), dimension (im, km), optional :: dqdt
 
    real :: eps10, ap1, tmin
 
    real, dimension (im, km) :: es
 
    integer :: i, k, it
 
    tmin = table_ice - 160.
    eps10 = 10. * eps
 
    if (.not. tables_are_initialized) then
        call qsmith_init
    endif
 
    do k = ks, km
        do i = 1, im
            ap1 = 10. * dim(t(i, k), tmin) + 1.
            ap1 = min(2621., ap1)
            it = ap1
            es(i, k) = table(it) + (ap1 - it) * des(it)
            qs(i, k) = eps * es(i, k) * (1. + zvir * q(i, k)) / p(i, k)
        enddo
    enddo
 
    if (present (dqdt)) then
        do k = ks, km
            do i = 1, im
                ap1 = 10. * dim(t(i, k), tmin) + 1.
                ap1 = min(2621., ap1) - 0.5
                it = ap1
                dqdt(i, k) = eps10 * (des(it) + (ap1 - it) * (des(it + 1) - des(it))) * (1. + zvir * q(i, k)) / p(i, k)
            enddo
        enddo
    endif
 
end subroutine qsmith
 
! =======================================================================
subroutine neg_adj (ktop, kbot, pt, dp, qv, ql, qr, qi, qs, qg)
 
    implicit none
 
    integer, intent (in) :: ktop, kbot
 
    real, intent (in), dimension (ktop:kbot) :: dp
 
    real, intent (inout), dimension (ktop:kbot) :: pt, qv, ql, qr, qi, qs, qg
 
    real, dimension (ktop:kbot) :: lcpk, icpk
 
    real :: dq, cvm
 
    integer :: k
 
    ! -----------------------------------------------------------------------
    ! define heat capacity and latent heat coefficient
    ! -----------------------------------------------------------------------
 
    do k = ktop, kbot
        cvm = c_air + qv(k) * c_vap + (qr(k) + ql(k)) * c_liq + (qi(k) + qs(k) + qg(k)) * c_ice
        lcpk(k) = (lv00 + d0_vap * pt(k)) / cvm
        icpk(k) = (li00 + dc_ice * pt(k)) / cvm
    enddo
 
    do k = ktop, kbot
 
        ! -----------------------------------------------------------------------
        ! ice phase:
        ! -----------------------------------------------------------------------
 
        ! if cloud ice < 0, borrow from snow
        if (qi(k) < 0.) then
            qs(k) = qs(k) + qi(k)
            qi(k) = 0.
        endif
        ! if snow < 0, borrow from graupel
        if (qs(k) < 0.) then
            qg(k) = qg(k) + qs(k)
            qs(k) = 0.
        endif
        ! if graupel < 0, borrow from rain
        if (qg(k) < 0.) then
            qr(k) = qr(k) + qg(k)
            pt(k) = pt(k) - qg(k) * icpk(k) ! heating
            qg(k) = 0.
        endif
 
        ! -----------------------------------------------------------------------
        ! liquid phase:
        ! -----------------------------------------------------------------------
 
        ! if rain < 0, borrow from cloud water
        if (qr(k) < 0.) then
            ql(k) = ql(k) + qr(k)
            qr(k) = 0.
        endif
        ! if cloud water < 0, borrow from water vapor
        if (ql(k) < 0.) then
            qv(k) = qv(k) + ql(k)
            pt(k) = pt(k) - ql(k) * lcpk(k) ! heating
            ql(k) = 0.
        endif
 
    enddo
 
    ! -----------------------------------------------------------------------
    ! fix water vapor; borrow from below
    ! -----------------------------------------------------------------------
 
    do k = ktop, kbot - 1
        if (qv(k) < 0.) then
            qv(k + 1) = qv(k + 1) + qv(k) * dp(k) / dp(k + 1)
            qv(k) = 0.
        endif
    enddo
 
    ! -----------------------------------------------------------------------
    ! bottom layer; borrow from above
    ! -----------------------------------------------------------------------
 
    if (qv(kbot) < 0. .and. qv(kbot - 1) > 0.) then
        dq = min(- qv(kbot) * dp(kbot), qv(kbot - 1) * dp(kbot - 1))
        qv(kbot - 1) = qv(kbot - 1) - dq / dp(kbot - 1)
        qv(kbot) = qv(kbot) + dq / dp(kbot)
    endif
 
end subroutine neg_adj
 
! =======================================================================
! compute global sum
! =======================================================================
 
!real function g_sum (p, ifirst, ilast, jfirst, jlast, area, mode)
!
! use mpp_mod, only: mpp_sum
!
! implicit none
!
! integer, intent (in) :: ifirst, ilast, jfirst, jlast
! integer, intent (in) :: mode ! if == 1 divided by area
!
! real, intent (in), dimension (ifirst:ilast, jfirst:jlast) :: p, area
!
! integer :: i, j
!
! real :: gsum
!
! if (global_area < 0.) then
! global_area = 0.
! do j = jfirst, jlast
! do i = ifirst, ilast
! global_area = global_area + area (i, j)
! enddo
! enddo
! call mpp_sum (global_area)
! endif
!
! gsum = 0.
! do j = jfirst, jlast
! do i = ifirst, ilast
! gsum = gsum + p (i, j) * area (i, j)
! enddo
! enddo
! call mpp_sum (gsum)
!
! if (mode == 1) then
! g_sum = gsum / global_area
! else
! g_sum = gsum
! endif
!
!end function g_sum
 
! ==========================================================================
subroutine interpolate_z (is, ie, js, je, km, zl, hgt, a3, a2)
 
    implicit none
 
    integer, intent (in) :: is, ie, js, je, km
 
    real, intent (in), dimension (is:ie, js:je, km) :: a3
 
    real, intent (in), dimension (is:ie, js:je, km + 1) :: hgt ! hgt (k) > hgt (k + 1)
 
    real, intent (in) :: zl
 
    real, intent (out), dimension (is:ie, js:je) :: a2
 
    real, dimension (km) :: zm
 
    integer :: i, j, k
 
    !$omp parallel do default (none) shared (is, ie, js, je, km, hgt, zl, a2, a3) private (zm)
 
    do j = js, je
        do i = is, ie
            do k = 1, km
                zm(k) = 0.5 * (hgt(i, j, k) + hgt(i, j, k + 1))
            enddo
            if (zl >= zm(1)) then
                a2(i, j) = a3(i, j, 1)
            elseif (zl <= zm(km)) then
                a2(i, j) = a3(i, j, km)
            else
                do k = 1, km - 1
                    if (zl <= zm(k) .and. zl >= zm(k + 1)) then
                        a2(i, j) = a3(i, j, k) + (a3(i, j, k + 1) - a3(i, j, k)) * (zm(k) - zl) / (zm(k) - zm(k + 1))
                        exit
                    endif
                enddo
            endif
        enddo
    enddo
 
end subroutine interpolate_z
 
! =======================================================================
! =======================================================================
subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms, qmg, t, &
        rew, rei, rer, res, reg)
 
    implicit none
 
    integer, intent (in) :: is, ie, ks, ke
    integer, intent (in), dimension (is:ie) :: lsm ! land sea mask, 0: ocean, 1: land, 2: sea ice
 
    real, intent (in), dimension (is:ie, ks:ke) :: den, delp, t
    real, intent (in), dimension (is:ie, ks:ke) :: qmw, qmi, qmr, qms, qmg
 
    real, intent (out), dimension (is:ie, ks:ke) :: rew, rei, rer, res, reg
 
    real, dimension (is:ie, ks:ke) :: qcw, qci, qcr, qcs, qcg
    
    integer :: i, k
 
    real :: lambdar, lambdas, lambdag
    real :: dpg, rei_fac, mask, ccn, bw
    real, parameter :: rho_0 = 50.e-3
 
    real :: rhow = 1.0e3, rhor = 1.0e3, rhos = 1.0e2, rhog = 4.0e2
    real :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6
    real :: alphar = 0.8, alphas = 0.25, alphag = 0.5
    real :: gammar = 17.837789, gammas = 8.2850630, gammag = 11.631769
    real :: qmin = 1.0e-12, beta = 1.22, qmin1 = 9.e-6
 
    do k = ks, ke
        do i = is, ie
            
            dpg = abs(delp(i, k)) / grav
            mask = min(max(real(lsm(i)), 0.0), 2.0)
 
            ! -----------------------------------------------------------------------
            ! cloud water (Martin et al., 1994)
            ! -----------------------------------------------------------------------
 
            ccn = 0.80 * (- 1.15e-3 * (ccn_o ** 2) + 0.963 * ccn_o + 5.30) * abs(mask - 1.0) + &
                0.67 * (- 2.10e-4 * (ccn_l ** 2) + 0.568 * ccn_l - 27.9) * (1.0 - abs(mask - 1.0))
 
            if (qmw(i, k) .gt. qmin) then
                qcw(i, k) = dpg * qmw(i, k) * 1.0e3
                rew(i, k) = exp(1.0 / 3.0 * log((3.0 * den(i, k) * qmw(i, k)) / (4.0 * pi * rhow * ccn))) * 1.0e4
                rew(i, k) = max(rewmin, min(rewmax, rew(i, k)))
            else
                qcw(i, k) = 0.0
                rew(i, k) = rewmin
            endif
            
            if (reiflag .eq. 1) then
 
            ! -----------------------------------------------------------------------
                ! cloud ice (Heymsfield and Mcfarquhar, 1996)
            ! -----------------------------------------------------------------------
 
                if (qmi(i, k) .gt. qmin1) then
                    qci(i, k) = dpg * qmi(i, k) * 1.0e3
                    rei_fac = log(1.0e3 * qmi(i, k) * den(i, k))
                    if (t(i, k) - tice .lt. - 50) then
                        rei(i, k) = beta / 9.917 * exp(0.109 * rei_fac) * 1.0e3
                    elseif (t(i, k) - tice .lt. - 40) then
                        rei(i, k) = beta / 9.337 * exp(0.080 * rei_fac) * 1.0e3
                    elseif (t(i, k) - tice .lt. - 30) then
                        rei(i, k) = beta / 9.208 * exp(0.055 * rei_fac) * 1.0e3
                else
                        rei(i, k) = beta / 9.387 * exp(0.031 * rei_fac) * 1.0e3
                endif
                    rei(i, k) = max(reimin, min(reimax, rei(i, k)))
            else
                    qci(i, k) = 0.0
                    rei(i, k) = reimin
            endif
 
            endif
            
            if (reiflag .eq. 2) then
 
            ! -----------------------------------------------------------------------
                ! cloud ice (Wyser, 1998)
            ! -----------------------------------------------------------------------
 
                if (qmi(i, k) .gt. qmin1) then
                    qci(i, k) = dpg * qmi(i, k) * 1.0e3
                    bw = - 2. + 1.e-3 * log10(den(i, k) * qmi(i, k) / rho_0) * max(0.0, tice - t(i, k)) ** 1.5
                    rei(i, k) = 377.4 + bw * (203.3 + bw * (37.91 + 2.3696 * bw))
                    rei(i, k) = max(reimin, min(reimax, rei(i, k)))
            else
                    qci(i, k) = 0.0
                    rei(i, k) = reimin
                endif
 
            endif
 
            ! -----------------------------------------------------------------------
            ! rain (Lin et al., 1983)
            ! -----------------------------------------------------------------------
 
            if (qmr(i, k) .gt. qmin) then
                qcr(i, k) = dpg * qmr(i, k) * 1.0e3
                lambdar = exp(0.25 * log(pi * rhor * n0r / qmr(i, k) / den(i, k)))
                rer(i, k) = 0.5 * exp(log(gammar / 6) / alphar) / lambdar * 1.0e6
                rer(i, k) = max(rermin, min(rermax, rer(i, k)))
            else
                qcr(i, k) = 0.0
                rer(i, k) = rermin
            endif
 
            ! -----------------------------------------------------------------------
            ! snow (Lin et al., 1983)
            ! -----------------------------------------------------------------------
 
            if (qms(i, k) .gt. qmin1) then
                qcs(i, k) = dpg * qms(i, k) * 1.0e3
                lambdas = exp(0.25 * log(pi * rhos * n0s / qms(i, k) / den(i, k)))
                res(i, k) = 0.5 * exp(log(gammas / 6) / alphas) / lambdas * 1.0e6
                res(i, k) = max(resmin, min(resmax, res(i, k)))
            else
                qcs(i, k) = 0.0
                res(i, k) = resmin
            endif
            
            ! -----------------------------------------------------------------------
            ! graupel (Lin et al., 1983)
            ! -----------------------------------------------------------------------
            
            if (qmg(i, k) .gt. qmin) then
                qcg(i, k) = dpg * qmg(i, k) * 1.0e3
                lambdag = exp(0.25 * log(pi * rhog * n0g / qmg(i, k) / den(i, k)))
                reg(i, k) = 0.5 * exp(log(gammag / 6) / alphag) / lambdag * 1.0e6
                reg(i, k) = max(regmin, min(regmax, reg(i, k)))
            else
                qcg(i, k) = 0.0
                reg(i, k) = regmin
            endif
 
        enddo
    enddo
 
end subroutine cloud_diagnosis
 
!+---+-----------------------------------------------------------------+
      subroutine refl10cm_gfdl (qv1d, qr1d, qs1d, qg1d,                 &
                       t1d, p1d, dBZ, kts, kte, ii, jj, melti)
 
      IMPLICIT NONE
 
!..Sub arguments
      INTEGER, INTENT(IN):: kts, kte, ii,jj
      REAL, DIMENSION(kts:kte), INTENT(IN)::                            &
                      qv1d, qr1d, qs1d, qg1d, t1d, p1d
      REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
 
!..Local variables
      REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
      REAL, DIMENSION(kts:kte):: rr, rs, rg
!      REAL:: temp_C
 
      DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams, ilamg
      DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s, N0_g
      DOUBLE PRECISION:: lamr, lams, lamg
      LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
 
      REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
      DOUBLE PRECISION:: fmelt_s, fmelt_g
 
      INTEGER:: i, k, k_0, kbot, n
      LOGICAL, INTENT(IN):: melti
      DOUBLE PRECISION:: cback, x, eta, f_d
!+---+
 
      do k = kts, kte
         dbz(k) = -35.0
      enddo
 
!+---+-----------------------------------------------------------------+
!..Put column of data into local arrays.
!+---+-----------------------------------------------------------------+
      do k = kts, kte
         temp(k) = t1d(k)
!         temp_C = min(-0.001, temp(K)-273.15)
         qv(k) = max(1.e-10, qv1d(k))
         pres(k) = p1d(k)
         rho(k) = 0.622*pres(k)/(rdgas*temp(k)*(qv(k)+0.622))
 
         if (qr1d(k) .gt. 1.e-9) then
            rr(k) = qr1d(k)*rho(k)
            n0_r(k) = n0r
            lamr = (xam_r*xcrg(3)*n0_r(k)/rr(k))**(1./xcre(1))
            ilamr(k) = 1./lamr
            l_qr(k) = .true.
         else
            rr(k) = 1.e-12
            l_qr(k) = .false.
         endif
 
         if (qs1d(k) .gt. 1.e-9) then
            rs(k) = qs1d(k)*rho(k)
            n0_s(k) = n0s
            lams = (xam_s*xcsg(3)*n0_s(k)/rs(k))**(1./xcse(1))
            ilams(k) = 1./lams
            l_qs(k) = .true.
         else
            rs(k) = 1.e-12
            l_qs(k) = .false.
         endif
 
         if (qg1d(k) .gt. 1.e-9) then
            rg(k) = qg1d(k)*rho(k)
            n0_g(k) = n0g
            lamg = (xam_g*xcgg(3)*n0_g(k)/rg(k))**(1./xcge(1))
            ilamg(k) = 1./lamg
            l_qg(k) = .true.
         else
            rg(k) = 1.e-12
            l_qg(k) = .false.
         endif
      enddo
 
!+---+-----------------------------------------------------------------+
!..Locate K-level of start of melting (k_0 is level above).
!+---+-----------------------------------------------------------------+
      k_0 = kts
      k_loop:do k = kte-1, kts, -1
         if ( melti .and. (temp(k).gt.273.15) .and. l_qr(k)             &
              .and. (l_qs(k+1).or.l_qg(k+1)) ) then
            k_0 = max(k+1, k_0)
            EXIT k_loop
         endif
      enddo k_loop
!+---+-----------------------------------------------------------------+
!..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
!.. and non-water-coated snow and graupel when below freezing are
!.. simple. Integrations of m(D)*m(D)*N(D)*dD.
!+---+-----------------------------------------------------------------+
      do k = kts, kte
         ze_rain(k) = 1.e-22
         ze_snow(k) = 1.e-22
         ze_graupel(k) = 1.e-22
         if (l_qr(k)) ze_rain(k) = n0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
         if (l_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/pi)*(6.0/pi)     &
                                 * (xam_s/900.0)*(xam_s/900.0)          &
                                 * n0_s(k)*xcsg(4)*ilams(k)**xcse(4)
         if (l_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/pi)*(6.0/pi)  &
                                    * (xam_g/900.0)*(xam_g/900.0)       &
                                    * n0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
      enddo
 
 
!+---+-----------------------------------------------------------------+
!..Special case of melting ice (snow/graupel) particles.  Assume the
!.. ice is surrounded by the liquid water.  Fraction of meltwater is
!.. extremely simple based on amount found above the melting level.
!.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
!.. routines).
!+---+-----------------------------------------------------------------+
 
      if (melti .and. k_0.ge.kts+1) then
       do k = k_0-1, kts, -1
 
!..Reflectivity contributed by melting snow
          if (l_qs(k) .and. l_qs(k_0) ) then
           fmelt_s = max(0.005d0, min(1.0d0-rs(k)/rs(k_0), 0.99d0))
           eta = 0.d0
           lams = 1./ilams(k)
           do n = 1, nrbins
              x = xam_s * xxds(n)**xbm_s
              call rayleigh_soak_wetgraupel (x,dble(xocms),dble(xobms), &
                    fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
                    cback, mixingrulestring_s, matrixstring_s,          &
                    inclusionstring_s, hoststring_s,                    &
                    hostmatrixstring_s, hostinclusionstring_s)
              f_d = n0_s(k)*xxds(n)**xmu_s * dexp(-lams*xxds(n))
              eta = eta + f_d * cback * simpson(n) * xdts(n)
           enddo
           ze_snow(k) = sngl(lamda4 / (pi5 * k_w) * eta)
          endif
 
 
!..Reflectivity contributed by melting graupel
 
          if (l_qg(k) .and. l_qg(k_0) ) then
           fmelt_g = max(0.005d0, min(1.0d0-rg(k)/rg(k_0), 0.99d0))
           eta = 0.d0
           lamg = 1./ilamg(k)
           do n = 1, nrbins
              x = xam_g * xxdg(n)**xbm_g
              call rayleigh_soak_wetgraupel (x,dble(xocmg),dble(xobmg), &
                    fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
                    cback, mixingrulestring_g, matrixstring_g,          &
                    inclusionstring_g, hoststring_g,                    &
                    hostmatrixstring_g, hostinclusionstring_g)
              f_d = n0_g(k)*xxdg(n)**xmu_g * dexp(-lamg*xxdg(n))
              eta = eta + f_d * cback * simpson(n) * xdtg(n)
           enddo
           ze_graupel(k) = sngl(lamda4 / (pi5 * k_w) * eta)
          endif
 
       enddo
      endif
 
      do k = kte, kts, -1
         dbz(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
      enddo
 
 
      end subroutine refl10cm_gfdl
!+---+-----------------------------------------------------------------+
 
end module gfdl_cloud_microphys_mod