Index: /LMDZ6/trunk/libf/phylmd/lmdz_cloudth.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_cloudth.F90	(revision 4909)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_cloudth.F90	(revision 4910)
@@ -1497,29 +1497,22 @@
 
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE cloudth_mpc(klon,klev,ind2,iflag_mpc_bl,mpc_bl_points,           &
-&           temp,ztv,po,zqta,fraca,zpspsk,paprs,pplay,ztla,zthl,            &
-&           ratqs,zqs,snowflux,qcloud,qincloud,icefrac,ctot,ctot_vol,       &
+SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points,                        &
+&           temp,qt,qt_th,frac_th,zpspsk,paprsup, paprsdn,play,thetal_th,   &
+&           ratqs,qcloud,qincloud,icefrac,ctot,ctot_vol,                    &
 &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! Author : Arnaud Octavio Jam (LMD/CNRS), Etienne Vignon (LMDZ/CNRS)
-! Date: Adapted from cloudth_vert_v3 in 2021
-! Aim : computes qc and rneb in thermals with cold microphysical considerations
-!       + for mixed phase boundary layer clouds, calculate ql and qi from 
-!       a stationary MPC model
-! IMPORTANT NOTE: we assume iflag_clouth_vert=3
+! Author : Etienne Vignon (LMDZ/CNRS)
+! Date: April 2024
+! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam
+! IMPORTANT NOTE: we assume iflag_cloudth_vert=7
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
       use lmdz_cloudth_ini, only: iflag_cloudth_vert,iflag_ratqs
-      use lmdz_cloudth_ini, only:  C_mpc ,Ni,C_cap,Ei,d_top ,vert_alpha, vert_alpha_th ,sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin ,iflag_cloudth_vert_noratqs
-      use lmdz_lscp_tools, ONLY: CALC_QSAT_ECMWF, FALLICE_VELOCITY
-
-      use phys_local_var_mod, ONLY : qlth, qith, qsith, wiceth
+      use lmdz_cloudth_ini, only: C_mpc ,Ni,C_cap,Ei,d_top ,vert_alpha, vert_alpha_th ,sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin ,iflag_cloudth_vert_noratqs
+      use lmdz_lscp_tools,  only: CALC_QSAT_ECMWF
+      use lmdz_lscp_ini,    only: RTT, RG, RPI, RD, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th
 
       IMPLICIT NONE
 
-      INCLUDE "YOMCST.h"
-      INCLUDE "YOETHF.h"
-      INCLUDE "FCTTRE.h"
-      
 
 !------------------------------------------------------------------------------
@@ -1532,28 +1525,27 @@
       
 
-      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  temp          ! Temperature [K]
-      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  ztv           ! Virtual potential temp [K]
-      REAL, DIMENSION(klon),      INTENT(IN)      ::  po            ! specific humidity [kg/kg]
-      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  zqta          ! specific humidity within thermals [kg/kg]
-      REAL, DIMENSION(klon,klev+1), INTENT(IN)    ::  fraca         ! Fraction of the mesh covered by thermals [0-1] 
-      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  zpspsk
-      REAL, DIMENSION(klon,klev+1), INTENT(IN)    ::  paprs         ! Pressure at layer interfaces [Pa]
-      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  pplay         ! Pressure at the center of layers [Pa]
-      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  ztla          ! Liquid temp [K]
-      REAL, DIMENSION(klon,klev), INTENT(INOUT)      ::  zthl       ! Liquid potential temp [K]
-      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  ratqs         ! Parameter that determines the width of the total water distrib.
-      REAL, DIMENSION(klon),      INTENT(IN)      ::  zqs           ! Saturation specific humidity in the mesh [kg/kg]
-      REAL, DIMENSION(klon,klev+1), INTENT(IN)    ::  snowflux      ! snow flux at the interface of the layer [kg/m2/s]
-      INTEGER,                      INTENT(IN)    ::  iflag_mpc_bl  ! option flag for mpc boundary layer clouds param.
+      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)      ::  paprsup       ! Pressure at top layer interface [Pa]
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  paprsdn       ! Pressure at bottom layer interface [Pa]
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  play          ! Pressure of layers [Pa]
+      REAL, DIMENSION(klon),      INTENT(IN)      ::  ratqs         ! Parameter that determines the width of the total water distrib.
 
 
       INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points  ! grid points with convective (thermals) mixed phase clouds
       
-      REAL, DIMENSION(klon,klev), INTENT(OUT)      ::  ctot         ! Cloud fraction [0-1]
-      REAL, DIMENSION(klon,klev), INTENT(OUT)      ::  ctot_vol     ! Volume cloud fraction [0-1]
+      REAL, DIMENSION(klon),      INTENT(OUT)      ::  ctot         ! Cloud fraction [0-1]
+      REAL, DIMENSION(klon),      INTENT(OUT)      ::  ctot_vol     ! Volume cloud fraction [0-1]
       REAL, DIMENSION(klon),      INTENT(OUT)      ::  qcloud       ! In cloud total water content [kg/kg] 
-      REAL, DIMENSION(klon),      INTENT(OUT)      ::  qincloud       ! In cloud condensed water content [kg/kg] 
-      REAL, DIMENSION(klon,klev), INTENT(OUT)      ::  icefrac      ! Fraction of ice in clouds [0-1]
-      real, dimension(klon,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
+      REAL, DIMENSION(klon),      INTENT(OUT)      ::  qincloud     ! In cloud condensed water content [kg/kg] 
+      REAL, DIMENSION(klon),      INTENT(OUT)      ::  icefrac      ! Fraction of ice in clouds [0-1]
+      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
 
 
@@ -1562,20 +1554,15 @@
       INTEGER itap,ind1,l,ig,iter,k 
       INTEGER iflag_topthermals, niter
-      LOGICAL falseklon(klon)
-
-      REAL zqsatth(klon), zqsatenv(klon), zqsatenvonly(klon)
-      REAL sigma1(klon,klev)                                                         
-      REAL sigma2(klon,klev)
-      REAL qcth(klon,klev)
-      REAL qcenv(klon,klev)
-      REAL qctot(klon,klev) 
-      REAL cth(klon,klev)
-      REAL cenv(klon,klev)   
-      REAL cth_vol(klon,klev)
-      REAL cenv_vol(klon,klev)
-      REAL rneb(klon,klev)
-      REAL zqenv(klon), zqth(klon), zqenvonly(klon) 
-      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
-      REAL rdd,cppd,Lv
+
+      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 icefracenv(klon), icefracth(klon)
+      REAL sqrtpi,sqrt2,sqrt2pi
       REAL alth,alenv,ath,aenv
       REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
@@ -1585,20 +1572,19 @@
       REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
       REAL zdelta,qsatbef,zcor
-      REAL Tbefth(klon), Tbefenv(klon), Tbefenvonly(klon), rhoth(klon)
+      REAL Tbefth(klon), Tbefenv(klon), zrho(klon), rhoth(klon)
       REAL qlbef
-      REAL dqsatenv(klon), dqsatth(klon), dqsatenvonly(klon)
+      REAL dqsatenv(klon), dqsatth(klon)
       REAL erf
       REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
       REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) 
-      REAL rhodz(klon,klev)
-      REAL zrho(klon,klev)
-      REAL dz(klon,klev)
-      REAL qslenv(klon), qslth(klon) 
+      REAL rhodz(klon)
+      REAL rho(klon)
+      REAL dz(klon)
+      REAL qslenv(klon), qslth(klon), qsienv(klon), qsith(klon)  
       REAL alenvl, aenvl
       REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
       REAL senvi, senvl, qbase, sbase, qliqth, qiceth
       REAL qmax, ttarget, stmp, cout, coutref, fraci
-      REAL maxi, mini, pas, temp_lim
-      REAL deltazlev_mpc(klev), temp_mpc(klev), pres_mpc(klev), fraca_mpc(klev+1), snowf_mpc(klev+1)
+      REAL maxi, mini, pas
 
       ! Modifty the saturation deficit PDF in thermals 
@@ -1615,53 +1601,19 @@
 ! Few initial checksS
 
-      IF (iflag_cloudth_vert.NE.3) THEN
-         abort_message = 'clouth_mpc cannot be used if iflag_cloudth_vert .NE. 3'
-         CALL abort_physic(routname,abort_message,1)
-      ENDIF
-
       DO k = 1,klev
       DO ind1 = 1, klon
-        rhodz(ind1,k) = (paprs(ind1,k)-paprs(ind1,k+1))/rg !kg/m2
-        zrho(ind1,k) = pplay(ind1,k)/temp(ind1,k)/rd !kg/m3
-        dz(ind1,k) = rhodz(ind1,k)/zrho(ind1,k) !m : epaisseur de la couche en metre
+        rhodz(ind1) = (paprsdn(ind1)-paprsup(ind1))/rg !kg/m2
+        zrho(ind1)  = play(ind1)/temp(ind1)/rd !kg/m3
+        dz(ind1)    = rhodz(ind1)/zrho(ind1) !m : epaisseur de la couche en metre
       END DO
       END DO
 
 
-      sigma1(:,:)=0.
-      sigma2(:,:)=0.
-      qcth(:,:)=0.
-      qcenv(:,:)=0.  
-      qctot(:,:)=0.
-      qlth(:,ind2)=0.
-      qith(:,ind2)=0.
-      wiceth(:,ind2)=0.
-      rneb(:,:)=0.
-      qcloud(:)=0.
-      cth(:,:)=0.
-      cenv(:,:)=0.
-      ctot(:,:)=0.
-      cth_vol(:,:)=0.
-      cenv_vol(:,:)=0.
-      ctot_vol(:,:)=0.
-      falseklon(:)=.false.
-      qsatmmussig1=0.
-      qsatmmussig2=0.
-      rdd=287.04
-      cppd=1005.7
-      pi=3.14159 
-      sqrt2pi=sqrt(2.*pi)
+      icefracth(:)=0.
+      icefracenv(:)=0.
+      sqrt2pi=sqrt(2.*rpi)
       sqrt2=sqrt(2.)
-      sqrtpi=sqrt(pi)
-      icefrac(:,ind2)=0.
-      althl=0.
-      althi=0.
-      athl=0.
-      athi=0.
-      senvl=0.
-      senvi=0.
-      sthi=0.
-      sthl=0.
-      sthil=0.
+      sqrtpi=sqrt(rpi)
+      icefrac(:)=0.
 
 !-------------------------------------------------------------------------------
@@ -1669,10 +1621,9 @@
 !-------------------------------------------------------------------------------
 
-      temp_lim=RTT-40.0
 
       DO ind1=1,klon
-            IF ((temp(ind1,ind2) .LT. RTT) .AND. (temp(ind1,ind2) .GT. temp_lim) &
-            .AND. (iflag_mpc_bl .GE. 1) .AND. (ind2<=klev-2)  &
-            .AND. (ztv(ind1,1).GT.ztv(ind1,2)) .AND.(fraca(ind1,ind2).GT.1.e-10)) THEN
+            IF ((temp(ind1) .LT. RTT) .AND. (temp(ind1) .GT. temp_nowater) &
+            .AND. (ind2<=klev-2)  &
+            .AND. (frac_th(ind1).GT.min_frac_th_cld)) THEN
                 mpc_bl_points(ind1,ind2)=1
             ELSE
@@ -1686,105 +1637,71 @@
 !-------------------------------------------------------------------------------  
 
-! calculation of temperature, humidity and saturation specific humidity
-
-Tbefenv(:)=zthl(:,ind2)*zpspsk(:,ind2)
-Tbefth(:)=ztla(:,ind2)*zpspsk(:,ind2)
-Tbefenvonly(:)=temp(:,ind2)
-rhoth(:)=paprs(:,ind2)/Tbefth(:)/RD
-
-zqenv(:)=(po(:)-fraca(:,ind2)*zqta(:,ind2))/(1.-fraca(:,ind2)) !qt = a*qtth + (1-a)*qtenv 
-zqth(:)=zqta(:,ind2)
-zqenvonly(:)=po(:)
-
-
-CALL CALC_QSAT_ECMWF(klon,Tbefenvonly,zqenvonly,paprs(:,ind2),RTT,0,.false.,zqsatenvonly,dqsatenv)
-
-CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,0,.false.,zqsatenv,dqsatenv)
-CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,1,.false.,qslenv,dqsatenv)
-
-CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,1,.false.,qslth,dqsatth)
-CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,2,.false.,qsith(:,ind2),dqsatth)
-CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,0,.false.,zqsatth,dqsatth)
-
-
-  DO ind1=1,klon
-
-
-    IF ((ztv(ind1,1).GT.ztv(ind1,2)).AND.(fraca(ind1,ind2).GT.1.e-10)) THEN !Thermal and environnement
-
+! initialisations and calculation of temperature, humidity and saturation specific humidity
+
+DO ind1=1,klon
+  
+   rhodz(ind1)   = (paprsdn(ind1)-paprsup(ind1))/rg  ! kg/m2
+   rho(ind1)     = play(ind1)/temp(ind1)/rd          ! kg/m3
+   dz(ind1)      = rhodz(ind1)/rho(ind1)             ! m : epaisseur de la couche en metre
+   Tbefenv(ind1) = temp(ind1)
+   thetal_env(ind1) = Tbefenv(ind1)/zpspsk(ind1)
+   Tbefth(ind1)  = thetal_th(ind1)*zpspsk(ind1)
+   rhoth(ind1)   = play(ind1)/Tbefth(ind1)/RD
+   qt_env(ind1)  = (qt(ind1)-frac_th(ind1)*qt_th(ind1))/(1.-frac_th(ind1)) !qt = a*qtth + (1-a)*qtenv 
+
+ENDDO
+
+! Calculation of saturation specific humidity
+CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,1,.false.,qslenv,dqsatenv)
+CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,2,.false.,qsienv,dqsatenv)
+CALL CALC_QSAT_ECMWF(klon,Tbefth,qt_th,play,RTT,1,.false.,qslth,dqsatth)
+
+
+DO ind1=1,klon
+
+
+    IF (frac_th(ind1).GT.min_frac_th_cld) THEN !Thermal and environnement
+
+! unlike in the other cloudth routine, 
+! We consider distributions of the saturation deficit WRT liquid 
+! at positive AND negative celcius temperature
+! subsequently, cloud fraction corresponds to the part of the pdf corresponding
+! to superstauration with respect to liquid WHATEVER the temperature
 
 ! Environment:
 
 
-        IF (Tbefenv(ind1) .GE. RTT) THEN
-               Lv=RLVTT
-        ELSE
-               Lv=RLSTT
-        ENDIF
-        
-
-        alenv=(0.622*Lv*zqsatenv(ind1))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
-        aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
-        senv=aenv*(po(ind1)-zqsatenv(ind1))                          !s, p84
+        alenv=(0.622*RLVTT*qslenv(ind1))/(rd*thetal_env(ind1)**2)     
+        aenv=1./(1.+(alenv*RLVTT/rcpd))                             
+        senv=aenv*(qt_env(ind1)-qslenv(ind1))                           
      
-        ! For MPCs:
-        IF (mpc_bl_points(ind1,ind2) .EQ. 1) THEN
-        alenvl=(0.622*RLVTT*qslenv(ind1))/(rdd*zthl(ind1,ind2)**2)     
-        aenvl=1./(1.+(alenv*Lv/cppd))         
-        senvl=aenvl*(po(ind1)-qslenv(ind1))    
-        senvi=senv
-        ENDIF
-
 
 ! Thermals:
 
 
-        IF (Tbefth(ind1) .GE. RTT) THEN
-            Lv=RLVTT
-        ELSE
-            Lv=RLSTT
-        ENDIF
-
-        
-        alth=(0.622*Lv*zqsatth(ind1))/(rdd*ztla(ind1,ind2)**2)       
-        ath=1./(1.+(alth*Lv/cppd))                                                         
-        sth=ath*(zqta(ind1,ind2)-zqsatth(ind1))                      
-
-       ! For MPCs:
-        IF (mpc_bl_points(ind1,ind2) .GT. 0) THEN
-         althi=alth
-         althl=(0.622*RLVTT*qslth(ind1))/(rdd*ztla(ind1,ind2)**2)                    
-         athl=1./(1.+(alth*RLVTT/cppd)) 
-         athi=alth     
-         sthl=athl*(zqta(ind1,ind2)-qslth(ind1))   
-         sthi=athi*(zqta(ind1,ind2)-qsith(ind1,ind2))  
-         sthil=athi*(zqta(ind1,ind2)-qslth(ind1))
-        ENDIF     
-
-
-
-!-------------------------------------------------------------------------------
-!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
-!  Rq: in this subroutine, we assume iflag_clouth_vert .EQ. 3 
-!-------------------------------------------------------------------------------
-
-        IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC
-
-       ! Standard deviation of the distributions
-
-           sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
-           &                (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
+        alth=(0.622*RLVTT*qslth(ind1))/(rd*thetal_th(ind1)**2)       
+        ath=1./(1.+(alth*RLVTT/rcpd))                                                         
+        sth=ath*(qt_th(ind1)-qslth(ind1))                      
+
+
+!       IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC
+
+
+! Standard deviation of the distributions
+
+           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*po(ind1)
+             sigma1s_ratqs = cloudth_ratqsmin*qt(ind1)
            ELSE
-             sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
+             sigma1s_ratqs = ratqs(ind1)*qt(ind1)
            ENDIF
  
            sigma1s = sigma1s_fraca + sigma1s_ratqs
            IF (iflag_ratqs.eq.11) then
-              sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
+              sigma1s = ratqs(ind1)*qt(ind1)*aenv
            ENDIF
-           sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
+           sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((frac_th(ind1)+0.02)**sigma2s_power))+0.002*qt_th(ind1)
 
 
@@ -1801,12 +1718,12 @@
            exp_xth2 = exp(-1.*xth2**2)
       
-      !surface CF
-
-           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) 
-
-
-      !volume CF and condensed water
+!surface CF
+
+           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) 
+
+
+!volume CF, condensed water and ice fraction
 
             !environnement
@@ -1816,6 +1733,6 @@
 
             IF (deltasenv .LT. 1.e-10) THEN 
-              qcenv(ind1,ind2)=IntJ
-              cenv_vol(ind1,ind2)=IntJ_CF
+              qcenv(ind1)=IntJ
+              cenv_vol(ind1)=IntJ_CF
             ELSE
               IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
@@ -1824,6 +1741,10 @@
               IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
               IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
-              qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
-              cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
+              qcenv(ind1)=IntJ+IntI1+IntI2+IntI3
+              cenv_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
+              IF (Tbefenv(ind1) .LT. temp_nowater) THEN
+                  ! freeze all droplets in cirrus temperature regime
+                  icefracenv(ind1)=1. 
+              ENDIF
             ENDIF
              
@@ -1836,6 +1757,6 @@
      
             IF (deltasth .LT. 1.e-10) THEN
-              qcth(ind1,ind2)=IntJ
-              cth_vol(ind1,ind2)=IntJ_CF
+              qcth(ind1)=IntJ
+              cth_vol(ind1)=IntJ_CF
             ELSE
               IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
@@ -1844,288 +1765,58 @@
               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
+              qcth(ind1)=IntJ+IntI1+IntI2+IntI3
+              cth_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
+              IF (Tbefth(ind1) .LT. temp_nowater) THEN
+                  ! freeze all droplets in cirrus temperature regime
+                  icefracth(ind1)=1.
+              ENDIF
+
             ENDIF
 
-              qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)
-              ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
-              
-
-            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)
+            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)=qslenv(ind1)
                 qincloud(ind1)=0.
+                icefrac(ind1)=0.
             ELSE               
-                qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
-                qincloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)
+                qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
+                qincloud(ind1)=qctot(ind1)/ctot(ind1)
+                IF (qctot(ind1) .GT. 0) THEN
+                  icefrac(ind1)=(frac_th(ind1)*qcth(ind1)*icefracth(ind1)+(1.-1.*frac_th(ind1))*qcenv(ind1)*icefracenv(ind1))/qctot(ind1)
+                  icefrac(ind1)=max(min(1.,icefrac(ind1)), 0.)
+                ENDIF
             ENDIF 
 
 
-        ELSE ! mpc_bl_points>0
-
-            ! Treat boundary layer mixed phase clouds
-            
-            ! thermals
-            !=========
-
-            ! ice phase
-            !...........
-
-            qiceth=0.
-            deltazlev_mpc=dz(ind1,:)
-            temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
-            pres_mpc=pplay(ind1,:)
-            fraca_mpc=fraca(ind1,:)
-            snowf_mpc=snowflux(ind1,:)
-            iflag_topthermals=0
-            IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0))  THEN
-                iflag_topthermals = 1
-            ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) & 
-                    .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
-                iflag_topthermals = 2
-            ELSE 
-                iflag_topthermals = 0
-            ENDIF
-
-            CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), &
-                                   qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:))
-
-            ! qmax calculation
-            sigma2s=(sigma2s_factor*((MAX((sthl-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
-            deltasth=athl*vert_alpha_th*sigma2s
-            xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
-            xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
-            exp_xth1 = exp(-1.*xth1**2)
-            exp_xth2 = exp(-1.*xth2**2)
-            IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
-            IntJ_CF=0.5*(1.-1.*erf(xth2))
-            IntI1=(((sthl+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*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
-            IntI1_CF=((sthl+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
-            IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
-            qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)
-
-
-            ! Liquid phase
-            !................
-            ! We account for the effect of ice crystals in thermals on sthl
-            ! and on the width of the distribution
-
-            sthlc=sthl*1./(1.+C_mpc*qith(ind1,ind2))  &
-                + (1.-1./(1.+C_mpc*qith(ind1,ind2))) * athl*(qsith(ind1,ind2)-qslth(ind1))  
-
-            sigma2sc=(sigma2s_factor*((MAX((sthlc-senvl),0.)**2)**0.5) &
-                 /((fraca(ind1,ind2)+0.02)**sigma2s_power)) &
-                 +0.002*zqta(ind1,ind2)
-            deltasthc=athl*vert_alpha_th*sigma2sc
-     
-            
-            xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
-            xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)           
-            exp_xth1 = exp(-1.*xth1**2)
-            exp_xth2 = exp(-1.*xth2**2)
-            IntJ=0.5*sthlc*(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
-            IntJ_CF=0.5*(1.-1.*erf(xth2))
-            IntI1=(((sthlc+deltasthc)**2+(sigma2sc)**2)/(8*deltasthc))*(erf(xth2)-erf(xth1))
-            IntI2=(sigma2sc**2/(4*deltasthc*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
-            IntI3=((sqrt2*sigma2sc*(sthlc+deltasthc))/(4*sqrtpi*deltasthc))*(exp_xth1-exp_xth2)
-            IntI1_CF=((sthlc+deltasthc)*(erf(xth2)-erf(xth1)))/(4*deltasthc)
-            IntI3_CF=(sqrt2*sigma2sc*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasthc)
-            qliqth=IntJ+IntI1+IntI2+IntI3
-            
-            qlth(ind1,ind2)=MAX(0.,qliqth)
-
-            ! Condensed water
-            
-            qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
-
-
-            ! consistency with subgrid distribution
-           
-             IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
-                 fraci=qith(ind1,ind2)/qcth(ind1,ind2)
-                 qcth(ind1,ind2)=qmax
-                 qith(ind1,ind2)=fraci*qmax
-                 qlth(ind1,ind2)=(1.-fraci)*qmax
-             ENDIF
-
-            ! Cloud Fraction
-            !...............
-            ! calculation of qbase which is the value of the water vapor within mixed phase clouds
-            ! such that the total water in cloud = qbase+qliqth+qiceth
-            ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
-            ! sbase and qbase calculation (note that sbase is wrt liq so negative)
-            ! look for an approximate solution with iteration
-            
-            ttarget=qcth(ind1,ind2)
-            mini= athl*(qsith(ind1,ind2)-qslth(ind1))
-            maxi= 0. !athl*(3.*sqrt(sigma2s))
-            niter=20
-            pas=(maxi-mini)/niter
-            stmp=mini
-            sbase=stmp
-            coutref=1.E6
-            DO iter=1,niter
-                cout=ABS(sigma2s/SQRT(2.*RPI)*EXP(-((sthl-stmp)/sigma2s)**2)+(sthl-stmp)/SQRT(2.)*(1.-erf(-(sthl-stmp)/sigma2s)) &
-                     + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
-               IF (cout .LT. coutref) THEN
-                     sbase=stmp
-                     coutref=cout
-                ELSE
-                     stmp=stmp+pas
-                ENDIF
-            ENDDO
-            qbase=MAX(0., sbase/athl+qslth(ind1))
-
-            ! surface cloud fraction in thermals
-            cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
-            cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)
-
-
-            !volume cloud fraction in thermals
-            !to be checked
-            xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
-            xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)           
-            exp_xth1 = exp(-1.*xth1**2)
-            exp_xth2 = exp(-1.*xth2**2)
-
-            IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
-            IntJ_CF=0.5*(1.-1.*erf(xth2))
-     
-            IF (deltasth .LT. 1.e-10) THEN
-              cth_vol(ind1,ind2)=IntJ_CF
-            ELSE
-              IntI1=(((sthl+deltasth-sbase)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
-              IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
-              IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
-              IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
-              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
-              cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
-            ENDIF
-              cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
-
-
-
-            ! Environment
-            !=============
-            ! In the environment/downdrafts, only liquid clouds
-            ! See Shupe et al. 2008, JAS
-
-            ! standard deviation of the distribution in the environment
-            sth=sthl
-            senv=senvl
-            sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
-                          &                (1-fraca(ind1,ind2))*(MAX((sth-senv),0.)**2)**0.5
-            ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
-            ! in the environement
-
-            sigma1s_ratqs=1E-10
-            IF (cloudth_ratqsmin>0.) THEN
-                sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
-            ENDIF
-
-            sigma1s = sigma1s_fraca + sigma1s_ratqs
-            IF (iflag_ratqs.eq.11) then
-               sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
-            ENDIF
-            IF (iflag_ratqs.eq.11) then
-              sigma1s = ratqs(ind1,ind2)*po(ind1)*aenvl
-            ENDIF
-            deltasenv=aenvl*vert_alpha*sigma1s
-            xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
-            xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
-            exp_xenv1 = exp(-1.*xenv1**2)
-            exp_xenv2 = exp(-1.*xenv2**2)
-
-            !surface CF
-            cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
-
-            !volume CF and condensed water
-            IntJ=0.5*senvl*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
-            IntJ_CF=0.5*(1.-1.*erf(xenv2))
-
-            IF (deltasenv .LT. 1.e-10) THEN
-              qcenv(ind1,ind2)=IntJ
-              cenv_vol(ind1,ind2)=IntJ_CF
-            ELSE
-              IntI1=(((senvl+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=((senvl+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
-              IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
-              qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
-              cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
-            ENDIF
-
-            qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
-            cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)
-
-
-            
-            ! Thermals + environment
-            ! ======================= 
-            ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) 
-            qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)
-            ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
-            IF (qcth(ind1,ind2) .GT. 0) THEN
-               icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) &
-                    /(fraca(ind1,ind2)*qcth(ind1,ind2) &
-                    +(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2))
-                icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
-            ELSE
-                icefrac(ind1,ind2)=0.
-            ENDIF
-
-            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.
-                qincloud(ind1)=0.
-                qcloud(ind1)=zqsatenv(ind1)
-            ELSE               
-                qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
-                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
-                qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
-                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
-            ENDIF 
-
-        ENDIF ! mpc_bl_points
-
+!        ELSE ! mpc_bl_points>0
 
     ELSE  ! gaussian for environment only
 
       
-
-
-        IF (Tbefenvonly(ind1) .GE. RTT) THEN
-                Lv=RLVTT
-        ELSE
-                Lv=RLSTT
-        ENDIF
-        
-
-        zthl(ind1,ind2)=temp(ind1,ind2)*(101325./paprs(ind1,ind2))**(rdd/cppd)
-        alenv=(0.622*Lv*zqsatenvonly(ind1))/(rdd*zthl(ind1,ind2)**2)
-        aenv=1./(1.+(alenv*Lv/cppd))
-        senv=aenv*(po(ind1)-zqsatenvonly(ind1))
+        alenv=(0.622*RLVTT*qslenv(ind1))/(rd*thetal_env(ind1)**2)
+        aenv=1./(1.+(alenv*RLVTT/rcpd))
+        senv=aenv*(qt_env(ind1)-qslenv(ind1))
         sth=0.
-      
-        sigma1s=ratqs(ind1,ind2)*zqenvonly(ind1)
+        sigma1s=ratqs(ind1)*qt_env(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)
-        qctot(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)=zqsatenvonly(ind1)
+        cenv(ind1)=0.5*(1.-1.*erf(xenv))
+        ctot(ind1)=0.5*(1.+1.*erf(xenv))
+        ctot_vol(ind1)=ctot(ind1)
+        qctot(ind1)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1))
+
+
+        IF (ctot(ind1).LT.min_neb_th) THEN
+          ctot(ind1)=0.
+          qcloud(ind1)=qslenv(ind1)
           qincloud(ind1)=0.
         ELSE 
-          qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenvonly(ind1)
-          qincloud(ind1)=MAX(qctot(ind1,ind2)/ctot(ind1,ind2),0.)
+          qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
+          qincloud(ind1)=MAX(qctot(ind1)/ctot(ind1),0.)
         ENDIF
  
@@ -2134,15 +1825,12 @@
 
     ! 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
+    cloudth_senv(ind1) = senv
+    cloudth_sth(ind1) = sth
+    cloudth_sigmaenv(ind1) = sigma1s
+    cloudth_sigmath(ind1) = sigma2s
 
 
     ENDDO       !loop on klon
 
-    ! Calcule ice fall velocity in thermals
-
-    CALL FALLICE_VELOCITY(klon,qith(:,ind2),Tbefth(:),rhoth(:),paprs(:,ind2),falseklon(:),wiceth(:,ind2))
 
 RETURN
@@ -2150,4 +1838,238 @@
 
 END SUBROUTINE cloudth_mpc
+
+
+
+!        ELSE ! mpc_bl_points>0
+!
+!            ! Treat boundary layer mixed phase clouds
+!            
+!            ! thermals
+!            !=========
+!
+!           ! ice phase
+!           !...........
+!
+!            qiceth=0.
+!            deltazlev_mpc=dz(ind1,:)
+!            temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
+!            pres_mpc=pplay(ind1,:)
+!            fraca_mpc=fraca(ind1,:)
+!            snowf_mpc=snowflux(ind1,:)
+!            iflag_topthermals=0
+!            IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0))  THEN
+!                iflag_topthermals = 1
+!            ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) & 
+!                    .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
+!                iflag_topthermals = 2
+!            ELSE 
+!                iflag_topthermals = 0
+!            ENDIF
+!
+!            CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), &
+!                                   qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:))
+!
+!            ! qmax calculation
+!            sigma2s=(sigma2s_factor*((MAX((sthl-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
+!            deltasth=athl*vert_alpha_th*sigma2s
+!            xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
+!            xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
+!            exp_xth1 = exp(-1.*xth1**2)
+!            exp_xth2 = exp(-1.*xth2**2)
+!            IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
+!            IntJ_CF=0.5*(1.-1.*erf(xth2))
+!            IntI1=(((sthl+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*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
+!            IntI1_CF=((sthl+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
+!            IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
+!            qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)
+!
+!
+!            ! Liquid phase
+!            !................
+!            ! We account for the effect of ice crystals in thermals on sthl
+!            ! and on the width of the distribution
+!
+!            sthlc=sthl*1./(1.+C_mpc*qith(ind1,ind2))  &
+!                + (1.-1./(1.+C_mpc*qith(ind1,ind2))) * athl*(qsith(ind1,ind2)-qslth(ind1))  
+!
+!            sigma2sc=(sigma2s_factor*((MAX((sthlc-senvl),0.)**2)**0.5) &
+!                 /((fraca(ind1,ind2)+0.02)**sigma2s_power)) &
+!                 +0.002*zqta(ind1,ind2)
+!            deltasthc=athl*vert_alpha_th*sigma2sc
+!     
+!            
+!            xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
+!            xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)           
+!            exp_xth1 = exp(-1.*xth1**2)
+!            exp_xth2 = exp(-1.*xth2**2)
+!            IntJ=0.5*sthlc*(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
+!            IntJ_CF=0.5*(1.-1.*erf(xth2))
+!            IntI1=(((sthlc+deltasthc)**2+(sigma2sc)**2)/(8*deltasthc))*(erf(xth2)-erf(xth1))
+!            IntI2=(sigma2sc**2/(4*deltasthc*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
+!            IntI3=((sqrt2*sigma2sc*(sthlc+deltasthc))/(4*sqrtpi*deltasthc))*(exp_xth1-exp_xth2)
+!            IntI1_CF=((sthlc+deltasthc)*(erf(xth2)-erf(xth1)))/(4*deltasthc)
+!            IntI3_CF=(sqrt2*sigma2sc*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasthc)
+!            qliqth=IntJ+IntI1+IntI2+IntI3
+!            
+!            qlth(ind1,ind2)=MAX(0.,qliqth)
+!
+!            ! Condensed water
+!            
+!            qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
+!
+!
+!            ! consistency with subgrid distribution
+!           
+!             IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
+!                 fraci=qith(ind1,ind2)/qcth(ind1,ind2)
+!                 qcth(ind1,ind2)=qmax
+!                 qith(ind1,ind2)=fraci*qmax
+!                 qlth(ind1,ind2)=(1.-fraci)*qmax
+!             ENDIF
+!
+!            ! Cloud Fraction
+!            !...............
+!            ! calculation of qbase which is the value of the water vapor within mixed phase clouds
+!            ! such that the total water in cloud = qbase+qliqth+qiceth
+!            ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
+!            ! sbase and qbase calculation (note that sbase is wrt liq so negative)
+!            ! look for an approximate solution with iteration
+!            
+!            ttarget=qcth(ind1,ind2)
+!            mini= athl*(qsith(ind1,ind2)-qslth(ind1))
+!            maxi= 0. !athl*(3.*sqrt(sigma2s))
+!            niter=20
+!            pas=(maxi-mini)/niter
+!            stmp=mini
+!            sbase=stmp
+!            coutref=1.E6
+!            DO iter=1,niter
+!                cout=ABS(sigma2s/SQRT(2.*RPI)*EXP(-((sthl-stmp)/sigma2s)**2)+(sthl-stmp)/SQRT(2.)*(1.-erf(-(sthl-stmp)/sigma2s)) &
+!                     + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
+!               IF (cout .LT. coutref) THEN
+!                     sbase=stmp
+!                     coutref=cout
+!                ELSE
+!                     stmp=stmp+pas
+!                ENDIF
+!            ENDDO
+!            qbase=MAX(0., sbase/athl+qslth(ind1))
+!
+!            ! surface cloud fraction in thermals
+!            cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
+!            cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)
+!
+!
+!            !volume cloud fraction in thermals
+!            !to be checked
+!            xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
+!            xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)           
+!            exp_xth1 = exp(-1.*xth1**2)
+!            exp_xth2 = exp(-1.*xth2**2)
+!
+!            IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
+!            IntJ_CF=0.5*(1.-1.*erf(xth2))
+!     
+!            IF (deltasth .LT. 1.e-10) THEN
+!              cth_vol(ind1,ind2)=IntJ_CF
+!            ELSE
+!              IntI1=(((sthl+deltasth-sbase)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
+!              IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
+!              IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
+!              IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
+!              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
+!              cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
+!            ENDIF
+!              cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
+!
+!
+!
+!            ! Environment
+!            !=============
+!            ! In the environment/downdrafts, only liquid clouds
+!            ! See Shupe et al. 2008, JAS
+!
+!            ! standard deviation of the distribution in the environment
+!            sth=sthl
+!            senv=senvl
+!            sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
+!                          &                (1-fraca(ind1,ind2))*(MAX((sth-senv),0.)**2)**0.5
+!            ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
+!            ! in the environement
+!
+!            sigma1s_ratqs=1E-10
+!            IF (cloudth_ratqsmin>0.) THEN
+!                sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
+!            ENDIF
+!
+!            sigma1s = sigma1s_fraca + sigma1s_ratqs
+!            IF (iflag_ratqs.eq.11) then
+!               sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
+!            ENDIF
+!            IF (iflag_ratqs.eq.11) then
+!              sigma1s = ratqs(ind1,ind2)*po(ind1)*aenvl
+!            ENDIF
+!            deltasenv=aenvl*vert_alpha*sigma1s
+!            xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
+!            xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
+!            exp_xenv1 = exp(-1.*xenv1**2)
+!            exp_xenv2 = exp(-1.*xenv2**2)
+!
+!            !surface CF
+!            cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
+!
+!            !volume CF and condensed water
+!            IntJ=0.5*senvl*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
+!            IntJ_CF=0.5*(1.-1.*erf(xenv2))
+!
+!            IF (deltasenv .LT. 1.e-10) THEN
+!              qcenv(ind1,ind2)=IntJ
+!              cenv_vol(ind1,ind2)=IntJ_CF
+!            ELSE
+!              IntI1=(((senvl+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=((senvl+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
+!              IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
+!              qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
+!              cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
+!            ENDIF
+!
+!            qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
+!            cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)
+!
+!
+!            
+!            ! Thermals + environment
+!            ! ======================= 
+!            ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) 
+!            qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)
+!            ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
+!            IF (qcth(ind1,ind2) .GT. 0) THEN
+!               icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) &
+!                    /(fraca(ind1,ind2)*qcth(ind1,ind2) &
+!                    +(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2))
+!                icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
+!            ELSE
+!                icefrac(ind1,ind2)=0.
+!            ENDIF
+!
+!            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.
+!                qincloud(ind1)=0.
+!                qcloud(ind1)=zqsatenv(ind1)
+!            ELSE               
+!                qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
+!                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
+!                qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
+!                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
+!            ENDIF 
+!
+!        ENDIF ! mpc_bl_points
+
+
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Index: /LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90	(revision 4909)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90	(revision 4910)
@@ -12,6 +12,6 @@
      radocond, radicefrac, rain, snow,                  &
      frac_impa, frac_nucl, beta,                        &
-     prfl, psfl, rhcl, zqta, fraca,                     &
-     ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
+     prfl, psfl, rhcl, qta, fraca,                     &
+     tv, pspsk, tla, thl, iflag_cld_th,             &
      iflag_ice_thermo, ok_ice_sursat, qsatl, qsats,     &
      distcltop,temp_cltop,                              &
@@ -104,9 +104,9 @@
 USE lmdz_lscp_ini, ONLY : prt_level, lunout
 USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
-USE lmdz_lscp_ini, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
+USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
 USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min 
 USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
 USE lmdz_lscp_ini, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo
-USE lmdz_lscp_ini, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc
+USE lmdz_lscp_ini, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc, min_frac_th_cld
 USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
 USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_fonte_lscp
@@ -138,14 +138,14 @@
   !Inputs associated with thermal plumes
 
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztv             ! virtual potential temperature [K]
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zqta            ! specific humidity within thermals [kg/kg]
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca           ! fraction of thermals within the mesh [-]
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zpspsk          ! exner potential (p/100000)**(R/cp) 
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztla            ! liquid temperature within thermals [K]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv             ! virtual potential temperature [K]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta            ! specific humidity within thermals [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca          ! fraction of thermals within the mesh [-]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk          ! exner potential (p/100000)**(R/cp) 
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla            ! liquid temperature within thermals [K]
 
   ! INPUT/OUTPUT variables
   !------------------------
   
-  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: zthl         ! liquid potential temperature [K]
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl          ! liquid potential temperature [K]
   REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs        ! function of pressure that sets the large-scale
 
@@ -236,5 +236,5 @@
   REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2
   REAL :: erf
-  REAL, DIMENSION(klon,klev) :: icefrac_mpc
+  REAL, DIMENSION(klon) :: zfice_th
   REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
   REAL, DIMENSION(klon) :: zrfl, zrfln
@@ -713,12 +713,14 @@
 
         IF (iflag_cld_th.GE.5) THEN
-
-             IF (iflag_mpc_bl .LT. 1) THEN 
-
-                IF (iflag_cloudth_vert.LE.2) THEN
-
-                    CALL cloudth(klon,klev,k,ztv,                  &
-                         zq,zqta,fraca,                            &
-                         qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
+               ! Cloud cover and content in meshes affected by shallow convection, 
+               ! are retrieved from a bi-gaussian distribution of the saturation deficit
+               ! following Jam et al. 2013
+
+               IF (iflag_cloudth_vert.LE.2) THEN
+                  ! Old version of Arnaud Jam
+
+                    CALL cloudth(klon,klev,k,tv,                  &
+                         zq,qta,fraca,                            &
+                         qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
                          ratqs,zqs,temp,                              &
                          cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
@@ -726,35 +728,33 @@
 
                 ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
-
-                    CALL cloudth_v3(klon,klev,k,ztv,                        &
-                         zq,zqta,fraca,                                     &
-                         qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
+                   ! Default version of Arnaud Jam
+                        
+                    CALL cloudth_v3(klon,klev,k,tv,                        &
+                         zq,qta,fraca,                                     &
+                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
                          ratqs,zqs,temp, &
                          cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
 
-                !Jean Jouhaud's version, Decembre 2018 
-                !-------------------------------------
 
                 ELSEIF (iflag_cloudth_vert.EQ.6) THEN
-
-                    CALL cloudth_v6(klon,klev,k,ztv,                        &
-                         zq,zqta,fraca,                                     &
-                         qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
+                   ! Jean Jouhaud's version, with specific separation between surface and volume
+                   ! cloud fraction Decembre 2018 
+
+                    CALL cloudth_v6(klon,klev,k,tv,                        &
+                         zq,qta,fraca,                                     &
+                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
                          ratqs,zqs,temp, &
                          cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
 
-                ENDIF
-
-         ELSE 
-                ! cloudth_v3 + cold microphysical considerations 
-                ! + stationary mixed-phase cloud model
-
-                    CALL cloudth_mpc(klon,klev,k,iflag_mpc_bl,mpc_bl_points, &
-                         temp,ztv,zq,zqta,fraca, zpspsk,                      &
-                         paprs,pplay,ztla,zthl,ratqs,zqs,psfl,             &
-                         qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol,    &
-                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
-
-         ENDIF ! iflag_mpc_bl
+                ELSEIF (iflag_cloudth_vert .EQ. 7) THEN
+                   ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
+                   ! for boundary-layer mixed phase clouds following Vignon et al.  
+                    CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
+                                     pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
+                                     ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
+                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
+
+                ENDIF
+
 
                 DO i=1,klon
@@ -774,5 +774,5 @@
             
                 ! lognormal distribution when no thermals 
-            lognormale = fraca(:,k) < 1e-10
+            lognormale = fraca(:,k) < min_frac_th_cld
 
         ELSE
@@ -934,7 +934,6 @@
             temp_cltop(:,k)=temp_cltop1D(:)
         ENDIF   
-
         ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs)
-        CALL icefrac_lscp(klon,zt,iflag_ice_thermo,distcltop1D,temp_cltop1D,zfice, dzfice)
+        CALL icefrac_lscp(klon,zt,iflag_ice_thermo,distcltop1D,temp_cltop1D,zfice,dzfice)
 
 
@@ -943,4 +942,5 @@
 
             IF (mpc_bl_points(i,k) .GT. 0) THEN
+               
                 zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
                 ! following line is very strange and probably wrong
@@ -948,5 +948,5 @@
                 ! water vapor update and partition function if thermals
                 zq(i) = zq(i) - zcond(i)        
-                zfice(i)=icefrac_mpc(i,k)
+                zfice(i)=zfice_th(i)
 
             ELSE
@@ -1314,4 +1314,5 @@
     ENDIF
 
+
     ENDIF ! ok_poprecip
 
Index: /LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90	(revision 4909)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90	(revision 4910)
@@ -12,6 +12,12 @@
   !$OMP THREADPRIVATE(ok_bug_fonte_lscp)
 
-  REAL, SAVE, PROTECTED :: seuil_neb=0.001      ! cloud fraction threshold: a cloud really exists when exceeded
+  REAL, SAVE, PROTECTED :: seuil_neb=0.001      ! cloud fraction threshold: a cloud can precipitate when exceeded
   !$OMP THREADPRIVATE(seuil_neb)
+
+  REAL, SAVE, PROTECTED :: min_neb_th=1e-10      ! a cloud produced by bi-gaussian really exists when exceeded
+  !$OMP THREADPRIVATE(min_neb_th)
+
+  REAL, SAVE, PROTECTED :: min_frac_thermals=1.e-10   ! minimum thermals fraction for use of bigaussian distribution
+  !$OMP THREADPRIVATE(min_frac_thermals)
 
   INTEGER, SAVE :: lunout, prt_level            ! Logical unit number and level for standard output
@@ -43,7 +49,7 @@
   !$OMP THREADPRIVATE(a_tr_sca)
   
-  INTEGER, SAVE, PROTECTED ::  iflag_mpc_bl=0   ! flag to activate boundary layer mixed phase cloud param
-  !$OMP THREADPRIVATE(iflag_mpc_bl)
-  
+  REAL, SAVE, PROTECTED ::  min_frac_th_cld=1.e-10   ! minimum thermal fraction to compute a thermal cloud fraction
+  !$OMP THREADPRIVATE(min_frac_th_cld)
+
   LOGICAL, SAVE, PROTECTED :: ok_radocond_snow=.false. ! take into account the mass of ice precip in the cloud ice content seen by radiation
   !$OMP THREADPRIVATE(ok_radocond_snow)
@@ -260,5 +266,4 @@
     CALL getin_p('iflag_evap_prec',iflag_evap_prec)
     CALL getin_p('seuil_neb',seuil_neb)
-    CALL getin_p('iflag_mpc_bl',iflag_mpc_bl)
     CALL getin_p('ok_radocond_snow',ok_radocond_snow) 
     CALL getin_p('t_glace_max',t_glace_max)
@@ -320,5 +325,4 @@
     WRITE(lunout,*) 'lscp_ini, iflag_evap_prec:', iflag_evap_prec
     WRITE(lunout,*) 'lscp_ini, seuil_neb:', seuil_neb
-    WRITE(lunout,*) 'lscp_ini, iflag_mpc_bl:', iflag_mpc_bl
     WRITE(lunout,*) 'lscp_ini, ok_radocond_snow:', ok_radocond_snow
     WRITE(lunout,*) 'lscp_ini, t_glace_max:', t_glace_max
