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, Povera, 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

    !--Initialisation
    issrfra(i)  = 0.
    qissr(i)    = 0.

    !--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.

      !--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 in contact with clear sky, and can be non-integer.
        !--In particular, it is 0 if cldfra = 1.
        !--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))
        !--P / a is a function of b / a only, that we can calculate
        !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
        Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )
        !--The clouds perimeter is imposed using the formula from Morcrette 2012,
        !--based on observations.
        !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra )
        !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have:
        !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra )
        !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) )
        a_mix = Povera / coef_mixing_lscp / RPI / ( 1. - cldfra(i) ) / bovera
        !--and finally,
        !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a )
        N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - cldfra(i) ) * cell_area(i) &
                  / Povera / a_mix

        !--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


!**********************************************************************************

!**********************************************************************************
SUBROUTINE  condensation_cloudth(klon,                     &
&           temp,qt,qt_th,frac_th,zpspsk,play,thetal_th,   &
&           ratqs,sigma_qtherm,qsth,qsenv,qcloud,ctot,ctotth,ctot_vol,  &
&           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
! This routine computes the condensation of clouds in convective boundary layers
! with thermals assuming two separate distribution of the saturation deficit in 
! the thermal plumes and in the environment
! It is based on the work of Arnaud Jam (Jam et al. 2013, BLM)
! Author : Etienne Vignon (LMDZ/CNRS)
! Date: February 2025
! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam
! IMPORTANT NOTE: we assume iflag_cloudth_vert=7
!-----------------------------------------------------------------------------------

      use lmdz_lscp_ini,    only: iflag_cloudth_vert,iflag_ratqs,iflag_cloudth_vert_noratqs
      use lmdz_lscp_ini,    only: vert_alpha, vert_alpha_th ,sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin
      use lmdz_lscp_ini,    only: RTT, RG, RPI, RD, RV, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th

      IMPLICIT NONE


!------------------------------------------------------------------------------
! Declarations
!------------------------------------------------------------------------------

! INPUT/OUTPUT

      INTEGER, INTENT(IN)                         :: klon
      

      REAL, DIMENSION(klon),      INTENT(IN)      ::  temp          ! Temperature (liquid temperature) in the mesh [K] : has seen evap of precip
      REAL, DIMENSION(klon),      INTENT(IN)      ::  qt            ! total water specific humidity in the mesh [kg/kg]: has seen evap of precip
      REAL, DIMENSION(klon),      INTENT(IN)      ::  qt_th         ! total water specific humidity in thermals [kg/kg]: has not seen evap of precip
      REAL, DIMENSION(klon),      INTENT(IN)      ::  thetal_th     ! Liquid potential temperature in thermals [K]: has not seen the evap of precip
      REAL, DIMENSION(klon),      INTENT(IN)      ::  frac_th       ! Fraction of the mesh covered by thermals [0-1] 
      REAL, DIMENSION(klon),      INTENT(IN)      ::  zpspsk        ! Exner potential
      REAL, DIMENSION(klon),      INTENT(IN)      ::  play          ! Pressure of layers [Pa]
      REAL, DIMENSION(klon),      INTENT(IN)      ::  ratqs         ! Parameter that determines the width of the water distrib     [-]
      REAL, DIMENSION(klon),      INTENT(IN)      ::  sigma_qtherm  ! Parameter determining the width of the distrib in thermals   [-]
      REAL, DIMENSION(klon),      INTENT(IN)      ::  qsth          ! Saturation specific humidity in thermals
      REAL, DIMENSION(klon),      INTENT(IN)      ::  qsenv         ! Saturation specific humidity in environment
      
      REAL, DIMENSION(klon),      INTENT(INOUT)      ::  ctot         ! Cloud fraction [0-1]
      REAL, DIMENSION(klon),      INTENT(INOUT)      ::  ctotth       ! Cloud fraction [0-1] in thermals
      REAL, DIMENSION(klon),      INTENT(INOUT)      ::  ctot_vol     ! Volume cloud fraction [0-1]
      REAL, DIMENSION(klon),      INTENT(INOUT)      ::  qcloud       ! In cloud total water content [kg/kg] 
      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sth    ! mean saturation deficit in thermals
      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_senv   ! mean saturation deficit in environment
      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sigmath  ! std of saturation deficit in thermals
      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sigmaenv ! std of saturation deficit in environment


! LOCAL VARIABLES

      INTEGER itap,ind1,l,ig,iter,k 
      INTEGER iflag_topthermals, niter

      REAL qcth(klon)
      REAL qcenv(klon)
      REAL qctot(klon) 
      REAL cth(klon)
      REAL cenv(klon)   
      REAL cth_vol(klon)
      REAL cenv_vol(klon)
      REAL qt_env(klon), thetal_env(klon)
      REAL sqrtpi,sqrt2,sqrt2pi
      REAL alth,alenv,ath,aenv
      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
      REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
      REAL zdelta,qsatbef,zcor
      REAL Tbefth(klon), Tbefenv(klon)
      REAL qlbef
      REAL dqsatenv(klon), dqsatth(klon)
      REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
      REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) 
      REAL qincloud(klon)
      REAL alenvl, aenvl
      REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc


!------------------------------------------------------------------------------
! Initialisation
!------------------------------------------------------------------------------


      sqrt2pi=sqrt(2.*rpi)
      sqrt2=sqrt(2.)
      sqrtpi=sqrt(rpi)

!-------------------------------------------------------------------------------
! Thermal fraction calculation and standard deviation of the distribution
!-------------------------------------------------------------------------------  

! initialisations and calculation of temperature, humidity and saturation specific humidity

cloudth_senv(:) = 0.
cloudth_sth(:) = 0.
cloudth_sigmaenv(:) = 0.
cloudth_sigmath(:) = 0.


DO ind1=1,klon
  
   Tbefenv(ind1) = temp(ind1)
   thetal_env(ind1) = Tbefenv(ind1)/zpspsk(ind1)
   Tbefth(ind1)  = thetal_th(ind1)*zpspsk(ind1)
   qt_env(ind1)  = (qt(ind1)-frac_th(ind1)*qt_th(ind1))/(1.-frac_th(ind1)) !qt = a*qtth + (1-a)*qtenv 

ENDDO



DO ind1=1,klon


    IF (frac_th(ind1).GT.min_frac_th_cld) THEN !Thermal and environnement

! Environment:


        alenv=(RD/RV*RLVTT*qsenv(ind1))/(rd*thetal_env(ind1)**2)     
        aenv=1./(1.+(alenv*RLVTT/rcpd))                             
        senv=aenv*(qt_env(ind1)-qsenv(ind1))                           
     

! Thermals:


        alth=(RD/RV*RLVTT*qsth(ind1))/(rd*thetal_th(ind1)**2)       
        ath=1./(1.+(alth*RLVTT/rcpd))                                                         
        sth=ath*(qt_th(ind1)-qsth(ind1))                      


! Standard deviation of the distributions

           ! environment
           sigma1s_fraca = (sigma1s_factor**0.5)*(frac_th(ind1)**sigma1s_power) / &
           &                (1-frac_th(ind1))*((sth-senv)**2)**0.5

           IF (cloudth_ratqsmin>0.) THEN
             sigma1s_ratqs = cloudth_ratqsmin*qt(ind1)
           ELSE
             sigma1s_ratqs = ratqs(ind1)*qt(ind1)
           ENDIF 
           sigma1s = sigma1s_fraca + sigma1s_ratqs

           IF (iflag_ratqs.eq.10.or.iflag_ratqs.eq.11) then
              sigma1s = ratqs(ind1)*qt(ind1)*aenv
           ENDIF

           ! thermals
           sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((frac_th(ind1)+0.02)**sigma2s_power))+0.002*qt_th(ind1)

          IF (iflag_ratqs.eq.10.and.sigma_qtherm(ind1).ne.0) then
             sigma2s = sigma_qtherm(ind1)*ath
          ENDIF

 
! surface cloud fraction 

           deltasenv=aenv*vert_alpha*sigma1s
           deltasth=ath*vert_alpha_th*sigma2s

           xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
           xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
           exp_xenv1 = exp(-1.*xenv1**2)
           exp_xenv2 = exp(-1.*xenv2**2)
           xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
           xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
           exp_xth1 = exp(-1.*xth1**2)
           exp_xth2 = exp(-1.*xth2**2)
           cth(ind1)=0.5*(1.-1.*erf(xth1))
           cenv(ind1)=0.5*(1.-1.*erf(xenv1))
           ctot(ind1)=frac_th(ind1)*cth(ind1)+(1.-1.*frac_th(ind1))*cenv(ind1) 
           ctotth(ind1)=frac_th(ind1)*cth(ind1)
        

!volume cloud fraction and condensed water 

            !environnement

            IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
            IntJ_CF=0.5*(1.-1.*erf(xenv2))

            IF (deltasenv .LT. 1.e-10) THEN 
              qcenv(ind1)=IntJ
              cenv_vol(ind1)=IntJ_CF
            ELSE
              IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
              IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
              IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
              IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
              IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
              qcenv(ind1)=IntJ+IntI1+IntI2+IntI3
              cenv_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
            ENDIF
             


            !thermals

            IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
            IntJ_CF=0.5*(1.-1.*erf(xth2))
     
            IF (deltasth .LT. 1.e-10) THEN
              qcth(ind1)=IntJ
              cth_vol(ind1)=IntJ_CF
            ELSE
              IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
              IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
              IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
              IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
              qcth(ind1)=IntJ+IntI1+IntI2+IntI3
              cth_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
            ENDIF

            ! total 

            qctot(ind1)=frac_th(ind1)*qcth(ind1)+(1.-1.*frac_th(ind1))*qcenv(ind1)
            ctot_vol(ind1)=frac_th(ind1)*cth_vol(ind1)+(1.-1.*frac_th(ind1))*cenv_vol(ind1)

            IF (cenv(ind1).LT.min_neb_th.and.cth(ind1).LT.min_neb_th) THEN
                ctot(ind1)=0.
                ctot_vol(ind1)=0.
                qcloud(ind1)=qsenv(ind1)
                qincloud(ind1)=0.
            ELSE             
                qincloud(ind1)=qctot(ind1)/ctot(ind1)
                !to prevent situations with cloud condensed water greater than available total water
                qincloud(ind1)=min(qincloud(ind1),qt(ind1)/ctot(ind1))
                ! we assume that water vapor in cloud is qsenv
                qcloud(ind1)=qincloud(ind1)+qsenv(ind1)
            ENDIF 



           ! Outputs used to check the PDFs
           cloudth_senv(ind1) = senv
           cloudth_sth(ind1) = sth
           cloudth_sigmaenv(ind1) = sigma1s
           cloudth_sigmath(ind1) = sigma2s

      ENDIF       ! selection of grid points concerned by thermals 


    ENDDO       !loop on klon


RETURN


END SUBROUTINE condensation_cloudth


!*****************************************************************************************
!*****************************************************************************************
! pre-cmip7 routines are below and are becoming obsolete
!*****************************************************************************************
!*****************************************************************************************


       SUBROUTINE cloudth(ngrid,klev,ind2,  &
     &           ztv,po,zqta,fraca, & 
     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
     &           ratqs,zqs,t, &
     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)


      use lmdz_lscp_ini, only: iflag_cloudth_vert,iflag_ratqs

      USE yomcst_mod_h
      USE yoethf_mod_h
IMPLICIT NONE


!===========================================================================
! Auteur : Arnaud Octavio Jam (LMD/CNRS)
! Date : 25 Mai 2010
! Objet : calcule les valeurs de qc et rneb dans les thermiques
!===========================================================================

      INCLUDE "FCTTRE.h"

      INTEGER itap,ind1,ind2
      INTEGER ngrid,klev,klon,l,ig 
      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
      
      REAL ztv(ngrid,klev)
      REAL po(ngrid)
      REAL zqenv(ngrid)   
      REAL zqta(ngrid,klev)
          
      REAL fraca(ngrid,klev+1)
      REAL zpspsk(ngrid,klev)
      REAL paprs(ngrid,klev+1)
      REAL pplay(ngrid,klev)
      REAL ztla(ngrid,klev)
      REAL zthl(ngrid,klev)

      REAL zqsatth(ngrid,klev)
      REAL zqsatenv(ngrid,klev)
      
      
      REAL sigma1(ngrid,klev)
      REAL sigma2(ngrid,klev)
      REAL qlth(ngrid,klev)
      REAL qlenv(ngrid,klev)
      REAL qltot(ngrid,klev) 
      REAL cth(ngrid,klev)  
      REAL cenv(ngrid,klev)   
      REAL ctot(ngrid,klev)
      REAL rneb(ngrid,klev)
      REAL t(ngrid,klev)
      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
      REAL rdd,cppd,Lv
      REAL alth,alenv,ath,aenv
      REAL sth,senv,sigma1s,sigma2s,xth,xenv
      REAL Tbef,zdelta,qsatbef,zcor
      REAL qlbef  
      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
      
      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
      REAL zqs(ngrid), qcloud(ngrid)




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Gestion de deux versions de cloudth
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      IF (iflag_cloudth_vert.GE.1) THEN
      CALL cloudth_vert(ngrid,klev,ind2,  &
     &           ztv,po,zqta,fraca, & 
     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
     &           ratqs,zqs,t)
      RETURN
      ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!-------------------------------------------------------------------------------
! Initialisation des variables r?elles
!-------------------------------------------------------------------------------
      sigma1(:,ind2)=0.
      sigma2(:,ind2)=0.
      qlth(:,ind2)=0.
      qlenv(:,ind2)=0.  
      qltot(:,ind2)=0.
      rneb(:,ind2)=0.
      qcloud(:)=0.
      cth(:,ind2)=0.
      cenv(:,ind2)=0.
      ctot(:,ind2)=0.
      qsatmmussig1=0.
      qsatmmussig2=0.
      rdd=287.04
      cppd=1005.7
      pi=3.14159 
      Lv=2.5e6
      sqrt2pi=sqrt(2.*pi)



!-------------------------------------------------------------------------------
! Calcul de la fraction du thermique et des ?cart-types des distributions
!-------------------------------------------------------------------------------                 
      do ind1=1,ngrid

      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then 

      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))


!      zqenv(ind1)=po(ind1)
      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef




      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
      aenv=1./(1.+(alenv*Lv/cppd))
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 




      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatth(ind1,ind2)=qsatbef
            
      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
      ath=1./(1.+(alth*Lv/cppd))
      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) 
      
      

!------------------------------------------------------------------------------
! Calcul des ?cart-types pour s
!------------------------------------------------------------------------------

!      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2) 
!       if (paprs(ind1,ind2).gt.90000) then
!       ratqs(ind1,ind2)=0.002
!       else 
!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
!       endif
       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 
!       sigma1s=ratqs(ind1,ind2)*po(ind1)
!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003  
 
!------------------------------------------------------------------------------
! Calcul de l'eau condens?e et de la couverture nuageuse
!------------------------------------------------------------------------------
      sqrt2pi=sqrt(2.*pi)
      xth=sth/(sqrt(2.)*sigma2s)
      xenv=senv/(sqrt(2.)*sigma1s)
      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) 
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   

      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      if (ctot(ind1,ind2).lt.1.e-10) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2) 

      else   
                
      ctot(ind1,ind2)=ctot(ind1,ind2) 
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)

      endif                           
      
          


      else  ! gaussienne environnement seule
      
      zqenv(ind1)=po(ind1)
      Tbef=t(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef
      

!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
      aenv=1./(1.+(alenv*Lv/cppd))
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 
      

      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)

      sqrt2pi=sqrt(2.*pi)
      xenv=senv/(sqrt(2.)*sigma1s)
      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
      
      if (ctot(ind1,ind2).lt.1.e-3) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2) 

      else   
                
      ctot(ind1,ind2)=ctot(ind1,ind2) 
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)

      endif    
 
 
 
 
 
 
      endif   
      enddo
     
      return
!     end
END SUBROUTINE cloudth



!===========================================================================
     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
     &           ztv,po,zqta,fraca, & 
     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
     &           ratqs,zqs,t)

!===========================================================================
! Auteur : Arnaud Octavio Jam (LMD/CNRS)
! Date : 25 Mai 2010
! Objet : calcule les valeurs de qc et rneb dans les thermiques
!===========================================================================


USE yoethf_mod_h
            use lmdz_lscp_ini, only: iflag_cloudth_vert, vert_alpha

      USE yomcst_mod_h
IMPLICIT NONE


      INCLUDE "FCTTRE.h"
      
      INTEGER itap,ind1,ind2
      INTEGER ngrid,klev,klon,l,ig 
      
      REAL ztv(ngrid,klev)
      REAL po(ngrid)
      REAL zqenv(ngrid)   
      REAL zqta(ngrid,klev)
          
      REAL fraca(ngrid,klev+1)
      REAL zpspsk(ngrid,klev)
      REAL paprs(ngrid,klev+1)
      REAL pplay(ngrid,klev)
      REAL ztla(ngrid,klev)
      REAL zthl(ngrid,klev)

      REAL zqsatth(ngrid,klev)
      REAL zqsatenv(ngrid,klev)
      
      
      REAL sigma1(ngrid,klev)                                                         
      REAL sigma2(ngrid,klev)
      REAL qlth(ngrid,klev)
      REAL qlenv(ngrid,klev)
      REAL qltot(ngrid,klev) 
      REAL cth(ngrid,klev)  
      REAL cenv(ngrid,klev)   
      REAL ctot(ngrid,klev)
      REAL rneb(ngrid,klev)
      REAL t(ngrid,klev)                                                                  
      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
      REAL rdd,cppd,Lv,sqrt2,sqrtpi
      REAL alth,alenv,ath,aenv
      REAL sth,senv,sigma1s,sigma2s,xth,xenv
      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
      REAL Tbef,zdelta,qsatbef,zcor
      REAL qlbef  
      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
      ! (J Jouhaud, JL Dufresne, JB Madeleine)
      
      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
      REAL zqs(ngrid), qcloud(ngrid)

!------------------------------------------------------------------------------
! Initialisation des variables r?elles
!------------------------------------------------------------------------------
      sigma1(:,ind2)=0.
      sigma2(:,ind2)=0.
      qlth(:,ind2)=0.
      qlenv(:,ind2)=0.  
      qltot(:,ind2)=0.
      rneb(:,ind2)=0.
      qcloud(:)=0.
      cth(:,ind2)=0.
      cenv(:,ind2)=0.
      ctot(:,ind2)=0.
      qsatmmussig1=0.
      qsatmmussig2=0.
      rdd=287.04
      cppd=1005.7
      pi=3.14159 
      Lv=2.5e6
      sqrt2pi=sqrt(2.*pi)
      sqrt2=sqrt(2.)
      sqrtpi=sqrt(pi)

!-------------------------------------------------------------------------------
! Calcul de la fraction du thermique et des ?cart-types des distributions
!-------------------------------------------------------------------------------                 
      do ind1=1,ngrid

      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then 

      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))


!      zqenv(ind1)=po(ind1)
      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef




      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
      aenv=1./(1.+(alenv*Lv/cppd))
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 




      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatth(ind1,ind2)=qsatbef
            
      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
      ath=1./(1.+(alth*Lv/cppd))
      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) 
      
      

!------------------------------------------------------------------------------
! Calcul des ?cart-types pour s
!------------------------------------------------------------------------------

      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) 
!       if (paprs(ind1,ind2).gt.90000) then
!       ratqs(ind1,ind2)=0.002
!       else 
!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
!       endif
!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 
!       sigma1s=ratqs(ind1,ind2)*po(ind1)
!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003  
 
!------------------------------------------------------------------------------
! Calcul de l'eau condens?e et de la couverture nuageuse
!------------------------------------------------------------------------------
      sqrt2pi=sqrt(2.*pi)
      xth=sth/(sqrt(2.)*sigma2s)
      xenv=senv/(sqrt(2.)*sigma1s)
      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) 
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   

      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
     
       IF (iflag_cloudth_vert == 1) THEN
!-------------------------------------------------------------------------------
!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
!-------------------------------------------------------------------------------
!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
!      deltasenv=aenv*0.01*po(ind1)
!     deltasth=ath*0.01*zqta(ind1,ind2)    
      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
      
      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) 
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   

      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))

      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
!      qlenv(ind1,ind2)=IntJ 
!      print*, qlenv(ind1,ind2),'VERIF EAU'


      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
!      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
!      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
!      qlth(ind1,ind2)=IntJ
!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)

      ELSE IF (iflag_cloudth_vert == 2) THEN

!-------------------------------------------------------------------------------
!  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
!-------------------------------------------------------------------------------
!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
      deltasenv=aenv*vert_alpha*sigma1s
      deltasth=ath*vert_alpha*sigma2s
      
      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
!     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
!     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
      
      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)

      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))

!      IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
!      IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
!      IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))

      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
!      qlenv(ind1,ind2)=IntJ 
!      print*, qlenv(ind1,ind2),'VERIF EAU'

      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
      
      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
!      qlth(ind1,ind2)=IntJ
!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
      



      ENDIF ! of if (iflag_cloudth_vert==1 or 2)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2) 

      else 
                
      ctot(ind1,ind2)=ctot(ind1,ind2) 
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) 

      endif  
                       
      
          
!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'


      else  ! gaussienne environnement seule
      
      zqenv(ind1)=po(ind1)
      Tbef=t(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef
      

!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
      aenv=1./(1.+(alenv*Lv/cppd))
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 
      

      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)

      sqrt2pi=sqrt(2.*pi)
      xenv=senv/(sqrt(2.)*sigma1s)
      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
      
      if (ctot(ind1,ind2).lt.1.e-3) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2) 

      else   
                
      ctot(ind1,ind2)=ctot(ind1,ind2) 
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)

      endif    
 
 
 
 
 
 
      endif   
      enddo
     
      return
!     end
END SUBROUTINE cloudth_vert




       SUBROUTINE cloudth_v3(ngrid,klev,ind2,  &
     &           ztv,po,zqta,fraca, & 
     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
     &           ratqs,sigma_qtherm,zqs,t, &
     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)

      use lmdz_lscp_ini, only: iflag_cloudth_vert

      USE yomcst_mod_h
      USE yoethf_mod_h
IMPLICIT NONE


!===========================================================================
! Author : Arnaud Octavio Jam (LMD/CNRS)
! Date : 25 Mai 2010
! Objet : calcule les valeurs de qc et rneb dans les thermiques
!===========================================================================
      INCLUDE "FCTTRE.h"

      integer, intent(in) :: ind2
      integer, intent(in) :: ngrid,klev
      
      real, dimension(ngrid,klev), intent(in) :: ztv
      real, dimension(ngrid), intent(in) :: po
      real, dimension(ngrid,klev), intent(in) :: zqta
      real, dimension(ngrid,klev+1), intent(in) :: fraca
      real, dimension(ngrid), intent(out) :: qcloud
      real, dimension(ngrid,klev), intent(out) :: ctot
      real, dimension(ngrid,klev), intent(out) :: ctot_vol
      real, dimension(ngrid,klev), intent(in) :: zpspsk
      real, dimension(ngrid,klev+1), intent(in) :: paprs
      real, dimension(ngrid,klev), intent(in) :: pplay
      real, dimension(ngrid,klev), intent(in) :: ztla
      real, dimension(ngrid,klev), intent(inout) :: zthl
      real, dimension(ngrid,klev), intent(in) :: ratqs,sigma_qtherm
      real, dimension(ngrid), intent(in) :: zqs
      real, dimension(ngrid,klev), intent(in) :: t
      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv


      REAL zqenv(ngrid)   
      REAL zqsatth(ngrid,klev)
      REAL zqsatenv(ngrid,klev)
      
      REAL sigma1(ngrid,klev)                                                         
      REAL sigma2(ngrid,klev)
      REAL qlth(ngrid,klev)
      REAL qlenv(ngrid,klev)
      REAL qltot(ngrid,klev) 
      REAL cth(ngrid,klev)
      REAL cenv(ngrid,klev)   
      REAL cth_vol(ngrid,klev)
      REAL cenv_vol(ngrid,klev)
      REAL rneb(ngrid,klev)      
      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
      REAL rdd,cppd,Lv
      REAL alth,alenv,ath,aenv
      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
      REAL Tbef,zdelta,qsatbef,zcor
      REAL qlbef  
      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 


      INTEGER :: ind1,l, ig

      IF (iflag_cloudth_vert.GE.1) THEN
      CALL cloudth_vert_v3(ngrid,klev,ind2,  &
     &           ztv,po,zqta,fraca, & 
     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
     &           ratqs,sigma_qtherm,zqs,t, &
     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
      RETURN
      ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!-------------------------------------------------------------------------------
! Initialisation des variables r?elles
!-------------------------------------------------------------------------------
      sigma1(:,ind2)=0.
      sigma2(:,ind2)=0.
      qlth(:,ind2)=0.
      qlenv(:,ind2)=0.  
      qltot(:,ind2)=0.
      rneb(:,ind2)=0.
      qcloud(:)=0.
      cth(:,ind2)=0.
      cenv(:,ind2)=0.
      ctot(:,ind2)=0.
      cth_vol(:,ind2)=0.
      cenv_vol(:,ind2)=0.
      ctot_vol(:,ind2)=0.
      qsatmmussig1=0.
      qsatmmussig2=0.
      rdd=287.04
      cppd=1005.7
      pi=3.14159 
      Lv=2.5e6
      sqrt2pi=sqrt(2.*pi)
      sqrt2=sqrt(2.)
      sqrtpi=sqrt(pi)


!-------------------------------------------------------------------------------
! Cloud fraction in the thermals and standard deviation of the PDFs
!-------------------------------------------------------------------------------                 
      do ind1=1,ngrid

      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then 

      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))

      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef


      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84

!po = qt de l'environnement ET des thermique
!zqenv = qt environnement
!zqsatenv = qsat environnement
!zthl = Tl environnement


      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatth(ind1,ind2)=qsatbef
            
      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
      
!zqta = qt thermals
!zqsatth = qsat thermals
!ztla = Tl thermals

!------------------------------------------------------------------------------
! s standard deviations
!------------------------------------------------------------------------------

!     tests
!     sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
!     sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
!     sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
!     final option
      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
 
!------------------------------------------------------------------------------
! Condensed water and cloud cover
!------------------------------------------------------------------------------
      xth=sth/(sqrt2*sigma2s)
      xenv=senv/(sqrt2*sigma1s)
      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))       !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
      ctot_vol(ind1,ind2)=ctot(ind1,ind2)

      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)

      if (ctot(ind1,ind2).lt.1.e-10) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2) 
      else
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
      endif

      else  ! Environnement only, follow the if l.110
      
      zqenv(ind1)=po(ind1)
      Tbef=t(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef

!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
      aenv=1./(1.+(alenv*Lv/cppd))
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
      
      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)

      xenv=senv/(sqrt2*sigma1s)
      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))

      if (ctot(ind1,ind2).lt.1.e-3) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2)
      else   
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
      endif


      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
      enddo       ! from the loop on ngrid l.108
      return
!     end
END SUBROUTINE cloudth_v3



!===========================================================================
     SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2,  &
     &           ztv,po,zqta,fraca, & 
     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
     &           ratqs,sigma_qtherm,zqs,t, &
     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)

!===========================================================================
! Auteur : Arnaud Octavio Jam (LMD/CNRS)
! Date : 25 Mai 2010
! Objet : calcule les valeurs de qc et rneb dans les thermiques
!===========================================================================

      use yoethf_mod_h
      use lmdz_lscp_ini, only : iflag_cloudth_vert,iflag_ratqs
      use lmdz_lscp_ini, only : vert_alpha,vert_alpha_th, sigma1s_factor, sigma1s_power , sigma2s_factor , sigma2s_power , cloudth_ratqsmin , iflag_cloudth_vert_noratqs

      USE yomcst_mod_h
IMPLICIT NONE




      INCLUDE "FCTTRE.h"
      
      INTEGER itap,ind1,ind2
      INTEGER ngrid,klev,klon,l,ig 
      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
      
      REAL ztv(ngrid,klev)
      REAL po(ngrid)
      REAL zqenv(ngrid)   
      REAL zqta(ngrid,klev)
          
      REAL fraca(ngrid,klev+1)
      REAL zpspsk(ngrid,klev)
      REAL paprs(ngrid,klev+1)
      REAL pplay(ngrid,klev)
      REAL ztla(ngrid,klev)
      REAL zthl(ngrid,klev)

      REAL zqsatth(ngrid,klev)
      REAL zqsatenv(ngrid,klev)
      
      REAL sigma1(ngrid,klev)                                                         
      REAL sigma2(ngrid,klev)
      REAL qlth(ngrid,klev)
      REAL qlenv(ngrid,klev)
      REAL qltot(ngrid,klev) 
      REAL cth(ngrid,klev)
      REAL cenv(ngrid,klev)   
      REAL ctot(ngrid,klev)
      REAL cth_vol(ngrid,klev)
      REAL cenv_vol(ngrid,klev)
      REAL ctot_vol(ngrid,klev)
      REAL rneb(ngrid,klev)
      REAL t(ngrid,klev)                                                                  
      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
      REAL rdd,cppd,Lv
      REAL alth,alenv,ath,aenv
      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
      REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
      REAL Tbef,zdelta,qsatbef,zcor
      REAL qlbef  
      REAL ratqs(ngrid,klev),sigma_qtherm(ngrid,klev) ! determine la largeur de distribution de vapeur
      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
      ! (J Jouhaud, JL Dufresne, JB Madeleine)

      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
      REAL zqs(ngrid), qcloud(ngrid)

      REAL rhodz(ngrid,klev)
      REAL zrho(ngrid,klev)
      REAL dz(ngrid,klev)

      DO ind1 = 1, ngrid
        !Layer calculation
        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2
        zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3
        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre
      END DO

!------------------------------------------------------------------------------
! Initialize
!------------------------------------------------------------------------------

      sigma1(:,ind2)=0.
      sigma2(:,ind2)=0.
      qlth(:,ind2)=0.
      qlenv(:,ind2)=0.  
      qltot(:,ind2)=0.
      rneb(:,ind2)=0.
      qcloud(:)=0.
      cth(:,ind2)=0.
      cenv(:,ind2)=0.
      ctot(:,ind2)=0.
      cth_vol(:,ind2)=0.
      cenv_vol(:,ind2)=0.
      ctot_vol(:,ind2)=0.
      qsatmmussig1=0.
      qsatmmussig2=0.
      rdd=287.04
      cppd=1005.7
      pi=3.14159 
      Lv=2.5e6
      sqrt2pi=sqrt(2.*pi)
      sqrt2=sqrt(2.)
      sqrtpi=sqrt(pi)



!-------------------------------------------------------------------------------
! Calcul de la fraction du thermique et des ecart-types des distributions
!-------------------------------------------------------------------------------                 
      do ind1=1,ngrid

      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement

      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv


      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef


      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84

!zqenv = qt environnement
!zqsatenv = qsat environnement
!zthl = Tl environnement


      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatth(ind1,ind2)=qsatbef
            
      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
      
      
!zqta = qt thermals
!zqsatth = qsat thermals
!ztla = Tl thermals
!------------------------------------------------------------------------------
! s standard deviation
!------------------------------------------------------------------------------

      sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
     &                (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
!     sigma1s_fraca = (1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
      IF (cloudth_ratqsmin>0.) THEN
         sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
      ELSE
         sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
      ENDIF
      sigma1s = sigma1s_fraca + sigma1s_ratqs
      sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
      IF (iflag_ratqs.eq.10.or.iflag_ratqs.eq.11) then
         sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
         IF (iflag_ratqs.eq.10.and.sigma_qtherm(ind1,ind2).ne.0) then
            sigma2s = sigma_qtherm(ind1,ind2)*ath
         ENDIF
      ENDIF
      
!      tests
!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
!      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
!      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) 
!      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
!       if (paprs(ind1,ind2).gt.90000) then
!       ratqs(ind1,ind2)=0.002
!       else 
!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
!       endif
!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 
!       sigma1s=ratqs(ind1,ind2)*po(ind1)
!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003  
 
       IF (iflag_cloudth_vert == 1) THEN
!-------------------------------------------------------------------------------
!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
!-------------------------------------------------------------------------------

      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)

      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
      
      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) 
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   

      ! Environment
      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))

      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3

      ! Thermal
      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)

      ELSE IF (iflag_cloudth_vert >= 3) THEN
      IF (iflag_cloudth_vert < 5) THEN
!-------------------------------------------------------------------------------
!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
!-------------------------------------------------------------------------------
!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
      IF (iflag_cloudth_vert == 3) THEN
        deltasenv=aenv*vert_alpha*sigma1s
        deltasth=ath*vert_alpha_th*sigma2s
      ELSE IF (iflag_cloudth_vert == 4) THEN
        IF (iflag_cloudth_vert_noratqs == 1) THEN
          deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
          deltasth=vert_alpha_th*sigma2s
        ELSE
          deltasenv=vert_alpha*sigma1s
          deltasth=vert_alpha_th*sigma2s
        ENDIF
      ENDIF
      
      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
      exp_xenv1 = exp(-1.*xenv1**2)
      exp_xenv2 = exp(-1.*xenv2**2)
      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
      exp_xth1 = exp(-1.*xth1**2)
      exp_xth2 = exp(-1.*xth2**2)
      
      !CF_surfacique
      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) 


      !CF_volumique & eau condense
      !environnement
      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
      IntJ_CF=0.5*(1.-1.*erf(xenv2))
      if (deltasenv .lt. 1.e-10) then 
      qlenv(ind1,ind2)=IntJ
      cenv_vol(ind1,ind2)=IntJ_CF
      else
      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
      IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
      IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
      cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
      endif

      !thermique
      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
      IntJ_CF=0.5*(1.-1.*erf(xth2))
      if (deltasth .lt. 1.e-10) then
      qlth(ind1,ind2)=IntJ
      cth_vol(ind1,ind2)=IntJ_CF
      else
      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
      IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
      IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
      cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
      endif

      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)

      ELSE IF (iflag_cloudth_vert == 5) THEN
         sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) &
              /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) &
              +ratqs(ind1,ind2)*po(ind1) !Environment
      sigma2s=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.02)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2)                   !Thermals
      !sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
      !sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 
      xth=sth/(sqrt(2.)*sigma2s)
      xenv=senv/(sqrt(2.)*sigma1s)

      !Volumique
      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
      !print *,'jeanjean_CV=',ctot_vol(ind1,ind2)

      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth_vol(ind1,ind2))
      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv_vol(ind1,ind2))  
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)

      !Surfacique
      !Neggers
      !beta=0.0044
      !inverse_rho=1.+beta*dz(ind1,ind2)
      !print *,'jeanjean : beta=',beta
      !cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
      !cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
      !ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)

      !Brooks
      a_Brooks=0.6694
      b_Brooks=0.1882
      A_Maj_Brooks=0.1635 !-- sans shear
      !A_Maj_Brooks=0.17   !-- ARM LES
      !A_Maj_Brooks=0.18   !-- RICO LES
      !A_Maj_Brooks=0.19   !-- BOMEX LES
      Dx_Brooks=200000.
      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
      !print *,'jeanjean_f=',f_Brooks

      cth(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
      cenv(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
      ctot(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
      !print *,'JJ_ctot_1',ctot(ind1,ind2)





      ENDIF ! of if (iflag_cloudth_vert<5)
      ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)

!      if (ctot(ind1,ind2).lt.1.e-10) then
      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
      ctot(ind1,ind2)=0.
      ctot_vol(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2)

      else
                
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) 

      endif  

      else  ! gaussienne environnement seule
      

      zqenv(ind1)=po(ind1)
      Tbef=t(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef
      

!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
      aenv=1./(1.+(alenv*Lv/cppd))
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
      sth=0.
      

      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
      sigma2s=0.

      sqrt2pi=sqrt(2.*pi)
      xenv=senv/(sqrt(2.)*sigma1s)
      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
      
      if (ctot(ind1,ind2).lt.1.e-3) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2)

      else   
                
!      ctot(ind1,ind2)=ctot(ind1,ind2) 
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)

      endif  
 



      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
      ! Outputs used to check the PDFs
      cloudth_senv(ind1,ind2) = senv
      cloudth_sth(ind1,ind2) = sth
      cloudth_sigmaenv(ind1,ind2) = sigma1s
      cloudth_sigmath(ind1,ind2) = sigma2s

      enddo       ! from the loop on ngrid l.333
      return
!     end
END SUBROUTINE cloudth_vert_v3
!











       SUBROUTINE cloudth_v6(ngrid,klev,ind2,  &
     &           ztv,po,zqta,fraca, & 
     &           qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
     &           ratqs,zqs,T, &
     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)

      USE yoethf_mod_h
      USE lmdz_lscp_ini, only: iflag_cloudth_vert

      USE yomcst_mod_h
IMPLICIT NONE



      INCLUDE "FCTTRE.h"


        !Domain variables
      INTEGER ngrid !indice Max lat-lon
      INTEGER klev  !indice Max alt
      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
      INTEGER ind1  !indice in [1:ngrid]
      INTEGER ind2  !indice in [1:klev]
        !thermal plume fraction
      REAL fraca(ngrid,klev+1)   !thermal plumes fraction in the gridbox
        !temperatures
      REAL T(ngrid,klev)       !temperature
      REAL zpspsk(ngrid,klev)  !factor (p/p0)**kappa (used for potential variables)
      REAL ztv(ngrid,klev)     !potential temperature (voir thermcell_env.F90)
      REAL ztla(ngrid,klev)    !liquid temperature in the thermals (Tl_th)
      REAL zthl(ngrid,klev)    !liquid temperature in the environment (Tl_env)
        !pressure
      REAL paprs(ngrid,klev+1)   !pressure at the interface of levels
      REAL pplay(ngrid,klev)     !pressure at the middle of the level
        !humidity
      REAL ratqs(ngrid,klev)   !width of the total water subgrid-scale distribution
      REAL po(ngrid)           !total water (qt)
      REAL zqenv(ngrid)        !total water in the environment (qt_env)
      REAL zqta(ngrid,klev)    !total water in the thermals (qt_th)
      REAL zqsatth(ngrid,klev)   !water saturation level in the thermals (q_sat_th)
      REAL zqsatenv(ngrid,klev)  !water saturation level in the environment (q_sat_env)
      REAL qlth(ngrid,klev)    !condensed water in the thermals
      REAL qlenv(ngrid,klev)   !condensed water in the environment
      REAL qltot(ngrid,klev)   !condensed water in the gridbox
        !cloud fractions
      REAL cth_vol(ngrid,klev)   !cloud fraction by volume in the thermals
      REAL cenv_vol(ngrid,klev)  !cloud fraction by volume in the environment
      REAL ctot_vol(ngrid,klev)  !cloud fraction by volume in the gridbox
      REAL cth_surf(ngrid,klev)  !cloud fraction by surface in the thermals
      REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment  
      REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
        !PDF of saturation deficit variables
      REAL rdd,cppd,Lv
      REAL Tbef,zdelta,qsatbef,zcor
      REAL alth,alenv,ath,aenv
      REAL sth,senv              !saturation deficits in the thermals and environment
      REAL sigma_env,sigma_th    !standard deviations of the biGaussian PDF
        !cloud fraction variables
      REAL xth,xenv
      REAL inverse_rho,beta                                  !Neggers et al. (2011) method
      REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
        !Incloud total water variables
      REAL zqs(ngrid)    !q_sat
      REAL qcloud(ngrid) !eau totale dans le nuage
        !Some arithmetic variables
      REAL  pi,sqrt2,sqrt2pi
        !Depth of the layer
      REAL dz(ngrid,klev)    !epaisseur de la couche en metre
      REAL rhodz(ngrid,klev)
      REAL zrho(ngrid,klev)
      DO ind1 = 1, ngrid
        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
        zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd          ![kg/m3]
        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2)            ![m]
      END DO

!------------------------------------------------------------------------------
! Initialization
!------------------------------------------------------------------------------
      qlth(:,ind2)=0.
      qlenv(:,ind2)=0.  
      qltot(:,ind2)=0.
      cth_vol(:,ind2)=0.
      cenv_vol(:,ind2)=0.
      ctot_vol(:,ind2)=0.
      cth_surf(:,ind2)=0.
      cenv_surf(:,ind2)=0.
      ctot_surf(:,ind2)=0.
      qcloud(:)=0.
      rdd=287.04
      cppd=1005.7
      pi=3.14159 
      Lv=2.5e6
      sqrt2=sqrt(2.)
      sqrt2pi=sqrt(2.*pi)


      DO ind1=1,ngrid
!-------------------------------------------------------------------------------
!Both thermal and environment in the gridbox
!-------------------------------------------------------------------------------
      IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN
        !--------------------------------------------
        !calcul de qsat_env
        !--------------------------------------------
      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef
        !--------------------------------------------
        !calcul de s_env
        !--------------------------------------------
      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam
        !--------------------------------------------
        !calcul de qsat_th
        !--------------------------------------------
      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatth(ind1,ind2)=qsatbef
        !--------------------------------------------
        !calcul de s_th  
        !--------------------------------------------
      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84 these Arnaud Jam
      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84 these Arnaud Jam
      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84 these Arnaud Jam
        !--------------------------------------------
        !calcul standard deviations bi-Gaussian PDF
        !--------------------------------------------
      sigma_th=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.01)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2)
      sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) &
           /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) &
           +ratqs(ind1,ind2)*po(ind1)
      xth=sth/(sqrt2*sigma_th)
      xenv=senv/(sqrt2*sigma_env)
        !--------------------------------------------
        !Cloud fraction by volume CF_vol
        !--------------------------------------------
      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
        !-------------------------------------------- 
        !Condensed water qc
        !--------------------------------------------
      qlth(ind1,ind2)=sigma_th*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth_vol(ind1,ind2))
      qlenv(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv_vol(ind1,ind2))  
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
        !--------------------------------------------
        !Cloud fraction by surface CF_surf
        !--------------------------------------------
        !Method Neggers et al. (2011) : ok for cumulus clouds only
      !beta=0.0044 (Jouhaud et al.2018)
      !inverse_rho=1.+beta*dz(ind1,ind2)
      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
        !Method Brooks et al. (2005) : ok for all types of clouds
      a_Brooks=0.6694
      b_Brooks=0.1882
      A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
      Dx_Brooks=200000.   !-- si l'on considere des mailles de 200km de cote
      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
        !--------------------------------------------
        !Incloud Condensed water qcloud
        !--------------------------------------------
      if (ctot_surf(ind1,ind2) .lt. 1.e-10) then
      ctot_vol(ind1,ind2)=0.
      ctot_surf(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2)
      else
      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
      endif



!-------------------------------------------------------------------------------
!Environment only in the gridbox
!-------------------------------------------------------------------------------
      ELSE
        !--------------------------------------------
        !calcul de qsat_env
        !--------------------------------------------
      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef
        !--------------------------------------------
        !calcul de s_env
        !--------------------------------------------
      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam 
        !--------------------------------------------
        !calcul standard deviations Gaussian PDF
        !--------------------------------------------
      zqenv(ind1)=po(ind1)
      sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
      xenv=senv/(sqrt2*sigma_env)
        !--------------------------------------------
        !Cloud fraction by volume CF_vol
        !--------------------------------------------
      ctot_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
        !--------------------------------------------
        !Condensed water qc
        !--------------------------------------------
      qltot(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*ctot_vol(ind1,ind2))
        !--------------------------------------------
        !Cloud fraction by surface CF_surf
        !--------------------------------------------
        !Method Neggers et al. (2011) : ok for cumulus clouds only
      !beta=0.0044 (Jouhaud et al.2018)
      !inverse_rho=1.+beta*dz(ind1,ind2)
      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
        !Method Brooks et al. (2005) : ok for all types of clouds
      a_Brooks=0.6694
      b_Brooks=0.1882
      A_Maj_Brooks=0.1635 !-- sans dependence au shear
      Dx_Brooks=200000.
      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
        !--------------------------------------------
        !Incloud Condensed water qcloud
        !--------------------------------------------
      if (ctot_surf(ind1,ind2) .lt. 1.e-8) then
      ctot_vol(ind1,ind2)=0.
      ctot_surf(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2)
      else
      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
      endif


      END IF  ! From the separation (thermal/envrionnement) et (environnement only)

      ! Outputs used to check the PDFs
      cloudth_senv(ind1,ind2) = senv
      cloudth_sth(ind1,ind2) = sth
      cloudth_sigmaenv(ind1,ind2) = sigma_env
      cloudth_sigmath(ind1,ind2) = sigma_th

      END DO  ! From the loop on ngrid
      return

END SUBROUTINE cloudth_v6



END MODULE lmdz_lscp_condensation
