Index: LMDZ6/trunk/libf/phylmd/lscp_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 3999)
+++ LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 4059)
@@ -6,12 +6,13 @@
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE LSCP(dtime,paprs,pplay,t,q,ptconv,ratqs,     &
-     d_t, d_q, d_ql, d_qi, rneb, radliq, radicefrac,    &
-     rain, snow,                                        &
+SUBROUTINE LSCP(dtime,missing_val,                      & 
+     paprs,pplay,t,q,ptconv,ratqs,                      &
+     d_t, d_q, d_ql, d_qi, rneb, rneb_seri,             & 
+     radliq, radicefrac, rain, snow,                    &
      pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
      frac_impa, frac_nucl, beta,                        &
      prfl, psfl, rhcl, zqta, fraca,                     &
      ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
-     iflag_ice_thermo)
+     iflag_ice_thermo, iflag_ice_sursat)
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -93,9 +94,11 @@
 USE phys_local_var_mod, ONLY: rneblsvol
 USE lscp_tools_mod, ONLY : CALC_QSAT_ECMWF, ICEFRAC_LSCP, CALC_GAMMASAT, FALLICE_VELOCITY
-
+USE ice_sursat_mod
+!--ice supersaturation
+USE phys_local_var_mod, ONLY: zqsats, zqsatl
+USE phys_local_var_mod, ONLY: qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss
+USE phys_local_var_mod, ONLY: Tcontr, qcontr, qcontr2, fcontrN, fcontrP
 
 IMPLICIT NONE
-
-
 
 !===============================================================================
@@ -114,4 +117,6 @@
 
   REAL,                            INTENT(IN)   :: dtime           ! time step [s]
+  REAL, INTENT(IN)                              :: missing_val     ! missing value for output
+
   REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
@@ -121,4 +126,6 @@
   INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
                                                                    ! CR: if iflag_ice_thermo=2, only convection is active    
+  INTEGER,                         INTENT(IN)   :: iflag_ice_sursat ! 0 = sursat desativee, 1 = sursat activee
+
   LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
 
@@ -138,4 +145,7 @@
   REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratqs            ! function of pressure that sets the large-scale
                                                                     ! cloud PDF (sigma=ratqs*qt)
+
+  ! Input sursaturation en glace
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rneb_seri        ! fraction nuageuse en memoire
   
   ! OUTPUT variables
@@ -388,4 +398,13 @@
 d_tot_zneb(:) = 0.0   
 
+!--ice sursaturation
+gamma_ss(:,:) = 1.
+qss(:,:) = 0.
+rnebss(:,:) = 0.
+Tcontr(:,:) = missing_val
+qcontr(:,:) = missing_val
+qcontr2(:,:) = missing_val
+fcontrN(:,:) = 0.0
+fcontrP(:,:) = 0.0
 
 !===============================================================================
@@ -645,6 +664,5 @@
             qtot=zq(i)+zmqc(i)
             CALL CALC_QSAT_ECMWF(zt(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
-            zdqsdT_raw(i) = zdqs(i)*  &
-            & RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
+            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
          
             IF (zq(i) .LT. 1.e-15) THEN
@@ -677,6 +695,4 @@
         qincloud_mpc(:)=0.
 
-
-
         IF (iflag_cld_th.GE.5) THEN
 
@@ -778,5 +794,4 @@
                         ! new temperature:
                         Tbef(i)=Tbef(i)+DT(i)
-
 
                         ! Rneb, qzn and zcond for lognormal PDFs
@@ -800,11 +815,31 @@
                         zpdf_e2(i)=1.-erf(zpdf_e2(i))
              
-                        IF (zpdf_e1(i).LT.1.e-10) THEN
-                            rneb(i,k)=0.
-                            zqn(i)=gammasat(i)*zqs(i)
+                        !--ice sursaturation by Audran
+                        IF ((iflag_ice_sursat.EQ.0).OR.(Tbef(i).GT.t_glace_min)) THEN
+
+                          IF (zpdf_e1(i).LT.1.e-10) THEN
+                              rneb(i,k)=0.
+                              zqn(i)=gammasat(i)*zqs(i)
+                          ELSE
+                              rneb(i,k)=0.5*zpdf_e1(i)
+                              zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
+                          ENDIF
+
+                          rnebss(i,k)=0.0   !--ajout OB (necessaire car boucle de convergence sur le temps)
+                          fcontrN(i,k)=0.0  !--idem
+                          fcontrP(i,k)=0.0  !--idem
+                          qss(i,k)=0.0      !--idem
+
                         ELSE
-                            rneb(i,k)=0.5*zpdf_e1(i)
-                            zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
-                        ENDIF
+                        !------------------------------------
+                        ! SURSATURATION EN GLACE
+                        !------------------------------------
+
+                        CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), &
+                             gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),               &
+                             rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                 &
+                             Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )
+
+                        ENDIF ! ((flag_ice_sursat.eq.0).or.(Tbef(i).gt.t_glace_min))
 
                         ! If vertical heterogeneity, change fraction by volume as well
@@ -823,5 +858,4 @@
                        ! EV: calculation of icefrac in one sole function
                         CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
-
                 
                         IF (zfice(i).LT.1) THEN
@@ -985,5 +1019,4 @@
     ! remaining water in the cloud during the time step that is seen by the radiation
     ! -------------------------------------------------------------------------------
-
      
     DO n = 1, ninter
@@ -1128,7 +1161,4 @@
    ENDDO
 
-
-
-                                                                                                         
     ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min                                                                              
     ! if iflag_evap_pre=4
@@ -1137,5 +1167,4 @@
         DO i=1, klon                                        
                 
-
             IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN     
                 znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
@@ -1144,5 +1173,4 @@
                 znebprecipclr(i)=0.                                                                   
             ENDIF  
-
                                                                                                          
             IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN                                                  
@@ -1152,5 +1180,4 @@
                 znebprecipcld(i)=0.                                                                    
             ENDIF     
-
 
         ENDDO       
@@ -1190,10 +1217,7 @@
         ENDIF
 
-        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
-        * (paprs(i,k)-paprs(i,k+1))/RG
-
+        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
 
         IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
-
 
             IF (t(i,k) .GE. t_glace_min) THEN
@@ -1202,5 +1226,4 @@
                 zalpha_tr = a_tr_sca(4)
             ENDIF
-
 
             zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
@@ -1241,11 +1264,28 @@
     ENDDO
      
-
-END DO
+    !--save some variables for ice sursaturation
+    !
+    DO i = 1, klon
+        ! pour la mémoire
+        rneb_seri(i,k) = rneb(i,k)
+
+        ! pour les diagnostics
+        rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)
+
+        qvc(i,k) = zqs(i) * rneb(i,k)
+        qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k))  !--ajout OB a cause de cas pathologiques avec lognormale=F
+        qcld(i,k) = qvc(i,k) + zcond(i)
+
+        !q_sat
+        CALL CALC_QSAT_ECMWF(Tbef(i),0.,pplay(i,k),RTT,1,.false.,zqsatl(i,k),zdqs(i))
+        CALL CALC_QSAT_ECMWF(Tbef(i),0.,pplay(i,k),RTT,2,.false.,zqsats(i,k),zdqs(i))
+
+     ENDDO
+
+ENDDO
 
 !======================================================================
 !                      END OF VERTICAL LOOP
 !======================================================================
-  
 
   ! Rain or snow at the surface (depending on the first layer temperature)
@@ -1254,6 +1294,4 @@
       rain(i) = zrfl(i)
   ENDDO
-   
-  
 
   IF (ncoreczq>0) THEN
@@ -1261,7 +1299,4 @@
   ENDIF
 
-
-
-
 END SUBROUTINE LSCP
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
