Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90	(revision 5633)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90	(revision 5634)
@@ -174,5 +174,5 @@
       pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
       pdf_x = qcritcont(i) / qsatl(i) * 100.
-      pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+      pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
       IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
         pdf_fra_above_qcritcont = 0.
@@ -191,5 +191,5 @@
     
       pdf_x = qsat(i) / qsatl(i) * 100.
-      pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+      pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
       IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
         pdf_fra_above_qsat = 0.
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90	(revision 5633)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90	(revision 5634)
@@ -556,5 +556,5 @@
               !--the old one, it is comprised between 0 and 1
               tauinv_depsub = qiceincld**(1./3.) * kappa_depsub
-              qice_ratio = tauinv_depsub * dtime / 1.5 / qiceinmix * ( qsat(i) - qvapincld )
+              qice_ratio = tauinv_depsub * dtime / 1.5 / qiceincld * ( qsat(i) - qvapincld )
               !--The new vapor in the cloud is increased with the
               !--sublimated ice
@@ -683,6 +683,6 @@
         pdf_gamma(i) = GAMMA(1. + 1. / pdf_alpha(i))
         !--Barrier to avoid overflows
-        pdf_scale(i) = MAX(0.01, pdf_std / SQRT( &
-                                    GAMMA(1. + 2. / pdf_alpha(i)) - pdf_gamma(i)**2 ))
+        pdf_scale(i) = MAX(eps, MIN(rhl_clr / pdf_gamma(i), pdf_std / SQRT( &
+                                    GAMMA(1. + 2. / pdf_alpha(i)) - pdf_gamma(i)**2 )))
         pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
 
@@ -692,5 +692,5 @@
 
         pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
-        pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+        pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
         IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
           pdf_fra_above_nuc = 0.
@@ -747,5 +747,5 @@
         pdf_alpha(i) = 1.
         pdf_gamma(i) = 1.
-        pdf_scale(i) = 0.01
+        pdf_scale(i) = eps
       ENDIF ! clrfra(i) .GT. eps
 
@@ -786,5 +786,6 @@
         !--are spherical.
         !-- b / a = bovera = MAX(0.1, cldfra)
-        bovera = MAX(0.1, cldfra(i))
+        !bovera = MAX(0.1, cldfra(i))
+        bovera = 0.111111111
         !--P / a is a function of b / a only, that we can calculate
         !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
@@ -877,7 +878,8 @@
 
           rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
+          pdf_x = qvapinclr_lim / qsatl(i) * 100.
           pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
-          pdf_x = qvapinclr_lim / qsatl(i) * 100.
-          pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+          pdf_x = qsat(i) / qsatl(i) * 100.
+          pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
           IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
             pdf_fra_above_lim = 0.
@@ -1034,7 +1036,8 @@
 
           rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
+          pdf_x = qvapinclr_lim / qsatl(i) * 100.
           pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
-          pdf_x = qvapinclr_lim / qsatl(i) * 100.
-          pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+          pdf_x = qsat(i) / qsatl(i) * 100.
+          pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
           IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
             pdf_fra_above_lim = 0.
@@ -1138,7 +1141,8 @@
 
           rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
+          pdf_x = qvapinclr_lim / qsatl(i) * 100.
           pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
-          pdf_x = qvapinclr_lim / qsatl(i) * 100.
-          pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+          pdf_x = qsat(i) / qsatl(i) * 100.
+          pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
           IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
             pdf_fra_above_lim = 0.
@@ -1193,5 +1197,5 @@
         pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
         pdf_x = qsat(i) / qsatl(i) * 100.
-        pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+        pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
         IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
           issrfra(i) = 0.
@@ -1206,4 +1210,8 @@
           issrfra(i) = EXP( - pdf_y ) * clrfra(i)
           qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
+        ENDIF
+        IF ( issrfra(i) .LT. eps ) THEN
+          issrfra(i) = 0.
+          qissr(i) = 0.
         ENDIF
       ELSE