Index: /LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90	(revision 5209)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90	(revision 5210)
@@ -975,6 +975,6 @@
                           
                    CALL condensation_lognormal( &
-                       klon, Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), &
-                       keepgoing(:), rneb(:,k), zqn(:), qvc(:))
+                       klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
+                       keepgoing, rneb(:,k), zqn, qvc)
 
 
Index: /LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.F90	(revision 5209)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.F90	(revision 5210)
@@ -124,5 +124,5 @@
                                 mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, &
                                 rhlmid_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, rhl0_pdf_lscp, &
-                                cond_thresh_pdf_lscp, coef_mixing_lscp, coef_shear_lscp, & 
+                                coef_mixing_lscp, coef_shear_lscp, & 
                                 chi_mixing_lscp, rho_ice
 
@@ -212,5 +212,9 @@
 REAL :: pres_sat, kappa
 REAL :: air_thermal_conduct, water_vapor_diff
-REAL :: r_ice
+REAL :: iwc
+REAL :: iwc_log_inf100, iwc_log_sup100
+REAL :: iwc_inf100, alpha_inf100, coef_inf100
+REAL :: mu_sup100, sigma_sup100, coef_sup100
+REAL :: Dm_ice, rm_ice
 !
 ! for sublimation
@@ -404,12 +408,12 @@
         !--during deposition, and alpha = 1. during sublimation.
         !--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
-        !-- C = capa_cond_cirrus * r_ice
+        !-- C = capa_cond_cirrus * rm_ice
         !
         !--We have qice = Nice * mi, where Nice is the ice crystal
         !--number concentration per kg of moist air
         !--HYPOTHESIS 1: the ice crystals are spherical, therefore
-        !-- mi = 4/3 * pi * r_ice**3 * rho_ice
+        !-- mi = 4/3 * pi * rm_ice**3 * rho_ice
         !--HYPOTHESIS 2: the ice crystals are monodisperse with the
-        !--initial radius r_ice_0.
+        !--initial radius rm_ice_0.
         !--NB. this is notably different than the assumption
         !--of a distributed qice in the cloud made in the sublimation process;
@@ -417,20 +421,20 @@
         !
         !--As the deposition process does not create new ice crystals,
-        !--and because we assume a same r_ice value for all crystals
+        !--and because we assume a same rm_ice value for all crystals
         !--therefore the sublimation process does not destroy ice crystals
         !--(or, in a limit case, it destroys all ice crystals), then
         !--Nice is a constant during the sublimation/deposition process.
-        !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI r_ice_0**3 rho_ice )
+        !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
         !
         !--The deposition equation then reads:
-        !-- dqi/dt = alpha*4pi*capa_cond_cirrus*r_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice
-        !-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *r_ice_0*(qvc-qsat)/qsat &
+        !-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice
+        !-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
         !--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
-        !--                         * qi_0 / ( 4/3 RPI r_ice_0**3 rho_ice )
+        !--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
         !-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
-        !--  *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 r_ice_0**2 rho_ice )
+        !--  *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 rm_ice_0**2 rho_ice )
         !--and we have
-        !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / r_ice_0**2
-        !-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / r_ice_0**2
+        !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
+        !-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
         !--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
         !
@@ -489,7 +493,22 @@
             !--the new vapor qvapincld_new
 
-            !--r_ice formula from Sun and Rikus (1999)
-            r_ice = 1.e-6 * ( 45.8966 * ( qiceincld * rho * 1e3 )**0.2214 &
-                  + 0.7957 * ( qiceincld * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2.
+            !--rm_ice formula from McFarquhar and Heymsfield (1997)
+            iwc = qiceincld * rho * 1e3
+            iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837)
+            iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) )
+            iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) )
+
+            alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100
+            coef_inf100 = iwc_inf100 * alpha_inf100**3. / 120.
+
+            mu_sup100 = ( 5.2 + 0.0013 * ( temp(i) - RTT ) ) &
+                      + ( 0.026 - 1.2E-3 * ( temp(i) - RTT ) ) * iwc_log_sup100
+            sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp(i) - RTT ) ) &
+                         + ( 0.018 - 2.1E-4 * ( temp(i) - RTT ) ) * iwc_log_sup100
+            coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3 * mu_sup100 + 4.5 * sigma_sup100**2. )
+
+            Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2. ) &
+                   * coef_sup100 ) / ( coef_inf100 + coef_sup100 )
+            rm_ice = Dm_ice / 2. * 1.E-6
 
             IF ( qvapincld .GE. qsat(i) ) THEN
@@ -497,5 +516,5 @@
               !--Exact explicit integration (qvc exact, qice explicit)
               qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) &
-                            * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / r_ice**2. )
+                            * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2. )
             ELSE
               !--If the cloud is initially subsaturated
@@ -505,5 +524,5 @@
               !--qice_ratio is the ratio between the new ice content and
               !--the old one, it is comprised between 0 and 1
-              qice_ratio = ( 1. - 2. / 3. / kappa / r_ice**2. * dtime * ( qsat(i) - qvapincld ) )
+              qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2. * dtime * ( qsat(i) - qvapincld ) )
 
               IF ( qice_ratio .LT. 0. ) THEN
@@ -538,4 +557,19 @@
             qvapincld_new = qsat(i)
           ENDIF ! ok_unadjusted_clouds
+          
+          !--Adjustment of the IWC to the new vapor in cloud
+          !--(this can be either positive or negative)
+          dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )
+          dqi_adj(i) = - dqvc_adj(i)
+
+          !--Add tendencies
+          !--The vapor in the cloud is updated, but not qcld as it is constant
+          !--through this process, as well as cldfra which is unmodified
+          qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))
+          
+
+          !------------------------------------
+          !--    DISSIPATION OF THE CLOUD    --
+          !------------------------------------
 
           !--If the vapor in cloud is below vapor needed for the cloud to survive
@@ -577,6 +611,5 @@
 
             !--Tendencies and diagnostics
-            dqvc_sub(i) = dcf_sub(i) * qvapincld
-            dqi_sub(i) = dqt_sub - dqvc_sub(i)
+            dqvc_sub(i) = dqt_sub
 
             !--Add tendencies
@@ -586,14 +619,4 @@
 
           ENDIF ! qvapincld .LT. qvapincld_new
-
-          !--Adjustment of the IWC of the remaining cloud
-          !--to the new vapor in cloud (this can be either positive or negative)
-          dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )
-          dqi_adj(i) = - dqvc_adj(i)
-
-          !--Add tendencies
-          !--The vapor in the cloud is updated, but not qcld as it is constant
-          !--through this process, as well as cldfra which is unmodified
-          qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))
 
         ENDIF   ! qiceincld .LT. eps
@@ -653,31 +676,15 @@
 
         IF ( ok_warm_cloud ) THEN
-          !--If the statistical scheme is activated, the first
-          !--calculation of the tendencies is kept for the diagnosis of
-          !--the cloud properties, and the condensation process stops
+          !--If the statistical scheme is activated, the calculated increase is equal
+          !--to the cloud fraction, we assume saturation adjustment, and the
+          !--condensation process stops
           cldfra(i) = cf_cond
           qcld(i) = qt_cond
           qvc(i) = cldfra(i) * qsat(i)
 
-        ELSEIF ( cf_cond .GT. cond_thresh_pdf_lscp ) THEN
-          !--We check that there is something to condense, if not the
-          !--condensation process stops
-
-          !--Calculation of the limit qclr in cloud value (stored in pdf_e4)
-
-          pdf_e3 = ( LOG( 1. / cond_thresh_pdf_lscp ) ** ( 1. / pdf_alpha ) - pdf_e2 ) &
-                 / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. )
-          pdf_e4 = ( 2. * pdf_e3 * pdf_e1 - pdf_x ) &
-                 / ( pdf_e3 * pdf_e1 / rhlmid_pdf_lscp - 1. )
-          IF ( ( MIN(pdf_e4, rhl_clr) .LT. rhlmid_pdf_lscp ) .OR. ( pdf_e4 .GT. rhl_clr ) ) THEN
-            pdf_e4 = pdf_x / ( pdf_e3 * pdf_e1 / rhlmid_pdf_lscp + 1. )
-          ENDIF
-          pdf_e4 = pdf_e4 * qsatl(i) / 100.
-
-          !--Calculation of the final tendencies
-          dcf_con(i) = ( 1. - cldfra(i) ) * ( pdf_e4 - qclr / ( 1. - cldfra(i) ) ) &
-                     / ( pdf_e4 - qt_cond / cf_cond )
-          dqt_con = qclr - pdf_e4 * ( 1. - cldfra(i) - dcf_con(i) )
-
+        ELSEIF ( cf_cond .GT. eps ) THEN
+
+          dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond
+          dqt_con = ( 1. - cldfra(i) ) * qt_cond
           
           !--Barriers
@@ -692,12 +699,28 @@
             qvapincld = gamma_cond(i) * qsat(i)
             qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
-            !--r_ice formula from Sun and Rikus (1999)
-            r_ice = 1.e-6 * ( 45.8966 * ( qiceincld * rho * 1e3 )**0.2214 &
-                  + 0.7957 * ( qiceincld * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2.
+
+            !--rm_ice formula from McFarquhar and Heymsfield (1997)
+            iwc = qiceincld * rho * 1e3
+            iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837)
+            iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) )
+            iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) )
+
+            alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100
+            coef_inf100 = iwc_inf100 * alpha_inf100**3. / 120.
+
+            mu_sup100 = ( 5.2 + 0.0013 * ( temp(i) - RTT ) ) &
+                      + ( 0.026 - 1.2E-3 * ( temp(i) - RTT ) ) * iwc_log_sup100
+            sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp(i) - RTT ) ) &
+                         + ( 0.018 - 2.1E-4 * ( temp(i) - RTT ) ) * iwc_log_sup100
+            coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3. * mu_sup100 + 4.5 * sigma_sup100**2. )
+
+            Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2. ) &
+                   * coef_sup100 ) / ( coef_inf100 + coef_sup100 )
+            rm_ice = Dm_ice / 2. * 1.E-6
             !--As qvapincld is necessarily greater than qsat, we only
             !--use the exact explicit formulation
             !--Exact explicit version
             qvapincld = qsat(i) + ( qvapincld - qsat(i) ) &
-                      * EXP( - depo_coef_cirrus * dtime / 2. * qiceincld / kappa / r_ice**2. )
+                      * EXP( - depo_coef_cirrus * dtime / 2. * qiceincld / kappa / rm_ice**2. )
           ELSE
             !--We keep the saturation adjustment hypothesis, and the vapor in the
@@ -715,5 +738,5 @@
           !-- qtot(i) - dqi_con(i) - dqvc_con(i)
 
-        ENDIF ! ok_warm_cloud, cf_cond .GT. cond_thresh_pdf_lscp
+        ENDIF ! ok_warm_cloud, cf_cond .GT. eps
       ENDIF ! ( 1. - cldfra(i) ) .GT. eps
 
@@ -854,5 +877,6 @@
         L_shear = coef_shear_lscp * shear(i) * dz * dtime
         !--therefore, the total increase in fraction is
-        shear_fra = RPI * L_shear * a_mix * bovera * N_cld_mix &
+        !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
+        shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix &
                   / cell_area(i) / ( 1. - dcf_con(i) )
         !--and the environment and cloud mixed fractions are the same,
Index: /LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90	(revision 5209)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90	(revision 5210)
@@ -185,7 +185,4 @@
   REAL, SAVE, PROTECTED :: rhl0_pdf_lscp=88.7                ! [%] factor for the PDF fit of water vapor in UTLS
   !$OMP THREADPRIVATE(rhl0_pdf_lscp)
-  
-  REAL, SAVE, PROTECTED :: cond_thresh_pdf_lscp=1.E-10       ! [-] threshold for the formation of new cloud
-  !$OMP THREADPRIVATE(cond_thresh_pdf_lscp)
   
   REAL, SAVE, PROTECTED :: a_homofreez=2.349                 ! [-] factor for the Koop homogeneous freezing fit
@@ -433,5 +430,4 @@
     CALL getin_p('kappa_pdf_lscp',kappa_pdf_lscp)
     CALL getin_p('rhl0_pdf_lscp',rhl0_pdf_lscp)
-    CALL getin_p('cond_thresh_pdf_lscp',cond_thresh_pdf_lscp)
     CALL getin_p('a_homofreez',a_homofreez)
     CALL getin_p('b_homofreez',b_homofreez)
