Index: LMDZ6/trunk/libf/phylmd/calcratqs_multi_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/calcratqs_multi_mod.F90	(revision 4071)
+++ LMDZ6/trunk/libf/phylmd/calcratqs_multi_mod.F90	(revision 4072)
@@ -230,6 +230,5 @@
 
 INTEGER :: i,k,nsrf
-REAL :: ratiom, qsat2m, dqsatdT
-REAL, DIMENSION(klon) :: xsi0
+REAL, DIMENSION(klon) :: xsi0, ratiom, qsat2m, dqsatdT
 REAL, DIMENSION (klon,klev) :: zlay
 
@@ -239,17 +238,15 @@
 !-------------------------------------------
 
-DO i=1,klon
     
-    ratiom=0.
-    xsi0(i)=0.
+    ratiom(:)=0.
+    xsi0(:)=0.
     
     DO nsrf=1,nbsrf
-    CALL CALC_QSAT_ECMWF(t2m(i,nsrf),q2m(i,nsrf),paprs(i,1),RTT,0,.false.,qsat2m,dqsatdT)
-    ratiom=ratiom+pctsrf(i,nsrf)*(q2m(i,nsrf)/qsat2m)
-    xsi0(i)=xsi0(i)+pctsrf(i,nsrf)*((q2m(i,nsrf)/qsat2m-ratiom)**2)
+    CALL CALC_QSAT_ECMWF(klon,t2m(:,nsrf),q2m(:,nsrf),paprs(:,1),RTT,0,.false.,qsat2m,dqsatdT)
+    ratiom(:)=ratiom(:)+pctsrf(:,nsrf)*(q2m(:,nsrf)/qsat2m(:))
+    xsi0(:)=xsi0(:)+pctsrf(:,nsrf)*((q2m(:,nsrf)/qsat2m(:)-ratiom(:))**2)
     END DO
     
-    xsi0(i)=sqrt(xsi0(i))/(ratiom+1E-6)
-END DO
+    xsi0(:)=sqrt(xsi0(:))/(ratiom(:)+1E-6)
 
 
Index: LMDZ6/trunk/libf/phylmd/cloudth_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/cloudth_mod.F90	(revision 4071)
+++ LMDZ6/trunk/libf/phylmd/cloudth_mod.F90	(revision 4072)
@@ -1539,5 +1539,5 @@
 
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points,                     &
+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)
@@ -1554,6 +1554,6 @@
       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
+      USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF
+      USE phys_local_var_mod, ONLY : qlth, qith, qsith
 
       IMPLICIT NONE
@@ -1573,5 +1573,5 @@
 
       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]
@@ -1589,5 +1589,6 @@
       REAL, DIMENSION(klon,klev+1), INTENT(IN)    ::  snowflux      ! snow flux at the interface of the layer [kg/m2/s]
 
-
+      INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points  ! grid points with convective (thermals) mixed phase clouds
+      
       REAL, DIMENSION(klon,klev), INTENT(OUT)      ::  ctot         ! Cloud fraction [0-1]
       REAL, DIMENSION(klon,klev), INTENT(OUT)      ::  ctot_vol     ! Volume cloud fraction [0-1]
@@ -1600,8 +1601,8 @@
 
       INTEGER itap,ind1,l,ig,iter,k 
-      LOGICAL flag_topthermals
-
-
-      REAL zqsatth(klon,klev), zqsatenv(klon,klev)
+      INTEGER iflag_topthermals, niter
+
+
+      REAL zqsatth(klon), zqsatenv(klon), zqsatenvonly(klon)
       REAL sigma1(klon,klev)                                                         
       REAL sigma2(klon,klev)
@@ -1614,6 +1615,5 @@
       REAL cenv_vol(klon,klev)
       REAL rneb(klon,klev)
-      REAL zqenv(klon)   
-
+      REAL zqenv(klon), zqth(klon), zqenvonly(klon) 
       REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
       REAL rdd,cppd,Lv
@@ -1624,6 +1624,8 @@
       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 zdelta,qsatbef,zcor
+      REAL Tbefth(klon), Tbefenv(klon), Tbefenvonly(klon)
+      REAL qlbef
+      REAL dqsatenv(klon), dqsatth(klon), dqsatenvonly(klon)
       REAL erf
       REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
@@ -1632,12 +1634,11 @@
       REAL zrho(klon,klev)
       REAL dz(klon,klev)
-      REAL qslth, qsith, qslenv, alenvl, aenvl
-      REAL sthi, sthl, althl, athl
+      REAL qslenv(klon), qslth(klon) 
+      REAL alenvl, aenvl
+      REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
       REAL senvi, senvl, qbase, sbase, qliqth, qiceth
-      REAL qimax, ttarget, stmp, cout, coutref
+      REAL qmax, ttarget, stmp, cout, coutref, fraci
       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
+      REAL deltazlev_mpc(klev), temp_mpc(klev), pres_mpc(klev), fraca_mpc(klev+1), snowf_mpc(klev+1)
 
       ! Modifty the saturation deficit PDF in thermals 
@@ -1645,4 +1646,6 @@
       REAL,SAVE :: C_mpc
       !$OMP THREADPRIVATE(C_mpc)
+      REAL, SAVE    :: Ni,C_cap,Ei,d_top
+      !$OMP THREADPRIVATE(Ni,C_cap,Ei,d_top)
       ! Change the width of the PDF used for vertical subgrid scale heterogeneity
       ! (J Jouhaud, JL Dufresne, JB Madeleine)
@@ -1709,4 +1712,13 @@
       sqrtpi=sqrt(pi)
       icefrac(:,ind2)=0.
+      althl=0.
+      althi=0.
+      athl=0.
+      athi=0.
+      senvl=0.
+      senvi=0.
+      sthi=0.
+      sthl=0.
+      sthil=0.
 
 
@@ -1741,7 +1753,19 @@
         WRITE(*,*) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs
         ! Modifies the PDF in thermals when ice crystals are present
-        C_mpc=1.e2
+        C_mpc=1.e4
         CALL getin_p('C_mpc',C_mpc)
         WRITE(*,*) 'C_mpc = ', C_mpc
+        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=1.2
+        CALL getin_p('d_top', d_top)
+        WRITE(*,*) 'd_top = ', d_top
 
         firstcall=.FALSE.
@@ -1759,5 +1783,5 @@
       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. (iflag_mpc_bl .GE. 1) .AND. (ind2<=klev-2)  &
             .AND. (ztv(ind1,1).GT.ztv(ind1,2)) .AND.(fraca(ind1,ind2).GT.1.e-10)) THEN
                 mpc_bl_points(ind1,ind2)=1
@@ -1772,5 +1796,26 @@
 !-------------------------------------------------------------------------------  
 
-    DO ind1=1,klon
+! calculation of temperature, humidity and saturation specific humidity
+
+Tbefenv(:)=zthl(:,ind2)*zpspsk(:,ind2)
+Tbefth(:)=ztla(:,ind2)*zpspsk(:,ind2)
+Tbefenvonly(:)=temp(:,ind2)
+
+zqenv(:)=(po(:)-fraca(:,ind2)*zqta(:,ind2))/(1.-fraca(:,ind2)) !qt = a*qtth + (1-a)*qtenv 
+zqth(:)=zqta(:,ind2)
+zqenvonly(:)=po(:)
+
+
+CALL CALC_QSAT_ECMWF(klon,Tbefenvonly,zqenvonly,paprs(:,ind2),RTT,0,.false.,zqsatenvonly,dqsatenv)
+
+CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,0,.false.,zqsatenv,dqsatenv)
+CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,1,.false.,qslenv,dqsatenv)
+
+CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,1,.false.,qslth,dqsatth)
+CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,2,.false.,qsith(:,ind2),dqsatth)
+CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,0,.false.,zqsatth,dqsatth)
+
+
+  DO ind1=1,klon
 
 
@@ -1780,11 +1825,6 @@
 ! 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
+
+        IF (Tbefenv(ind1) .GE. RTT) THEN
                Lv=RLVTT
         ELSE
@@ -1793,14 +1833,14 @@
         
 
-        alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
+        alenv=(0.622*Lv*zqsatenv(ind1))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
         aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
-        senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
+        senv=aenv*(po(ind1)-zqsatenv(ind1))                          !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)     
+        alenvl=(0.622*RLVTT*qslenv(ind1))/(rdd*zthl(ind1,ind2)**2)     
         aenvl=1./(1.+(alenv*Lv/cppd))         
-        senvl=aenvl*(po(ind1)-qslenv)    
+        senvl=aenvl*(po(ind1)-qslenv(ind1))    
+        senvi=senv
         ENDIF
 
@@ -1808,9 +1848,6 @@
 ! 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
+
+        IF (Tbefth(ind1) .GE. RTT) THEN
             Lv=RLVTT
         ELSE
@@ -1819,16 +1856,17 @@
 
         
-        alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       
+        alth=(0.622*Lv*zqsatth(ind1))/(rdd*ztla(ind1,ind2)**2)       
         ath=1./(1.+(alth*Lv/cppd))                                                         
-        sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      
+        sth=ath*(zqta(ind1,ind2)-zqsatth(ind1))                      
 
        ! 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)  
+         althi=alth
+         althl=(0.622*RLVTT*qslth(ind1))/(rdd*ztla(ind1,ind2)**2)                    
+         athl=1./(1.+(alth*RLVTT/cppd)) 
+         athi=alth     
+         sthl=athl*(zqta(ind1,ind2)-qslth(ind1))   
+         sthi=athi*(zqta(ind1,ind2)-qsith(ind1,ind2))  
+         sthil=athi*(zqta(ind1,ind2)-qslth(ind1))
         ENDIF     
 
@@ -1923,5 +1961,5 @@
                 ctot(ind1,ind2)=0.
                 ctot_vol(ind1,ind2)=0.
-                qcloud(ind1)=zqsatenv(ind1,ind2)
+                qcloud(ind1)=zqsatenv(ind1)
                 qincloud(ind1)=0.
             ELSE               
@@ -1941,4 +1979,5 @@
             !...........
 
+            qiceth=0.
             deltazlev_mpc=dz(ind1,:)
             temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
@@ -1946,33 +1985,21 @@
             fraca_mpc=fraca(ind1,:)
             snowf_mpc=snowflux(ind1,:)
-            qth_mpc=zqta(ind1,:)
-            flag_topthermals=.FALSE.
+            iflag_topthermals=0
             IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0))  THEN
-                flag_topthermals = .TRUE.
+                iflag_topthermals = 1
+            ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) & 
+                    .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
+                iflag_topthermals = 2
+            ELSE 
+                iflag_topthermals = 0
             ENDIF
 
-            CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,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)
+            CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:),qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,snowf_mpc,fraca_mpc,qith(ind1,:))
+
+            ! qmax calculation
+            sigma2s=(sigma2s_factor*((MAX((sthl-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
             deltasth=athl*vert_alpha_th*sigma2s
-      
-            ! Liquid phase
-            !.............
             xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
-            xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)           
+            xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
             exp_xth1 = exp(-1.*xth1**2)
             exp_xth2 = exp(-1.*xth2**2)
@@ -1984,28 +2011,50 @@
             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)           
+            qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)
+
+
+            ! Liquid phase
+            !................
+            ! We account for the effect of ice crystals in thermals on sthl
+            ! and on the width of the distribution
+
+            sthlc=sthl*1./(1.+C_mpc*qith(ind1,ind2))  &
+                + (1.-1./(1.+C_mpc*qith(ind1,ind2))) * athl*(qsith(ind1,ind2)-qslth(ind1))  
+
+            sigma2sc=(sigma2s_factor*((MAX((sthlc-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
+            deltasthc=athl*vert_alpha_th*sigma2sc
+     
+            
+            xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
+            xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)           
             exp_xth1 = exp(-1.*xth1**2)
             exp_xth2 = exp(-1.*xth2**2)
-            IntJ=0.5*sthi*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
+            IntJ=0.5*sthlc*(1-erf(xth2))+(sigma2sc/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
-
+            IntI1=(((sthlc+deltasthc)**2+(sigma2sc)**2)/(8*deltasthc))*(erf(xth2)-erf(xth1))
+            IntI2=(sigma2sc**2/(4*deltasthc*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
+            IntI3=((sqrt2*sigma2sc*(sthlc+deltasthc))/(4*sqrtpi*deltasthc))*(exp_xth1-exp_xth2)
+            IntI1_CF=((sthlc+deltasthc)*(erf(xth2)-erf(xth1)))/(4*deltasthc)
+            IntI3_CF=(sqrt2*sigma2sc*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasthc)
+            qliqth=IntJ+IntI1+IntI2+IntI3
+            
+            qlth(ind1,ind2)=MAX(0.,qliqth)
 
             ! Condensed water
-            ! 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)
 
+
+            ! consistency with subgrid distribution
+           
+             IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
+                 fraci=qith(ind1,ind2)/qcth(ind1,ind2)
+                 qcth(ind1,ind2)=qmax
+                 qith(ind1,ind2)=fraci*qmax
+                 qlth(ind1,ind2)=(1.-fraci)*qmax
+             ENDIF
+
+            ! Cloud Fraction
+            !...............
             ! calculation of qbase which is the value of the water vapor within mixed phase clouds
             ! such that the total water in cloud = qbase+qliqth+qiceth
@@ -2015,6 +2064,7 @@
             
             ttarget=qcth(ind1,ind2)
-            mini=athl*(qsith-qslth)
-            maxi=0. 
+            mini= athl*(qsith(ind1,ind2)-qslth(ind1))
+            maxi= 0. !athl*(3.*sqrt(sigma2s))
+            niter=20
             pas=(maxi-mini)/niter
             stmp=mini
@@ -2031,5 +2081,5 @@
                 ENDIF
             ENDDO
-            qbase=MAX(0., sbase/athl+qslth)
+            qbase=MAX(0., sbase/athl+qslth(ind1))
 
             ! surface cloud fraction in thermals
@@ -2053,5 +2103,5 @@
               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)
+              IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
               IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
               IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
@@ -2059,4 +2109,6 @@
             ENDIF
               cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
+
+
 
             ! Environment
@@ -2110,5 +2162,6 @@
 
             
-            ! Thermals + environment 
+            ! 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)
@@ -2125,8 +2178,8 @@
                 ctot_vol(ind1,ind2)=0.
                 qincloud(ind1)=0.
-                qcloud(ind1)=zqsatenv(ind1,ind2)
+                qcloud(ind1)=zqsatenv(ind1)
             ELSE               
                 qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
-                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv)
+                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
                 qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
                             +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
@@ -2139,11 +2192,7 @@
 
       
-        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
+
+
+        IF (Tbefenvonly(ind1) .GE. RTT) THEN
                 Lv=RLVTT
         ELSE
@@ -2153,10 +2202,10 @@
 
         zthl(ind1,ind2)=temp(ind1,ind2)*(101325./paprs(ind1,ind2))**(rdd/cppd)
-        alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
+        alenv=(0.622*Lv*zqsatenvonly(ind1))/(rdd*zthl(ind1,ind2)**2)
         aenv=1./(1.+(alenv*Lv/cppd))
-        senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
+        senv=aenv*(po(ind1)-zqsatenvonly(ind1))
         sth=0.
       
-        sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
+        sigma1s=ratqs(ind1,ind2)*zqenvonly(ind1)
         sigma2s=0.
 
@@ -2169,8 +2218,8 @@
         IF (ctot(ind1,ind2).LT.1.e-3) THEN
           ctot(ind1,ind2)=0.
-          qcloud(ind1)=zqsatenv(ind1,ind2)
+          qcloud(ind1)=zqsatenvonly(ind1)
           qincloud(ind1)=0.
         ELSE 
-          qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
+          qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenvonly(ind1)
           qincloud(ind1)=MAX(qctot(ind1,ind2)/ctot(ind1,ind2),0.)
         ENDIF
@@ -2194,5 +2243,5 @@
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev,flag_topthermals,temp,pres,qth,qlth,qith,deltazlev,snowf,fraca,qi)
+SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp,pres,qth,qsith,qlth,deltazlev,snowf,fraca,qith)
 
 ! parameterization of ice for boundary
@@ -2260,26 +2309,22 @@
     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 ?
+    INTEGER, INTENT(IN) :: ind1,ind2, klev           ! horizontal and vertical indices and dimensions
+    INTEGER, INTENT(IN) :: iflag_topthermals         ! uppermost layer of thermals ?
+    REAL, INTENT(IN)    :: Ni                        ! ice number concentration [m-3]
+    REAL, INTENT(IN)    :: Ei                        ! ice-droplet collision efficiency
+    REAL, INTENT(IN)    :: C_cap                     ! ice crystal capacitance
+    REAL, INTENT(IN)    :: d_top                     ! cloud-top ice crystal mixing parameter
     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) :: qsith        ! saturation specific humidity wrt ice 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)
+    REAL,  DIMENSION(klev), INTENT(INOUT) :: qith       ! condensed ice water , thermals [kg/kg]
 
 
@@ -2287,41 +2332,12 @@
     REAL rho(klev)
     REAL unsurtaudet, unsurtaustardep, unsurtaurim
-    REAL qsl, qsi, dqs, AA, BB, Ka, Dv, rhoi
-    REAL  p0, t0
+    REAL qi, AA, BB, Ka, Dv, rhoi
+    REAL p0, t0, fp1, fp2
     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
@@ -2334,26 +2350,23 @@
 
 
-    IF (flag_topthermals) THEN ! uppermost thermals level, solve a third order polynomial with Cardan's method
+    IF (iflag_topthermals .GT. 0) THEN ! uppermost thermals levels
 
     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)
+    BB=1./(rho(ind2)*Dv*qsith(ind2))
+    unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2)-qsith(ind2))/qsith(ind2) &
+                    *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=fraca(ind2)*unsurtaustardep/MAX((d_top*unsurtaudet),1E-12)
     qi=MAX(qi,0.)**(3./2.)
 
@@ -2370,35 +2383,40 @@
     ! 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
+    BB=1./(rho(ind2p1)*Dv*qsith(ind2p1))
+    unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2p1)-qsith(ind2p1))/qsith(ind2p1)*4.*RPI/(AA+BB)*(6.*rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33)
+    dep_term=rho(ind2p1)*fraca(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
+    rim_term=rho(ind2p1)*fraca(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))
+    ! We assume that there is no solid precipitation outside thermals 
+    ! so the precipitation flux within thermals is equal to the precipitation flux
+    ! at mesh-scale divided by  thermals fraction
+    
+    fp2=0.
+    fp1=0.
+    IF (fraca(ind2p1) .GT. 0.) THEN
+    fp2=-snowf(ind2p2)/fraca(ind2p1) ! flux defined positive upward
+    fp1=-snowf(ind2p1)/fraca(ind2p1)
+    ENDIF
+
+    precip_term=-1./deltazlev(ind2p1)*(fp2-fp1)
 
     ! 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)
+          qi= 1./fm_therm(ind1,ind2p1)* &
+              (deltazlev(ind2p1)*(-rim_term-dep_term-det_term-precip_term) + &
+              fm_therm(ind1,ind2p2)*(qith(ind2p1)))
+    END IF
+
+    ENDIF ! top thermals
+
+    qith(ind2)=MAX(0.,qi)
 
     RETURN
Index: LMDZ6/trunk/libf/phylmd/conf_phys_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 4071)
+++ LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 4072)
@@ -175,5 +175,5 @@
     INTEGER,SAVE :: iflag_clw_omp
     REAL,SAVE :: cld_lc_lsc_omp,cld_lc_con_omp,cld_tau_lsc_omp,cld_tau_con_omp
-    REAL,SAVE :: ffallv_lsc_omp, ffallv_con_omp,coef_eva_omp
+    REAL,SAVE :: ffallv_lsc_omp, ffallv_con_omp,coef_eva_omp,coef_eva_i_omp
     LOGICAL,SAVE :: reevap_ice_omp
     INTEGER,SAVE :: iflag_pdf_omp
@@ -1131,4 +1131,12 @@
     CALL getin('coef_eva',coef_eva_omp)
     !
+    !Config Key  = coef_eva_i
+    !Config Desc =  
+    !Config Def  = 2.e-5
+    !Config Help = 
+    !
+    coef_eva_i_omp = coef_eva_omp
+    CALL getin('coef_eva_i',coef_eva_i_omp)
+    !
     !Config Key  = reevap_ice
     !Config Desc =  
@@ -2492,4 +2500,5 @@
     ffallv_con = ffallv_con_omp
     coef_eva = coef_eva_omp
+    coef_eva_i = coef_eva_i_omp
     reevap_ice = reevap_ice_omp
     iflag_pdf = iflag_pdf_omp
@@ -2914,4 +2923,5 @@
     WRITE(lunout,*) ' ffallv_con = ', ffallv_con
     WRITE(lunout,*) ' coef_eva = ', coef_eva
+    WRITE(lunout,*) ' coef_eva_i = ', coef_eva_i
     WRITE(lunout,*) ' reevap_ice = ', reevap_ice
     WRITE(lunout,*) ' iflag_pdf = ', iflag_pdf
Index: LMDZ6/trunk/libf/phylmd/fisrtilp.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/fisrtilp.h	(revision 4071)
+++ LMDZ6/trunk/libf/phylmd/fisrtilp.h	(revision 4072)
@@ -11,5 +11,5 @@
       REAL cld_tau_lsc,cld_tau_con
       REAL ffallv_lsc,ffallv_con
-      REAL coef_eva
+      REAL coef_eva,coef_eva_i
       LOGICAL reevap_ice
       INTEGER iflag_pdf
@@ -25,4 +25,5 @@
      &     ,ffallv_con                                                  &
      &     ,coef_eva                                                    &
+     &     ,coef_eva_i                                                  &
      &     ,reevap_ice                                                  &
      &     ,iflag_fisrtilp_qsat                                         &
Index: LMDZ6/trunk/libf/phylmd/lscp_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 4071)
+++ LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 4072)
@@ -6,7 +6,7 @@
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE LSCP(dtime,missing_val,                      & 
+SUBROUTINE LSCP(dtime,missing_val,                      &
      paprs,pplay,t,q,ptconv,ratqs,                      &
-     d_t, d_q, d_ql, d_qi, rneb, rneb_seri,             & 
+     d_t, d_q, d_ql, d_qi, rneb, rneb_seri,             &
      radliq, radicefrac, rain, snow,                    &
      pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
@@ -25,5 +25,5 @@
 !      
 ! This code is a new version of the fisrtilp.F90 routine, which is itself a 
-! fusion of 'first' (superrsaturation physics, P. LeVan K. Laval)
+! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
 ! and 'ilp' (il pleut, L. Li)
 !
@@ -126,5 +126,5 @@
   INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
                                                                    ! CR: if iflag_ice_thermo=2, only convection is active    
-  LOGICAL,                         INTENT(IN)   :: ok_ice_sursat   ! flag to determine if ice sursaturation is activated 
+  LOGICAL,                         INTENT(IN)   :: ok_ice_sursat   ! flag to determine if ice sursaturation is activated
 
   LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
@@ -138,5 +138,5 @@
   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]
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: zthl         ! liquid potential temperature [K]
 
   ! INPUT/OUTPUT variables
@@ -145,4 +145,5 @@
   REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratqs            ! function of pressure that sets the large-scale
                                                                     ! cloud PDF (sigma=ratqs*qt)
+  
 
   ! Input sursaturation en glace
@@ -204,18 +205,21 @@
   PARAMETER (ztfondue=278.15)
 
-  REAL, SAVE    :: rain_int_min=0.001                               ! Minimum local rain intensity [mm/s] before the decrease in  associates precipitation fraction 
+  REAL, SAVE    :: rain_int_min=0.001                               ! Minimum local rain intensity [mm/s] before the decrease in associated precipitation fraction 
   !$OMP THREADPRIVATE(rain_int_min)
 
-
-  
+  LOGICAL, SAVE :: ok_radliq_precip=.false.                         ! Inclusion of all precipitation in radliq (cloud water seen by radiation)
+  !$OMP THREADPRIVATE(ok_radliq_precip)
+
+
+
+
   ! LOCAL VARIABLES:
   !----------------
 
 
-  REAL qsl, qsi
+  REAL qsl(klon), qsi(klon)
   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)
@@ -230,11 +234,13 @@
   REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
   REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon)
-  REAL zoliqp(klon), zoliqi(klon)
+  REAL zoliql(klon), zoliqi(klon)
   REAL zt(klon)
-  REAL zdz(klon),zrho(klon),ztot, zrhol(klon)
+  REAL zdz(klon),zrho(klon),iwc(klon)
   REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon)
-  REAL zmelt,zpluie,zice
+  REAL zmelt,zrain,zsnow,zprecip
   REAL dzfice(klon)
-  REAL zsolid,qtot,dqsl,dqsi
+  REAL zsolid
+  REAL qtot(klon), qzero(klon)
+  REAL dqsl(klon),dqsi(klon)
   REAL smallestreal
   !  Variables for Bergeron process
@@ -262,11 +268,12 @@
   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 velo(klon,klev), vr
   REAL qlmpc, qimpc, rnebmpc
-  REAL radliqi(klon,klev)
+  REAL radliqi(klon,klev), radliql(klon,klev)
 
   INTEGER i, k, n, kk
   INTEGER n_i(klon), iter
   INTEGER ncoreczq  
+  INTEGER mpc_bl_points(klon,klev)
   INTEGER,SAVE :: itap=0
   !$OMP THREADPRIVATE(itap)
@@ -334,8 +341,10 @@
     CALL getin_p('seuil_neb',seuil_neb)
     CALL getin_p('rain_int_min',rain_int_min) 
+    CALL getin_p('ok_radliq_precip',ok_radliq_precip)
     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
+    WRITE(lunout,*) 'lscp, ok_radliq_precip:', ok_radliq_precip
 
     ! check for precipitation sub-time steps
@@ -397,4 +406,5 @@
 tot_znebn(:) = 0.0                                                                               
 d_tot_zneb(:) = 0.0   
+qzero(:)=0.0
 
 !--ice sursaturation
@@ -468,5 +478,4 @@
     ! --------------------------------------------------------------------
     ! A part of the precipitation coming from above is evaporated/sublimated.
-    !
     ! --------------------------------------------------------------------
 
@@ -474,4 +483,11 @@
     IF (iflag_evap_prec.GE.1) THEN
 
+        ! Calculation of saturation specific humidity
+        ! depending on temperature:
+        CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
+        ! wrt liquid water
+        CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
+        ! wrt ice
+        CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
 
         DO i = 1, klon
@@ -480,6 +496,4 @@
             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                                                          
@@ -503,10 +517,5 @@
                 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)
@@ -515,13 +524,13 @@
 
                 IF (iflag_evap_prec.EQ.3) THEN
-                    zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl) &
+                    zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
                     *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) &
+                     zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
                     *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)) &
+                    zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
                     *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                 ENDIF
@@ -531,18 +540,15 @@
                 *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) &
+                    zqevti = znebprecip(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) &
                     *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) &
+                     zqevti = znebprecipclr(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) &
                     *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)) &
+                    zqevti = 1.*coef_eva_i*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
                     *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                 ENDIF
@@ -566,19 +572,4 @@
 
 
-                ! 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))  &
@@ -659,11 +650,10 @@
     ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
     !-------------------------------------------------------
-       
-        DO i = 1, klon
+      
+     qtot(:)=zq(:)+zmqc(:)
+     CALL CALC_QSAT_ECMWF(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
+     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
@@ -671,5 +661,5 @@
             ENDIF
 
-        ENDDO
+     ENDDO
     
 
@@ -777,5 +767,5 @@
                 DO i=1,klon
 
-                    ! convergence = .true. since convergence is not satisfied
+                    ! convergence = .true. until when convergence is not satisfied
                     convergence(i)=ABS(DT(i)).GT.DDT0
                    
@@ -796,10 +786,22 @@
 
                         ! 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) 
+                        qtot(i)=zq(i)+zmqc(i)
+                    
+                      ENDIF
+
+                  ENDDO
+
+                  ! Calculation of saturation specific humidity and ce fraction
+                  CALL CALC_QSAT_ECMWF(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
+                  CALL CALC_GAMMASAT(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))
+                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
+                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
+                  CALL ICEFRAC_LSCP(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
+
+
+
+                  DO i=1,klon
+
+                      IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
 
                         zpdf_sig(i)=ratqs(i,k)*zq(i)
@@ -826,17 +828,18 @@
                           ENDIF
 
-                          rnebss(i,k)=0.0   !--ajout OB (necessaire car boucle de convergence sur le temps)
+                          rnebss(i,k)=0.0   !--added by OB (needed because of the convergence loop on time)
                           fcontrN(i,k)=0.0  !--idem
                           fcontrP(i,k)=0.0  !--idem
                           qss(i,k)=0.0      !--idem
-
+                       
                         ELSE
+               
                         !------------------------------------
-                        ! SURSATURATION EN GLACE
+                        ! ICE SUPERSATURATION
                         !------------------------------------
 
                         CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), &
-                             gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),               &
-                             rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                 &
+                             gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),               &  
+                             rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                 & 
                              Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )
 
@@ -850,12 +853,10 @@
 
 
-            ! 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(:))
+                       ! P2.2.2> Approximative calculation of temperature variation DT 
+                       !           due to condensation.
+                       ! Calculated variables: 
+                       ! dT : temperature change due to condensation
+                       !---------------------------------------------------------------
+
                 
                         IF (zfice(i).LT.1) THEN
@@ -896,5 +897,5 @@
 
         ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs)
-        CALL icefrac_lscp(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))
+        CALL ICEFRAC_LSCP(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))
 
         DO i=1, klon
@@ -1002,14 +1003,15 @@
 
     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)
+            zoliqi(i)= zoliq(i)*zfice(i)
+            zoliql(i)= zoliq(i)*(1.-zfice(i))
+            zrho(i)  = pplay(i,k) / zt(i) / RD
+            zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
+            iwc(i)   = 0. 
+            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
+            radliql(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ninter+1)
     ENDDO
 
@@ -1024,9 +1026,9 @@
         DO i=1,klon
             IF (rneb(i,k).GT.0.0) THEN
-                zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
+                iwc(i) = zrho(i) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
             ENDIF
         ENDDO
 
-        CALL FALLICE_VELOCITY(klon,zrhol(:),zt(:),zrho(:),pplay(:,k),ptconv(:,k),velo(:))
+        CALL FALLICE_VELOCITY(klon,iwc(:),zt(:),zrho(:),pplay(:,k),ptconv(:,k),velo(:,k))
 
         DO i = 1, klon
@@ -1034,13 +1036,13 @@
             IF (rneb(i,k).GT.0.0) THEN
 
-                ! Initialization of zpluie, zice and ztot:
-                zpluie=0.
-                zice=0.
-                ztot=0.
+                ! Initialization of zrain, zsnow and zprecip:
+                zrain=0.
+                zsnow=0.
+                zprecip=0.
 
                 IF (zneb(i).EQ.seuil_neb) THEN
-                    ztot = 0.0
-                    zice = 0.0
-                    zpluie= 0.0
+                    zprecip = 0.0
+                    zsnow = 0.0
+                    zrain= 0.0
                 ELSE
                     ! water quantity to remove: zchau (Sundqvist, 1978)
@@ -1049,10 +1051,10 @@
                         zcl=cld_lc_con
                         zct=1./cld_tau_con
-                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i)*zfice(i)
+                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i,k)*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)
+                            *velo(i,k) * zfice(i)
                     ENDIF
 
@@ -1061,4 +1063,5 @@
                     ! 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) &
@@ -1069,22 +1072,19 @@
                     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)
+                    zrain   = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
+                    zsnow   = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i)) 
+                    zprecip = MAX(zrain + zsnow,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)
+                zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
+                zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
+                zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
 
                 radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
+                radliql(i,k) = radliql(i,k) + zoliql(i)/REAL(ninter+1)
                 radliqi(i,k) = radliqi(i,k) + zoliqi(i)/REAL(ninter+1)
            
-            ENDIF
+            ENDIF ! rneb >0
 
         ENDDO  ! i = 1,klon
@@ -1093,7 +1093,17 @@
 
     
-    ! Caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme
+    ! Include the contribution to non-evaporated precipitation to radliq 
+    IF ((ok_radliq_precip) .AND. (k .LT. klev)) THEN
+        ! for rain drops, we assume a fallspeed of 5m/s
+        vr=5.0
+        radliql(:,k)=radliql(:,k)+zrfl(:)/vr/zrho(:)
+        ! for ice crystals, we take the fallspeed from level above 
+        radliqi(:,k)=radliqi(:,k)+zifl(:)/zrho(:)/velo(:,k+1)
+        radliq(:,k)=radliql(:,k)+radliqi(:,k)
+    ENDIF
+
+    ! 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
+        IF (radliq(i,k) .GT. 0.) THEN
             radicefrac(i,k)=MIN(MAX(radliqi(i,k)/radliq(i,k),0.),1.)
         ENDIF
@@ -1101,45 +1111,44 @@
 
 
-    ! 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
-
-
+    ! Precipitation flux calculation
 
     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.)
+                
+
+               ! 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
+
+                    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 converted to ice precip 
+                    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) &                                                   
@@ -1165,5 +1174,5 @@
     IF (iflag_evap_prec.EQ.4) THEN  
 
-        DO i=1, klon                                        
+        DO i=1,klon                                        
                 
             IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN     
@@ -1192,8 +1201,10 @@
     ! Outputs:
     ! Precipitation fluxes at layer interfaces
-    ! and temperature and vapor tendencies
+    ! and temperature and water species tendencies
     DO i = 1, klon
         psfl(i,k)=zifl(i) 
         prfl(i,k)=zrfl(i)
+        d_ql(i,k) = (1-zfice(i))*zoliq(i)
+        d_qi(i,k) = zfice(i)*zoliq(i)
         d_q(i,k) = zq(i) - q(i,k)
         d_t(i,k) = zt(i) - t(i,k)
@@ -1264,24 +1275,22 @@
     ENDDO
      
-    !--save some variables for ice sursaturation
+    !--save some variables for ice supersaturation
     !
     DO i = 1, klon
-        ! pour la mémoire
+        ! for memory 
         rneb_seri(i,k) = rneb(i,k)
 
-        ! pour les diagnostics
+        ! for diagnostics
         rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)
 
         qvc(i,k) = zqs(i) * rneb(i,k)
-        qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k))  !--ajout OB a cause de cas pathologiques avec lognormale=F
+        qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k))  !--added by OB because of pathologic cases with the lognormal
         qcld(i,k) = qvc(i,k) + zcond(i)
-
-        !q_sat
-        CALL CALC_QSAT_ECMWF(Tbef(i),0.,pplay(i,k),RTT,1,.false.,zqsatl(i,k),zdqs(i))
-        CALL CALC_QSAT_ECMWF(Tbef(i),0.,pplay(i,k),RTT,2,.false.,zqsats(i,k),zdqs(i))
-
-     ENDDO
-
-ENDDO
+   END DO
+   !q_sat
+   CALL CALC_QSAT_ECMWF(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,zqsatl(:,k),zdqs(:))
+   CALL CALC_QSAT_ECMWF(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,zqsats(:,k),zdqs(:))     
+
+END DO
 
 !======================================================================
Index: LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90	(revision 4071)
+++ LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90	(revision 4072)
@@ -47,5 +47,5 @@
 
         tempc=temp(i)-273.15 ! celcius temp
-        iwcg=iwc(i)*1000.    ! iwc in g/m3
+        iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0
         phpa=pres(i)/100.    ! pressure in hPa
 
@@ -76,7 +76,7 @@
                 **(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=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)
 
        
@@ -93,16 +93,14 @@
 
         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)
+        cvel=velo(i)/((iwcg/1000.*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)
+        velo(i) = fallv_tun*3.29/2.0*((iwcg/1000.)**0.16)
         dvel=0.16
         cvel=fallv_tun*3.29/2.0*(rho(i)**0.16)
     ENDIF
-
     ENDDO
 
@@ -202,4 +200,5 @@
     ENDDO
 
+
     RETURN
 
@@ -209,8 +208,9 @@
 
 
-SUBROUTINE CALC_QSAT_ECMWF(temp,qtot,pressure,tref,phase,flagth,qs,dqs)
+SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     ! Calculate qsat following ECMWF method
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
 
     IMPLICIT NONE
@@ -220,10 +220,10 @@
     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
+    INTEGER, INTENT(IN) :: klon  ! number of horizontal grid points
+    REAL, INTENT(IN), DIMENSION(klon) :: temp     ! temperature in K
+    REAL, INTENT(IN), DIMENSION(klon) :: qtot     ! total specific water in kg/kg
+    REAL, INTENT(IN), DIMENSION(klon) :: 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)
@@ -231,9 +231,12 @@
     !        2=solid
 
-    REAL, INTENT(OUT) :: qs      ! saturation specific humidity [kg/kg]
-    REAL, INTENT(OUT) :: dqs     ! derivation of saturation specific humidity wrt T
+    REAL, INTENT(OUT), DIMENSION(klon) :: qs      ! saturation specific humidity [kg/kg]
+    REAL, INTENT(OUT), DIMENSION(klon) :: dqs     ! derivation of saturation specific humidity wrt T
 
     REAL delta, cor, cvm5
-   
+    INTEGER i
+
+    DO i=1,klon
+
     IF (phase .EQ. 1) THEN
         delta=0.
@@ -241,5 +244,5 @@
         delta=1.
     ELSE
-        delta=MAX(0.,SIGN(1.,tref-temp))
+        delta=MAX(0.,SIGN(1.,tref-temp(i)))
     ENDIF
 
@@ -248,12 +251,14 @@
     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)
+    cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
+    ENDIF
+
+    qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
+    qs(i)=MIN(0.5,qs(i))
+    cor=1./(1.-RETV*qs(i))
+    qs(i)=qs(i)*cor
+    dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
+
+    END DO
 
 END SUBROUTINE CALC_QSAT_ECMWF
@@ -262,5 +267,5 @@
 
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE CALC_GAMMASAT(temp,qtot,pressure,gammasat,dgammasatdt)
+SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
 
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -279,34 +284,38 @@
     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
+    INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
+    REAL, INTENT(IN), DIMENSION(klon) :: temp         ! temperature in K
+    REAL, INTENT(IN), DIMENSION(klon) :: qtot         ! total specific water in kg/kg
+
+    REAL, INTENT(IN), DIMENSION(klon) :: pressure     ! pressure in Pa
+
+    REAL, INTENT(OUT), DIMENSION(klon) :: gammasat    ! coefficient to multiply qsat with to calculate saturation
+    REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
+
+    REAL, DIMENSION(klon) ::  qsi,qsl,dqsl,dqsi
+    REAL  fcirrus, fac
     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
+    INTEGER i
+    
+        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
+        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
+
+    DO i=1,klon
+
+        IF (temp(i) .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
+            gammasat(i)=1.
+            dgammasatdt(i)=0.
+
+        ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN
             
             IF (iflag_gammasat .GE. 2) THEN          
-               gammasat=qsl/qsi
-               dgammasatdt=(dqsl*qsi-dqsi*qsl)/qsi/qsi
+               gammasat(i)=qsl(i)/qsi(i)
+               dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
             ELSE
-               gammasat=1.
-               dgammasatdt=0.
+               gammasat(i)=1.
+               dgammasatdt(i)=0.
             ENDIF
 
@@ -317,17 +326,17 @@
             ! 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
+               fcirrus=acirrus-temp(i)/bcirrus
+               IF (fcirrus .LT. qsl(i)/qsi(i)) THEN
+                  gammasat(i)=qsl(i)/qsi(i)
+                  dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
                ELSE
-                  gammasat=fcirrus
-                  dgammasatdt=-1.0/bcirrus
+                  gammasat(i)=fcirrus
+                  dgammasatdt(i)=-1.0/bcirrus
                ENDIF
            
             ELSE
 
-               gammasat=1.
-               dgammasatdt=0.
+               gammasat(i)=1.
+               dgammasatdt(i)=0.
 
             ENDIF
@@ -335,4 +344,7 @@
         ENDIF
    
+    END DO
+
+
 END SUBROUTINE CALC_GAMMASAT
 
Index: LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 4071)
+++ LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 4072)
@@ -438,6 +438,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(:,:) :: qlth, qith, qsith
+!$OMP THREADPRIVATE(qlth, qith, qsith)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: ref_liq, ref_ice, theta, zphi
 !$OMP THREADPRIVATE(ref_liq, ref_ice, theta, zphi)
@@ -458,4 +458,5 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: phiwriteSTD, qwriteSTD, twriteSTD, rhwriteSTD
 !$OMP THREADPRIVATE(phiwriteSTD, qwriteSTD, twriteSTD, rhwriteSTD)
+
 
 ! ug et d'autres encore:
@@ -777,5 +778,5 @@
       ALLOCATE(rain_lsc(klon))
       ALLOCATE(rain_num(klon))
-      ALLOCATE(qlth(klon,klev), qith(klon,klev))
+      ALLOCATE(qlth(klon,klev), qith(klon,klev), qsith(klon,klev))
       !
       ALLOCATE(sens_x(klon), sens_w(klon))
@@ -823,4 +824,5 @@
       ALLOCATE(ref_liq_pi(klon, klev), ref_ice_pi(klon, klev))
       ALLOCATE(zphi(klon, klev), zx_rh(klon, klev), zx_rhl(klon,klev), zx_rhi(klon,klev))
+      zx_rhl(:,:)=0.; zx_rhi(:,:)=0. ! because not always defined
       ALLOCATE(pmfd(klon, klev), pmfu(klon, klev))
 
@@ -1087,5 +1089,5 @@
       DEALLOCATE(rain_lsc)
       DEALLOCATE(rain_num)
-      DEALLOCATE(qlth, qith)
+      DEALLOCATE(qlth, qith, qsith)
 !
       DEALLOCATE(sens_x, sens_w)
Index: LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90	(revision 4071)
+++ LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90	(revision 4072)
@@ -559,6 +559,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(:,:) :: qlth, qith, qsith
+!$OMP THREADPRIVATE(qlth, qith, qsith)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: ref_liq, ref_ice, theta, zphi
 !$OMP THREADPRIVATE(ref_liq, ref_ice, theta, zphi)
@@ -943,5 +943,5 @@
       ALLOCATE(rain_lsc(klon))
       ALLOCATE(rain_num(klon))
-      ALLOCATE(qlth(klon,klev), qith(klon,klev))
+      ALLOCATE(qlth(klon,klev), qith(klon,klev), qsith(klon,klev))
       !
 #ifdef ISO
@@ -1004,4 +1004,5 @@
       ALLOCATE(ref_liq_pi(klon, klev), ref_ice_pi(klon, klev))
       ALLOCATE(zphi(klon, klev), zx_rh(klon, klev), zx_rhl(klon,klev), zx_rhi(klon,klev))
+      zx_rhl(:,:)=0.; zx_rhi(:,:)=0. ! because not always defined
       ALLOCATE(pmfd(klon, klev), pmfu(klon, klev))
 
@@ -1324,5 +1325,5 @@
       DEALLOCATE(rain_lsc)
       DEALLOCATE(rain_num)
-      DEALLOCATE(qlth, qith)
+      DEALLOCATE(qlth, qith, qsith)
 !
       DEALLOCATE(sens_x, sens_w)
