source: LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.f90 @ 5408

MODULE lmdz_lscp_condensation
!----------------------------------------------------------------
! Module for condensation of clouds routines
! that are called in LSCP


IMPLICIT NONE

CONTAINS

!**********************************************************************************
SUBROUTINE condensation_lognormal( &
      klon, temp, qtot, qsat, gamma_cond, ratqs, &
      keepgoing, cldfra, qincld, qvc)

!----------------------------------------------------------------------
! This subroutine calculates the formation of clouds, using a
! statistical cloud scheme. It uses a generalised lognormal distribution
! of total water in the gridbox
! See Bony and Emanuel, 2001
!----------------------------------------------------------------------

USE lmdz_lscp_ini, ONLY: eps

IMPLICIT NONE

!
! Input
!
INTEGER,  INTENT(IN)                     :: klon          ! number of horizontal grid points
!
REAL,     INTENT(IN)   , DIMENSION(klon) :: temp          ! temperature [K]
REAL,     INTENT(IN)   , DIMENSION(klon) :: qtot          ! total specific humidity (without precip) [kg/kg]
REAL,     INTENT(IN)   , DIMENSION(klon) :: qsat          ! saturation specific humidity [kg/kg]
REAL,     INTENT(IN)   , DIMENSION(klon) :: gamma_cond    ! condensation threshold w.r.t. qsat [-]
REAL,     INTENT(IN)   , DIMENSION(klon) :: ratqs         ! ratio between the variance of the total water distribution and its average [-]
LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: keepgoing     ! .TRUE. if a new condensation loop should be computed
!
!  Output
!  NB. those are in INOUT because of the convergence loop on temperature
!      (in some cases, the values are not re-computed) but the values
!      are never used explicitely
!
REAL,     INTENT(INOUT), DIMENSION(klon) :: cldfra   ! cloud fraction [-]
REAL,     INTENT(INOUT), DIMENSION(klon) :: qincld   ! cloud-mean in-cloud total specific water [kg/kg]
REAL,     INTENT(INOUT), DIMENSION(klon) :: qvc      ! gridbox-mean vapor in the cloud [kg/kg]
!
! Local
!
INTEGER :: i
REAL :: pdf_std, pdf_k, pdf_delta
REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
!
!--Loop on klon
!
DO i = 1, klon

  !--If a new calculation of the condensation is needed,
  !--i.e., temperature has not yet converged (or the cloud is
  !--formed elsewhere)
  IF (keepgoing(i)) THEN

    pdf_std   = ratqs(i) * qtot(i)
    pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2 ) )
    pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
    pdf_a     = pdf_delta / ( pdf_k * SQRT(2.) )
    pdf_b     = pdf_k / (2. * SQRT(2.))
    pdf_e1    = pdf_a - pdf_b
    pdf_e1    = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
    pdf_e1    = 1. - ERF(pdf_e1)
    pdf_e2    = pdf_a + pdf_b
    pdf_e2    = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
    pdf_e2    = 1. - ERF(pdf_e2)


    IF ( pdf_e1 .LT. eps ) THEN
      cldfra(i) = 0.
      qincld(i) = qsat(i)
      !--AB grid-mean vapor in the cloud - we assume saturation adjustment
      qvc(i) = 0.
    ELSE
      cldfra(i) = 0.5 * pdf_e1
      qincld(i) = qtot(i) * pdf_e2 / pdf_e1
      !--AB grid-mean vapor in the cloud - we assume saturation adjustment
      qvc(i) = qsat(i) * cldfra(i)
    ENDIF

  ENDIF   ! end keepgoing

ENDDO     ! end loop on i
END SUBROUTINE condensation_lognormal
!**********************************************************************************

!**********************************************************************************
SUBROUTINE condensation_ice_supersat( &
      klon, dtime, missing_val, pplay, paprsdn, paprsup, &
      cf_seri, rvc_seri, ql_seri, qi_seri, shear, pbl_eps, cell_area, &
      temp, qtot, qsat, gamma_cond, ratqs, keepgoing, &
      cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, &
      dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, &
      Tcontr, qcontr, qcontr2, fcontrN, fcontrP, flight_dist, flight_h2o, &
      dcf_avi, dqi_avi, dqvc_avi)

!----------------------------------------------------------------------
! This subroutine calculates the formation, evolution and dissipation
! of clouds, using a process-oriented treatment of the cloud properties
! (cloud fraction, vapor in the cloud, condensed water in the cloud).
! It allows for ice supersaturation in cold regions, in clear sky.
! If ok_unadjusted_clouds, it also allows for sub- and supersaturated
! cloud water vapors.
! It also allows for the formation and evolution of condensation trails
! (contrails) from aviation.
! Authors: Audran Borella, Etienne Vignon, Olivier Boucher
! April 2024
!----------------------------------------------------------------------

USE lmdz_lscp_tools,      ONLY: calc_qsat_ecmwf, calc_gammasat, GAMMAINC
USE lmdz_lscp_ini,        ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI, EPS_W
USE lmdz_lscp_ini,        ONLY: eps, temp_nowater, ok_weibull_warm_clouds
USE lmdz_lscp_ini,        ONLY: ok_unadjusted_clouds, iflag_cloud_sublim_pdf
USE lmdz_lscp_ini,        ONLY: lunout

USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, std_subl_pdf_lscp, &
                                mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, &
                                std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, &
                                coef_mixing_lscp, coef_shear_lscp, & 
                                chi_mixing_lscp, rho_ice

IMPLICIT NONE

!
! Input
!
INTEGER,  INTENT(IN)                     :: klon          ! number of horizontal grid points
REAL,     INTENT(IN)                     :: dtime         ! time step [s]
REAL,     INTENT(IN)                     :: missing_val   ! missing value for output
!
REAL,     INTENT(IN)   , DIMENSION(klon) :: pplay         ! layer pressure [Pa]
REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsdn       ! pressure at the lower interface [Pa]
REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsup       ! pressure at the upper interface [Pa]
REAL,     INTENT(IN)   , DIMENSION(klon) :: cf_seri       ! cloud fraction [-]
REAL,     INTENT(IN)   , DIMENSION(klon) :: rvc_seri      ! gridbox-mean water vapor in cloud [kg/kg]
REAL,     INTENT(IN)   , DIMENSION(klon) :: ql_seri       ! specific liquid water content [kg/kg]
REAL,     INTENT(IN)   , DIMENSION(klon) :: qi_seri       ! specific ice water content [kg/kg]
REAL,     INTENT(IN)   , DIMENSION(klon) :: shear         ! vertical shear [s-1]
REAL,     INTENT(IN)   , DIMENSION(klon) :: pbl_eps       ! 
REAL,     INTENT(IN)   , DIMENSION(klon) :: cell_area     ! 
REAL,     INTENT(IN)   , DIMENSION(klon) :: temp          ! temperature [K]
REAL,     INTENT(IN)   , DIMENSION(klon) :: qtot          ! total specific humidity (without precip) [kg/kg]
REAL,     INTENT(IN)   , DIMENSION(klon) :: qsat          ! saturation specific humidity [kg/kg]
REAL,     INTENT(IN)   , DIMENSION(klon) :: gamma_cond    ! condensation threshold w.r.t. qsat [-]
REAL,     INTENT(IN)   , DIMENSION(klon) :: ratqs         ! ratio between the variance of the total water distribution and its average [-]
LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: keepgoing     ! .TRUE. if a new condensation loop should be computed
!
! Input for aviation
!
REAL,     INTENT(IN),    DIMENSION(klon) :: flight_dist   ! 
REAL,     INTENT(IN),    DIMENSION(klon) :: flight_h2o    ! 
!
!  Output
!  NB. cldfra and qincld should be outputed as cf_seri and qi_seri,
!      or as tendencies (maybe in the future)
!  NB. those are in INOUT because of the convergence loop on temperature
!      (in some cases, the values are not re-computed) but the values
!      are never used explicitely
!
REAL,     INTENT(INOUT), DIMENSION(klon) :: cldfra   ! cloud fraction [-]
REAL,     INTENT(INOUT), DIMENSION(klon) :: qincld   ! cloud-mean in-cloud total specific water [kg/kg]
REAL,     INTENT(INOUT), DIMENSION(klon) :: qvc      ! gridbox-mean vapor in the cloud [kg/kg]
REAL,     INTENT(INOUT), DIMENSION(klon) :: issrfra  ! ISSR fraction [-]
REAL,     INTENT(INOUT), DIMENSION(klon) :: qissr    ! gridbox-mean ISSR specific water [kg/kg]
!
!  Diagnostics for condensation and ice supersaturation
!  NB. idem for the INOUT
!
REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_sub  ! cloud fraction tendency because of sublimation [s-1]
REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_con  ! cloud fraction tendency because of condensation [s-1]
REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_mix  ! cloud fraction tendency because of cloud mixing [s-1]
REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_adj  ! specific ice content tendency because of temperature adjustment [kg/kg/s]
REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_sub  ! specific ice content tendency because of sublimation [kg/kg/s]
REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_con  ! specific ice content tendency because of condensation [kg/kg/s]
REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_mix  ! specific ice content tendency because of cloud mixing [kg/kg/s]
REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_adj ! specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_sub ! specific cloud water vapor tendency because of sublimation [kg/kg/s]
REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_con ! specific cloud water vapor tendency because of condensation [kg/kg/s]
REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
!
!  Diagnostics for aviation
!  NB. idem for the INOUT
!
REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcontr   ! critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) [K]
REAL,     INTENT(INOUT), DIMENSION(klon) :: qcontr   ! 
REAL,     INTENT(INOUT), DIMENSION(klon) :: qcontr2  ! 
REAL,     INTENT(INOUT), DIMENSION(klon) :: fcontrN  ! 
REAL,     INTENT(INOUT), DIMENSION(klon) :: fcontrP  ! 
REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_avi  ! cloud fraction tendency because of aviation [s-1]
REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_avi  ! specific ice content tendency because of aviation [kg/kg/s]
REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_avi ! specific cloud water vapor tendency because of aviation [kg/kg/s]
!
! Local
!
INTEGER :: i
LOGICAL :: ok_warm_cloud
REAL, DIMENSION(klon) :: qcld, qzero
!
! for lognormal
REAL :: pdf_std, pdf_k, pdf_delta
REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
!
! for unadjusted clouds
REAL :: qvapincld, qvapincld_new
REAL :: qiceincld
!
! for sublimation
REAL :: pdf_alpha
REAL :: dqt_sub
!
! for condensation
REAL, DIMENSION(klon) :: qsatl, dqsatl
REAL :: clrfra, qclr, sl_clr, rhl_clr
REAL :: pdf_ratqs, pdf_skew, pdf_scale, pdf_loc
REAL :: pdf_x, pdf_y, pdf_T
REAL :: pdf_e3, pdf_e4
REAL :: cf_cond, qt_cond, dqt_con
!
! for mixing
REAL, DIMENSION(klon) :: subfra, qsub
REAL :: dqt_mix_sub, dqt_mix_issr
REAL :: dcf_mix_sub, dcf_mix_issr
REAL :: dqvc_mix_sub, dqvc_mix_issr
REAL :: dqt_mix
REAL :: a_mix, bovera, N_cld_mix, L_mix
REAL :: envfra_mix, cldfra_mix
REAL :: L_shear, shear_fra
REAL :: sigma_mix, issrfra_mix, subfra_mix, qvapinmix
!
! for cell properties
REAL :: rho, rhodz, dz
!REAL :: V_cell, M_cell
!
! for aviation and cell properties
!REAL :: dqt_avi
!REAL :: contrail_fra
!
!
!--more local variables for diagnostics
!--imported from YOMCST.h 
!--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1)
!--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air
!--values from Schumann, Meteorol Zeitschrift, 1996
!--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen
!--Qheat = 43.  /  50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen
!REAL, PARAMETER :: EiH2O=1.25  !--emission index of water vapour for kerosene (kg kg-1)
!REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1) 
!REAL, PARAMETER :: eta=0.3     !--average propulsion efficiency of the aircraft
!--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute 
!--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010.
!REAL :: Gcontr
!--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) 
!--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1)
!REAL :: qsatliqcontr


!-----------------------------------------------
! Initialisations
!-----------------------------------------------

! Ajout des émissions de H2O dues à l'aviation
! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O )  
!      = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) ) 
! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q) 
! flight_h2O is in kg H2O / s / cell
!
!IF (ok_plane_h2o) THEN 
!  q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) ) 
!ENDIF


qzero(:) = 0.

!--Calculation of qsat w.r.t. liquid
CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 1, .FALSE., qsatl, dqsatl)

!
!--Loop on klon
!
DO i = 1, klon

  !--If a new calculation of the condensation is needed,
  !--i.e., temperature has not yet converged (or the cloud is
  !--formed elsewhere)
  IF (keepgoing(i)) THEN

    !--If the temperature is higher than the threshold below which
    !--there is no liquid in the gridbox, we activate the usual scheme
    !--(generalised lognormal from Bony and Emanuel 2001)
    !--If ok_weibull_warm_clouds = .TRUE., the Weibull law is used for
    !--all clouds, and the lognormal scheme is not activated
    IF ( ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds ) THEN

      pdf_std   = ratqs(i) * qtot(i)
      pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2 ) )
      pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
      pdf_a     = pdf_delta / ( pdf_k * SQRT(2.) )
      pdf_b     = pdf_k / (2. * SQRT(2.))
      pdf_e1    = pdf_a - pdf_b
      pdf_e1    = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
      pdf_e1    = 1. - ERF(pdf_e1)
      pdf_e2    = pdf_a + pdf_b
      pdf_e2    = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
      pdf_e2    = 1. - ERF(pdf_e2)


      IF ( pdf_e1 .LT. eps ) THEN
        cldfra(i) = 0.
        qincld(i) = qsat(i)
        qvc(i) = 0.
      ELSE
        cldfra(i) = 0.5 * pdf_e1
        qincld(i) = qtot(i) * pdf_e2 / pdf_e1
        qvc(i) = qsat(i) * cldfra(i)
      ENDIF

    !--If the temperature is lower than temp_nowater, we use the new
    !--condensation scheme that allows for ice supersaturation
    ELSE

      !--Initialisation
      IF ( temp(i) .GT. temp_nowater ) THEN
        !--If the air mass is warm (liquid water can exist),
        !--all the memory is lost and the scheme becomes statistical,
        !--i.e., the sublimation and mixing processes are deactivated,
        !--and the condensation process is slightly adapted
        !--This can happen only if ok_weibull_warm_clouds = .TRUE.
        ! AB WARM CLOUD
        !cldfra(i) = 0.
        !qcld(i)   = 0.
        !qvc(i)    = 0.
        cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
        qcld(i)   = MAX(0., MIN(qtot(i), ql_seri(i) + qi_seri(i) + rvc_seri(i) * qtot(i)))
        qvc(i)    = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
        ok_warm_cloud = .TRUE.
      ELSE
        !--The following barriers ensure that the traced cloud properties
        !--are consistent. In some rare cases, i.e. the cloud water vapor
        !--can be greater than the total water in the gridbox
        cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
        qcld(i)   = MAX(0., MIN(qtot(i), rvc_seri(i) * qtot(i) + qi_seri(i)))
        qvc(i)    = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
        ok_warm_cloud = .FALSE.
      ENDIF

      dcf_sub(i)  = 0.
      dqi_sub(i)  = 0.
      dqvc_sub(i) = 0.
      dqi_adj(i)  = 0.
      dqvc_adj(i) = 0.
      dcf_con(i)  = 0.
      dqi_con(i)  = 0.
      dqvc_con(i) = 0.
      dcf_mix(i)  = 0.
      dqi_mix(i)  = 0.
      dqvc_mix(i) = 0.

      issrfra(i)  = 0.
      qissr(i)    = 0.
      subfra(i)   = 0.
      qsub(i)     = 0.

      !--Initialisation of the cell properties
      !--Dry density [kg/m3]
      rho = pplay(i) / temp(i) / RD
      !--Dry air mass [kg/m2]
      rhodz = ( paprsdn(i) - paprsup(i) ) / RG
      !--Cell thickness [m]
      dz = rhodz / rho
      !--Cell volume [m3]
      !V_cell = dz * cell_area(i)
      !--Cell dry air mass [kg]
      !M_cell = rhodz * cell_area(i)


      !-------------------------------------------------------------------
      !--    SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CLOUD    --
      !-------------------------------------------------------------------
      
      !--If there is a cloud
      IF ( cldfra(i) .GT. eps ) THEN
        
        qvapincld = qvc(i) / cldfra(i)
        qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
        
        !--If the ice water content is too low, the cloud is purely sublimated
        !--Most probably, we advected a cloud with no ice water content (possible
        !--if the entire cloud precipited for example)
        IF ( qiceincld .LT. eps ) THEN
          dcf_sub(i) = - cldfra(i)
          dqvc_sub(i) = - qvc(i)
          dqi_sub(i) = - ( qcld(i) - qvc(i) )

          cldfra(i) = 0.
          qcld(i) = 0.
          qvc(i) = 0.

        !--Else, the cloud is adjusted and sublimated
        ELSE

          !--The vapor in cloud cannot be higher than the
          !--condensation threshold
          qvapincld = MIN(qvapincld, gamma_cond(i) * qsat(i))
          qiceincld = ( qcld(i) / cldfra(i) - qvapincld )

          ! AB WARM CLOUD
          !IF ( ok_unadjusted_clouds ) THEN
          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
            CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
                                        pplay(i), dtime, qvapincld_new)
            IF ( qvapincld_new .EQ. 0. ) THEN
              !--If all the ice has been sublimated, we sublimate
              !--completely the cloud and do not activate the sublimation
              !--process
              !--Tendencies and diagnostics
              dcf_sub(i) = - cldfra(i)
              dqvc_sub(i) = - qvc(i)
              dqi_sub(i) = - ( qcld(i) - qvc(i) )

              cldfra(i) = 0.
              qcld(i) = 0.
              qvc(i) = 0.
            ENDIF
          ELSE
            !--We keep the saturation adjustment hypothesis, and the vapor in the
            !--cloud is set equal to the saturation vapor
            qvapincld_new = qsat(i)
          ENDIF ! ok_unadjusted_clouds
          
          !--Adjustment of the IWC to the new vapor in cloud
          !--(this can be either positive or negative)
          dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )
          dqi_adj(i) = - dqvc_adj(i)

          !--Add tendencies
          !--The vapor in the cloud is updated, but not qcld as it is constant
          !--through this process, as well as cldfra which is unmodified
          qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))
          

          !------------------------------------
          !--    DISSIPATION OF THE CLOUD    --
          !------------------------------------

          !--If the vapor in cloud is below vapor needed for the cloud to survive
          IF ( ( qvapincld .LT. qvapincld_new ) .OR. ( iflag_cloud_sublim_pdf .GE. 4 ) ) THEN
            !--Sublimation of the subsaturated cloud
            !--iflag_cloud_sublim_pdf selects the PDF of the ice water content
            !--to use.
            !--iflag = 1 --> uniform distribution
            !--iflag = 2 --> exponential distribution
            !--iflag = 3 --> gamma distribution (Karcher et al 2018)

            IF ( iflag_cloud_sublim_pdf .EQ. 1 ) THEN
              !--Uniform distribution starting at qvapincld
              pdf_e1 = 1. / ( 2. * qiceincld )

              dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1
              dqt_sub    = - cldfra(i) * ( qvapincld_new**2 - qvapincld**2 ) &
                                       * pdf_e1 / 2.

            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 2 ) THEN
              !--Exponential distribution starting at qvapincld
              pdf_alpha = 1. / qiceincld
              pdf_e1 = EXP( - pdf_alpha * ( qvapincld_new - qvapincld ) )

              dcf_sub(i) = - cldfra(i) * ( 1. - pdf_e1 )
              dqt_sub    = - cldfra(i) * ( ( 1. - pdf_e1 ) / pdf_alpha &
                                       + qvapincld - qvapincld_new * pdf_e1 )

            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 3 ) THEN
              !--Gamma distribution starting at qvapincld
              pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld
              pdf_y = pdf_alpha * ( qvapincld_new - qvapincld )
              pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y )
              pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y )

              dcf_sub(i) = - cldfra(i) * pdf_e1
              dqt_sub    = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 )

            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 4 ) THEN
              !--Normal distribution around qvapincld + qiceincld with width
              !--constant in the RHi space
              pdf_y = ( qvapincld_new - qvapincld - qiceincld ) &
                    / ( std_subl_pdf_lscp / 100. * qsat(i))
              pdf_e1 = 0.5 * ( 1. + ERF( pdf_y / SQRT(2.) ) )
              pdf_e2 = EXP( - pdf_y**2 / 2. ) / SQRT( 2. * RPI )

              dcf_sub(i) = - cldfra(i) * pdf_e1
              dqt_sub    = - cldfra(i) * ( ( qvapincld + qiceincld ) * pdf_e1 &
                                         - std_subl_pdf_lscp / 100. * qsat(i) * pdf_e2 )
                                                                                       
            ENDIF

            !--Tendencies and diagnostics
            dqvc_sub(i) = dqt_sub

            !--Add tendencies
            cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i))
            qcld(i)   = MAX(0., qcld(i) + dqt_sub)
            qvc(i)    = MAX(0., qvc(i) + dqvc_sub(i))

          ENDIF ! qvapincld .LT. qvapincld_new

        ENDIF   ! qiceincld .LT. eps
      ENDIF ! cldfra(i) .GT. eps


      !--------------------------------------------------------------------------
      !--    CONDENSATION AND DIAGNOTICS OF SUB- AND SUPERSATURATED REGIONS    --
      !--------------------------------------------------------------------------
      !--This section relies on a distribution of water in the clear-sky region of
      !--the mesh.

      !--If there is a clear-sky region
      IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN

        !--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
        qclr = qtot(i) - qcld(i)

        !--New PDF
        rhl_clr = qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100.

        !--Calculation of the properties of the PDF
        !--Parameterization from IAGOS observations
        !--pdf_e1 and pdf_e2 will be reused below

        !--Coefficient for standard deviation:
        !--  tuning coef * (clear sky area**0.25) * (function of temperature)
        !pdf_e1 = beta_pdf_lscp &
        !       * ( ( 1. - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) &
        !       * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
        !--  tuning coef * (cell area**0.25) * (function of temperature)
        pdf_e1 = beta_pdf_lscp * ( ( 1. - cldfra(i) ) * cell_area(i) )**0.25 &
                * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
        IF ( rhl_clr .GT. 50. ) THEN
          pdf_std = ( pdf_e1 - std100_pdf_lscp ) * ( 100. - rhl_clr ) / 50. + std100_pdf_lscp
        ELSE
          pdf_std = pdf_e1 * rhl_clr / 50.
        ENDIF
        pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
        pdf_alpha = EXP( rhl_clr / 100. ) * pdf_e3
        pdf_alpha = MIN(10., pdf_alpha)
        
        IF ( ok_warm_cloud ) THEN
          !--If the statistical scheme is activated, the standard deviation is adapted
          !--to depend on the pressure level. It is multiplied by ratqs, so that near the
          !--surface std is almost 0, and upper than about 450 hPa the std is left untouched
          pdf_std = pdf_std * ratqs(i)
        ENDIF

        pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
        pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2 ))
        pdf_loc = rhl_clr - pdf_scale * pdf_e2

        !--Calculation of the newly condensed water and fraction (pronostic)
        !--Integration of the clear sky PDF between gamma_cond*qsat and +inf
        !--NB. the calculated values are clear-sky averaged

        pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
        pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
        cf_cond = EXP( - pdf_y )
        qt_cond = ( pdf_e3 + pdf_loc * cf_cond ) * qsatl(i) / 100.

        ! AB WARM CLOUD
        !IF ( ok_warm_cloud ) THEN
        !  !--If the statistical scheme is activated, the calculated increase is equal
        !  !--to the cloud fraction, we assume saturation adjustment, and the
        !  !--condensation process stops
        !  cldfra(i) = cf_cond
        !  qcld(i) = qt_cond
        !  qvc(i) = cldfra(i) * qsat(i)

        !ELSEIF ( cf_cond .GT. eps ) THEN
        IF ( cf_cond .GT. eps ) THEN

          dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond
          dqt_con = ( 1. - cldfra(i) ) * qt_cond
          
          !--Barriers
          dcf_con(i) = MIN(dcf_con(i), 1. - cldfra(i))
          dqt_con = MIN(dqt_con, qclr)


          ! AB WARM CLOUD
          !IF ( ok_unadjusted_clouds ) THEN
          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
            !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute
            !--the new vapor qvapincld. The timestep is divided by two because we do not
            !--know when the condensation occurs
            qvapincld = gamma_cond(i) * qsat(i)
            qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
            CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
                                        pplay(i), dtime / 2., qvapincld_new)
            qvapincld = qvapincld_new
          ELSE
            !--We keep the saturation adjustment hypothesis, and the vapor in the
            !--newly formed cloud is set equal to the saturation vapor.
            qvapincld = qsat(i)
          ENDIF

          !--Tendency on cloud vapor and diagnostic
          dqvc_con(i) = qvapincld * dcf_con(i)
          dqi_con(i) = dqt_con - dqvc_con(i)

          !--Add tendencies
          cldfra(i) = MIN(1., cldfra(i) + dcf_con(i))
          qcld(i)   = MIN(qtot(i), qcld(i) + dqt_con)
          qvc(i)    = MIN(qcld(i), qvc(i) + dqvc_con(i))

        ENDIF ! ok_warm_cloud, cf_cond .GT. eps
        
        !--We then calculate the part that is greater than qsat
        !--and lower than gamma_cond * qsat, which is the ice supersaturated region
        pdf_x = qsat(i) / qsatl(i) * 100.
        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
        pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
        issrfra(i) = EXP( - pdf_y ) * ( 1. - cldfra(i) )
        qissr(i) = ( pdf_e3 * ( 1. - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.

        issrfra(i) = issrfra(i) - dcf_con(i)
        qissr(i) = qissr(i) - dqvc_con(i) - dqi_con(i)

      ENDIF ! ( 1. - cldfra(i) ) .GT. eps

      !--Calculation of the subsaturated clear sky fraction and water
      subfra(i) = 1. - cldfra(i) - issrfra(i)
      qsub(i) = qtot(i) - qcld(i) - qissr(i)


      !--------------------------------------
      !--           CLOUD MIXING           --
      !--------------------------------------
      !--This process mixes the cloud with its surroundings: the subsaturated clear sky,
      !--and the supersaturated clear sky. It is activated if the cloud is big enough,
      !--but does not cover the entire mesh.
      !
      ! AB WARM CLOUD
      !IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) &
      !        .AND. .NOT. ok_warm_cloud ) THEN
      IF ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) ) THEN

        !--Initialisation
        dcf_mix_sub   = 0.
        dqt_mix_sub   = 0.
        dqvc_mix_sub  = 0.
        dcf_mix_issr  = 0.
        dqt_mix_issr  = 0.
        dqvc_mix_issr = 0.


        !-- PART 1 - TURBULENT DIFFUSION

        !--Clouds within the mesh are assumed to be ellipses. The length of the
        !--semi-major axis is a and the length of the semi-minor axis is b.
        !--N_cld_mix is the number of clouds within the mesh, and
        !--clouds_perim is the total perimeter of the clouds within the mesh,
        !--not considering interfaces with other meshes (only the interfaces with clear
        !--sky are taken into account).
        !--
        !--The area of each cloud is A = a * b * RPI,
        !--and the perimeter of each cloud is
        !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
        !--
        !--With cell_area the area of the cell, we have:
        !-- cldfra = A * N_cld_mix / cell_area
        !-- clouds_perim = P * N_cld_mix
        !--
        !--We assume that the ratio between b and a is a function of
        !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
        !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
        !--are spherical.
        !-- b / a = bovera = MAX(0.1, cldfra)
        bovera = MAX(0.1, cldfra(i))
        !--The clouds perimeter is imposed using the formula from Morcrette 2012,
        !--based on observations.
        !-- clouds_perim_normalized = alpha * cldfra * ( 1. - cldfra ) = clouds_perim / ( norm_constant * cell_area )
        !--
        !--With all this, we have
        !-- a = SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) )
        !-- P = RPI * a * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
        !--and therefore, using the perimeter
        !-- alpha * cldfra * ( 1. - cldfra ) * norm_constant * cell_area
        !--             = N_cld_mix * RPI &
        !--             * SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) ) &
        !--             * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
        !--and finally
        N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - cldfra(i) )**2 &
                  * cell_area(i) * bovera / RPI &
                  / ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )**2
        !--where coef_mix_lscp = ( alpha * norm_constant )**2
        !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer
        !--In particular, it is 0 if cldfra = 1
        a_mix = SQRT( cell_area(i) * cldfra(i) / bovera / N_cld_mix / RPI )

        !--The time required for turbulent diffusion to homogenize a region of size
        !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
        !--We compute L_mix and assume that the cloud is mixed over this length
        L_mix = SQRT( dtime**3 * pbl_eps(i) )
        !--The mixing length cannot be greater than the semi-minor axis. In this case,
        !--the entire cloud is mixed.
        L_mix = MIN(L_mix, a_mix * bovera)

        !--The fraction of clear sky mixed is
        !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
        envfra_mix = N_cld_mix * RPI / cell_area(i) &
                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
        !--The fraction of cloudy sky mixed is
        !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
        cldfra_mix = N_cld_mix * RPI / cell_area(i) &
                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )


        !-- PART 2 - SHEARING

        !--The clouds are then sheared. We keep the shape and number
        !--assumptions from before. The clouds are sheared along their
        !--semi-major axis (a_mix), on the entire cell heigh dz.
        !--The increase in size is
        L_shear = coef_shear_lscp * shear(i) * dz * dtime
        !--therefore, the fraction of clear sky mixed is
        !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
        !--and the fraction of cloud mixed is
        !-- N_cld_mix * ( (a * b) - (a - L_shear) * b ) * RPI / 2. / cell_area
        !--(note that they are equal)
        shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i)
        !--and the environment and cloud mixed fractions are the same,
        !--which we add to the previous calculated mixed fractions.
        !--We therefore assume that the sheared clouds and the turbulent
        !--mixed clouds are different.
        envfra_mix = envfra_mix + shear_fra
        cldfra_mix = cldfra_mix + shear_fra


        !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES

        !--The environment fraction is allocated to subsaturated sky or supersaturated sky,
        !--according to the factor sigma_mix. This is computed as the ratio of the
        !--subsaturated sky fraction to the environment fraction, corrected by a factor
        !--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the
        !--supersaturated sky is favoured. Physically, this means that it is more likely
        !--to have supersaturated sky around the cloud than subsaturated sky.
        sigma_mix = subfra(i) / ( subfra(i) + chi_mixing_lscp * issrfra(i) )
        subfra_mix = MIN( sigma_mix * envfra_mix, subfra(i) )
        issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra(i) )
        cldfra_mix = MIN( cldfra_mix, cldfra(i) )

        !--First, we mix the subsaturated sky (subfra_mix) and the cloud close
        !--to this fraction (sigma_mix * cldfra_mix).
        IF ( subfra(i) .GT. eps ) THEN
          !--We compute the total humidity in the mixed air, which
          !--can be either sub- or supersaturated.
          qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &
                    + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &
                    / ( subfra_mix + cldfra_mix * sigma_mix )

          ! AB WARM CLOUD
          !IF ( ok_unadjusted_clouds ) THEN
          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
            qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * sigma_mix / cldfra(i) &
                      / ( subfra_mix + cldfra_mix * sigma_mix )
            CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
                                        pplay(i), dtime, qvapincld_new)
            IF ( qvapincld_new .EQ. 0. ) THEN
              !--If all the ice has been sublimated, we sublimate
              !--completely the cloud and do not activate the sublimation
              !--process
              !--Tendencies and diagnostics
              dcf_mix_sub = - sigma_mix * cldfra_mix
              dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
              dqvc_mix_sub = dcf_mix_sub * qvc(i) / cldfra(i)
            ELSE
              dcf_mix_sub = subfra_mix
              dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
              dqvc_mix_sub = dcf_mix_sub * qvapincld_new
            ENDIF
          ELSE
            IF ( qvapinmix .GT. qsat(i) ) THEN
              !--If the mixed air is supersaturated, we condense the subsaturated
              !--region which was mixed.
              dcf_mix_sub = subfra_mix
              dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
              dqvc_mix_sub = dcf_mix_sub * qsat(i)
            ELSE
              !--Else, we sublimate the cloud which was mixed.
              dcf_mix_sub = - sigma_mix * cldfra_mix
              dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
              dqvc_mix_sub = dcf_mix_sub * qsat(i)
            ENDIF
          ENDIF ! ok_unadjusted_clouds
        ENDIF ! subfra .GT. eps
   
        !--We then mix the supersaturated sky (issrfra_mix) and the cloud,
        !--for which the mixed air is always supersatured, therefore
        !--the cloud necessarily expands
        IF ( issrfra(i) .GT. eps ) THEN

          ! AB WARM CLOUD
          !IF ( ok_unadjusted_clouds ) THEN
          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
            qvapinmix = ( qissr(i) * issrfra_mix / issrfra(i) &
                      + qcld(i) * cldfra_mix * (1. - sigma_mix) / cldfra(i) ) &
                      / ( issrfra_mix + cldfra_mix * (1. -  sigma_mix) )
            qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * (1. - sigma_mix) / cldfra(i) &
                      / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) )
            CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
                                        pplay(i), dtime, qvapincld_new)
            dcf_mix_issr = issrfra_mix
            dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
            dqvc_mix_issr = dcf_mix_issr * qvapincld_new
          ELSE
            !--In this case, the additionnal vapor condenses
            dcf_mix_issr = issrfra_mix
            dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
            dqvc_mix_issr = dcf_mix_issr * qsat(i)
          ENDIF ! ok_unadjusted_clouds


        ENDIF ! issrfra .GT. eps

        !--Sum up the tendencies from subsaturated sky and supersaturated sky
        dcf_mix(i)  = dcf_mix_sub  + dcf_mix_issr
        dqt_mix     = dqt_mix_sub  + dqt_mix_issr
        dqvc_mix(i) = dqvc_mix_sub + dqvc_mix_issr
        dqi_mix(i)  = dqt_mix - dqvc_mix(i)

        !--Add tendencies
        issrfra(i) = MAX(0., issrfra(i) - dcf_mix_issr)
        qissr(i)   = MAX(0., qissr(i) - dqt_mix_issr)
        cldfra(i)  = MAX(0., MIN(1., cldfra(i) + dcf_mix(i)))
        qcld(i)    = MAX(0., MIN(qtot(i), qcld(i) + dqt_mix))
        qvc(i)     = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i)))

      ENDIF ! ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) )


      !----------------------------------------
      !--       CONTRAILS AND AVIATION       --
      !----------------------------------------

      !--Add a source of cirrus from aviation contrails
      !IF ( ok_plane_contrail ) THEN
      !  dcf_avi(i) = 0.
      !  dqi_avi(i) = 0.
      !  dqvc_avi(i) = 0.
      !  ! TODO implement ok_unadjusted_clouds
      !  IF ( issrfra(i) .GT. eps ) THEN
      !    contrail_fra = MIN(1., flight_m(i,k) * dtime * contrail_cross_section / V_cell)
      !    dcf_avi(i) = issrfra(i) * contrail_fra
      !    dqt_avi = dcf_avi(i) * qissr(i) / issrfra(i)
      !    dqvc_avi(i) = qsat(i) * dcf_avi(i)
      !    
      !    !--Add tendencies
      !    cldfra(i) = cldfra(i) + dcf_avi(i)
      !    issrfra(i) = issrfra(i) - dcf_avi(i)
      !    qcld(i) = qcld(i) + dqt_avi
      !    qvc(i) = qvc(i) + dqvc_avi(i)
      !    qissr(i) = qissr(i) - dqt_avi

      !    !--Diagnostics
      !    dqi_avi(i) = dqt_avi - qsat(i) * dcf_avi(i)
      !  ENDIF
      !  dcf_avi(i)  = dcf_avi(i)  / dtime
      !  dqi_avi(i)  = dqi_avi(i)  / dtime
      !  dqvc_avi(i) = dqvc_avi(i) / dtime
      !ENDIF



      !-------------------------------------------
      !--       FINAL BARRIERS AND OUTPUTS      --
      !-------------------------------------------

      IF ( cldfra(i) .LT. eps ) THEN
        !--If the cloud is too small, it is sublimated.
        cldfra(i) = 0.
        qcld(i)   = 0.
        qvc(i)    = 0.
        qincld(i) = qsat(i)
      ELSE
        qincld(i) = qcld(i) / cldfra(i)
      ENDIF ! cldfra .LT. eps

      !--Diagnostics
      dcf_sub(i)  = dcf_sub(i)  / dtime
      dcf_con(i)  = dcf_con(i)  / dtime 
      dcf_mix(i)  = dcf_mix(i)  / dtime 
      dqi_adj(i)  = dqi_adj(i)  / dtime 
      dqi_sub(i)  = dqi_sub(i)  / dtime 
      dqi_con(i)  = dqi_con(i)  / dtime 
      dqi_mix(i)  = dqi_mix(i)  / dtime 
      dqvc_adj(i) = dqvc_adj(i) / dtime
      dqvc_sub(i) = dqvc_sub(i) / dtime
      dqvc_con(i) = dqvc_con(i) / dtime
      dqvc_mix(i) = dqvc_mix(i) / dtime

    ENDIF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds

  ENDIF   ! end keepgoing

ENDDO     ! end loop on i

END SUBROUTINE condensation_ice_supersat
!**********************************************************************************

SUBROUTINE deposition_sublimation( &
    qvapincld, qiceincld, temp, qsat, pplay, dtime, &
    qvapincld_new)

USE lmdz_lscp_ini,        ONLY: RV, RLSTT, RTT, RD, EPS_W
USE lmdz_lscp_ini,        ONLY: eps
USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice

REAL,     INTENT(IN) :: qvapincld     ! 
REAL,     INTENT(IN) :: qiceincld     ! 
REAL,     INTENT(IN) :: temp          ! 
REAL,     INTENT(IN) :: qsat          ! 
REAL,     INTENT(IN) :: pplay         ! 
REAL,     INTENT(IN) :: dtime         ! time step [s]
REAL,     INTENT(OUT):: qvapincld_new ! 

!
! for unadjusted clouds
REAL :: qice_ratio
REAL :: pres_sat, rho, kappa
REAL :: air_thermal_conduct, water_vapor_diff
REAL :: iwc
REAL :: iwc_log_inf100, iwc_log_sup100
REAL :: iwc_inf100, alpha_inf100, coef_inf100
REAL :: mu_sup100, sigma_sup100, coef_sup100
REAL :: Dm_ice, rm_ice

!--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
!--hypothesis is lost, and the vapor in the cloud is purely prognostic.
!
!--The deposition equation is
!-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
!--from Lohmann et al. (2016), where
!--alpha is the deposition coefficient [-]
!--mi is the mass of one ice crystal [kg]
!--C is the capacitance of an ice crystal [m]
!--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
!--R_v is the specific gas constant for humid air [J/kg/K]
!--T is the temperature [K]
!--esi is the saturation pressure w.r.t. ice [Pa]
!--Dv is the diffusivity of water vapor [m2/s]
!--Ls is the specific latent heat of sublimation [J/kg/K]
!--ka is the thermal conductivity of dry air [J/m/s/K]
!
!--alpha is a coefficient to take into account the fact that during deposition, a water
!--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
!--coherent (with the same structure). It has no impact for sublimation.
!--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
!--during deposition, and alpha = 1. during sublimation.
!--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
!-- C = capa_cond_cirrus * rm_ice
!
!--We have qice = Nice * mi, where Nice is the ice crystal
!--number concentration per kg of moist air
!--HYPOTHESIS 1: the ice crystals are spherical, therefore
!-- mi = 4/3 * pi * rm_ice**3 * rho_ice
!--HYPOTHESIS 2: the ice crystals are monodisperse with the
!--initial radius rm_ice_0.
!--NB. this is notably different than the assumption
!--of a distributed qice in the cloud made in the sublimation process;
!--should it be consistent?
!
!--As the deposition process does not create new ice crystals,
!--and because we assume a same rm_ice value for all crystals
!--therefore the sublimation process does not destroy ice crystals
!--(or, in a limit case, it destroys all ice crystals), then
!--Nice is a constant during the sublimation/deposition process.
!-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
!
!--The deposition equation then reads:
!-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice
!-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
!--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
!--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
!-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
!--  *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 rm_ice_0**2 rho_ice )
!--and we have
!-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
!-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
!--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
!
!--This system of equations can be resolved with an exact
!--explicit numerical integration, having one variable resolved
!--explicitly, the other exactly. The exactly resolved variable is
!--the one decreasing, so it is qvc if the process is
!--condensation, qi if it is sublimation.
!
!--kappa is computed as an initialisation constant, as it depends only
!--on temperature and other pre-computed values
pres_sat = qsat / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay
!--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184
!--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4
kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat &
      * ( RV * temp / water_vapor_diff / pres_sat &
        + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) )
!--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process

!--Dry density [kg/m3]
rho = pplay / temp / RD

!--Here, the initial vapor in the cloud is qvapincld, and we compute
!--the new vapor qvapincld_new

!--rm_ice formula from McFarquhar and Heymsfield (1997)
iwc = qiceincld * rho * 1e3
iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837)
iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) )
iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) )

alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100
coef_inf100 = iwc_inf100 * alpha_inf100**3 / 120.

mu_sup100 = ( 5.2 + 0.0013 * ( temp - RTT ) ) &
          + ( 0.026 - 1.2E-3 * ( temp - RTT ) ) * iwc_log_sup100
sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp - RTT ) ) &
             + ( 0.018 - 2.1E-4 * ( temp - RTT ) ) * iwc_log_sup100
coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3. * mu_sup100 + 4.5 * sigma_sup100**2 )

Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2 ) &
       * coef_sup100 ) / ( coef_inf100 + coef_sup100 )
rm_ice = Dm_ice / 2. * 1.E-6

IF ( qvapincld .GE. qsat ) THEN
  !--If the cloud is initially supersaturated
  !--Exact explicit integration (qvc exact, qice explicit)
  qvapincld_new = qsat + ( qvapincld - qsat ) &
                * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
ELSE
  !--If the cloud is initially subsaturated
  !--Exact explicit integration (qice exact, qvc explicit)
  !--The barrier is set so that the resulting vapor in cloud
  !--cannot be greater than qsat
  !--qice_ratio is the ratio between the new ice content and
  !--the old one, it is comprised between 0 and 1
  qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2 * dtime * ( qsat - qvapincld ) )

  IF ( qice_ratio .LT. 0. ) THEN
    !--The new vapor in cloud is set to 0 so that the
    !--sublimation process does not activate
    qvapincld_new = 0.
  ELSE
    !--Else, the sublimation process is activated with the
    !--diagnosed new cloud water vapor
    !--The new vapor in the cloud is increased with the
    !--sublimated ice
    qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**1.5 )
    !--The new vapor in the cloud cannot be greater than qsat
    qvapincld_new = MIN(qvapincld_new, qsat)
  ENDIF ! qice_ratio .LT. 0.
ENDIF ! qvapincld .GT. qsat

END SUBROUTINE deposition_sublimation

END MODULE lmdz_lscp_condensation