!>\file fv_sat_adj.F90
!! This file contains the GFDL in-core fast saturation adjustment.
!! and it is an "intermediate physics" implemented in the remapping Lagrangian to 
!! Eulerian loop of FV3 solver.
!***********************************************************************
!*                   GNU Lesser General Public License
!*
!* This file is part of the GFDL Cloud Microphysics.
!*
!* The GFDL Cloud Microphysics is free software: you can 
!8 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/>.
!***********************************************************************

!> This module contains the GFDL in-core fast saturation adjustment  
!! called in FV3 dynamics solver.
module fv_sat_adj
! Modules Included:
! <table>
! <tr>
!     <th>Module Name</th>
!     <th>Functions Included</th>
!   </tr>
!   <tr>
!     <td>constants_mod</td>
!     <td>rvgas, rdgas, grav, hlv, hlf, cp_air</td>
!   </tr>
!   <tr>
!     <td>fv_arrays_mod</td>
!     <td> r_grid</td>
!   </tr>
!   <tr>
!   <tr>
!     <td>fv_mp_mod</td>
!     <td>is_master</td>
!   </tr>
!   <tr>
!     <td>gfdl_cloud_microphys_mod</td>
!     <td>ql_gen, qi_gen, qi0_max, ql_mlt, ql0_max, qi_lim, qs_mlt,
!         tau_r2g, tau_smlt, tau_i2s, tau_v2l, tau_l2v, tau_imlt, tau_l2r,
!         rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land, tintqs</td>
!   </tr>
! </table>
    ! DH* TODO - MAKE THIS INPUT ARGUMENTS *DH
    use physcons, only : rdgas => con_rd_dyn, &
                         rvgas => con_rv_dyn, &
                         grav => con_g_dyn,   &
                         hlv => con_hvap_dyn, &
                         hlf => con_hfus_dyn, &
                         cp_air => con_cp_dyn
    ! *DH
    use machine,                  only: kind_grid, kind_dyn
    use gfdl_cloud_microphys_mod, only: ql_gen, qi_gen, qi0_max, ql_mlt, ql0_max, qi_lim, qs_mlt
    use gfdl_cloud_microphys_mod, only: icloud_f, sat_adj0, t_sub, cld_min
    use gfdl_cloud_microphys_mod, only: tau_r2g, tau_smlt, tau_i2s, tau_v2l, tau_l2v, tau_imlt, tau_l2r
    use gfdl_cloud_microphys_mod, only: rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land, tintqs
#ifdef MULTI_GASES
    use ccpp_multi_gases_mod, only: multi_gases_init,     &
                                    multi_gases_finalize, &
                                    virq_qpz, vicpqd_qpz, &
                                    vicvqd_qpz, num_gas
#endif

    implicit none
    
    private

    public fv_sat_adj_init, fv_sat_adj_run, fv_sat_adj_finalize

    logical :: is_initialized = .false.

    real(kind=kind_dyn), parameter :: rrg = -rdgas/grav
    ! real, parameter :: cp_air = cp_air           ! 1004.6, heat capacity of dry air at constant pressure, come from constants_mod
    real(kind=kind_dyn), parameter :: cp_vap = 4.0 * rvgas        !< 1846.0, heat capacity of water vapor at constant pressure
    real(kind=kind_dyn), parameter :: cv_air = cp_air - rdgas     !< 717.55, heat capacity of dry air at constant volume
    real(kind=kind_dyn), parameter :: cv_vap = 3.0 * rvgas        !< 1384.5, heat capacity of water vapor at constant volume
    ! http: / / www.engineeringtoolbox.com / ice - thermal - properties - d_576.html
    ! c_ice = 2050.0 at 0 deg c
    ! c_ice = 1972.0 at - 15 deg c
    ! c_ice = 1818.0 at - 40 deg c
    ! http: / / www.engineeringtoolbox.com / water - thermal - properties - d_162.html
    ! c_liq = 4205.0 at 4 deg c
    ! c_liq = 4185.5 at 15 deg c
    ! c_liq = 4178.0 at 30 deg c
    ! real, parameter :: c_ice = 2106.0            ! ifs: heat capacity of ice at 0 deg c
    ! real, parameter :: c_liq = 4218.0            ! ifs: heat capacity of liquid at 0 deg c
    real(kind=kind_dyn), parameter :: c_ice = 1972.0              !< gfdl: heat capacity of ice at - 15 deg c
    real(kind=kind_dyn), parameter :: c_liq = 4185.5              !< gfdl: heat capacity of liquid at 15 deg c
    real(kind=kind_dyn), parameter :: dc_vap = cp_vap - c_liq     !< - 2339.5, isobaric heating / cooling
    real(kind=kind_dyn), parameter :: dc_ice = c_liq - c_ice      !< 2213.5, isobaric heating / colling
    real(kind=kind_dyn), parameter :: tice = 273.16               !< freezing temperature
    real(kind=kind_dyn), parameter :: t_wfr = tice - 40.          !< homogeneous freezing temperature
    real(kind=kind_dyn), parameter :: lv0 = hlv - dc_vap * tice   !< 3.13905782e6, evaporation latent heat coefficient at 0 deg k
    real(kind=kind_dyn), parameter :: li00 = hlf - dc_ice * tice  !< - 2.7105966e5, fusion latent heat coefficient at 0 deg k
    ! real (kind_grid), parameter :: e00 = 610.71  ! gfdl: saturation vapor pressure at 0 deg c
    real (kind_grid), parameter  :: e00 = 611.21    !< ifs: saturation vapor pressure at 0 deg c
    real (kind_grid), parameter :: d2ice = dc_vap + dc_ice !< - 126, isobaric heating / cooling
    real (kind_grid), parameter :: li2 = lv0 + li00        !< 2.86799816e6, sublimation latent heat coefficient at 0 deg k
    real(kind=kind_dyn), parameter :: lat2 = (hlv + hlf) ** 2     !< used in bigg mechanism
    real(kind=kind_dyn) :: d0_vap                                 !< the same as dc_vap, except that cp_vap can be cp_vap or cv_vap
    real(kind=kind_dyn) :: lv00                                   !< the same as lv0, except that cp_vap can be cp_vap or cv_vap
    real(kind=kind_dyn), allocatable :: table (:), table2 (:), tablew (:), des2 (:), desw (:)

contains

!>\brief The subroutine 'fv_sat_adj_init' initializes lookup tables for the saturation mixing ratio.
!! \section arg_table_fv_sat_adj_init Argument Table
!! \htmlinclude fv_sat_adj_init.html
!!
subroutine fv_sat_adj_init(do_sat_adj, kmp, nwat, ngas, rilist, cpilist, &
                           mpirank, mpiroot, errmsg, errflg)

    implicit none

    ! Interface variables
    logical,          intent(in   ) :: do_sat_adj
    integer,          intent(in   ) :: kmp
    integer,          intent(in   ) :: nwat
    integer,          intent(in   ) :: ngas
    real(kind_dyn),   intent(in   ) :: rilist(0:ngas)
    real(kind_dyn),   intent(in   ) :: cpilist(0:ngas)
    integer,          intent(in   ) :: mpirank
    integer,          intent(in   ) :: mpiroot
    character(len=*), intent(  out) :: errmsg
    integer,          intent(  out) :: errflg

    ! Local variables
    integer, parameter :: length = 2621
    integer :: i

    ! Initialize the CCPP error handling variables
    errmsg = ''
    errflg = 0

    ! If saturation adjustment is not used, return immediately
    if (.not.do_sat_adj) then
      write(errmsg,'(a)') 'Logic error: fv_sat_adj_init is called but do_sat_adj is set to false'
      errflg = 1
      return
    end if

    if (.not.nwat==6) then
      write(errmsg,'(a)') 'Logic error: fv_sat_adj requires six water species (nwat=6)'
      errflg = 1
      return
    end if

    if (is_initialized) return

    ! generate es table (dt = 0.1 deg c)

    allocate (table (length))
    allocate (table2 (length))
    allocate (tablew (length))
    allocate (des2 (length))
    allocate (desw (length))

    call qs_table (length)
    call qs_table2 (length)
    call qs_tablew (length)

    do i = 1, length - 1
        des2 (i) = max (0., table2 (i + 1) - table2 (i))
        desw (i) = max (0., tablew (i + 1) - tablew (i))
    enddo
    des2 (length) = des2 (length - 1)
    desw (length) = desw (length - 1)

#ifdef MULTI_GASES
    call multi_gases_init(ngas,nwat,rilist,cpilist,mpirank==mpiroot)
#endif

    is_initialized = .true.

end subroutine fv_sat_adj_init

!\ingroup fast_sat_adj
!>\brief The subroutine 'fv_sat_adj_finalize' deallocates lookup tables for the saturation mixing ratio.
!! \section arg_table_fv_sat_adj_finalize Argument Table
!! \htmlinclude fv_sat_adj_finalize.html
!!
subroutine fv_sat_adj_finalize (errmsg, errflg)

    implicit none

    character(len=*), intent(out) :: errmsg
    integer,          intent(out) :: errflg

    ! Initialize the CCPP error handling variables
    errmsg = ''
    errflg = 0

    if (.not.is_initialized) return

    if (allocated(table )) deallocate(table )
    if (allocated(table2)) deallocate(table2)
    if (allocated(tablew)) deallocate(tablew)
    if (allocated(des2  )) deallocate(des2  )
    if (allocated(desw  )) deallocate(desw  )

#ifdef MULTI_GASES
    call multi_gases_finalize()
#endif

    is_initialized = .false.

end subroutine fv_sat_adj_finalize

!>\defgroup fast_sat_adj GFDL In-Core Fast Saturation Adjustment Module
!> @{
!! The subroutine 'fv_sat_adj' implements the fast processes in the GFDL
!! Cloud MP. It is part of the GFDL Cloud MP.
!>\author Shian-Jiann Lin, Linjiong Zhou
!!
!>\brief The subroutine 'fv_sat_adj' performs the fast processes in the GFDL microphysics.
!>\details This is designed for single-moment 6-class cloud microphysics schemes.
!! It handles the heat release due to in situ phase changes.
!!
!! \section arg_table_fv_sat_adj_run Argument Table
!! \htmlinclude fv_sat_adj_run.html
!!
subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, kmp, km, kmdelz, js, je, jsd, jed, &
                 ng, hydrostatic, fast_mp_consv, te0_2d, te0, ngas, qvi, qv, ql, qi, qr,  &
                 qs, qg, hs, peln, delz, delp, pt, pkz, q_con, akap, cappa, area, dtdt,   &
                 out_dt, last_step, do_qa, qa,                                            &
                 nthreads, errmsg, errflg)

    implicit none

    ! Interface variables
    real(kind=kind_dyn), intent(in)    :: mdt
    real(kind=kind_dyn), intent(in)    :: zvir
    integer,             intent(in)    :: is
    integer,             intent(in)    :: ie
    integer,             intent(in)    :: isd
    integer,             intent(in)    :: ied
    integer,             intent(in)    :: kmp
    integer,             intent(in)    :: km
    integer,             intent(in)    :: kmdelz
    integer,             intent(in)    :: js
    integer,             intent(in)    :: je
    integer,             intent(in)    :: jsd
    integer,             intent(in)    :: jed
    integer,             intent(in)    :: ng
    logical,             intent(in)    :: hydrostatic
    logical,             intent(in)    :: fast_mp_consv
    real(kind=kind_dyn), intent(inout) :: te0_2d(is:ie, js:je)
    real(kind=kind_dyn), intent(  out) :: te0(isd:ied, jsd:jed, 1:km)
    ! If multi-gases physics are not used, ngas is one and qvi identical to qv
    integer,             intent(in)    :: ngas
    real(kind=kind_dyn), intent(inout) :: qvi(isd:ied, jsd:jed, 1:km, 1:ngas)
    real(kind=kind_dyn), intent(inout) :: qv(isd:ied, jsd:jed, 1:km)
    real(kind=kind_dyn), intent(inout) :: ql(isd:ied, jsd:jed, 1:km)
    real(kind=kind_dyn), intent(inout) :: qi(isd:ied, jsd:jed, 1:km)
    real(kind=kind_dyn), intent(inout) :: qr(isd:ied, jsd:jed, 1:km)
    real(kind=kind_dyn), intent(inout) :: qs(isd:ied, jsd:jed, 1:km)
    real(kind=kind_dyn), intent(inout) :: qg(isd:ied, jsd:jed, 1:km)
    real(kind=kind_dyn), intent(in)    :: hs(isd:ied, jsd:jed)
    real(kind=kind_dyn), intent(in)    :: peln(is:ie, 1:km+1, js:je)
    ! For hydrostatic build, kmdelz=1, otherwise kmdelz=km (see fv_arrays.F90)
    real(kind=kind_dyn), intent(in)    :: delz(is:ie, js:je, 1:kmdelz)
    real(kind=kind_dyn), intent(in)    :: delp(isd:ied, jsd:jed, 1:km)
    real(kind=kind_dyn), intent(inout) :: pt(isd:ied, jsd:jed, 1:km)
    real(kind=kind_dyn), intent(inout) :: pkz(is:ie, js:je, 1:km)
#ifdef USE_COND
    real(kind=kind_dyn), intent(inout) :: q_con(isd:ied, jsd:jed, 1:km)
#else
    real(kind=kind_dyn), intent(inout) :: q_con(isd:isd, jsd:jsd, 1)
#endif
    real(kind=kind_dyn), intent(in)    :: akap
#ifdef MOIST_CAPPA
    real(kind=kind_dyn), intent(inout) :: cappa(isd:ied, jsd:jed, 1:km)
#else
    real(kind=kind_dyn), intent(inout) :: cappa(isd:ied, jsd:jed, 1)
#endif
    ! DH* WARNING, allocation in fv_arrays.F90 is area(isd_2d:ied_2d, jsd_2d:jed_2d),
    ! where normally isd_2d = isd etc, but for memory optimization, these can be set
    ! to isd_2d=0, ied_2d=-1 etc. I don't believe this optimization is actually used,
    ! as it would break a whole lot of code (including the one below)!
    ! Assume thus that isd_2d = isd etc.
    real(kind_grid),     intent(in)    :: area(isd:ied, jsd:jed)
    real(kind=kind_dyn), intent(inout) :: dtdt(is:ie, js:je, 1:km)
    logical,             intent(in)    :: out_dt
    logical,             intent(in)    :: last_step
    logical,             intent(in)    :: do_qa
    real(kind=kind_dyn), intent(  out) :: qa(isd:ied, jsd:jed, 1:km)
    integer,             intent(in)    :: nthreads
    character(len=*),    intent(  out) :: errmsg
    integer,             intent(  out) :: errflg

    ! Local variables
    real(kind=kind_dyn), dimension(is:ie,js:je) :: dpln
    integer :: kdelz
    integer :: k, j, i

    ! Initialize the CCPP error handling variables
    errmsg = ''
    errflg = 0

#ifndef FV3
! Open parallel region if not already opened by host model
!$OMP parallel num_threads(nthreads) default(none) &
!$OMP          shared(kmp,km,js,je,is,ie,peln,mdt, &
!$OMP                 isd,jsd,delz,q_con,cappa,qa, &
!$OMP                 do_qa,last_step,out_dt,dtdt, &
!$OMP                 area,delp,pt,hs,qg,qs,qr,qi, &
!$OMP                 ql,qv,te0,fast_mp_consv,     &
!$OMP                 hydrostatic,ng,zvir,pkz,     &
!$OMP                 akap,te0_2d,ngas,qvi)        &
!$OMP          private(k,j,i,kdelz,dpln)
#endif

!$OMP do
    do k=kmp,km
       do j=js,je
          do i=is,ie
             dpln(i,j) = peln(i,k+1,j) - peln(i,k,j)
          enddo
       enddo
       if (hydrostatic) then
          kdelz = 1
       else
          kdelz = k
       end if
       call fv_sat_adj_work(abs(mdt), zvir, is, ie, js, je, ng, hydrostatic, fast_mp_consv, &
                            te0(isd,jsd,k),                                                 &
#ifdef MULTI_GASES
                            qvi(isd,jsd,k,1:ngas),                                          &
#else
                            qv(isd,jsd,k),                                                  &
#endif
                            ql(isd,jsd,k), qi(isd,jsd,k),                                   &
                            qr(isd,jsd,k), qs(isd,jsd,k), qg(isd,jsd,k),                    &
                            hs, dpln, delz(is:,js:,kdelz), pt(isd,jsd,k), delp(isd,jsd,k),&
                            q_con(isd:,jsd:,k), cappa(isd:,jsd:,k), area, dtdt(is,js,k),    &
                            out_dt, last_step, do_qa, qa(isd,jsd,k))
       if ( .not. hydrostatic  ) then
          do j=js,je
             do i=is,ie
#ifdef MOIST_CAPPA
               pkz(i,j,k) = exp(cappa(i,j,k)*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
#else
#ifdef MULTI_GASES
               pkz(i,j,k) = exp(akap*(virqd(q(i,j,k,1:num_gas))/vicpqd(q(i,j,k,1:num_gas))*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
#else
               pkz(i,j,k) = exp(akap*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
#endif
#endif
             enddo
          enddo
       endif
    enddo
!$OMP end do

    if ( fast_mp_consv ) then
!$OMP do
       do j=js,je
          do i=is,ie
             do k=kmp,km
                te0_2d(i,j) = te0_2d(i,j) + te0(i,j,k)
             enddo
          enddo
       enddo
!$OMP end do
    endif

#ifndef FV3
!$OMP end parallel
#endif

    return

end subroutine fv_sat_adj_run

!>\ingroup fast_sat_adj
!> This subroutine includes the entity of the fast saturation adjustment processes.
!>\section fast_gen GFDL Cloud Fast Physics General Algorithm
!> @{
subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, te0, &
#ifdef MULTI_GASES
        qvi, &
#else 
        qv, &
#endif 
        ql, qi, qr, qs, qg, hs, dpln, delz, pt, dp, q_con, cappa, &
        area, dtdt, out_dt, last_step, do_qa, qa)

    implicit none

    ! Interface variables
    integer, intent (in) :: is, ie, js, je, ng
    logical, intent (in) :: hydrostatic, consv_te, out_dt, last_step, do_qa
    real(kind=kind_dyn), intent (in) :: zvir, mdt ! remapping time step
    real(kind=kind_dyn), intent (in), dimension (is - ng:ie + ng, js - ng:je + ng) :: dp,  hs
    real(kind=kind_dyn), intent (in), dimension (is:ie, js:je) :: dpln, delz
    real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: pt
#ifdef MULTI_GASES
    real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng, 1:1, 1:num_gas) :: qvi
#else
    real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: qv
#endif
    real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: ql, qi, qr, qs, qg
    real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: q_con, cappa
    real(kind=kind_dyn), intent (inout), dimension (is:ie, js:je) :: dtdt
    real(kind=kind_dyn), intent (out), dimension (is - ng:ie + ng, js - ng:je + ng) :: qa, te0
    real (kind_grid), intent (in), dimension (is - ng:ie + ng, js - ng:je + ng) :: area

    ! Local variables
#ifdef MULTI_GASES
    real, dimension (is - ng:ie + ng, js - ng:je + ng) :: qv
#endif
    real(kind=kind_dyn), dimension (is:ie) :: wqsat, dq2dt, qpz, cvm, t0, pt1, qstar
    real(kind=kind_dyn), dimension (is:ie) :: icp2, lcp2, tcp2, tcp3
    real(kind=kind_dyn), dimension (is:ie) :: den, q_liq, q_sol, q_cond, src, sink, hvar
    real(kind=kind_dyn), dimension (is:ie) :: mc_air, lhl, lhi
    real(kind=kind_dyn) :: qsw, rh
    real(kind=kind_dyn) :: tc, qsi, dqsdt, dq, dq0, pidep, qi_crt, tmp, dtmp
    real(kind=kind_dyn) :: tin, rqi, q_plus, q_minus
    real(kind=kind_dyn) :: sdt, dt_bigg, adj_fac
    real(kind=kind_dyn) :: fac_smlt, fac_r2g, fac_i2s, fac_imlt, fac_l2r, fac_v2l, fac_l2v
    real(kind=kind_dyn) :: factor, qim, tice0, c_air, c_vap, dw
    integer :: i, j

#ifdef MULTI_GASES
    qv(:,:) = qvi(:,:,1,1)
#endif
    sdt = 0.5 * mdt                   ! half remapping time step
    dt_bigg = mdt                     ! bigg mechinism time step
    
    tice0 = tice - 0.01               ! 273.15, standard freezing temperature
    
    ! -----------------------------------------------------------------------
    !> - Define conversion scalar / factor.
    ! -----------------------------------------------------------------------
    
    fac_i2s = 1. - exp (- mdt / tau_i2s)
    fac_v2l = 1. - exp (- sdt / tau_v2l)
    fac_r2g = 1. - exp (- mdt / tau_r2g)
    fac_l2r = 1. - exp (- mdt / tau_l2r)
    
    fac_l2v = 1. - exp (- sdt / tau_l2v)
    fac_l2v = min (sat_adj0, fac_l2v)
    
    fac_imlt = 1. - exp (- sdt / tau_imlt)
    fac_smlt = 1. - exp (- mdt / tau_smlt)
    
    ! -----------------------------------------------------------------------
    !> - Define heat capacity of dry air and water vapor based on hydrostatical property.
    ! -----------------------------------------------------------------------
    
    if (hydrostatic) then
        c_air = cp_air
        c_vap = cp_vap
    else
        c_air = cv_air
        c_vap = cv_vap
    endif
    d0_vap = c_vap - c_liq
    lv00 = hlv - d0_vap * tice
    ! dc_vap = cp_vap - c_liq ! - 2339.5
    ! d0_vap = cv_vap - c_liq ! - 2801.0
    
    do j = js, je ! start j loop
        
        do i = is, ie
            q_liq (i) = ql (i, j) + qr (i, j)
            q_sol (i) = qi (i, j) + qs (i, j) + qg (i, j)
            qpz (i) = q_liq (i) + q_sol (i)
#ifdef MULTI_GASES
            pt1 (i) = pt (i, j) / virq_qpz(qvi(i,j,1,1:num_gas),qpz(i))
#else
#ifdef USE_COND
            pt1 (i) = pt (i, j) / ((1 + zvir * qv (i, j)) * (1 - qpz (i)))
#else
            pt1 (i) = pt (i, j) / (1 + zvir * qv (i, j))
#endif
#endif
            t0 (i) = pt1 (i) ! true temperature
            qpz (i) = qpz (i) + qv (i, j) ! total_wat conserved in this routine
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Define air density based on hydrostatical property.
        ! -----------------------------------------------------------------------
        
        if (hydrostatic) then
            do i = is, ie
                den (i) = dp (i, j) / (dpln (i, j) * rdgas * pt (i, j))
            enddo
        else
            do i = is, ie
                den (i) = - dp (i, j) / (grav * delz (i, j)) ! moist_air density
            enddo
        endif
        
        ! -----------------------------------------------------------------------
        !> - Define heat capacity and latend heat coefficient.
        ! -----------------------------------------------------------------------
        
        do i = is, ie
#ifdef MULTI_GASES
            if (hydrostatic) then
                c_air = cp_air * vicpqd_qpz(qvi(i,j,1,1:num_gas),qpz(i))
            else
                c_air = cv_air * vicvqd_qpz(qvi(i,j,1,1:num_gas),qpz(i))
            endif
#endif
            mc_air (i) = (1. - qpz (i)) * c_air ! constant
            cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
            lhi (i) = li00 + dc_ice * pt1 (i)
            icp2 (i) = lhi (i) / cvm (i)
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Fix energy conservation.
        ! -----------------------------------------------------------------------
        
        if (consv_te) then
            if (hydrostatic) then
                do i = is, ie
#ifdef MULTI_GASES
                    c_air = cp_air * vicpqd_qpz(qvi(i,j,1,1:num_gas),qpz(i))
#endif
                    te0 (i, j) = - c_air * t0 (i)
                enddo
            else
                do i = is, ie
#ifdef USE_COND
                    te0 (i, j) = - cvm (i) * t0 (i)
#else
#ifdef MULTI_GASES
                    c_air = cv_air * vicvqd_qpz(qvi(i,j,1,1:num_gas),qpz(i))
#endif
                    te0 (i, j) = - c_air * t0 (i)
#endif
                enddo
            endif
        endif
        
        ! -----------------------------------------------------------------------
        !> - Fix negative cloud ice with snow.
        ! -----------------------------------------------------------------------
        
            do i = is, ie
            if (qi (i, j) < 0.) then
                qs (i, j) = qs (i, j) + qi (i, j)
                qi (i, j) = 0.
            endif
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Melting of cloud ice to cloud water and rain.
        ! -----------------------------------------------------------------------
        
        do i = is, ie
            if (qi (i, j) > 1.e-8 .and. pt1 (i) > tice) then
                sink (i) = min (qi (i, j), fac_imlt * (pt1 (i) - tice) / icp2 (i))
                qi (i, j) = qi (i, j) - sink (i)
                ! sjl, may 17, 2017
                ! tmp = min (sink (i), dim (ql_mlt, ql (i, j))) ! max ql amount
                ! ql (i, j) = ql (i, j) + tmp
                ! qr (i, j) = qr (i, j) + sink (i) - tmp
                ! sjl, may 17, 2017
                ql (i, j) = ql (i, j) + sink (i)
                q_liq (i) = q_liq (i) + sink (i)
                q_sol (i) = q_sol (i) - sink (i)
                cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
                pt1 (i) = pt1 (i) - sink (i) * lhi (i) / cvm (i)
            endif
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Update latend heat coefficient.
        ! -----------------------------------------------------------------------
        
            do i = is, ie
            lhi (i) = li00 + dc_ice * pt1 (i)
            icp2 (i) = lhi (i) / cvm (i)
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Fix negative snow with graupel or graupel with available snow.
        ! -----------------------------------------------------------------------
        
        do i = is, ie
            if (qs (i, j) < 0.) then
                qg (i, j) = qg (i, j) + qs (i, j)
                qs (i, j) = 0.
            elseif (qg (i, j) < 0.) then
                tmp = min (- qg (i, j), max (0., qs (i, j)))
                qg (i, j) = qg (i, j) + tmp
                qs (i, j) = qs (i, j) - tmp
            endif
        enddo
        
        ! after this point cloud ice & snow are positive definite
        
        ! -----------------------------------------------------------------------
        !> - Fix negative cloud water with rain or rain with available cloud water.
        ! -----------------------------------------------------------------------
        
            do i = is, ie
            if (ql (i, j) < 0.) then
                tmp = min (- ql (i, j), max (0., qr (i, j)))
                ql (i, j) = ql (i, j) + tmp
                qr (i, j) = qr (i, j) - tmp
            elseif (qr (i, j) < 0.) then
                tmp = min (- qr (i, j), max (0., ql (i, j)))
                ql (i, j) = ql (i, j) - tmp
                qr (i, j) = qr (i, j) + tmp
            endif
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Enforce complete freezing of cloud water to cloud ice below - 48 c.
        ! -----------------------------------------------------------------------

        do i = is, ie
            dtmp = tice - 48. - pt1 (i)
            if (ql (i, j) > 0. .and. dtmp > 0.) then
                sink (i) = min (ql (i, j), dtmp / icp2 (i))
                ql (i, j) = ql (i, j) - sink (i)
                qi (i, j) = qi (i, j) + sink (i)
                q_liq (i) = q_liq (i) - sink (i)
                q_sol (i) = q_sol (i) + sink (i)
                cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
                pt1 (i) = pt1 (i) + sink (i) * lhi (i) / cvm (i)
            endif
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Update latend heat coefficient.
        ! -----------------------------------------------------------------------
        
           do i = is, ie
            lhl (i) = lv00 + d0_vap * pt1 (i)
            lhi (i) = li00 + dc_ice * pt1 (i)
            lcp2 (i) = lhl (i) / cvm (i)
            icp2 (i) = lhi (i) / cvm (i)
            tcp3 (i) = lcp2 (i) + icp2 (i) * min (1., dim (tice, pt1 (i)) /48.)
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Condensation/evaporation between water vapor and cloud water.
        ! -----------------------------------------------------------------------
        
            call wqs2_vect (is, ie, pt1, den, wqsat, dq2dt)
        
            adj_fac = sat_adj0
        do i = is, ie
            dq0 = (qv (i, j) - wqsat (i)) / (1. + tcp3 (i) * dq2dt (i))
            if (dq0 > 0.) then ! whole grid - box saturated
                src (i) = min (adj_fac * dq0, max (ql_gen - ql (i, j), fac_v2l * dq0))
            else ! evaporation of ql
                ! sjl 20170703 added ql factor to prevent the situation of high ql and rh<1
                ! factor = - min (1., fac_l2v * sqrt (max (0., ql (i, j)) / 1.e-5) * 10. * (1. - qv (i, j) / wqsat (i)))
                ! factor = - fac_l2v
                ! factor = - 1
                factor = - min (1., fac_l2v * 10. * (1. - qv (i, j) / wqsat (i))) ! the rh dependent factor = 1 at 90%
                src (i) = - min (ql (i, j), factor * dq0)
            endif
            qv (i, j) = qv (i, j) - src (i)
#ifdef MULTI_GASES
            qvi(i,j,1,1) = qv (i, j)
#endif
            ql (i, j) = ql (i, j) + src (i)
            q_liq (i) = q_liq (i) + src (i)
            cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
            pt1 (i) = pt1 (i) + src (i) * lhl (i) / cvm (i)
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Update latend heat coefficient.
        ! -----------------------------------------------------------------------
        
        do i = is, ie
            lhl (i) = lv00 + d0_vap * pt1 (i)
            lhi (i) = li00 + dc_ice * pt1 (i)
            lcp2 (i) = lhl (i) / cvm (i)
            icp2 (i) = lhi (i) / cvm (i)
            tcp3 (i) = lcp2 (i) + icp2 (i) * min (1., dim (tice, pt1 (i)) / 48.)
        enddo
        
        if (last_step) then
            
            ! -----------------------------------------------------------------------
            !> - condensation/evaporation between water vapor and cloud water, last time step
            !! enforce upper (no super_sat) & lower (critical rh) bounds.
            ! final iteration:
            ! -----------------------------------------------------------------------
            
            call wqs2_vect (is, ie, pt1, den, wqsat, dq2dt)
            
            do i = is, ie
                dq0 = (qv (i, j) - wqsat (i)) / (1. + tcp3 (i) * dq2dt (i))
                if (dq0 > 0.) then ! remove super - saturation, prevent super saturation over water
                    src (i) = dq0
                else ! evaporation of ql
                    ! factor = - min (1., fac_l2v * sqrt (max (0., ql (i, j)) / 1.e-5) * 10. * (1. - qv (i, j) / wqsat (i))) ! the rh dependent factor = 1 at 90%
                    ! factor = - fac_l2v
                    ! factor = - 1
                    factor = - min (1., fac_l2v * 10. * (1. - qv (i, j) / wqsat (i))) ! the rh dependent factor = 1 at 90%
                    src (i) = - min (ql (i, j), factor * dq0)
                endif
                adj_fac = 1.
                qv (i, j) = qv (i, j) - src (i)
#ifdef MULTI_GASES
                qvi(i,j,1,1) = qv(i,j)
#endif
                ql (i, j) = ql (i, j) + src (i)
                q_liq (i) = q_liq (i) + src (i)
                cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
                pt1 (i) = pt1 (i) + src (i) * lhl (i) / cvm (i)
            enddo
            
            ! -----------------------------------------------------------------------
            !> - Update latend heat coefficient.
            ! -----------------------------------------------------------------------
            
            do i = is, ie
                lhl (i) = lv00 + d0_vap * pt1 (i)
                lhi (i) = li00 + dc_ice * pt1 (i)
                lcp2 (i) = lhl (i) / cvm (i)
                icp2 (i) = lhi (i) / cvm (i)
            enddo
            
        endif
        
        ! -----------------------------------------------------------------------
        !> - Homogeneous freezing of cloud water to cloud ice.
        ! -----------------------------------------------------------------------
        
        do i = is, ie
            dtmp = t_wfr - pt1 (i) ! [ - 40, - 48]
            if (ql (i, j) > 0. .and. dtmp > 0.) then
                sink (i) = min (ql (i, j), ql (i, j) * dtmp * 0.125, dtmp / icp2 (i))
                ql (i, j) = ql (i, j) - sink (i)
                qi (i, j) = qi (i, j) + sink (i)
                q_liq (i) = q_liq (i) - sink (i)
                q_sol (i) = q_sol (i) + sink (i)
                cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
                pt1 (i) = pt1 (i) + sink (i) * lhi (i) / cvm (i)
            endif
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Update latend heat coefficient.
        ! -----------------------------------------------------------------------
        
        do i = is, ie
            lhi (i) = li00 + dc_ice * pt1 (i)
            icp2 (i) = lhi (i) / cvm (i)
        enddo
        
        ! -----------------------------------------------------------------------
        !> - bigg mechanism (heterogeneous freezing of cloud water to cloud ice).
        ! -----------------------------------------------------------------------
        
        do i = is, ie
            tc = tice0 - pt1 (i)
            if (ql (i, j) > 0.0 .and. tc > 0.) then
                sink (i) = 3.3333e-10 * dt_bigg * (exp (0.66 * tc) - 1.) * den (i) * ql (i, j) ** 2
                sink (i) = min (ql (i, j), tc / icp2 (i), sink (i))
                ql (i, j) = ql (i, j) - sink (i)
                qi (i, j) = qi (i, j) + sink (i)
                q_liq (i) = q_liq (i) - sink (i)
                q_sol (i) = q_sol (i) + sink (i)
                cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
                pt1 (i) = pt1 (i) + sink (i) * lhi (i) / cvm (i)
            endif
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Update latend heat coefficient.
        ! -----------------------------------------------------------------------
        
        do i = is, ie
            lhi (i) = li00 + dc_ice * pt1 (i)
            icp2 (i) = lhi (i) / cvm (i)
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Freezing of rain to graupel.
        ! -----------------------------------------------------------------------
        
            do i = is, ie
            dtmp = (tice - 0.1) - pt1 (i)
            if (qr (i, j) > 1.e-7 .and. dtmp > 0.) then
                tmp = min (1., (dtmp * 0.025) ** 2) * qr (i, j) ! no limit on freezing below - 40 deg c
                sink (i) = min (tmp, fac_r2g * dtmp / icp2 (i))
                qr (i, j) = qr (i, j) - sink (i)
                qg (i, j) = qg (i, j) + sink (i)
                q_liq (i) = q_liq (i) - sink (i)
                q_sol (i) = q_sol (i) + sink (i)
                cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
                pt1 (i) = pt1 (i) + sink (i) * lhi (i) / cvm (i)
            endif
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Update latend heat coefficient.
        ! -----------------------------------------------------------------------
        
            do i = is, ie
            lhi (i) = li00 + dc_ice * pt1 (i)
            icp2 (i) = lhi (i) / cvm (i)
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Melting of snow to rain or cloud water.
        ! -----------------------------------------------------------------------
        
            do i = is, ie
            dtmp = pt1 (i) - (tice + 0.1)
            if (qs (i, j) > 1.e-7 .and. dtmp > 0.) then
                tmp = min (1., (dtmp * 0.1) ** 2) * qs (i, j) ! no limter on melting above 10 deg c
                sink (i) = min (tmp, fac_smlt * dtmp / icp2 (i))
                tmp = min (sink (i), dim (qs_mlt, ql (i, j))) ! max ql due to snow melt
                qs (i, j) = qs (i, j) - sink (i)
                ql (i, j) = ql (i, j) + tmp
                qr (i, j) = qr (i, j) + sink (i) - tmp
                ! qr (i, j) = qr (i, j) + sink (i)
                q_liq (i) = q_liq (i) + sink (i)
                q_sol (i) = q_sol (i) - sink (i)
                cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
                pt1 (i) = pt1 (i) - sink (i) * lhi (i) / cvm (i)
            endif
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Autoconversion from cloud water to rain.
        ! -----------------------------------------------------------------------
        
            do i = is, ie
            if (ql (i, j) > ql0_max) then
                sink (i) = fac_l2r * (ql (i, j) - ql0_max)
                qr (i, j) = qr (i, j) + sink (i)
                ql (i, j) = ql (i, j) - sink (i)
            endif
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Update latend heat coefficient.
        ! -----------------------------------------------------------------------
        
            do i = is, ie
            lhi (i) = li00 + dc_ice * pt1 (i)
            lhl (i) = lv00 + d0_vap * pt1 (i)
            lcp2 (i) = lhl (i) / cvm (i)
            icp2 (i) = lhi (i) / cvm (i)
            tcp2 (i) = lcp2 (i) + icp2 (i)
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Sublimation/deposition between water vapor and cloud ice.
        ! -----------------------------------------------------------------------
        
            do i = is, ie
            src (i) = 0.
            if (pt1 (i) < t_sub) then ! too cold to be accurate; freeze qv as a fix
                src (i) = dim (qv (i, j), 1.e-6)
            elseif (pt1 (i) < tice0) then
                qsi = iqs2 (pt1 (i), den (i), dqsdt)
                dq = qv (i, j) - qsi
                sink (i) = adj_fac * dq / (1. + tcp2 (i) * dqsdt)
                if (qi (i, j) > 1.e-8) then
                    pidep = sdt * dq * 349138.78 * exp (0.875 * log (qi (i, j) * den (i))) &
                         / (qsi * den (i) * lat2 / (0.0243 * rvgas * pt1 (i) ** 2) + 4.42478e4)
                else
                    pidep = 0.
                endif
                if (dq > 0.) then ! vapor - > ice
                    tmp = tice - pt1 (i)
                    qi_crt = qi_gen * min (qi_lim, 0.1 * tmp) / den (i)
                    src (i) = min (sink (i), max (qi_crt - qi (i, j), pidep), tmp / tcp2 (i))
                else
                    pidep = pidep * min (1., dim (pt1 (i), t_sub) * 0.2)
                    src (i) = max (pidep, sink (i), - qi (i, j))
                endif
            endif
            qv (i, j) = qv (i, j) - src (i)
#ifdef MULTI_GASES
            qvi(i,j,1,1) = qv(i,j)
#endif
            qi (i, j) = qi (i, j) + src (i)
            q_sol (i) = q_sol (i) + src (i)
            cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
            pt1 (i) = pt1 (i) + src (i) * (lhl (i) + lhi (i)) / cvm (i)
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Virtual temperature updated.
        ! -----------------------------------------------------------------------
        
            do i = is, ie
#ifdef USE_COND
            q_con (i, j) = q_liq (i) + q_sol (i)
#ifdef MULTI_GASES
            pt (i, j) = pt1 (i) * virq_qpz(qvi(i,j,1,1:num_gas),q_con(i,j))
#else
            tmp = 1. + zvir * qv (i, j)
            pt (i, j) = pt1 (i) * tmp * (1. - q_con (i, j))
#endif
            tmp = rdgas * tmp
            cappa (i, j) = tmp / (tmp + cvm (i))
#else
#ifdef MULTI_GASES
            q_con (i, j) = q_liq (i) + q_sol (i)
            pt (i, j) = pt1 (i) * virq_qpz(qvi(i,j,1,1:num_gas),q_con(i,j)) * (1. - q_con(i,j))
#else
            pt (i, j) = pt1 (i) * (1. + zvir * qv (i, j))
#endif
#endif
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Fix negative graupel with available cloud ice.
        ! -----------------------------------------------------------------------
        
            do i = is, ie
            if (qg (i, j) < 0.) then
                tmp = min (- qg (i, j), max (0., qi (i, j)))
                qg (i, j) = qg (i, j) + tmp
                qi (i, j) = qi (i, j) - tmp
            endif
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Autoconversion from cloud ice to snow.
        ! -----------------------------------------------------------------------
        
            do i = is, ie
            qim = qi0_max / den (i)
            if (qi (i, j) > qim) then
                sink (i) = fac_i2s * (qi (i, j) - qim)
                qi (i, j) = qi (i, j) - sink (i)
                qs (i, j) = qs (i, j) + sink (i)
            endif
        enddo
        
            if (out_dt) then
            do i = is, ie
                dtdt (i, j) = dtdt (i, j) + pt1 (i) - t0 (i)
            enddo
        endif
        
        ! -----------------------------------------------------------------------
        !> - Fix energy conservation.
        ! -----------------------------------------------------------------------
        
            if (consv_te) then
            do i = is, ie
                if (hydrostatic) then
#ifdef MULTI_GASES
                    c_air = cp_air * vicpqd_qpz(qvi(i,j,1,1:num_gas),qpz(i))
#endif
                    te0 (i, j) = dp (i, j) * (te0 (i, j) + c_air * pt1 (i))
                else
#ifdef USE_COND
                    te0 (i, j) = dp (i, j) * (te0 (i, j) + cvm (i) * pt1 (i))
#else
#ifdef MULTI_GASES
                    c_air = cv_air * vicvqd_qpz(qvi(i,j,1,1:num_gas),qpz(i))
#endif
                    te0 (i, j) = dp (i, j) * (te0 (i, j) + c_air * pt1 (i))
#endif
                endif
            enddo
        endif
        
        ! -----------------------------------------------------------------------
        !> - Update latend heat coefficient.
        ! -----------------------------------------------------------------------
        
            do i = is, ie
            lhi (i) = li00 + dc_ice * pt1 (i)
            lhl (i) = lv00 + d0_vap * pt1 (i)
            cvm (i) = mc_air (i) + (qv (i, j) + q_liq (i) + q_sol (i)) * c_vap
            lcp2 (i) = lhl (i) / cvm (i)
            icp2 (i) = lhi (i) / cvm (i)
        enddo
        
        ! -----------------------------------------------------------------------
        !> - Compute cloud fraction.
        ! -----------------------------------------------------------------------
        
        if (do_qa .and. last_step) then
            
            ! -----------------------------------------------------------------------
            !>  - If it is the last step, combine water species.
            ! -----------------------------------------------------------------------
            
            if (rad_snow) then
                if (rad_graupel) then
                    do i = is, ie
                        q_sol (i) = qi (i, j) + qs (i, j) + qg (i, j)
                    enddo
                else
                    do i = is, ie
                        q_sol (i) = qi (i, j) + qs (i, j)
                    enddo
                endif
            else
                do i = is, ie
                    q_sol (i) = qi (i, j)
                enddo
            endif
            if (rad_rain) then
                do i = is, ie
                    q_liq (i) = ql (i, j) + qr (i, j)
                enddo
            else
                do i = is, ie
                    q_liq (i) = ql (i, j)
                enddo
            endif
            do i = is, ie
                q_cond (i) = q_sol (i) + q_liq (i)
            enddo
            
            ! -----------------------------------------------------------------------
            !>  - Use the "liquid - frozen water temperature" (tin) to compute saturated specific humidity.
            ! -----------------------------------------------------------------------
            
            do i = is, ie
                
                if(tintqs) then 
                  tin = pt1(i)
                else 
                  tin = pt1 (i) - (lcp2 (i) * q_cond (i) + icp2 (i) * q_sol (i)) ! minimum temperature
                ! tin = pt1 (i) - ((lv00 + d0_vap * pt1 (i)) * q_cond (i) + &
                ! (li00 + dc_ice * pt1 (i)) * q_sol (i)) / (mc_air (i) + qpz (i) * c_vap)
                endif 
                
                ! -----------------------------------------------------------------------
                ! determine saturated specific humidity
                ! -----------------------------------------------------------------------
                
                if (tin <= t_wfr) then
                    ! ice phase:
                    qstar (i) = iqs1 (tin, den (i))
                elseif (tin >= tice) then
                    ! liquid phase:
                    qstar (i) = wqs1 (tin, den (i))
                else
                    ! mixed phase:
                    qsi = iqs1 (tin, den (i))
                    qsw = wqs1 (tin, den (i))
                    if (q_cond (i) > 1.e-6) then
                        rqi = q_sol (i) / q_cond (i)
                    else
                        ! mostly liquid water clouds at initial cloud development stage
                        rqi = ((tice - tin) / (tice - t_wfr))
                    endif
                    qstar (i) = rqi * qsi + (1. - rqi) * qsw
                endif
                !>   - higher than 10 m is considered "land" and will have higher subgrid variability
                dw = dw_ocean + (dw_land - dw_ocean) * min (1., abs (hs (i, j)) / (10. * grav))
                !>   - "scale - aware" subgrid variability: 100 - km as the base
                hvar (i) = min (0.2, max (0.01, dw * sqrt (sqrt (area (i, j)) / 100.e3)))
                
                ! -----------------------------------------------------------------------
                !>   - calculate partial cloudiness by pdf;
                !! assuming subgrid linear distribution in horizontal; this is effectively a smoother for the
                !! binary cloud scheme; qa = 0.5 if qstar (i) == qpz
                ! -----------------------------------------------------------------------
                
                rh = qpz (i) / qstar (i)
                
                ! -----------------------------------------------------------------------
                ! icloud_f = 0: bug - fixed
                ! icloud_f = 1: old fvgfs gfdl) mp implementation
                ! icloud_f = 2: binary cloud scheme (0 / 1)
                ! -----------------------------------------------------------------------
                
                if (rh > 0.75 .and. qpz (i) > 1.e-8) then
                    dq = hvar (i) * qpz (i)
                    q_plus = qpz (i) + dq
                    q_minus = qpz (i) - dq
                    if (icloud_f == 2) then
                        if (qpz (i) > qstar (i)) then
                            qa (i, j) = 1.
                        elseif (qstar (i) < q_plus .and. q_cond (i) > 1.e-8) then
                            qa (i, j) = ((q_plus - qstar (i)) / dq) ** 2
                            qa (i, j) = min (1., qa (i, j))
                        else
                            qa (i, j) = 0.
                        endif
                    else
                        if (qstar (i) < q_minus) then
                            qa (i, j) = 1.
                        else
                            if (qstar (i) < q_plus) then
                                if (icloud_f == 0) then
                                    qa (i, j) = (q_plus - qstar (i)) / (dq + dq)
                                else
                                    qa (i, j) = (q_plus - qstar (i)) / (2. * dq * (1. - q_cond (i)))
                                endif
                            else
                                qa (i, j) = 0.
                            endif
                            ! impose minimum cloudiness if substantial q_cond (i) exist
                            if (q_cond (i) > 1.e-8) then
                                qa (i, j) = max (cld_min, qa (i, j))
                            endif
                            qa (i, j) = min (1., qa (i, j))
                        endif
                    endif
                else
                    qa (i, j) = 0.
                endif
                
            enddo
            
        endif
        
    enddo ! end j loop
    
end subroutine fv_sat_adj_work
!> @}

! =======================================================================
!>\ingroup fast_sat_adj
!>\brief the function 'wqs1' computes the 
!! saturated specific humidity for table ii.
! =======================================================================
real(kind=kind_dyn) function wqs1 (ta, den)

    implicit none

    ! pure water phase; universal dry / moist formular using air density
    ! input "den" can be either dry or moist air density

    real(kind=kind_dyn), intent (in) :: ta, den

    real(kind=kind_dyn) :: es, ap1, tmin

    integer :: it

    tmin = tice - 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

! =======================================================================
!>\ingroup fast_sat_adj
!>\brief the function 'wqs1' computes the  saturated specific humidity 
!! for table iii
! =======================================================================
real(kind=kind_dyn) function iqs1 (ta, den)

    implicit none

    ! water - ice phase; universal dry / moist formular using air density
    ! input "den" can be either dry or moist air density

    real(kind=kind_dyn), intent (in) :: ta, den

    real(kind=kind_dyn) :: es, ap1, tmin

    integer :: it

    tmin = tice - 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

! =======================================================================
!>\ingroup fast_sat_adj
!>\brief The function 'wqs2'computes the gradient of saturated specific 
!! humidity for table ii
! =======================================================================
real(kind=kind_dyn) 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(kind=kind_dyn), intent (in) :: ta, den

    real(kind=kind_dyn), intent (out) :: dqdt

    real(kind=kind_dyn) :: es, ap1, tmin

    integer :: it

    tmin = tice - 160.
    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

! =======================================================================
!>\ingroup fast_sat_adj
!>\brief The function wqs2_vect computes the gradient of saturated 
!! specific humidity for table ii.
!! It is the same as "wqs2", but written as vector function.
! =======================================================================
subroutine wqs2_vect (is, ie, ta, den, wqsat, dqdt)

    implicit none

    ! pure water phase; universal dry / moist formular using air density
    ! input "den" can be either dry or moist air density

    integer, intent (in) :: is, ie

    real(kind=kind_dyn), intent (in), dimension (is:ie) :: ta, den

    real(kind=kind_dyn), intent (out), dimension (is:ie) :: wqsat, dqdt

    real(kind=kind_dyn) :: es, ap1, tmin

    integer :: i, it

    tmin = tice - 160.

    do i = is, ie
        ap1 = 10. * dim (ta (i), tmin) + 1.
        ap1 = min (2621., ap1)
        it = ap1
        es = tablew (it) + (ap1 - it) * desw (it)
        wqsat (i) = es / (rvgas * ta (i) * den (i))
        it = ap1 - 0.5
        ! finite diff, del_t = 0.1:
        dqdt (i) = 10. * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) / (rvgas * ta (i) * den (i))
    enddo

end subroutine wqs2_vect

! =======================================================================
!>\ingroup fast_sat_adj
!>\brief The function 'iqs2' computes the gradient of saturated specific 
!! humidity for table iii.
! =======================================================================
real(kind=kind_dyn) 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(kind=kind_dyn), intent (in) :: ta, den

    real(kind=kind_dyn), intent (out) :: dqdt

    real(kind=kind_dyn) :: es, ap1, tmin

    integer :: it

    tmin = tice - 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
    ! finite diff, del_t = 0.1:
    dqdt = 10. * (des2 (it) + (ap1 - it) * (des2 (it + 1) - des2 (it))) / (rvgas * ta * den)

end function iqs2

! =======================================================================
!>\ingroup fast_sat_adj
!! saturation water vapor pressure table i
! 3 - phase table
! =======================================================================

subroutine qs_table (n)
    
    implicit none
    
    integer, intent (in) :: n
    real (kind_grid) :: delt = 0.1
    real (kind_grid) :: tmin, tem, esh20
    real (kind_grid) :: wice, wh2o, fac0, fac1, fac2
    real (kind_grid) :: esupc (200)
    integer :: i
    
    tmin = tice - 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 - tice) / (tem * tice)
        fac1 = fac0 * li2
        fac2 = (d2ice * log (tem / tice) + 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 - tice) / (tem * tice)
        fac1 = fac0 * lv0
        fac2 = (dc_vap * log (tem / tice) + 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 * (tice - tem)
        wh2o = 0.05 * (tem - 253.16)
        table (i + 1400) = wice * table (i + 1400) + wh2o * esupc (i)
    enddo
    
end subroutine qs_table

! =======================================================================
!>\ingroup fast_sat_adj
!! saturation water vapor pressure table ii.
! 1 - phase table
! =======================================================================

subroutine qs_tablew (n)
    
    implicit none
    
    integer, intent (in) :: n
    real (kind_grid) :: delt = 0.1
    real (kind_grid) :: tmin, tem, fac0, fac1, fac2
    integer :: i
    
    tmin = tice - 160.
    
    ! -----------------------------------------------------------------------
    ! compute es over water
    ! -----------------------------------------------------------------------
    
    do i = 1, n
        tem = tmin + delt * real (i - 1)
        fac0 = (tem - tice) / (tem * tice)
        fac1 = fac0 * lv0
        fac2 = (dc_vap * log (tem / tice) + fac1) / rvgas
        tablew (i) = e00 * exp (fac2)
    enddo
    
end subroutine qs_tablew

! =======================================================================
!>\ingroup fast_sat_adj
!! saturation water vapor pressure table iii.
! 2 - phase table
! =======================================================================

subroutine qs_table2 (n)
    
    implicit none
    
    integer, intent (in) :: n
    real (kind_grid) :: delt = 0.1
    real (kind_grid) :: tmin, tem0, tem1, fac0, fac1, fac2
    integer :: i, i0, i1
    
    tmin = tice - 160.
    
    do i = 1, n
        tem0 = tmin + delt * real (i - 1)
        fac0 = (tem0 - tice) / (tem0 * tice)
        if (i <= 1600) then
            ! -----------------------------------------------------------------------
            ! compute es over ice between - 160 deg c and 0 deg c.
            ! -----------------------------------------------------------------------
            fac1 = fac0 * li2
            fac2 = (d2ice * log (tem0 / tice) + fac1) / rvgas
        else
            ! -----------------------------------------------------------------------
            ! compute es over water between 0 deg c and 102 deg c.
            ! -----------------------------------------------------------------------
            fac1 = fac0 * lv0
            fac2 = (dc_vap * log (tem0 / tice) + 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

end module fv_sat_adj
!> @}