Index: LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_ini.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_ini.F90	(revision 4974)
+++ LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_ini.F90	(revision 5019)
@@ -40,5 +40,5 @@
   !$OMP THREADPRIVATE(ztfondue)
 
-  REAL, SAVE, PROTECTED :: temp_nowater=233.15  ! temperature below which liquid water no longer exists
+  REAL, SAVE, PROTECTED :: temp_nowater=235.15  ! temperature below which liquid water no longer exists
   !$OMP THREADPRIVATE(temp_nowater)
 
@@ -363,4 +363,5 @@
     CALL getin_p('tresh_cl',tresh_cl)
     CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
+    CALL getin_p('temp_nowater',temp_nowater)
     ! for poprecip
     CALL getin_p('ok_poprecip',ok_poprecip)
@@ -441,4 +442,5 @@
     WRITE(lunout,*) 'lscp_ini, iflag_oldbug_fisrtilp', iflag_oldbug_fisrtilp
     WRITE(lunout,*) 'lscp_ini, fl_cor_ebil', fl_cor_ebil
+    WRITE(lunout,*) 'lscp_ini, temp_nowater', temp_nowater
     ! for poprecip
     WRITE(lunout,*) 'lscp_ini, ok_poprecip', ok_poprecip
Index: LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_tools.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_tools.F90	(revision 4974)
+++ LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_tools.F90	(revision 5019)
@@ -334,42 +334,77 @@
         CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
 
-    DO i=1,klon
-
-        IF (temp(i) .GE. RTT) THEN
+    DO i = 1, klon
+
+        IF ( temp(i) .GE. RTT ) THEN
             ! warm clouds: condensation at saturation wrt liquid
-            gammasat(i)=1.
-            dgammasatdt(i)=0.
-
-        ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. temp_nowater)) THEN
+            gammasat(i) = 1.
+            dgammasatdt(i) = 0.
+
+        ELSE
+            ! cold clouds: qsi > qsl
             
-            IF (iflag_gammasat .GE. 2) THEN          
-               gammasat(i)=qsl(i)/qsi(i)
-               dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
+            ! homogeneous freezing of aerosols, according to
+            ! Koop, 2000 and Ren and MacKenzie, 2005 (QJRMS)
+            ! 'Cirrus regime'
+            ! if f_homofreez > qsl / qsi, liquid nucleation
+            ! if f_homofreez < qsl / qsi, homogeneous freezing of aerosols
+            ! Note: f_homofreez = qsl / qsi for temp ~= -38degC
+            f_homofreez = a_homofreez - temp(i) / b_homofreez
+            
+            IF ( iflag_gammasat .GE. 3 ) THEN
+              ! condensation at homogeneous freezing threshold for temp < -38 degC
+              ! condensation at liquid saturation for temp > -38 degC
+              IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
+                gammasat(i) = f_homofreez
+                dgammasatdt(i) = - 1. / b_homofreez
+              ELSE
+                gammasat(i) = qsl(i) / qsi(i)
+                dgammasatdt(i) = ( dqsl(i) * qsi(i) - dqsi(i) * qsl(i) ) / qsi(i) / qsi(i)
+              ENDIF
+
+            ELSEIF ( iflag_gammasat .EQ. 2 ) THEN
+              ! condensation at homogeneous freezing threshold for temp < -38 degC
+              ! condensation at a threshold linearly decreasing between homogeneous
+              ! freezing and ice saturation for -38 degC < temp < temp_nowater
+              ! condensation at ice saturation for temp > temp_nowater
+              ! If temp_nowater = 235.15 K, this is equivalent to iflag_gammasat = 1
+              IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
+                gammasat(i) = f_homofreez
+                dgammasatdt(i) = - 1. / b_homofreez
+              ELSEIF ( temp(i) .LE. temp_nowater ) THEN
+                ! Here, we assume that f_homofreez = qsl / qsi for temp = -38 degC = 235.15 K
+                dgammasatdt(i) = ( a_homofreez - 235.15 / b_homofreez - 1. ) &
+                               / ( 235.15 - temp_nowater )
+                gammasat(i) = dgammasatdt(i) * ( temp(i) - temp_nowater ) + 1.
+              ELSE
+                gammasat(i) = 1.
+                dgammasatdt(i) = 0.
+              ENDIF
+
+            ELSEIF ( iflag_gammasat .EQ. 1 ) THEN
+              ! condensation at homogeneous freezing threshold for temp < -38 degC
+              ! condensation at ice saturation for temp > -38 degC
+              IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
+                gammasat(i) = f_homofreez
+                dgammasatdt(i) = - 1. / b_homofreez
+              ELSE
+                gammasat(i) = 1.
+                dgammasatdt(i) = 0.
+              ENDIF
+
             ELSE
-               gammasat(i)=1.
-               dgammasatdt(i)=0.
+              ! condensation at ice saturation for temp < -38 degC
+              ! condensation at ice saturation for temp > -38 degC
+              gammasat(i) = 1.
+              dgammasatdt(i) = 0.
+
             ENDIF
 
-        ELSE
-
-            IF (iflag_gammasat .GE.1) THEN
-               ! homogeneous freezing of aerosols, according to
-               ! Koop, 2000 and Ren and MacKenzie, 2005 (QJRMS)
-               ! 'Cirrus regime'
-               f_homofreez=a_homofreez-temp(i)/b_homofreez
-               IF (f_homofreez .GT. qsl(i)/qsi(i)) THEN
-                  gammasat(i)=qsl(i)/qsi(i)
-                  dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
-               ELSE
-                  gammasat(i)= (1.-delta_hetfreez) + delta_hetfreez * f_homofreez
-                  dgammasatdt(i)=-delta_hetfreez/b_homofreez
-               ENDIF
-           
-            ELSE
-
-               gammasat(i)=1.
-               dgammasatdt(i)=0.
-
-            ENDIF
+            ! Note that the delta_hetfreez parameter allows to linearly decrease the
+            ! condensation threshold between the calculated threshold and the ice saturation
+            ! for delta_hetfreez = 1, the threshold is the calculated condensation threshold
+            ! for delta_hetfreez = 0, the threshold is the ice saturation
+            gammasat(i) = ( 1. - delta_hetfreez ) + delta_hetfreez * gammasat(i)
+            dgammasatdt(i) = delta_hetfreez * dgammasatdt(i)
 
         ENDIF
