Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90	(revision 4830)
+++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90	(revision 4832)
@@ -136,5 +136,4 @@
     !--multiplying by dz = - dP / g / rho (line 3-4)
     !--formula from Sundqvist 1988, Klemp & Wilhemson 1978
-
     draineva = - precipfracclr(i) * coef_eva * (1. - qvap(i) / qsatl(i)) &
              * ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva &
@@ -197,6 +196,7 @@
 
     !--Add tendencies
-    rainclr(i) = rainclr(i) + draineva
-    snowclr(i) = snowclr(i) + dsnowsub
+
+    rainclr(i) = MAX(0., rainclr(i) + draineva)
+    snowclr(i) = MAX(0., snowclr(i) + dsnowsub)
 
     !--If there is no more precip fluxes, the precipitation fraction in clear
@@ -248,5 +248,5 @@
                           cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, & 
                           thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
-                          rho_rain, rho_snow, r_rain, r_snow,                  &
+                          rho_rain, rho_snow, r_rain, r_snow, rho_ice, r_ice,  &
                           tau_auto_snow_min, tau_auto_snow_max,                &
                           thresh_precip_frac, eps, air_thermal_conduct,        &
@@ -340,5 +340,5 @@
 REAL :: dqrfreez_max
 REAL :: tau_freez
-REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez
+REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2
 REAL :: coef_freez
 REAL :: dqifreez        !--loss of ice cloud content due to collection of ice from rain [kg/kg/s]
@@ -371,5 +371,5 @@
   dqifreez= 0.
 
-  !--hum_to_flux: coef to convert a specific quantity to a flux 
+  !--hum_to_flux: coef to convert a delta in specific quantity to a flux 
   !-- hum_to_flux = rho * dz/dt = 1 / g * dP/dt
   hum_to_flux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime
@@ -398,8 +398,16 @@
   !--evaporated (see barrier at the end of poprecip_precld)
   !--NB. precipfraccld(i) is here the cloud fraction of the layer above
-  precipfractot = 1. - ( 1. - precipfractot ) * &
+  !precipfractot = 1. - ( 1. - precipfractot ) * &
+  !               ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
+  !             / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
+
+
+  IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN
+    precipfractot = 1.
+  ELSE
+    precipfractot = 1. - ( 1. - precipfractot ) * &
                  ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
-               / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
-
+               / ( 1. - precipfraccld(i) )
+  ENDIF
 
   !--precipfraccld(i) is here the cloud fraction of the layer above
@@ -451,10 +459,10 @@
   precipfraccld(i) = precipfraccld(i) + dprecipfraccld
   precipfracclr(i) = precipfracclr(i) + dprecipfracclr
-  rainclr(i) = rainclr(i) + drainclr
-  snowclr(i) = snowclr(i) + dsnowclr
-  raincld(i) = raincld(i) - drainclr
-  snowcld(i) = snowcld(i) - dsnowclr
+
+  rainclr(i) = MAX( 0. , rainclr(i) + drainclr )
+  snowclr(i) = MAX( 0. , snowclr(i) + dsnowclr )
+  raincld(i) = MAX( 0. , raincld(i) - drainclr )
+  snowcld(i) = MAX( 0. , snowcld(i) - dsnowclr )
   
-
   !--If vertical heterogeneity is taken into account, we use
   !--the "true" volume fraction instead of a modified 
@@ -513,5 +521,5 @@
       !--Add tendencies
       qliq(i) = qliq(i) + dqlcol
-      raincld(i) = raincld(i) - dqlcol * hum_to_flux(i)
+      raincld(i) =  MAX( 0. , raincld(i) - dqlcol * hum_to_flux(i) )
 
       !--Diagnostic tendencies
@@ -535,5 +543,5 @@
       !--Add tendencies
       qice(i) = qice(i) + dqiagg
-      snowcld(i) = snowcld(i) - dqiagg * hum_to_flux(i)
+      snowcld(i) = MAX( 0. , snowcld(i) - dqiagg * hum_to_flux(i) )
 
       !--Diagnostic tendencies
@@ -599,6 +607,6 @@
     qliq(i) = qliq(i) + dqlauto
     qice(i) = qice(i) + dqiauto
-    raincld(i) = raincld(i) - dqlauto * hum_to_flux(i)
-    snowcld(i) = snowcld(i) - dqiauto * hum_to_flux(i)
+    raincld(i) = MAX( 0. , raincld(i) - dqlauto * hum_to_flux(i) )
+    snowcld(i) = MAX( 0. , snowcld(i) - dqiauto * hum_to_flux(i) )
 
     !--Diagnostic tendencies
@@ -631,5 +639,5 @@
       !--Add tendencies
       qliq(i) = qliq(i) + dqlrim
-      snowcld(i) = snowcld(i) - dqlrim * hum_to_flux(i)
+      snowcld(i) = MAX( 0. , snowcld(i) - dqlrim * hum_to_flux(i) )
 
       !--Temperature adjustment with the release of latent
@@ -699,5 +707,5 @@
       !--Barrier to limit the melting flux to the clr snow flux in the mesh
       dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / hum_to_flux(i))
-   ENDIF
+    ENDIF
 
     !--In cloudy air
@@ -726,11 +734,14 @@
       dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt
       dqstotmelt = dqsmelt_max
+
     ENDIF
 
     !--Add tendencies
-    rainclr(i) = rainclr(i) - dqsclrmelt * hum_to_flux(i)
-    raincld(i) = raincld(i) - dqscldmelt * hum_to_flux(i)
-    snowclr(i) = snowclr(i) + dqsclrmelt * hum_to_flux(i)
-    snowcld(i) = snowcld(i) + dqscldmelt * hum_to_flux(i)
+
+    rainclr(i) = MAX( 0. , rainclr(i) - dqsclrmelt * hum_to_flux(i) )
+    raincld(i) = MAX( 0. , raincld(i) - dqscldmelt * hum_to_flux(i) )
+    snowclr(i) = MAX( 0. , snowclr(i) + dqsclrmelt * hum_to_flux(i) )
+    snowcld(i) = MAX( 0. , snowcld(i) + dqscldmelt * hum_to_flux(i) )
+
 
     !--Temperature adjustment with the release of latent
@@ -749,8 +760,10 @@
   !--                  FREEZING
   !---------------------------------------------------------
-  !--Process through which rain freezes into snow. This is
-  !--parameterized as an exponential decrease of the rain
-  !--water content.
-  !--The formula is homemade.
+  !--Process through which rain freezes into snow. 
+  !-- We parameterize it as a 2 step process:
+  !-- first: freezing following collision with ice crystals
+  !-- second: immersion freezing following (inspired by Bigg 1953)
+  !-- the latter is parameterized as an exponential decrease of the rain
+  !-- water content with a homemade formulya
   !--This is based on a caracteritic time of freezing, which
   !--exponentially depends on temperature so that it is
@@ -758,5 +771,4 @@
   !--0 for RTT (=273.15 K).
   !--NB.: this process needs a temperature adjustment
-
   !--dqrfreez_max: maximum rain freezing so that temperature
   !--              stays lower than 273 K [kg/kg]
@@ -768,30 +780,25 @@
   IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
 
-    !--Computed according to
-    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
-    dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
+      
+    !--1st step: freezing following collision with ice crystals
+    !--Sub-process of freezing which quantifies the collision between
+    !--ice crystals in suspension and falling rain droplets.
+    !--The rain droplets freeze, becoming graupel, and carrying
+    !--the ice crystal (which acted as an ice nucleating particle).
+    !--The formula is adapted from the riming formula.
+    !-- it works only in the cloudy part
+
+    dqifreez = 0.
+    dqrtotfreez_step1 = 0.
+
+    IF ((qice(i) .GT. 0.) .AND. (raincld(i) .GT. 0.)) THEN
+      dqrclrfreez = 0.
+      dqrcldfreez = 0.
+
+      !--Computed according to
+      !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
+      dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
                          * ( 1. + RVTMP2 * qtot(i) ))
  
-    tau_freez = 1. / ( beta_freez &
-              * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
-
-    !--Initialisation 
-    dqrclrfreez = 0.
-    dqrcldfreez = 0.
-
-    !--In clear air
-    IF ( rainclr(i) .GT. 0. ) THEN
-      !--Exact solution of dqrain/dt = -qrain/tau_freez
-      dqrclrfreez = rainclr(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
-    ENDIF
-
-    !--In cloudy air
-    IF ( raincld(i) .GT. 0. ) THEN
-      !--Exact solution of dqrain/dt = -qrain/tau_freez
-      dqrcldfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
-
-      !---------------------------------------------------------
-      !--    FREEZING OF RAIN WITH CONTACT OF ICE CRYSTALS
-      !---------------------------------------------------------
       !--Sub-process of freezing which quantifies the collision between
       !--ice crystals in suspension and falling rain droplets.
@@ -800,43 +807,86 @@
       !--The formula is adapted from the riming formula.
 
-      !--The sticking efficacy is perfect.
+      !--The collision efficiency is assumed unity 
       Eff_rain_ice = 1.
       coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice
-      !-- ATTENTION Double implicit version? + barriers if needed
       !--Exact version, which does not need a barrier because of
-      !--the exponential decrease
+      !--the exponential decrease.
       dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. )
 
-      !--We add the two freezing processes in cloud
-      dqrcldfreez = dqrcldfreez - dqifreez
-    ENDIF
-
-    !--Barriers
+      !--We add the part of rain water that freezes, limited by a temperature barrier
+      !-- this quantity is calculated assuming that the number of drop that freeze correspond to the number
+      !-- of crystals collected (and assuming uniform distributions of ice crystals and rain drops)
+      dqrcldfreez = MAX(dqifreez*rho_rain*r_rain/(rho_ice*r_ice),-raincld(i)/hum_to_flux(i))
+      dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max) 
+      dqrtotfreez_step1 = dqrcldfreez
+
+      !--Add tendencies
+      qice(i) = qice(i) + dqifreez
+      raincld(i) = MAX( 0. , raincld(i) + dqrcldfreez * hum_to_flux(i) )
+      snowcld(i) = MAX( 0. , snowcld(i) - dqrcldfreez * hum_to_flux(i) - dqifreez * hum_to_flux(i) )
+      temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD &
+                      / ( 1. + RVTMP2 * qtot(i) )
+
+    ENDIF
+    
+    !-- Second step immersion freezing of rain drops
+    !-- with a homemade timeconstant depending on temperature
+
+    dqrclrfreez = 0.
+    dqrcldfreez = 0.
+    dqrtotfreez_step2 = 0.
+    !--Computed according to
+    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
+
+    dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
+                         * ( 1. + RVTMP2 * qtot(i) ))
+
+    
+    tau_freez = 1. / ( beta_freez &
+              * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
+
+
+    !--In clear air
+    IF ( rainclr(i) .GT. 0. ) THEN
+      !--Exact solution of dqrain/dt = -qrain/tau_freez
+      dqrclrfreez = rainclr(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
+    ENDIF
+
+    !--In cloudy air
+    IF ( raincld(i) .GT. 0. ) THEN
+      !--Exact solution of dqrain/dt = -qrain/tau_freez
+      dqrcldfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
+    ENDIF
+
+    !--temperature barrie step 2
     !--It is activated if the total is LOWER than the max
     !--because everything is negative
-    dqrtotfreez = dqrclrfreez + dqrcldfreez
-    IF ( dqrtotfreez .LT. dqrfreez_max ) THEN
+    dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez
+  
+    IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN
       !--We redistribute the max freezed rain keeping
       !--the clear/cloud partition of the freezing rain
-      dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez
-      dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez
-      dqrtotfreez = dqrfreez_max
+      dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2
+      dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2
+      dqrtotfreez_step2 = dqrfreez_max
     ENDIF
 
 
     !--Add tendencies
-    rainclr(i) = rainclr(i) + dqrclrfreez * hum_to_flux(i)
-    raincld(i) = raincld(i) + dqrcldfreez * hum_to_flux(i)
-    snowclr(i) = snowclr(i) - dqrclrfreez * hum_to_flux(i)
-    snowcld(i) = snowcld(i) - dqrcldfreez * hum_to_flux(i)
+    rainclr(i) = MAX( 0. , rainclr(i) + dqrclrfreez * hum_to_flux(i) )
+    raincld(i) = MAX( 0. , raincld(i) + dqrcldfreez * hum_to_flux(i) )
+    snowclr(i) = MAX( 0. , snowclr(i) - dqrclrfreez * hum_to_flux(i) )
+    snowcld(i) = MAX( 0. , snowcld(i) - dqrcldfreez * hum_to_flux(i) )
+
 
     !--Temperature adjustment with the uptake of latent
     !--heat because of freezing
-    temp(i) = temp(i) - dqrtotfreez * RLMLT / RCPD &
+    temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD &
                       / ( 1. + RVTMP2 * qtot(i) )
 
+    dqrtotfreez=dqrtotfreez_step1+dqrtotfreez_step2          
     !--Diagnostic tendencies
     dqrfreez(i) = dqrtotfreez / dtime
-    dqsfreez(i) = - dqrtotfreez / dtime
+    dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime
 
   ENDIF
