Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90	(revision 5380)
+++ LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90	(revision 5383)
@@ -7,5 +7,5 @@
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
-     paprs,pplay,temp,qt,qice_save,ptconv,ratqs,sigma_qtherm, &
+     paprs,pplay,omega,temp,qt,ptconv,ratqs,sigma_qtherm, &
      d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
      pfraclr, pfracld,                                  &
@@ -106,5 +106,4 @@
 USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb
 USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top
-!USE lmdz_lscp_condensation, ONLY : condensation_lognormal_test, condensation_ice_supersat
 USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
 USE lmdz_lscp_poprecip, ONLY : poprecip_precld, poprecip_postcld
@@ -138,6 +137,6 @@
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: omega           ! vertical velocity [Pa/s]
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg] 
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qice_save       ! ice specific from previous time step [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
@@ -267,4 +266,5 @@
   !----------------
   REAL,DIMENSION(klon) :: qsl, qsi                                ! saturation threshold at current vertical level
+  REAL,DIMENSION(klon) :: qice_ini
   REAL :: zct, zcl,zexpo
   REAL, DIMENSION(klon,klev) :: ctot
@@ -569,4 +569,5 @@
     ENDIF
 
+    qice_ini = zq(:)*ratio_qi_qtot(:,k)
     ! --------------------------------------------------------------------
     ! P2> Precipitation evaporation/sublimation/melting
@@ -579,9 +580,9 @@
         ! Calculation of saturation specific humidity
         ! depending on temperature:
-        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
+        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(:))
+        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(:))
+        CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsi,dqsi)
 
         DO i = 1, klon
@@ -763,5 +764,5 @@
 
      qtot(:)=zq(:)+zmqc(:)
-     CALL calc_qsat_ecmwf(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
+     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)))
@@ -874,10 +875,11 @@
 
         ! Treatment of non-boundary layer clouds (lognormale)
-        ! condensation with qsat(T) variation (adaptation)
-        ! Iterative resolution to converge towards qsat
-        ! with update of temperature, ice fraction and qsat at 
-        ! each iteration
-
-        ! todoan -> sensitivity to iflag_fisrtilp_qsat
+        ! We iterate here to take into account the change in qsat(T) and ice fraction
+        ! during the condensation process 
+        ! the increment in temperature is calculated using the first principle of 
+        ! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction
+        ! and a clear sky fraction)
+        ! note that we assume that the vapor in the cloud is at saturation for this calculation     
+
         DO iter=1,iflag_fisrtilp_qsat+1
 
@@ -913,6 +915,6 @@
 
                   ! Calculation of saturation specific humidity and ice 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(:))
+                  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(:)
@@ -923,5 +925,5 @@
                   ENDIF
 
-                  CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, zdistcltop(:),ztemp_cltop(:),zfice(:),dzfice(:))
+                  CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop,ztemp_cltop,zfice,dzfice)
 
                   !--AB Activates a condensation scheme that allows for
@@ -975,7 +977,4 @@
                        klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
                        keepgoing, rneb(:,k), zqn, qvc)
-!                   CALL condensation_lognormal_test( &
-!                       klon, klev, k, keepgoing, ratqs, gammasat, zq, zqs, &
-!                       zqn, qvc, rneb)
 
 
@@ -1057,8 +1056,6 @@
         ! Partition function depending on tke for non shallow-convective clouds 
         IF (iflag_icefrac .GE. 1) THEN
-
-           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), qice_save(:,k), ziflcld, zqn, &
+           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_ini, ziflcld, zqn, &
            rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
-
         ENDIF
 
Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.f90	(revision 5380)
+++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.f90	(revision 5383)
@@ -8,56 +8,4 @@
 
 CONTAINS
-
-!*********************************************************************************
-SUBROUTINE condensation_lognormal_test(klon,klev,k,keepgoing,ratqs, gammasat,zq, zqs, zqn, qvc, rneb)
-
-
-IMPLICIT NONE
-
-INTEGER, INTENT(IN)   :: klon,klev,k
-LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: keepgoing     ! .TRUE. if a new condensation loop should be computed
-REAL, INTENT(IN), DIMENSION(klon,klev)   :: ratqs
-REAL, INTENT(IN), DIMENSION(klon)        :: zq, zqs,gammasat
-REAL, INTENT(INOUT), DIMENSION(klon)     :: zqn, qvc
-REAL, INTENT(INOUT), DIMENSION(klon,klev) :: rneb
-
-REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta             ! lognormal parameters
-REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2          ! lognormal intermediate variables
-INTEGER               :: i
-
-
-DO i=1,klon !todoan : check if loop in i is needed
-
-    IF (keepgoing(i)) THEN
-
-    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)=zqs(i)
-        !--AB grid-mean vapor in the cloud - we assume saturation adjustment
-        qvc(i) = 0.
-    ELSE
-        rneb(i,k)=0.5*zpdf_e1(i)
-        zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
-        !--AB grid-mean vapor in the cloud - we assume saturation adjustment
-        qvc(i) = rneb(i,k) * zqs(i)
-    ENDIF
-
-    ENDIF ! keepgoing
-ENDDO ! loop on klon
-
-
-END SUBROUTINE condensation_lognormal_test
 
 
Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.f90	(revision 5380)
+++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.f90	(revision 5383)
@@ -222,33 +222,33 @@
   !$OMP THREADPRIVATE(rain_int_min)
 
-  REAL, SAVE, PROTECTED :: thresh_precip_frac=1.E-6         ! precipitation fraction threshold TODO [-]
+  REAL, SAVE, PROTECTED :: thresh_precip_frac=1.E-6         ! precipitation fraction threshold [-]
   !$OMP THREADPRIVATE(thresh_precip_frac)
 
-  REAL, SAVE, PROTECTED :: tau_mixenv=100000                ! Homogeneization time of mixed phase clouds [s]
-  !$OMP THREADPRIVATE(tau_mixenv)
-
-    REAL, SAVE, PROTECTED :: capa_crystal=1.                ! Sursaturation of ice part in mixed phase clouds [-]
+  REAL, SAVE, PROTECTED :: capa_crystal=1.                  ! Crystal capacitance (shape factor) for lscp_icefrac_turb [-]
   !$OMP THREADPRIVATE(capa_crystal)
-
-  REAL, SAVE, PROTECTED :: lmix_mpc=1000                    ! Length of turbulent zones in Mixed Phase Clouds [m]
-  !$OMP THREADPRIVATE(lmix_mpc)
 
   REAL, SAVE, PROTECTED :: naero5=0.5                       ! Number concentration of aerosol larger than 0.5 microns [scm-3]
   !$OMP THREADPRIVATE(naero5)
 
-  REAL, SAVE, PROTECTED :: gamma_snwretro = 0.              ! Proportion of snow taken into account in ice retroaction in icefrac_turb [-]
+  REAL, SAVE, PROTECTED :: gamma_snwretro = 0.              ! Proportion of snow taken into account in ice retroaction in lscp_icefrac_turb [-]
   !$OMP THREADPRIVATE(gamma_snwretro)
 
-  REAL, SAVE, PROTECTED :: gamma_taud = 1.                  ! Tuning coeff for tau_dissipturb [-]
+  REAL, SAVE, PROTECTED :: gamma_taud = 1.                  ! Tuning coeff for Lagrangian decorrelation timescale in lscp_icefrac_turb [-] 
   !$OMP THREADPRIVATE(gamma_taud)
 
-  REAL, SAVE, PROTECTED :: gamma_col=1.                     ! A COMMENTER TODO [-]
+  REAL, SAVE, PROTECTED :: gamma_col=1.                     ! Tuning coefficient for rain collection efficiency (poprecip) [-]
   !$OMP THREADPRIVATE(gamma_col)
 
-  REAL, SAVE, PROTECTED :: gamma_agg=1.                     ! A COMMENTER TODO [-]
+  REAL, SAVE, PROTECTED :: gamma_agg=1.                     ! Tuning coefficient for snow aggregation efficiency (poprecip) [-]
   !$OMP THREADPRIVATE(gamma_agg)
 
-  REAL, SAVE, PROTECTED :: gamma_rim=1.                     ! A COMMENTER TODO [-]
+  REAL, SAVE, PROTECTED :: gamma_rim=1.                     ! Tuning coefficient for riming efficiency (poprecip) [-]
   !$OMP THREADPRIVATE(gamma_rim)
+  
+  REAL, SAVE, PROTECTED :: gamma_melt=1.                    ! Tuning coefficient for snow melting efficiency (poprecip) [-]
+  !$OMP THREADPRIVATE(gamma_melt)
+  
+  REAL, SAVE, PROTECTED :: gamma_freez=1.                   ! Tuning coefficient for rain collision freezing efficiency (poprecip) [-]
+  !$OMP THREADPRIVATE(gamma_freez)
 
   REAL, SAVE, PROTECTED :: rho_rain=1000.                   ! Rain density [kg/m3]
@@ -258,47 +258,41 @@
   !$OMP THREADPRIVATE(rho_ice)
 
-  REAL, SAVE, PROTECTED :: r_rain=500.E-6                   ! Rain droplets radius for POPRECIP [m]
+  REAL, SAVE, PROTECTED :: r_rain=500.E-6                   ! Rain droplets radius (poprecip) [m]
   !$OMP THREADPRIVATE(r_rain)
 
-  REAL, SAVE, PROTECTED :: r_snow=1.E-3                     ! Ice crystals radius for POPRECIP [m]
+  REAL, SAVE, PROTECTED :: r_snow=1.E-3                     ! Ice crystals radius (poprecip) [m]
   !$OMP THREADPRIVATE(r_snow)
 
-  REAL, SAVE, PROTECTED :: tau_auto_snow_min=100.           ! A COMMENTER TODO [s]
+  REAL, SAVE, PROTECTED :: tau_auto_snow_min=100.           ! Snow autoconversion minimal timescale (when liquid) [s]
   !$OMP THREADPRIVATE(tau_auto_snow_min)
 
-  REAL, SAVE, PROTECTED :: tau_auto_snow_max=1000.          ! A COMMENTER TODO [s]
+  REAL, SAVE, PROTECTED :: tau_auto_snow_max=1000.          ! Snow autoconversion minimal timescale (when only ice) [s]
   !$OMP THREADPRIVATE(tau_auto_snow_max)
 
-  REAL, SAVE, PROTECTED :: eps=1.E-10                       ! A COMMENTER TODO [-]
+  REAL, SAVE, PROTECTED :: eps=1.E-10                       ! Treshold 0 [-]
   !$OMP THREADPRIVATE(eps)
 
-  REAL, SAVE, PROTECTED :: gamma_melt=1.                    ! A COMMENTER TODO [-]
-  !$OMP THREADPRIVATE(gamma_melt)
-
-  REAL, SAVE, PROTECTED :: alpha_freez=4.                   ! A COMMENTER TODO [-]
+  REAL, SAVE, PROTECTED :: alpha_freez=4.                   ! Slope of exponential for immersion freezing timescale [-]
   !$OMP THREADPRIVATE(alpha_freez)
 
-  REAL, SAVE, PROTECTED :: beta_freez=0.1                   ! A COMMENTER TODO [m-3.s-1]
+  REAL, SAVE, PROTECTED :: beta_freez=0.1                   ! Inv.time immersion freezing [s-1]
   !$OMP THREADPRIVATE(beta_freez)
 
-  REAL, SAVE, PROTECTED :: gamma_freez=1.                   ! A COMMENTER TODO [-]
-  !$OMP THREADPRIVATE(gamma_freez)
-
-  REAL, SAVE, PROTECTED :: rain_fallspeed=4.                ! A COMMENTER TODO [m/s]
+  REAL, SAVE, PROTECTED :: rain_fallspeed=4.                ! Rain fall velocity [m/s]
   !$OMP THREADPRIVATE(rain_fallspeed)
 
-  REAL, SAVE, PROTECTED :: rain_fallspeed_clr                ! A COMMENTER TODO [m/s]
+  REAL, SAVE, PROTECTED :: rain_fallspeed_clr               ! Rain fall velocity in clear sky [m/s]
   !$OMP THREADPRIVATE(rain_fallspeed_clr)
 
-  REAL, SAVE, PROTECTED :: rain_fallspeed_cld               ! A COMMENTER TODO [m/s]
+  REAL, SAVE, PROTECTED :: rain_fallspeed_cld               ! Rain fall velocity in cloudy sky [m/s]
   !$OMP THREADPRIVATE(rain_fallspeed_cld)
 
-  REAL, SAVE, PROTECTED :: snow_fallspeed=1.               ! A COMMENTER TODO [m/s]
+  REAL, SAVE, PROTECTED :: snow_fallspeed=1.                ! Snow fall velocity [m/s]
   !$OMP THREADPRIVATE(snow_fallspeed)
 
-  REAL, SAVE, PROTECTED :: snow_fallspeed_clr               ! A COMMENTER TODO [m/s]
+  REAL, SAVE, PROTECTED :: snow_fallspeed_clr               ! Snow fall velocity in clear sky [m/s]
   !$OMP THREADPRIVATE(snow_fallspeed_clr)
 
-  REAL, SAVE, PROTECTED :: snow_fallspeed_cld               ! A COMMENTER TODO [m/s]
+  REAL, SAVE, PROTECTED :: snow_fallspeed_cld               ! Snow fall velocity in cloudy sky [m/s]
   !$OMP THREADPRIVATE(snow_fallspeed_cld)
   !--End of the parameters for poprecip
@@ -386,7 +380,5 @@
     CALL getin_p('dist_liq',dist_liq)
     CALL getin_p('tresh_cl',tresh_cl)
-    CALL getin_p('tau_mixenv',tau_mixenv)
     CALL getin_p('capa_crystal',capa_crystal)
-    CALL getin_p('lmix_mpc',lmix_mpc)
     CALL getin_p('naero5',naero5)
     CALL getin_p('gamma_snwretro',gamma_snwretro)
@@ -403,4 +395,6 @@
     CALL getin_p('gamma_freez',gamma_freez)
     CALL getin_p('gamma_melt',gamma_melt)
+    CALL getin_p('tau_auto_snow_max',tau_auto_snow_max)
+    CALL getin_p('tau_auto_snow_min',tau_auto_snow_min)
     CALL getin_p('r_snow',r_snow)
     CALL getin_p('rain_fallspeed',rain_fallspeed)
@@ -471,7 +465,5 @@
     WRITE(lunout,*) 'lscp_ini, dist_liq', dist_liq
     WRITE(lunout,*) 'lscp_ini, tresh_cl', tresh_cl
-    WRITE(lunout,*) 'lscp_ini, tau_mixenv', tau_mixenv
     WRITE(lunout,*) 'lscp_ini, capa_crystal', capa_crystal
-    WRITE(lunout,*) 'lscp_ini, lmix_mpc', lmix_mpc
     WRITE(lunout,*) 'lscp_ini, naero5', naero5
     WRITE(lunout,*) 'lscp_ini, gamma_snwretro', gamma_snwretro
@@ -489,4 +481,6 @@
     WRITE(lunout,*) 'lscp_ini, gamma_freez:', gamma_freez
     WRITE(lunout,*) 'lscp_ini, gamma_melt:', gamma_melt
+    WRITE(lunout,*) 'lscp_ini, tau_auto_snow_max:',tau_auto_snow_max
+    WRITE(lunout,*) 'lscp_ini, tau_auto_snow_min:',tau_auto_snow_min
     WRITE(lunout,*) 'lscp_ini, r_snow:', r_snow
     WRITE(lunout,*) 'lscp_ini, rain_fallspeed_clr:', rain_fallspeed_clr
Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.f90	(revision 5380)
+++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.f90	(revision 5383)
@@ -240,9 +240,11 @@
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
-SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, temp, pplay, paprsdn, paprsup, qice_ini, snowcld, qtot_incl, cldfra, tke, tke_dissip, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb)
+
+SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, temp, pplay, paprsdn, paprsup, omega, qice_ini, snowcld, qtot_incl, cldfra, tke,   &
+                             tke_dissip, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb)
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   ! Compute the liquid, ice and vapour content (+ice fraction) based 
   ! on turbulence (see Fields 2014, Furtado 2016, Raillard 2025)
-  ! L.Raillard (30/08/24)
+  ! L.Raillard (23/09/24)
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
@@ -251,64 +253,63 @@
    USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI
    USE lmdz_lscp_ini, ONLY : seuil_neb, temp_nowater
-   USE lmdz_lscp_ini, ONLY : tau_mixenv, lmix_mpc, naero5, gamma_snwretro, gamma_taud, capa_crystal
+   USE lmdz_lscp_ini, ONLY : naero5, gamma_snwretro, gamma_taud, capa_crystal
    USE lmdz_lscp_ini, ONLY : eps
 
    IMPLICIT NONE
 
-   INTEGER,   INTENT(IN)                           :: klon              !--number of horizontal grid points
-   REAL,      INTENT(IN)                           :: dtime             !--time step [s]
-
-   REAL,      INTENT(IN),       DIMENSION(klon)    :: temp              !--temperature
-   REAL,      INTENT(IN),       DIMENSION(klon)    :: pplay             !--pressure in the middle of the layer       [Pa]
-   REAL,      INTENT(IN),       DIMENSION(klon)    :: paprsdn           !--pressure at the bottom interface of the layer [Pa]
-   REAL,      INTENT(IN),       DIMENSION(klon)    :: paprsup           !--pressure at the top interface of the layer [Pa]
-   REAL,      INTENT(IN),       DIMENSION(klon)    :: qtot_incl         !--specific total cloud water content, in-cloud content [kg/kg]
-   REAL,      INTENT(IN),       DIMENSION(klon)    :: cldfra            !--cloud fraction in gridbox [-]
-   REAL,      INTENT(IN),       DIMENSION(klon)    :: tke               !--turbulent kinetic energy [m2/s2]
-   REAL,      INTENT(IN),       DIMENSION(klon)    :: tke_dissip        !--TKE dissipation [m2/s3]
-
-   REAL,      INTENT(IN),       DIMENSION(klon)    :: qice_ini          !--initial specific ice content gridbox-mean [kg/kg]
+   INTEGER,   INTENT(IN)                           :: klon                !--number of horizontal grid points
+   REAL,      INTENT(IN)                           :: dtime               !--time step [s]
+
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: temp                !--temperature
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: pplay               !--pressure in the middle of the layer           [Pa]
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: paprsdn             !--pressure at the bottom interface of the layer [Pa]
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: paprsup             !--pressure at the top interface of the layer    [Pa]
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: omega               !--resolved vertical velocity                    [Pa/s]
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: qtot_incl           !--specific total cloud water in-cloud content   [kg/kg]
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: cldfra              !--cloud fraction in gridbox                     [-]
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: tke                 !--turbulent kinetic energy                      [m2/s2]
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: tke_dissip          !--TKE dissipation                               [m2/s3]
+
+   REAL,      INTENT(IN),       DIMENSION(klon)    :: qice_ini            !--initial specific ice content gridbox-mean     [kg/kg]
    REAL,      INTENT(IN),       DIMENSION(klon)    :: snowcld
-   REAL,      INTENT(OUT),      DIMENSION(klon)    :: qliq              !--specific liquid content gridbox-mean [kg/kg]
-   REAL,      INTENT(OUT),      DIMENSION(klon)    :: qvap_cld          !--specific cloud vapor content, gridbox-mean [kg/kg]
-   REAL,      INTENT(OUT),      DIMENSION(klon)    :: qice              !--specific ice content gridbox-mean [kg/kg]
-   REAL,      INTENT(OUT),      DIMENSION(klon)    :: icefrac           !--fraction of ice in condensed water [-]
+   REAL,      INTENT(OUT),      DIMENSION(klon)    :: qliq                !--specific liquid content gridbox-mean          [kg/kg]
+   REAL,      INTENT(OUT),      DIMENSION(klon)    :: qvap_cld            !--specific cloud vapor content, gridbox-mean    [kg/kg]
+   REAL,      INTENT(OUT),      DIMENSION(klon)    :: qice                !--specific ice content gridbox-mean             [kg/kg]
+   REAL,      INTENT(OUT),      DIMENSION(klon)    :: icefrac             !--fraction of ice in condensed water            [-]
    REAL,      INTENT(OUT),      DIMENSION(klon)    :: dicefracdT
 
-   REAL,      INTENT(OUT),      DIMENSION(klon)    :: cldfraliq         !--fraction of cldfra which is liquid only
-   REAL,      INTENT(OUT),      DIMENSION(klon)    :: sigma2_icefracturb     !--Temporary 
-   REAL,      INTENT(OUT),      DIMENSION(klon)    :: mean_icefracturb      !--Temporary 
-
-   REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati         !--specific humidity saturation values
+   REAL,      INTENT(OUT),      DIMENSION(klon)    :: cldfraliq           !--fraction of cldfra where liquid               [-]
+   REAL,      INTENT(OUT),      DIMENSION(klon)    :: sigma2_icefracturb  !--Sigma2 of the ice supersaturation PDF         [-]
+   REAL,      INTENT(OUT),      DIMENSION(klon)    :: mean_icefracturb    !--Mean of the ice supersaturation PDF           [-]
+
+   REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati           !--specific humidity saturation values
    INTEGER :: i
 
-   REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl                !--In-cloud specific quantities [kg/kg]
-   REAL :: qsnowcld_incl
-   !REAL :: capa_crystal                                                 !--Capacitance of ice crystals  [-]
-   REAL :: water_vapor_diff                                             !--Water-vapour diffusion coefficient in air [m2/s] (function of T&P)
-   REAL :: air_thermal_conduct                                          !--Thermal conductivity of air [J/m/K/s] (function of T)
-   REAL :: C0                                                           !--Lagrangian structure function [-]
-   REAL :: tau_mixingenv
+   REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl                  !--In-cloud specific quantities                  [kg/kg]
+   REAL :: water_vapor_diff                                               !--Water-vapour diffusion coeff in air (f(T,P))  [m2/s] 
+   REAL :: air_thermal_conduct                                            !--Thermal conductivity of air (f(T))            [J/m/K/s] 
+   REAL :: C0                                                             !--Lagrangian structure function                 [-]
    REAL :: tau_dissipturb
-   REAL :: invtau_phaserelax
+   REAL :: tau_phaserelax
    REAL :: sigma2_pdf, mean_pdf
    REAL :: ai, bi, B0
    REAL :: sursat_iceliq
-   REAL :: sursat_env
+   REAL :: sursat_equ
    REAL :: liqfra_max
    REAL :: sursat_iceext
-   REAL :: nb_crystals                                                  !--number concentration of ice crystals [#/m3]
-   REAL :: moment1_PSD                                                  !--1st moment of ice PSD 
-   REAL :: N0_PSD, lambda_PSD                                           !--parameters of the exponential PSD
-
-   REAL :: rho_ice                                                      !--ice density [kg/m3]
+   REAL :: nb_crystals                                                    !--number concentration of ice crystals          [#/m3]
+   REAL :: moment1_PSD                                                    !--1st moment of ice PSD 
+   REAL :: N0_PSD, lambda_PSD                                             !--parameters of the exponential PSD
+
+   REAL :: rho_ice                                                        !--ice density                                   [kg/m3]
    REAL :: cldfra1D
-   REAL :: deltaz, rho_air
-   REAL :: psati                                                        !--saturation vapor pressure wrt i [Pa]
+   REAL :: rho_air
+   REAL :: psati                                                          !--saturation vapor pressure wrt ice             [Pa]
    
-   C0            = 10.                                                  !--value assumed in Field2014            
+   REAL :: vitw                                                           !--vertical velocity                             [m/s]
+                                                                       
+   C0            = 10.                                                    !--value assumed in Field2014            
    rho_ice       = 950.
    sursat_iceext = -0.1
-   !capa_crystal  = 1. !r_ice                                       
    qzero(:)      = 0.
    cldfraliq(:)  = 0.
@@ -317,7 +318,7 @@
 
    sigma2_icefracturb(:) = 0.
-   mean_icefracturb(:)  = 0.
-
-   !--wrt liquid water
+   mean_icefracturb(:)   = 0.
+
+   !--wrt liquid 
    CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
    !--wrt ice
@@ -327,7 +328,6 @@
     DO i=1,klon
 
-
      rho_air  = pplay(i) / temp(i) / RD
-     !deltaz   = ( paprsdn(i) - paprsup(i) ) / RG / rho_air(i)
+
      ! because cldfra is intent in, but can be locally modified due to test 
      cldfra1D = cldfra(i)
@@ -364,31 +364,47 @@
            dicefracdT(i) = 0.
 
-        ! MPC temperature
+        !---------------------------------------------------------
+        !--             MIXED PHASE TEMPERATURE REGIME          
+        !--------------------------------------------------------- 
+        !--In the mixed phase regime (-38°C< T <0°C) we distinguish
+        !--3 possible subcases. 
+        !--1.  No pre-existing ice  
+        !--2A. Pre-existing ice and no turbulence
+        !--2B. Pre-existing ice and turbulence
         ELSE
-           ! Not enough TKE     
-           IF ( tke_dissip(i) .LE. eps )  THEN
-              qvap_cld(i)   = qsati(i) * cldfra1D
-              qliq(i)       = 0.
-              qice(i)       = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D   
-              cldfraliq(i)  = 0.
-              icefrac(i)    = 1.
-              dicefracdT(i) = 0.
+
+           vitw = -omega(i) / RG / rho_air
+           qiceini_incl  = qice_ini(i) / cldfra1D + snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) ) / cldfra1D
+
+           !--1. No preexisting ice : if vertical motion, fully liquid
+           !--cloud else fully iced cloud
+           IF ( qiceini_incl .LT. eps ) THEN
+              IF ( (vitw .GT. eps) .OR. (tke(i) .GT. eps) ) THEN
+                 qvap_cld(i)   = qsatl(i) * cldfra1D
+                 qliq(i)       = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D
+                 qice(i)       = 0.
+                 cldfraliq(i)  = 1.
+                 icefrac(i)    = 0.
+                 dicefracdT(i) = 0.
+              ELSE
+                 qvap_cld(i)   = qsati(i) * cldfra1D
+                 qliq(i)       = 0.
+                 qice(i)       = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D
+                 cldfraliq(i)  = 0.
+                 icefrac(i)    = 1.
+                 dicefracdT(i) = 0.
+              ENDIF
            
-           ! Enough TKE   
-           ELSE  
-              print*,"MOUCHOIRACTIVE"
-              !---------------------------------------------------------
-              !--               ICE SUPERSATURATION PDF   
-              !--------------------------------------------------------- 
-              !--If -38°C< T <0°C and there is enough turbulence, 
-              !--we compute the cloud liquid properties with a Gaussian PDF 
-              !--of ice supersaturation F(Si) (Field2014, Furtado2016). 
-              !--Parameters of the PDF are function of turbulence and 
-              !--microphysics/existing ice.
+
+           !--2. Pre-existing ice :computation of ice properties for
+           !--feedback
+           ELSE
+              ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. )
+
+              sursat_equ    = ai * vitw * tau_phaserelax
 
               sursat_iceliq = qsatl(i)/qsati(i) - 1.
               psati         = qsati(i) * pplay(i) / (RD/RV)
-
-              !-------------- MICROPHYSICAL TERMS --------------
+              
               !--We assume an exponential ice PSD whose parameters 
               !--are computed following Morrison&Gettelman 2008
@@ -399,16 +415,9 @@
               !--tau_phase_relax is the typical time of vapor deposition
               !--onto ice crystals
-              
-              qiceini_incl  = qice_ini(i) / cldfra1D
-              qsnowcld_incl = snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) ) / cldfra1D
-              sursat_env    = max(0., (qtot_incl(i) - qiceini_incl)/qsati(i) - 1.)
-              IF ( qiceini_incl .GT. eps ) THEN
-                nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )**3.33 * naero5**(0.0264*(RTT-temp(i))+0.0033)
-                lambda_PSD  = ( (RPI*rho_ice*nb_crystals) / (rho_air*(qiceini_incl + gamma_snwretro * qsnowcld_incl)) ) ** (1./3.)
-                N0_PSD      = nb_crystals * lambda_PSD
-                moment1_PSD = N0_PSD/lambda_PSD**2
-              ELSE
-                moment1_PSD = 0.
-              ENDIF
+
+              nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )**3.33 * naero5**(0.0264*(RTT-temp(i))+0.0033)
+              lambda_PSD  = ( (RPI*rho_ice*nb_crystals) / (rho_air * qiceini_incl ) ) ** (1./3.)
+              N0_PSD      = nb_crystals * lambda_PSD
+              moment1_PSD = N0_PSD/lambda_PSD**2
 
               !--Formulae for air thermal conductivity and water vapor diffusivity 
@@ -422,85 +431,97 @@
               B0 = 4. * RPI * capa_crystal * 1. / (  RLSTT**2 / air_thermal_conduct / RV / temp(i)**2  &
                                                   +  RV * temp(i) / psati / water_vapor_diff  )
-
-              invtau_phaserelax  = (bi * B0 * moment1_PSD )
-
-!             Old way of estimating moment1 : spherical crystals + monodisperse PSD              
-!             nb_crystals = rho_air * qiceini_incl / ( 4. / 3. * RPI * r_ice**3. * rho_ice )
-!             moment1_PSD = nb_crystals * r_ice 
-
-              !----------------- TURBULENT SOURCE/SINK TERMS -----------------
-              !--Tau_mixingenv is the time needed to homogeneize the parcel 
-              !--with its environment by turbulent diffusion over the parcel
-              !--length scale
-              !--if lmix_mpc <0, tau_mixigenv value is prescribed
-              !--else tau_mixigenv value is derived from tke_dissip and lmix_mpc
-              !--Tau_dissipturb is the time needed turbulence to decay due to
-              !--viscosity
-
+              tau_phaserelax = 1. / (bi * B0 * moment1_PSD )
+              
               ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. )
-              IF ( lmix_mpc .GT. 0 ) THEN 
-                 tau_mixingenv = ( lmix_mpc**2. / tke_dissip(i) )**(1./3.)
-              ELSE 
-                 tau_mixingenv = tau_mixenv
-              ENDIF 
-              
-              tau_dissipturb = gamma_taud * 2. * 2./3. * tke(i) / tke_dissip(i) / C0
-
-              !--------------------- PDF COMPUTATIONS ---------------------
-              !--Formulae for sigma2_pdf (variance), mean of PDF in Furtado2016
-              !--cloud liquid fraction and in-cloud liquid content are given 
-              !--by integrating resp. F(Si) and Si*F(Si)
-              !--Liquid is limited by the available water vapor trough a 
-              !--maximal liquid fraction
-
-              liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - qiceini_incl - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) )
-              sigma2_pdf = 1./2. * ( ai**2 ) *  2./3. * tke(i) * tau_dissipturb / ( invtau_phaserelax + 1./tau_mixingenv )
-              mean_pdf   = sursat_env * 1./tau_mixingenv / ( invtau_phaserelax + 1./tau_mixingenv )
-              cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - mean_pdf) / (SQRT(2.* sigma2_pdf) ) ) )
-              IF (cldfraliq(i) .GT. liqfra_max) THEN
-                  cldfraliq(i) = liqfra_max
-              ENDIF 
-              
-              qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.*RPI) * EXP( -1.*(sursat_iceliq - mean_pdf)**2. / (2.*sigma2_pdf) )  &
-                        - qsati(i) * cldfraliq(i) * (sursat_iceliq - mean_pdf )
-              
-              sigma2_icefracturb(i)= sigma2_pdf
-              mean_icefracturb(i)  = mean_pdf 
+
+              !--2A. No TKE : stationnary binary solution depending on omega
+              ! If Sequ > Siw liquid cloud, else ice cloud
+              IF ( tke_dissip(i) .LE. eps )  THEN
+                 IF (sursat_equ .GT. sursat_iceliq) THEN
+                    qvap_cld(i)   = qsatl(i) * cldfra1D
+                    qliq(i)       = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D
+                    qice(i)       = 0.
+                    cldfraliq(i)  = 1.
+                    icefrac(i)    = 0.
+                    dicefracdT(i) = 0.
+                 ELSE
+                    qvap_cld(i)   = qsati(i) * cldfra1D
+                    qliq(i)       = 0.
+                    qice(i)       = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D
+                    cldfraliq(i)  = 0.
+                    icefrac(i)    = 1.
+                    dicefracdT(i) = 0.
+                 ENDIF
+                 
+              !--2B. TKE and ice : ice supersaturation PDF 
+              !--we compute the cloud liquid properties with a Gaussian PDF 
+              !--of ice supersaturation F(Si) (Field2014, Furtado2016). 
+              !--Parameters of the PDF are function of turbulence and 
+              !--microphysics/existing ice.
+              ELSE  
+                      
+                 !--Tau_dissipturb is the time needed for turbulence to decay 
+                 !--due to viscosity
+                 tau_dissipturb = gamma_taud * 2. * 2./3. * tke(i) / tke_dissip(i) / C0
+
+                 !--------------------- PDF COMPUTATIONS ---------------------
+                 !--Formulae for sigma2_pdf (variance), mean of PDF in Raillard2025
+                 !--cloud liquid fraction and in-cloud liquid content are given 
+                 !--by integrating resp. F(Si) and Si*F(Si)
+                 !--Liquid is limited by the available water vapor trough a 
+                 !--maximal liquid fraction
+                 !--qice_ini(i) / cldfra1D = qiceincld without precip
+
+                 liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - (qice_ini(i) / cldfra1D) - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) )
+                 sigma2_pdf = 1./2. * ( ai**2 ) *  2./3. * tke(i) * tau_dissipturb * tau_phaserelax
+                 
+                 mean_pdf = ai * vitw * tau_phaserelax 
+                 
+                 cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - mean_pdf) / (SQRT(2.* sigma2_pdf) ) ) )
+                 IF (cldfraliq(i) .GT. liqfra_max) THEN
+                     cldfraliq(i) = liqfra_max
+                 ENDIF 
+                 
+                 qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.*RPI) * EXP( -1.*(sursat_iceliq - mean_pdf)**2. / (2.*sigma2_pdf) )  &
+                           - qsati(i) * cldfraliq(i) * (sursat_iceliq - mean_pdf )
+                 
+                 sigma2_icefracturb(i)= sigma2_pdf
+                 mean_icefracturb(i)  = mean_pdf 
      
-              !------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION  ------------
-
-              IF ( (qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps) ) THEN
-                  qliq_incl    = 0.
-                  cldfraliq(i) = 0.
-              END IF
-               
-              !--Specific humidity is the max between qsati and the weighted mean between 
-              !--qv in MPC patches and qv in ice-only parts. We assume that MPC parts are
-              !--always at qsatl and ice-only parts slightly subsaturated (qsati*sursat_iceext+1)
-              !--The whole cloud can therefore be supersaturated but never subsaturated.
-
-              qvap_incl = MAX(qsati(i), ( 1. - cldfraliq(i) ) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i) )
-
-
-              IF ( qvap_incl  .GE. qtot_incl(i) ) THEN
-                 qvap_incl = qsati(i)
-                 qliq_incl = qtot_incl(i) - qvap_incl
-                 qice_incl = 0.
-
-              ELSEIF ( (qvap_incl + qliq_incl) .GE. qtot_incl(i) ) THEN
-                 qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)
-                 qice_incl = 0.
-              ELSE
-                 qice_incl = qtot_incl(i) - qvap_incl - qliq_incl
-              END IF
-
-              qvap_cld(i)   = qvap_incl * cldfra1D
-              qliq(i)       = qliq_incl * cldfra1D
-              qice(i)       = qice_incl * cldfra1D
-              icefrac(i)    = qice(i) / ( qice(i) + qliq(i) )
-              dicefracdT(i) = 0.
-              !print*,'MPC turb'
-
-           END IF ! Enough TKE
+                 !------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION  ------------
+
+                 IF ( (qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps) ) THEN
+                     qliq_incl    = 0.
+                     cldfraliq(i) = 0.
+                 END IF
+                  
+                 !--Specific humidity is the max between qsati and the weighted mean between 
+                 !--qv in MPC patches and qv in ice-only parts. We assume that MPC parts are
+                 !--always at qsatl and ice-only parts slightly subsaturated (qsati*sursat_iceext+1)
+                 !--The whole cloud can therefore be supersaturated but never subsaturated.
+
+                 qvap_incl = MAX(qsati(i), ( 1. - cldfraliq(i) ) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i) )
+
+                 IF ( qvap_incl  .GE. qtot_incl(i) ) THEN
+                    qvap_incl = qsati(i)
+                    qliq_incl = qtot_incl(i) - qvap_incl
+                    qice_incl = 0.
+
+                 ELSEIF ( (qvap_incl + qliq_incl) .GE. qtot_incl(i) ) THEN
+                    qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)
+                    qice_incl = 0.
+                 ELSE
+                    qice_incl = qtot_incl(i) - qvap_incl - qliq_incl
+                 END IF
+
+                 qvap_cld(i)   = qvap_incl * cldfra1D
+                 qliq(i)       = qliq_incl * cldfra1D
+                 qice(i)       = qice_incl * cldfra1D
+                 icefrac(i)    = qice(i) / ( qice(i) + qliq(i) )
+                 dicefracdT(i) = 0.
+
+              END IF ! Enough TKE
+           
+           END IF ! End qini
 
         END IF ! ! MPC temperature
@@ -510,4 +531,5 @@
    ENDDO ! klon
 END SUBROUTINE ICEFRAC_LSCP_TURB
+!
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
Index: LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5380)
+++ LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5383)
@@ -1540,4 +1540,6 @@
   TYPE(ctrl_out), SAVE :: o_ocond = ctrl_out((/ 2, 3, 4, 10, 10, 10, 11, 11, 11, 11/), &
     'ocond', 'Condensed water', 'kg/kg', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_oice = ctrl_out((/ 2, 3, 4, 10, 10, 10, 11, 11, 11, 11/), &
+    'oice', 'Ice water', 'kg/kg', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_qbs = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'qbs', 'Specific content of blowing snow', 'kg/kg', (/ ('', i=1, 10) /))
Index: LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 5380)
+++ LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 5383)
@@ -138,5 +138,5 @@
          o_cldicemxrat, o_cldwatmxrat, o_reffclwtop, o_ec550aer, &
          o_lwcon, o_iwcon, o_temp, o_theta, &
-         o_ovapinit, o_ovap, o_oliq, o_ocond, o_geop,o_qbs, &
+         o_ovapinit, o_ovap, o_oliq, o_ocond, o_oice, o_geop,o_qbs, &
          o_vitu, o_vitv, o_vitw, o_pres, o_paprs, &
          o_zfull, o_zhalf, o_rneb, o_rnebjn, o_rnebcon, &
@@ -1997,4 +1997,7 @@
        IF (vars_defined) zx_tmp_fi3d = ql_seri+qs_seri
        CALL histwrite_phy(o_ocond, zx_tmp_fi3d)
+      
+       IF (vars_defined) zx_tmp_fi3d = qs_seri
+       CALL histwrite_phy(o_oice, zx_tmp_fi3d)
 
        CALL histwrite_phy(o_geop, zphi)
Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 5380)
+++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 5383)
@@ -1243,7 +1243,4 @@
     !albedo SB <<<
 
-    !--Lea Raillard qs_ini
-    REAL, dimension(klon,klev) :: qs_ini
-
     !--OB variables for mass fixer (hard coded for now)
     REAL qql1(klon),qql2(klon),corrqql
@@ -2434,6 +2431,4 @@
        ENDDO
     ENDDO
-    ! Lea Raillard qs_ini for cloud phase param.
-    qs_ini(:,:)=qs_seri(:,:)
     !
     !--OB water mass fixer
@@ -2654,4 +2649,14 @@
             ratio_qi_qtot(i,k) = 0.
             rvc_seri(i,k) = 0.
+          ENDIF
+        ENDDO
+      ENDDO
+    ELSE
+      DO k = 1, klev
+        DO i = 1, klon
+          IF ( ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) ) .GT. 0. ) THEN
+            ratio_qi_qtot(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
+          ELSE
+            ratio_qi_qtot(i,k) = 0.
           ENDIF
         ENDDO
@@ -3832,6 +3837,6 @@
     !ENDIF
 
-    CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay, &
-         t_seri, q_seri,qs_ini,ptconv,ratqs,sigma_qtherm, &
+    CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay,omega, &
+         t_seri, q_seri, ptconv, ratqs, sigma_qtherm, &
          d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, &
          pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb,  &
Index: LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 5380)
+++ LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 5383)
@@ -5065,7 +5065,7 @@
     !ENDIF
 
-    CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay, &
-         t_seri, q_seri,qs_ini,ptconv,ratqs,sigma_qtherm, &
-         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, & 
+    CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay,omega, &
+         t_seri, q_seri, ptconv, ratqs, sigma_qtherm, &
+         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, &
          pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb,  &
          radocond, picefra, rain_lsc, snow_lsc, &
@@ -5073,5 +5073,5 @@
          prfl, psfl, rhcl,  &
          zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
-         iflag_ice_thermo, distcltop, temp_cltop,  &
+         iflag_ice_thermo, distcltop, temp_cltop,   &
          pbl_tke(:,:,is_ave), pbl_eps(:,:,is_ave), &
          cell_area, &
