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
! reference: Bony and Emanuel (2001), doi:10.1175/1520-0469(2001)058
!----------------------------------------------------------------------

      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)
            END IF

         END IF   ! end keepgoing

      END DO     ! 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
! contact: E. Vignon, etienne.vignon@lmd.ipsl.fr
! reference: Borella et al. (2025), doi:10.1029/2024MS004918
!----------------------------------------------------------------------

      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)
               END IF

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

               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.
                        END IF
                     ELSE
                        !--We keep the saturation adjustment hypothesis, and the vapor in the
                        !--cloud is set equal to the saturation vapor
                        qvapincld_new = qsat(i)
                     END IF ! 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)

                        END IF

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

                     END IF ! qvapincld .LT. qvapincld_new

                  END IF   ! qiceincld .LT. eps
               END IF ! 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.
                  END IF
                  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)
                  END IF

                  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)
                     END IF

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

                  END IF ! 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)

               END IF ! ( 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
                        END IF
                     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)
                        END IF
                     END IF ! ok_unadjusted_clouds
                  END IF ! 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)
                     END IF ! ok_unadjusted_clouds

                  END IF ! 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)))

               END IF ! ( ( 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)
               END IF ! 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

            END IF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds

         END IF   ! end keepgoing

      END DO     ! 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     ! in-cloud specific humidity [kg/kg]
      REAL, INTENT(IN) :: qiceincld     ! in-cloud specific ice content [kg/kg]
      REAL, INTENT(IN) :: temp          ! temperature [K]
      REAL, INTENT(IN) :: qsat          ! saturation specific humidity [kg/kg]
      REAL, INTENT(IN) :: pplay         ! pressure [Pa]
      REAL, INTENT(IN) :: dtime         ! time step [s]
      REAL, INTENT(OUT):: qvapincld_new ! update qvapincld [kg/kg]

!
! 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)
         END IF ! qice_ratio .LT. 0.
      END IF ! 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
! Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam
! IMPORTANT NOTE: we assume iflag_cloudth_vert=7
! contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr
! Date: February 2025
! reference: Jam et al. (2013), doi: 10.1007/s10546-012-9789-3
!-----------------------------------------------------------------------------------

      use lmdz_lscp_ini, only: iflag_cloudth_vert, iflag_ratqs, iflag_cloudth_vert_noratqs
      use lmdz_lscp_ini, only: vert_alpha, vert_alpha_th
      use lmdz_lscp_ini, only: 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

      END DO

      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*Tbefenv(ind1)**2)
            aenv = 1./(1.+(alenv*RLVTT/RCPD))
            senv = aenv*(qt_env(ind1) - qsenv(ind1))

! Thermals:

            alth = (RD/RV*RLVTT*qsth(ind1))/(RD*Tbefth(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)
            END IF
            sigma1s = sigma1s_fraca + sigma1s_ratqs

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

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

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

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

            ! 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)
            END IF

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

         END IF       ! selection of grid points concerned by thermals

      END DO       !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, ok_bug_deficit

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

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

            IF (ok_bug_deficit) THEN
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*zthl(ind1, ind2)**2)
            ELSE
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*Tbef**2)
            ENDIF   
            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

            IF (ok_bug_deficit) THEN
               alth = (0.622*Lv*zqsatth(ind1, ind2))/(rdd*ztla(ind1, ind2)**2)
            ELSE
               alth = (0.622*Lv*zqsatth(ind1, ind2))/(rdd*Tbef**2)
            END IF
            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)

            end if

         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)
            IF (ok_bug_deficit) THEN
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*zthl(ind1, ind2)**2)
            ELSE
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*Tbef**2)
            ENDIF
            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)

            end if

         end if
      end do

      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, ok_bug_deficit
      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
            IF (ok_bug_deficit) THEN
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*zthl(ind1, ind2)**2)
            ELSE
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*Tbef**2)
            ENDIF
            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
            IF (ok_bug_deficit) THEN
               alth = (0.622*Lv*zqsatth(ind1, ind2))/(rdd*ztla(ind1, ind2)**2)
            ELSE
               alth = (0.622*Lv*zqsatth(ind1, ind2))/(rdd*Tbef**2)
            ENDIF
            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)

            END IF ! 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)

            end if

!     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)
            
            IF (ok_bug_deficit) THEN
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*zthl(ind1, ind2)**2)
            ELSE
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*Tbef**2)
            ENDIF
            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)

            end if

         end if
      end do

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

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

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

            IF (ok_bug_deficit) THEN
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*zthl(ind1, ind2)**2)     !qsl, p84
            ELSE
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*Tbef**2)     !qsl, p84
            ENDIF
            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

            IF (ok_bug_deficit) THEN
               alth = (0.622*Lv*zqsatth(ind1, ind2))/(rdd*ztla(ind1, ind2)**2)       !qsl, p84
            ELSE
               alth = (0.622*Lv*zqsatth(ind1, ind2))/(rdd*Tbef**2)       !qsl, p84
            ENDIF    
            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)
            end if

         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)
            IF (ok_bug_deficit) THEN
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*zthl(ind1, ind2)**2)
            ELSE
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*Tbef**2)
            ENDIF
            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)
            end if

         end if       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
      end do       ! 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, ok_bug_deficit
      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
            IF (ok_bug_deficit) THEN
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*zthl(ind1, ind2)**2)     !qsl, p84
            ELSE
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*Tbef**2)     !qsl, p84
            ENDIF
            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
            IF (ok_bug_deficit) THEN
               alth = (0.622*Lv*zqsatth(ind1, ind2))/(rdd*ztla(ind1, ind2)**2)       !qsl, p84
            ELSE
               alth = (0.622*Lv*zqsatth(ind1, ind2))/(rdd*Tbef**2)       !qsl, p84
            ENDIF
            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)
            END IF
            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
               END IF
            END IF

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

                  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
                  end if

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

                  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)

               END IF ! of if (iflag_cloudth_vert<5)
            END IF ! 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)

            end if

         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)
            
            IF (ok_bug_deficit) THEN
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*zthl(ind1, ind2)**2)
            ELSE
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*Tbef**2)
            ENDIF
            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)

            end if

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

      end do       ! 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, ok_bug_deficit

      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
            !--------------------------------------------
            IF (ok_bug_deficit) THEN
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*zthl(ind1, ind2)**2)     !qsl, p84 these Arnaud Jam
            ELSE
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*Tbef**2)     !qsl, p84 these Arnaud Jam
            ENDIF
            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
            !--------------------------------------------
            IF (ok_bug_deficit) THEN
               alth = (0.622*Lv*zqsatth(ind1, ind2))/(rdd*ztla(ind1, ind2)**2)       !qsl, p84 these Arnaud Jam
            ELSE
               alth = (0.622*Lv*zqsatth(ind1, ind2))/(rdd*Tbef**2)       !qsl, p84 these Arnaud Jam
            ENDIF
            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)
            end if

!-------------------------------------------------------------------------------
!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
            !--------------------------------------------
            IF (ok_bug_deficit) THEN
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*zthl(ind1, ind2)**2)     !qsl, p84 these Arnaud Jam
            ELSE
               alenv = (0.622*Lv*zqsatenv(ind1, ind2))/(rdd*Tbef**2)     !qsl, p84 these Arnaud Jam
            ENDIF
            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)
            end if

         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
