Index: LMDZ6/trunk/libf/phylmd/cloudth_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/cloudth_mod.F90	(revision 3006)
+++ LMDZ6/trunk/libf/phylmd/cloudth_mod.F90	(revision 3493)
@@ -7,5 +7,5 @@
        SUBROUTINE cloudth(ngrid,klev,ind2,  &
      &           ztv,po,zqta,fraca, & 
-     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
+     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
      &           ratqs,zqs,t)
 
@@ -38,4 +38,5 @@
       REAL zpspsk(ngrid,klev)
       REAL paprs(ngrid,klev+1)
+      REAL pplay(ngrid,klev)
       REAL ztla(ngrid,klev)
       REAL zthl(ngrid,klev)
@@ -78,5 +79,5 @@
       CALL cloudth_vert(ngrid,klev,ind2,  &
      &           ztv,po,zqta,fraca, & 
-     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
+     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
      &           ratqs,zqs,t)
       RETURN
@@ -251,5 +252,5 @@
      SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
      &           ztv,po,zqta,fraca, & 
-     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
+     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
      &           ratqs,zqs,t)
 
@@ -282,4 +283,5 @@
       REAL zpspsk(ngrid,klev)
       REAL paprs(ngrid,klev+1)
+      REAL pplay(ngrid,klev)
       REAL ztla(ngrid,klev)
       REAL zthl(ngrid,klev)
@@ -585,7 +587,10 @@
 END SUBROUTINE cloudth_vert
 
+
+
+
        SUBROUTINE cloudth_v3(ngrid,klev,ind2,  &
      &           ztv,po,zqta,fraca, & 
-     &           qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
+     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
      &           ratqs,zqs,t)
 
@@ -618,4 +623,5 @@
       REAL zpspsk(ngrid,klev)
       REAL paprs(ngrid,klev+1)
+      REAL pplay(ngrid,klev)
       REAL ztla(ngrid,klev)
       REAL zthl(ngrid,klev)
@@ -641,4 +647,5 @@
       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  
@@ -654,5 +661,5 @@
       CALL cloudth_vert_v3(ngrid,klev,ind2,  &
      &           ztv,po,zqta,fraca, & 
-     &           qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
+     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
      &           ratqs,zqs,t)
       RETURN
@@ -808,5 +815,5 @@
      SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2,  &
      &           ztv,po,zqta,fraca, & 
-     &           qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
+     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
      &           ratqs,zqs,t)
 
@@ -841,4 +848,5 @@
       REAL zpspsk(ngrid,klev)
       REAL paprs(ngrid,klev+1)
+      REAL pplay(ngrid,klev)
       REAL ztla(ngrid,klev)
       REAL zthl(ngrid,klev)
@@ -864,4 +872,5 @@
       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
@@ -888,4 +897,16 @@
       REAL zqs(ngrid), qcloud(ngrid)
       REAL erf
+
+      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
+
 
 !------------------------------------------------------------------------------
@@ -1050,5 +1071,5 @@
 
       ELSE IF (iflag_cloudth_vert >= 3) THEN
-
+      IF (iflag_cloudth_vert < 5) THEN
 !-------------------------------------------------------------------------------
 !  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
@@ -1119,17 +1140,68 @@
       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 (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
+      if (ctot(ind1,ind2).lt.1.e-10) then
       ctot(ind1,ind2)=0.
       ctot_vol(ind1,ind2)=0.
-      qcloud(ind1)=zqsatenv(ind1,ind2) 
-
-      else 
+      qcloud(ind1)=zqsatenv(ind1,ind2)
+
+      else
                 
       qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
@@ -1139,5 +1211,5 @@
       endif  
 
-      else  ! Environment only
+      else  ! gaussienne environnement seule
       
       zqenv(ind1)=po(ind1)
@@ -1151,18 +1223,17 @@
       
 
-!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
+!      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.
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 
+      
 
       sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
-      sigma2s=0.
-
-      xenv=senv/(sqrt2*sigma1s)
+
+      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*sqrt2*cenv(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
@@ -1172,8 +1243,12 @@
       else   
                 
+      ctot(ind1,ind2)=ctot(ind1,ind2) 
       qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
 
-      endif    
- 
+      endif  
+ 
+
+
+
       endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
       ! Outputs used to check the PDFs
@@ -1184,8 +1259,275 @@
 
       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)
+
+
+      USE ioipsl_getin_p_mod, ONLY : getin_p
+      USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
+     &                                cloudth_sigmath,cloudth_sigmaenv
+
+      IMPLICIT NONE
+
+#include "YOMCST.h"
+#include "YOETHF.h"
+#include "FCTTRE.h"
+#include "thermcell.h"
+#include "nuage.h"
+
+
+        !Domain variables
+      INTEGER ngrid !indice Max lat-lon
+      INTEGER klev  !indice Max alt
+      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 erf,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(:,:)=0.
+      qlenv(:,:)=0.  
+      qltot(:,:)=0.
+      cth_vol(:,:)=0.
+      cenv_vol(:,:)=0.
+      ctot_vol(:,:)=0.
+      cth_surf(:,:)=0.
+      cenv_surf(:,:)=0.
+      ctot_surf(:,:)=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 cloudth_mod
+
+
+
+
