Changeset 5614 for LMDZ6/trunk

Timestamp:

Apr 14, 2025, 9:21:07 PM (8 weeks ago)

Author:

evignon

Message:

Commission liée à un update majeur de la routine de condensation grande echelle suite au travail
de Lea, Audran et Etienne
Elle inclue une restructuration des routines pour clarifier le role "moniteur" de la routine lscp_main,
une mise à jour de la parametrisation de partitionnement de phase de Lea pour inclure les nuages de couche limite,
ainsi que des corrections des routines de precipitations "poprecip".
Convergence numerique verifiee en prod et debug pour les physiques NPv6.3 et 7.0.1c

Location:

LMDZ6/trunk

Files:

• DefLists/field_def_lmdz.xml (modified) (1 diff)
• libf/phylmd/lmdz_cloud_optics_prop.f90 (modified) (2 diffs)
• libf/phylmd/lmdz_cloudth.f90 (deleted)
• libf/phylmd/lmdz_cloudth_ini.f90 (deleted)
• libf/phylmd/lmdz_lscp.f90 (deleted)
• libf/phylmd/lmdz_lscp_condensation.f90 (modified) (2 diffs)
• libf/phylmd/lmdz_lscp_ini.f90 (modified) (25 diffs)
• libf/phylmd/lmdz_lscp_main.f90 (added)
• libf/phylmd/lmdz_lscp_old.f90 (modified) (1 diff)
• libf/phylmd/lmdz_lscp_precip.f90 (modified) (24 diffs)
• libf/phylmd/lmdz_lscp_tools.f90 (modified) (18 diffs)
• libf/phylmd/phys_local_var_mod.F90 (modified) (3 diffs)
• libf/phylmd/phys_output_ctrlout_mod.F90 (modified) (1 diff)
• libf/phylmd/phys_output_write_mod.F90 (modified) (4 diffs)
• libf/phylmd/physiq_mod.F90 (modified) (3 diffs)
• libf/phylmdiso/lmdz_cloudth.f90 (deleted)
• libf/phylmdiso/lmdz_cloudth_ini.f90 (deleted)
• libf/phylmdiso/lmdz_lscp.f90 (deleted)
• libf/phylmdiso/lmdz_lscp_old.F90 (modified) (1 diff)
• libf/phylmdiso/physiq_mod.F90 (modified) (3 diffs)

LMDZ6/trunk/DefLists/field_def_lmdz.xml

@@ -731 +731 @@
         <field id="sigma2_icefracturb"    long_name="Variance of the diagnostic supersaturation distribution (icefrac_turb)"    unit="-" />
         <field id="mean_icefracturb"     long_name="Mean of the diagnostic supersaturation distribution (icefrac_turb)"    unit="-" />
+        <field id="cldfraliqth"          long_name="Fraction of liquid cloud in thermals"    unit="-" />
+        <field id="sigma2_icefracturbth"    long_name="Variance of the diagnostic supersaturation distribution in thermals (icefrac_turb)"    unit="-" />
+        <field id="mean_icefracturbth"     long_name="Mean of the diagnostic supersaturation distribution in thermals (icefrac_turb)"    unit="-" />
         <field id="rnebcon"    long_name="Convective Cloud Fraction"    unit="-" />
         <field id="rnebls"    long_name="LS Cloud fraction"    unit="-" />

LMDZ6/trunk/libf/phylmd/lmdz_cloud_optics_prop.f90

@@ -176 +176 @@
   reice_pi = 0.
 
-  IF (iflag_t_glace.EQ.0) THEN
+  IF ((.NOT. ok_new_lscp) .AND. iflag_t_glace.EQ.0) THEN
     DO k = 1, klev
       DO i = 1, klon
@@ -202 +202 @@
 
       DO i = 1, klon
-        
+
         IF ((.NOT. ptconv(i,k)) .AND. ok_new_lscp .AND. ok_icefra_lscp) THEN
         ! EV: take the ice fraction directly from the lscp code

LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.f90

@@ -909 +909 @@
 !**********************************************************************************
 
+!**********************************************************************************
 SUBROUTINE deposition_sublimation( &
     qvapincld, qiceincld, temp, qsat, pplay, dtime, &
@@ -1064 +1065 @@
 END SUBROUTINE deposition_sublimation
 
+
+!**********************************************************************************
+
+!**********************************************************************************
+SUBROUTINE  condensation_cloudth(klon,                     &
+&           temp,qt,qt_th,frac_th,zpspsk,play,thetal_th,   &
+&           ratqs,sigma_qtherm,qsth,qsenv,qcloud,ctot,ctotth,ctot_vol,  &
+&           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
+! This routine computes the condensation of clouds in convective boundary layers
+! with thermals assuming two separate distribution of the saturation deficit in 
+! the thermal plumes and in the environment
+! It is based on the work of Arnaud Jam (Jam et al. 2013, BLM)
+! Author : Etienne Vignon (LMDZ/CNRS)
+! Date: February 2025
+! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam
+! IMPORTANT NOTE: we assume iflag_cloudth_vert=7
+!-----------------------------------------------------------------------------------
+
+      use lmdz_lscp_ini,    only: iflag_cloudth_vert,iflag_ratqs,iflag_cloudth_vert_noratqs
+      use lmdz_lscp_ini,    only: vert_alpha, vert_alpha_th ,sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin
+      use lmdz_lscp_ini,    only: RTT, RG, RPI, RD, RV, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th
+
+      IMPLICIT NONE
+
+
+!------------------------------------------------------------------------------
+! Declarations
+!------------------------------------------------------------------------------
+
+! INPUT/OUTPUT
+
+      INTEGER, INTENT(IN)                         :: klon
+      
+
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  temp          ! Temperature (liquid temperature) in the mesh [K] : has seen evap of precip
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  qt            ! total water specific humidity in the mesh [kg/kg]: has seen evap of precip
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  qt_th         ! total water specific humidity in thermals [kg/kg]: has not seen evap of precip
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  thetal_th     ! Liquid potential temperature in thermals [K]: has not seen the evap of precip
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  frac_th       ! Fraction of the mesh covered by thermals [0-1] 
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  zpspsk        ! Exner potential
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  play          ! Pressure of layers [Pa]
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  ratqs         ! Parameter that determines the width of the water distrib     [-]
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  sigma_qtherm  ! Parameter determining the width of the distrib in thermals   [-]
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  qsth          ! Saturation specific humidity in thermals
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  qsenv         ! Saturation specific humidity in environment
+      
+      REAL, DIMENSION(klon),      INTENT(INOUT)      ::  ctot         ! Cloud fraction [0-1]
+      REAL, DIMENSION(klon),      INTENT(INOUT)      ::  ctotth       ! Cloud fraction [0-1] in thermals
+      REAL, DIMENSION(klon),      INTENT(INOUT)      ::  ctot_vol     ! Volume cloud fraction [0-1]
+      REAL, DIMENSION(klon),      INTENT(INOUT)      ::  qcloud       ! In cloud total water content [kg/kg] 
+      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sth    ! mean saturation deficit in thermals
+      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_senv   ! mean saturation deficit in environment
+      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sigmath  ! std of saturation deficit in thermals
+      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sigmaenv ! std of saturation deficit in environment
+
+
+! LOCAL VARIABLES
+
+      INTEGER itap,ind1,l,ig,iter,k 
+      INTEGER iflag_topthermals, niter
+
+      REAL qcth(klon)
+      REAL qcenv(klon)
+      REAL qctot(klon) 
+      REAL cth(klon)
+      REAL cenv(klon)   
+      REAL cth_vol(klon)
+      REAL cenv_vol(klon)
+      REAL qt_env(klon), thetal_env(klon)
+      REAL sqrtpi,sqrt2,sqrt2pi
+      REAL alth,alenv,ath,aenv
+      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
+      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
+      REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
+      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
+      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
+      REAL zdelta,qsatbef,zcor
+      REAL Tbefth(klon), Tbefenv(klon)
+      REAL qlbef
+      REAL dqsatenv(klon), dqsatth(klon)
+      REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
+      REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) 
+      REAL qincloud(klon)
+      REAL alenvl, aenvl
+      REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
+
+
+!------------------------------------------------------------------------------
+! Initialisation
+!------------------------------------------------------------------------------
+
+
+      sqrt2pi=sqrt(2.*rpi)
+      sqrt2=sqrt(2.)
+      sqrtpi=sqrt(rpi)
+
+!-------------------------------------------------------------------------------
+! Thermal fraction calculation and standard deviation of the distribution
+!-------------------------------------------------------------------------------  
+
+! initialisations and calculation of temperature, humidity and saturation specific humidity
+
+cloudth_senv(:) = 0.
+cloudth_sth(:) = 0.
+cloudth_sigmaenv(:) = 0.
+cloudth_sigmath(:) = 0.
+
+
+DO ind1=1,klon
+  
+   Tbefenv(ind1) = temp(ind1)
+   thetal_env(ind1) = Tbefenv(ind1)/zpspsk(ind1)
+   Tbefth(ind1)  = thetal_th(ind1)*zpspsk(ind1)
+   qt_env(ind1)  = (qt(ind1)-frac_th(ind1)*qt_th(ind1))/(1.-frac_th(ind1)) !qt = a*qtth + (1-a)*qtenv 
+
+ENDDO
+
+
+
+DO ind1=1,klon
+
+
+    IF (frac_th(ind1).GT.min_frac_th_cld) THEN !Thermal and environnement
+
+! Environment:
+
+
+        alenv=(RD/RV*RLVTT*qsenv(ind1))/(rd*thetal_env(ind1)**2)     
+        aenv=1./(1.+(alenv*RLVTT/rcpd))                             
+        senv=aenv*(qt_env(ind1)-qsenv(ind1))                           
+     
+
+! Thermals:
+
+
+        alth=(RD/RV*RLVTT*qsth(ind1))/(rd*thetal_th(ind1)**2)       
+        ath=1./(1.+(alth*RLVTT/rcpd))                                                         
+        sth=ath*(qt_th(ind1)-qsth(ind1))                      
+
+
+! Standard deviation of the distributions
+
+           ! environment
+           sigma1s_fraca = (sigma1s_factor**0.5)*(frac_th(ind1)**sigma1s_power) / &
+           &                (1-frac_th(ind1))*((sth-senv)**2)**0.5
+
+           IF (cloudth_ratqsmin>0.) THEN
+             sigma1s_ratqs = cloudth_ratqsmin*qt(ind1)
+           ELSE
+             sigma1s_ratqs = ratqs(ind1)*qt(ind1)
+           ENDIF 
+           sigma1s = sigma1s_fraca + sigma1s_ratqs
+
+           IF (iflag_ratqs.eq.10.or.iflag_ratqs.eq.11) then
+              sigma1s = ratqs(ind1)*qt(ind1)*aenv
+           ENDIF
+
+           ! thermals
+           sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((frac_th(ind1)+0.02)**sigma2s_power))+0.002*qt_th(ind1)
+
+          IF (iflag_ratqs.eq.10.and.sigma_qtherm(ind1).ne.0) then
+             sigma2s = sigma_qtherm(ind1)*ath
+          ENDIF
+
+ 
+! surface cloud fraction 
+
+           deltasenv=aenv*vert_alpha*sigma1s
+           deltasth=ath*vert_alpha_th*sigma2s
+
+           xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
+           xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
+           exp_xenv1 = exp(-1.*xenv1**2)
+           exp_xenv2 = exp(-1.*xenv2**2)
+           xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
+           xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
+           exp_xth1 = exp(-1.*xth1**2)
+           exp_xth2 = exp(-1.*xth2**2)
+           cth(ind1)=0.5*(1.-1.*erf(xth1))
+           cenv(ind1)=0.5*(1.-1.*erf(xenv1))
+           ctot(ind1)=frac_th(ind1)*cth(ind1)+(1.-1.*frac_th(ind1))*cenv(ind1) 
+           ctotth(ind1)=frac_th(ind1)*cth(ind1)
+        
+
+!volume cloud fraction and condensed water 
+
+            !environnement
+
+            IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
+            IntJ_CF=0.5*(1.-1.*erf(xenv2))
+
+            IF (deltasenv .LT. 1.e-10) THEN 
+              qcenv(ind1)=IntJ
+              cenv_vol(ind1)=IntJ_CF
+            ELSE
+              IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
+              IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
+              IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
+              IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
+              IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
+              qcenv(ind1)=IntJ+IntI1+IntI2+IntI3
+              cenv_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
+            ENDIF
+             
+
+
+            !thermals
+
+            IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
+            IntJ_CF=0.5*(1.-1.*erf(xth2))
+     
+            IF (deltasth .LT. 1.e-10) THEN
+              qcth(ind1)=IntJ
+              cth_vol(ind1)=IntJ_CF
+            ELSE
+              IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
+              IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
+              IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
+              IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
+              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
+              qcth(ind1)=IntJ+IntI1+IntI2+IntI3
+              cth_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
+            ENDIF
+
+            ! total 
+
+            qctot(ind1)=frac_th(ind1)*qcth(ind1)+(1.-1.*frac_th(ind1))*qcenv(ind1)
+            ctot_vol(ind1)=frac_th(ind1)*cth_vol(ind1)+(1.-1.*frac_th(ind1))*cenv_vol(ind1)
+
+            IF (cenv(ind1).LT.min_neb_th.and.cth(ind1).LT.min_neb_th) THEN
+                ctot(ind1)=0.
+                ctot_vol(ind1)=0.
+                qcloud(ind1)=qsenv(ind1)
+                qincloud(ind1)=0.
+            ELSE             
+                qincloud(ind1)=qctot(ind1)/ctot(ind1)
+                !to prevent situations with cloud condensed water greater than available total water
+                qincloud(ind1)=min(qincloud(ind1),qt(ind1)/ctot(ind1))
+                ! we assume that water vapor in cloud is qsenv
+                qcloud(ind1)=qincloud(ind1)+qsenv(ind1)
+            ENDIF 
+
+
+
+           ! Outputs used to check the PDFs
+           cloudth_senv(ind1) = senv
+           cloudth_sth(ind1) = sth
+           cloudth_sigmaenv(ind1) = sigma1s
+           cloudth_sigmath(ind1) = sigma2s
+
+      ENDIF       ! selection of grid points concerned by thermals 
+
+
+    ENDDO       !loop on klon
+
+
+RETURN
+
+
+END SUBROUTINE condensation_cloudth
+
+
+!*****************************************************************************************
+!*****************************************************************************************
+! pre-cmip7 routines are below and are becoming obsolete
+!*****************************************************************************************
+!*****************************************************************************************
+
+
+       SUBROUTINE cloudth(ngrid,klev,ind2,  &
+     &           ztv,po,zqta,fraca, & 
+     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
+     &           ratqs,zqs,t, &
+     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
+
+
+      use lmdz_lscp_ini, only: iflag_cloudth_vert,iflag_ratqs
+
+      USE yomcst_mod_h
+      USE yoethf_mod_h
+IMPLICIT NONE
+
+
+!===========================================================================
+! Auteur : Arnaud Octavio Jam (LMD/CNRS)
+! Date : 25 Mai 2010
+! Objet : calcule les valeurs de qc et rneb dans les thermiques
+!===========================================================================
+
+      INCLUDE "FCTTRE.h"
+
+      INTEGER itap,ind1,ind2
+      INTEGER ngrid,klev,klon,l,ig 
+      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
+      
+      REAL ztv(ngrid,klev)
+      REAL po(ngrid)
+      REAL zqenv(ngrid)   
+      REAL zqta(ngrid,klev)
+          
+      REAL fraca(ngrid,klev+1)
+      REAL zpspsk(ngrid,klev)
+      REAL paprs(ngrid,klev+1)
+      REAL pplay(ngrid,klev)
+      REAL ztla(ngrid,klev)
+      REAL zthl(ngrid,klev)
+
+      REAL zqsatth(ngrid,klev)
+      REAL zqsatenv(ngrid,klev)
+      
+      
+      REAL sigma1(ngrid,klev)
+      REAL sigma2(ngrid,klev)
+      REAL qlth(ngrid,klev)
+      REAL qlenv(ngrid,klev)
+      REAL qltot(ngrid,klev) 
+      REAL cth(ngrid,klev)  
+      REAL cenv(ngrid,klev)   
+      REAL ctot(ngrid,klev)
+      REAL rneb(ngrid,klev)
+      REAL t(ngrid,klev)
+      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
+      REAL rdd,cppd,Lv
+      REAL alth,alenv,ath,aenv
+      REAL sth,senv,sigma1s,sigma2s,xth,xenv
+      REAL Tbef,zdelta,qsatbef,zcor
+      REAL qlbef  
+      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
+      
+      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
+      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
+      REAL zqs(ngrid), qcloud(ngrid)
+
+
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Gestion de deux versions de cloudth
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      IF (iflag_cloudth_vert.GE.1) THEN
+      CALL cloudth_vert(ngrid,klev,ind2,  &
+     &           ztv,po,zqta,fraca, & 
+     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
+     &           ratqs,zqs,t)
+      RETURN
+      ENDIF
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+!-------------------------------------------------------------------------------
+! Initialisation des variables r?elles
+!-------------------------------------------------------------------------------
+      sigma1(:,ind2)=0.
+      sigma2(:,ind2)=0.
+      qlth(:,ind2)=0.
+      qlenv(:,ind2)=0.  
+      qltot(:,ind2)=0.
+      rneb(:,ind2)=0.
+      qcloud(:)=0.
+      cth(:,ind2)=0.
+      cenv(:,ind2)=0.
+      ctot(:,ind2)=0.
+      qsatmmussig1=0.
+      qsatmmussig2=0.
+      rdd=287.04
+      cppd=1005.7
+      pi=3.14159 
+      Lv=2.5e6
+      sqrt2pi=sqrt(2.*pi)
+
+
+
+!-------------------------------------------------------------------------------
+! Calcul de la fraction du thermique et des ?cart-types des distributions
+!-------------------------------------------------------------------------------                 
+      do ind1=1,ngrid
+
+      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then 
+
+      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
+
+
+!      zqenv(ind1)=po(ind1)
+      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatenv(ind1,ind2)=qsatbef
+
+
+
+
+      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
+      aenv=1./(1.+(alenv*Lv/cppd))
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 
+
+
+
+
+      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatth(ind1,ind2)=qsatbef
+            
+      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
+      ath=1./(1.+(alth*Lv/cppd))
+      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) 
+      
+      
+
+!------------------------------------------------------------------------------
+! Calcul des ?cart-types pour s
+!------------------------------------------------------------------------------
+
+!      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
+!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2) 
+!       if (paprs(ind1,ind2).gt.90000) then
+!       ratqs(ind1,ind2)=0.002
+!       else 
+!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
+!       endif
+       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
+       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 
+!       sigma1s=ratqs(ind1,ind2)*po(ind1)
+!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003  
+ 
+!------------------------------------------------------------------------------
+! Calcul de l'eau condens?e et de la couverture nuageuse
+!------------------------------------------------------------------------------
+      sqrt2pi=sqrt(2.*pi)
+      xth=sth/(sqrt(2.)*sigma2s)
+      xenv=senv/(sqrt(2.)*sigma1s)
+      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
+      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) 
+      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
+
+      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
+      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
+      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      if (ctot(ind1,ind2).lt.1.e-10) then
+      ctot(ind1,ind2)=0.
+      qcloud(ind1)=zqsatenv(ind1,ind2) 
+
+      else   
+                
+      ctot(ind1,ind2)=ctot(ind1,ind2) 
+      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
+
+      endif                           
+      
+          
+
+
+      else  ! gaussienne environnement seule
+      
+      zqenv(ind1)=po(ind1)
+      Tbef=t(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatenv(ind1,ind2)=qsatbef
+      
+
+!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
+      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
+      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
+      aenv=1./(1.+(alenv*Lv/cppd))
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 
+      
+
+      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
+
+      sqrt2pi=sqrt(2.*pi)
+      xenv=senv/(sqrt(2.)*sigma1s)
+      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
+      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
+      
+      if (ctot(ind1,ind2).lt.1.e-3) then
+      ctot(ind1,ind2)=0.
+      qcloud(ind1)=zqsatenv(ind1,ind2) 
+
+      else   
+                
+      ctot(ind1,ind2)=ctot(ind1,ind2) 
+      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
+
+      endif    
+ 
+ 
+ 
+ 
+ 
+ 
+      endif   
+      enddo
+     
+      return
+!     end
+END SUBROUTINE cloudth
+
+
+
+!===========================================================================
+     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
+     &           ztv,po,zqta,fraca, & 
+     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
+     &           ratqs,zqs,t)
+
+!===========================================================================
+! Auteur : Arnaud Octavio Jam (LMD/CNRS)
+! Date : 25 Mai 2010
+! Objet : calcule les valeurs de qc et rneb dans les thermiques
+!===========================================================================
+
+
+USE yoethf_mod_h
+            use lmdz_lscp_ini, only: iflag_cloudth_vert, vert_alpha
+
+      USE yomcst_mod_h
+IMPLICIT NONE
+
+
+      INCLUDE "FCTTRE.h"
+      
+      INTEGER itap,ind1,ind2
+      INTEGER ngrid,klev,klon,l,ig 
+      
+      REAL ztv(ngrid,klev)
+      REAL po(ngrid)
+      REAL zqenv(ngrid)   
+      REAL zqta(ngrid,klev)
+          
+      REAL fraca(ngrid,klev+1)
+      REAL zpspsk(ngrid,klev)
+      REAL paprs(ngrid,klev+1)
+      REAL pplay(ngrid,klev)
+      REAL ztla(ngrid,klev)
+      REAL zthl(ngrid,klev)
+
+      REAL zqsatth(ngrid,klev)
+      REAL zqsatenv(ngrid,klev)
+      
+      
+      REAL sigma1(ngrid,klev)                                                         
+      REAL sigma2(ngrid,klev)
+      REAL qlth(ngrid,klev)
+      REAL qlenv(ngrid,klev)
+      REAL qltot(ngrid,klev) 
+      REAL cth(ngrid,klev)  
+      REAL cenv(ngrid,klev)   
+      REAL ctot(ngrid,klev)
+      REAL rneb(ngrid,klev)
+      REAL t(ngrid,klev)                                                                  
+      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
+      REAL rdd,cppd,Lv,sqrt2,sqrtpi
+      REAL alth,alenv,ath,aenv
+      REAL sth,senv,sigma1s,sigma2s,xth,xenv
+      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
+      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
+      REAL Tbef,zdelta,qsatbef,zcor
+      REAL qlbef  
+      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
+      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
+      ! (J Jouhaud, JL Dufresne, JB Madeleine)
+      
+      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
+      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
+      REAL zqs(ngrid), qcloud(ngrid)
+
+!------------------------------------------------------------------------------
+! Initialisation des variables r?elles
+!------------------------------------------------------------------------------
+      sigma1(:,ind2)=0.
+      sigma2(:,ind2)=0.
+      qlth(:,ind2)=0.
+      qlenv(:,ind2)=0.  
+      qltot(:,ind2)=0.
+      rneb(:,ind2)=0.
+      qcloud(:)=0.
+      cth(:,ind2)=0.
+      cenv(:,ind2)=0.
+      ctot(:,ind2)=0.
+      qsatmmussig1=0.
+      qsatmmussig2=0.
+      rdd=287.04
+      cppd=1005.7
+      pi=3.14159 
+      Lv=2.5e6
+      sqrt2pi=sqrt(2.*pi)
+      sqrt2=sqrt(2.)
+      sqrtpi=sqrt(pi)
+
+!-------------------------------------------------------------------------------
+! Calcul de la fraction du thermique et des ?cart-types des distributions
+!-------------------------------------------------------------------------------                 
+      do ind1=1,ngrid
+
+      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then 
+
+      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
+
+
+!      zqenv(ind1)=po(ind1)
+      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatenv(ind1,ind2)=qsatbef
+
+
+
+
+      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
+      aenv=1./(1.+(alenv*Lv/cppd))
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 
+
+
+
+
+      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatth(ind1,ind2)=qsatbef
+            
+      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
+      ath=1./(1.+(alth*Lv/cppd))
+      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) 
+      
+      
+
+!------------------------------------------------------------------------------
+! Calcul des ?cart-types pour s
+!------------------------------------------------------------------------------
+
+      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
+      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) 
+!       if (paprs(ind1,ind2).gt.90000) then
+!       ratqs(ind1,ind2)=0.002
+!       else 
+!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
+!       endif
+!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
+!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 
+!       sigma1s=ratqs(ind1,ind2)*po(ind1)
+!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003  
+ 
+!------------------------------------------------------------------------------
+! Calcul de l'eau condens?e et de la couverture nuageuse
+!------------------------------------------------------------------------------
+      sqrt2pi=sqrt(2.*pi)
+      xth=sth/(sqrt(2.)*sigma2s)
+      xenv=senv/(sqrt(2.)*sigma1s)
+      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
+      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) 
+      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
+
+      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
+      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
+      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
+     
+       IF (iflag_cloudth_vert == 1) THEN
+!-------------------------------------------------------------------------------
+!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
+!-------------------------------------------------------------------------------
+!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
+!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
+      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
+      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
+!      deltasenv=aenv*0.01*po(ind1)
+!     deltasth=ath*0.01*zqta(ind1,ind2)    
+      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
+      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
+      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
+      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
+      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
+      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
+      
+      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
+      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) 
+      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
+
+      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
+      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
+      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
+      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
+
+      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
+!      qlenv(ind1,ind2)=IntJ 
+!      print*, qlenv(ind1,ind2),'VERIF EAU'
+
+
+      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
+!      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
+!      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
+      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
+      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
+      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
+      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
+!      qlth(ind1,ind2)=IntJ
+!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
+      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
+
+      ELSE IF (iflag_cloudth_vert == 2) THEN
+
+!-------------------------------------------------------------------------------
+!  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
+!-------------------------------------------------------------------------------
+!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
+!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
+!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
+!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
+      deltasenv=aenv*vert_alpha*sigma1s
+      deltasth=ath*vert_alpha*sigma2s
+      
+      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
+      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
+      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
+      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
+!     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
+!     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
+      
+      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
+      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
+      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
+
+      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
+      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
+      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
+      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))
+
+!      IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
+!      IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
+!      IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
+
+      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
+!      qlenv(ind1,ind2)=IntJ 
+!      print*, qlenv(ind1,ind2),'VERIF EAU'
+
+      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
+      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
+      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
+      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
+      
+      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
+!      qlth(ind1,ind2)=IntJ
+!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
+      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
+      
+
+
+
+      ENDIF ! of if (iflag_cloudth_vert==1 or 2)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
+      ctot(ind1,ind2)=0.
+      qcloud(ind1)=zqsatenv(ind1,ind2) 
+
+      else 
+                
+      ctot(ind1,ind2)=ctot(ind1,ind2) 
+      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
+!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
+!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) 
+
+      endif  
+                       
+      
+          
+!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
+
+
+      else  ! gaussienne environnement seule
+      
+      zqenv(ind1)=po(ind1)
+      Tbef=t(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatenv(ind1,ind2)=qsatbef
+      
+
+!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
+      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
+      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
+      aenv=1./(1.+(alenv*Lv/cppd))
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 
+      
+
+      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
+
+      sqrt2pi=sqrt(2.*pi)
+      xenv=senv/(sqrt(2.)*sigma1s)
+      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
+      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
+      
+      if (ctot(ind1,ind2).lt.1.e-3) then
+      ctot(ind1,ind2)=0.
+      qcloud(ind1)=zqsatenv(ind1,ind2) 
+
+      else   
+                
+      ctot(ind1,ind2)=ctot(ind1,ind2) 
+      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
+
+      endif    
+ 
+ 
+ 
+ 
+ 
+ 
+      endif   
+      enddo
+     
+      return
+!     end
+END SUBROUTINE cloudth_vert
+
+
+
+
+       SUBROUTINE cloudth_v3(ngrid,klev,ind2,  &
+     &           ztv,po,zqta,fraca, & 
+     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
+     &           ratqs,sigma_qtherm,zqs,t, &
+     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
+
+      use lmdz_lscp_ini, only: iflag_cloudth_vert
+
+      USE yomcst_mod_h
+      USE yoethf_mod_h
+IMPLICIT NONE
+
+
+!===========================================================================
+! Author : Arnaud Octavio Jam (LMD/CNRS)
+! Date : 25 Mai 2010
+! Objet : calcule les valeurs de qc et rneb dans les thermiques
+!===========================================================================
+      INCLUDE "FCTTRE.h"
+
+      integer, intent(in) :: ind2
+      integer, intent(in) :: ngrid,klev
+      
+      real, dimension(ngrid,klev), intent(in) :: ztv
+      real, dimension(ngrid), intent(in) :: po
+      real, dimension(ngrid,klev), intent(in) :: zqta
+      real, dimension(ngrid,klev+1), intent(in) :: fraca
+      real, dimension(ngrid), intent(out) :: qcloud
+      real, dimension(ngrid,klev), intent(out) :: ctot
+      real, dimension(ngrid,klev), intent(out) :: ctot_vol
+      real, dimension(ngrid,klev), intent(in) :: zpspsk
+      real, dimension(ngrid,klev+1), intent(in) :: paprs
+      real, dimension(ngrid,klev), intent(in) :: pplay
+      real, dimension(ngrid,klev), intent(in) :: ztla
+      real, dimension(ngrid,klev), intent(inout) :: zthl
+      real, dimension(ngrid,klev), intent(in) :: ratqs,sigma_qtherm
+      real, dimension(ngrid), intent(in) :: zqs
+      real, dimension(ngrid,klev), intent(in) :: t
+      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
+
+
+      REAL zqenv(ngrid)   
+      REAL zqsatth(ngrid,klev)
+      REAL zqsatenv(ngrid,klev)
+      
+      REAL sigma1(ngrid,klev)                                                         
+      REAL sigma2(ngrid,klev)
+      REAL qlth(ngrid,klev)
+      REAL qlenv(ngrid,klev)
+      REAL qltot(ngrid,klev) 
+      REAL cth(ngrid,klev)
+      REAL cenv(ngrid,klev)   
+      REAL cth_vol(ngrid,klev)
+      REAL cenv_vol(ngrid,klev)
+      REAL rneb(ngrid,klev)      
+      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
+      REAL rdd,cppd,Lv
+      REAL alth,alenv,ath,aenv
+      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
+      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
+      REAL Tbef,zdelta,qsatbef,zcor
+      REAL qlbef  
+      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
+      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
+
+
+      INTEGER :: ind1,l, ig
+
+      IF (iflag_cloudth_vert.GE.1) THEN
+      CALL cloudth_vert_v3(ngrid,klev,ind2,  &
+     &           ztv,po,zqta,fraca, & 
+     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
+     &           ratqs,sigma_qtherm,zqs,t, &
+     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
+      RETURN
+      ENDIF
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+!-------------------------------------------------------------------------------
+! Initialisation des variables r?elles
+!-------------------------------------------------------------------------------
+      sigma1(:,ind2)=0.
+      sigma2(:,ind2)=0.
+      qlth(:,ind2)=0.
+      qlenv(:,ind2)=0.  
+      qltot(:,ind2)=0.
+      rneb(:,ind2)=0.
+      qcloud(:)=0.
+      cth(:,ind2)=0.
+      cenv(:,ind2)=0.
+      ctot(:,ind2)=0.
+      cth_vol(:,ind2)=0.
+      cenv_vol(:,ind2)=0.
+      ctot_vol(:,ind2)=0.
+      qsatmmussig1=0.
+      qsatmmussig2=0.
+      rdd=287.04
+      cppd=1005.7
+      pi=3.14159 
+      Lv=2.5e6
+      sqrt2pi=sqrt(2.*pi)
+      sqrt2=sqrt(2.)
+      sqrtpi=sqrt(pi)
+
+
+!-------------------------------------------------------------------------------
+! Cloud fraction in the thermals and standard deviation of the PDFs
+!-------------------------------------------------------------------------------                 
+      do ind1=1,ngrid
+
+      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then 
+
+      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
+
+      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatenv(ind1,ind2)=qsatbef
+
+
+      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
+      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
+
+!po = qt de l'environnement ET des thermique
+!zqenv = qt environnement
+!zqsatenv = qsat environnement
+!zthl = Tl environnement
+
+
+      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatth(ind1,ind2)=qsatbef
+            
+      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
+      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
+      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
+      
+!zqta = qt thermals
+!zqsatth = qsat thermals
+!ztla = Tl thermals
+
+!------------------------------------------------------------------------------
+! s standard deviations
+!------------------------------------------------------------------------------
+
+!     tests
+!     sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
+!     sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
+!     sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
+!     final option
+      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
+      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
+ 
+!------------------------------------------------------------------------------
+! Condensed water and cloud cover
+!------------------------------------------------------------------------------
+      xth=sth/(sqrt2*sigma2s)
+      xenv=senv/(sqrt2*sigma1s)
+      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))       !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
+      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
+      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
+      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
+
+      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
+      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
+      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
+
+      if (ctot(ind1,ind2).lt.1.e-10) then
+      ctot(ind1,ind2)=0.
+      qcloud(ind1)=zqsatenv(ind1,ind2) 
+      else
+      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
+      endif
+
+      else  ! Environnement only, follow the if l.110
+      
+      zqenv(ind1)=po(ind1)
+      Tbef=t(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatenv(ind1,ind2)=qsatbef
+
+!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
+      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
+      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
+      aenv=1./(1.+(alenv*Lv/cppd))
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
+      
+      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
+
+      xenv=senv/(sqrt2*sigma1s)
+      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
+      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
+      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
+
+      if (ctot(ind1,ind2).lt.1.e-3) then
+      ctot(ind1,ind2)=0.
+      qcloud(ind1)=zqsatenv(ind1,ind2)
+      else   
+      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
+      endif
+
+
+      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
+      enddo       ! from the loop on ngrid l.108
+      return
+!     end
+END SUBROUTINE cloudth_v3
+
+
+
+!===========================================================================
+     SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2,  &
+     &           ztv,po,zqta,fraca, & 
+     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
+     &           ratqs,sigma_qtherm,zqs,t, &
+     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
+
+!===========================================================================
+! Auteur : Arnaud Octavio Jam (LMD/CNRS)
+! Date : 25 Mai 2010
+! Objet : calcule les valeurs de qc et rneb dans les thermiques
+!===========================================================================
+
+      use yoethf_mod_h
+      use lmdz_lscp_ini, only : iflag_cloudth_vert,iflag_ratqs
+      use lmdz_lscp_ini, only : vert_alpha,vert_alpha_th, sigma1s_factor, sigma1s_power , sigma2s_factor , sigma2s_power , cloudth_ratqsmin , iflag_cloudth_vert_noratqs
+
+      USE yomcst_mod_h
+IMPLICIT NONE
+
+
+
+
+      INCLUDE "FCTTRE.h"
+      
+      INTEGER itap,ind1,ind2
+      INTEGER ngrid,klev,klon,l,ig 
+      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
+      
+      REAL ztv(ngrid,klev)
+      REAL po(ngrid)
+      REAL zqenv(ngrid)   
+      REAL zqta(ngrid,klev)
+          
+      REAL fraca(ngrid,klev+1)
+      REAL zpspsk(ngrid,klev)
+      REAL paprs(ngrid,klev+1)
+      REAL pplay(ngrid,klev)
+      REAL ztla(ngrid,klev)
+      REAL zthl(ngrid,klev)
+
+      REAL zqsatth(ngrid,klev)
+      REAL zqsatenv(ngrid,klev)
+      
+      REAL sigma1(ngrid,klev)                                                         
+      REAL sigma2(ngrid,klev)
+      REAL qlth(ngrid,klev)
+      REAL qlenv(ngrid,klev)
+      REAL qltot(ngrid,klev) 
+      REAL cth(ngrid,klev)
+      REAL cenv(ngrid,klev)   
+      REAL ctot(ngrid,klev)
+      REAL cth_vol(ngrid,klev)
+      REAL cenv_vol(ngrid,klev)
+      REAL ctot_vol(ngrid,klev)
+      REAL rneb(ngrid,klev)
+      REAL t(ngrid,klev)                                                                  
+      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
+      REAL rdd,cppd,Lv
+      REAL alth,alenv,ath,aenv
+      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
+      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
+      REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
+      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
+      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
+      REAL Tbef,zdelta,qsatbef,zcor
+      REAL qlbef  
+      REAL ratqs(ngrid,klev),sigma_qtherm(ngrid,klev) ! determine la largeur de distribution de vapeur
+      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
+      ! (J Jouhaud, JL Dufresne, JB Madeleine)
+
+      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
+      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
+      REAL zqs(ngrid), qcloud(ngrid)
+
+      REAL rhodz(ngrid,klev)
+      REAL zrho(ngrid,klev)
+      REAL dz(ngrid,klev)
+
+      DO ind1 = 1, ngrid
+        !Layer calculation
+        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2
+        zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3
+        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre
+      END DO
+
+!------------------------------------------------------------------------------
+! Initialize
+!------------------------------------------------------------------------------
+
+      sigma1(:,ind2)=0.
+      sigma2(:,ind2)=0.
+      qlth(:,ind2)=0.
+      qlenv(:,ind2)=0.  
+      qltot(:,ind2)=0.
+      rneb(:,ind2)=0.
+      qcloud(:)=0.
+      cth(:,ind2)=0.
+      cenv(:,ind2)=0.
+      ctot(:,ind2)=0.
+      cth_vol(:,ind2)=0.
+      cenv_vol(:,ind2)=0.
+      ctot_vol(:,ind2)=0.
+      qsatmmussig1=0.
+      qsatmmussig2=0.
+      rdd=287.04
+      cppd=1005.7
+      pi=3.14159 
+      Lv=2.5e6
+      sqrt2pi=sqrt(2.*pi)
+      sqrt2=sqrt(2.)
+      sqrtpi=sqrt(pi)
+
+
+
+!-------------------------------------------------------------------------------
+! Calcul de la fraction du thermique et des ecart-types des distributions
+!-------------------------------------------------------------------------------                 
+      do ind1=1,ngrid
+
+      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
+
+      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
+
+
+      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatenv(ind1,ind2)=qsatbef
+
+
+      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
+      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
+
+!zqenv = qt environnement
+!zqsatenv = qsat environnement
+!zthl = Tl environnement
+
+
+      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatth(ind1,ind2)=qsatbef
+            
+      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
+      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
+      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
+      
+      
+!zqta = qt thermals
+!zqsatth = qsat thermals
+!ztla = Tl thermals
+!------------------------------------------------------------------------------
+! s standard deviation
+!------------------------------------------------------------------------------
+
+      sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
+     &                (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
+!     sigma1s_fraca = (1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
+      IF (cloudth_ratqsmin>0.) THEN
+         sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
+      ELSE
+         sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
+      ENDIF
+      sigma1s = sigma1s_fraca + sigma1s_ratqs
+      sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
+      IF (iflag_ratqs.eq.10.or.iflag_ratqs.eq.11) then
+         sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
+         IF (iflag_ratqs.eq.10.and.sigma_qtherm(ind1,ind2).ne.0) then
+            sigma2s = sigma_qtherm(ind1,ind2)*ath
+         ENDIF
+      ENDIF
+      
+!      tests
+!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
+!      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
+!      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) 
+!      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
+!       if (paprs(ind1,ind2).gt.90000) then
+!       ratqs(ind1,ind2)=0.002
+!       else 
+!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
+!       endif
+!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
+!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 
+!       sigma1s=ratqs(ind1,ind2)*po(ind1)
+!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003  
+ 
+       IF (iflag_cloudth_vert == 1) THEN
+!-------------------------------------------------------------------------------
+!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
+!-------------------------------------------------------------------------------
+
+      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
+      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
+
+      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
+      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
+      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
+      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
+      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
+      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
+      
+      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
+      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) 
+      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
+
+      ! Environment
+      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
+      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
+      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
+      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
+
+      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
+
+      ! Thermal
+      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
+      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
+      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
+      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
+      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
+      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
+
+      ELSE IF (iflag_cloudth_vert >= 3) THEN
+      IF (iflag_cloudth_vert < 5) THEN
+!-------------------------------------------------------------------------------
+!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
+!-------------------------------------------------------------------------------
+!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
+!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
+!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
+!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
+      IF (iflag_cloudth_vert == 3) THEN
+        deltasenv=aenv*vert_alpha*sigma1s
+        deltasth=ath*vert_alpha_th*sigma2s
+      ELSE IF (iflag_cloudth_vert == 4) THEN
+        IF (iflag_cloudth_vert_noratqs == 1) THEN
+          deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
+          deltasth=vert_alpha_th*sigma2s
+        ELSE
+          deltasenv=vert_alpha*sigma1s
+          deltasth=vert_alpha_th*sigma2s
+        ENDIF
+      ENDIF
+      
+      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
+      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
+      exp_xenv1 = exp(-1.*xenv1**2)
+      exp_xenv2 = exp(-1.*xenv2**2)
+      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
+      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
+      exp_xth1 = exp(-1.*xth1**2)
+      exp_xth2 = exp(-1.*xth2**2)
+      
+      !CF_surfacique
+      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
+      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
+      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) 
+
+
+      !CF_volumique & eau condense
+      !environnement
+      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
+      IntJ_CF=0.5*(1.-1.*erf(xenv2))
+      if (deltasenv .lt. 1.e-10) then 
+      qlenv(ind1,ind2)=IntJ
+      cenv_vol(ind1,ind2)=IntJ_CF
+      else
+      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
+      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
+      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
+      IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
+      IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
+      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
+      cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
+      endif
+
+      !thermique
+      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
+      IntJ_CF=0.5*(1.-1.*erf(xth2))
+      if (deltasth .lt. 1.e-10) then
+      qlth(ind1,ind2)=IntJ
+      cth_vol(ind1,ind2)=IntJ_CF
+      else
+      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
+      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
+      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
+      IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
+      IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
+      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
+      cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
+      endif
+
+      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
+      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
+
+      ELSE IF (iflag_cloudth_vert == 5) THEN
+         sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) &
+              /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) &
+              +ratqs(ind1,ind2)*po(ind1) !Environment
+      sigma2s=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.02)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2)                   !Thermals
+      !sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
+      !sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 
+      xth=sth/(sqrt(2.)*sigma2s)
+      xenv=senv/(sqrt(2.)*sigma1s)
+
+      !Volumique
+      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
+      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
+      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
+      !print *,'jeanjean_CV=',ctot_vol(ind1,ind2)
+
+      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth_vol(ind1,ind2))
+      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv_vol(ind1,ind2))  
+      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
+
+      !Surfacique
+      !Neggers
+      !beta=0.0044
+      !inverse_rho=1.+beta*dz(ind1,ind2)
+      !print *,'jeanjean : beta=',beta
+      !cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
+      !cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
+      !ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
+
+      !Brooks
+      a_Brooks=0.6694
+      b_Brooks=0.1882
+      A_Maj_Brooks=0.1635 !-- sans shear
+      !A_Maj_Brooks=0.17   !-- ARM LES
+      !A_Maj_Brooks=0.18   !-- RICO LES
+      !A_Maj_Brooks=0.19   !-- BOMEX LES
+      Dx_Brooks=200000.
+      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
+      !print *,'jeanjean_f=',f_Brooks
+
+      cth(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
+      cenv(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
+      ctot(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
+      !print *,'JJ_ctot_1',ctot(ind1,ind2)
+
+
+
+
+
+      ENDIF ! of if (iflag_cloudth_vert<5)
+      ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
+
+!      if (ctot(ind1,ind2).lt.1.e-10) then
+      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
+      ctot(ind1,ind2)=0.
+      ctot_vol(ind1,ind2)=0.
+      qcloud(ind1)=zqsatenv(ind1,ind2)
+
+      else
+                
+      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
+!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
+!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) 
+
+      endif  
+
+      else  ! gaussienne environnement seule
+      
+
+      zqenv(ind1)=po(ind1)
+      Tbef=t(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatenv(ind1,ind2)=qsatbef
+      
+
+!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
+      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
+      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
+      aenv=1./(1.+(alenv*Lv/cppd))
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
+      sth=0.
+      
+
+      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
+      sigma2s=0.
+
+      sqrt2pi=sqrt(2.*pi)
+      xenv=senv/(sqrt(2.)*sigma1s)
+      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
+      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
+      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
+      
+      if (ctot(ind1,ind2).lt.1.e-3) then
+      ctot(ind1,ind2)=0.
+      qcloud(ind1)=zqsatenv(ind1,ind2)
+
+      else   
+                
+!      ctot(ind1,ind2)=ctot(ind1,ind2) 
+      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
+
+      endif  
+ 
+
+
+
+      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
+      ! Outputs used to check the PDFs
+      cloudth_senv(ind1,ind2) = senv
+      cloudth_sth(ind1,ind2) = sth
+      cloudth_sigmaenv(ind1,ind2) = sigma1s
+      cloudth_sigmath(ind1,ind2) = sigma2s
+
+      enddo       ! from the loop on ngrid l.333
+      return
+!     end
+END SUBROUTINE cloudth_vert_v3
+!
+
+
+
+
+
+
+
+
+
+
+
+       SUBROUTINE cloudth_v6(ngrid,klev,ind2,  &
+     &           ztv,po,zqta,fraca, & 
+     &           qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
+     &           ratqs,zqs,T, &
+     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
+
+      USE yoethf_mod_h
+      USE lmdz_lscp_ini, only: iflag_cloudth_vert
+
+      USE yomcst_mod_h
+IMPLICIT NONE
+
+
+
+      INCLUDE "FCTTRE.h"
+
+
+        !Domain variables
+      INTEGER ngrid !indice Max lat-lon
+      INTEGER klev  !indice Max alt
+      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
+      INTEGER ind1  !indice in [1:ngrid]
+      INTEGER ind2  !indice in [1:klev]
+        !thermal plume fraction
+      REAL fraca(ngrid,klev+1)   !thermal plumes fraction in the gridbox
+        !temperatures
+      REAL T(ngrid,klev)       !temperature
+      REAL zpspsk(ngrid,klev)  !factor (p/p0)**kappa (used for potential variables)
+      REAL ztv(ngrid,klev)     !potential temperature (voir thermcell_env.F90)
+      REAL ztla(ngrid,klev)    !liquid temperature in the thermals (Tl_th)
+      REAL zthl(ngrid,klev)    !liquid temperature in the environment (Tl_env)
+        !pressure
+      REAL paprs(ngrid,klev+1)   !pressure at the interface of levels
+      REAL pplay(ngrid,klev)     !pressure at the middle of the level
+        !humidity
+      REAL ratqs(ngrid,klev)   !width of the total water subgrid-scale distribution
+      REAL po(ngrid)           !total water (qt)
+      REAL zqenv(ngrid)        !total water in the environment (qt_env)
+      REAL zqta(ngrid,klev)    !total water in the thermals (qt_th)
+      REAL zqsatth(ngrid,klev)   !water saturation level in the thermals (q_sat_th)
+      REAL zqsatenv(ngrid,klev)  !water saturation level in the environment (q_sat_env)
+      REAL qlth(ngrid,klev)    !condensed water in the thermals
+      REAL qlenv(ngrid,klev)   !condensed water in the environment
+      REAL qltot(ngrid,klev)   !condensed water in the gridbox
+        !cloud fractions
+      REAL cth_vol(ngrid,klev)   !cloud fraction by volume in the thermals
+      REAL cenv_vol(ngrid,klev)  !cloud fraction by volume in the environment
+      REAL ctot_vol(ngrid,klev)  !cloud fraction by volume in the gridbox
+      REAL cth_surf(ngrid,klev)  !cloud fraction by surface in the thermals
+      REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment  
+      REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
+        !PDF of saturation deficit variables
+      REAL rdd,cppd,Lv
+      REAL Tbef,zdelta,qsatbef,zcor
+      REAL alth,alenv,ath,aenv
+      REAL sth,senv              !saturation deficits in the thermals and environment
+      REAL sigma_env,sigma_th    !standard deviations of the biGaussian PDF
+        !cloud fraction variables
+      REAL xth,xenv
+      REAL inverse_rho,beta                                  !Neggers et al. (2011) method
+      REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
+        !Incloud total water variables
+      REAL zqs(ngrid)    !q_sat
+      REAL qcloud(ngrid) !eau totale dans le nuage
+        !Some arithmetic variables
+      REAL  pi,sqrt2,sqrt2pi
+        !Depth of the layer
+      REAL dz(ngrid,klev)    !epaisseur de la couche en metre
+      REAL rhodz(ngrid,klev)
+      REAL zrho(ngrid,klev)
+      DO ind1 = 1, ngrid
+        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
+        zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd          ![kg/m3]
+        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2)            ![m]
+      END DO
+
+!------------------------------------------------------------------------------
+! Initialization
+!------------------------------------------------------------------------------
+      qlth(:,ind2)=0.
+      qlenv(:,ind2)=0.  
+      qltot(:,ind2)=0.
+      cth_vol(:,ind2)=0.
+      cenv_vol(:,ind2)=0.
+      ctot_vol(:,ind2)=0.
+      cth_surf(:,ind2)=0.
+      cenv_surf(:,ind2)=0.
+      ctot_surf(:,ind2)=0.
+      qcloud(:)=0.
+      rdd=287.04
+      cppd=1005.7
+      pi=3.14159 
+      Lv=2.5e6
+      sqrt2=sqrt(2.)
+      sqrt2pi=sqrt(2.*pi)
+
+
+      DO ind1=1,ngrid
+!-------------------------------------------------------------------------------
+!Both thermal and environment in the gridbox
+!-------------------------------------------------------------------------------
+      IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN
+        !--------------------------------------------
+        !calcul de qsat_env
+        !--------------------------------------------
+      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatenv(ind1,ind2)=qsatbef
+        !--------------------------------------------
+        !calcul de s_env
+        !--------------------------------------------
+      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
+      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam
+        !--------------------------------------------
+        !calcul de qsat_th
+        !--------------------------------------------
+      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatth(ind1,ind2)=qsatbef
+        !--------------------------------------------
+        !calcul de s_th  
+        !--------------------------------------------
+      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84 these Arnaud Jam
+      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84 these Arnaud Jam
+      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84 these Arnaud Jam
+        !--------------------------------------------
+        !calcul standard deviations bi-Gaussian PDF
+        !--------------------------------------------
+      sigma_th=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.01)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2)
+      sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) &
+           /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) &
+           +ratqs(ind1,ind2)*po(ind1)
+      xth=sth/(sqrt2*sigma_th)
+      xenv=senv/(sqrt2*sigma_env)
+        !--------------------------------------------
+        !Cloud fraction by volume CF_vol
+        !--------------------------------------------
+      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
+      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
+      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
+        !-------------------------------------------- 
+        !Condensed water qc
+        !--------------------------------------------
+      qlth(ind1,ind2)=sigma_th*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth_vol(ind1,ind2))
+      qlenv(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv_vol(ind1,ind2))  
+      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
+        !--------------------------------------------
+        !Cloud fraction by surface CF_surf
+        !--------------------------------------------
+        !Method Neggers et al. (2011) : ok for cumulus clouds only
+      !beta=0.0044 (Jouhaud et al.2018)
+      !inverse_rho=1.+beta*dz(ind1,ind2)
+      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
+        !Method Brooks et al. (2005) : ok for all types of clouds
+      a_Brooks=0.6694
+      b_Brooks=0.1882
+      A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
+      Dx_Brooks=200000.   !-- si l'on considere des mailles de 200km de cote
+      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
+      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
+        !--------------------------------------------
+        !Incloud Condensed water qcloud
+        !--------------------------------------------
+      if (ctot_surf(ind1,ind2) .lt. 1.e-10) then
+      ctot_vol(ind1,ind2)=0.
+      ctot_surf(ind1,ind2)=0.
+      qcloud(ind1)=zqsatenv(ind1,ind2)
+      else
+      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
+      endif
+
+
+
+!-------------------------------------------------------------------------------
+!Environment only in the gridbox
+!-------------------------------------------------------------------------------
+      ELSE
+        !--------------------------------------------
+        !calcul de qsat_env
+        !--------------------------------------------
+      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatenv(ind1,ind2)=qsatbef
+        !--------------------------------------------
+        !calcul de s_env
+        !--------------------------------------------
+      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
+      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam 
+        !--------------------------------------------
+        !calcul standard deviations Gaussian PDF
+        !--------------------------------------------
+      zqenv(ind1)=po(ind1)
+      sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
+      xenv=senv/(sqrt2*sigma_env)
+        !--------------------------------------------
+        !Cloud fraction by volume CF_vol
+        !--------------------------------------------
+      ctot_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
+        !--------------------------------------------
+        !Condensed water qc
+        !--------------------------------------------
+      qltot(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*ctot_vol(ind1,ind2))
+        !--------------------------------------------
+        !Cloud fraction by surface CF_surf
+        !--------------------------------------------
+        !Method Neggers et al. (2011) : ok for cumulus clouds only
+      !beta=0.0044 (Jouhaud et al.2018)
+      !inverse_rho=1.+beta*dz(ind1,ind2)
+      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
+        !Method Brooks et al. (2005) : ok for all types of clouds
+      a_Brooks=0.6694
+      b_Brooks=0.1882
+      A_Maj_Brooks=0.1635 !-- sans dependence au shear
+      Dx_Brooks=200000.
+      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
+      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
+        !--------------------------------------------
+        !Incloud Condensed water qcloud
+        !--------------------------------------------
+      if (ctot_surf(ind1,ind2) .lt. 1.e-8) then
+      ctot_vol(ind1,ind2)=0.
+      ctot_surf(ind1,ind2)=0.
+      qcloud(ind1)=zqsatenv(ind1,ind2)
+      else
+      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
+      endif
+
+
+      END IF  ! From the separation (thermal/envrionnement) et (environnement only)
+
+      ! Outputs used to check the PDFs
+      cloudth_senv(ind1,ind2) = senv
+      cloudth_sth(ind1,ind2) = sth
+      cloudth_sigmaenv(ind1,ind2) = sigma_env
+      cloudth_sigmath(ind1,ind2) = sigma_th
+
+      END DO  ! From the loop on ngrid
+      return
+
+END SUBROUTINE cloudth_v6
+
+
+
 END MODULE lmdz_lscp_condensation

LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.f90

@@ -9 +9 @@
   !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RV, RG, RPI, EPS_W)
   
-  REAL, SAVE, PROTECTED :: seuil_neb=0.001      ! cloud fraction threshold: a cloud can precipitate when exceeded
+  INTEGER, SAVE, PROTECTED :: iflag_ratqs        ! control of ratqs option
+  !$OMP THREADPRIVATE(iflag_ratqs)
+  
+  REAL, SAVE, PROTECTED :: seuil_neb=0.001       ! cloud fraction threshold: a cloud can precipitate when exceeded
   !$OMP THREADPRIVATE(seuil_neb)
 
@@ -67 +70 @@
   !$OMP THREADPRIVATE(iflag_t_glace)
 
-  INTEGER, SAVE, PROTECTED :: iflag_cloudth_vert=0          ! option for determining cloud fraction and content in convective boundary layers
-  !$OMP THREADPRIVATE(iflag_cloudth_vert)
-
   INTEGER, SAVE, PROTECTED :: iflag_gammasat=0              ! which threshold for homogeneous nucleation below -40oC
   !$OMP THREADPRIVATE(iflag_gammasat)
@@ -133 +133 @@
   !$OMP THREADPRIVATE(expo_sub)
 
-  REAL, SAVE, PROTECTED :: cice_velo=1.645                  ! factor in the ice fall velocity formulation
+  REAL, SAVE, PROTECTED :: cice_velo=1.645                  ! factor in the ice fall velocity formulation. It is half the value of
+                                                            ! Heymsfield and Donner 1990 to concur with previous LMDZ versions
   !$OMP THREADPRIVATE(cice_velo)
 
@@ -205 +206 @@
   !--End of the parameters for condensation and ice supersaturation
 
-  !--Parameters for poprecip
+  !--Parameters for poprecip and cloud phase
   LOGICAL, SAVE, PROTECTED :: ok_poprecip=.FALSE.           ! use the processes-oriented formulation of precipitations
   !$OMP THREADPRIVATE(ok_poprecip)
@@ -212 +213 @@
   !$OMP THREADPRIVATE(ok_corr_vap_evasub)
 
+  LOGICAL, SAVE, PROTECTED :: ok_growth_precip_deposition=.FALSE. ! allows growth of snowfall through vapor deposition in supersat. regions
+  !$OMP THREADPRIVATE(ok_growth_precip_deposition)
+
   REAL, SAVE, PROTECTED :: cld_lc_lsc_snow=2.e-5            ! snow autoconversion coefficient, stratiform. default from  Chaboureau and PInty 2006
   !$OMP THREADPRIVATE(cld_lc_lsc_snow)
@@ -233 +237 @@
   !$OMP THREADPRIVATE(gamma_snwretro)
 
+  REAL, SAVE, PROTECTED :: gamma_mixth = 1.                 ! Tuning coeff for mixing with thermals/env in lscp_icefrac_turb [-]
+  !$OMP THREADPRIVATE(gamma_mixth)
+
   REAL, SAVE, PROTECTED :: gamma_taud = 1.                  ! Tuning coeff for Lagrangian decorrelation timescale in lscp_icefrac_turb [-]
   !$OMP THREADPRIVATE(gamma_taud)
@@ -254 +261 @@
   !$OMP THREADPRIVATE(rho_rain)
 
-  REAL, SAVE, PROTECTED :: rho_ice=920.                     ! Ice density  [kg/m3]
+  REAL, SAVE, PROTECTED :: rho_ice=920.                     ! Ice crystal density (assuming spherical geometry) [kg/m3]
   !$OMP THREADPRIVATE(rho_ice)
 
@@ -268 +275 @@
   REAL, SAVE, PROTECTED :: tau_auto_snow_max=1000.          ! Snow autoconversion minimal timescale (when only ice) [s]
   !$OMP THREADPRIVATE(tau_auto_snow_max)
+
+  REAL, SAVE, PROTECTED :: expo_tau_auto_snow=0.1          ! Snow autoconversion timescale exponent for icefrac dependency
+  !$OMP THREADPRIVATE(expo_tau_auto_snow)
 
   REAL, SAVE, PROTECTED :: eps=1.E-10                       ! Treshold 0 [-]
@@ -297 +307 @@
   !--End of the parameters for poprecip
 
-! Two parameters used for lmdz_lscp_old only
+  ! Parameters for cloudth routines
+  LOGICAL, SAVE, PROTECTED :: ok_lscp_mergecond=.false.     ! more consistent condensation stratiform and shallow convective clouds
+  !$OMP THREADPRIVATE(ok_lscp_mergecond)
+  
+  INTEGER, SAVE, PROTECTED :: iflag_cloudth_vert=0          ! option for determining cloud fraction and content in convective boundary layers
+  !$OMP THREADPRIVATE(iflag_cloudth_vert)
+
+  INTEGER, SAVE, PROTECTED :: iflag_cloudth_vert_noratqs=0  ! option to control the width of gaussian distrib in a specific case
+  !$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs)
+
+  REAL, SAVE, PROTECTED :: cloudth_ratqsmin=-1.             ! minimum ratqs in cloudth
+  !$OMP THREADPRIVATE(cloudth_ratqsmin)
+
+  REAL, SAVE, PROTECTED :: sigma1s_factor=1.1               ! factor for standard deviation of gaussian distribution of environment
+  !$OMP THREADPRIVATE(sigma1s_factor)
+
+  REAL, SAVE, PROTECTED :: sigma2s_factor=0.09              ! factor for standard deviation of gaussian distribution of thermals
+  !$OMP THREADPRIVATE(sigma2s_factor)
+
+
+  REAL, SAVE, PROTECTED :: sigma1s_power=0.6                ! exponent for standard deviation of gaussian distribution of environment
+  !$OMP THREADPRIVATE(sigma1s_power)
+    
+  REAL, SAVE, PROTECTED :: sigma2s_power=0.5                ! exponent for standard deviation of gaussian distribution of thermals
+  !$OMP THREADPRIVATE(sigma2s_power)
+
+  REAL, SAVE, PROTECTED :: vert_alpha=0.5                   ! tuning coefficient for standard deviation of gaussian distribution of thermals
+  !$OMP THREADPRIVATE(vert_alpha)
+
+  REAL, SAVE, PROTECTED :: vert_alpha_th=0.5                ! tuning coefficient for standard deviation of gaussian distribution of thermals
+  !$OMP THREADPRIVATE(vert_alpha_th)
+  ! End of parameters for cloudth routines
+
+  ! Two parameters used for lmdz_lscp_old only
   INTEGER, SAVE, PROTECTED :: iflag_oldbug_fisrtilp=0, fl_cor_ebil 
   !$OMP THREADPRIVATE(iflag_oldbug_fisrtilp,fl_cor_ebil)
@@ -303 +346 @@
 CONTAINS
 
-SUBROUTINE lscp_ini(dtime,lunout_in,prt_level_in,ok_ice_supersat_in, iflag_ratqs, fl_cor_ebil_in, &
+SUBROUTINE lscp_ini(dtime,lunout_in,prt_level_in,ok_ice_supersat_in, iflag_ratqs_in, fl_cor_ebil_in, &
                     RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, RVTMP2_in, &
                     RTT_in, RD_in, RV_in, RG_in, RPI_in, EPS_W_in)
@@ -309 +352 @@
 
    USE ioipsl_getin_p_mod, ONLY : getin_p
-   USE lmdz_cloudth_ini, ONLY : cloudth_ini
 
    REAL, INTENT(IN)      :: dtime
-   INTEGER, INTENT(IN)   :: lunout_in,prt_level_in,iflag_ratqs,fl_cor_ebil_in
+   INTEGER, INTENT(IN)   :: lunout_in,prt_level_in,iflag_ratqs_in,fl_cor_ebil_in
    LOGICAL, INTENT(IN)   :: ok_ice_supersat_in
 
@@ -324 +366 @@
     prt_level=prt_level_in
     fl_cor_ebil=fl_cor_ebil_in
-
+    iflag_ratqs=iflag_ratqs_in
     ok_ice_supersat=ok_ice_supersat_in
 
@@ -351 +393 @@
     CALL getin_p('iflag_vice',iflag_vice)
     CALL getin_p('iflag_t_glace',iflag_t_glace)
-    CALL getin_p('iflag_cloudth_vert',iflag_cloudth_vert)
     CALL getin_p('iflag_gammasat',iflag_gammasat)
     CALL getin_p('iflag_rain_incloud_vol',iflag_rain_incloud_vol)
@@ -369 +410 @@
     CALL getin_p('ffallv_lsc',ffallv_lsc)
     CALL getin_p('ffallv_lsc',ffallv_con)
+    ! for poprecip and cloud phase
     CALL getin_p('coef_eva',coef_eva)
     coef_sub=coef_eva
@@ -383 +425 @@
     CALL getin_p('gamma_snwretro',gamma_snwretro)
     CALL getin_p('gamma_taud',gamma_taud)
+    CALL getin_p('gamma_mixth',gamma_mixth)
     CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
     CALL getin_p('temp_nowater',temp_nowater)
     CALL getin_p('ok_bug_phase_lscp',ok_bug_phase_lscp)
-    ! for poprecip
     CALL getin_p('ok_poprecip',ok_poprecip)
     CALL getin_p('ok_corr_vap_evasub',ok_corr_vap_evasub)
+    CALL getin_p('ok_growth_precip_deposition',ok_growth_precip_deposition)
     CALL getin_p('rain_int_min',rain_int_min) 
     CALL getin_p('gamma_agg',gamma_agg)
@@ -397 +440 @@
     CALL getin_p('tau_auto_snow_max',tau_auto_snow_max)
     CALL getin_p('tau_auto_snow_min',tau_auto_snow_min)
+    CALL getin_p('expo_tau_auto_snow', expo_tau_auto_snow)
+    CALL getin_p('alpha_freez',alpha_freez)
+    CALL getin_p('beta_freez',beta_freez)
     CALL getin_p('r_snow',r_snow)
     CALL getin_p('rain_fallspeed',rain_fallspeed)
@@ -427 +473 @@
     CALL getin_p('coef_shear_lscp',coef_shear_lscp)
     CALL getin_p('chi_mixing_lscp',chi_mixing_lscp)
-
-
+    ! for cloudth routines
+    CALL getin_p('ok_lscp_mergecond',ok_lscp_mergecond)
+    CALL getin_p('iflag_cloudth_vert',iflag_cloudth_vert)
+    CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin)
+    CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor)
+    CALL getin_p('cloudth_sigma1s_power',sigma1s_power)
+    CALL getin_p('cloudth_sigma2s_factor',sigma2s_factor)
+    CALL getin_p('cloudth_sigma2s_power',sigma2s_power)
+    CALL getin_p('cloudth_vert_alpha',vert_alpha)
+    vert_alpha_th=vert_alpha
+    CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
+    CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs)
 
     WRITE(lunout,*) 'lscp_ini, niter_lscp:', niter_lscp
@@ -439 +495 @@
     WRITE(lunout,*) 'lscp_ini, iflag_vice:', iflag_vice
     WRITE(lunout,*) 'lscp_ini, iflag_t_glace:', iflag_t_glace
-    WRITE(lunout,*) 'lscp_ini, iflag_cloudth_vert:', iflag_cloudth_vert
     WRITE(lunout,*) 'lscp_ini, iflag_gammasat:', iflag_gammasat
     WRITE(lunout,*) 'lscp_ini, iflag_rain_incloud_vol:', iflag_rain_incloud_vol
@@ -467 +522 @@
     WRITE(lunout,*) 'lscp_ini, naero5', naero5
     WRITE(lunout,*) 'lscp_ini, gamma_snwretro', gamma_snwretro
+    WRITE(lunout,*) 'lscp_ini, gamma_mixth', gamma_mixth
     WRITE(lunout,*) 'lscp_ini, gamma_taud', gamma_taud
     WRITE(lunout,*) 'lscp_ini, iflag_oldbug_fisrtilp', iflag_oldbug_fisrtilp
@@ -475 +531 @@
     WRITE(lunout,*) 'lscp_ini, ok_poprecip', ok_poprecip
     WRITE(lunout,*) 'lscp_ini, ok_corr_vap_evasub', ok_corr_vap_evasub
+    WRITE(lunout,*) 'lscp_ini, ok_growth_precip_deposition', ok_growth_precip_deposition
     WRITE(lunout,*) 'lscp_ini, rain_int_min:', rain_int_min
     WRITE(lunout,*) 'lscp_ini, gamma_agg:', gamma_agg
@@ -483 +540 @@
     WRITE(lunout,*) 'lscp_ini, tau_auto_snow_max:',tau_auto_snow_max
     WRITE(lunout,*) 'lscp_ini, tau_auto_snow_min:',tau_auto_snow_min
+    WRITE(lunout,*) 'lscp_ini, expo_tau_auto_snow:',expo_tau_auto_snow
     WRITE(lunout,*) 'lscp_ini, r_snow:', r_snow
+    WRITE(lunout,*) 'lscp_ini, alpha_freez:', alpha_freez
+    WRITE(lunout,*) 'lscp_ini, beta_freez:', beta_freez
     WRITE(lunout,*) 'lscp_ini, rain_fallspeed_clr:', rain_fallspeed_clr
     WRITE(lunout,*) 'lscp_ini, rain_fallspeed_cld:', rain_fallspeed_cld
@@ -508 +568 @@
     WRITE(lunout,*) 'lscp_ini, coef_shear_lscp:', coef_shear_lscp
     WRITE(lunout,*) 'lscp_ini, chi_mixing_lscp:', chi_mixing_lscp
-
-
-
-
+    ! for cloudth routines
+    WRITE(lunout,*) 'lscp_ini, ok_lscp_mergecond:', ok_lscp_mergecond
+    WRITE(lunout,*) 'lscp_ini, iflag_cloudth_vert:', iflag_cloudth_vert
+    WRITE(lunout,*) 'lscp_ini, cloudth_ratqsmin:', cloudth_ratqsmin
+    WRITE(lunout,*) 'lscp_ini, cloudth_sigma1s_factor:', sigma1s_factor
+    WRITE(lunout,*) 'lscp_ini, cloudth_sigma1s_power:', sigma1s_power
+    WRITE(lunout,*) 'lscp_ini, cloudth_sigma2s_factor:', sigma2s_factor
+    WRITE(lunout,*) 'lscp_ini, cloudth_sigma2s_power:', sigma2s_power
+    WRITE(lunout,*) 'lscp_ini, cloudth_vert_alpha:', vert_alpha
+    WRITE(lunout,*) 'lscp_ini, cloudth_vert_alpha_th:', vert_alpha_th
+    WRITE(lunout,*) 'lscp_ini, iflag_cloudth_vert_noratqs:', iflag_cloudth_vert_noratqs
+
+
+    ! check consistency for cloud phase partitioning options
+
+    IF ((iflag_icefrac .GE. 2) .AND. (.NOT. ok_lscp_mergecond)) THEN
+      abort_message = 'in lscp, iflag_icefrac .GE. 2 works only if ok_lscp_mergecond=.TRUE.'
+      CALL abort_physic (modname,abort_message,1)
+    ENDIF
 
     ! check for precipitation sub-time steps
@@ -522 +597 @@
     ! and other options
    
-    IF (iflag_autoconversion .EQ. 2) THEN
+    IF ((iflag_autoconversion .EQ. 2) .AND. .NOT. ok_poprecip) THEN
         IF ((iflag_vice .NE. 0) .OR. (niter_lscp .GT. 1)) THEN
            abort_message = 'in lscp, iflag_autoconversion=2 requires iflag_vice=0 and niter_lscp=1'
@@ -539 +614 @@
       CALL abort_physic (modname,abort_message,1)
     ENDIF
+
+    IF ( (iflag_icefrac .GE. 1) .AND. (.NOT. ok_poprecip .AND. (iflag_evap_prec .LT. 4)) ) THEN
+      abort_message = 'in lscp, icefracturb works with poprecip or with precip evap option >=4'
+      CALL abort_physic (modname,abort_message,1)
+    ENDIF
+
 
 
@@ -547 +628 @@
     a_tr_sca(4) = -0.5
     
-    CALL cloudth_ini(iflag_cloudth_vert,iflag_ratqs)
 
 RETURN

LMDZ6/trunk/libf/phylmd/lmdz_lscp_old.f90

@@ -70 +70 @@
   USE yomcst_mod_h
   USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
-  USE lmdz_cloudth, only : cloudth, cloudth_v3, cloudth_v6
+  USE lmdz_lscp_condensation, only : cloudth, cloudth_v3, cloudth_v6
 
   USE lmdz_lscp_ini, ONLY: prt_level, lunout

LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90

@@ -330 +330 @@
 LOGICAL, INTENT(IN)                     :: iftop          !--if top of the column
 
+
 REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
 REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
@@ -540 +541 @@
           IF (zoliqi(i) .GT. 0.) THEN
             zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) & 
-              +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo))
+              +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)/zneb(i)*zrho(i)*ffallv)**(-1./dice_velo))
           ELSE
             zfroi=0.
@@ -618 +619 @@
     ! Computation of DT if all the liquid precip freezes
     DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
+    
+
     ! T should not exceed the freezing point
     ! that is Delta > RTT-zt(i)
@@ -722 +725 @@
 USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
 USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub, ok_ice_supersat, ok_unadjusted_clouds
-USE lmdz_lscp_ini, ONLY : eps, temp_nowater
+USE lmdz_lscp_ini, ONLY : eps, temp_nowater, ok_growth_precip_deposition
 USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
 
@@ -792 +795 @@
 dqreva(:) = 0.
 dqssub(:) = 0.
-dqrevap   = 0.
-dqssubl   = 0.
 
 !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
@@ -880 +881 @@
 DO i = 1, klon
 
+  dqrevap   = 0.
+  dqssubl   = 0.
   !--If there is precipitation from the layer above
   IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN
@@ -913 +916 @@
       IF ( ( 1. - precipfraccld(i) ) .GT. eps ) THEN
         qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) )
+      ELSE
+        qvapclr = 0.
+      ENDIF
+      IF (  precipfraccld(i)  .GT. eps ) THEN
+        qvapcld = MAX(qtotupnew(i)-qvapclrup(i) , 0.) / qtotupnew(i) * qvap(i) /  precipfraccld(i) 
+      ELSE
+        qvapcld = 0.
       ENDIF
     ELSE
@@ -918 +928 @@
       !--for the evap / subl process
       qvapclr = qvap(i)
+      qvapcld = qvap(i)
     ENDIF
 
@@ -940 +951 @@
       !--NB. with ok_ice_supersat activated, this barrier should be useless
       drainclreva = MIN(0., drainclreva)
-
+     
+      ! we set it to 0 as not sufficiently tested
+      drainclreva = 0.
 
       !--Sublimation of the solid precipitation coming from above
@@ -986 +999 @@
       ENDIF
 
+    ELSE
+            
+    !--All the precipitation is sublimated if the fraction is zero
+       drainclreva = - rainclr_tmp(i)
+       dsnowclrsub = - snowclr_tmp(i)
+
     ENDIF ! precipfracclr_tmp .GT. eps
 
@@ -993 +1012 @@
     !---------------------------
 
-    IF ( ok_unadjusted_clouds .AND. ( temp(i) .LE. temp_nowater ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN
+    IF ( ok_growth_precip_deposition .AND. ( temp(i) .LE. RTT ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN
       !--Evaporation of liquid precipitation coming from above
       !--in the cloud only
@@ -1000 +1019 @@
       !--Exact explicit formulation (raincld is resolved exactly, qvap explicitly)
       !--which does not need a barrier on raincld, because included in the formula
-      !draincldeva = precipfraccld_tmp(i) * MAX(0., &
-      !            - coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) &
-      !            + ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) &
-      !            )**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i)
-               
+      
+      draincldeva = precipfraccld_tmp(i) * MAX(0., &
+                  - coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) &
+                  + ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) &
+                  )**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i)
+
       !--Evaporation is limited by 0
       !--NB. with ok_ice_supersat activated, this barrier should be useless
-      !draincldeva = MIN(0., draincldeva)
-      draincldeva = 0.
+      draincldeva = MIN(0., draincldeva)
 
 
@@ -1277 +1296 @@
 
 USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb,    &
-                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, & 
+                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc,               & 
                           thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
                           rho_rain, r_rain, r_snow, rho_ice,                   &
+                          expo_tau_auto_snow,                                  &
                           tau_auto_snow_min, tau_auto_snow_max,                &
-                          thresh_precip_frac, eps,                             &
+                          thresh_precip_frac, eps, rain_int_min,               &
                           gamma_melt, alpha_freez, beta_freez, temp_nowater,   &
                           iflag_cloudth_vert, iflag_rain_incloud_vol,          &
@@ -1340 +1360 @@
 REAL, DIMENSION(klon) :: dhum_to_dflux
 REAL, DIMENSION(klon) :: qtot                             !--includes vap, liq, ice and precip
+REAL                  :: min_precip                       !--minimum precip flux below which precip fraction decreases
 
 !--Collection, aggregation and riming
@@ -1502 +1523 @@
         !--tau for snow depends on the ice fraction in mixed-phase clouds     
         tau_auto_snow = tau_auto_snow_max &
-                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
+                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ((1. - icefrac(i))**expo_tau_auto_snow)
 
         expo_auto_rain = cld_expo_con
@@ -1513 +1534 @@
         !--tau for snow depends on the ice fraction in mixed-phase clouds     
         tau_auto_snow = tau_auto_snow_max &
-                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
+                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ((1. - icefrac(i))**expo_tau_auto_snow)
 
         expo_auto_rain = cld_expo_lsc
@@ -1534 +1555 @@
                - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
 
-
     !--Barriers so that we don't create more rain/snow
     !--than there is liquid/ice
@@ -1543 +1563 @@
     qliq(i) = qliq(i) + dqlauto
     qice(i) = qice(i) + dqiauto
+
     raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i)
     snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i)
@@ -1710 +1731 @@
   !--second: immersion freezing following (inspired by Bigg 1953)
   !--the latter is parameterized as an exponential decrease of the rain
-  !--water content with a homemade formulya
+  !--water content with a homemade formula
   !--This is based on a caracteritic time of freezing, which
   !--exponentially depends on temperature so that it is
@@ -1717 +1738 @@
   !--NB.: this process needs a temperature adjustment
   !--dqrfreez_max : maximum rain freezing so that temperature
-  !--              stays lower than 273 K [kg/kg]
+  !--               stays lower than 273 K [kg/kg]
   !--tau_freez    : caracteristic time of freezing [s]
   !--gamma_freez  : tuning parameter [s-1]
@@ -1798 +1819 @@
               * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
 
-
     !--In clear air
     IF ( rainclr(i) .GT. 0. ) THEN
@@ -1827 +1847 @@
     !--Add tendencies
     !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
+    
     rainclr(i) = MAX(0., rainclr(i) + dqrclrfreez * dhum_to_dflux(i))
     raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
@@ -1833 +1854 @@
 
 
+
     !--Temperature adjustment with the uptake of latent
     !--heat because of freezing
+
     temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD &
                       / ( 1. + RVTMP2 * qtot(i) )
-
     !--Diagnostic tendencies
     dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2          
@@ -1847 +1869 @@
 
 
-  !--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
-  !--we reduce the precipiration fraction in the clear/cloudy air so that the new
-  !--local flux of rain+snow is equal to rain_int_min.
+  !--If the local flux of rain+snow in clear air is lower than min_precip,
+  !--we reduce the precipiration fraction in the clear air so that the new
+  !--local flux of rain+snow is equal to min_precip.
+  !--we apply the minimum only on the clear-sky fraction because the cloudy precip fraction
+  !--already decreases out of clouds
   !--Here, rain+snow is the gridbox-mean flux of precip.
   !--Therefore, (rain+snow)/precipfrac is the local flux of precip.
-  !--If the local flux of precip is lower than rain_int_min, i.e.,
-  !-- (rain+snow)/precipfrac < rain_int_min , i.e., 
-  !-- (rain+snow)/rain_int_min < precipfrac , then we want to reduce
-  !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/rain_int_min.
+  !--If the local flux of precip is lower than min_precip, i.e.,
+  !-- (rain+snow)/precipfrac < min_precip , i.e., 
+  !-- (rain+snow)/min_precip < precipfrac , then we want to reduce
+  !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/min_precip.
   !--Note that this is physically different than what is proposed in LTP thesis.
-  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
-  precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
+  !--min_precip is either equal to rain_int_min or calculated as a very small fraction
+  !--of the minimum precip flux estimated as the flux associated with the 
+  !--autoconversion threshold mass content
+  !min_precip=1.e-6*(pplay(i)/RD/temp(i))*MIN(rain_fallspeed_clr*cld_lc_lsc,snow_fallspeed_clr*cld_lc_lsc_snow)
+  min_precip=rain_int_min
+  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / min_precip )
 
   !--Calculate outputs

LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.f90

@@ -225 +225 @@
         ENDIF
 
-        ! if temperature of cloud top <-40°C, 
+        ! if temperature or temperature of cloud top <-40°C, 
         IF (iflag_t_glace .GE. 4) THEN
                 IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN 
@@ -241 +241 @@
 
 
-SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, temp, pplay, paprsdn, paprsup, omega, qice_ini, snowcld, qtot_incl, cldfra, tke,   &
-                             tke_dissip, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb)
+SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, pticefracturb, temp, pplay, paprsdn, paprsup, wvel, qice_ini, snowcld, qtot_incl, cldfra, tke,   &
+                             tke_dissip, sursat_e, invtau_e, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb)
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   ! Compute the liquid, ice and vapour content (+ice fraction) based 
   ! on turbulence (see Fields 2014, Furtado 2016, Raillard 2025)
   ! L.Raillard (23/09/24)
+  ! E.Vignon (03/2025) : additional elements for treatment of convective 
+  !                      boundary layer clouds
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
@@ -253 +255 @@
    USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI
    USE lmdz_lscp_ini, ONLY : seuil_neb, temp_nowater
-   USE lmdz_lscp_ini, ONLY : naero5, gamma_snwretro, gamma_taud, capa_crystal
+   USE lmdz_lscp_ini, ONLY : naero5, gamma_snwretro, gamma_taud, capa_crystal, rho_ice
    USE lmdz_lscp_ini, ONLY : eps
 
@@ -260 +262 @@
    INTEGER,   INTENT(IN)                           :: klon                !--number of horizontal grid points
    REAL,      INTENT(IN)                           :: dtime               !--time step [s]
-
+   LOGICAL,   INTENT(IN),       DIMENSION(klon)    :: pticefracturb       !--grid points concerned by this routine  
    REAL,      INTENT(IN),       DIMENSION(klon)    :: temp                !--temperature
    REAL,      INTENT(IN),       DIMENSION(klon)    :: pplay               !--pressure in the middle of the layer           [Pa]
    REAL,      INTENT(IN),       DIMENSION(klon)    :: paprsdn             !--pressure at the bottom interface of the layer [Pa]
    REAL,      INTENT(IN),       DIMENSION(klon)    :: paprsup             !--pressure at the top interface of the layer    [Pa]
-   REAL,      INTENT(IN),       DIMENSION(klon)    :: omega               !--resolved vertical velocity                    [Pa/s]
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: wvel                !--vertical velocity                             [m/s]
    REAL,      INTENT(IN),       DIMENSION(klon)    :: qtot_incl           !--specific total cloud water in-cloud content   [kg/kg]
    REAL,      INTENT(IN),       DIMENSION(klon)    :: cldfra              !--cloud fraction in gridbox                     [-]
@@ -272 +274 @@
 
    REAL,      INTENT(IN),       DIMENSION(klon)    :: qice_ini            !--initial specific ice content gridbox-mean     [kg/kg]
-   REAL,      INTENT(IN),       DIMENSION(klon)    :: snowcld
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: snowcld             !--in-cloud snowfall flux                        [kg/m2/s]
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: sursat_e            !--environment supersaturation                   [-]
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: invtau_e            !--inverse time-scale of mixing with environment [s-1]
    REAL,      INTENT(OUT),      DIMENSION(klon)    :: qliq                !--specific liquid content gridbox-mean          [kg/kg]
    REAL,      INTENT(OUT),      DIMENSION(klon)    :: qvap_cld            !--specific cloud vapor content, gridbox-mean    [kg/kg]
    REAL,      INTENT(OUT),      DIMENSION(klon)    :: qice                !--specific ice content gridbox-mean             [kg/kg]
-   REAL,      INTENT(OUT),      DIMENSION(klon)    :: icefrac             !--fraction of ice in condensed water            [-]
-   REAL,      INTENT(OUT),      DIMENSION(klon)    :: dicefracdT
+
+   REAL,      INTENT(INOUT),    DIMENSION(klon)    :: icefrac             !--fraction of ice in condensed water            [-]
+   REAL,      INTENT(INOUT),    DIMENSION(klon)    :: dicefracdT
 
    REAL,      INTENT(OUT),      DIMENSION(klon)    :: cldfraliq           !--fraction of cldfra where liquid               [-]
@@ -291 +296 @@
    REAL :: C0                                                             !--Lagrangian structure function                 [-]
    REAL :: tau_dissipturb
-   REAL :: tau_phaserelax
-   REAL :: sigma2_pdf, mean_pdf
+   REAL :: invtau_phaserelax
+   REAL :: sigma2_pdf
    REAL :: ai, bi, B0
    REAL :: sursat_iceliq
@@ -302 +307 @@
    REAL :: N0_PSD, lambda_PSD                                             !--parameters of the exponential PSD
 
-   REAL :: rho_ice                                                        !--ice density                                   [kg/m3]
    REAL :: cldfra1D
    REAL :: rho_air
    REAL :: psati                                                          !--saturation vapor pressure wrt ice             [Pa]
-   
-   REAL :: vitw                                                           !--vertical velocity                             [m/s]
+
                                                                        
+    REAL :: tempvig1, tempvig2
+
+   tempvig1    = -21.06 + RTT
+   tempvig2    = -30.35 + RTT
    C0            = 10.                                                    !--value assumed in Field2014            
-   rho_ice       = 950.
    sursat_iceext = -0.1
    qzero(:)      = 0.
    cldfraliq(:)  = 0.
-   icefrac(:)    = 0.
-   dicefracdT(:) = 0.
-
+   qliq(:)       = 0.
+   qice(:)       = 0.
+   qvap_cld(:)   = 0.
    sigma2_icefracturb(:) = 0.
    mean_icefracturb(:)   = 0.
@@ -327 +333 @@
 
     DO i=1,klon
-
      rho_air  = pplay(i) / temp(i) / RD
-
      ! because cldfra is intent in, but can be locally modified due to test 
      cldfra1D = cldfra(i)
-     IF (cldfra(i) .LE. 0.) THEN
-        qvap_cld(i)   = 0.
-        qliq(i)       = 0.
-        qice(i)       = 0.
-        cldfraliq(i)  = 0.
-        icefrac(i)    = 0.
-        dicefracdT(i) = 0.
-
-     ! If there is a cloud
-     ELSE
+     ! activate param for concerned grid points and for cloudy conditions
+     IF ((pticefracturb(i)) .AND. (cldfra(i) .GT. 0.)) THEN
         IF (cldfra(i) .GE. 1.0) THEN
            cldfra1D = 1.0
@@ -364 +360 @@
            dicefracdT(i) = 0.
 
+
         !---------------------------------------------------------
         !--             MIXED PHASE TEMPERATURE REGIME          
@@ -374 +371 @@
         ELSE
 
-           vitw = -omega(i) / RG / rho_air
-           qiceini_incl  = qice_ini(i) / cldfra1D + snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) ) / cldfra1D
-
-           !--1. No preexisting ice : if vertical motion, fully liquid
+           qiceini_incl  = qice_ini(i) / cldfra1D + gamma_snwretro * snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) ) / cldfra1D
+           ! assuming snowfall velocity = 1 m/s:
+           !qiceini_incl  = qice_ini(i) / cldfra1D + gamma_snwretro * snowcld(i) / pplay(i) * RD * temp(i) / 1. / cldfra1D
+
+           !--1. No preexisting ice and no mixing with environment: if vertical motion, fully liquid
            !--cloud else fully iced cloud
-           IF ( qiceini_incl .LT. eps ) THEN
-              IF ( (vitw .GT. eps) .OR. (tke(i) .GT. eps) ) THEN
+           IF ( (qiceini_incl .LT. eps) .AND. (invtau_e(i) .LT. eps) ) THEN
+              IF ( (wvel(i) .GT. eps) .OR. (tke(i) .GT. eps) ) THEN
                  qvap_cld(i)   = qsatl(i) * cldfra1D
                  qliq(i)       = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D
@@ -397 +395 @@
            
 
-           !--2. Pre-existing ice :computation of ice properties for
+           !--2. Pre-existing ice and/or mixing with environment:computation of ice properties for
            !--feedback
            ELSE
-              ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. )
-
-              sursat_equ    = ai * vitw * tau_phaserelax
 
               sursat_iceliq = qsatl(i)/qsati(i) - 1.
@@ -410 +405 @@
               !--are computed following Morrison&Gettelman 2008
               !--Ice number density is assumed equals to INP density 
-              !--which is a function of temperature (DeMott 2010)  
+              !--which is for naero5>0 a function of temperature (DeMott 2010)    
               !--bi and B0 are microphysical function characterizing 
               !--vapor/ice interactions
@@ -416 +411 @@
               !--onto ice crystals
 
-              nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )**3.33 * naero5**(0.0264*(RTT-temp(i))+0.0033)
-              lambda_PSD  = ( (RPI*rho_ice*nb_crystals) / (rho_air * qiceini_incl ) ) ** (1./3.)
+              !--For naero5<=0 INP density is derived from the empirical fit 
+              !--from MARCUS campaign from Vignon 2021
+              !--/!\ Note that option is very specific and should be use for
+              !--the Southern Ocean and the Antarctic
+
+              IF (naero5 .LE. 0) THEN
+                 IF ( temp(i) .GT. tempvig1 ) THEN
+                      nb_crystals = 1.e3 * 10**(-0.14*(temp(i)-tempvig1) - 2.88)
+                 ELSE IF ( temp(i) .GT. tempvig2 ) THEN 
+                      nb_crystals = 1.e3 * 10**(-0.31*(temp(i)-tempvig1) - 2.88)
+                 ELSE
+                      nb_crystals = 1.e3 * 10**(0.)
+                 ENDIF
+              ELSE
+                 nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )**3.33 * naero5**(0.0264*(RTT-temp(i))+0.0033)
+              ENDIF
+              lambda_PSD  = ( (RPI*rho_ice*nb_crystals) / (rho_air * MAX(qiceini_incl , eps) ) ) ** (1./3.)
               N0_PSD      = nb_crystals * lambda_PSD
               moment1_PSD = N0_PSD/lambda_PSD**2
@@ -431 +441 @@
               B0 = 4. * RPI * capa_crystal * 1. / (  RLSTT**2 / air_thermal_conduct / RV / temp(i)**2  &
                                                   +  RV * temp(i) / psati / water_vapor_diff  )
-              tau_phaserelax = 1. / (bi * B0 * moment1_PSD )
+              invtau_phaserelax = bi * B0 * moment1_PSD 
               
               ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. )
-
-              !--2A. No TKE : stationnary binary solution depending on omega
+              sursat_equ    = (ai * wvel(i) + sursat_e(i)*invtau_e(i)) / (invtau_phaserelax + invtau_e(i))
+              sursat_equ    = MAX(-1., MIN(sursat_equ , 1.)) 
+              ! as sursaturation is by definition lower than -1 and 
+              ! because local supersaturation > 1 are never found in the atmosphere
+
+              !--2A. No TKE : stationnary binary solution depending on vertical velocity and mixing with env.
               ! If Sequ > Siw liquid cloud, else ice cloud
               IF ( tke_dissip(i) .LE. eps )  THEN
+                 sigma2_icefracturb(i)= 0.
+                 mean_icefracturb(i)  = sursat_equ
                  IF (sursat_equ .GT. sursat_iceliq) THEN
                     qvap_cld(i)   = qsatl(i) * cldfra1D
@@ -474 +490 @@
 
                  liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - (qice_ini(i) / cldfra1D) - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) )
-                 sigma2_pdf = 1./2. * ( ai**2 ) *  2./3. * tke(i) * tau_dissipturb * tau_phaserelax
-                 
-                 mean_pdf = ai * vitw * tau_phaserelax 
-                 
-                 cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - mean_pdf) / (SQRT(2.* sigma2_pdf) ) ) )
+                 sigma2_pdf = MIN(1., 1./2. * ( ai**2 ) *  2./3. * tke(i) * tau_dissipturb / (invtau_phaserelax + invtau_e(i))) 
+                 ! sursat ranges between -1 and 1, so we prevent sigma2 so exceed 1
+                 cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - sursat_equ) / (SQRT(2.* sigma2_pdf) ) ) )
                  IF (cldfraliq(i) .GT. liqfra_max) THEN
                      cldfraliq(i) = liqfra_max
                  ENDIF 
                  
-                 qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.*RPI) * EXP( -1.*(sursat_iceliq - mean_pdf)**2. / (2.*sigma2_pdf) )  &
-                           - qsati(i) * cldfraliq(i) * (sursat_iceliq - mean_pdf )
+                 qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.*RPI) * EXP( -1.*(sursat_iceliq - sursat_equ)**2. / (2.*sigma2_pdf) )  &
+                           - qsati(i) * cldfraliq(i) * (sursat_iceliq - sursat_equ )
                  
                  sigma2_icefracturb(i)= sigma2_pdf
-                 mean_icefracturb(i)  = mean_pdf 
+                 mean_icefracturb(i)  = sursat_equ 
      
                  !------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION  ------------
@@ -505 +519 @@
                  IF ( qvap_incl  .GE. qtot_incl(i) ) THEN
                     qvap_incl = qsati(i)
-                    qliq_incl = qtot_incl(i) - qvap_incl
+                    qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)
                     qice_incl = 0.
 
@@ -518 +532 @@
                  qliq(i)       = qliq_incl * cldfra1D
                  qice(i)       = qice_incl * cldfra1D
-                 icefrac(i)    = qice(i) / ( qice(i) + qliq(i) )
+                 IF ((qice(i)+qliq(i)) .GT. 0.) THEN
+                    icefrac(i)    = qice(i) / ( qice(i) + qliq(i) )
+                 ELSE
+                    icefrac(i)    = 1. ! to keep computation of qsat wrt ice in condensation loop in lmdz_lscp_main
+                 ENDIF
                  dicefracdT(i) = 0.
 
@@ -527 +545 @@
         END IF ! ! MPC temperature
 
-     END IF ! cldfra
+     END IF ! pticefracturb and cldfra
    
    ENDDO ! klon

LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90

@@ -623 +623 @@
       REAL, SAVE, ALLOCATABLE :: sigma2_icefracturb(:,:)
 !$OMP THREADPRIVATE(sigma2_icefracturb)
+      REAL, SAVE, ALLOCATABLE :: cldfraliqth(:,:)
+!$OMP THREADPRIVATE(cldfraliqth)
+      REAL, SAVE, ALLOCATABLE ::mean_icefracturbth(:,:)
+!$OMP THREADPRIVATE(mean_icefracturbth)
+      REAL, SAVE, ALLOCATABLE :: sigma2_icefracturbth(:,:)
+!$OMP THREADPRIVATE(sigma2_icefracturbth)
 
 ! variables de sorties MM
@@ -1201 +1207 @@
       ALLOCATE(sigma2_icefracturb(klon,klev))
       ALLOCATE(mean_icefracturb(klon,klev))
+      ALLOCATE(cldfraliqth(klon,klev))
+      ALLOCATE(sigma2_icefracturbth(klon,klev))
+      ALLOCATE(mean_icefracturbth(klon,klev))
       ALLOCATE(distcltop(klon,klev))
       ALLOCATE(temp_cltop(klon,klev))
@@ -1619 +1628 @@
       DEALLOCATE(sigma2_icefracturb)
       DEALLOCATE(mean_icefracturb)
+      DEALLOCATE(cldfraliqth)
+      DEALLOCATE(sigma2_icefracturbth)
+      DEALLOCATE(mean_icefracturbth)
       DEALLOCATE (zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic)
       DEALLOCATE(distcltop)

LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90

@@ -1592 +1592 @@
   TYPE(ctrl_out), SAVE :: o_rneb = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'rneb', 'Cloud fraction', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_distcltop = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'distcltop', 'Distance from cloud top', 'm', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_tempcltop = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'tempcltop', 'Cloud top temperature', 'K', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_cldfraliq = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
-    'cldfraliq', 'Liquid fraction of the cloud', '-', (/ ('', i=1, 10) /))
+    'cldfraliq', 'Liquid fraction of the cloud part of the mesh', '-', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_sigma2_icefracturb = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     'sigma2_icefracturb', 'Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]', '-', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_mean_icefracturb = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     'mean_icefracturb', 'Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]', '-', (/ ('', i=1, 10) /))
- 
+   TYPE(ctrl_out), SAVE :: o_cldfraliqth = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'cldfraliqth', 'Liquid fraction of clouds in thermals', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_sigma2_icefracturbth = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'sigma2_icefracturbth', 'Variance of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-]', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_mean_icefracturbth = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'mean_icefracturbth', 'Mean of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-]', '-', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_rnebjn = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11,11, 11/), &      
     'rnebjn', 'Cloud fraction in day', '-', (/ ('', i=1, 10) /))

LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90

@@ -145 +145 @@
          o_zfull, o_zhalf, o_rneb, o_rnebjn, o_rnebcon, &
          o_rnebls, o_rneblsvol, o_rhum, o_rhl, o_rhi, o_ozone, o_ozone_light, &
+         o_distcltop, o_tempcltop,  &
          o_pfraclr, o_pfracld, o_cldfraliq, o_sigma2_icefracturb, o_mean_icefracturb,  &
+         o_cldfraliqth, o_sigma2_icefracturbth, o_mean_icefracturbth,  & 
          o_qrainlsc, o_qsnowlsc, o_dqreva, o_dqrauto, o_dqrcol, o_dqrmelt, o_dqrfreez, &
          o_dqssub, o_dqsauto, o_dqsagg, o_dqsrim, o_dqsmelt, o_dqsfreez, &
@@ -269 +271 @@
          o_SAD_sulfate, o_reff_sulfate, o_sulfmmr, o_nd_mode, o_sulfmmr_mode,o_SO2_chlm
 
-    USE lmdz_lscp_ini, ONLY: ok_poprecip
+    USE lmdz_lscp_ini, ONLY: ok_poprecip, iflag_icefrac
 
     USE phys_output_ctrlout_mod, ONLY: o_heat_volc, o_cool_volc !NL
@@ -377 +379 @@
          zphi, u_seri, v_seri, omega, cldfra, &
          rneb, rnebjn, rneblsvol,  &
-         zx_rh, zx_rhl, zx_rhi, &
+         zx_rh, zx_rhl, zx_rhi, distcltop, temp_cltop, &
          pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb, &
+         cldfraliqth, sigma2_icefracturbth, mean_icefracturbth, &
          qraindiag, qsnowdiag, dqreva, dqssub, &
          dqrauto,dqrcol,dqrmelt,dqrfreez, &
@@ -2087 +2090 @@
            CALL histwrite_phy(o_pfraclr, pfraclr)
            CALL histwrite_phy(o_pfracld, pfracld)
+           IF (iflag_icefrac .GT. 0) THEN
            CALL histwrite_phy(o_cldfraliq, cldfraliq)
            CALL histwrite_phy(o_sigma2_icefracturb, sigma2_icefracturb)
            CALL histwrite_phy(o_mean_icefracturb, mean_icefracturb)
+           CALL histwrite_phy(o_cldfraliqth, cldfraliqth)
+           CALL histwrite_phy(o_sigma2_icefracturbth, sigma2_icefracturbth)
+           CALL histwrite_phy(o_mean_icefracturbth, mean_icefracturbth)
+           ELSE
+           CALL histwrite_phy(o_distcltop, distcltop)
+           CALL histwrite_phy(o_tempcltop, temp_cltop)
+           ENDIF
            IF (ok_poprecip) THEN
            CALL histwrite_phy(o_qrainlsc, qraindiag)

LMDZ6/trunk/libf/phylmd/physiq_mod.F90

@@ -75 +75 @@
     USE write_field_phy
     use wxios_mod, ONLY: g_ctx, wxios_set_context
-    USE lmdz_lscp, ONLY : lscp
+    USE lmdz_lscp_main, ONLY : lscp
     USE lmdz_call_cloud_optics_prop, ONLY : call_cloud_optics_prop
     USE lmdz_lscp_old, ONLY : fisrtilp, fisrtilp_first
@@ -323 +323 @@
        !
        rneblsvol, &
-       pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb,  &
+       pfraclr, pfracld, & 
+       cldfraliq, sigma2_icefracturb, mean_icefracturb,  &
+       cldfraliqth, sigma2_icefracturbth, mean_icefracturbth,  &
        distcltop, temp_cltop,  &
        !-- LSCP - condensation and ice supersaturation variables
@@ -3885 +3887 @@
          t_seri, q_seri, ql_seri_lscp, qi_seri_lscp, ptconv, ratqs, sigma_qtherm, &
          d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, &
-         pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb,  &
+         pfraclr, pfracld, cldfraliq, cldfraliqth,              &
+         sigma2_icefracturb, sigma2_icefracturbth,              &
+         mean_icefracturb,  mean_icefracturbth,                 &
          radocond, picefra, rain_lsc, snow_lsc, &
          frac_impa, frac_nucl, beta_prec_fisrt, &
          prfl, psfl, rhcl,  &
-         zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
+         zqasc, fraca,ztv,zpspsk,ztla,zthl,zw2,iflag_cld_th, &
          iflag_ice_thermo, distcltop, temp_cltop,   &
          pbl_tke(:,:,is_ave), pbl_eps(:,:,is_ave), &
+         entr_therm, detr_therm, &
          cell_area, &
          cf_seri, rvc_seri, u_seri, v_seri, &

LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_old.F90

@@ -22 +22 @@
   USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
   USE print_control_mod, ONLY: prt_level, lunout
-  USE lmdz_cloudth, only : cloudth, cloudth_v3, cloudth_v6
+  USE lmdz_lscp_condensation, only : cloudth, cloudth_v3, cloudth_v6
   USE ioipsl_getin_p_mod, ONLY : getin_p
   USE phys_local_var_mod, ONLY: ql_seri,qs_seri

LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

@@ -75 +75 @@
     USE write_field_phy
     use wxios_mod, ONLY: g_ctx, wxios_set_context
-    USE lmdz_lscp, ONLY : lscp
+    USE lmdz_lscp_main, ONLY : lscp
     USE lmdz_call_cloud_optics_prop, ONLY : call_cloud_optics_prop
     USE lmdz_lscp_old, ONLY : fisrtilp
@@ -365 +365 @@
        rneblsvol, &
        pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb,  &
+       cldfraliqth, sigma2_icefracturbth, mean_icefracturbth,  &
        distcltop, temp_cltop,  &
        !-- LSCP - condensation and ice supersaturation variables
@@ -5081 +5082 @@
          t_seri, q_seri, ql_seri_lscp, qi_seri_lscp, ptconv, ratqs, sigma_qtherm, &
          d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, &
-         pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb,  &
+         pfraclr, pfracld, cldfraliq, cldfraliqth, &
+         sigma2_icefracturb, sigma2_icefracturbth, &
+         mean_icefracturb, mean_icefracturbth,     &
          radocond, picefra, rain_lsc, snow_lsc, &
          frac_impa, frac_nucl, beta_prec_fisrt, &
          prfl, psfl, rhcl,  &
-         zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
+         zqasc, fraca,ztv,zpspsk,ztla,zthl,zw2,iflag_cld_th, &
          iflag_ice_thermo, distcltop, temp_cltop,   &
          pbl_tke(:,:,is_ave), pbl_eps(:,:,is_ave), &
+         entr_therm, detr_therm, &
          cell_area, &
          cf_seri, rvc_seri, u_seri, v_seri, &