Index: LMDZ6/trunk/libf/phylmd/clesphys.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/clesphys.h	(revision 3998)
+++ LMDZ6/trunk/libf/phylmd/clesphys.h	(revision 3999)
@@ -93,4 +93,5 @@
        LOGICAL :: adjust_tropopause
        LOGICAL :: ok_daily_climoz
+       LOGICAL :: ok_new_lscp
 ! flag to bypass or not the phytrac module
        INTEGER :: iflag_phytrac
@@ -141,5 +142,5 @@
      &     , ok_chlorophyll,ok_conserv_q, adjust_tropopause             &
      &     , ok_daily_climoz, ok_all_xml, ok_lwoff                      &
-     &     , iflag_phytrac
+     &     , iflag_phytrac, ok_new_lscp
      
        save /clesphys/
Index: LMDZ6/trunk/libf/phylmd/cloudth_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/cloudth_mod.F90	(revision 3998)
+++ LMDZ6/trunk/libf/phylmd/cloudth_mod.F90	(revision 3999)
@@ -655,5 +655,4 @@
       REAL zqs(ngrid), qcloud(ngrid)
       REAL erf
-
 
 
@@ -911,8 +910,8 @@
       END DO
 
-
 !------------------------------------------------------------------------------
 ! Initialize
 !------------------------------------------------------------------------------
+
       sigma1(:,:)=0.
       sigma2(:,:)=0.
@@ -1013,5 +1012,4 @@
 !zqsatth = qsat thermals
 !ztla = Tl thermals
-
 !------------------------------------------------------------------------------
 ! s standard deviation
@@ -1217,4 +1215,5 @@
       else  ! gaussienne environnement seule
       
+
       zqenv(ind1)=po(ind1)
       Tbef=t(ind1,ind2)
@@ -1534,4 +1533,880 @@
 
 END SUBROUTINE cloudth_v6
+
+
+
+
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points,                     &
+&           temp,ztv,po,zqta,fraca,zpspsk,paprs,pplay,ztla,zthl,            &
+&           ratqs,zqs,snowflux,qcloud,qincloud,icefrac,ctot,ctot_vol)
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+! 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
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+      USE ioipsl_getin_p_mod, ONLY : getin_p
+      USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
+      USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF, ICEFRAC_LSCP
+      USE phys_local_var_mod, ONLY : qlth, qith
+
+      IMPLICIT NONE
+
+#include "YOMCST.h"
+#include "YOETHF.h"
+#include "FCTTRE.h"
+#include "thermcell.h"
+#include "nuage.h"
+      
+
+!------------------------------------------------------------------------------
+! Declaration
+!------------------------------------------------------------------------------
+
+! INPUT/OUTPUT
+
+      INTEGER, INTENT(IN)                         :: klon,klev,ind2
+      INTEGER, DIMENSION(klon,klev), INTENT(INOUT)   ::  mpc_bl_points ! 1 where BL MPC, 0 otherwise
+
+      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]
+
+
+      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)      ::  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]
+
+
+! LOCAL VARIABLES
+
+      INTEGER itap,ind1,l,ig,iter,k 
+      LOGICAL flag_topthermals
+
+
+      REAL zqsatth(klon,klev), zqsatenv(klon,klev)
+      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)   
+
+      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,dqsatdt
+      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 qslth, qsith, qslenv, alenvl, aenvl
+      REAL sthi, sthl, althl, athl
+      REAL senvi, senvl, qbase, sbase, qliqth, qiceth
+      REAL qimax, ttarget, stmp, cout, coutref
+      REAL maxi, mini, pas, temp_lim
+      REAL deltazlev_mpc(klev),qth_mpc(klev), temp_mpc(klev), pres_mpc(klev), fraca_mpc(klev+1), snowf_mpc(klev+1)
+
+      INTEGER, SAVE :: niter=20
+
+      ! Modifty the saturation deficit PDF in thermals 
+      ! in the presence of ice crystals
+      REAL,SAVE :: C_mpc
+      !$OMP THREADPRIVATE(C_mpc)
+      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
+      ! (J Jouhaud, JL Dufresne, JB Madeleine)
+      REAL,SAVE :: vert_alpha, vert_alpha_th
+      !$OMP THREADPRIVATE(vert_alpha, vert_alpha_th)
+      REAL,SAVE :: sigma1s_factor=1.1
+      REAL,SAVE :: sigma1s_power=0.6
+      REAL,SAVE :: sigma2s_factor=0.09
+      REAL,SAVE :: sigma2s_power=0.5
+      REAL,SAVE :: cloudth_ratqsmin=-1.
+      !$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin)
+      INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0
+      !$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs)
+      LOGICAL, SAVE :: firstcall = .TRUE.
+      !$OMP THREADPRIVATE(firstcall)
+
+      CHARACTER (len = 80) :: abort_message
+      CHARACTER (len = 20) :: routname = 'cloudth_mpc'
+
+
+!------------------------------------------------------------------------------
+! Initialisation
+!------------------------------------------------------------------------------
+
+
+! 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
+      END DO
+      END DO
+
+
+      sigma1(:,:)=0.
+      sigma2(:,:)=0.
+      qcth(:,:)=0.
+      qcenv(:,:)=0.  
+      qctot(:,:)=0.
+      qlth(:,ind2)=0.
+      qith(:,ind2)=0.
+      rneb(:,:)=0.
+      qcloud(:)=0.
+      cth(:,:)=0.
+      cenv(:,:)=0.
+      ctot(:,:)=0.
+      cth_vol(:,:)=0.
+      cenv_vol(:,:)=0.
+      ctot_vol(:,:)=0.
+      qsatmmussig1=0.
+      qsatmmussig2=0.
+      rdd=287.04
+      cppd=1005.7
+      pi=3.14159 
+      sqrt2pi=sqrt(2.*pi)
+      sqrt2=sqrt(2.)
+      sqrtpi=sqrt(pi)
+      icefrac(:,ind2)=0.
+
+
+
+      IF (firstcall) THEN
+
+        vert_alpha=0.5
+        CALL getin_p('cloudth_vert_alpha',vert_alpha)
+        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
+        ! The factor used for the thermal is equal to that of the environment
+        !   if nothing is explicitly specified in the def file
+        vert_alpha_th=vert_alpha
+        CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
+        WRITE(*,*) 'cloudth_vert_alpha_th = ', vert_alpha_th
+        ! Factor used in the calculation of sigma1s
+        CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor)
+        WRITE(*,*) 'cloudth_sigma1s_factor = ', sigma1s_factor
+        ! Power used in the calculation of sigma1s
+        CALL getin_p('cloudth_sigma1s_power',sigma1s_power)
+        WRITE(*,*) 'cloudth_sigma1s_power = ', sigma1s_power
+        ! Factor used in the calculation of sigma2s
+        CALL getin_p('cloudth_sigma2s_factor',sigma2s_factor)
+        WRITE(*,*) 'cloudth_sigma2s_factor = ', sigma2s_factor
+        ! Power used in the calculation of sigma2s
+        CALL getin_p('cloudth_sigma2s_power',sigma2s_power)
+        WRITE(*,*) 'cloudth_sigma2s_power = ', sigma2s_power
+        ! Minimum value for the environmental air subgrid water distrib
+        CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin)
+        WRITE(*,*) 'cloudth_ratqsmin = ', cloudth_ratqsmin
+        ! Remove the dependency to ratqs from the variance of the vertical PDF
+        CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs)
+        WRITE(*,*) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs
+        ! Modifies the PDF in thermals when ice crystals are present
+        C_mpc=1.e2
+        CALL getin_p('C_mpc',C_mpc)
+        WRITE(*,*) 'C_mpc = ', C_mpc
+
+        firstcall=.FALSE.
+
+      ENDIF
+
+
+
+!-------------------------------------------------------------------------------
+! Identify grid points with potential mixed-phase conditions
+!-------------------------------------------------------------------------------
+
+      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. 2) .AND. (ind2<=klev-2)  &
+            .AND. (ztv(ind1,1).GT.ztv(ind1,2)) .AND.(fraca(ind1,ind2).GT.1.e-10)) THEN
+                mpc_bl_points(ind1,ind2)=1
+            ELSE
+                mpc_bl_points(ind1,ind2)=0
+            ENDIF
+      ENDDO
+
+
+!-------------------------------------------------------------------------------
+! Thermal fraction calculation and standard deviation of the distribution
+!-------------------------------------------------------------------------------  
+
+    DO ind1=1,klon
+
+
+    IF ((ztv(ind1,1).GT.ztv(ind1,2)).AND.(fraca(ind1,ind2).GT.1.e-10)) THEN !Thermal and environnement
+
+
+! Environment:
+
+        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)
+            
+        CALL CALC_QSAT_ECMWF(Tbef,0.,paprs(ind1,ind2),RTT,0,.false.,qsatbef,dqsatdt)
+        zqsatenv(ind1,ind2)=qsatbef
+
+        IF (Tbef .GE. RTT) THEN
+               Lv=RLVTT
+        ELSE
+               Lv=RLSTT
+        ENDIF
+        
+
+        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
+     
+        ! For MPCs:
+        IF (mpc_bl_points(ind1,ind2) .EQ. 1) THEN
+        CALL CALC_QSAT_ECMWF(Tbef,0.,paprs(ind1,ind2),RTT,1,.false.,qslenv,dqsatdt)
+        alenvl=(0.622*RLVTT*qslenv)/(rdd*zthl(ind1,ind2)**2)     
+        aenvl=1./(1.+(alenv*Lv/cppd))         
+        senvl=aenvl*(po(ind1)-qslenv)    
+        ENDIF
+
+
+! Thermals:
+
+        Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
+        CALL CALC_QSAT_ECMWF(Tbef,0.,paprs(ind1,ind2),RTT,0,.false.,qsatbef,dqsatdt)
+        zqsatth(ind1,ind2)=qsatbef
+
+        IF (Tbef .GE. RTT) THEN
+            Lv=RLVTT
+        ELSE
+            Lv=RLSTT
+        ENDIF
+
+        
+        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))                      
+
+       ! For MPCs:
+        IF (mpc_bl_points(ind1,ind2) .GT. 0) THEN
+         CALL CALC_QSAT_ECMWF(Tbef,0.,paprs(ind1,ind2),RTT,1,.false.,qslth,dqsatdt)
+         CALL CALC_QSAT_ECMWF(Tbef,0.,paprs(ind1,ind2),RTT,2,.false.,qsith,dqsatdt)
+         althl=(0.622*RLVTT*qslth)/(rdd*ztla(ind1,ind2)**2)                    
+         athl=1./(1.+(alth*RLVTT/cppd))                                        
+         sthl=athl*(zqta(ind1,ind2)-qslth)   
+         sthi=athl*(zqta(ind1,ind2)-qsith)  
+        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
+
+           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)
+
+
+           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)
+      
+      !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
+
+            !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,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)
+              qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
+              cenv_vol(ind1,ind2)=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,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
+
+              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,ind2)
+                qincloud(ind1)=0.
+            ELSE               
+                qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
+                qincloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)
+            ENDIF 
+
+
+        ELSE ! mpc_bl_points>0
+
+            ! Treat boundary layer mixed phase clouds
+            
+            ! thermals
+            !=========
+
+            ! ice phase
+            !...........
+
+            deltazlev_mpc=dz(ind1,:)
+            temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
+            pres_mpc=pplay(ind1,:)
+            fraca_mpc=fraca(ind1,:)
+            snowf_mpc=snowflux(ind1,:)
+            qth_mpc=zqta(ind1,:)
+            flag_topthermals=.FALSE.
+            IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0))  THEN
+                flag_topthermals = .TRUE.
+            ENDIF
+
+            CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,flag_topthermals,temp_mpc,pres_mpc,qth_mpc,qlth(ind1,:),qith(ind1,:),deltazlev_mpc,snowf_mpc,fraca_mpc,qiceth)
+
+
+
+            ! We account for the effect of ice crystals in thermals on sthl
+            ! and on the width of the distribution
+
+            sthl=sthl*1./(1.+C_mpc*qiceth)  &
+                + (1.-1./(1.+C_mpc*qiceth)) * athl*(zqta(ind1,ind2)-(qsith+qiceth))  
+
+            sthi=sthi*1./(1.+C_mpc*qiceth)  &
+                + (1.-1./(1.+C_mpc*qiceth)) * athl*(zqta(ind1,ind2)-(qsith+qiceth))  
+
+           ! standard deviation of the water distribution in thermals
+            sth=sthl
+            senv=senvl
+            sigma2s=(sigma2s_factor*((MAX((sth-senv),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
+            deltasth=athl*vert_alpha_th*sigma2s
+      
+            ! Liquid phase
+            !.............
+            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)
+            qliqth=IntJ+IntI1+IntI2+IntI3
+
+            ! qimax calculation
+            xth1=-(sthi+deltasth)/(sqrt(2.)*sigma2s)
+            xth2=-(sthi-deltasth)/(sqrt(2.)*sigma2s)           
+            exp_xth1 = exp(-1.*xth1**2)
+            exp_xth2 = exp(-1.*xth2**2)
+            IntJ=0.5*sthi*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
+            IntJ_CF=0.5*(1.-1.*erf(xth2))
+            IntI1=(((sthi+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*(sthi+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
+            IntI1_CF=((sthi+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
+            IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
+            qimax=IntJ+IntI1+IntI2+IntI3
+            qimax=qimax-qliqth
+
+
+            ! Condensed water
+            ! Guarantee the consistency between qiceth and the subgrid scale PDF of total water
+            qlth(ind1,ind2)=MAX(0.,qliqth)
+            qith(ind1,ind2)=MAX(0.,MIN(qiceth,qimax))
+            qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
+
+            ! 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-qslth)
+            maxi=0. 
+            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)
+
+            ! 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*(sth+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
+            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,ind2)
+            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)
+                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  ! gaussian for environment only
+
+      
+        zqenv(ind1)=po(ind1)
+        Tbef=temp(ind1,ind2)
+
+        CALL CALC_QSAT_ECMWF(Tbef,0.,paprs(ind1,ind2),RTT,0,.false.,qsatbef,dqsatdt)
+        zqsatenv(ind1,ind2)=qsatbef
+
+        IF (Tbef .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*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)
+        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)=zqsatenv(ind1,ind2)
+          qincloud(ind1)=0.
+        ELSE 
+          qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
+          qincloud(ind1)=MAX(qctot(ind1,ind2)/ctot(ind1,ind2),0.)
+        ENDIF
+ 
+
+    ENDIF       ! From the separation (thermal/envrionnement) and (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       !loop on klon 
+
+RETURN
+
+
+END SUBROUTINE cloudth_mpc
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev,flag_topthermals,temp,pres,qth,qlth,qith,deltazlev,snowf,fraca,qi)
+
+! parameterization of ice for boundary
+! layer mixed-phase clouds assuming a stationary system
+!
+! Note that vapor deposition on ice crystals and riming of liquid droplets
+! depend on the ice number concentration Ni
+! One could assume that Ni depends on qi, e.g.,  Ni=beta*(qi*rho)**xi
+! and use values from Hong et al. 2004, MWR for instance
+! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962 
+! One could also think of a more complex expression of Ni;
+! function of qi, T, the concentration in aerosols or INP ..
+! Here we prefer fixing Ni to a tuning parameter
+! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard
+! in Mioche et al. 2017
+!
+!
+! References:
+!------------
+! This parameterization is thoroughly described in Vignon et al. 
+!
+! More specifically
+! for the Water vapor deposition process:
+!
+! Rotstayn, L. D., 1997: A physically based scheme for the treat-
+! ment of stratiform cloudfs and precipitation in large-scale
+! models. I: Description and evaluation of the microphysical
+! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282.
+!
+!  Morrison, H., and A. Gettelman, 2008: A new two-moment bulk
+!  stratiform cloud microphysics scheme in the NCAR Com-
+!  munity Atmosphere Model (CAM3). Part I: Description and
+!  numerical tests. J. Climate, 21, 3642–3659
+!
+! for the Riming process:
+!
+! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro-
+! scale structure and organization of clouds and precipitation in
+! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’
+! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206
+!
+! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit
+! forecasts of winter precipitation using an improved bulk
+! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit
+! forecasts of winter precipitation using an improved bulk
+! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542
+!
+! For the formation of clouds by thermals:
+!
+! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of
+! the Atmospheric Sciences, 65, 407–425.
+!
+! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a
+! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3
+!
+!
+!
+! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr
+!=============================================================================
+
+    USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF
+    USE ioipsl_getin_p_mod, ONLY : getin_p
+    USE phys_state_var_mod, ONLY : fm_therm, detr_therm, entr_therm
+
+    IMPLICIT none
+
+
+    INCLUDE "YOMCST.h"
+    INCLUDE "nuage.h"
+
+    INTEGER, INTENT(IN) :: ind1,ind2, klev ! horizontal and vertical indices and dimensions
+    LOGICAL, INTENT(IN) :: flag_topthermals ! uppermost layer of thermals ?
+    REAL,  DIMENSION(klev), INTENT(IN) :: temp       ! temperature [K] within thermals
+    REAL,  DIMENSION(klev), INTENT(IN) :: pres       ! pressure [Pa]
+    REAL,  DIMENSION(klev), INTENT(IN) :: qth        ! mean specific water content in thermals [kg/kg]
+    REAL,  DIMENSION(klev), INTENT(IN) :: qlth       ! condensed liquid water in thermals, approximated value [kg/kg]
+    REAL,  DIMENSION(klev), INTENT(IN) :: qith       ! condensed ice water , thermals [kg/kg]
+    REAL,  DIMENSION(klev), INTENT(IN) :: deltazlev  ! layer thickness [m]
+    REAL,  DIMENSION(klev+1), INTENT(IN) :: snowf      ! snow flux at the upper inferface
+    REAL,  DIMENSION(klev+1), INTENT(IN) :: fraca      ! fraction of the mesh covered by thermals
+
+    REAL,  INTENT(OUT) :: qi        ! ice cloud specific content [kg/kg]
+
+
+    REAL, SAVE    :: Ni,  C_cap, Ei, d_top
+    !$OMP THREADPRIVATE(Ni, C_cap,Ei, d_top)
+    LOGICAL, SAVE :: firstcall = .TRUE.
+    !$OMP THREADPRIVATE(firstcall)
+
+
+    INTEGER ind2p1,ind2p2
+    REAL rho(klev)
+    REAL unsurtaudet, unsurtaustardep, unsurtaurim
+    REAL qsl, qsi, dqs, AA, BB, Ka, Dv, rhoi
+    REAL  p0, t0
+    REAL alpha, flux_term
+    REAL det_term, precip_term, rim_term, dep_term
+
+  
+    IF (firstcall) THEN
+        Ni=2.0e3
+        CALL getin_p('Ni', Ni)
+        WRITE(*,*) 'Ni = ', Ni
+
+        Ei=0.5
+        CALL getin_p('Ei', Ei)
+        WRITE(*,*) 'Ei = ', Ei
+
+        C_cap=0.5
+        CALL getin_p('C_cap', C_cap)
+        WRITE(*,*) 'C_cap = ', C_cap
+
+        d_top=0.8
+        CALL getin_p('d_top', d_top)
+        WRITE(*,*) 'd_top = ', d_top
+
+
+        firstcall=.FALSE.
+    ENDIF
+
+
+    ind2p1=ind2+1
+    ind2p2=ind2+2
+
+    ! Liquid water content:
+    !=====================
+    ! the liquid water content is not calculated in this routine
+
+    ! Ice water content
+    ! ==================
+
+    rho=pres/temp/RD  ! air density kg/m3
+
+    Ka=2.4e-2      ! thermal conductivity of the air, SI
+    p0=101325.0    ! ref pressure
+    T0=273.15      ! ref temp
+    rhoi=500.0     ! cloud ice density following Reisner et al. 1998
+    alpha=700.     ! fallvelocity param
+
+
+    IF (flag_topthermals) THEN ! uppermost thermals level, solve a third order polynomial with Cardan's method
+
+    Dv=0.0001*0.211*(p0/pres(ind2))*((temp(ind2)/T0)**1.94) ! water vapor diffusivity in air, SI
+
+    ! Detrainment term:
+    unsurtaudet=detr_therm(ind1,ind2)/rho(ind2)/deltazlev(ind2)
+ 
+    ! vertical flux
+    
+    flux_term=d_top*fm_therm(ind1,ind2)/deltazlev(ind2)
+
+    ! Deposition term
+    CALL CALC_QSAT_ECMWF(temp(ind2),0.,pres(ind2),RTT,2,.false.,qsi,dqs)
+    CALL CALC_QSAT_ECMWF(temp(ind2),0.,pres(ind2),RTT,1,.false.,qsl,dqs)
+    AA=RLSTT/Ka/temp(ind2)*(RLSTT/RV/temp(ind2)-1.)
+    BB=1./(rho(ind2)*Dv*qsi)
+    unsurtaustardep=C_cap*(Ni**0.66)*(qsl-qsi)/qsi*4.*RPI/(AA+BB)*(6.*rho(ind2)/rhoi/RPI/Gamma(4.))**(0.33)
+
+    ! Riming term  neglected at this level
+    !unsurtaurim=rho(ind2)*alpha*3./rhoi/2.*Ei*qlth(ind2)*((p0/pres(ind2))**0.4)
+
+    qi=rho(ind2)*unsurtaustardep/MAX((rho(ind2)*unsurtaudet-flux_term),1E-12)
+    qi=MAX(qi,0.)**(3./2.)
+
+    ELSE ! other levels, estimate qi(k) from variables at k+1 and k+2
+
+    Dv=0.0001*0.211*(p0/pres(ind2p1))*((temp(ind2p1)/T0)**1.94) ! water vapor diffusivity in air, SI
+
+    ! Detrainment term:
+
+    unsurtaudet=detr_therm(ind1,ind2p1)/rho(ind2p1)/deltazlev(ind2p1)
+    det_term=-unsurtaudet*qith(ind2p1)*rho(ind2p1)
+    
+    
+    ! Deposition term
+
+    CALL CALC_QSAT_ECMWF(temp(ind2p1),0.,pres(ind2p1),RTT,2,.false.,qsi,dqs)
+    CALL CALC_QSAT_ECMWF(temp(ind2p1),0.,pres(ind2p1),RTT,1,.false.,qsl,dqs)
+    AA=RLSTT/Ka/temp(ind2p1)*(RLSTT/RV/temp(ind2p1)-1.)
+    BB=1./(rho(ind2p1)*Dv*qsi)
+    unsurtaustardep=C_cap*(Ni**0.66)*(qsl-qsi)/qsi*4.*RPI/(AA+BB)*(6.*rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33)
+    dep_term=rho(ind2p1)*(qith(ind2p1)**0.33)*unsurtaustardep
+ 
+    ! Riming term
+
+    unsurtaurim=rho(ind2p1)*alpha*3./rhoi/2.*Ei*qlth(ind2p1)*((p0/pres(ind2p1))**0.4)
+    rim_term=rho(ind2p1)*qith(ind2p1)*unsurtaurim
+
+    ! Precip term
+
+    !precip_term=-1./deltazlev(ind2p1)*(fraca(ind2p2)*snowf(ind2p2)-fraca(ind2p1)*snowf(ind2p1))
+    ! We assume that there is no solid precipitation outside thermals (so no multiplication by fraca)
+    precip_term=-1./deltazlev(ind2p1)*(snowf(ind2p2)-snowf(ind2p1))
+
+    ! Calculation in a top-to-bottom loop
+
+    IF (fm_therm(ind1,ind2p1) .GT. 0.) THEN
+       qi= 1./fm_therm(ind1,ind2p1)* &
+          (deltazlev(ind2p1)*(-rim_term-dep_term-det_term-precip_term) + &
+          fm_therm(ind1,ind2p2)*(qith(ind2p1)))
+    ELSE
+       qi=0.
+    ENDIF
+
+    ENDIF ! flag_topthermals
+  
+    qi=MAX(0.,qi)
+
+    RETURN
+
+END SUBROUTINE ICE_MPC_BL_CLOUDS
+
+
+
+
 END MODULE cloudth_mod
 
Index: LMDZ6/trunk/libf/phylmd/conf_phys_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 3998)
+++ LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 3999)
@@ -182,7 +182,9 @@
     INTEGER,SAVE :: iflag_cloudth_vert_omp
     INTEGER,SAVE :: iflag_rain_incloud_vol_omp
+    INTEGER,SAVE :: iflag_vice_omp
     REAL,SAVE :: rad_froid_omp, rad_chau1_omp, rad_chau2_omp
     REAL,SAVE :: t_glace_min_omp, t_glace_max_omp
     REAL,SAVE :: exposant_glace_omp
+    INTEGER,SAVE :: iflag_gammasat_omp, iflag_mpc_bl_omp
     REAL,SAVE :: rei_min_omp, rei_max_omp
     INTEGER,SAVE :: iflag_sic_omp, iflag_inertie_omp
@@ -251,4 +253,6 @@
     LOGICAL, SAVE :: adjust_tropopause_omp
     LOGICAL, SAVE :: ok_daily_climoz_omp
+    LOGICAL, SAVE :: ok_new_lscp_omp
+    LOGICAL, SAVE :: ok_icefra_lscp_omp
 
     INTEGER, INTENT(OUT):: read_climoz ! read ozone climatology, OpenMP shared
@@ -1433,4 +1437,25 @@
 
     !
+    !Config Key  = iflag_gammasat
+    !Config Desc =  
+    !Config Def  = 0
+    !Config Help = 
+    !
+    iflag_gammasat_omp=0
+    CALL getin('iflag_gammasat',iflag_gammasat_omp)
+
+
+    !
+    !Config Key  = iflag_mpc_bl
+    !Config Desc =
+    !Config Def  = 0
+    !Config Help =
+    !
+    iflag_mpc_bl_omp=0
+    CALL getin('iflag_mpc_bl',iflag_mpc_bl_omp)
+
+
+
+    !
     !Config Key  = iflag_t_glace
     !Config Desc =  
@@ -1458,4 +1483,15 @@
     iflag_rain_incloud_vol_omp = 0
     CALL getin('iflag_rain_incloud_vol',iflag_rain_incloud_vol_omp)
+
+    !
+    !Config Key  = iflag_vice
+    !Config Desc =  
+    !Config Def  = 0
+    !Config Help = 
+    !
+    iflag_vice_omp = 0
+    CALL getin('iflag_vice',iflag_vice_omp)
+
+
 
     !
@@ -2313,5 +2349,26 @@
     !Config Help = .FALSE. ensure much fewer (no calendar dependency)
     !  and lighter monthly climoz files, inetrpolated in time at gcm run time.
-    !
+    
+    ok_new_lscp_omp = .FALSE.
+    CALL getin('ok_new_lscp', ok_new_lscp_omp)
+    !
+    !Config Key  = ok_new_lscp_omp
+    !Config Desc = new cloud scheme ith ice and mixed phase (Etienne and JB)
+    !Config Def  = .FALSE.
+    !Config Help = ...
+
+
+
+    ok_icefra_lscp_omp = .FALSE.
+    CALL getin('ok_icefra_lscp', ok_icefra_lscp_omp)
+    !
+    !Config Key  = ok_icefra_lscp_omp
+    !Config Desc = ice fraction in radiation from lscp
+    !Config Def  = .FALSE.
+    !Config Help = ...
+
+
+
+
     ecrit_LES_omp = 1./8.
     CALL getin('ecrit_LES', ecrit_LES_omp)
@@ -2427,7 +2484,10 @@
     t_glace_max = t_glace_max_omp
     exposant_glace = exposant_glace_omp
+    iflag_gammasat=iflag_gammasat_omp
+    iflag_mpc_bl=iflag_mpc_bl_omp
     iflag_t_glace = iflag_t_glace_omp
     iflag_cloudth_vert=iflag_cloudth_vert_omp
     iflag_rain_incloud_vol=iflag_rain_incloud_vol_omp
+    iflag_vice=iflag_vice_omp
     iflag_ice_thermo = iflag_ice_thermo_omp
     rei_min = rei_min_omp
@@ -2631,4 +2691,6 @@
     carbon_cycle_rad = carbon_cycle_rad_omp
     level_coupling_esm = level_coupling_esm_omp
+    ok_new_lscp = ok_new_lscp_omp
+    ok_icefra_lscp=ok_icefra_lscp_omp
     read_fco2_ocean_cor = read_fco2_ocean_cor_omp
     var_fco2_ocean_cor = var_fco2_ocean_cor_omp
@@ -2849,7 +2911,10 @@
     WRITE(lunout,*) ' t_glace_max = ',t_glace_max
     WRITE(lunout,*) ' exposant_glace = ',exposant_glace
+    WRITE(lunout,*) ' iflag_gammasat = ',iflag_gammasat
+    WRITE(lunout,*) ' iflag_mpc_bl = ',iflag_mpc_bl
     WRITE(lunout,*) ' iflag_t_glace = ',iflag_t_glace
     WRITE(lunout,*) ' iflag_cloudth_vert = ',iflag_cloudth_vert
     WRITE(lunout,*) ' iflag_rain_incloud_vol = ',iflag_rain_incloud_vol
+    WRITE(lunout,*) ' iflag_vice = ',iflag_vice
     WRITE(lunout,*) ' iflag_ice_thermo = ',iflag_ice_thermo
     WRITE(lunout,*) ' rei_min = ',rei_min
@@ -2964,4 +3029,6 @@
     WRITE(lunout,*) ' adjust_tropopause = ', adjust_tropopause
     WRITE(lunout,*) ' ok_daily_climoz = ',ok_daily_climoz
+    WRITE(lunout,*) ' ok_new_lscp = ', ok_new_lscp
+    WRITE(lunout,*) ' ok_icefra_lscp = ', ok_icefra_lscp
     WRITE(lunout,*) ' read_climoz = ', read_climoz
     WRITE(lunout,*) ' carbon_cycle_tr = ', carbon_cycle_tr
Index: LMDZ6/trunk/libf/phylmd/lscp_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 3999)
+++ LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 3999)
@@ -0,0 +1,1269 @@
+MODULE LSCP_mod
+
+IMPLICIT NONE
+
+CONTAINS
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+SUBROUTINE LSCP(dtime,paprs,pplay,t,q,ptconv,ratqs,     &
+     d_t, d_q, d_ql, d_qi, rneb, radliq, radicefrac,    &
+     rain, snow,                                        &
+     pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
+     frac_impa, frac_nucl, beta,                        &
+     prfl, psfl, rhcl, zqta, fraca,                     &
+     ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
+     iflag_ice_thermo)
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
+!          A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
+!--------------------------------------------------------------------------------
+! Date: 01/2021
+!--------------------------------------------------------------------------------
+! Aim: Large Scale Clouds and Precipitation (LSCP)
+!      
+! This code is a new version of the fisrtilp.F90 routine, which is itself a 
+! fusion of 'first' (superrsaturation physics, P. LeVan K. Laval)
+! and 'ilp' (il pleut, L. Li)
+!
+! Compared to fisrtilp, LSCP 
+! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
+! -> consider always precipitation thermalisation (fl_cor_ebil>0)
+! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T) 
+! -> option "oldbug" by JYG has been removed
+! -> iflag_t_glace >0 always
+! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
+! -> rectangular distribution from L. Li no longer available
+! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt) 
+!--------------------------------------------------------------------------------
+! References: 
+!
+! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
+! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
+! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
+! -------------------------------------------------------------------------------
+! Code structure:
+!
+! P0>     Thermalization of the  precipitation coming from the overlying layer
+! P1>     Evaporation of the precipitation (falling from the k+1 level)
+! P2>     Cloud formation (at the k level) 
+! P2.A.0> Cloud properties calculation from a rectangular pdf
+! P2.A.1> With the new PDFs, calculation of cloud properties using the inital
+!         values of T and Q
+! P2.A.2> Coupling between condensed water and temperature
+! P2.A.3> Calculation of final quantities associated with cloud formation
+! P2.B>   'All or nothing' cloud
+! P2.C>   Release of Latent heat after cloud formation
+! P3>     Autoconversion to precipitation (k-level)
+! P4>     Wet scavenging
+!------------------------------------------------------------------------------
+! Some preliminary comments (JBM) :
+!
+! The cloud water that the radiation scheme sees is not the same that the cloud
+! water used in the physics and the dynamics 
+! 
+! During the autoconversion to precipitation (P3 step), radliq (cloud water used
+! by the radiation scheme) is calculated as an average of the water that remains
+! in the cloud during the precipitation and not the water remaining at the end
+! of the time step. The latter is used in the rest of the physics and advected
+! by the dynamics. 
+!
+! In summary:
+!
+! Radiation:
+! xflwc(newmicro)+xfiwc(newmicro) =
+!   cldliq(physiq)=radliq(fisrt)=lwcon(nc)+iwcon(nc)
+!
+! Physics/Dynamics:
+! ql_seri(physiq)+qs_seri(physiq)=ocond(nc)
+!
+! Notetheless, be aware of:
+!
+! radliq(fisrt) .NE. ocond(nc)
+! i.e.:
+! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
+! (which is not trivial)
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  
+USE dimphy
+USE print_control_mod, ONLY: prt_level, lunout
+USE cloudth_mod 
+USE ioipsl_getin_p_mod, ONLY : getin_p
+USE phys_local_var_mod, ONLY: ql_seri,qs_seri
+USE phys_local_var_mod, ONLY: rneblsvol
+USE lscp_tools_mod, ONLY : CALC_QSAT_ECMWF, ICEFRAC_LSCP, CALC_GAMMASAT, FALLICE_VELOCITY
+
+
+IMPLICIT NONE
+
+
+
+!===============================================================================
+! VARIABLE DECLARATION
+!===============================================================================
+
+include "YOMCST.h"
+include "YOETHF.h"
+include "FCTTRE.h"
+include "fisrtilp.h"
+include "nuage.h" 
+
+  
+! INPUT VARIABLES:
+!-----------------
+
+  REAL,                            INTENT(IN)   :: dtime           ! time step [s]
+  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: t               ! temperature (K)
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: q               ! specific humidity [kg/kg]
+  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
+  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
+                                                                   ! CR: if iflag_ice_thermo=2, only convection is active    
+  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
+
+  
+  !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(INOUT)   :: zthl            ! liquid potential temperature [K]
+
+  ! INPUT/OUTPUT variables
+  !------------------------
+
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratqs            ! function of pressure that sets the large-scale
+                                                                    ! cloud PDF (sigma=ratqs*qt)
+  
+  ! OUTPUT variables
+  !-----------------
+
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t              ! temperature increment [K]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q              ! specific humidity increment [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql             ! liquid water increment [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi             ! cloud ice mass increment [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb             ! cloud fraction [-]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radliq           ! condensed water used in the radiation scheme [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac       ! ice fraction of condensed water for radiation scheme
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl             ! clear-sky relative humidity [-]
+  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! large-scale rainfall [kg/s/m2] 
+  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! large-scale snowfall [kg/s/m2]
+  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
+  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
+
+ 
+  ! cofficients of scavenging fraction (for off-line)
+ 
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_nucl       ! scavenging fraction due tu nucleation [-]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_1nucl      ! scavenging fraction due tu nucleation with a -1 factor [-]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_impa       ! scavening fraction due to impaction [-]
+  
+  ! fraction of aerosol scavenging through impaction and nucleation (for on-line)
+  
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa        ! scavenging fraction due tu impaction [-]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl        ! scavenging fraction due tu nucleation [-]
+  
+
+
+
+  ! PROGRAM PARAMETERS:
+  !--------------------
+  
+  
+  REAL, SAVE :: seuil_neb=0.001                                     ! cloud fraction threshold: a cloud really exists when exceeded
+  !$OMP THREADPRIVATE(seuil_neb)
+
+  INTEGER, SAVE :: ninter=5                                         ! number of iterations to calculate autoconversion to precipitation
+  !$OMP THREADPRIVATE(ninter)
+
+  INTEGER,SAVE :: iflag_evap_prec=1                                 ! precipitation evap. flag. 0: nothing, 1: "old way", 2: Max cloud fraction above to calculate the max of reevaporation
+                                                                    ! 4: LTP'method i.e. evaporation in the clear-sky fraction of the mesh only
+  !$OMP THREADPRIVATE(iflag_evap_prec)
+                                                                    
+
+  REAL t_coup                                                       ! temperature threshold which determines the phase 
+  PARAMETER (t_coup=234.0)                                          ! for which the saturation vapor pressure is calculated
+
+  REAL DDT0                                                         ! iteration parameter 
+  PARAMETER (DDT0=.01)
+
+  REAL ztfondue                                                     ! parameter to calculate melting fraction of precipitation
+  PARAMETER (ztfondue=278.15)
+
+  REAL, SAVE    :: rain_int_min=0.001                               ! Minimum local rain intensity [mm/s] before the decrease in  associates precipitation fraction 
+  !$OMP THREADPRIVATE(rain_int_min)
+
+
+  
+  ! LOCAL VARIABLES:
+  !----------------
+
+
+  REAL qsl, qsi
+  REAL zct, zcl
+  REAL ctot(klon,klev)
+  REAL ctot_vol(klon,klev)
+  INTEGER mpc_bl_points(klon,klev)
+  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5  
+  REAL zdqsdT_raw(klon)
+  REAL gammasat(klon),dgammasatdt(klon)                                 ! coefficient to make cold condensation at the correct RH and derivative wrt T
+  REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
+  REAL cste
+  REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
+  REAL Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
+  REAL erf   
+  REAL qcloud(klon), icefrac_mpc(klon,klev), qincloud_mpc(klon)
+  REAL zrfl(klon), zrfln(klon), zqev, zqevt 
+  REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
+  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon)
+  REAL zoliqp(klon), zoliqi(klon)
+  REAL zt(klon)
+  REAL zdz(klon),zrho(klon),ztot, zrhol(klon)
+  REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon)
+  REAL zmelt,zpluie,zice
+  REAL dzfice(klon)
+  REAL zsolid,qtot,dqsl,dqsi
+  REAL smallestreal
+  !  Variables for Bergeron process
+  REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
+  REAL zqpreci(klon), zqprecl(klon)
+  ! Variables precipitation energy conservation
+  REAL zmqc(klon)
+  ! Variables for tracers
+  ! temporary: alpha parameter for scavenging
+  ! 4 possible scavenging processes
+  REAL a_tr_sca(4)
+  save a_tr_sca
+  !$OMP THREADPRIVATE(a_tr_sca)  
+  REAL zalpha_tr
+  REAL zfrac_lessi
+  REAL zprec_cond(klon)
+  REAL beta(klon,klev)                   ! conversion rate of condensed water
+  REAL zmair(klon), zcpair, zcpeau
+  REAL zlh_solid(klon), zm_solid         ! for liquid -> solid conversion
+  REAL zrflclr(klon), zrflcld(klon)                                                                    
+  REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon)                                                   
+  REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon) 
+  REAL ziflclr(klon), ziflcld(klon)                                                                    
+  REAL znebprecipclr(klon), znebprecipcld(klon)                                                        
+  REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon)                                              
+  REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon)
+  REAL velo(klon)
+  REAL qlmpc, qimpc, rnebmpc
+  REAL radliqi(klon,klev)
+
+  INTEGER i, k, n, kk
+  INTEGER n_i(klon), iter
+  INTEGER ncoreczq  
+  INTEGER,SAVE :: itap=0
+  !$OMP THREADPRIVATE(itap)
+
+
+  LOGICAL lognormale(klon)
+  LOGICAL convergence(klon)
+  LOGICAL appel1er
+  SAVE appel1er
+  !$OMP THREADPRIVATE(appel1er)
+
+
+  CHARACTER (len = 20) :: modname = 'lscp'
+  CHARACTER (len = 80) :: abort_message
+
+  
+  DATA appel1er /.TRUE./
+
+
+
+
+!===============================================================================
+! INITIALISATION
+!===============================================================================
+
+! Few initial checksS
+
+IF (iflag_t_glace.EQ.0) THEN
+    abort_message = 'lscp cannot be used if iflag_t_glace=0'
+    CALL abort_physic(modname,abort_message,1)
+ENDIF
+
+IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN
+    abort_message = 'lscp cannot be used without ice thermodynamics'
+    CALL abort_physic(modname,abort_message,1)
+ENDIF
+
+IF (.NOT.thermcep) THEN
+     abort_message = 'lscp cannot be used when thermcep=false'
+     CALL abort_physic(modname,abort_message,1)
+ENDIF
+
+IF (iflag_fisrtilp_qsat .LT. 0) THEN
+        abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
+     CALL abort_physic(modname,abort_message,1)
+ENDIF
+
+
+! Few initialisations
+
+itap=itap+1
+znebprecip(:)=0.
+ctot_vol(1:klon,1:klev)=0.0
+rneblsvol(1:klon,1:klev)=0.0
+smallestreal=1.e-9                                                                                  
+znebprecipclr(:)=0.    
+znebprecipcld(:)=0.                                                                            
+mpc_bl_points(:,:)=0
+
+IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
+
+IF (appel1er) THEN
+    CALL getin_p('ninter',ninter)
+    CALL getin_p('iflag_evap_prec',iflag_evap_prec)
+    CALL getin_p('seuil_neb',seuil_neb)
+    CALL getin_p('rain_int_min',rain_int_min) 
+    WRITE(lunout,*) 'lscp, ninter:', ninter
+    WRITE(lunout,*) 'lscp, iflag_evap_prec:', iflag_evap_prec
+    WRITE(lunout,*) 'lscp, seuil_neb:', seuil_neb
+    WRITE(lunout,*) 'lscp, rain_int_min:', rain_int_min
+
+    ! check for precipitation sub-time steps
+    IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
+        WRITE(lunout,*) 'lscp: it is not expected, see Z.X.Li', dtime
+        WRITE(lunout,*) 'I would prefer a 6 min sub-timestep'
+    ENDIF
+     
+    !AA Temporary initialisation
+    a_tr_sca(1) = -0.5
+    a_tr_sca(2) = -0.5
+    a_tr_sca(3) = -0.5
+    a_tr_sca(4) = -0.5
+
+    !AA Set scavenged fractions to 1
+    DO k = 1, klev
+        DO i = 1, klon
+            pfrac_nucl(i,k)=1.
+            pfrac_1nucl(i,k)=1.
+            pfrac_impa(i,k)=1.
+            beta(i,k)=0.       
+        ENDDO
+    ENDDO
+
+    appel1er = .FALSE.
+
+ENDIF     !(appel1er)  
+
+! AA for 'safety' reasons
+zalpha_tr   = 0.
+zfrac_lessi = 0.
+
+  
+! Initialisation of output variables (JYG):
+prfl(:,:) = 0.0
+psfl(:,:) = 0.0
+
+d_t(:,:) = 0.0
+d_q(:,:) = 0.0
+d_ql(:,:) = 0.0
+d_qi(:,:) = 0.0
+rneb(:,:) = 0.0
+radliq(:,:) = 0.0
+radicefrac(:,:) = 0.0
+frac_nucl(:,:) = 1. 
+frac_impa(:,:) = 1. 
+
+rain(:) = 0.0
+snow(:) = 0.0
+zoliq(:)=0.
+zrfl(:) = 0.0
+zifl(:) = 0.0
+zneb(:) = seuil_neb
+zrflclr(:) = 0.0                                                                            
+ziflclr(:) = 0.0                                                                            
+zrflcld(:) = 0.0                                                                               
+ziflcld(:) = 0.0                                                                                  
+tot_zneb(:) = 0.0                                                                                 
+tot_znebn(:) = 0.0                                                                               
+d_tot_zneb(:) = 0.0   
+
+
+!===============================================================================
+!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
+!===============================================================================
+
+ncoreczq=0
+
+DO k = klev, 1, -1
+
+
+    ! Initialisation temperature and specific humidity
+    DO i = 1, klon
+        zt(i)=t(i,k)
+        zq(i)=q(i,k)
+    ENDDO
+    
+    ! --------------------------------------------------------------------
+    ! P0> Thermalization of precipitation falling from the overlying layer
+    ! --------------------------------------------------------------------
+    ! Computes air temperature variation due to latent heat transported by
+    ! precipitation. Precipitation is then thermalized with the air in the
+    ! layer. 
+    ! The precipitation should remain thermalized throughout the different
+    ! thermodynamical transformations. The corresponding water mass should
+    ! be added when calculating the layer's enthalpy change with 
+    ! temperature
+    ! ---------------------------------------------------------------------
+    
+    IF (k.LE.klevm1) THEN 
+
+        DO i = 1, klon
+ 
+            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
+            ! no condensed water so cp=cp(vapor+dry air)
+            ! RVTMP2=rcpv/rcpd-1
+            zcpair=RCPD*(1.0+RVTMP2*zq(i)) 
+            zcpeau=RCPD*RVTMP2             
+                                          
+            ! zmqc: precipitation mass that has to be thermalized with 
+            ! layer's air so that precipitation at the ground has the
+            ! same temperature as the lowermost layer
+            zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
+            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
+            zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
+             / (zcpair + zmqc(i)*zcpeau)
+           
+        ENDDO
+ 
+    ELSE  
+ 
+        DO i = 1, klon 
+            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG 
+            zmqc(i) = 0.
+        ENDDO
+ 
+    ENDIF
+
+    ! --------------------------------------------------------------------
+    ! P1> Precipitation evaporation/sublimation
+    ! --------------------------------------------------------------------
+    ! A part of the precipitation coming from above is evaporated/sublimated.
+    !
+    ! --------------------------------------------------------------------
+
+
+    IF (iflag_evap_prec.GE.1) THEN
+
+
+        DO i = 1, klon
+     
+            ! if precipitation
+            IF (zrfl(i)+zifl(i).GT.0.) THEN
+                    
+                CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
+
+            ! LTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4).
+                IF (iflag_evap_prec.EQ.4) THEN                                                          
+                    zrfl(i) = zrflclr(i)                                                          
+                    zifl(i) = ziflclr(i)                                                           
+                ENDIF                      
+
+                IF (iflag_evap_prec.EQ.1) THEN
+                    znebprecip(i)=zneb(i)
+                ELSE
+                    znebprecip(i)=MAX(zneb(i),znebprecip(i))
+                ENDIF
+     
+                                                                                                                                                                                                                  
+                IF (iflag_evap_prec.EQ.4) THEN                                                                 
+                ! Max evaporation not to saturate the whole mesh                                             
+                    zqev0 = MAX(0.0, zqs(i)-zq(i))                                                              
+                ELSE       
+                ! Max evap not to saturate the fraction below the cloud                                                   
+                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))                                                                                  
+                ENDIF   
+
+                ! A. JAM
+                ! We consider separately qsatice and qstal
+                ! qsat wrt liquid phase according to thermcep
+                CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,1,.false.,qsl,dqsl)
+
+         
+                ! Evaporation of liquid precipitation coming from above
+                ! dP/dz=beta*(1-q/qsat)*sqrt(P)
+                ! formula from Sundquist 1988, Klemp & Wilhemson 1978
+                ! LTP: evaporation only in the clear sky part (iflag_evap_prec=4)
+
+                IF (iflag_evap_prec.EQ.3) THEN
+                    zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl) &
+                    *SQRT(zrfl(i)/max(1.e-4,znebprecip(i)))        &
+                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+                ELSE IF (iflag_evap_prec.EQ.4) THEN
+                     zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl) &
+                    *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
+                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+                ELSE
+                    zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
+                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+                ENDIF
+
+
+                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
+                *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 
+         
+                ! qsat wrt ice according to thermcep
+                CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,2,.false.,qsi,dqsi)
+
+                ! sublimation of the solid precipitation coming from above
+                IF (iflag_evap_prec.EQ.3) THEN
+                    zqevti = znebprecip(i)*coef_eva*(1.0-zq(i)/qsi) &
+                    *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
+                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+                ELSE IF (iflag_evap_prec.EQ.4) THEN
+                     zqevti = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsi) &
+                    *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
+                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+                ELSE
+                    zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
+                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+                ENDIF
+            
+                zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
+                *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
+
+                ! A. JAM
+                ! Evaporation limit: we ensure that the layer's fraction below
+                ! the cloud or the whole mesh (depending on iflag_evap_prec)
+                ! does not reach saturation. In this case, we 
+                ! redistribute zqev0 conserving the ratio liquid/ice 
+   
+                IF (zqevt+zqevti.GT.zqev0) THEN
+                    zqev=zqev0*zqevt/(zqevt+zqevti)
+                    zqevi=zqev0*zqevti/(zqevt+zqevti)
+                ELSE
+                    zqev=zqevt
+                    zqevi=zqevti
+                ENDIF
+
+
+                ! EV: agrees with JL's comments below so correct and comment
+                ! JLD: I do not understand the lines below. We distribute the liquid 
+                ! and solid parts of the precipitation even if the layer does not get
+                ! saturated. To my opinion, we should all replace with:
+                ! zqev=zqevt
+                ! zqevi=zqevti
+                !    IF (zqevt+zqevti.GT.0.) THEN
+                !        zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
+                !        zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
+                !    ELSE
+                !        zqev=0.
+                !        zqevi=0.
+                !    ENDIF
+                ! ENDIF
+
+                ! New solid and liquid precipitation fluxes
+                zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))  &
+                         /RG/dtime)
+                zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
+                         /RG/dtime)
+
+                ! vapor, temperature, precip fluxes update
+                zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
+                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+                ! precip thermalization
+                zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
+                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+                zt(i) = zt(i) + (zrfln(i)-zrfl(i))        &
+                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
+                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
+                + (zifln(i)-zifl(i))                      &
+                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
+                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+
+                ! New values of liquid and solid precipitation 
+                zrfl(i) = zrfln(i)
+                zifl(i) = zifln(i)
+
+                IF (iflag_evap_prec.EQ.4) THEN                                                                 
+                    zrflclr(i) = zrfl(i)                                                                  
+                    ziflclr(i) = zifl(i)    
+                    IF(zrflclr(i) + ziflclr(i).LE.0) THEN                                                
+                        znebprecipclr(i) = 0.                                                          
+                    ENDIF                                                                                
+                    zrfl(i) = zrflclr(i) + zrflcld(i)                                                      
+                    zifl(i) = ziflclr(i) + ziflcld(i)                                                      
+                ENDIF      
+
+
+                ! CR Be careful: ice melting has been moved
+                zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! JYG
+                ! precip fraction that is melted
+                zmelt = MIN(MAX(zmelt,0.),1.)
+                ! Ice melting
+                IF (iflag_evap_prec.EQ.4) THEN                                                              
+                    zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)                                            
+                    zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)                                              
+                    zrfl(i)=zrflclr(i)+zrflcld(i)                                                       
+                ELSE                                                                                        
+                    zrfl(i)=zrfl(i)+zmelt*zifl(i)                                                      
+                ENDIF                  
+
+                ! Latent heat of melting with precipitation thermalization
+                zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
+                *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+                
+                IF (iflag_evap_prec.EQ.4) THEN                                                             
+                    ziflclr(i)=ziflclr(i)*(1.-zmelt)                                                    
+                    ziflcld(i)=ziflcld(i)*(1.-zmelt)                                                    
+                    zifl(i)=ziflclr(i)+ziflcld(i)                                                       
+                ELSE                                                                                      
+                    zifl(i)=zifl(i)*(1.-zmelt)                                                          
+                ENDIF
+
+            ELSE
+                ! if no precip, we reinitialize the cloud fraction used for the precip to 0 
+                znebprecip(i)=0.
+
+            ENDIF ! (zrfl(i)+zifl(i).GT.0.)
+
+        ENDDO
+         
+    ENDIF ! (iflag_evap_prec>=1)
+
+    ! --------------------------------------------------------------------
+    ! End precip evaporation
+    ! --------------------------------------------------------------------
+
+
+
+    
+    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
+    !-------------------------------------------------------
+       
+        DO i = 1, klon
+            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
+            qtot=zq(i)+zmqc(i)
+            CALL CALC_QSAT_ECMWF(zt(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
+            zdqsdT_raw(i) = zdqs(i)*  &
+            & RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
+         
+            IF (zq(i) .LT. 1.e-15) THEN
+                ncoreczq=ncoreczq+1
+                zq(i)=1.e-15
+            ENDIF
+
+        ENDDO
+    
+
+    ! --------------------------------------------------------------------
+    ! P2> Cloud formation
+    !---------------------------------------------------------------------
+    !
+    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach 
+    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
+    ! is prescribed and depends on large scale variables and boundary layer
+    ! properties)
+    ! The decrease in condensed part due tu latent heating is taken into
+    ! account
+    ! -------------------------------------------------------------------
+
+           
+        ! P2.1> With the PDFs (log-normal, bigaussian) 
+        ! cloud propertues calculation with the initial values of t and q
+        ! ----------------------------------------------------------------
+
+        ! initialise gammasat and qincloud_mpc
+        gammasat(:)=1.
+        qincloud_mpc(:)=0.
+
+
+
+        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, &
+                         ratqs,zqs,t)
+
+                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, &
+                         ratqs,zqs,t)
+
+
+                !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, &
+                         ratqs,zqs,t)
+
+                ENDIF
+
+         ELSE 
+                ! cloudth_v3 + cold microphysical considerations 
+                ! + stationary mixed-phase cloud model
+
+                    CALL cloudth_mpc(klon,klev,k,mpc_bl_points,            &
+                         t,ztv,zq,zqta,fraca, zpspsk,                      &
+                         paprs,pplay,ztla,zthl,ratqs,zqs,psfl,             &
+                         qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol)
+
+         ENDIF ! iflag_mpc_bl
+
+                DO i=1,klon
+                    rneb(i,k)=ctot(i,k)
+                    rneblsvol(i,k)=ctot_vol(i,k)
+                    zqn(i)=qcloud(i)
+                ENDDO
+
+        ENDIF
+
+        IF (iflag_cld_th .LE. 4) THEN
+            
+                ! lognormal 
+            lognormale = .TRUE.
+
+        ELSEIF (iflag_cld_th .GE. 6) THEN
+            
+                ! lognormal distribution when no thermals 
+            lognormale = fraca(:,k) < 1e-10
+
+        ELSE
+                ! When iflag_cld_th=5, we always assume  
+                ! bi-gaussian distribution
+            lognormale = .FALSE.
+        
+        ENDIF
+
+
+        DT(:) = 0.
+        n_i(:)=0
+        Tbef(:)=zt(:)
+        qlbef(:)=0.
+
+        ! Treatment of non-boundary layer clouds (lognormale)
+        ! condensation with qsat(T) variation (adaptation)
+        ! Iterative Loop to converge towards qsat
+        
+        
+        DO iter=1,iflag_fisrtilp_qsat+1
+                
+                DO i=1,klon
+
+                    ! convergence = .true. since convergence is not satisfied
+                    convergence(i)=ABS(DT(i)).GT.DDT0
+                   
+                    IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
+                          
+                        ! if not convergence:
+
+                        ! P2.2.1> cloud fraction and condensed water mass calculation
+                        ! Calculated variables: 
+                        ! rneb : cloud fraction
+                        ! zqn : total water within the cloud
+                        ! zcond: mean condensed water within the mesh
+                        ! rhcl: clear-sky relative humidity
+                        !---------------------------------------------------------------
+
+                        ! new temperature:
+                        Tbef(i)=Tbef(i)+DT(i)
+
+
+                        ! Rneb, qzn and zcond for lognormal PDFs
+                        qtot=zq(i)+zmqc(i)
+                        CALL CALC_QSAT_ECMWF(Tbef(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
+                        CALL CALC_GAMMASAT(Tbef(i),qtot,pplay(i,k),gammasat(i),dgammasatdt(i))
+                        
+                        ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
+                        zdqs(i) = gammasat(i)*zdqs(i)+zqs(i)*dgammasatdt(i) 
+
+                        zpdf_sig(i)=ratqs(i,k)*zq(i)
+                        zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
+                        zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i)))
+                        zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
+                        zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
+                        zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
+                        zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i))
+                        zpdf_e1(i)=1.-erf(zpdf_e1(i))
+                        zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
+                        zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i))
+                        zpdf_e2(i)=1.-erf(zpdf_e2(i))
+             
+                        IF (zpdf_e1(i).LT.1.e-10) THEN
+                            rneb(i,k)=0.
+                            zqn(i)=gammasat(i)*zqs(i)
+                        ELSE
+                            rneb(i,k)=0.5*zpdf_e1(i)
+                            zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
+                        ENDIF
+
+                        ! If vertical heterogeneity, change fraction by volume as well
+                        IF (iflag_cloudth_vert.GE.3) THEN
+                            ctot_vol(i,k)=rneb(i,k)
+                            rneblsvol(i,k)=ctot_vol(i,k)
+                        ENDIF
+
+
+            ! P2.2.2> Approximative calculation of temperature variation DT 
+            !           due to condensation.
+            ! Calculated variables: 
+            ! dT : temperature change due to condensation
+            ! ---------------------------------------------------------------
+
+                       ! EV: calculation of icefrac in one sole function
+                        CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
+
+                
+                        IF (zfice(i).LT.1) THEN
+                            cste=RLVTT
+                        ELSE
+                            cste=RLSTT
+                        ENDIF
+
+                        qlbef(i)=max(0.,zqn(i)-zqs(i))
+                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
+                        &   +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
+                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
+                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
+                        &     *qlbef(i)*dzfice(i)
+                        DT(i)=num/denom
+                        n_i(i)=n_i(i)+1
+
+                    ENDIF ! end convergence
+
+                ENDDO     ! end loop on i
+
+        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
+
+
+
+        ! P2.3> Final quantities calculation 
+        ! Calculated variables: 
+        !   rneb : cloud fraction
+        !   zcond: mean condensed water in the mesh
+        !   zqn  : mean water vapor in the mesh
+        !   zt   : temperature
+        !   rhcl : clear-sky relative humidity
+        ! ----------------------------------------------------------------
+
+
+
+        ! Water vapor update, Phase determination and subsequent latent heat exchange
+
+        ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs)
+        CALL icefrac_lscp(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))
+
+        DO i=1, klon
+
+
+            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
+                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
+                ! water vapor update and partition function if thermals
+                zq(i) = zq(i) - zcond(i)        
+                zfice(i)=icefrac_mpc(i,k)
+
+            ELSE
+
+                ! Checks on rneb, rhcl and zqn
+                IF (rneb(i,k) .LE. 0.0) THEN
+                    zqn(i) = 0.0
+                    rneb(i,k) = 0.0
+                    zcond(i) = 0.0
+                    rhcl(i,k)=zq(i)/zqs(i)
+                ELSE IF (rneb(i,k) .GE. 1.0) THEN
+                    zqn(i) = zq(i)
+                    rneb(i,k) = 1.0
+                    zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i))
+                    rhcl(i,k)=1.0
+                ELSE
+                    zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i))*rneb(i,k)
+                    ! following line is very strange and probably wrong:
+                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
+                ENDIF
+
+
+                ! water vapor update 
+                zq(i) = zq(i) - zcond(i)        
+            
+            ENDIF
+
+            ! temperature update due to phase change
+            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) & 
+            &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
+                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
+        ENDDO
+
+
+        ! If vertical heterogeneity, change volume fraction 
+        IF (iflag_cloudth_vert .GE. 3) THEN
+          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
+          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
+        ENDIF
+
+    ! ----------------------------------------------------------------
+    ! P3> Precipitation formation and calculation of condensed
+    !     water seen by the radiation scheme
+    ! ----------------------------------------------------------------
+    
+    ! LTP:
+    IF (iflag_evap_prec==4) THEN
+
+        !Partitionning between precipitation coming from clouds and that coming from CS
+
+        !0) Calculate tot_zneb, total cloud fraction above the cloud
+        !assuming a maximum-random overlap (voir Jakob and Klein, 2000)
+        
+        DO i=1, klon
+                tot_znebn(i) = 1 - (1-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
+                        /(1-min(zneb(i),1-smallestreal))
+                d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
+                tot_zneb(i) = tot_znebn(i)
+
+
+                !1) Cloudy to clear air
+                d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i)) 
+                IF (znebprecipcld(i) .GT. 0) THEN
+                        d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i) 
+                        d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i) 
+                ELSE 
+                        d_zrfl_cld_clr(i) = 0. 
+                        d_zifl_cld_clr(i) = 0. 
+                ENDIF
+
+                !2) Clear to cloudy air
+                d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) &
+                        - d_tot_zneb(i) - zneb(i)))
+                IF (znebprecipclr(i) .GT. 0) THEN
+                        d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i) 
+                        d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
+                ELSE
+                        d_zrfl_clr_cld(i) = 0. 
+                        d_zifl_clr_cld(i) = 0. 
+                ENDIF 
+
+                !Update variables
+                znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)  
+                znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i) 
+                zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i) 
+                ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
+                zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
+                ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
+
+        ENDDO
+    ENDIF
+
+    ! Initialisation of zoliq and radliq variables
+
+    DO i = 1, klon
+        IF (rneb(i,k).GT.0.0) THEN
+            zoliq(i) = zcond(i)
+            zrho(i) = pplay(i,k) / zt(i) / RD
+            zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
+    
+            zneb(i) = MAX(rneb(i,k), seuil_neb)
+            radliq(i,k) = zoliq(i)/REAL(ninter+1)
+            radicefrac(i,k) = zfice(i)
+            radliqi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1)
+        ENDIF
+    ENDDO
+
+    ! -------------------------------------------------------------------------------
+    ! Iterations to calculate a "mean" cloud water content during the precipitation
+    ! Comment: it is not the remaining water after precipitation but a mean
+    ! remaining water in the cloud during the time step that is seen by the radiation
+    ! -------------------------------------------------------------------------------
+
+     
+    DO n = 1, ninter
+
+        DO i=1,klon
+            IF (rneb(i,k).GT.0.0) THEN
+                zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
+            ENDIF
+        ENDDO
+
+        CALL FALLICE_VELOCITY(klon,zrhol(:),zt(:),zrho(:),pplay(:,k),ptconv(:,k),velo(:))
+
+        DO i = 1, klon
+
+            IF (rneb(i,k).GT.0.0) THEN
+
+                ! Initialization of zpluie, zice and ztot:
+                zpluie=0.
+                zice=0.
+                ztot=0.
+
+                IF (zneb(i).EQ.seuil_neb) THEN
+                    ztot = 0.0
+                    zice = 0.0
+                    zpluie= 0.0
+                ELSE
+                    ! water quantity to remove: zchau (Sundqvist, 1978)
+                    ! same thing for the ice: zfroi (Zender & Kiehl, 1997)
+                    IF (ptconv(i,k)) THEN ! if convective point
+                        zcl=cld_lc_con
+                        zct=1./cld_tau_con
+                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i)*zfice(i)
+                    ELSE
+                        zcl=cld_lc_lsc
+                        zct=1./cld_tau_lsc
+                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
+                            *velo(i) * zfice(i)
+                    ENDIF
+
+                    ! if vertical heterogeneity is taken into account, we use
+                    ! the "true" volume fraction instead of a modified 
+                    ! surface fraction (which is larger and artificially
+                    ! reduces the in-cloud water).
+                    IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
+                        zchau = zct   *dtime/REAL(ninter) * zoliq(i) &
+                        *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl)**2)) *(1.-zfice(i))
+                    ELSE
+                        zchau = zct   *dtime/REAL(ninter) * zoliq(i) &
+                        *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl)**2)) *(1.-zfice(i)) ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
+                    ENDIF
+
+                    zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
+                    zice   = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i)) 
+                    ztot   = MAX(zpluie + zice,0.0)
+
+                ENDIF
+            
+                ztot    = MIN(ztot,zoliq(i))
+                zice    = MIN(zice,ztot)
+                zpluie  = MIN(zpluie,ztot)
+
+                zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
+                zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
+                zoliq(i)  = MAX(zoliq(i)-ztot , 0.0)
+
+                radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
+                radliqi(i,k) = radliqi(i,k) + zoliqi(i)/REAL(ninter+1)
+           
+            ENDIF
+
+        ENDDO  ! i = 1,klon
+
+    ENDDO      ! n = 1,niter
+
+    
+    ! Caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme
+    DO i=1,klon
+        IF (radliq(i,k) .GT. 0) THEN
+            radicefrac(i,k)=MIN(MAX(radliqi(i,k)/radliq(i,k),0.),1.)
+        ENDIF
+    ENDDO
+
+
+    ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
+    ! If T<0C, liquid precip are converted into ice, which leads to an increase in
+    ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
+    ! taken into account through a linearization of the equations and by approximating
+    ! the liquid precipitation process with a "threshold" process. We assume that 
+    ! condensates are not modified during this operation. Liquid precipitation is
+    ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases. 
+    ! Water vapor increases as well
+    ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
+
+
+
+    DO i = 1, klon
+
+            IF (rneb(i,k) .GT. 0.0) THEN
+
+                zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
+                zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
+                zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
+                coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
+                ! 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)
+                DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
+                zt(i)  = zt(i) + DeltaT
+                ! water vaporization due to temp. increase
+                Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
+                ! we add this vaporized water to the total vapor and we remove it from the precip
+                zq(i) = zq(i) +  Deltaq
+                ! The three "max" lines herebelow protect from rounding errors
+                zcond(i) = max( zcond(i)- Deltaq, 0. )
+                ! liquid precipitation freezes oe evaporates
+                Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
+                zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
+                ! iced water budget
+                zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
+           
+                d_ql(i,k) = (1-zfice(i))*zoliq(i)
+                d_qi(i,k) = zfice(i)*zoliq(i)
+        
+                IF (iflag_evap_prec.EQ.4) THEN                                                            
+                    zrflcld(i) = zrflcld(i)+zqprecl(i) &                                                   
+                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                                
+                    ziflcld(i) = ziflcld(i)+ zqpreci(i) &                                                  
+                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                            
+                    znebprecipcld(i) = rneb(i,k)                                                          
+                    zrfl(i) = zrflcld(i) + zrflclr(i)                                                      
+                    zifl(i) = ziflcld(i) + ziflclr(i)                                                    
+                ELSE                                 
+                    zrfl(i) = zrfl(i)+ zqprecl(i)        &
+                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
+                    zifl(i) = zifl(i)+ zqpreci(i)        &
+                   *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)  
+                ENDIF
+ 
+           ENDIF ! rneb>0
+
+   ENDDO
+
+
+
+                                                                                                         
+    ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min                                                                              
+    ! if iflag_evap_pre=4
+    IF (iflag_evap_prec.EQ.4) THEN  
+
+        DO i=1, klon                                        
+                
+
+            IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN     
+                znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
+                (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
+            ELSE
+                znebprecipclr(i)=0.                                                                   
+            ENDIF  
+
+                                                                                                         
+            IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN                                                  
+                znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
+                (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
+            ELSE                                                                                     
+                znebprecipcld(i)=0.                                                                    
+            ENDIF     
+
+
+        ENDDO       
+
+    ENDIF  
+
+    ! --------------------------------
+    ! End of precipitation formation
+    ! --------------------------------
+
+
+    ! Outputs:
+    ! Precipitation fluxes at layer interfaces
+    ! and temperature and vapor tendencies
+    DO i = 1, klon
+        psfl(i,k)=zifl(i) 
+        prfl(i,k)=zrfl(i)
+        d_q(i,k) = zq(i) - q(i,k)
+        d_t(i,k) = zt(i) - t(i,k)
+    ENDDO
+     
+
+
+    ! ----------------------------------------------------------------
+    ! P4> Wet scavenging
+    ! ----------------------------------------------------------------
+
+
+    !Scavenging through nucleation in the layer
+
+    DO i = 1,klon
+        
+        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
+            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
+        ELSE
+            beta(i,k) = 0.
+        ENDIF
+
+        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
+        * (paprs(i,k)-paprs(i,k+1))/RG
+
+
+        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
+
+
+            IF (t(i,k) .GE. t_glace_min) THEN
+                zalpha_tr = a_tr_sca(3)
+            ELSE
+                zalpha_tr = a_tr_sca(4)
+            ENDIF
+
+
+            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
+            pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
+            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi 
+    
+            ! Nucleation with a factor of -1 instead of -0.5
+            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
+            pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
+
+        ENDIF
+        
+    ENDDO      
+
+     
+    ! Scavenging through impaction in the underlying layer
+
+    DO kk = k-1, 1, -1
+
+        DO i = 1, klon
+
+            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
+
+                IF (t(i,kk) .GE. t_glace_min) THEN
+                    zalpha_tr = a_tr_sca(1)
+                ELSE
+                    zalpha_tr = a_tr_sca(2)
+                ENDIF
+
+              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
+              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
+              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
+
+            ENDIF
+
+        ENDDO
+
+    ENDDO
+     
+
+END DO
+
+!======================================================================
+!                      END OF VERTICAL LOOP
+!======================================================================
+  
+
+  ! Rain or snow at the surface (depending on the first layer temperature)
+  DO i = 1, klon
+      snow(i) = zifl(i)
+      rain(i) = zrfl(i)
+  ENDDO
+   
+  
+
+  IF (ncoreczq>0) THEN
+      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
+  ENDIF
+
+
+
+
+END SUBROUTINE LSCP
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+END MODULE LSCP_MOD
Index: LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90	(revision 3999)
+++ LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90	(revision 3999)
@@ -0,0 +1,353 @@
+MODULE LSCP_TOOLS_MOD
+
+    IMPLICIT NONE
+
+CONTAINS
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,velo)
+
+    ! Ref:
+    ! Stubenrauch, C. J., Bonazzola, M.,
+    ! Protopapadaki, S. E., & Musat, I. (2019).
+    ! New cloud system metrics to assess bulk
+    ! ice cloud schemes in a GCM. Journal of
+    ! Advances in Modeling Earth Systems, 11,
+    ! 3212–3234. https://doi.org/10.1029/2019MS001642
+    
+    
+    IMPLICIT NONE
+
+    INCLUDE "nuage.h"
+    INCLUDE "fisrtilp.h"
+
+    INTEGER, INTENT(IN) :: klon
+    REAL, INTENT(IN), DIMENSION(klon) :: iwc       ! specific ice water content [kg/m3]
+    REAL, INTENT(IN), DIMENSION(klon) :: temp      ! temperature [K]
+    REAL, INTENT(IN), DIMENSION(klon) :: rho       ! dry air density [kg/m3]
+    REAL, INTENT(IN), DIMENSION(klon) :: pres      ! air pressure [Pa]
+    LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv    ! convective point  [-]
+
+    REAL, INTENT(OUT), DIMENSION(klon) :: velo    ! fallspeed velocity of crystals [m/s]
+
+
+    INTEGER i
+    REAL logvm,iwcg,tempc,phpa,cvel,dvel,fallv_tun
+    REAL m2ice, m2snow, vmice, vmsnow
+    REAL aice, bice, asnow, bsnow
+    
+
+    DO i=1,klon
+
+        IF (ptconv(i)) THEN
+            fallv_tun=ffallv_con
+        ELSE
+            fallv_tun=ffallv_lsc
+        ENDIF
+
+        tempc=temp(i)-273.15 ! celcius temp
+        iwcg=iwc(i)*1000.    ! iwc in g/m3
+        phpa=pres(i)/100.    ! pressure in hPa
+
+    IF (iflag_vice .EQ. 1) THEN
+        ! so-called 'empirical parameterization' in Stubenrauch et al. 2019
+        if (tempc .GE. -60.0) then
+
+            logvm= -0.0000414122*tempc*tempc*log(iwcg)-0.00538922*tempc*log(iwcg) &
+                    -0.0516344*log(iwcg)+0.00216078*tempc + 1.9714    
+            velo(i)=exp(logvm)
+        else
+            velo(i)=65.0*(iwcg**0.2)*(150./phpa)**0.15
+        endif
+        
+        velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s
+        dvel=0.2
+        cvel=fallv_tun*65.0*(rho(i)**0.2)*(150./phpa)**0.15
+
+    ELSE IF (iflag_vice .EQ. 2) THEN
+        ! so called  PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019
+        aice=0.587
+        bice=2.45
+        asnow=0.0444
+        bsnow=2.1
+        
+        m2ice=((iwcg*0.001/aice)/(exp(13.6-bice*7.76+0.479*bice**2)* & 
+                exp((-0.0361+bice*0.0151+0.00149*bice**2)*tempc)))   &
+                **(1./(0.807+bice*0.00581+0.0457*bice**2))
+
+        vmice=100.*1042.4*exp(13.6-(bice+1)*7.76+0.479*(bice+1.)**2)*exp((-0.0361+&
+                 (bice+1.)*0.0151+0.00149*(bice+1.)**2)*tempc)&
+                *(m2ice**(0.807+(bice+1.)*0.00581+0.0457*(bice+1.)**2))/(iwcg*0.001/aice)
+
+       
+        vmice=vmice*((1000./phpa)**0.2)
+      
+        m2snow=((iwcg*0.001/asnow)/(exp(13.6-bsnow*7.76+0.479*bsnow**2)* &
+               exp((-0.0361+bsnow*0.0151+0.00149*bsnow**2)*tempc)))         &
+               **(1./(0.807+bsnow*0.00581+0.0457*bsnow**2))
+
+
+        vmsnow=100.*14.3*exp(13.6-(bsnow+.416)*7.76+0.479*(bsnow+.416)**2)&
+                  *exp((-0.0361+(bsnow+.416)*0.0151+0.00149*(bsnow+.416)**2)*tempc)&
+                  *(m2snow**(0.807+(bsnow+.416)*0.00581+0.0457*(bsnow+.416)**2))/(iwcg*0.001/asnow)
+
+        vmsnow=vmsnow*((1000./phpa)**0.35)
+
+        velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
+        dvel=0.2
+        cvel=velo(i)/((iwc(i)*rho(i))**dvel)
+
+    ELSE
+        ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
+        velo(i) = fallv_tun*3.29/2.0 * ((iwc(i))**0.16)
+        dvel=0.16
+        cvel=fallv_tun*3.29/2.0*(rho(i)**0.16)
+    ENDIF
+
+    ENDDO
+
+END SUBROUTINE FALLICE_VELOCITY
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+SUBROUTINE ICEFRAC_LSCP(klon, temp, sig, icefrac, dicefracdT)
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  
+  ! Compute the ice fraction 1-xliq (see e.g.
+  ! Doutriaux-Boucher & Quaas 2004, section 2.2.)
+  ! as a function of temperature
+  ! see also Fig 3 of Madeleine et al. 2020, JAMES
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+    USE print_control_mod, ONLY: lunout, prt_level
+
+    IMPLICIT none
+
+
+    INCLUDE "YOMCST.h"
+    INCLUDE "nuage.h"
+    INCLUDE "clesphys.h"
+
+
+  ! nuage.h contains:
+  ! t_glace_min: if T < Tmin, the cloud is only made of water ice
+  ! t_glace_max: if T > Tmax, the cloud is only made of liquid water
+  ! exposant_glace: controls the sharpness of the transition
+
+    INTEGER, INTENT(IN)                 :: klon       ! number of horizontal grid points
+    REAL, INTENT(IN), DIMENSION(klon)   :: temp       ! temperature
+    REAL, INTENT(IN), DIMENSION(klon)   :: sig
+    REAL, INTENT(OUT), DIMENSION(klon)  :: icefrac
+    REAL, INTENT(OUT), DIMENSION(klon)  :: dicefracdT
+
+
+    INTEGER i
+    REAL    sig0,www,tmin_tmp,liqfrac_tmp
+    REAL    Dv, denomdep,beta,qsi,dqsidt
+    INTEGER exposant_glace_old
+    REAL t_glace_min_old
+    LOGICAL ice_thermo
+
+    sig0=0.8
+    t_glace_min_old = RTT - 15.0
+    ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
+    IF (ice_thermo) THEN
+      exposant_glace_old = 2
+    ELSE
+      exposant_glace_old = 6
+    ENDIF
+
+
+! calculation of icefrac and dicefrac/dT
+
+    DO i=1,klon
+
+    IF (iflag_t_glace.EQ.1) THEN
+            ! Transition to ice close to surface for T<Tmax
+            ! w=1 at the surface and 0 for sig < sig0
+            www=(max(sig(i)-sig0,0.))/(1.-sig0)
+    ELSEIF (iflag_t_glace.GE.2) THEN
+            ! No convertion to ice close to surface
+            www = 0.
+    ENDIF
+     
+    tmin_tmp=www*t_glace_max+(1.-www)*t_glace_min
+    liqfrac_tmp=  (temp(i)-tmin_tmp) / (t_glace_max-tmin_tmp)
+    liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
+       
+    IF (iflag_t_glace.GE.3) THEN
+            icefrac(i) = 1.0-liqfrac_tmp**exposant_glace
+            IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0)) THEN
+                 dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) &
+                           / (t_glace_min - t_glace_max)
+            ELSE
+
+                 dicefracdT(i)=0.
+            ENDIF
+
+    ELSE
+            icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace
+            IF (icefrac(i) .GT.0.) THEN
+                 dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) &
+                           / (t_glace_min - t_glace_max)
+            ENDIF
+
+            IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN
+                 dicefracdT(i)=0.
+            ENDIF
+
+    ENDIF
+    
+    ENDDO
+
+
+    RETURN
+
+END SUBROUTINE ICEFRAC_LSCP
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+
+SUBROUTINE CALC_QSAT_ECMWF(temp,qtot,pressure,tref,phase,flagth,qs,dqs)
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    ! Calculate qsat following ECMWF method
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+    IMPLICIT none
+
+    include "YOMCST.h"
+    include "YOETHF.h"
+    include "FCTTRE.h"
+
+    REAL, INTENT(IN) :: temp     ! temperature in K
+    REAL, INTENT(IN) :: qtot     ! total specific water in kg/kg
+    REAL, INTENT(IN) :: pressure ! pressure in Pa
+    REAL, INTENT(IN) :: tref     ! reference temperature in K
+    LOGICAL, INTENT(IN) :: flagth     ! flag for qsat calculation for thermals
+
+    INTEGER, INTENT(IN) :: phase 
+    ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
+    !        1=liquid
+    !        2=solid
+
+    REAL, INTENT(OUT) :: qs      ! saturation specific humidity [kg/kg]
+    REAL, INTENT(OUT) :: dqs     ! derivation of saturation specific humidity wrt T
+
+
+    REAL delta, cor, cvm5
+
+   
+    IF (phase .EQ. 1) THEN
+        delta=0.
+    ELSEIF (phase .EQ. 2) THEN
+        delta=1.
+    ELSE
+        delta=MAX(0.,SIGN(1.,tref-temp))
+    ENDIF
+
+    IF (flagth) THEN
+    cvm5=R5LES*(1.-delta) + R5IES*delta
+    ELSE
+    cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta
+    cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot))
+    ENDIF
+
+    qs= R2ES*FOEEW(temp,delta)/pressure
+    qs=MIN(0.5,qs)
+    cor=1./(1.-RETV*qs)
+    qs=qs*cor
+    dqs= FOEDE(temp,delta,cvm5,qs,cor)
+
+
+
+END SUBROUTINE CALC_QSAT_ECMWF
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+SUBROUTINE CALC_GAMMASAT(temp,qtot,pressure,gammasat,dgammasatdt)
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    ! programme that calculates the gammasat parameter that determines the
+    ! homogeneous condensation thresholds for cold (<0oC) clouds
+    ! condensation at q>gammasat*qsat
+    ! Etienne Vignon, March 2021
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+    IMPLICIT none
+
+    include "YOMCST.h"
+    include "YOETHF.h"
+    include "FCTTRE.h"
+    include "nuage.h" 
+
+
+    REAL, INTENT(IN) :: temp     ! temperature in K
+    REAL, INTENT(IN) :: qtot     ! total specific water in kg/kg
+
+    REAL, INTENT(IN) :: pressure ! pressure in Pa
+
+    REAL, INTENT(OUT) :: gammasat           ! coefficient to multiply qsat with to calculate saturation
+    REAL, INTENT(OUT) :: dgammasatdt       ! derivative of gammasat wrt temperature
+
+    REAL qsi,qsl,fac,dqsl,dqsi,fcirrus
+    REAL, PARAMETER :: acirrus=2.349
+    REAL, PARAMETER :: bcirrus=259.0
+
+
+        CALL CALC_QSAT_ECMWF(temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
+        CALL CALC_QSAT_ECMWF(temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
+
+        IF (temp .GE. RTT) THEN
+            ! warm clouds: condensation at saturation wrt liquid
+            gammasat=1.
+            dgammasatdt=0.
+
+        ELSEIF ((temp .LT. RTT) .AND. (temp .GT. t_glace_min)) THEN
+            
+            IF (iflag_gammasat .GE. 2) THEN          
+               gammasat=qsl/qsi
+               dgammasatdt=(dqsl*qsi-dqsi*qsl)/qsi/qsi
+            ELSE
+               gammasat=1.
+               dgammasatdt=0.
+            ENDIF
+
+        ELSE
+
+            IF (iflag_gammasat .GE.1) THEN
+            ! homogeneous freezing of aerosols, according to
+            ! Koop, 2000 and Karcher 2008, QJRMS
+            ! 'Cirrus regime'
+               fcirrus=acirrus-temp/bcirrus
+               IF (fcirrus .LT. qsl/qsi) THEN
+                  gammasat=qsl/qsi
+                  dgammasatdt=(dqsl*qsi-dqsi*qsl)/qsi/qsi
+               ELSE
+                  gammasat=fcirrus
+                  dgammasatdt=-1.0/bcirrus
+               ENDIF
+           
+            ELSE
+
+               gammasat=1.
+               dgammasatdt=0.
+
+            ENDIF
+
+        ENDIF
+   
+
+
+
+END SUBROUTINE CALC_GAMMASAT
+
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+END MODULE LSCP_TOOLS_MOD
+
+
Index: LMDZ6/trunk/libf/phylmd/newmicro.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/newmicro.F90	(revision 3998)
+++ LMDZ6/trunk/libf/phylmd/newmicro.F90	(revision 3999)
@@ -1,5 +1,5 @@
 ! $Id$
 
-SUBROUTINE newmicro(flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, pclc, &
+SUBROUTINE newmicro(flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, picefra, pclc, &
     pcltau, pclemi, pch, pcl, pcm, pct, pctlwp, xflwp, xfiwp, xflwc, xfiwc, &
     mass_solu_aero, mass_solu_aero_pi, pcldtaupi, re, fl, reliq, reice, &
@@ -9,9 +9,11 @@
   USE phys_local_var_mod, ONLY: scdnc, cldncl, reffclwtop, lcc, reffclws, &
       reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  & 
-      zfice, dNovrN
+      zfice, dNovrN, ptconv
   USE phys_state_var_mod, ONLY: rnebcon, clwcon
   USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
   USE ioipsl_getin_p_mod, ONLY : getin_p
   USE print_control_mod, ONLY: lunout
+  USE lscp_tools_mod, only: icefrac_lscp
+
 
 
@@ -31,4 +33,5 @@
   ! pqlwp---input-R-eau liquide nuageuse dans l'atmosphere dans la partie
   ! nuageuse (kg/kg)
+  ! picefra--input-R-fraction de glace dans les nuages
   ! pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
   ! mass_solu_aero-----input-R-total mass concentration for all soluble
@@ -58,4 +61,5 @@
   include "radepsi.h"
   include "radopt.h"
+  include "clesphys.h"
 
   ! choix de l'hypothese de recouvrement nuageuse via radopt.h (IM, 19.07.2016)
@@ -81,5 +85,5 @@
   REAL t(klon, klev)
   REAL pclc(klon, klev)
-  REAL pqlwp(klon, klev)
+  REAL pqlwp(klon, klev), picefra(klon,klev)
   REAL pcltau(klon, klev)
   REAL pclemi(klon, klev)
@@ -148,4 +152,5 @@
   ! jq-end
   ! IM cf. CR:parametres supplementaires
+  REAL dzfice(klon,klev)
   REAL zclear(klon)
   REAL zcloud(klon)
@@ -229,11 +234,24 @@
   ELSE ! of IF (iflag_t_glace.EQ.0)
     DO k = 1, klev
-        CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:,k))
- 
-
-        ! JBM: icefrac_lsc is now contained icefrac_lsc_mod
+
+! JBM: icefrac_lsc is now contained icefrac_lsc_mod
 !       zfice(i, k) = icefrac_lsc(t(i,k), t_glace_min, &
 !                                 t_glace_max, exposant_glace)
-      DO i = 1, klon
+
+      IF (ok_new_lscp) THEN
+          CALL icefrac_lscp(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:,k),dzfice(:,k))
+      ELSE
+          CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:,k))
+      ENDIF
+
+      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
+        ! consistent only for non convective grid points
+        ! critical for mixed phase clouds
+            zfice(i,k)=picefra(i,k)
+        ENDIF
+
         ! -layer calculation
         rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
Index: LMDZ6/trunk/libf/phylmd/nuage.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/nuage.F90	(revision 3998)
+++ LMDZ6/trunk/libf/phylmd/nuage.F90	(revision 3999)
@@ -1,9 +1,11 @@
 ! $Id$
 
-SUBROUTINE nuage(paprs, pplay, t, pqlwp, pclc, pcltau, pclemi, pch, pcl, pcm, &
+SUBROUTINE nuage(paprs, pplay, t, pqlwp,picefra, pclc, pcltau, pclemi, pch, pcl, pcm, &
     pct, pctlwp, ok_aie, mass_solu_aero, mass_solu_aero_pi, bl95_b0, bl95_b1, &
     cldtaupi, re, fl)
   USE dimphy
+  USE lscp_tools_mod, only: icefrac_lscp
   USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
+  USE phys_local_var_mod, ONLY: ptconv
   IMPLICIT NONE
   ! ======================================================================
@@ -14,4 +16,5 @@
   ! t-------input-R-temperature
   ! pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
+  ! picefra--inout-R-fraction de glace dans les nuages (-)
   ! pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
   ! ok_aie--input-L-apply aerosol indirect effect or not
@@ -36,4 +39,5 @@
   include "YOMCST.h"
   include "nuage.h" ! JBM 3/14
+  include "clesphys.h"
 
   REAL paprs(klon, klev+1), pplay(klon, klev)
@@ -41,5 +45,5 @@
 
   REAL pclc(klon, klev)
-  REAL pqlwp(klon, klev)
+  REAL pqlwp(klon, klev), picefra(klon,klev)
   REAL pcltau(klon, klev), pclemi(klon, klev)
 
@@ -89,4 +93,5 @@
 
   REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
+  REAl dzfice(klon)
   ! jq-end
 
@@ -106,5 +111,21 @@
 !       zfice(i) = icefrac_lsc(t(i,k), t_glace_min, &
 !                           t_glace_max, exposant_glace)
-         CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:))
+        IF (ok_new_lscp) THEN
+            CALL icefrac_lscp(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
+        ELSE
+            CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:))
+
+        ENDIF
+
+        IF ((.NOT. ptconv(i,k)) .AND. ok_new_lscp .AND. ok_icefra_lscp) THEN
+        ! EV: take the ice fraction directly from the lscp code
+        ! consistent only for non convective grid points
+        ! critical for mixed phase clouds
+            DO i=1,klon
+            zfice(i)=picefra(i,k)
+            ENDDO
+        ENDIF
+
+
      ENDIF
 
Index: LMDZ6/trunk/libf/phylmd/nuage.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/nuage.h	(revision 3998)
+++ LMDZ6/trunk/libf/phylmd/nuage.h	(revision 3999)
@@ -11,4 +11,7 @@
       INTEGER iflag_t_glace, iflag_cloudth_vert, iflag_cld_cv
       INTEGER iflag_rain_incloud_vol
+   
+      INTEGER iflag_mpc_bl, iflag_gammasat, iflag_vice
+      LOGICAL ok_icefra_lscp
 
       common /nuagecom/ rad_froid,rad_chau1, rad_chau2,t_glace_max,     &
@@ -17,4 +20,6 @@
      &                  tmax_fonte_cv,                                  &
      &                  iflag_t_glace,iflag_cloudth_vert,iflag_cld_cv,  &
-     &                  iflag_rain_incloud_vol
+     &                  iflag_rain_incloud_vol,                         &
+     &                  ok_icefra_lscp,                                 &
+     &                  iflag_mpc_bl, iflag_gammasat, iflag_vice    
 !$OMP THREADPRIVATE(/nuagecom/)
Index: LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 3998)
+++ LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 3999)
@@ -434,4 +434,6 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: cldemi, cldfra, cldtau, fiwc, fl, re, flwc
 !$OMP THREADPRIVATE(cldemi, cldfra, cldtau, fiwc, fl, re, flwc)
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: qlth, qith
+!$OMP THREADPRIVATE(qlth, qith)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: ref_liq, ref_ice, theta, zphi
 !$OMP THREADPRIVATE(ref_liq, ref_ice, theta, zphi)
@@ -730,5 +732,6 @@
       ALLOCATE(rain_lsc(klon))
       ALLOCATE(rain_num(klon))
-!
+      ALLOCATE(qlth(klon,klev), qith(klon,klev))
+      !
       ALLOCATE(sens_x(klon), sens_w(klon))
       ALLOCATE(zxfluxlat_x(klon), zxfluxlat_w(klon))
@@ -1029,4 +1032,5 @@
       DEALLOCATE(rain_lsc)
       DEALLOCATE(rain_num)
+      DEALLOCATE(qlth, qith)
 !
       DEALLOCATE(sens_x, sens_w)
Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 3998)
+++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 3999)
@@ -74,4 +74,5 @@
     USE VERTICAL_LAYERS_MOD, ONLY: aps,bps, ap, bp
     USE write_field_phy
+    USE lscp_mod, ONLY : lscp
 
     !USE cmp_seri_mod
@@ -960,4 +961,6 @@
     !IM cf. AM 081204 BEG
     LOGICAL ptconvth(klon,klev)
+
+    REAL picefra(klon,klev)
     !IM cf. AM 081204 END
     !
@@ -3510,5 +3513,20 @@
        print *,'itap, ->fisrtilp ',itap
     ENDIF
-    !
+
+    picefra(:,:)=0.
+
+    IF (ok_new_lscp) THEN
+
+    CALL lscp(phys_tstep,paprs,pplay, &
+         t_seri, q_seri,ptconv,ratqs, &
+         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, picefra, &
+         rain_lsc, snow_lsc, &
+         pfrac_impa, pfrac_nucl, pfrac_1nucl, &
+         frac_impa, frac_nucl, beta_prec_fisrt, &
+         prfl, psfl, rhcl,  &
+         zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
+         iflag_ice_thermo)
+
+    ELSE
     CALL fisrtilp(phys_tstep,paprs,pplay, &
          t_seri, q_seri,ptconv,ratqs, &
@@ -3520,5 +3538,5 @@
          zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
          iflag_ice_thermo)
-    !
+    ENDIF
     WHERE (rain_lsc < 0) rain_lsc = 0.
     WHERE (snow_lsc < 0) snow_lsc = 0.
@@ -4030,5 +4048,5 @@
           ENDIF
           CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
-               paprs, pplay, t_seri, cldliq, cldfra, &
+               paprs, pplay, t_seri, cldliq, picefra, cldfra, &
                cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
                flwp, fiwp, flwc, fiwc, &
@@ -4038,5 +4056,5 @@
        ELSE
           CALL nuage (paprs, pplay, &
-               t_seri, cldliq, cldfra, cldtau, cldemi, &
+               t_seri, cldliq, picefra, cldfra, cldtau, cldemi, &
                cldh, cldl, cldm, cldt, cldq, &
                ok_aie, &
