Index: LMDZ6/branches/cirrus/libf/phylmdiso/ice_sursat_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmdiso/ice_sursat_mod.F90	(revision 5054)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/ice_sursat_mod.F90
Index: LMDZ6/branches/cirrus/libf/phylmdiso/lmdz_lscp_condensation.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmdiso/lmdz_lscp_condensation.F90	(revision 5055)
+++ LMDZ6/branches/cirrus/libf/phylmdiso/lmdz_lscp_condensation.F90	(revision 5055)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_lscp_condensation.F90
Index: LMDZ6/branches/cirrus/libf/phylmdiso/phyetat0_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmdiso/phyetat0_mod.F90	(revision 5054)
+++ LMDZ6/branches/cirrus/libf/phylmdiso/phyetat0_mod.F90	(revision 5055)
@@ -25,5 +25,6 @@
        du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
        falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien, &
-       ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, qbs_ancien, rneb_ancien, radpas, radsol, rain_fall, ratqs, &
+       ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, qbs_ancien, &
+       cf_ancien, rvc_ancien, radpas, radsol, rain_fall, ratqs, &
        rnebcon, rugoro, sig1, snow_fall, bs_fall, solaire_etat0, sollw, sollwdown, &
        solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
@@ -412,5 +413,4 @@
   ancien_ok=ancien_ok.AND.phyetat0_get(ql_ancien,"QLANCIEN","QLANCIEN",0.)
   ancien_ok=ancien_ok.AND.phyetat0_get(qs_ancien,"QSANCIEN","QSANCIEN",0.)
-  ancien_ok=ancien_ok.AND.phyetat0_get(rneb_ancien,"RNEBANCIEN","RNEBANCIEN",0.)
   ancien_ok=ancien_ok.AND.phyetat0_get(u_ancien,"UANCIEN","UANCIEN",0.)
   ancien_ok=ancien_ok.AND.phyetat0_get(v_ancien,"VANCIEN","VANCIEN",0.)
@@ -427,5 +427,13 @@
      prbsw_ancien(:)=0.
   ENDIF
-
+  
+  ! cas specifique des variables de la sursaturation par rapport a la glace
+  IF ( ok_ice_supersat ) THEN
+    ancien_ok=ancien_ok.AND.phyetat0_get(cf_ancien,"CFANCIEN","CFANCIEN",0.)
+    ancien_ok=ancien_ok.AND.phyetat0_get(rvc_ancien,"RVCANCIEN","RVCANCIEN",0.)
+  ELSE
+    cf_ancien(:,:)=0.
+    rvc_ancien(:,:)=0.
+  ENDIF
 
   ! Ehouarn: addtional tests to check if t_ancien, q_ancien contain
@@ -435,5 +443,4 @@
        (maxval(ql_ancien).EQ.minval(ql_ancien))     .OR. &
        (maxval(qs_ancien).EQ.minval(qs_ancien))     .OR. &
-       (maxval(rneb_ancien).EQ.minval(rneb_ancien)) .OR. &
        (maxval(prw_ancien).EQ.minval(prw_ancien))   .OR. &
        (maxval(prlw_ancien).EQ.minval(prlw_ancien)) .OR. &
@@ -450,4 +457,10 @@
   ENDIF
 
+  IF ( ok_ice_supersat ) THEN
+    IF ( (maxval(cf_ancien).EQ.minval(cf_ancien))     .OR. &
+         (maxval(rvc_ancien).EQ.minval(rvc_ancien)) ) THEN
+       ancien_ok=.false.
+     ENDIF
+  ENDIF
 
   found=phyetat0_get(clwcon,"CLWCON","CLWCON",0.)
Index: LMDZ6/branches/cirrus/libf/phylmdiso/phyredem.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmdiso/phyredem.F90	(revision 5054)
+++ LMDZ6/branches/cirrus/libf/phylmdiso/phyredem.F90	(revision 5055)
@@ -19,6 +19,7 @@
                                 zval, rugoro, t_ancien, q_ancien,            &
                                 prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien,      &
-                                ql_ancien, qs_ancien, qbs_ancien,  u_ancien, &
-                                v_ancien, clwcon, rnebcon, ratqs, pbl_tke,   &
+                                ql_ancien, qs_ancien, qbs_ancien, cf_ancien, &
+                                rvc_ancien, u_ancien, v_ancien,              &
+                                clwcon, rnebcon, ratqs, pbl_tke,             &
                                 wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
                                 wake_deltat, wake_deltaq, wake_s, wake_dens, &
@@ -265,4 +266,9 @@
        CALL put_field(pass,"QBSANCIEN", "QBSANCIEN", qbs_ancien)
        CALL put_field(pass,"PRBSWANCIEN", "PRBSWANCIEN", prbsw_ancien)
+    ENDIF
+
+    IF ( ok_ice_supersat ) THEN
+      CALL put_field(pass,"CFANCIEN", "CFANCIEN", cf_ancien)
+      CALL put_field(pass,"RVCANCIEN", "RVCANCIEN", rvc_ancien)
     ENDIF
 
Index: LMDZ6/branches/cirrus/libf/phylmdiso/phys_local_var_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmdiso/phys_local_var_mod.F90	(revision 5054)
+++ LMDZ6/branches/cirrus/libf/phylmdiso/phys_local_var_mod.F90	(revision 5055)
@@ -18,8 +18,6 @@
       REAL, SAVE, ALLOCATABLE :: u_seri(:,:), v_seri(:,:)
       !$OMP THREADPRIVATE(u_seri, v_seri)
-      REAL, SAVE, ALLOCATABLE :: rneb_seri(:,:)
-      !$OMP THREADPRIVATE(rneb_seri)
-      REAL, SAVE, ALLOCATABLE :: d_rneb_dyn(:,:)
-      !$OMP THREADPRIVATE(d_rneb_dyn)
+      REAL, SAVE, ALLOCATABLE :: cf_seri(:,:), rvc_seri(:,:)
+      !$OMP THREADPRIVATE(cf_seri, rvc_seri)
       REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:),l_mix(:,:,:),wprime(:,:,:)
       !$OMP THREADPRIVATE(l_mixmin, l_mix, wprime)
@@ -38,4 +36,6 @@
       REAL, SAVE, ALLOCATABLE :: d_u_dyn(:,:), d_v_dyn(:,:)
       !$OMP THREADPRIVATE(d_u_dyn, d_v_dyn)
+      REAL, SAVE, ALLOCATABLE :: d_cf_dyn(:,:), d_rvc_dyn(:,:)
+      !$OMP THREADPRIVATE(d_cf_dyn, d_rvc_dyn)
       REAL, SAVE, ALLOCATABLE :: d_tr_dyn(:,:,:)
       !$OMP THREADPRIVATE(d_tr_dyn)
@@ -348,5 +348,5 @@
 !$OMP THREADPRIVATE(cdragm, cdragh)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: cldh, cldl, cldm, cldq, cldt, qsat2m
-!$OMP THREADPRIVATE(cldh, cldl, cldm, cldq, cldt, qsat2m )
+!$OMP THREADPRIVATE(cldh, cldl, cldm, cldq, cldt, qsat2m)
 !AS: cldhjn, cldljn, cldmjn,cldtjn pas utilisés en tant que variables, juste noms de diagnostics
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: JrNt
@@ -573,4 +573,5 @@
 !$OMP THREADPRIVATE(phiwriteSTD, qwriteSTD, twriteSTD, rhwriteSTD)
 
+
 ! ug et d'autres encore:
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: beta_prec
@@ -580,8 +581,4 @@
       REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: pfraclr,pfracld
 !$OMP THREADPRIVATE(pfraclr,pfracld)
-      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: distcltop
-!$OMP THREADPRIVATE(distcltop)
-      REAL, SAVE, ALLOCATABLE :: temp_cltop(:,:) 
-!$OMP THREADPRIVATE(temp_cltop)
 
 ! variables de sorties MM
@@ -605,46 +602,39 @@
 !$OMP THREADPRIVATE(zn2mout)
 
-      REAL, SAVE, ALLOCATABLE :: qclr(:,:)
-      !$OMP THREADPRIVATE(qclr)
-      REAL, SAVE, ALLOCATABLE :: qcld(:,:)
-      !$OMP THREADPRIVATE(qcld)
-      REAL, SAVE, ALLOCATABLE :: qss(:,:)
-      !$OMP THREADPRIVATE(qss)
-      REAL, SAVE, ALLOCATABLE :: qvc(:,:)
-      !$OMP THREADPRIVATE(qvc)
-      REAL, SAVE, ALLOCATABLE :: rnebclr(:,:)
-      !$OMP THREADPRIVATE(rnebclr)
-      REAL, SAVE, ALLOCATABLE :: rnebss(:,:)
-      !$OMP THREADPRIVATE(rnebss)
-      REAL, SAVE, ALLOCATABLE :: gamma_ss(:,:)
-      !$OMP THREADPRIVATE(gamma_ss)
-      REAL, SAVE, ALLOCATABLE :: N1_ss(:,:)
-      !$OMP THREADPRIVATE(N1_ss)
-      REAL, SAVE, ALLOCATABLE :: N2_ss(:,:)
-      !$OMP THREADPRIVATE(N2_ss)
-      REAL, SAVE, ALLOCATABLE :: drneb_sub(:,:)
-      !$OMP THREADPRIVATE(drneb_sub)
-      REAL, SAVE, ALLOCATABLE :: drneb_con(:,:)
-      !$OMP THREADPRIVATE(drneb_con)
-      REAL, SAVE, ALLOCATABLE :: drneb_tur(:,:)
-      !$OMP THREADPRIVATE(drneb_tur)
-      REAL, SAVE, ALLOCATABLE :: drneb_avi(:,:)
-      !$OMP THREADPRIVATE(drneb_avi)
-      REAL, SAVE, ALLOCATABLE :: zqsatl(:,:)
-      !$OMP THREADPRIVATE(zqsatl)
-      REAL, SAVE, ALLOCATABLE :: zqsats(:,:)
-      !$OMP THREADPRIVATE(zqsats)
-      REAL, SAVE, ALLOCATABLE :: Tcontr(:,:)
-      !$OMP THREADPRIVATE(Tcontr)
-      REAL, SAVE, ALLOCATABLE :: qcontr(:,:)
-      !$OMP THREADPRIVATE(qcontr)
-      REAL, SAVE, ALLOCATABLE :: qcontr2(:,:)
-      !$OMP THREADPRIVATE(qcontr2)
-      REAL, SAVE, ALLOCATABLE :: fcontrN(:,:)
-      !$OMP THREADPRIVATE(fcontrN)
-      REAL, SAVE, ALLOCATABLE :: fcontrP(:,:)
-      !$OMP THREADPRIVATE(fcontrP)
-
-!--POPRECIP variables
+!-- LSCP - condensation and ice supersaturation variables
+      REAL, SAVE, ALLOCATABLE :: qsub(:,:), qissr(:,:), qcld(:,:)
+      !$OMP THREADPRIVATE(qsub, qissr, qcld)
+      REAL, SAVE, ALLOCATABLE :: subfra(:,:), issrfra(:,:)
+      !$OMP THREADPRIVATE(subfra, issrfra)
+      REAL, SAVE, ALLOCATABLE :: gamma_cond(:,:)
+      !$OMP THREADPRIVATE(gamma_cond)
+      REAL, SAVE, ALLOCATABLE :: ratio_qi_qtot(:,:)
+      !$OMP THREADPRIVATE(ratio_qi_qtot)
+      REAL, SAVE, ALLOCATABLE :: dcf_sub(:,:), dcf_con(:,:), dcf_mix(:,:)
+      !$OMP THREADPRIVATE(dcf_sub, dcf_con, dcf_mix)
+      REAL, SAVE, ALLOCATABLE :: dqi_adj(:,:), dqi_sub(:,:), dqi_con(:,:), dqi_mix(:,:)
+      !$OMP THREADPRIVATE(dqi_adj, dqi_sub, dqi_con, dqi_mix)
+      REAL, SAVE, ALLOCATABLE :: dqvc_adj(:,:), dqvc_sub(:,:), dqvc_con(:,:), dqvc_mix(:,:)
+      !$OMP THREADPRIVATE(dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix)
+      REAL, SAVE, ALLOCATABLE :: qsatliq(:,:), qsatice(:,:)
+      !$OMP THREADPRIVATE(qsatliq, qsatice)
+
+!-- LSCP - aviation and contrails variables
+      REAL, SAVE, ALLOCATABLE :: Tcontr(:,:), qcontr(:,:), qcontr2(:,:)
+      !$OMP THREADPRIVATE(Tcontr, qcontr, qcontr2)
+      REAL, SAVE, ALLOCATABLE :: fcontrN(:,:), fcontrP(:,:)
+      !$OMP THREADPRIVATE(fcontrN, fcontrP)
+      REAL, SAVE, ALLOCATABLE :: dcf_avi(:,:), dqi_avi(:,:), dqvc_avi(:,:)
+      !$OMP THREADPRIVATE(dcf_avi, dqi_avi, dqvc_avi)
+      REAL, SAVE, ALLOCATABLE :: flight_dist(:,:), flight_h2o(:,:)
+      !$OMP THREADPRIVATE(flight_dist, flight_h2o)
+
+!-- LSCP - mixed phase clouds variables
+      REAL, SAVE, ALLOCATABLE :: distcltop(:,:)
+      !$OMP THREADPRIVATE(distcltop)
+      REAL, SAVE, ALLOCATABLE :: temp_cltop(:,:)
+      !$OMP THREADPRIVATE(temp_cltop)
+
+!-- LSCP - POPRECIP variables
       REAL, SAVE, ALLOCATABLE :: qraindiag(:,:)
       !$OMP THREADPRIVATE(qraindiag)
@@ -675,5 +665,5 @@
 
 
-
+      
 
 
@@ -771,7 +761,10 @@
       ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev), qbs_seri(klon,klev))
       ALLOCATE(u_seri(klon,klev),v_seri(klon,klev))
+      ALLOCATE(cf_seri(klon,klev),rvc_seri(klon,klev))
       ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
+      ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))
+      pbl_eps(:,:,:)=0.
       l_mix(:,:,:)=0.;l_mixmin(:,:,:)=0.;wprime(:,:,:)=0. ! doit etre initialse car pas toujours remplis
-      ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))      
+      ALLOCATE(rhcl(klon,klev))
       ALLOCATE(tr_seri(klon,klev,nbtr))
       ALLOCATE(d_t_dyn(klon,klev),d_q_dyn(klon,klev))
@@ -779,4 +772,5 @@
       ALLOCATE(d_q_dyn2d(klon),d_ql_dyn2d(klon),d_qs_dyn2d(klon), d_qbs_dyn2d(klon))
       ALLOCATE(d_u_dyn(klon,klev),d_v_dyn(klon,klev))
+      ALLOCATE(d_cf_dyn(klon,klev),d_rvc_dyn(klon,klev))
       ALLOCATE(d_tr_dyn(klon,klev,nbtr))                   !RomP
       ALLOCATE(d_t_con(klon,klev),d_q_con(klon,klev),d_q_con_zmasse(klon,klev))
@@ -1103,4 +1097,5 @@
       ALLOCATE(temp_cltop(klon,klev))
 
+
       ALLOCATE (zxsnow(klon),snowhgt(klon),qsnow(klon),to_ice(klon))
       ALLOCATE (sissnow(klon),runoff(klon),albsol3_lic(klon))
@@ -1117,17 +1112,22 @@
       ALLOCATE(zn2mout(klon,6))
 
-! Supersaturation
-      ALLOCATE(rneb_seri(klon,klev))
-      ALLOCATE(d_rneb_dyn(klon,klev))
-      ALLOCATE(qclr(klon,klev), qcld(klon,klev), qss(klon,klev), qvc(klon,klev))
-      ALLOCATE(rnebclr(klon,klev), rnebss(klon,klev), gamma_ss(klon,klev))
-      ALLOCATE(N1_ss(klon,klev), N2_ss(klon,klev))
-      ALLOCATE(drneb_sub(klon,klev), drneb_con(klon,klev), drneb_tur(klon,klev), drneb_avi(klon,klev))
-      ALLOCATE(zqsatl(klon,klev), zqsats(klon,klev))
-      ALLOCATE(Tcontr(klon,klev), qcontr(klon,klev), qcontr2(klon,klev), fcontrN(klon,klev), fcontrP(klon,klev))
-
-!--POPRECIP variables
-      ALLOCATE(dqreva(klon,klev),dqssub(klon,klev))
-      ALLOCATE(qraindiag(klon,klev), qsnowdiag(klon,klev)) 
+!-- LSCP - condensation and ice supersaturation variables
+      ALLOCATE(qsub(klon,klev), qissr(klon,klev), qcld(klon,klev))
+      ALLOCATE(subfra(klon,klev), issrfra(klon,klev))
+      ALLOCATE(gamma_cond(klon,klev), ratio_qi_qtot(klon,klev))
+      ALLOCATE(dcf_sub(klon,klev), dcf_con(klon,klev), dcf_mix(klon,klev))
+      ALLOCATE(dqi_adj(klon,klev), dqi_sub(klon,klev), dqi_con(klon,klev), dqi_mix(klon,klev))
+      ALLOCATE(dqvc_adj(klon,klev), dqvc_sub(klon,klev), dqvc_con(klon,klev), dqvc_mix(klon,klev))
+      ALLOCATE(qsatliq(klon,klev), qsatice(klon,klev))
+
+!-- LSCP - aviation and contrails variables
+      ALLOCATE(Tcontr(klon,klev), qcontr(klon,klev), qcontr2(klon,klev))
+      ALLOCATE(fcontrN(klon,klev), fcontrP(klon,klev))
+      ALLOCATE(dcf_avi(klon,klev), dqi_avi(klon,klev), dqvc_avi(klon,klev))
+      ALLOCATE(flight_dist(klon,klev), flight_h2o(klon,klev))
+
+!-- LSCP - POPRECIP variables
+      ALLOCATE(qraindiag(klon,klev), qsnowdiag(klon,klev))
+      ALLOCATE(dqreva(klon,klev), dqssub(klon,klev))
       ALLOCATE(dqrauto(klon,klev), dqrcol(klon,klev), dqrmelt(klon,klev), dqrfreez(klon,klev))
       ALLOCATE(dqsauto(klon,klev), dqsagg(klon,klev), dqsrim(klon,klev), dqsmelt(klon,klev), dqsfreez(klon,klev))
@@ -1182,7 +1182,8 @@
       DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri, qbs_seri)
       DEALLOCATE(u_seri,v_seri)
+      DEALLOCATE(cf_seri,rvc_seri)
       DEALLOCATE(l_mixmin,l_mix,wprime)
       DEALLOCATE(pbl_eps)
-
+      DEALLOCATE(rhcl)
       DEALLOCATE(tr_seri)
       DEALLOCATE(d_t_dyn,d_q_dyn)
@@ -1190,4 +1191,5 @@
       DEALLOCATE(d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, d_qbs_dyn2d)
       DEALLOCATE(d_u_dyn,d_v_dyn)
+      DEALLOCATE(d_cf_dyn,d_rvc_dyn)
       DEALLOCATE(d_tr_dyn)                      !RomP
       DEALLOCATE(d_t_con,d_q_con,d_q_con_zmasse)
@@ -1352,5 +1354,5 @@
       DEALLOCATE(zxfqcalving, zxfluxlat)
       DEALLOCATE(zxrunofflic)
-      DEALLOCATE(zxustartlic, zxrhoslic)
+      DEALLOCATE(zxustartlic, zxrhoslic, zxqsaltlic)
       DEALLOCATE(zxtsol, snow_lsc, zxfqfonte, zxqsurf)
       DEALLOCATE(rain_lsc)
@@ -1482,19 +1484,24 @@
       DEALLOCATE(zn2mout)
 
-! Supersaturation
-      DEALLOCATE(rneb_seri)
-      DEALLOCATE(d_rneb_dyn)
-      DEALLOCATE(qclr, qcld, qss, qvc)
-      DEALLOCATE(rnebclr, rnebss, gamma_ss)
-      DEALLOCATE(N1_ss, N2_ss)
-      DEALLOCATE(drneb_sub, drneb_con, drneb_tur, drneb_avi)
-      DEALLOCATE(zqsatl, zqsats)
-      DEALLOCATE(Tcontr, qcontr, qcontr2, fcontrN, fcontrP)
-
-!--POPRECIP variables
-      DEALLOCATE(qraindiag, qsnowdiag) 
-      DEALLOCATE(dqreva,dqssub)
-      DEALLOCATE(dqrauto,dqrcol,dqrmelt,dqrfreez)
-      DEALLOCATE(dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez)
+!-- LSCP - condensation and ice supersaturation variables
+      DEALLOCATE(qsub, qissr, qcld)
+      DEALLOCATE(subfra, issrfra)
+      DEALLOCATE(gamma_cond, ratio_qi_qtot)
+      DEALLOCATE(dcf_sub, dcf_con, dcf_mix)
+      DEALLOCATE(dqi_adj, dqi_sub, dqi_con, dqi_mix)
+      DEALLOCATE(dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix)
+      DEALLOCATE(qsatliq, qsatice)
+
+!-- LSCP - aviation and contrails variables
+      DEALLOCATE(Tcontr, qcontr, qcontr2)
+      DEALLOCATE(fcontrN, fcontrP)
+      DEALLOCATE(dcf_avi, dqi_avi, dqvc_avi)
+      DEALLOCATE(flight_dist, flight_h2o)
+
+!-- LSCP - POPRECIP variables
+      DEALLOCATE(qraindiag, qsnowdiag)
+      DEALLOCATE(dqreva, dqssub)
+      DEALLOCATE(dqrauto, dqrcol, dqrmelt, dqrfreez)
+      DEALLOCATE(dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
 
 #ifdef CPP_StratAer
Index: LMDZ6/branches/cirrus/libf/phylmdiso/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmdiso/phys_output_ctrlout_mod.F90	(revision 5054)
+++ LMDZ6/branches/cirrus/libf/phylmdiso/phys_output_ctrlout_mod.F90	(revision 5055)
@@ -450,10 +450,8 @@
   TYPE(ctrl_out), SAVE :: o_SWdn200clr = ctrl_out((/ 10, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'SWdn200clr', 'SWdn clear sky at 200mb', 'W/m2', (/ ('', i=1, 10) /))
-
-  ! arajouter
-  !  type(ctrl_out),save :: o_LWupTOA     = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'LWupTOA', &
-  !    (/ ('', i=1, 10) /))
-  !  type(ctrl_out),save :: o_LWupTOAclr  = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'LWupTOAclr', &
-  !    (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_LWupTOA = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'LWupTOA', 'LWup at TOA', 'W/m2', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_LWupTOAclr = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'LWupTOAclr', 'LWup clear sky at TOA', 'W/m2', (/ ('', i=1, 10) /))
   !  type(ctrl_out),save :: o_LWdnTOA     = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'LWdnTOA', &
   !    (/ ('', i=1, 10) /))
@@ -2093,51 +2091,73 @@
     'albslw3', 'Surface albedo LW3', '-', (/ ('', i=1, 10) /))
 
-!--aviation & supersaturation
-  TYPE(ctrl_out), SAVE :: o_oclr = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'oclr', 'Clear sky total water', 'kg/kg', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_ocld = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'ocld', 'Cloudy sky total water', 'kg/kg', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_oss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'oss', 'ISSR total water', 'kg/kg', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_ovc = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'ovc', 'In-cloup vapor', 'kg/kg', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_rnebclr = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'rnebclr', 'Clear sky fraction', '-', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_rnebss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'rnebss', 'ISSR fraction', '-', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_rnebseri = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'rnebseri', 'Cloud fraction', '-', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_gammass = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'gammass', 'Gamma supersaturation', '', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_N1_ss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'N1ss', 'N1', '', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_N2_ss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'N2ss', 'N2', '', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_drnebsub = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'drnebsub', 'Cloud fraction change because of sublimation', 's-1', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_drnebcon = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'drnebcon', 'Cloud fraction change because of condensation', 's-1', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_drnebtur = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'drnebtur', 'Cloud fraction change because of turbulence', 's-1', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_drnebavi = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'drnebavi', 'Cloud fraction change because of aviation', 's-1', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_qsatl = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'qsatl', 'Saturation with respect to liquid water', '', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_qsats = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'qsats', 'Saturation with respect to solid water', '', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_flight_m = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'flightm', 'Flown meters', 'm/s/mesh', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_flight_h2o = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
-    'flighth2o', 'H2O flight emission', 'kg H2O/s/mesh', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_Tcontr = ctrl_out((/ 1, 1, 1, 1, 11, 11, 11, 11, 11, 11/),&
+!-- LSCP - condensation and ice supersaturation variables
+  TYPE(ctrl_out), SAVE :: o_cfseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'cfseri', 'Cloud fraction', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dcfdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dcfdyn', 'Dynamics cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_rvcseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'rvcseri', 'Cloudy water vapor to total water vapor ratio', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_drvcdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'drvcdyn', 'Dynamics cloudy water vapor to total water vapor ratio tendency', 's-1', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_qsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'qsub', 'Subsaturated clear sky total water', 'kg/kg', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_qissr = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'qissr', 'Supersaturated clear sky total water', 'kg/kg', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_qcld = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'qcld', 'Cloudy sky total water', 'kg/kg', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_subfra = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'subfra', 'Subsaturated clear sky fraction', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_issrfra = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'issrfra', 'Supersaturated clear sky fraction', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_gammacond = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'gammacond', 'Condensation threshold w.r.t. saturation', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dcfsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dcfsub', 'Sublimation cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dcfcon = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dcfcon', 'Condensation cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dcfmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dcfmix', 'Cloud mixing cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqiadj = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqiadj', 'Temperature adjustment ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqisub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqisub', 'Sublimation ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqicon = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqicon', 'Condensation ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqimix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqimix', 'Cloud mixing ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqvcadj = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqvcadj', 'Temperature adjustment cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqvcsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqvcsub', 'Sublimation cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqvccon = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqvccon', 'Condensation cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqvcmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqvcmix', 'Cloud mixing cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_qsatl = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'qsatl', 'Saturation with respect to liquid', 'kg/kg', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_qsati = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'qsati', 'Saturation with respect to ice', 'kg/kg', (/ ('', i=1, 10) /))
+
+!-- LSCP - aviation variables
+  TYPE(ctrl_out), SAVE :: o_Tcontr = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     'Tcontr', 'Temperature threshold for contrail formation', 'K', (/ ('',i=1,10) /))
-  TYPE(ctrl_out), SAVE :: o_qcontr = ctrl_out((/ 1, 1, 1, 1, 11, 11, 11, 11, 11, 11/),&
-    'qcontr', 'Specific humidity threshold for contrail formation','Pa', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_qcontr2 = ctrl_out((/ 1, 1, 1, 1, 11, 11, 11, 11, 11, 11/),&
-    'qcontr2', 'Specific humidity threshold for contrail formation','kg/kg', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_fcontrN = ctrl_out((/ 2, 2, 2, 2, 2, 2, 11, 11, 11, 11/),&
+  TYPE(ctrl_out), SAVE :: o_qcontr = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
+    'qcontr', 'Specific humidity threshold for contrail formation', 'Pa', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_qcontr2 = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
+    'qcontr2', 'Specific humidity threshold for contrail formation', 'kg/kg', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_fcontrN = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     'fcontrN', 'Fraction with non-persistent contrail in clear-sky', '-', (/ ('', i=1,10)/))
-  TYPE(ctrl_out), SAVE :: o_fcontrP = ctrl_out((/ 2, 2, 2, 2, 2, 2, 11, 11, 11, 11/),&
+  TYPE(ctrl_out), SAVE :: o_fcontrP = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     'fcontrP', 'Fraction with persistent contrail in ISSR', '-', (/ ('', i=1,10)/))
+  TYPE(ctrl_out), SAVE :: o_dcfavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dcfavi', 'Aviation cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqiavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqiavi', 'Aviation ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqvcavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqvcavi', 'Aviation cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_flight_dist = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'flightdist', 'Aviation flown distance', 'm/s/mesh', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_flight_h2o = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'flighth2o', 'Aviation H2O flight emission', 'kg H2O/s/mesh', (/ ('', i=1, 10) /))
 
 !!!!!!!!!!!!! Sorties niveaux standards de pression NMC 
Index: LMDZ6/branches/cirrus/libf/phylmdiso/physiq_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmdiso/physiq_mod.F90	(revision 5054)
+++ LMDZ6/branches/cirrus/libf/phylmdiso/physiq_mod.F90	(revision 5055)
@@ -73,5 +73,4 @@
     USE tracinca_mod, ONLY: config_inca
     USE tropopause_m,     ONLY: dyn_tropopause
-    USE ice_sursat_mod,  ONLY: flight_init, airplane
     USE vampir
     USE write_field_phy
@@ -192,7 +191,10 @@
        ! [Variables internes non sauvegardees de la physique]
        ! Variables locales pour effectuer les appels en serie
-       t_seri,q_seri,ql_seri,qs_seri,qbs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
+       t_seri,q_seri,ql_seri,qs_seri,qbs_seri, &
+       u_seri,v_seri,cf_seri,rvc_seri,tr_seri, &
+       rhcl, &        
        ! Dynamic tendencies (diagnostics)
-       d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_qbs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn,d_rneb_dyn, &
+       d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_qbs_dyn, &
+       d_u_dyn,d_v_dyn,d_cf_dyn,d_rvc_dyn,d_tr_dyn, &
        d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d,d_qbs_dyn2d, &
        ! Physic tendencies
@@ -271,5 +273,5 @@
        JrNt,                             &
        dthmin, evap, snowerosion,fder, plcl, plfc,   &
-       prw, prlw, prsw, prbsw,                  &
+       prw, prlw, prsw, prbsw, water_budget,         &
        s_lcl, s_pblh, s_pblt, s_therm,   &
        cdragm, cdragh,                   &
@@ -346,7 +348,12 @@
        pfraclr,pfracld, &
        distcltop,temp_cltop, &
-       zqsatl, zqsats, &
-       qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
+       !-- LSCP - condensation and ice supersaturation variables
+       qsub, qissr, qcld, subfra, issrfra, gamma_cond, ratio_qi_qtot, &
+       dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
+       dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
+       !-- LSCP - aviation and contrails variables
        Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
+       dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
+       !
        cldemi,  &
        cldfra, cldtau, fiwc,  &
@@ -359,5 +366,5 @@
        t2m, fluxlat,  &
        fsollw, evap_pot,  &
-       fsolsw, wfbils, wfevap,  &
+       fsolsw, wfbils, wfevap, &  
        prfl, psfl,bsfl, fraca, Vprecip,  &
        zw2,  &
@@ -373,5 +380,5 @@
        rneb,  &
        zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
-       zxfluxt,zxfluxq
+       zxfluxt,zxfluxq 
 
 
@@ -549,6 +556,6 @@
     !
     ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
-    INTEGER,SAVE :: ivap, iliq, isol, irneb, ibs
-!$OMP THREADPRIVATE(ivap, iliq, isol, irneb, ibs)
+    INTEGER,SAVE :: ivap, iliq, isol, ibs, icf, irvc
+!$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, irvc)
     !
     !
@@ -921,5 +928,5 @@
     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !
-    REAL rhcl(klon,klev)    ! humiditi relative ciel clair
+!    REAL rhcl(klon,klev)    ! humiditi relative ciel clair
     REAL dialiq(klon,klev)  ! eau liquide nuageuse
     REAL diafra(klon,klev)  ! fraction nuageuse
@@ -1407,7 +1414,8 @@
        iliq = strIdx(tracers(:)%name, addPhase('H2O', 'l'))
        isol = strIdx(tracers(:)%name, addPhase('H2O', 's'))
-       irneb= strIdx(tracers(:)%name, addPhase('H2O', 'r'))
        ibs  = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
-       CALL init_etat0_limit_unstruct
+       icf  = strIdx(tracers(:)%name, addPhase('H2O', 'f'))
+       irvc = strIdx(tracers(:)%name, addPhase('H2O', 'c'))
+!       CALL init_etat0_limit_unstruct
        IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
        !CR:nvelles variables convection/poches froides
@@ -1454,25 +1462,25 @@
        ENDIF
 
-       IF (ok_ice_sursat.AND.(iflag_ice_thermo.EQ.0)) THEN
-          WRITE (lunout, *) ' ok_ice_sursat=y requires iflag_ice_thermo=1 as well'
+       IF (ok_ice_supersat.AND.(iflag_ice_thermo.EQ.0)) THEN
+          WRITE (lunout, *) ' ok_ice_supersat=y requires iflag_ice_thermo=1 as well'
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
        ENDIF
 
-       IF (ok_ice_sursat.AND.(nqo.LT.4)) THEN
-          WRITE (lunout, *) ' ok_ice_sursat=y requires 4 H2O tracers ', &
-               '(H2O_g, H2O_l, H2O_s, H2O_r) but nqo=', nqo, '. Might as well stop here.'
+       IF (ok_ice_supersat.AND.(nqo.LT.5)) THEN
+          WRITE (lunout, *) ' ok_ice_supersat=y requires 5 H2O tracers ', &
+               '(H2O_g, H2O_l, H2O_s, H2O_f, H2O_c) but nqo=', nqo, '. Might as well stop here.'
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
        ENDIF
 
-       IF (ok_plane_h2o.AND..NOT.ok_ice_sursat) THEN
-          WRITE (lunout, *) ' ok_plane_h2o=y requires ok_ice_sursat=y '
+       IF (ok_plane_h2o.AND..NOT.ok_ice_supersat) THEN
+          WRITE (lunout, *) ' ok_plane_h2o=y requires ok_ice_supersat=y '
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
        ENDIF
 
-       IF (ok_plane_contrail.AND..NOT.ok_ice_sursat) THEN
-          WRITE (lunout, *) ' ok_plane_contrail=y requires ok_ice_sursat=y '
+       IF (ok_plane_contrail.AND..NOT.ok_ice_supersat) THEN
+          WRITE (lunout, *) ' ok_plane_contrail=y requires ok_ice_supersat=y '
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
@@ -1482,5 +1490,5 @@
           abort_message='blowing snow cannot be activated with water isotopes yet'
           CALL abort_physic(modname,abort_message, 1)
-         IF ((ok_ice_sursat.AND.nqo .LT.5).OR.(.NOT.ok_ice_sursat.AND.nqo.LT.4)) THEN
+         IF ((ok_ice_supersat.AND.nqo .LT.6).OR.(.NOT.ok_ice_supersat.AND.nqo.LT.4)) THEN
              WRITE (lunout, *) 'activation of blowing snow needs a specific H2O tracer', &
                                'but nqo=', nqo
@@ -1950,5 +1958,6 @@
    &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
        CALL ratqs_ini(klon,klev,iflag_thermals,lunout,nbsrf,is_lic,is_ter,RG,RV,RD,RCPD,RLSTT,RLVTT,RTT)
-       CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RPI)
+       CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_supersat,iflag_ratqs,fl_cor_ebil, &
+                     RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RV,RG,RPI,EPS_W)
        CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, &
                              RVTMP2, RTT,RD,RG, RV, RPI)
@@ -2001,5 +2010,5 @@
                               ptconv, read_climoz, clevSTD,                   &
                               ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
-                              flag_aerosol, flag_aerosol_strat, ok_cdnc,t,u1,v1)
+                              flag_aerosol, flag_aerosol_strat, ok_cdnc, t, u1, v1)
 #endif
 
@@ -2434,4 +2443,6 @@
                     sollwdown(:))
 
+      !--Init for LSCP - condensation
+      ratio_qi_qtot(:,:) = 0.
 
 
@@ -2540,24 +2551,22 @@
           q_seri(i,k)  = qx(i,k,ivap)
           ql_seri(i,k) = qx(i,k,iliq)
-          qbs_seri(i,k) = 0.
+          qbs_seri(i,k)= 0.
+          cf_seri(i,k) = 0.
+          rvc_seri(i,k)= 0.
           !CR: ATTENTION, on rajoute la variable glace
           IF (nqo.EQ.2) THEN             !--vapour and liquid only
              qs_seri(i,k) = 0.
-             rneb_seri(i,k) = 0.
           ELSE IF (nqo.EQ.3) THEN        !--vapour, liquid and ice
              qs_seri(i,k) = qx(i,k,isol)
-             rneb_seri(i,k) = 0.
-          ELSE IF (nqo.GE.4) THEN        !--vapour, liquid, ice and rneb and blowing snow
+          ELSE IF (nqo.GE.4) THEN        !--vapour, liquid, ice, blowing snow, cloud fraction and cloudy water vapor to total water vapor ratio
              qs_seri(i,k) = qx(i,k,isol)
-             IF (ok_ice_sursat) THEN
-             rneb_seri(i,k) = qx(i,k,irneb)
+             IF (ok_ice_supersat) THEN
+               cf_seri(i,k) = qx(i,k,icf)
+               rvc_seri(i,k) = qx(i,k,irvc)
              ENDIF
              IF (ok_bs) THEN
-             qbs_seri(i,k)= qx(i,k,ibs)
+               qbs_seri(i,k)= qx(i,k,ibs)
              ENDIF
-
           ENDIF
-
-
        ENDDO
     ENDDO
@@ -2717,5 +2726,7 @@
        d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
        d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
-       d_qbs_dyn(:,:) = (qbs_seri(:,:)-qbs_ancien(:,:))/phys_tstep
+       d_qbs_dyn(:,:)= (qbs_seri(:,:)-qbs_ancien(:,:))/phys_tstep
+       d_cf_dyn(:,:) = (cf_seri(:,:)-cf_ancien(:,:))/phys_tstep
+       d_rvc_dyn(:,:)= (rvc_seri(:,:)-rvc_ancien(:,:))/phys_tstep
        CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
        d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
@@ -2729,6 +2740,4 @@
        IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
        ! !! RomP <<<
-       !!d_rneb_dyn(:,:)=(rneb_seri(:,:)-rneb_ancien(:,:))/phys_tstep
-       d_rneb_dyn(:,:)=0.0
 
 #ifdef ISO
@@ -2809,4 +2818,7 @@
        d_ql_dyn(:,:) = 0.0
        d_qs_dyn(:,:) = 0.0
+       d_qbs_dyn(:,:)= 0.0
+       d_cf_dyn(:,:) = 0.0
+       d_rvc_dyn(:,:)= 0.0
        d_q_dyn2d(:)  = 0.0
        d_ql_dyn2d(:) = 0.0
@@ -2835,6 +2847,4 @@
        IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
        ! !! RomP <<<
-       d_rneb_dyn(:,:)=0.0
-       d_qbs_dyn(:,:)=0.0
        ancien_ok = .TRUE.
     ENDIF
@@ -2944,4 +2954,19 @@
           ! "zmasse" changes a little.)
        ENDIF
+    ENDIF
+
+    !-- Needed for LSCP - condensation and ice supersaturation
+    IF (ok_ice_supersat) THEN
+      DO k = 1, klev
+        DO i = 1, klon
+          IF ( ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) ) .GT. 0. ) THEN
+            ratio_qi_qtot(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
+            rvc_seri(i,k) = rvc_seri(i,k) * q_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
+          ELSE
+            ratio_qi_qtot(i,k) = 0.
+            rvc_seri(i,k) = 0.
+          ENDIF
+        ENDDO
+      ENDDO
     ENDIF
 
@@ -4872,11 +4897,11 @@
 
     !--mise à jour de flight_m et flight_h2o dans leur module
-    IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
-      CALL airplane(debut,pphis,pplay,paprs,t_seri)
-    ENDIF
+    !IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
+    !  CALL airplane(debut,pphis,pplay,paprs,t_seri)
+    !ENDIF
 
     CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay, &
          t_seri, q_seri,ptconv,ratqs, &
-         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, & 
+         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, & 
          pfraclr,pfracld, &
          radocond, picefra, rain_lsc, snow_lsc, &
@@ -4884,7 +4909,11 @@
          prfl, psfl, rhcl,  &
          zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
-         iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, distcltop, temp_cltop, &
-         qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
-         Tcontr, qcontr, qcontr2, fcontrN, fcontrP , &
+         iflag_ice_thermo, distcltop, temp_cltop, cell_area, &
+         cf_seri, rvc_seri, u_seri, v_seri, pbl_eps(:,:,is_ave), &
+         qsub, qissr, qcld, subfra, issrfra, gamma_cond, ratio_qi_qtot, &
+         dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
+         dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
+         Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
+         dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
          cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
          qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
@@ -6823,7 +6852,8 @@
              d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
           ENDIF
-          !--ice_sursat: nqo=4, on ajoute rneb
-          IF (nqo.ge.4 .and. ok_ice_sursat) THEN
-             d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
+          !--ice_supersat: nqo=5, we add cloud fraction and cloudy water vapor to total water vapor ratio
+          IF (nqo.ge.5 .and. ok_ice_supersat) THEN
+             d_qx(i,k,icf) = ( cf_seri(i,k) - qx(i,k,icf) ) / phys_tstep
+             d_qx(i,k,irvc) = ( rvc_seri(i,k) - qx(i,k,irvc) ) / phys_tstep
           ENDIF
 
@@ -6831,5 +6861,4 @@
              d_qx(i,k,ibs) = ( qbs_seri(i,k) - qx(i,k,ibs) ) / phys_tstep
           ENDIF
-
 
        ENDDO
@@ -6926,6 +6955,7 @@
     ql_ancien(:,:) = ql_seri(:,:)
     qs_ancien(:,:) = qs_seri(:,:)
-    qbs_ancien(:,:) = qbs_seri(:,:)
-    rneb_ancien(:,:) = rneb_seri(:,:)
+    qbs_ancien(:,:)= qbs_seri(:,:)
+    cf_ancien(:,:) = cf_seri(:,:)
+    rvc_ancien(:,:)= rvc_seri(:,:)
 #ifdef ISO
     xt_ancien(:,:,:)=xt_seri(:,:,:)
