Index: LMDZ6/trunk/libf/phylmd/ice_sursat_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ice_sursat_mod.F90	(revision 4073)
+++ LMDZ6/trunk/libf/phylmd/ice_sursat_mod.F90	(revision 4074)
@@ -693,33 +693,35 @@
      !--critical T_LM below which no liquid contrail can form in exhaust
      !Tcontr(i,k) = 226.69+9.43*log(Gcontr-0.053)+0.72*(log(Gcontr-0.053))**2 !--K
-     Tcontr = 226.69+9.43*log(Gcontr-0.053)+0.72*(log(Gcontr-0.053))**2 !--K
-     !print *,'Tcontr=',iter,i,k,eps_w,pplay,Gcontr,Tcontr(i,k)
-     !--Psat with index 0 in FOEEW to get saturation wrt liquid water corresponding to Tcontr
-     !qsatliqcontr = RESTT*FOEEW(Tcontr(i,k),0.)                               !--Pa
-     qsatliqcontr = RESTT*FOEEW(Tcontr,0.)                               !--Pa
-     !--Critical water vapour above which there is contrail formation for ambiant temperature 
-     !qcontr(i,k) = Gcontr*(t-Tcontr(i,k)) + qsatliqcontr                      !--Pa
-     qcontr = Gcontr*(t-Tcontr) + qsatliqcontr                      !--Pa
-     !--Conversion of qcontr in specific humidity - method 1
-     !qcontr(i,k) = RD/RV*qcontr(i,k)/pplay      !--so as to return to something similar to R2ES*FOEEW/pplay
-     qcontr2 = RD/RV*qcontr/pplay      !--so as to return to something similar to R2ES*FOEEW/pplay
-     !qcontr(i,k) = min(0.5,qcontr(i,k))         !--and then we apply the same correction term as for qsat
-     qcontr2 = min(0.5,qcontr2)         !--and then we apply the same correction term as for qsat
-     !zcor = 1./(1.-RETV*qcontr(i,k))            !--for consistency with qsat but is it correct at all?
-     zcor = 1./(1.-RETV*qcontr2)            !--for consistency with qsat but is it correct at all as p is dry?
-     !zcor = 1./(1.+qcontr2)                 !--for consistency with qsat
-     !qcontr(i,k) = qcontr(i,k)*zcor
-     qcontr2 = qcontr2*zcor
-     qcontr2=MAX(1.e-10,qcontr2)            !--eliminate negative values due to extrapolation on dilution curve
-     !--Conversion of qcontr in specific humidity - method 2
-     !qcontr(i,k) = eps_w*qcontr(i,k) / (pplay+eps_w*qcontr(i,k))
-     !qcontr=MAX(1.E-10,qcontr)
-     !qcontr2 = eps_w*qcontr / (pplay+eps_w*qcontr)
-     !
-     IF (t .LT. Tcontr) THEN !--contrail formation is possible
-     !
-     !--compute fractions of persistent (P) and non-persistent(N) contrail potential regions
-     !!IF (qcontr(i,k).GE.qsat) THEN
-     IF (qcontr2.GE.qsat) THEN
+     IF (Gcontr .GT. 0.1) THEN
+     !
+       Tcontr = 226.69+9.43*log(Gcontr-0.053)+0.72*(log(Gcontr-0.053))**2 !--K
+       !print *,'Tcontr=',iter,i,k,eps_w,pplay,Gcontr,Tcontr(i,k)
+       !--Psat with index 0 in FOEEW to get saturation wrt liquid water corresponding to Tcontr
+       !qsatliqcontr = RESTT*FOEEW(Tcontr(i,k),0.)                               !--Pa
+       qsatliqcontr = RESTT*FOEEW(Tcontr,0.)                               !--Pa
+       !--Critical water vapour above which there is contrail formation for ambiant temperature 
+       !qcontr(i,k) = Gcontr*(t-Tcontr(i,k)) + qsatliqcontr                      !--Pa
+       qcontr = Gcontr*(t-Tcontr) + qsatliqcontr                      !--Pa
+       !--Conversion of qcontr in specific humidity - method 1
+       !qcontr(i,k) = RD/RV*qcontr(i,k)/pplay      !--so as to return to something similar to R2ES*FOEEW/pplay
+       qcontr2 = RD/RV*qcontr/pplay      !--so as to return to something similar to R2ES*FOEEW/pplay
+       !qcontr(i,k) = min(0.5,qcontr(i,k))         !--and then we apply the same correction term as for qsat
+       qcontr2 = min(0.5,qcontr2)         !--and then we apply the same correction term as for qsat
+       !zcor = 1./(1.-RETV*qcontr(i,k))            !--for consistency with qsat but is it correct at all?
+       zcor = 1./(1.-RETV*qcontr2)            !--for consistency with qsat but is it correct at all as p is dry?
+       !zcor = 1./(1.+qcontr2)                 !--for consistency with qsat
+       !qcontr(i,k) = qcontr(i,k)*zcor
+       qcontr2 = qcontr2*zcor
+       qcontr2=MAX(1.e-10,qcontr2)            !--eliminate negative values due to extrapolation on dilution curve
+       !--Conversion of qcontr in specific humidity - method 2
+       !qcontr(i,k) = eps_w*qcontr(i,k) / (pplay+eps_w*qcontr(i,k))
+       !qcontr=MAX(1.E-10,qcontr)
+       !qcontr2 = eps_w*qcontr / (pplay+eps_w*qcontr)
+       !
+       IF (t .LT. Tcontr) THEN !--contrail formation is possible
+       !
+       !--compute fractions of persistent (P) and non-persistent(N) contrail potential regions
+       !!IF (qcontr(i,k).GE.qsat) THEN
+       IF (qcontr2.GE.qsat) THEN
          !--none of the unsaturated clear sky is prone for contrail formation
          !!fcontrN(i,k) = 0.0
@@ -807,5 +809,5 @@
          fcontrP = fcontrP + MAX(0., 0.5*(1-pdf_N2)*(pdf_e1-pdf_e2))
          !
-     ELSE  !--qcontr LT qsat
+       ELSE  !--qcontr LT qsat
          !
          !--all of ISSR is prone for contrail formation
@@ -835,7 +837,9 @@
          !!fcontrN=2.0
          !
-     ENDIF
-     !
-     ENDIF !-- t < Tcontr
+       ENDIF
+       !
+       ENDIF !-- t < Tcontr
+     !
+     ENDIF !-- Gcontr > 0.1
 
   RETURN