Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.F90	(revision 5169)
+++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.F90	(revision 5170)
@@ -243,5 +243,6 @@
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   ! Compute the liquid, ice and vapour content (+ice fraction) based 
-  ! on turbulence (see Fields 2014, Furtado 2016)
+  ! on turbulence (see Fields 2014, Furtado 2016, Raillard 2025)
+  ! L.Raillard (30/08/24)
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
@@ -375,5 +376,6 @@
            
            ! Enough TKE   
-           ELSE   
+           ELSE  
+              print*,"MOUCHOIRACTIVE"
               !---------------------------------------------------------
               !--               ICE SUPERSATURATION PDF   
@@ -403,7 +405,7 @@
               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*24.) / (6.*(qiceini_incl + gamma_snwretro * qsnowcld_incl)) ) ** (1./3.)
+                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/2./lambda_PSD**2
+                moment1_PSD = N0_PSD/lambda_PSD**2
               ELSE
                 moment1_PSD = 0.
@@ -456,7 +458,7 @@
               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 
+              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) )  &
@@ -464,6 +466,7 @@
               
               sigma2_icefracturb(i)= sigma2_pdf
-              mean_icefracturb(i)  = mean_pdf      
-              !------------ ICE AMOUNT AND WATER CONSERVATION  ------------
+              mean_icefracturb(i)  = mean_pdf 
+     
+              !------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION  ------------
 
               IF ( (qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps) ) THEN
@@ -471,10 +474,13 @@
                   cldfraliq(i) = 0.
               END IF
-              
-              !--Choice for in-cloud vapor : 
-              !--1.Weighted mean between qvap in MPC parts and in ice-only parts
-              !--2.Always at ice saturation
+               
+              !--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)