Index: /LMDZ6/branches/cirrus/DefLists/field_def_lmdz.xml
===================================================================
--- /LMDZ6/branches/cirrus/DefLists/field_def_lmdz.xml	(revision 4950)
+++ /LMDZ6/branches/cirrus/DefLists/field_def_lmdz.xml	(revision 4951)
@@ -156,4 +156,6 @@
         <field id="LWupSFC"    long_name="Upwd. IR rad. at surface"    unit="W/m2" />
         <field id="LWupSFCclr"    long_name="CS Upwd. IR rad. at surface"    unit="W/m2" />
+        <field id="LWupTOA"    long_name="LWup at TOA"    unit="W/m2" />
+        <field id="LWupTOAclr"    long_name="LWup clear sky at TOA"    unit="W/m2" />
         <field id="LWupTOAcleanclr"    long_name="CS clean (no aerosol) Upwd. IR rad. at TOA"    unit="W/m2" />
         <field id="LWdnSFC"    long_name="Down. IR rad. at surface"    unit="W/m2" />
@@ -798,27 +800,37 @@
 	<field id="dqch4"        long_name="H2O due to CH4 oxidation and photolysis" unit="(kg/kg)/s" />
         <field id="pvap"    long_name="pvap intermediary variable" unit="-">pres*ovap*461.5 / (287.04*(1.+ (10.9491/18.0153)*ovap)) </field>
-        <field id="oclr"    long_name="Clear sky total water"    unit="kg/kg" />
-        <field id="ocld"    long_name="Cloudy sky total water"   unit="kg/kg" />
-        <field id="oss"     long_name="ISSR total water"    unit="kg/kg" />
-        <field id="ovc"     long_name="In-cloud vapor"    unit="kg/kg" />
-        <field id="rnebclr"    long_name="Clear sky fraction"    unit="-" />
-        <field id="rnebss"     long_name="ISSR fraction"     unit="-" />
-        <field id="rnebseri"   long_name="Cloud fraction"    unit="-" />
-        <field id="gammass"    long_name="Supersaturation ratio"    unit="-" />
-        <field id="N1ss"       long_name="N1ss"    unit="-" />
-        <field id="N2ss"       long_name="N2ss"    unit="-" />
-        <field id="drnebsub"   long_name="Cloud fraction change because of sublimation"    unit="s-1" />
-        <field id="drnebcon"   long_name="Cloud fraction change because of condensation"   unit="s-1" />
-        <field id="drnebtur"   long_name="Cloud fraction change because of turbulence"     unit="s-1" />
-        <field id="drnebavi"   long_name="Cloud fraction change because of aviation"       unit="s-1" />
-        <field id="qsatl"      long_name="Saturation with respect to liquid water"    unit="-" />
-        <field id="qsats"      long_name="Saturation with respect to liquid ice"      unit="-" />
-        <field id="flightm"    long_name="Aviation meters flown"      unit="m/s" />
-        <field id="flighth2o"  long_name="Aviation H2O emitted"       unit="kg H2O/s" />
-        <field id="fcontrP"    long_name="Gridbox fraction with potential persistent contrails"     unit="-" />
-        <field id="fcontrN"    long_name="Gridbox fraction with potential non-persistent contrails"     unit="-" />
-        <field id="qcontr"     long_name="Contrail qcontr"     unit="Pa" />
-        <field id="qcontr2"     long_name="Contrail qcontr2"     unit="kg/kg" />
-        <field id="Tcontr"     long_name="Contrail Tcontr"     unit="K" />
+        <field id="cfseri"     long_name="Cloud fraction"    unit="-" />
+        <field id="dcfdyn"     long_name="Dynamics cloud fraction tendency"    unit="s-1" />
+        <field id="rvcseri"    long_name="Cloudy water vapor to total water vapor ratio"    unit="-" />
+        <field id="drvcdyn"    long_name="Dynamics cloudy water vapor to total water vapor ratio tendency"    unit="s-1" />
+        <field id="qsub"       long_name="Subsaturated clear sky total water"    unit="kg/kg" />
+        <field id="qissr"      long_name="Supersaturated clear sky total water"    unit="kg/kg" />
+        <field id="qcld"       long_name="Cloudy sky total water"    unit="kg/kg" />
+        <field id="subfra"     long_name="Subsaturated clear sky fraction"    unit="kg/kg" />
+        <field id="issrfra"    long_name="Supersaturated clear sky fraction"    unit="kg/kg" />
+        <field id="gammacond"  long_name="Condensation threshold w.r.t. saturation"    unit="-" />
+        <field id="dcfsub"     long_name="Sublimation cloud fraction tendency"    unit="s-1" />
+        <field id="dcfcon"     long_name="Condensation cloud fraction tendency"    unit="s-1" />
+        <field id="dcfmix"     long_name="Cloud mixing cloud fraction tendency"    unit="s-1" />
+        <field id="dqiadj"     long_name="Temperature adjustment ice tendency"    unit="kg/kg/s" />
+        <field id="dqisub"     long_name="Sublimation ice tendency"    unit="kg/kg/s" />
+        <field id="dqicon"     long_name="Condensation ice tendency"    unit="kg/kg/s" />
+        <field id="dqimix"     long_name="Cloud mixing ice tendency"    unit="kg/kg/s" />
+        <field id="dqvcadj"    long_name="Temperature adjustment cloudy water vapor tendency"    unit="kg/kg/s" />
+        <field id="dqvcsub"    long_name="Sublimation cloudy water vapor tendency"    unit="kg/kg/s" />
+        <field id="dqvccon"    long_name="Condensation cloudy water vapor tendency"    unit="kg/kg/s" />
+        <field id="dqvcmix"    long_name="Cloud mixing cloudy water vapor tendency"    unit="kg/kg/s" />
+        <field id="qsatl"      long_name="Saturation with respect to liquid"    unit="kg/kg" />
+        <field id="qsati"      long_name="Saturation with respect to ice"    unit="kg/kg" />
+        <field id="Tcontr"     long_name="Temperature threshold for contrail formation"    unit="K" />
+        <field id="qcontr"     long_name="Specific humidity threshold for contrail formation"    unit="Pa" />
+        <field id="qcontr2"    long_name="Specific humidity threshold for contrail formation"    unit="kg/kg" />
+        <field id="fcontrN"    long_name="Fraction with non-persistent contrail in clear-sky"    unit="-" />
+        <field id="fcontrP"    long_name="Fraction with persistent contrail in ISSR"    unit="-" />
+        <field id="dcfavi"     long_name="Aviation cloud fraction tendency"    unit="s-1" />
+        <field id="dqiavi"     long_name="Aviation ice tendency"    unit="kg/kg/s" />
+        <field id="dqvcavi"    long_name="Aviation cloudy water vapor tendency"    unit="kg/kg/s" />
+        <field id="flightdist" long_name="Aviation flown distance"    unit="m/s/mesh" />
+        <field id="flighth2o"  long_name="Aviation H2O flight emission"    unit="kg H2O/s/mesh" />
 	<field id="fluxt"     long_name="flux h"     unit="W/m2" />
         <field id="fluxq"     long_name="flux q"     unit="-" />
Index: /LMDZ6/branches/cirrus/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 4951)
@@ -49,5 +49,6 @@
     prw_ancien, u10m,v10m, treedrg, u_ancien, v_ancien, wake_delta_pbl_TKE, wake_dens, &
     ale_bl, ale_bl_trig, alp_bl, &
-    ale_wake, ale_bl_stat, AWAKE_S
+    ale_wake, ale_bl_stat, AWAKE_S, &
+    cf_ancien, rvc_ancien
 
   USE comconst_mod, ONLY: pi, dtvr
@@ -94,5 +95,5 @@
   USE init_ssrf_m, ONLY: start_init_subsurf
   USE phys_state_var_mod, ONLY: beta_aridity, delta_tsurf, awake_dens, cv_gen, &
-       ratqs_inter_, rneb_ancien
+       ratqs_inter_
   !use ioipsl_getincom
   IMPLICIT NONE
@@ -242,4 +243,7 @@
   ale_wake=0.
   ale_bl_stat=0.
+
+  cf_ancien = 0.
+  rvc_ancien = 0.
   
   z0m(:,:)=0 ! ym missing 5th subsurface initialization
@@ -287,5 +291,4 @@
 
   ratqs_inter_ = 0.002
-  rneb_ancien = 0.
   CALL phyredem( "startphy.nc" )
 
Index: /LMDZ6/branches/cirrus/libf/misc/readTracFiles_mod.f90
===================================================================
--- /LMDZ6/branches/cirrus/libf/misc/readTracFiles_mod.f90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/misc/readTracFiles_mod.f90	(revision 4951)
@@ -114,9 +114,9 @@
   !--- SOME PARAMETERS THAT ARE NOT LIKELY TO CHANGE OFTEN
   CHARACTER(LEN=maxlen), SAVE      :: tran0        = 'air'      !--- Default transporting fluid
-  CHARACTER(LEN=maxlen), PARAMETER :: old_phases   = 'vlirb'     !--- Old phases for water (no separator)
-  CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsrb'     !--- Known phases initials
+  CHARACTER(LEN=maxlen), PARAMETER :: old_phases   = 'vlifbc'     !--- Old phases for water (no separator)
+  CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsfbc'     !--- Known phases initials
   INTEGER, PARAMETER :: nphases = LEN_TRIM(known_phases)        !--- Number of phases
   CHARACTER(LEN=maxlen), SAVE      :: phases_names(nphases) &   !--- Known phases names
-                                = ['gaseous', 'liquid ', 'solid  ', 'cloud  ','blosno ']
+                                = ['gaseous', 'liquid ', 'solid  ', 'fracld ', 'blosnow', 'cldvapr']
   CHARACTER(LEN=1), SAVE :: phases_sep  =  '_'                  !--- Phase separator
   LOGICAL,          SAVE :: tracs_merge = .TRUE.                !--- Merge/stack tracers lists
Index: /LMDZ6/branches/cirrus/libf/phylmd/clesphys.h
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/clesphys.h	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/phylmd/clesphys.h	(revision 4951)
@@ -94,5 +94,5 @@
        LOGICAL :: ok_airs
        INTEGER :: ip_ebil_phy, iflag_rrtm, iflag_ice_thermo, NSW, iflag_albedo
-       LOGICAL :: ok_ice_sursat, ok_plane_h2o, ok_plane_contrail
+       LOGICAL :: ok_ice_supersat, ok_plane_h2o, ok_plane_contrail
        LOGICAL :: ok_chlorophyll
        LOGICAL :: ok_strato
@@ -154,5 +154,5 @@
      &     , ok_lic_melt, ok_lic_cond, aer_type                         &
      &     , iflag_rrtm, ok_strato,ok_hines, ok_qch4                    &
-     &     , iflag_ice_thermo, ok_ice_sursat                            & 
+     &     , iflag_ice_thermo, ok_ice_supersat                          & 
      &     , ok_plane_h2o, ok_plane_contrail                            & 
      &     , ok_gwd_rando, NSW, iflag_albedo                            &
Index: /LMDZ6/branches/cirrus/libf/phylmd/conf_phys_m.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/conf_phys_m.F90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/phylmd/conf_phys_m.F90	(revision 4951)
@@ -173,5 +173,5 @@
     INTEGER,SAVE :: iflag_clw_omp
     INTEGER,SAVE :: iflag_ice_thermo_omp
-    LOGICAL,SAVE :: ok_ice_sursat_omp
+    LOGICAL,SAVE :: ok_ice_supersat_omp
     LOGICAL,SAVE :: ok_plane_h2o_omp, ok_plane_contrail_omp
     REAL,SAVE :: rad_froid_omp, rad_chau1_omp, rad_chau2_omp
@@ -1347,15 +1347,15 @@
 
     !
-    !Config Key  = ok_ice_sursat
-    !Config Desc =
-    !Config Def  = 0
-    !Config Help =
-    !
-    ok_ice_sursat_omp = 0
-    CALL getin('ok_ice_sursat',ok_ice_sursat_omp)
+    !Config Key  = ok_ice_supersat
+    !Config Desc = include ice supersaturation for cold clouds
+    !Config Def  = .FALSE.
+    !Config Help =
+    !
+    ok_ice_supersat_omp = .FALSE.
+    CALL getin('ok_ice_supersat',ok_ice_supersat_omp)
 
     !Config Key  = ok_plane_h2o
-    !Config Desc =
-    !Config Def  = 0
+    !Config Desc = include H2O emissions from aviation
+    !Config Def  = .FALSE.
     !Config Help =
     !
@@ -1364,6 +1364,6 @@
 
     !Config Key  = ok_plane_contrail
-    !Config Desc =
-    !Config Def  = 0
+    !Config Desc = include the formation of contrail cirrus clouds
+    !Config Def  = .FALSE.
     !Config Help =
     !
@@ -2345,5 +2345,5 @@
     rad_chau2 = rad_chau2_omp
     iflag_ice_thermo = iflag_ice_thermo_omp
-    ok_ice_sursat = ok_ice_sursat_omp
+    ok_ice_supersat = ok_ice_supersat_omp
     ok_plane_h2o = ok_plane_h2o_omp
     ok_plane_contrail = ok_plane_contrail_omp
@@ -2770,5 +2770,5 @@
     WRITE(lunout,*) ' rad_chau2 = ',rad_chau2
     WRITE(lunout,*) ' iflag_ice_thermo = ',iflag_ice_thermo
-    WRITE(lunout,*) ' ok_ice_sursat = ',ok_ice_sursat
+    WRITE(lunout,*) ' ok_ice_supersat = ',ok_ice_supersat
     WRITE(lunout,*) ' ok_plane_h2o = ',ok_plane_h2o
     WRITE(lunout,*) ' ok_plane_contrail = ',ok_plane_contrail
Index: /LMDZ6/branches/cirrus/libf/phylmd/create_etat0_unstruct_mod.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/create_etat0_unstruct_mod.F90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/phylmd/create_etat0_unstruct_mod.F90	(revision 4951)
@@ -258,4 +258,8 @@
     prw_ancien = 0.
     agesno = 0.
+    
+    ! LSCP condensation and ice supersaturation
+    cf_ancien = 0.
+    rvc_ancien = 0.
 
     wake_delta_pbl_TKE(:,:,:)=0
@@ -310,5 +314,4 @@
     END IF
     ratqs_inter_ = 0.002
-    rneb_ancien = 0.
 
     CALL gather_omp(cell_area,cell_area_mpi)
Index: DZ6/branches/cirrus/libf/phylmd/ice_sursat_mod.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/ice_sursat_mod.F90	(revision 4950)
+++ 	(revision )
@@ -1,906 +1,0 @@
-MODULE ice_sursat_mod
-
-IMPLICIT NONE
-
-!--flight inventories
-!
-REAL, SAVE, ALLOCATABLE :: flight_m(:,:)    !--flown distance m s-1 per cell
-!$OMP THREADPRIVATE(flight_m)
-REAL, SAVE, ALLOCATABLE :: flight_h2o(:,:)  !--emitted kg H2O s-1 per cell
-!$OMP THREADPRIVATE(flight_h2o)
-!
-!--Fixed Parameters
-!
-!--safety parameters for ERF function
-REAL, PARAMETER :: erf_lim = 5., eps = 1.e-10
-!
-!--Tuning parameters (and their default values)
-!
-!--chi gère la répartition statistique de la longueur des frontières
-!  entre les zones nuages et ISSR/ciel clair sous-saturé. Gamme de valeur :
-!  chi > 1, je n'ai pas regardé de limite max (pour chi = 1, la longueur de
-!  la frontière entre ne nuage et l'ISSR est proportionnelle à la
-!  répartition ISSR/ciel clair sous-sat dans la maille, i.e. il n'y a pas
-!  de favorisation de la localisation de l'ISSR près de nuage. Pour chi = inf, 
-!  le nuage n'est en contact qu'avec de l'ISSR, quelle que soit la taille 
-!  de l'ISSR dans la maille.)
-!
-!--l_turb est la longueur de mélange pour la turbulence. 
-!  dans les tests, ça n'a jamais été modifié pour l'instant. 
-!
-!--tun_N est le paramètre qui contrôle l'importance relative de N_2 par rapport à N_1. 
-!  La valeur est comprise entre 1 et 2 (tun_N = 1 => N_1 = N_2)
-!
-!--tun_ratqs : paramètre qui modifie ratqs en fonction de la valeur de
-!  alpha_cld selon la formule ratqs_new = ratqs_old / ( 1 + tun_ratqs *
-!  alpha_cld ). Dans le rapport il est appelé beta. Il varie entre 0 et 5
-!  (tun_ratqs = 0 => pas de modification de ratqs).
-!
-!--gamma0 and Tgamma: define RHcrit limit above which heterogeneous freezing occurs as a function of T
-!--Karcher and Lohmann (2002)
-!--gamma = 2.583 - t / 207.83
-!--Ren and MacKenzie (2005) reused by Kärcher
-!--gamma = 2.349 - t / 259.0
-!
-!--N_cld: number of clouds in cell (needs to be parametrized at some point)
-!
-!--contrail cross section: typical value found in Freudenthaler et al, GRL, 22, 3501-3504, 1995
-!--in m2, 1000x200 = 200 000 m2 after 15 min
-!
-REAL, SAVE :: chi=1.1, l_turb=50.0, tun_N=1.3, tun_ratqs=3.0
-REAL, SAVE :: gamma0=2.349, Tgamma=259.0, N_cld=100, contrail_cross_section=200000.0
-!$OMP THREADPRIVATE(chi,l_turb,tun_N,tun_ratqs,gamma0,Tgamma,N_cld,contrail_cross_section)
-
-CONTAINS
-
-!*******************************************************************
-!
-SUBROUTINE ice_sursat_init()
-
-  USE print_control_mod, ONLY: lunout
-  USE ioipsl_getin_p_mod, ONLY : getin_p
-
-  IMPLICIT NONE
-
-  CALL getin_p('flag_chi',chi)
-  CALL getin_p('flag_l_turb',l_turb)
-  CALL getin_p('flag_tun_N',tun_N)
-  CALL getin_p('flag_tun_ratqs',tun_ratqs)
-  CALL getin_p('gamma0',gamma0)
-  CALL getin_p('Tgamma',Tgamma)
-  CALL getin_p('N_cld',N_cld)
-  CALL getin_p('contrail_cross_section',contrail_cross_section)
-
-  WRITE(lunout,*) 'Parameters for ice_sursat param'
-  WRITE(lunout,*) 'flag_chi = ', chi
-  WRITE(lunout,*) 'flag_l_turb = ', l_turb
-  WRITE(lunout,*) 'flag_tun_N = ', tun_N
-  WRITE(lunout,*) 'flag_tun_ratqs = ', tun_ratqs
-  WRITE(lunout,*) 'gamma0 = ', gamma0
-  WRITE(lunout,*) 'Tgamma = ', Tgamma
-  WRITE(lunout,*) 'N_cld = ', N_cld
-  WRITE(lunout,*) 'contrail_cross_section = ', contrail_cross_section
-
-END SUBROUTINE ice_sursat_init
-
-!*******************************************************************
-!
-SUBROUTINE airplane(debut,pphis,pplay,paprs,t_seri)
-
-  USE dimphy
-  USE mod_grid_phy_lmdz,  ONLY: klon_glo
-  USE geometry_mod, ONLY: cell_area
-  USE phys_cal_mod, ONLY : mth_cur
-  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
-  USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
-  USE mod_phys_lmdz_para, ONLY: scatter, bcast
-  USE print_control_mod, ONLY: lunout
-
-  IMPLICIT NONE
-
-  INCLUDE "YOMCST.h"
-  INCLUDE 'netcdf.inc'
-
-  !--------------------------------------------------------
-  !--input variables 
-  !--------------------------------------------------------
-  LOGICAL, INTENT(IN) :: debut
-  REAL, INTENT(IN)    :: pphis(klon), pplay(klon,klev), paprs(klon,klev+1), t_seri(klon,klev)
-
-  !--------------------------------------------------------
-  !	... Local variables
-  !--------------------------------------------------------
-
-  CHARACTER (LEN=20) :: modname='airplane_mod'
-  INTEGER :: i, k, kori, iret, varid, error, ncida, klona
-  INTEGER,SAVE :: nleva, ntimea
-!$OMP THREADPRIVATE(nleva,ntimea)
-  REAL, ALLOCATABLE :: pkm_airpl_glo(:,:,:)    !--km/s
-  REAL, ALLOCATABLE :: ph2o_airpl_glo(:,:,:)   !--molec H2O/cm3/s
-  REAL, ALLOCATABLE, SAVE :: zmida(:), zinta(:)
-  REAL, ALLOCATABLE, SAVE :: pkm_airpl(:,:,:)
-  REAL, ALLOCATABLE, SAVE :: ph2o_airpl(:,:,:)
-!$OMP THREADPRIVATE(pkm_airpl,ph2o_airpl,zmida,zinta)
-  REAL :: zalt(klon,klev+1)
-  REAL :: zrho, zdz(klon,klev), zfrac
-
-  !
-  IF (debut) THEN
-  !--------------------------------------------------------------------------------
-  !       ... Open the file and read airplane emissions
-  !--------------------------------------------------------------------------------
-  !
-  IF (is_mpi_root .AND. is_omp_root) THEN
-      !
-      iret = nf_open('aircraft_phy.nc', 0, ncida)
-      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to open aircraft_phy.nc file',1)
-      ! ... Get lengths
-      iret = nf_inq_dimid(ncida, 'time', varid)
-      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimid in aircraft_phy.nc file',1)
-      iret = nf_inq_dimlen(ncida, varid, ntimea)
-      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimlen aircraft_phy.nc file',1)
-      iret = nf_inq_dimid(ncida, 'vector', varid)
-      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimid aircraft_phy.nc file',1)
-      iret = nf_inq_dimlen(ncida, varid, klona)
-      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimlen aircraft_phy.nc file',1)
-      iret = nf_inq_dimid(ncida, 'lev', varid)
-      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1)
-      iret = nf_inq_dimlen(ncida, varid, nleva)
-      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimlen aircraft_phy.nc file',1)
-      !
-      IF ( klona /= klon_glo ) THEN
-        WRITE(lunout,*) 'klona & klon_glo =', klona, klon_glo
-        CALL abort_physic(modname,'problem klon in aircraft_phy.nc file',1)
-      ENDIF
-      !
-      IF ( ntimea /= 12 ) THEN
-        WRITE(lunout,*) 'ntimea=', ntimea
-        CALL abort_physic(modname,'problem ntime<>12 in aircraft_phy.nc file',1)
-      ENDIF
-      !
-      ALLOCATE(zmida(nleva), STAT=error)
-      IF (error /= 0) CALL abort_physic(modname,'problem to allocate zmida',1)
-      ALLOCATE(pkm_airpl_glo(klona,nleva,ntimea), STAT=error)
-      IF (error /= 0) CALL abort_physic(modname,'problem to allocate pkm_airpl_glo',1)
-      ALLOCATE(ph2o_airpl_glo(klona,nleva,ntimea), STAT=error)
-      IF (error /= 0) CALL abort_physic(modname,'problem to allocate ph2o_airpl_glo',1)
-      !
-      iret = nf_inq_varid(ncida, 'lev', varid)
-      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1)
-      iret = nf_get_var_double(ncida, varid, zmida)
-      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read zmida file',1)
-      !
-      iret = nf_inq_varid(ncida, 'emi_co2_aircraft', varid)  !--CO2 as a proxy for m flown - 
-      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_distance dimid aircraft_phy.nc file',1)
-      iret = nf_get_var_double(ncida, varid, pkm_airpl_glo)
-      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read pkm_airpl file',1)
-      !
-      iret = nf_inq_varid(ncida, 'emi_h2o_aircraft', varid)
-      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_h2o_aircraft dimid aircraft_phy.nc file',1)
-      iret = nf_get_var_double(ncida, varid, ph2o_airpl_glo)
-      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read ph2o_airpl file',1)
-      ! 
-     ENDIF    !--is_mpi_root and is_omp_root
-     !
-     CALL bcast(nleva)
-     CALL bcast(ntimea)
-     ! 
-     IF (.NOT.ALLOCATED(zmida)) ALLOCATE(zmida(nleva), STAT=error)
-     IF (.NOT.ALLOCATED(zinta)) ALLOCATE(zinta(nleva+1), STAT=error)
-     !
-     ALLOCATE(pkm_airpl(klon,nleva,ntimea))
-     ALLOCATE(ph2o_airpl(klon,nleva,ntimea))
-     !
-     ALLOCATE(flight_m(klon,klev))
-     ALLOCATE(flight_h2o(klon,klev))
-     !
-     CALL bcast(zmida)
-     zinta(1)=0.0                                   !--surface
-     DO k=2, nleva
-       zinta(k) = (zmida(k-1)+zmida(k))/2.0*1000.0  !--conversion from km to m
-     ENDDO
-     zinta(nleva+1)=zinta(nleva)+(zmida(nleva)-zmida(nleva-1))*1000.0 !--extrapolation for last interface
-     !print *,'zinta=', zinta
-     !
-     CALL scatter(pkm_airpl_glo,pkm_airpl)
-     CALL scatter(ph2o_airpl_glo,ph2o_airpl)
-     !
-!$OMP MASTER
-     IF (is_mpi_root .AND. is_omp_root) THEN
-        DEALLOCATE(pkm_airpl_glo)
-        DEALLOCATE(ph2o_airpl_glo)
-     ENDIF   !--is_mpi_root
-!$OMP END MASTER
-
-  ENDIF !--debut
-!
-!--compute altitude of model level interfaces
-!
-  DO i = 1, klon
-    zalt(i,1)=pphis(i)/RG         !--in m
-  ENDDO
-!
-  DO k=1, klev
-    DO i = 1, klon
-      zrho=pplay(i,k)/t_seri(i,k)/RD
-      zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho/RG
-      zalt(i,k+1)=zalt(i,k)+zdz(i,k)   !--in m
-    ENDDO
-  ENDDO
-!
-!--vertical reprojection 
-!
-  flight_m(:,:)=0.0 
-  flight_h2o(:,:)=0.0 
-!
-  DO k=1, klev
-    DO kori=1, nleva 
-      DO i=1, klon
-        !--fraction of layer kori included in layer k
-        zfrac=max(0.0,min(zalt(i,k+1),zinta(kori+1))-max(zalt(i,k),zinta(kori)))/(zinta(kori+1)-zinta(kori))
-        !--reproject
-        flight_m(i,k)=flight_m(i,k) + pkm_airpl(i,kori,mth_cur)*zfrac 
-        !--reproject
-        flight_h2o(i,k)=flight_h2o(i,k) + ph2o_airpl(i,kori,mth_cur)*zfrac   
-      ENDDO
-    ENDDO
-  ENDDO
-!
-  DO k=1, klev
-    DO i=1, klon
-      !--molec.cm-3.s-1 / (molec/mol) * kg CO2/mol * m2 * m * cm3/m3 / (kg CO2/m) => m s-1 per cell
-      flight_m(i,k)=flight_m(i,k)/RNAVO*44.e-3*cell_area(i)*zdz(i,k)*1.e6/16.37e-3
-      flight_m(i,k)=flight_m(i,k)*100.0  !--x100 to augment signal to noise 
-      !--molec.cm-3.s-1 / (molec/mol) * kg H2O/mol * m2 * m * cm3/m3 => kg H2O s-1 per cell
-      flight_h2o(i,k)=flight_h2o(i,k)/RNAVO*18.e-3*cell_area(i)*zdz(i,k)*1.e6
-    ENDDO
-  ENDDO
-!
-END SUBROUTINE airplane
-
-!********************************************************************
-! simple routine to initialise flight_m and test a flight corridor
-!--Olivier Boucher - 2021 
-!
-SUBROUTINE flight_init()
-  USE dimphy
-  USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
-  IMPLICIT NONE
-  INTEGER :: i
-
-  ALLOCATE(flight_m(klon,klev))
-  ALLOCATE(flight_h2o(klon,klev))
-  !
-  flight_m(:,:) = 0.0    !--initialisation
-  flight_h2o(:,:) = 0.0  !--initialisation
-  !
-  DO i=1, klon
-   IF (latitude_deg(i).GE.42.0.AND.latitude_deg(i).LE.48.0) THEN 
-     flight_m(i,38) = 50000.0  !--5000 m of flight/second in grid cell x 10 scaling
-   ENDIF
-  ENDDO
-  
-  RETURN
-END SUBROUTINE flight_init
-
-!*******************************************************************
-!--Routine to deal with ice supersaturation
-!--Determines the respective fractions of unsaturated clear sky, ice supersaturated clear sky and cloudy sky
-!--Diagnoses regions prone for non-persistent and persistent contrail formation
-!
-!--Audran Borella - 2021 
-!
-SUBROUTINE ice_sursat(pplay, dpaprs, dtime, i, k, t, q, gamma_ss, &
-                      qsat, t_actuel, rneb_seri, ratqs, rneb, qincld,   &
-                      rnebss, qss, Tcontr, qcontr, qcontr2, fcontrN, fcontrP)
-  !
-  USE dimphy
-  USE print_control_mod,    ONLY: prt_level, lunout
-  USE phys_state_var_mod,   ONLY: pbl_tke, t_ancien 
-  USE phys_local_var_mod,   ONLY: N1_ss, N2_ss
-  USE phys_local_var_mod,   ONLY: drneb_sub, drneb_con, drneb_tur, drneb_avi
-!!  USE phys_local_var_mod,   ONLY: Tcontr, qcontr, fcontrN, fcontrP
-  USE indice_sol_mod,       ONLY: is_ave 
-  USE geometry_mod,         ONLY: cell_area
-  !
-  IMPLICIT NONE
-  INCLUDE "YOMCST.h"
-  INCLUDE "YOETHF.h"
-  INCLUDE "FCTTRE.h"
-  INCLUDE "clesphys.h"
-
-  !
-  ! Input
-  ! Beware: this routine works on a gridpoint!
-  !
-  REAL,     INTENT(IN)    :: pplay     ! layer pressure (Pa)
-  REAL,     INTENT(IN)    :: dpaprs    ! layer delta pressure (Pa)
-  REAL,     INTENT(IN)    :: dtime     ! intervalle du temps (s)
-  REAL,     INTENT(IN)    :: t         ! température advectée (K)
-  REAL,     INTENT(IN)    :: qsat      ! vapeur de saturation
-  REAL,     INTENT(IN)    :: t_actuel  ! temperature actuelle de la maille (K)
-  REAL,     INTENT(IN)    :: rneb_seri ! fraction nuageuse en memoire
-  INTEGER,  INTENT(IN)    :: i, k
-  !
-  !  Input/output
-  !
-  REAL,     INTENT(INOUT) :: q         ! vapeur de la maille (=zq)
-  REAL,     INTENT(INOUT) :: ratqs     ! determine la largeur de distribution de vapeur
-  REAL,     INTENT(INOUT) :: Tcontr, qcontr, qcontr2, fcontrN, fcontrP
-  !
-  !  Output
-  !
-  REAL,     INTENT(OUT)   :: gamma_ss  !
-  REAL,     INTENT(OUT)   :: rneb      !  cloud fraction
-  REAL,     INTENT(OUT)   :: qincld    !  in-cloud total water
-  REAL,     INTENT(OUT)   :: rnebss    !  ISSR fraction
-  REAL,     INTENT(OUT)   :: qss       !  in-ISSR total water
-  !
-  ! Local
-  !
-  REAL PI
-  PARAMETER (PI=4.*ATAN(1.))
-  REAL rnebclr, gamma_prec
-  REAL qclr, qvc, qcld, qi
-  REAL zrho, zdz, zrhodz
-  REAL pdf_N, pdf_N1, pdf_N2
-  REAL pdf_a, pdf_b
-  REAL pdf_e1, pdf_e2, pdf_k
-  REAL drnebss, drnebclr, dqss, dqclr, sum_rneb_rnebss, dqss_avi
-  REAL V_cell !--volume of the cell 
-  REAL M_cell !--dry mass of the cell 
-  REAL tke, sig, L_tur, b_tur, q_eq
-  REAL V_env, V_cld, V_ss, V_clr
-  REAL zcor
-  !
-  !--more local variables for diagnostics
-  !--imported from YOMCST.h 
-  !--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1)
-  !--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air
-  !--values from Schumann, Meteorol Zeitschrift, 1996
-  !--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen
-  !--Qheat = 43.  /  50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen
-  REAL, PARAMETER :: EiH2O=1.25  !--emission index of water vapour for kerosene (kg kg-1)
-  REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1) 
-  REAL, PARAMETER :: eta=0.3     !--average propulsion efficiency of the aircraft
-  !--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute 
-  !--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010.
-  REAL :: Gcontr
-  !--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) 
-  !--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1)
-  REAL :: qsatliqcontr
-
-     ! Initialisations
-     zrho = pplay / t / RD            !--dry density kg m-3
-     zrhodz = dpaprs / RG             !--dry air mass kg m-2
-     zdz = zrhodz / zrho              !--cell thickness m
-     V_cell = zdz * cell_area(i)      !--cell volume m3
-     M_cell = zrhodz * cell_area(i)   !--cell dry air mass kg
-     !
-     ! Recuperation de la memoire sur la couverture nuageuse
-     rneb = rneb_seri
-     !
-     ! Ajout des émissions de H2O dues à l'aviation
-     ! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
-     ! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O )  
-     !      = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) ) 
-     ! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q) 
-     ! flight_h2O is in kg H2O / s / cell
-     !
-     IF (ok_plane_h2o) THEN 
-       q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) ) 
-     ENDIF
-     !
-     !--Estimating gamma 
-     gamma_ss = MAX(1.0, gamma0 - t_actuel/Tgamma)
-     !gamma_prec = MAX(1.0, gamma0 - t_ancien(i,k)/Tgamma)      !--formulation initiale d Audran
-     gamma_prec = MAX(1.0, gamma0 - t/Tgamma)                   !--autre formulation possible basée sur le T du pas de temps
-     !
-     ! Initialisation de qvc : q_sat du pas de temps precedent
-     !qvc = R2ES*FOEEW(t_ancien(i,k),1.)/pplay      !--formulation initiale d Audran
-     qvc = R2ES*FOEEW(t,1.)/pplay                   !--autre formulation possible basée sur le T du pas de temps
-     qvc = min(0.5,qvc)
-     zcor = 1./(1.-RETV*qvc)
-     qvc = qvc*zcor
-     !
-     ! Modification de ratqs selon formule proposee : ksi_new = ksi_old/(1+beta*alpha_cld)
-     ratqs = ratqs / (tun_ratqs*rneb_seri + 1.)
-     !
-     ! Calcul de N
-     pdf_k = -sqrt(log(1.+ratqs**2.))
-     pdf_a = log(qvc/q)/(pdf_k*sqrt(2.))
-     pdf_b = pdf_k/(2.*sqrt(2.))
-     pdf_e1 = pdf_a+pdf_b
-     IF (abs(pdf_e1).GE.erf_lim) THEN
-        pdf_e1 = sign(1.,pdf_e1)
-        pdf_N = max(0.,sign(rneb,pdf_e1))
-     ELSE
-        pdf_e1 = erf(pdf_e1)
-        pdf_e1 = 0.5*(1.+pdf_e1)
-        pdf_N = max(0.,rneb/pdf_e1)
-     ENDIF
-     !
-     ! On calcule ensuite N1 et N2. Il y a deux cas : N > 1 et N <= 1
-     ! Cas 1 : N > 1. N'arrive en theorie jamais, c'est une barriere
-     ! On perd la memoire sur la temperature (sur qvc) pour garder
-     ! celle sur alpha_cld
-     IF (pdf_N.GT.1.) THEN
-        ! On inverse alpha_cld = int_qvc^infty P(q) dq
-        ! pour determiner qvc = f(alpha_cld)
-        ! On approxime en serie entiere erf-1(x)
-        qvc = 2.*rneb-1.
-        qvc = qvc + PI/12.*qvc**3 + 7.*PI**2/480.*qvc**5 &
-             + 127.*PI**3/40320.*qvc**7 + 4369.*PI**4/5806080.*qvc**9 &
-             + 34807.*PI**5/182476800.*qvc**11
-        qvc = sqrt(PI)/2.*qvc
-        qvc = (qvc-pdf_b)*pdf_k*sqrt(2.)
-        qvc = q*exp(qvc)
-
-        ! On met a jour rneb avec la nouvelle valeur de qvc
-        ! La maj est necessaire a cause de la serie entiere approximative
-        pdf_a = log(qvc/q)/(pdf_k*sqrt(2.))
-        pdf_e1 = pdf_a+pdf_b
-        IF (abs(pdf_e1).GE.erf_lim) THEN
-           pdf_e1 = sign(1.,pdf_e1)
-        ELSE
-           pdf_e1 = erf(pdf_e1)
-        ENDIF
-        pdf_e1 = 0.5*(1.+pdf_e1)
-        rneb = pdf_e1
-        
-        ! Si N > 1, N1 et N2 = 1
-        pdf_N1 = 1.
-        pdf_N2 = 1.
-        
-     ! Cas 2 : N <= 1
-     ELSE
-        ! D'abord on calcule N2 avec le tuning
-        pdf_N2 = min(1.,pdf_N*tun_N)
-        
-        ! Puis N1 pour assurer la conservation de alpha_cld
-        pdf_a = log(qvc*gamma_prec/q)/(pdf_k*sqrt(2.))
-        pdf_e2 = pdf_a+pdf_b
-        IF (abs(pdf_e2).GE.erf_lim) THEN
-           pdf_e2 = sign(1.,pdf_e2)
-        ELSE
-           pdf_e2 = erf(pdf_e2)
-        ENDIF
-        pdf_e2 = 0.5*(1.+pdf_e2) ! integrale sous P pour q > gamma qsat
-
-        IF (abs(pdf_e1-pdf_e2).LT.eps) THEN
-           pdf_N1 = pdf_N2
-        ELSE
-           pdf_N1 = (rneb-pdf_N2*pdf_e2)/(pdf_e1-pdf_e2)
-        ENDIF
-
-        ! Barriere qui traite le cas gamma_prec = 1.
-        IF (pdf_N1.LE.0.) THEN
-           pdf_N1 = 0.
-           IF (pdf_e2.GT.eps) THEN
-              pdf_N2 = rneb/pdf_e2
-           ELSE
-              pdf_N2 = 0.
-           ENDIF
-        ENDIF
-     ENDIF
-
-     ! Physique 1
-     ! Sublimation
-     IF (qvc.LT.qsat) THEN
-        pdf_a = log(qvc/q)/(pdf_k*sqrt(2.))
-        pdf_e1 = pdf_a+pdf_b
-        IF (abs(pdf_e1).GE.erf_lim) THEN
-           pdf_e1 = sign(1.,pdf_e1)
-        ELSE
-           pdf_e1 = erf(pdf_e1)
-        ENDIF
-
-        pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
-        pdf_e2 = pdf_a+pdf_b
-        IF (abs(pdf_e2).GE.erf_lim) THEN
-           pdf_e2 = sign(1.,pdf_e2)
-        ELSE
-           pdf_e2 = erf(pdf_e2)
-        ENDIF
-
-        pdf_e1 = 0.5*pdf_N1*(pdf_e1-pdf_e2)
-        
-        ! Calcul et ajout de la tendance
-        drneb_sub(i,k) = - pdf_e1/dtime    !--unit [s-1]
-        rneb = rneb + drneb_sub(i,k)*dtime
-     ELSE
-        drneb_sub(i,k) = 0.
-     ENDIF
-     
-     ! NOTE : verifier si ca marche bien pour gamma proche de 1.
-
-     ! Condensation
-     IF (gamma_ss*qsat.LT.gamma_prec*qvc) THEN
-     
-        pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.))
-        pdf_e1 = pdf_a+pdf_b
-        IF (abs(pdf_e1).GE.erf_lim) THEN
-           pdf_e1 = sign(1.,pdf_e1)
-        ELSE
-           pdf_e1 = erf(pdf_e1)
-        ENDIF
-
-        pdf_a = log(gamma_prec*qvc/q)/(pdf_k*sqrt(2.))
-        pdf_e2 = pdf_a+pdf_b
-        IF (abs(pdf_e2).GE.erf_lim) THEN
-           pdf_e2 = sign(1.,pdf_e2)
-        ELSE
-           pdf_e2 = erf(pdf_e2)
-        ENDIF
-
-        pdf_e1 = 0.5*(1.-pdf_N1)*(pdf_e1-pdf_e2)
-        pdf_e2 = 0.5*(1.-pdf_N2)*(1.+pdf_e2)
-
-        ! Calcul et ajout de la tendance
-        drneb_con(i,k) = (pdf_e1 + pdf_e2)/dtime         !--unit [s-1]
-        rneb = rneb + drneb_con(i,k)*dtime
-        
-     ELSE
-     
-        pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.))
-        pdf_e1 = pdf_a+pdf_b
-        IF (abs(pdf_e1).GE.erf_lim) THEN
-           pdf_e1 = sign(1.,pdf_e1)
-        ELSE
-           pdf_e1 = erf(pdf_e1)
-        ENDIF
-        pdf_e1 = 0.5*(1.-pdf_N2)*(1.+pdf_e1)
-
-        ! Calcul et ajout de la tendance
-        drneb_con(i,k) = pdf_e1/dtime         !--unit [s-1]
-        rneb = rneb + drneb_con(i,k)*dtime
-        
-     ENDIF
-
-     ! Calcul des grandeurs diagnostiques
-     ! Determination des grandeurs ciel clair
-     pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
-     pdf_e1 = pdf_a+pdf_b
-     IF (abs(pdf_e1).GE.erf_lim) THEN
-        pdf_e1 = sign(1.,pdf_e1)
-     ELSE
-        pdf_e1 = erf(pdf_e1)
-     ENDIF
-     pdf_e1 = 0.5*(1.-pdf_e1)
-
-     pdf_e2 = pdf_a-pdf_b
-     IF (abs(pdf_e2).GE.erf_lim) THEN
-        pdf_e2 = sign(1.,pdf_e2)
-     ELSE
-        pdf_e2 = erf(pdf_e2)
-     ENDIF
-     pdf_e2 = 0.5*q*(1.-pdf_e2)
-
-     rnebclr = pdf_e1
-     qclr = pdf_e2
-
-     ! Determination de q_cld
-     ! Partie 1
-     pdf_a = log(max(qsat,qvc)/q)/(pdf_k*sqrt(2.))
-     pdf_e1 = pdf_a-pdf_b
-     IF (abs(pdf_e1).GE.erf_lim) THEN
-        pdf_e1 = sign(1.,pdf_e1)
-     ELSE
-        pdf_e1 = erf(pdf_e1)
-     ENDIF
-
-     pdf_a = log(min(gamma_ss*qsat,gamma_prec*qvc)/q)/(pdf_k*sqrt(2.))
-     pdf_e2 = pdf_a-pdf_b
-     IF (abs(pdf_e2).GE.erf_lim) THEN
-        pdf_e2 = sign(1.,pdf_e2)
-     ELSE
-        pdf_e2 = erf(pdf_e2)
-     ENDIF
-
-     pdf_e1 = 0.5*q*pdf_N1*(pdf_e1-pdf_e2)
-     
-     qcld = pdf_e1
-
-     ! Partie 2 (sous condition)
-     IF (gamma_ss*qsat.GT.gamma_prec*qvc) THEN
-        pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.))
-        pdf_e1 = pdf_a-pdf_b
-        IF (abs(pdf_e1).GE.erf_lim) THEN
-           pdf_e1 = sign(1.,pdf_e1)
-        ELSE
-           pdf_e1 = erf(pdf_e1)
-        ENDIF
-
-        pdf_e2 = 0.5*q*pdf_N2*(pdf_e2-pdf_e1)
-
-        qcld = qcld + pdf_e2
-
-        pdf_e2 = pdf_e1  
-     ENDIF
-
-     ! Partie 3
-     pdf_e2 = 0.5*q*(1.+pdf_e2)
-     
-     qcld = qcld + pdf_e2
-
-     ! Fin du calcul de q_cld
-     
-     ! Determination des grandeurs ISSR via les equations de conservation
-     rneb=MIN(rneb, 1. - rnebclr - eps)      !--ajout OB - barrière
-     rnebss = MAX(0.0, 1. - rnebclr - rneb)  !--ajout OB
-     qss = MAX(0.0, q - qclr - qcld)         !--ajout OB
-
-     ! Physique 2 : Turbulence
-     IF (rneb.GT.eps.AND.rneb.LT.1.-eps) THEN ! rneb != 0 and != 1
-       !
-       tke = pbl_tke(i,k,is_ave)
-       ! A MODIFIER formule a revoir
-       L_tur = min(l_turb, sqrt(tke)*dtime)
-
-       ! On fait pour l'instant l'hypothese a = 3b. V = 4/3 pi a b**2 = alpha_cld
-       ! donc b = alpha_cld/4pi **1/3. 
-       b_tur = (rneb*V_cell/4./PI/N_cld)**(1./3.)
-       ! On verifie que la longeur de melange n'est pas trop grande
-       IF (L_tur.GT.b_tur) THEN
-          L_tur = b_tur
-       ENDIF
-       
-       V_env = N_cld*4.*PI*(3.*(b_tur**2.)*L_tur + L_tur**3. + 3.*b_tur*(L_tur**2.))
-       V_cld = N_cld*4.*PI*(3.*(b_tur**2.)*L_tur + L_tur**3. - 3.*b_tur*(L_tur**2.))
-       V_cld = abs(V_cld)
-
-       ! Repartition statistique des zones de contact nuage-ISSR et nuage-ciel clair
-       sig = rnebss/(chi*rnebclr+rnebss)
-       V_ss = MIN(sig*V_env,rnebss*V_cell)
-       V_clr = MIN((1.-sig)*V_env,rnebclr*V_cell)
-       V_cld = MIN(V_cld,rneb*V_cell)
-       V_env = V_ss + V_clr
-
-       ! ISSR => rneb
-       drnebss = -1.*V_ss/V_cell
-       dqss = drnebss*qss/MAX(eps,rnebss)
-
-       ! Clear sky <=> rneb
-       q_eq = V_env*qclr/MAX(eps,rnebclr) + V_cld*qcld/MAX(eps,rneb)
-       q_eq = q_eq/(V_env + V_cld)
-
-       IF (q_eq.GT.qsat) THEN
-          drnebclr = - V_clr/V_cell
-          dqclr = drnebclr*qclr/MAX(eps,rnebclr)
-       ELSE
-          drnebclr = V_cld/V_cell
-          dqclr = drnebclr*qcld/MAX(eps,rneb)
-       ENDIF
-
-       ! Maj des variables avec les tendances
-       rnebclr = MAX(0.0,rnebclr + drnebclr)   !--OB ajout d'un max avec eps (il faudrait modified drnebclr pour le diag)
-       qclr = MAX(0.0, qclr + dqclr)           !--OB ajout d'un max avec 0
-
-       rneb = rneb - drnebclr - drnebss
-       qcld = qcld - dqclr - dqss
-
-       rnebss = MAX(0.0,rnebss + drnebss)     !--OB ajout d'un max avec eps (il faudrait modifier drnebss pour le diag)
-       qss = MAX(0.0, qss + dqss)             !--OB ajout d'un max avec 0
-
-       ! Tendances pour le diagnostic
-       drneb_tur(i,k) = (drnebclr + drnebss)/dtime  !--unit [s-1]
-
-     ENDIF ! rneb
-
-     !--add a source of cirrus from aviation contrails
-     IF (ok_plane_contrail) THEN
-       drneb_avi(i,k) = rnebss*flight_m(i,k)*contrail_cross_section/V_cell    !--tendency rneb due to aviation [s-1]
-       drneb_avi(i,k) = MIN(drneb_avi(i,k), rnebss/dtime)                     !--majoration
-       dqss_avi = qss*drneb_avi(i,k)/MAX(eps,rnebss)                          !--tendency q aviation [kg kg-1 s-1]
-       rneb = rneb + drneb_avi(i,k)*dtime                                     !--add tendency to rneb
-       qcld = qcld + dqss_avi*dtime                                           !--add tendency to qcld
-       rnebss = rnebss - drneb_avi(i,k)*dtime                                 !--add tendency to rnebss
-       qss = qss - dqss_avi*dtime                                             !--add tendency to qss
-     ELSE 
-       drneb_avi(i,k)=0.0 
-     ENDIF
-
-     ! Barrieres
-     ! ISSR trop petite
-     IF (rnebss.LT.eps) THEN
-        rneb = MIN(rneb + rnebss,1.0-eps) !--ajout OB barriere
-        qcld = qcld + qss
-        rnebss = 0.
-        qss = 0.
-     ENDIF
-
-     ! le nuage est trop petit
-     IF (rneb.LT.eps) THEN
-        ! s'il y a une ISSR on met tout dans l'ISSR, sinon dans le
-        ! clear sky
-        IF (rnebss.LT.eps) THEN
-           rnebclr = 1.
-           rnebss = 0. !--ajout OB
-           qclr = q
-        ELSE
-           rnebss = MIN(rnebss + rneb,1.0-eps) !--ajout OB barriere
-           qss = qss + qcld
-        ENDIF
-        rneb = 0.
-        qcld = 0.
-        qincld = qsat * gamma_ss
-     ELSE
-        qincld = qcld / rneb
-     ENDIF
-
-     !--OB ajout borne superieure
-     sum_rneb_rnebss=rneb+rnebss
-     rneb=rneb*MIN(1.-eps,sum_rneb_rnebss)/MAX(eps,sum_rneb_rnebss)
-     rnebss=rnebss*MIN(1.-eps,sum_rneb_rnebss)/MAX(eps,sum_rneb_rnebss)
-
-     ! On ecrit dans la memoire
-     N1_ss(i,k) = pdf_N1
-     N2_ss(i,k) = pdf_N2
-   
-     !--Diagnostics only used from last iteration
-     !--test
-     !!Tcontr(i,k)=200.
-     !!fcontrN(i,k)=1.0
-     !!fcontrP(i,k)=0.5
-     !
-     !--slope of dilution line in exhaust
-     !--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1) 
-     Gcontr = EiH2O * RCPD * pplay / (eps_w*Qheat*(1.-eta))             !--Pa K-1
-     !--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
-     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
-         fcontrN = 0.0
-         !
-         !--integral of P(q) from qsat to qcontr in ISSR 
-         pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
-         pdf_e1 = pdf_a+pdf_b
-         IF (abs(pdf_e1).GE.erf_lim) THEN
-            pdf_e1 = sign(1.,pdf_e1)
-         ELSE
-            pdf_e1 = erf(pdf_e1)
-         ENDIF
-         !
-         !!pdf_a = log(MIN(qcontr(i,k),qvc)/q)/(pdf_k*sqrt(2.))
-         pdf_a = log(MIN(qcontr2,qvc)/q)/(pdf_k*sqrt(2.))
-         pdf_e2 = pdf_a+pdf_b
-         IF (abs(pdf_e2).GE.erf_lim) THEN
-            pdf_e2 = sign(1.,pdf_e2)
-         ELSE
-            pdf_e2 = erf(pdf_e2)
-         ENDIF
-         !
-         !!fcontrP(i,k) = MAX(0., 0.5*(pdf_e1-pdf_e2))
-         fcontrP = MAX(0., 0.5*(pdf_e1-pdf_e2))
-         !
-         pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
-         pdf_e1 = pdf_a+pdf_b
-         IF (abs(pdf_e1).GE.erf_lim) THEN
-            pdf_e1 = sign(1.,pdf_e1)
-         ELSE
-            pdf_e1 = erf(pdf_e1)
-         ENDIF
-         !
-         !!pdf_a = log(MIN(qcontr(i,k),qvc)/q)/(pdf_k*sqrt(2.))
-         pdf_a = log(MIN(qcontr2,qvc)/q)/(pdf_k*sqrt(2.))
-         pdf_e2 = pdf_a+pdf_b
-         IF (abs(pdf_e2).GE.erf_lim) THEN
-            pdf_e2 = sign(1.,pdf_e2)
-         ELSE
-            pdf_e2 = erf(pdf_e2)
-         ENDIF
-         !
-         !!fcontrP(i,k) = MAX(0., 0.5*(pdf_e1-pdf_e2))
-         fcontrP = MAX(0., 0.5*(pdf_e1-pdf_e2))
-         !
-         pdf_a = log(MAX(qsat,qvc)/q)/(pdf_k*sqrt(2.))
-         pdf_e1 = pdf_a+pdf_b
-         IF (abs(pdf_e1).GE.erf_lim) THEN
-            pdf_e1 = sign(1.,pdf_e1)
-         ELSE
-            pdf_e1 = erf(pdf_e1)
-         ENDIF
-         !
-         !!pdf_a = log(MIN(qcontr(i,k),MIN(gamma_prec*qvc,gamma_ss*qsat))/q)/(pdf_k*sqrt(2.))
-         pdf_a = log(MIN(qcontr2,MIN(gamma_prec*qvc,gamma_ss*qsat))/q)/(pdf_k*sqrt(2.))
-         pdf_e2 = pdf_a+pdf_b
-         IF (abs(pdf_e2).GE.erf_lim) THEN
-            pdf_e2 = sign(1.,pdf_e2)
-         ELSE
-            pdf_e2 = erf(pdf_e2)
-         ENDIF
-         !
-         !!fcontrP(i,k) = fcontrP(i,k) + MAX(0., 0.5*(1-pdf_N1)*(pdf_e1-pdf_e2))
-         fcontrP = fcontrP + MAX(0., 0.5*(1-pdf_N1)*(pdf_e1-pdf_e2))
-         !
-         pdf_a = log(gamma_prec*qvc/q)/(pdf_k*sqrt(2.))
-         pdf_e1 = pdf_a+pdf_b
-         IF (abs(pdf_e1).GE.erf_lim) THEN
-            pdf_e1 = sign(1.,pdf_e1)
-         ELSE
-            pdf_e1 = erf(pdf_e1)
-         ENDIF
-         !
-         !!pdf_a = log(MIN(qcontr(i,k),gamma_ss*qsat)/q)/(pdf_k*sqrt(2.))
-         pdf_a = log(MIN(qcontr2,gamma_ss*qsat)/q)/(pdf_k*sqrt(2.))
-         pdf_e2 = pdf_a+pdf_b
-         IF (abs(pdf_e2).GE.erf_lim) THEN
-            pdf_e2 = sign(1.,pdf_e2)
-         ELSE
-            pdf_e2 = erf(pdf_e2)
-         ENDIF
-         !
-         !!fcontrP(i,k) = fcontrP(i,k) + MAX(0., 0.5*(1-pdf_N2)*(pdf_e1-pdf_e2))
-         fcontrP = fcontrP + MAX(0., 0.5*(1-pdf_N2)*(pdf_e1-pdf_e2))
-         !
-       ELSE  !--qcontr LT qsat
-         !
-         !--all of ISSR is prone for contrail formation
-         !!fcontrP(i,k) = rnebss
-         fcontrP = rnebss
-         !
-         !--integral of zq from qcontr to qsat in unsaturated clear-sky region
-         !!pdf_a = log(qcontr(i,k)/q)/(pdf_k*sqrt(2.))
-         pdf_a = log(qcontr2/q)/(pdf_k*sqrt(2.))
-         pdf_e1 = pdf_a+pdf_b   !--normalement pdf_b est deja defini
-         IF (abs(pdf_e1).GE.erf_lim) THEN
-            pdf_e1 = sign(1.,pdf_e1)
-         ELSE
-            pdf_e1 = erf(pdf_e1)
-         ENDIF
-         !
-         pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
-         pdf_e2 = pdf_a+pdf_b
-         IF (abs(pdf_e2).GE.erf_lim) THEN
-            pdf_e2 = sign(1.,pdf_e2)
-         ELSE
-            pdf_e2 = erf(pdf_e2)
-         ENDIF
-         !
-         !!fcontrN(i,k) = 0.5*(pdf_e1-pdf_e2)
-         fcontrN = 0.5*(pdf_e1-pdf_e2)
-         !!fcontrN=2.0
-         !
-       ENDIF
-       !
-       ENDIF !-- t < Tcontr
-     !
-     ENDIF !-- Gcontr > 0.1
-
-  RETURN
-END SUBROUTINE ice_sursat
-!
-!*******************************************************************
-!
-END MODULE ice_sursat_mod
Index: /LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp.F90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp.F90	(revision 4951)
@@ -8,5 +8,5 @@
 SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
      paprs,pplay,temp,qt,ptconv,ratqs,                  &
-     d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri,  &
+     d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
      pfraclr,pfracld,                                   &
      radocond, radicefrac, rain, snow,                  &
@@ -14,8 +14,12 @@
      prfl, psfl, rhcl, qta, fraca,                     &
      tv, pspsk, tla, thl, iflag_cld_th,             &
-     iflag_ice_thermo, ok_ice_sursat, qsatl, qsats,     &
-     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,        &
+     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, qsatl, qsati,        &
+     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,     &
@@ -98,9 +102,9 @@
 USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, icefrac_lscp, calc_gammasat
 USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top
-USE ice_sursat_mod, ONLY : ice_sursat
+USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
 USE lmdz_lscp_poprecip, ONLY : poprecip_precld, poprecip_postcld
 
 ! Use du module lmdz_lscp_ini contenant les constantes
-USE lmdz_lscp_ini, ONLY : prt_level, lunout
+USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
 USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
 USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
@@ -111,4 +115,5 @@
 USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
 USE lmdz_lscp_ini, ONLY : ok_poprecip
+USE lmdz_lscp_ini, ONLY : ok_external_lognormal, ok_ice_supersat, ok_unadjusted_clouds
 
 IMPLICIT NONE
@@ -132,6 +137,4 @@
   INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
                                                                    ! CR: if iflag_ice_thermo=2, only convection is active
-  LOGICAL,                         INTENT(IN)   :: ok_ice_sursat   ! flag to determine if ice sursaturation is activated
-
   LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
 
@@ -150,7 +153,18 @@
   REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs        ! function of pressure that sets the large-scale
 
-
-  ! Input sursaturation en glace
-  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rneb_seri        ! fraction nuageuse en memoire
+  ! INPUT/OUTPUT condensation and ice supersaturation
+  !--------------------------------------------------
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cf_seri          ! cloud fraction [-]
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratio_qi_qtot    ! solid specific water to total specific water ratio [-]
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rvc_seri         ! cloudy water vapor to total water vapor ratio [-]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
+  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: pbl_eps          ! TKE dissipation [?]
+  REAL, DIMENSION(klon),           INTENT(IN)   :: cell_area        ! area of each cell [m2]
+
+  ! INPUT/OUTPUT aviation
+  !--------------------------------------------------
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! Aviation distance flown within the mesh [m/s/mesh]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! Aviation H2O emitted within the mesh [kg H2O/s/mesh]
   
   ! OUTPUT variables
@@ -170,6 +184,4 @@
   REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! surface large-scale rainfall [kg/s/m2] 
   REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! surface large-scale snowfall [kg/s/m2]
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl            ! saturation specific humidity wrt liquid [kg/kg]
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsats            ! saturation specific humidity wrt ice [kg/kg]  
   REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
   REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
@@ -183,22 +195,37 @@
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl        ! scavenging fraction due tu nucleation [-]
   
-  ! for supersaturation and contrails parameterisation
-  
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qclr             ! specific total water content in clear sky region [kg/kg]
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld             ! specific total water content in cloudy region [kg/kg]
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qss              ! specific total water content in supersat region [kg/kg]
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qvc              ! specific vapor content in clouds [kg/kg] 
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rnebclr          ! mesh fraction of clear sky [-]   
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rnebss           ! mesh fraction of ISSR [-]  
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_ss         ! coefficient governing the ice nucleation RHi threshold [-]      
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcontr           ! threshold temperature for contrail formation [K]
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr           ! threshold humidity for contrail formation [kg/kg]
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2          ! // (2nd expression more consistent with LMDZ expression of q)          
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN          ! fraction of grid favourable to non-persistent contrails
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP          ! fraction of grid favourable to persistent contrails
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth      ! mean saturation deficit in thermals
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv     ! mean saturation deficit in environment
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath  ! std of saturation deficit in thermals
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv ! std of saturation deficit in environment
+  ! for condensation and ice supersaturation
+
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsub           !--specific total water content in sub-saturated clear sky region [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qissr          !--specific total water content in supersat region [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld           !--specific total water content in cloudy region [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: subfra         !--mesh fraction of subsaturated clear sky [-]   
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: issrfra        !--mesh fraction of ISSR [-]  
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_cond     !--coefficient governing the ice nucleation RHi threshold [-]      
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_sub        !--cloud fraction tendency because of sublimation [s-1]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_con        !--cloud fraction tendency because of condensation [s-1]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_mix        !--cloud fraction tendency because of cloud mixing [s-1]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_adj        !--specific ice content tendency because of temperature adjustment [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_sub        !--specific ice content tendency because of sublimation [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_con        !--specific ice content tendency because of condensation [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_mix        !--specific ice content tendency because of cloud mixing [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_adj       !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_sub       !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_con       !--specific cloud water vapor tendency because of condensation [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_mix       !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl          !--saturation specific humidity wrt liquid [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsati          !--saturation specific humidity wrt ice [kg/kg]  
+
+  ! for contrails and aviation
+
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcontr         !--threshold temperature for contrail formation [K]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr         !--threshold humidity for contrail formation [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2        !--// (2nd expression more consistent with LMDZ expression of q)
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN        !--fraction of grid favourable to non-persistent contrails
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP        !--fraction of grid favourable to persistent contrails
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_avi        !--cloud fraction tendency because of aviation [s-1]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_avi        !--specific ice content tendency because of aviation [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_avi       !--specific cloud water vapor tendency because of aviation [kg/kg/s]
+
 
   ! for POPRECIP
@@ -217,4 +244,11 @@
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
+
+  ! for thermals
+
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth      !--mean saturation deficit in thermals
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv     !--mean saturation deficit in environment
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath  !--std of saturation deficit in thermals
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv !--std of saturation deficit in environment
 
 
@@ -282,4 +316,15 @@
   REAL :: effective_zneb
   REAL, DIMENSION(klon) :: distcltop1D, temp_cltop1D
+  
+  ! for condensation and ice supersaturation
+  REAL, DIMENSION(klon) :: qvc, shear
+  REAL :: delta_z
+  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrails)
+  ! Constants used for calculating ratios that are advected (using a parent-child
+  ! formalism). This is not done in the dynamical core because at this moment,
+  ! only isotopes can use this parent-child formalism. Note that the two constants
+  ! are the same as the one use in the dynamical core, being also defined in
+  ! dyn3d_common/infotrac.F90
+  REAL :: min_qParent, min_ratio
 
 
@@ -363,15 +408,36 @@
 temp_cltop1D(:) = 0.0
 ztupnew(:)=0.0
-!--ice supersaturation
-gamma_ss(:,:) = 1.0
-qss(:,:) = 0.0
-rnebss(:,:) = 0.0
-Tcontr(:,:) = missing_val
-qcontr(:,:) = missing_val
-qcontr2(:,:) = missing_val
-fcontrN(:,:) = 0.0
-fcontrP(:,:) = 0.0
+
 distcltop(:,:)=0.
 temp_cltop(:,:)=0.
+
+!--Ice supersaturation
+gamma_cond(:,:) = 1.
+qissr(:,:)      = 0.
+issrfra(:,:)    = 0.
+dcf_sub(:,:)    = 0.
+dcf_con(:,:)    = 0.
+dcf_mix(:,:)    = 0.
+dqi_adj(:,:)    = 0.
+dqi_sub(:,:)    = 0.
+dqi_con(:,:)    = 0.
+dqi_mix(:,:)    = 0.
+dqvc_adj(:,:)   = 0.
+dqvc_sub(:,:)   = 0.
+dqvc_con(:,:)   = 0.
+dqvc_mix(:,:)   = 0.
+fcontrN(:,:)    = 0.
+fcontrP(:,:)    = 0.
+Tcontr(:,:)     = missing_val
+qcontr(:,:)     = missing_val
+qcontr2(:,:)    = missing_val
+dcf_avi(:,:)    = 0.
+dqi_avi(:,:)    = 0.
+dqvc_avi(:,:)   = 0.
+qvc(:)          = 0.
+shear(:)        = 0.
+min_qParent     = 1.e-30
+min_ratio       = 1.e-16
+
 !-- poprecip
 qraindiag(:,:)= 0.
@@ -759,4 +825,6 @@
                     rneblsvol(i,k)=ctot_vol(i,k)
                     zqn(i)=qcloud(i)
+                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
+                    qvc(i) = rneb(i,k) * zqs(i)
                 ENDDO
 
@@ -794,12 +862,15 @@
         DO iter=1,iflag_fisrtilp_qsat+1
 
+                keepgoing(:) = .FALSE.
+
                 DO i=1,klon
 
                     ! keepgoing = .true. while convergence is not satisfied
-                    keepgoing(i)=ABS(DT(i)).GT.DDT0
-
-                    IF ((keepgoing(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
+
+                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
 
                         ! if not convergence:
+                        ! we calculate a new iteration
+                        keepgoing(i) = .TRUE.
 
                         ! P2.2.1> cloud fraction and condensed water mass calculation
@@ -833,7 +904,65 @@
                   CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, distcltop1D(:),temp_cltop1D(:),zfice(:),dzfice(:))
 
+
+                  !--AB Activates a condensation scheme that allows for
+                  !--ice supersaturation and contrails evolution from aviation
+                  IF (ok_ice_supersat) THEN
+
+                    !--Calculate the shear value (input for condensation and ice supersat)
+                    DO i = 1, klon
+                      !--Cell thickness [m]
+                      delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
+                      IF ( iftop ) THEN
+                        ! top
+                        shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
+                                       + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
+                      ELSEIF ( k .EQ. 1 ) THEN
+                        ! surface
+                        shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
+                                       + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
+                      ELSE
+                        ! other layers
+                        shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
+                                           - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
+                                       + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
+                                           - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
+                      ENDIF
+                    ENDDO
+
+                    !---------------------------------------------
+                    !--   CONDENSATION AND ICE SUPERSATURATION  --
+                    !---------------------------------------------
+
+                    CALL condensation_ice_supersat( &
+                        klon, dtime, missing_val, &
+                        pplay(:,k), paprs(:,k), paprs(:,k+1), &
+                        cf_seri(:,k), rvc_seri(:,k), ratio_qi_qtot(:,k), &
+                        shear(:), pbl_eps(:,k), cell_area(:), &
+                        Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), keepgoing(:), &
+                        rneb(:,k), zqn(:), qvc(:), issrfra(:,k), qissr(:,k), &
+                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
+                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
+                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
+                        Tcontr(:,k), qcontr(:,k), qcontr2(:,k), fcontrN(:,k), fcontrP(:,k), &
+                        flight_dist(:,k), flight_h2o(:,k), &
+                        dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k))
+
+
+                  !--If .TRUE., calls an externalised version of the generalised
+                  !--lognormal condensation scheme (Bony and Emanuel 2001)
+                  !--Numerically, convergence is conserved with this option
+                  !--The objective is to simplify LSCP
+                  ELSEIF ( ok_external_lognormal ) THEN
+                          
+                   CALL condensation_lognormal( &
+                       klon, Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), &
+                       keepgoing(:), rneb(:,k), zqn(:), qvc(:))
+
+
+                 ELSE !--Generalised lognormal (Bony and Emanuel 2001)
+
                   DO i=1,klon !todoan : check if loop in i is needed
 
-                      IF ((keepgoing(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
+                      IF (keepgoing(i)) THEN
 
                         zpdf_sig(i)=ratqs(i,k)*zq(i)
@@ -849,31 +978,23 @@
                         zpdf_e2(i)=1.-erf(zpdf_e2(i))
 
-                        IF ((.NOT.ok_ice_sursat).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)
+                              zqn(i)=zqs(i)
+                              !--AB grid-mean vapor in the cloud - we assume saturation adjustment
+                              qvc(i) = 0.
                           ELSE
                               rneb(i,k)=0.5*zpdf_e1(i)
                               zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
+                              !--AB grid-mean vapor in the cloud - we assume saturation adjustment
+                              qvc(i) = rneb(i,k) * zqs(i)
                           ENDIF
 
-                          rnebss(i,k)=0.0   !--added by OB (needed because of the convergence loop on time)
-                          fcontrN(i,k)=0.0  !--idem
-                          fcontrP(i,k)=0.0  !--idem
-                          qss(i,k)=0.0      !--idem
-
-                       ELSE ! in case of ice supersaturation by Audran
-
-                        !------------------------------------
-                        ! ICE SUPERSATURATION
-                        !------------------------------------
-
-                        CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, temp(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))
+                      ENDIF ! keepgoing
+                  ENDDO ! loop on klon
+
+                  ENDIF ! .NOT. ok_ice_supersat
+
+                  DO i=1,klon
+                      IF (keepgoing(i)) THEN
 
                         ! If vertical heterogeneity, change fraction by volume as well
@@ -958,8 +1079,20 @@
                     zqn(i) = zq(i)
                     rneb(i,k) = 1.0
-                    zcond(i) = MAX(0.0,zqn(i)-zqs(i))
+                    IF ( ok_unadjusted_clouds ) THEN
+                      !--AB We relax the saturation adjustment assumption
+                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
+                      zcond(i) = MAX(0., zqn(i) - qvc(i))
+                    ELSE
+                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
+                    ENDIF
                     rhcl(i,k)=1.0
                 ELSE
-                    zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
+                    IF ( ok_unadjusted_clouds ) THEN
+                      !--AB We relax the saturation adjustment assumption
+                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
+                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
+                    ELSE
+                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
+                    ENDIF
                     ! following line is very strange and probably wrong:
                     rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
@@ -985,4 +1118,47 @@
           rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
         ENDIF
+    
+    !--AB Write diagnostics and tracers for ice supersaturation
+    IF ( ok_ice_supersat ) THEN
+      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
+      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
+
+      DO i = 1, klon
+
+        cf_seri(i,k) = rneb(i,k)
+
+        IF ( .NOT. ok_unadjusted_clouds ) THEN
+          qvc(i) = zqs(i) * rneb(i,k)
+        ENDIF
+        IF ( zq(i) .GT. min_qParent ) THEN
+          rvc_seri(i,k) = qvc(i) / zq(i)
+        ELSE
+          rvc_seri(i,k) = min_ratio
+        ENDIF
+        !--The MIN barrier is NEEDED because of:
+        !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
+        !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
+        !--The MAX barrier is a safeguard that should not be activated
+        rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
+
+        !--Diagnostics
+        gamma_cond(i,k) = gammasat(i)
+        IF ( issrfra(i,k) .LT. eps ) THEN
+          issrfra(i,k) = 0.
+          qissr(i,k) = 0.
+        ENDIF
+        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
+        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
+        IF ( subfra(i,k) .LT. eps ) THEN
+          subfra(i,k) = 0.
+          qsub(i,k) = 0.
+        ENDIF
+        qcld(i,k) = qvc(i) + zcond(i)
+        IF ( cf_seri(i,k) .LT. eps ) THEN
+          qcld(i,k) = 0.
+        ENDIF
+      ENDDO
+    ENDIF
+
 
     ! ----------------------------------------------------------------
@@ -1407,23 +1583,4 @@
     ENDDO
 
-    !--save some variables for ice supersaturation
-    !
-    DO i = 1, klon
-        ! for memory 
-        rneb_seri(i,k) = rneb(i,k)
-
-        ! for 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))  !--added by OB because of pathologic cases with the lognormal
-        qcld(i,k) = qvc(i,k) + zcond(i)
-   ENDDO
-   !q_sat
-   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs(:))
-   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,qsats(:,k),zdqs(:))
-
-
-
     ! Outputs:
     !-------------------------------
Index: /LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_condensation.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_condensation.F90	(revision 4951)
+++ /LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_condensation.F90	(revision 4951)
@@ -0,0 +1,1009 @@
+MODULE lmdz_lscp_condensation
+!----------------------------------------------------------------
+! Module for condensation of clouds routines
+! that are called in LSCP
+
+
+IMPLICIT NONE
+
+CONTAINS
+
+!**********************************************************************************
+SUBROUTINE condensation_lognormal( &
+      klon, temp, qtot, qsat, gamma_cond, ratqs, &
+      keepgoing, cldfra, qincld, qvc)
+
+!----------------------------------------------------------------------
+! This subroutine calculates the formation of clouds, using a
+! statistical cloud scheme. It uses a generalised lognormal distribution
+! of total water in the gridbox
+! See Bony and Emanuel, 2001
+!----------------------------------------------------------------------
+
+USE lmdz_lscp_ini, ONLY: eps
+
+IMPLICIT NONE
+
+!
+! Input
+!
+INTEGER,  INTENT(IN)                     :: klon          ! number of horizontal grid points
+!
+REAL,     INTENT(IN)   , DIMENSION(klon) :: temp          ! temperature [K]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: qtot          ! total specific humidity (without precip) [kg/kg]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: qsat          ! saturation specific humidity [kg/kg]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: gamma_cond    ! condensation threshold w.r.t. qsat [-]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: ratqs         ! ratio between the variance of the total water distribution and its average [-]
+LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: keepgoing     ! .TRUE. if a new condensation loop should be computed
+!
+!  Output
+!  NB. those are in INOUT because of the convergence loop on temperature
+!      (in some cases, the values are not re-computed) but the values
+!      are never used explicitely
+!
+REAL,     INTENT(INOUT), DIMENSION(klon) :: cldfra   ! cloud fraction [-]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: qincld   ! cloud-mean in-cloud total specific water [kg/kg]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: qvc      ! gridbox-mean vapor in the cloud [kg/kg]
+!
+! Local
+!
+INTEGER :: i
+REAL :: pdf_std, pdf_k, pdf_delta
+REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
+!
+!--Loop on klon
+!
+DO i = 1, klon
+
+  !--If a new calculation of the condensation is needed,
+  !--i.e., temperature has not yet converged (or the cloud is
+  !--formed elsewhere)
+  IF (keepgoing(i)) THEN
+
+    pdf_std   = ratqs(i) * qtot(i)
+    pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2. ) )
+    pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
+    pdf_a     = pdf_delta / ( pdf_k * SQRT(2.) )
+    pdf_b     = pdf_k / (2. * SQRT(2.))
+    pdf_e1    = pdf_a - pdf_b
+    pdf_e1    = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
+    pdf_e1    = 1. - ERF(pdf_e1)
+    pdf_e2    = pdf_a + pdf_b
+    pdf_e2    = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
+    pdf_e2    = 1. - ERF(pdf_e2)
+
+
+    IF ( pdf_e1 .LT. eps ) THEN
+      cldfra(i) = 0.
+      qincld(i) = qsat(i)
+      !--AB grid-mean vapor in the cloud - we assume saturation adjustment
+      qvc(i) = 0.
+    ELSE
+      cldfra(i) = 0.5 * pdf_e1
+      qincld(i) = qtot(i) * pdf_e2 / pdf_e1
+      !--AB grid-mean vapor in the cloud - we assume saturation adjustment
+      qvc(i) = qsat(i) * cldfra(i)
+    ENDIF
+
+  ENDIF   ! end keepgoing
+
+ENDDO     ! end loop on i
+END SUBROUTINE condensation_lognormal
+!**********************************************************************************
+
+!**********************************************************************************
+SUBROUTINE condensation_ice_supersat( &
+      klon, dtime, missing_val, pplay, paprsdn, paprsup, &
+      cf_seri, rvc_seri, ratio_qi_qtot, shear, pbl_eps, cell_area, &
+      temp, qtot, qsat, gamma_cond, ratqs, keepgoing, &
+      cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, &
+      dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, &
+      Tcontr, qcontr, qcontr2, fcontrN, fcontrP, flight_dist, flight_h2o, &
+      dcf_avi, dqi_avi, dqvc_avi)
+
+!----------------------------------------------------------------------
+! This subroutine calculates the formation, evolution and dissipation
+! of clouds, using a process-oriented treatment of the cloud properties
+! (cloud fraction, vapor in the cloud, condensed water in the cloud).
+! It allows for ice supersaturation in cold regions, in clear sky.
+! If ok_unadjusted_clouds, it also allows for sub- and supersaturated
+! cloud water vapors.
+! It also allows for the formation and evolution of condensation trails
+! (contrails) from aviation.
+! Authors: Audran Borella, Etienne Vignon, Olivier Boucher
+! April 2024
+!----------------------------------------------------------------------
+
+USE lmdz_lscp_tools,      ONLY: calc_qsat_ecmwf, calc_gammasat, GAMMAINC
+USE lmdz_lscp_ini,        ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI, EPS_W
+USE lmdz_lscp_ini,        ONLY: eps, temp_nowater, ok_weibull_warm_clouds
+USE lmdz_lscp_ini,        ONLY: ok_unadjusted_clouds, iflag_cloud_sublim_pdf
+USE lmdz_lscp_ini,        ONLY: lunout
+
+USE lmdz_lscp_ini,        ONLY: mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, &
+                                rhlmid_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, rhl0_pdf_lscp, &
+                                cond_thresh_pdf_lscp, coef_mixing_lscp, coef_shear_lscp, & 
+                                chi_mixing_lscp, rho_ice
+
+IMPLICIT NONE
+
+!
+! Input
+!
+INTEGER,  INTENT(IN)                     :: klon          ! number of horizontal grid points
+REAL,     INTENT(IN)                     :: dtime         ! time step [s]
+REAL,     INTENT(IN)                     :: missing_val   ! missing value for output
+!
+REAL,     INTENT(IN)   , DIMENSION(klon) :: pplay         ! layer pressure [Pa]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsdn       ! pressure at the lower interface [Pa]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsup       ! pressure at the upper interface [Pa]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: cf_seri       ! cloud fraction [-]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: rvc_seri      ! gridbox-mean water vapor in cloud [kg/kg]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: ratio_qi_qtot ! specific ice water content to total specific water ratio [-]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: shear         ! vertical shear [s-1]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: pbl_eps       ! 
+REAL,     INTENT(IN)   , DIMENSION(klon) :: cell_area     ! 
+REAL,     INTENT(IN)   , DIMENSION(klon) :: temp          ! temperature [K]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: qtot          ! total specific humidity (without precip) [kg/kg]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: qsat          ! saturation specific humidity [kg/kg]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: gamma_cond    ! condensation threshold w.r.t. qsat [-]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: ratqs         ! ratio between the variance of the total water distribution and its average [-]
+LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: keepgoing     ! .TRUE. if a new condensation loop should be computed
+!
+! Input for aviation
+!
+REAL,     INTENT(IN),    DIMENSION(klon) :: flight_dist   ! 
+REAL,     INTENT(IN),    DIMENSION(klon) :: flight_h2o    ! 
+!
+!  Output
+!  NB. cldfra and qincld should be outputed as cf_seri and qi_seri,
+!      or as tendencies (maybe in the future)
+!  NB. those are in INOUT because of the convergence loop on temperature
+!      (in some cases, the values are not re-computed) but the values
+!      are never used explicitely
+!
+REAL,     INTENT(INOUT), DIMENSION(klon) :: cldfra   ! cloud fraction [-]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: qincld   ! cloud-mean in-cloud total specific water [kg/kg]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: qvc      ! gridbox-mean vapor in the cloud [kg/kg]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: issrfra  ! ISSR fraction [-]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: qissr    ! gridbox-mean ISSR specific water [kg/kg]
+!
+!  Diagnostics for condensation and ice supersaturation
+!  NB. idem for the INOUT
+!
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_sub  ! cloud fraction tendency because of sublimation [s-1]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_con  ! cloud fraction tendency because of condensation [s-1]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_mix  ! cloud fraction tendency because of cloud mixing [s-1]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_adj  ! specific ice content tendency because of temperature adjustment [kg/kg/s]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_sub  ! specific ice content tendency because of sublimation [kg/kg/s]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_con  ! specific ice content tendency because of condensation [kg/kg/s]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_mix  ! specific ice content tendency because of cloud mixing [kg/kg/s]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_adj ! specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_sub ! specific cloud water vapor tendency because of sublimation [kg/kg/s]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_con ! specific cloud water vapor tendency because of condensation [kg/kg/s]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
+!
+!  Diagnostics for aviation
+!  NB. idem for the INOUT
+!
+REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcontr   ! critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) [K]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: qcontr   ! 
+REAL,     INTENT(INOUT), DIMENSION(klon) :: qcontr2  ! 
+REAL,     INTENT(INOUT), DIMENSION(klon) :: fcontrN  ! 
+REAL,     INTENT(INOUT), DIMENSION(klon) :: fcontrP  ! 
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_avi  ! cloud fraction tendency because of aviation [s-1]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_avi  ! specific ice content tendency because of aviation [kg/kg/s]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_avi ! specific cloud water vapor tendency because of aviation [kg/kg/s]
+!
+! Local
+!
+INTEGER :: i
+LOGICAL :: ok_warm_cloud
+REAL, DIMENSION(klon) :: qcld, qzero
+!
+! for lognormal
+REAL :: pdf_std, pdf_k, pdf_delta
+REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
+!
+! for unadjusted clouds
+REAL :: qvapincld, qvapincld_new
+REAL :: qiceincld, qice_ratio
+REAL :: pres_sat, kappa
+REAL :: air_thermal_conduct, water_vapor_diff
+REAL :: r_ice
+!
+! for sublimation
+REAL :: pdf_alpha
+REAL :: dqt_sub
+!
+! for condensation
+REAL, DIMENSION(klon) :: qsatl, dqsatl
+REAL :: clrfra, qclr, sl_clr, rhl_clr
+REAL :: pdf_ratqs, pdf_skew, pdf_scale, pdf_loc
+REAL :: pdf_x, pdf_y, pdf_T
+REAL :: pdf_e3, pdf_e4
+REAL :: cf_cond, qt_cond, dqt_con
+!
+! for mixing
+REAL, DIMENSION(klon) :: subfra, qsub
+REAL :: dqt_mix_sub, dqt_mix_issr
+REAL :: dcf_mix_sub, dcf_mix_issr
+REAL :: dqvc_mix_sub, dqvc_mix_issr
+REAL :: dqt_mix
+REAL :: a_mix, bovera, N_cld_mix, L_mix
+REAL :: envfra_mix, cldfra_mix
+REAL :: L_shear, shear_fra
+REAL :: sigma_mix, issrfra_mix, subfra_mix, qvapinmix
+!
+! for cell properties
+REAL :: rho, rhodz, dz
+!REAL :: V_cell, M_cell
+!
+! for aviation and cell properties
+!REAL :: dqt_avi
+!REAL :: contrail_fra
+!
+!
+!--more local variables for diagnostics
+!--imported from YOMCST.h 
+!--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1)
+!--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air
+!--values from Schumann, Meteorol Zeitschrift, 1996
+!--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen
+!--Qheat = 43.  /  50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen
+!REAL, PARAMETER :: EiH2O=1.25  !--emission index of water vapour for kerosene (kg kg-1)
+!REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1) 
+!REAL, PARAMETER :: eta=0.3     !--average propulsion efficiency of the aircraft
+!--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute 
+!--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010.
+!REAL :: Gcontr
+!--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) 
+!--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1)
+!REAL :: qsatliqcontr
+
+
+!-----------------------------------------------
+! Initialisations
+!-----------------------------------------------
+
+! Ajout des émissions de H2O dues à l'aviation
+! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
+! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O )  
+!      = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) ) 
+! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q) 
+! flight_h2O is in kg H2O / s / cell
+!
+!IF (ok_plane_h2o) THEN 
+!  q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) ) 
+!ENDIF
+
+
+qzero(:) = 0.
+
+!--Calculation of qsat w.r.t. liquid
+CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 1, .FALSE., qsatl, dqsatl)
+
+!
+!--Loop on klon
+!
+DO i = 1, klon
+
+  !--If a new calculation of the condensation is needed,
+  !--i.e., temperature has not yet converged (or the cloud is
+  !--formed elsewhere)
+  IF (keepgoing(i)) THEN
+
+    !--If the temperature is higher than the threshold below which
+    !--there is no liquid in the gridbox, we activate the usual scheme
+    !--(generalised lognormal from Bony and Emanuel 2001)
+    !--If ok_weibull_warm_clouds = .TRUE., the Weibull law is used for
+    !--all clouds, and the lognormal scheme is not activated
+    IF ( ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds ) THEN
+
+      pdf_std   = ratqs(i) * qtot(i)
+      pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2. ) )
+      pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
+      pdf_a     = pdf_delta / ( pdf_k * SQRT(2.) )
+      pdf_b     = pdf_k / (2. * SQRT(2.))
+      pdf_e1    = pdf_a - pdf_b
+      pdf_e1    = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
+      pdf_e1    = 1. - ERF(pdf_e1)
+      pdf_e2    = pdf_a + pdf_b
+      pdf_e2    = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
+      pdf_e2    = 1. - ERF(pdf_e2)
+
+
+      IF ( pdf_e1 .LT. eps ) THEN
+        cldfra(i) = 0.
+        qincld(i) = qsat(i)
+        qvc(i) = 0.
+      ELSE
+        cldfra(i) = 0.5 * pdf_e1
+        qincld(i) = qtot(i) * pdf_e2 / pdf_e1
+        qvc(i) = qsat(i) * cldfra(i)
+      ENDIF
+
+    !--If the temperature is lower than temp_nowater, we use the new
+    !--condensation scheme that allows for ice supersaturation
+    ELSE
+
+      !--Initialisation
+      IF ( temp(i) .GT. temp_nowater ) THEN
+        !--If the air mass is warm (liquid water can exist),
+        !--all the memory is lost and the scheme becomes statistical,
+        !--i.e., the sublimation and mixing processes are deactivated,
+        !--and the condensation process is slightly adapted
+        !--This can happen only if ok_weibull_warm_clouds = .TRUE.
+        cldfra(i) = 0.
+        qvc(i)    = 0.
+        qcld(i)   = 0.
+        ok_warm_cloud = .TRUE.
+      ELSE
+        !--The following barriers ensure that the traced cloud properties
+        !--are consistent. In some rare cases, i.e. the cloud water vapor
+        !--can be greater than the total water in the gridbox
+        cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
+        qcld(i)   = MAX(0., MIN(qtot(i), ( ratio_qi_qtot(i) + rvc_seri(i) ) * qtot(i))) 
+        qvc(i)    = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
+        ok_warm_cloud = .FALSE.
+      ENDIF
+
+      dcf_sub(i)  = 0.
+      dqi_sub(i)  = 0.
+      dqvc_sub(i) = 0.
+      dqi_adj(i)  = 0.
+      dqvc_adj(i) = 0.
+      dcf_con(i)  = 0.
+      dqi_con(i)  = 0.
+      dqvc_con(i) = 0.
+      dcf_mix(i)  = 0.
+      dqi_mix(i)  = 0.
+      dqvc_mix(i) = 0.
+
+      issrfra(i)  = 0.
+      qissr(i)    = 0.
+      subfra(i)   = 0.
+      qsub(i)     = 0.
+
+      !--Initialisation of the cell properties
+      !--Dry density [kg/m3]
+      rho = pplay(i) / temp(i) / RD
+      !--Dry air mass [kg/m2]
+      rhodz = ( paprsdn(i) - paprsup(i) ) / RG
+      !--Cell thickness [m]
+      dz = rhodz / rho
+      !--Cell volume [m3]
+      !V_cell = dz * cell_area(i)
+      !--Cell dry air mass [kg]
+      !M_cell = rhodz * cell_area(i)
+
+
+      IF ( ok_unadjusted_clouds ) THEN
+        !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
+        !--hypothesis is lost, and the vapor in the cloud is purely prognostic.
+        !
+        !--The deposition equation is
+        !-- dmi/dt = alpha*4pi*C*Svi / ( R*T/esi/Mw/Dv + Ls/ka/T * (Ls*Mw/R/T - 1) )
+        !--from Pruppacher and Klett (2010), where
+        !--mi is the mass of one ice crystal [kg]
+        !--C is the capacitance of an ice crystal [m]
+        !--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
+        !--R is the perfect gas constant [J/mol/K]
+        !--T is the temperature [K]
+        !--esi is the saturation pressure w.r.t. ice [Pa]
+        !--Mw is the molar mass of water [kg/mol]
+        !--Dv is the diffusivity of water vapor [m2/s]
+        !--Ls is the specific latent heat of sublimation [J/kg/K]
+        !--ka is the thermal conductivity of dry air [J/m/s/K]
+        !
+        !--We have qice = Nice * mi, where Nice is the ice crystal
+        !--number concentration per kg of moist air
+        !--HYPOTHESIS 1: the ice crystals are spherical, therefore
+        !-- C = r_ice
+        !-- mi = 4/3 * pi * r_ice**3 * rho_ice
+        !--HYPOTHESIS 2: the ice crystals are monodisperse with the
+        !--initial radius r_ice_0.
+        !--NB. this is notably different than the assumption
+        !--of a distributed qice in the cloud made in the sublimation process;
+        !--should it be consistent?
+        !
+        !--As the deposition process does not create new ice crystals,
+        !--and because we assume a same r_ice value for all crystals
+        !--therefore the sublimation process does not destroy ice crystals
+        !--(or, in a limit case, it destroys all ice crystals), then
+        !--Nice is a constant during the sublimation/deposition process.
+        !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI r_ice_0**3 rho_ice )
+        !
+        !--The deposition equation then reads:
+        !-- dqi/dt = alpha*4pi*r_ice*(qvc-qsat)/qsat / ( R*T/esi/Mw/Dv + Ls/ka/T * (Ls*Mw/R/T - 1) ) * Nice
+        !-- dqi/dt = alpha*4pi* (qi / qi_0)**(1/3) *r_ice_0*(qvc-qsat)/qsat &
+        !--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
+        !--                         * qi_0 / ( 4/3 RPI r_ice_0**3 rho_ice )
+        !-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
+        !--  *alpha/qsat/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 r_ice_0**2 rho_ice )
+        !--and we have
+        !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * qi_0**(2/3) / r_ice_0**2
+        !-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * qi_0**(2/3) / r_ice_0**2
+        !--where kappa = (2/3*rho_ice)*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
+        !
+        !--This system of equations can be resolved with an exact
+        !--explicit numerical integration, having one variable resolved
+        !--explicitly, the other exactly. The exactly resolved variable is
+        !--the one decreasing, so it is qvc if the process is
+        !--condensation, qi if it is sublimation.
+        !
+        !--kappa is computed as an initialisation constant, as it depends only
+        !--on temperature and other pre-computed values
+        pres_sat = qsat(i) / ( EPS_W + ( 1. - EPS_W ) * qsat(i) ) * pplay(i)
+        !--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
+        air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
+        !--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
+        water_vapor_diff = 0.211 * ( temp(i) / RTT )**1.94 * ( 101325. / pplay(i) ) * 1.e-4
+        kappa = 2. / 3. * rho_ice * qsat(i) &
+              * ( RV * temp(i) / water_vapor_diff / pres_sat &
+                + RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) )
+      ENDIF
+
+
+      !-------------------------------------------------------------------
+      !--    SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CLOUD    --
+      !-------------------------------------------------------------------
+      
+      !--If there is a cloud
+      IF ( cldfra(i) .GT. eps ) THEN
+        
+        qvapincld = qvc(i) / cldfra(i)
+        !--The vapor in cloud cannot be higher than the
+        !--condensation threshold
+        qvapincld = MIN(qvapincld, gamma_cond(i) * qsat(i))
+        qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
+        
+        !--If the ice water content is too low, the cloud is purely sublimated
+        IF ( qiceincld .LT. eps ) THEN
+          dcf_sub(i) = - cldfra(i)
+          dqvc_sub(i) = - qvc(i)
+          dqi_sub(i) = - ( qcld(i) - qvc(i) )
+
+          cldfra(i) = 0.
+          qcld(i) = 0.
+          qvc(i) = 0.
+
+        !--Else, the cloud is adjusted and sublimated
+        ELSE
+
+          IF ( ok_unadjusted_clouds ) THEN
+            !--Here, the initial vapor in the cloud is qvapincld, and we compute
+            !--the new vapor qvapincld_new
+
+            !--r_ice formula from Sun and Rikus (1999)
+            r_ice = 1.e-6 * ( 45.8966 * ( qiceincld * rho * 1e3 )**0.2214 &
+                  + 0.7957 * ( qiceincld * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2.
+
+            IF ( qvapincld .GE. qsat(i) ) THEN
+              !--If the cloud is initially supersaturated
+              !--Exact explicit integration (qvc exact, qice explicit)
+              qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) &
+                            * EXP( - dtime * qiceincld / kappa / r_ice**2. )
+            ELSE
+              !--If the cloud is initially subsaturated
+              !--Exact explicit integration (qice exact, qvc explicit)
+              !--The barrier is set so that the resulting vapor in cloud
+              !--cannot be greater than qsat
+              !--qice_ratio is the ratio between the new ice content and
+              !--the old one, it is comprised between 0 and 1
+              qice_ratio = ( 1. - 2. / 3. / kappa / r_ice**2. * dtime * ( qsat(i) - qvapincld ) )
+
+              IF ( qice_ratio .LT. 0. ) THEN
+                !--If all the ice has been sublimated, we sublimate
+                !--completely the cloud and do not activate the sublimation
+                !--process
+                !--Tendencies and diagnostics
+                dcf_sub(i) = - cldfra(i)
+                dqvc_sub(i) = - qvc(i)
+                dqi_sub(i) = - ( qcld(i) - qvc(i) )
+
+                cldfra(i) = 0.
+                qcld(i) = 0.
+                qvc(i) = 0.
+
+                !--The new vapor in cloud is set to 0 so that the
+                !--sublimation process does not activate
+                qvapincld_new = 0.
+              ELSE
+                !--Else, the sublimation process is activated with the
+                !--diagnosed new cloud water vapor
+                !--The new vapor in the cloud is increased with the
+                !--sublimated ice
+                qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**(3./2.) )
+                !--The new vapor in the cloud cannot be greater than qsat
+                qvapincld_new = MIN(qvapincld_new, qsat(i))
+              ENDIF ! qice_ratio .LT. 0.
+            ENDIF ! qvapincld .GT. qsat(i)
+          ELSE
+            !--We keep the saturation adjustment hypothesis, and the vapor in the
+            !--cloud is set equal to the saturation vapor
+            qvapincld_new = qsat(i)
+          ENDIF ! ok_unadjusted_clouds
+
+          !--If the vapor in cloud is below vapor needed for the cloud to survive
+          IF ( qvapincld .LT. qvapincld_new ) THEN
+            !--Sublimation of the subsaturated cloud
+            !--iflag_cloud_sublim_pdf selects the PDF of the ice water content
+            !--to use.
+            !--iflag = 1 --> uniform distribution
+            !--iflag = 2 --> exponential distribution
+            !--iflag = 3 --> gamma distribution (Karcher et al 2018)
+
+            IF ( iflag_cloud_sublim_pdf .EQ. 1 ) THEN
+              !--Uniform distribution starting at qvapincld
+              pdf_e1 = 1. / ( 2. * qiceincld )
+
+              dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1
+              dqt_sub    = - cldfra(i) * ( qvapincld_new**2. - qvapincld**2. ) &
+                                       * pdf_e1 / 2.
+
+            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 2 ) THEN
+              !--Exponential distribution starting at qvapincld
+              pdf_alpha = 1. / qiceincld
+              pdf_e1 = EXP( - pdf_alpha * ( qvapincld_new - qvapincld ) )
+
+              dcf_sub(i) = - cldfra(i) * ( 1. - pdf_e1 )
+              dqt_sub    = - cldfra(i) * ( ( 1. - pdf_e1 ) / pdf_alpha &
+                                       + qvapincld - qvapincld_new * pdf_e1 )
+
+            ELSEIF ( iflag_cloud_sublim_pdf .GE. 3 ) THEN
+              !--Gamma distribution starting at qvapincld
+              pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld
+              pdf_y = pdf_alpha * ( qvapincld_new - qvapincld )
+              pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y )
+              pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y )
+
+              dcf_sub(i) = - cldfra(i) * pdf_e1
+              dqt_sub    = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 )
+            ENDIF
+
+            !--Tendencies and diagnostics
+            dqvc_sub(i) = dcf_sub(i) * qvapincld
+            dqi_sub(i) = dqt_sub - dqvc_sub(i)
+
+            !--Add tendencies
+            cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i))
+            qcld(i)   = MAX(0., qcld(i) + dqt_sub)
+            qvc(i)    = MAX(0., qvc(i) + dqvc_sub(i))
+
+          ENDIF ! qvapincld .LT. qvapincld_new
+
+          !--Adjustment of the IWC of the remaining cloud
+          !--to the new vapor in cloud (this can be either positive or negative)
+          dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )
+          dqi_adj(i) = - dqvc_adj(i)
+
+          !--Add tendencies
+          !--The vapor in the cloud is updated, but not qcld as it is constant
+          !--through this process, as well as cldfra which is unmodified
+          qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))
+
+        ENDIF   ! qiceincld .LT. eps
+      ENDIF ! cldfra(i) .GT. eps
+
+
+      !--------------------------------------------------------------------------
+      !--    CONDENSATION AND DIAGNOTICS OF SUB- AND SUPERSATURATED REGIONS    --
+      !--------------------------------------------------------------------------
+      !--This section relies on a distribution of water in the clear-sky region of
+      !--the mesh.
+
+      !--If there is a clear-sky region
+      IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN
+
+        !--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
+        qclr = qtot(i) - qcld(i)
+
+        !--New PDF
+        rhl_clr = qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100.
+        rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp)
+
+        !--Calculation of the properties of the PDF
+        !--Parameterization from IAGOS observations
+        !--pdf_e1 and pdf_e2 will be reused below
+
+        pdf_e1 = beta_pdf_lscp * cell_area(i)**( 1. / 4. ) &
+               * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
+        IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN
+          pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp
+        ELSE
+          pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp
+        ENDIF
+        pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
+        pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3
+
+        pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
+        pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. ))
+        pdf_loc = rhl_clr - pdf_scale * pdf_e2
+
+        !--Diagnostics of ratqs
+        ratqs(i) = pdf_std / ( qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100. )
+
+        !--Calculation of the newly condensed water and fraction (pronostic)
+        !--Integration of the clear sky PDF between gamma_cond*qsat and +inf
+        !--NB. the calculated values are clear-sky averaged
+
+        pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
+        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
+        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
+        pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
+        cf_cond = EXP( - pdf_y )
+        qt_cond = ( pdf_e3 + pdf_loc * cf_cond ) * qsatl(i) / 100.
+
+        IF ( ok_warm_cloud ) THEN
+          !--If the statistical scheme is activated, the first
+          !--calculation of the tendencies is kept for the diagnosis of
+          !--the cloud properties, and the condensation process stops
+          cldfra(i) = cf_cond
+          qcld(i) = qt_cond
+          qvc(i) = cldfra(i) * qsat(i)
+
+        ELSEIF ( cf_cond .GT. cond_thresh_pdf_lscp ) THEN
+          !--We check that there is something to condense, if not the
+          !--condensation process stops
+
+          !--Calculation of the limit qclr in cloud value (stored in pdf_e4)
+
+          pdf_e3 = ( LOG( 1. / cond_thresh_pdf_lscp ) ** ( 1. / pdf_alpha ) - pdf_e2 ) &
+                 / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. )
+          pdf_e4 = ( 2. * pdf_e3 * pdf_e1 - pdf_x ) &
+                 / ( pdf_e3 * pdf_e1 / rhlmid_pdf_lscp - 1. )
+          IF ( ( MIN(pdf_e4, rhl_clr) .LT. rhlmid_pdf_lscp ) .OR. ( pdf_e4 .GT. rhl_clr ) ) THEN
+            pdf_e4 = pdf_x / ( pdf_e3 * pdf_e1 / rhlmid_pdf_lscp + 1. )
+          ENDIF
+          pdf_e4 = pdf_e4 * qsatl(i) / 100.
+
+          !--Calculation of the final tendencies
+          dcf_con(i) = ( 1. - cldfra(i) ) * ( pdf_e4 - qclr / ( 1. - cldfra(i) ) ) &
+                     / ( pdf_e4 - qt_cond / cf_cond )
+          dqt_con = qclr - pdf_e4 * ( 1. - cldfra(i) - dcf_con(i) )
+
+          
+          !--Barriers
+          dcf_con(i) = MIN(dcf_con(i), 1. - cldfra(i))
+          dqt_con = MIN(dqt_con, qclr)
+
+
+          IF ( ok_unadjusted_clouds ) THEN
+            !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute
+            !--the new vapor qvapincld. The timestep is divided by two because we do not
+            !--know when the condensation occurs
+            qvapincld = gamma_cond(i) * qsat(i)
+            qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
+            !--r_ice formula from Sun and Rikus (1999)
+            r_ice = 1.e-6 * ( 45.8966 * ( qiceincld * rho * 1e3 )**0.2214 &
+                  + 0.7957 * ( qiceincld * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2.
+            !--As qvapincld is necessarily greater than qsat, we only
+            !--use the exact explicit formulation
+            !--Exact explicit version
+            qvapincld = qsat(i) + ( qvapincld - qsat(i) ) &
+                      * EXP( - dtime / 2. * qiceincld / kappa / r_ice**2. )
+          ELSE
+            !--We keep the saturation adjustment hypothesis, and the vapor in the
+            !--newly formed cloud is set equal to the saturation vapor.
+            qvapincld = qsat(i)
+          ENDIF
+
+          !--Tendency on cloud vapor and diagnostic
+          dqvc_con(i) = qvapincld * dcf_con(i)
+          dqi_con(i) = dqt_con - dqvc_con(i)
+
+          !--Note that the tendencies are NOT added because they are
+          !--added after the mixing process. In the following, the gridbox fraction is
+          !-- 1. - dcf_con(i), and the total water in the gridbox is
+          !-- qtot(i) - dqi_con(i) - dqvc_con(i)
+
+        ENDIF ! ok_warm_cloud, cf_cond .GT. cond_thresh_pdf_lscp
+      ENDIF ! ( 1. - cldfra(i) ) .GT. eps
+
+      !--If there is still clear sky, we diagnose the ISSR
+      !--We recalculte the PDF properties (after the condensation process)
+      IF ( ( ( 1. - dcf_con(i) - cldfra(i) ) .GT. eps ) .AND. .NOT. ok_warm_cloud ) THEN
+        !--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
+        qclr = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i)
+
+        !--New PDF
+        rhl_clr = qclr / ( 1. - dcf_con(i) - cldfra(i) ) / qsatl(i) * 100.
+        rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp)
+
+        !--Calculation of the properties of the PDF
+        !--Parameterization from IAGOS observations
+        !--pdf_e1 and pdf_e2 will be reused below
+
+        pdf_e1 = beta_pdf_lscp * cell_area(i)**( 1. / 4. ) &
+               * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
+        IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN
+          pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp
+        ELSE
+          pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp
+        ENDIF
+        pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
+        pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3
+
+        pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
+        pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. ))
+        pdf_loc = rhl_clr - pdf_scale * pdf_e2
+
+        !--We then calculate the part that is greater than qsat
+        !--and consider it supersaturated
+
+        pdf_x = qsat(i) / qsatl(i) * 100.
+        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
+        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
+        pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
+        issrfra(i) = EXP( - pdf_y ) * ( 1. - dcf_con(i) - cldfra(i) )
+        qissr(i) = ( pdf_e3 * ( 1. - dcf_con(i) - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
+      ENDIF
+
+      !--Calculation of the subsaturated clear sky fraction and water
+      subfra(i) = 1. - dcf_con(i) - cldfra(i) - issrfra(i)
+      qsub(i) = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i) - qissr(i)
+
+
+      !--------------------------------------
+      !--           CLOUD MIXING           --
+      !--------------------------------------
+      !--This process mixes the cloud with its surroundings: the subsaturated clear sky,
+      !--and the supersaturated clear sky. It is activated if the cloud is big enough,
+      !--but does not cover the entire mesh.
+      !
+      IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) &
+           .AND. .NOT. ok_warm_cloud ) THEN
+
+        !--Initialisation
+        dcf_mix_sub   = 0.
+        dqt_mix_sub   = 0.
+        dqvc_mix_sub  = 0.
+        dcf_mix_issr  = 0.
+        dqt_mix_issr  = 0.
+        dqvc_mix_issr = 0.
+
+
+        !-- PART 1 - TURBULENT DIFFUSION
+
+        !--Clouds within the mesh are assumed to be ellipses. The length of the
+        !--semi-major axis is a and the length of the semi-minor axis is b.
+        !--N_cld_mix is the number of clouds within the mesh, and
+        !--clouds_perim is the total perimeter of the clouds within the mesh,
+        !--not considering interfaces with other meshes (only the interfaces with clear
+        !--sky are taken into account).
+        !--
+        !--The area of each cloud is A = a * b * RPI,
+        !--and the perimeter of each cloud is
+        !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
+        !--
+        !--With cell_area the area of the cell, we have:
+        !-- cldfra = A * N_cld_mix / cell_area
+        !-- clouds_perim = P * N_cld_mix
+        !--
+        !--We assume that the ratio between b and a is a function of
+        !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
+        !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
+        !--are spherical.
+        !-- b / a = bovera = MAX(0.1, cldfra)
+        bovera = MAX(0.1, cldfra(i))
+        !--The clouds perimeter is imposed using the formula from Morcrette 2012,
+        !--based on observations.
+        !-- clouds_perim_normalized = alpha * cldfra * ( 1. - cldfra ) = clouds_perim / ( norm_constant * cell_area )
+        !--
+        !--With all this, we have
+        !-- a = SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) )
+        !-- P = RPI * a * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
+        !--and therefore, using the perimeter
+        !-- alpha * cldfra * ( 1. - cldfra ) * norm_constant * cell_area
+        !--             = N_cld_mix * RPI &
+        !--             * SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) ) &
+        !--             * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
+        !--and finally
+        N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - dcf_con(i) - cldfra(i) )**2. &
+                  * cell_area(i) * ( 1. - dcf_con(i) ) * bovera / RPI &
+                  / ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )**2.
+        !--where coef_mix_lscp = ( alpha * norm_constant )**2.
+        !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer
+        !--In particular, it is 0 if cldfra = 1
+        a_mix = SQRT( cell_area(i) * ( 1. - dcf_con(i) ) * cldfra(i) / bovera / N_cld_mix / RPI )
+
+        !--The time required for turbulent diffusion to homogenize a region of size
+        !--L_mix is defined as (L_mix**2./tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
+        !--We compute L_mix and assume that the cloud is mixed over this length
+        L_mix = SQRT( dtime**3. * pbl_eps(i) )
+        !--The mixing length cannot be greater than the semi-minor axis. In this case,
+        !--the entire cloud is mixed.
+        L_mix = MIN(L_mix, a_mix * bovera)
+
+        !--The fraction of clear sky mixed is
+        !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
+        envfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) &
+                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2. )
+        !--The fraction of cloudy sky mixed is
+        !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
+        cldfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) &
+                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2. )
+
+
+        !-- PART 2 - SHEARING
+
+        !--The clouds are then sheared. We keep the shape and number
+        !--assumptions from before. The clouds are sheared along their
+        !--semi-major axis (a_mix), on the entire cell heigh dz.
+        !--The increase in size is
+        L_shear = coef_shear_lscp * shear(i) * dz * dtime
+        !--therefore, the total increase in fraction is
+        shear_fra = RPI * L_shear * a_mix * bovera * N_cld_mix &
+                  / cell_area(i) / ( 1. - dcf_con(i) )
+        !--and the environment and cloud mixed fractions are the same,
+        !--which we add to the previous calculated mixed fractions.
+        !--We therefore assume that the sheared clouds and the turbulent
+        !--mixed clouds are different.
+        envfra_mix = envfra_mix + shear_fra
+        cldfra_mix = cldfra_mix + shear_fra
+
+
+        !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
+
+        !--The environment fraction is allocated to subsaturated sky or supersaturated sky,
+        !--according to the factor sigma_mix. This is computed as the ratio of the
+        !--subsaturated sky fraction to the environment fraction, corrected by a factor
+        !--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the
+        !--supersaturated sky is favoured. Physically, this means that it is more likely
+        !--to have supersaturated sky around the cloud than subsaturated sky.
+        sigma_mix = subfra(i) / ( subfra(i) + chi_mixing_lscp * issrfra(i) )
+        subfra_mix = MIN( sigma_mix * envfra_mix, subfra(i) )
+        issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra(i) )
+        cldfra_mix = MIN( cldfra_mix, cldfra(i) )
+
+        !--First, we mix the subsaturated sky (subfra_mix) and the cloud close
+        !--to this fraction (sigma_mix * cldfra_mix).
+        IF ( subfra(i) .GT. eps ) THEN
+
+          IF ( ok_unadjusted_clouds ) THEN
+            !--The subsaturated air is simply added to the cloud,
+            !--with the corresponding cloud fraction
+            !--If the cloud is too subsaturated, the sublimation process
+            !--activated in the following timestep will reduce the cloud
+            !--fraction
+            dcf_mix_sub = subfra_mix
+            dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
+            dqvc_mix_sub = dqt_mix_sub
+
+          ELSE
+            !--We compute the total humidity in the mixed air, which
+            !--can be either sub- or supersaturated.
+            qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &
+                      + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &
+                      / ( subfra_mix + cldfra_mix * sigma_mix )
+
+            IF ( qvapinmix .GT. qsat(i) ) THEN
+              !--If the mixed air is supersaturated, we condense the subsaturated
+              !--region which was mixed.
+              dcf_mix_sub = subfra_mix
+              dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
+              dqvc_mix_sub = dcf_mix_sub * qsat(i)
+            ELSE
+              !--Else, we sublimate the cloud which was mixed.
+              dcf_mix_sub = - sigma_mix * cldfra_mix
+              dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
+              dqvc_mix_sub = dcf_mix_sub * qsat(i)
+            ENDIF
+          ENDIF ! ok_unadjusted_clouds
+        ENDIF ! subfra .GT. eps
+   
+        !--We then mix the supersaturated sky (issrfra_mix) and the cloud,
+        !--for which the mixed air is always supersatured, therefore
+        !--the cloud necessarily expands
+        IF ( issrfra(i) .GT. eps ) THEN
+
+          IF ( ok_unadjusted_clouds ) THEN
+            !--The ice supersaturated air is simply added to the
+            !--cloud, and supersaturated vapor will be deposited on the
+            !--cloud ice crystals by the deposition process in the
+            !--following timestep
+            dcf_mix_issr = issrfra_mix
+            dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
+            dqvc_mix_issr = dqt_mix_issr
+          ELSE
+            !--In this case, the additionnal vapor condenses
+            dcf_mix_issr = issrfra_mix
+            dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
+            dqvc_mix_issr = dcf_mix_issr * qsat(i)
+          ENDIF ! ok_unadjusted_clouds
+
+
+        ENDIF ! issrfra .GT. eps
+
+        !--Sum up the tendencies from subsaturated sky and supersaturated sky
+        dcf_mix(i)  = dcf_mix_sub  + dcf_mix_issr
+        dqt_mix     = dqt_mix_sub  + dqt_mix_issr
+        dqvc_mix(i) = dqvc_mix_sub + dqvc_mix_issr
+        dqi_mix(i)  = dqt_mix - dqvc_mix(i)
+
+        !--Add tendencies
+        issrfra(i) = MAX(0., issrfra(i) - dcf_mix_issr)
+        qissr(i)   = MAX(0., qissr(i) - dqt_mix_issr)
+        cldfra(i)  = MAX(0., MIN(1. - dcf_con(i), cldfra(i) + dcf_mix(i)))
+        qcld(i)    = MAX(0., MIN(qtot(i) - dqi_con(i) - dqvc_con(i), qcld(i) + dqt_mix))
+        qvc(i)     = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i)))
+
+      ENDIF ! ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) )
+
+      !--Finally, we add the tendencies of condensation
+      cldfra(i) = MIN(1., cldfra(i) + dcf_con(i))
+      qcld(i)   = MIN(qtot(i), qcld(i) + dqvc_con(i) + dqi_con(i))
+      qvc(i)    = MIN(qcld(i), qvc(i) + dqvc_con(i))
+
+
+      !----------------------------------------
+      !--       CONTRAILS AND AVIATION       --
+      !----------------------------------------
+
+      !--Add a source of cirrus from aviation contrails
+      !IF ( ok_plane_contrail ) THEN
+      !  dcf_avi(i) = 0.
+      !  dqi_avi(i) = 0.
+      !  dqvc_avi(i) = 0.
+      !  ! TODO implement ok_unadjusted_clouds
+      !  IF ( issrfra(i) .GT. eps ) THEN
+      !    contrail_fra = MIN(1., flight_m(i,k) * dtime * contrail_cross_section / V_cell)
+      !    dcf_avi(i) = issrfra(i) * contrail_fra
+      !    dqt_avi = dcf_avi(i) * qissr(i) / issrfra(i)
+      !    dqvc_avi(i) = qsat(i) * dcf_avi(i)
+      !    
+      !    !--Add tendencies
+      !    cldfra(i) = cldfra(i) + dcf_avi(i)
+      !    issrfra(i) = issrfra(i) - dcf_avi(i)
+      !    qcld(i) = qcld(i) + dqt_avi
+      !    qvc(i) = qvc(i) + dqvc_avi(i)
+      !    qissr(i) = qissr(i) - dqt_avi
+
+      !    !--Diagnostics
+      !    dqi_avi(i) = dqt_avi - qsat(i) * dcf_avi(i)
+      !  ENDIF
+      !  dcf_avi(i)  = dcf_avi(i)  / dtime
+      !  dqi_avi(i)  = dqi_avi(i)  / dtime
+      !  dqvc_avi(i) = dqvc_avi(i) / dtime
+      !ENDIF
+
+
+
+      !-------------------------------------------
+      !--       FINAL BARRIERS AND OUTPUTS      --
+      !-------------------------------------------
+
+      IF ( cldfra(i) .LT. eps ) THEN
+        !--If the cloud is too small, it is sublimated.
+        cldfra(i) = 0.
+        qcld(i)   = 0.
+        qvc(i)    = 0.
+        qincld(i) = qsat(i)
+      ELSE
+        qincld(i) = qcld(i) / cldfra(i)
+      ENDIF ! cldfra .LT. eps
+
+      !--Diagnostics
+      dcf_sub(i)  = dcf_sub(i)  / dtime
+      dcf_con(i)  = dcf_con(i)  / dtime 
+      dcf_mix(i)  = dcf_mix(i)  / dtime 
+      dqi_adj(i)  = dqi_adj(i)  / dtime 
+      dqi_sub(i)  = dqi_sub(i)  / dtime 
+      dqi_con(i)  = dqi_con(i)  / dtime 
+      dqi_mix(i)  = dqi_mix(i)  / dtime 
+      dqvc_adj(i) = dqvc_adj(i) / dtime
+      dqvc_sub(i) = dqvc_sub(i) / dtime
+      dqvc_con(i) = dqvc_con(i) / dtime
+      dqvc_mix(i) = dqvc_mix(i) / dtime
+
+    ENDIF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds
+
+  ENDIF   ! end keepgoing
+
+ENDDO     ! end loop on i
+
+END SUBROUTINE condensation_ice_supersat
+!**********************************************************************************
+
+END MODULE lmdz_lscp_condensation
Index: /LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_ini.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_ini.F90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_ini.F90	(revision 4951)
@@ -6,6 +6,6 @@
   !--------------------
  
-  REAL RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
-  !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI)
+  REAL RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RV, RG, RPI, EPS_W
+  !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RV, RG, RPI, EPS_W)
   
   REAL, SAVE, PROTECTED :: seuil_neb=0.001      ! cloud fraction threshold: a cloud can precipitate when exceeded
@@ -136,6 +136,68 @@
   !$OMP THREADPRIVATE(dist_liq)
 
-  REAL, SAVE, PROTECTED    :: tresh_cl=0.0                  ! cloud fraction threshold for cloud top search
+  REAL, SAVE, PROTECTED  :: tresh_cl=0.0                  ! cloud fraction threshold for cloud top search
   !$OMP THREADPRIVATE(tresh_cl)
+
+  !--Parameters for condensation and ice supersaturation
+  LOGICAL, SAVE, PROTECTED :: ok_external_lognormal=.FALSE.  ! if True, the lognormal condensation scheme is calculated in the lmdz_lscp_condensation routine
+  !$OMP THREADPRIVATE(ok_external_lognormal)
+
+  LOGICAL, SAVE, PROTECTED :: ok_ice_supersat=.FALSE.        ! activates the condensation scheme that allows for ice supersaturation
+  !$OMP THREADPRIVATE(ok_ice_supersat)
+
+  LOGICAL, SAVE, PROTECTED :: ok_unadjusted_clouds=.FALSE.   ! if True, relax the saturation adjustment assumption for ice clouds
+  !$OMP THREADPRIVATE(ok_unadjusted_clouds)
+
+  LOGICAL, SAVE, PROTECTED :: ok_weibull_warm_clouds=.FALSE. ! if True, the weibull condensation scheme replaces the lognormal condensation scheme at positive temperatures
+  !$OMP THREADPRIVATE(ok_weibull_warm_clouds)
+
+  INTEGER, SAVE, PROTECTED :: iflag_cloud_sublim_pdf=3       ! iflag for the distribution of water inside ice clouds
+  !$OMP THREADPRIVATE(iflag_cloud_sublim_pdf)
+
+  REAL, SAVE, PROTECTED :: mu_subl_pdf_lscp=1./3.            ! [-] shape factor of the gamma distribution of water inside ice clouds
+  !$OMP THREADPRIVATE(mu_subl_pdf_lscp)
+  
+  REAL, SAVE, PROTECTED :: beta_pdf_lscp=2.8e-3              ! [] 
+  !$OMP THREADPRIVATE(beta_pdf_lscp)
+  
+  REAL, SAVE, PROTECTED :: temp_thresh_pdf_lscp=188.         ! [K] below this temperature, water vapor is homoegeneously distributed in the clear sky region
+  !$OMP THREADPRIVATE(temp_thresh_pdf_lscp)
+  
+  REAL, SAVE, PROTECTED :: rhlmid_pdf_lscp=57.6              ! [%] 
+  !$OMP THREADPRIVATE(rhlmid_pdf_lscp)
+  
+  REAL, SAVE, PROTECTED :: k0_pdf_lscp=1.75                  ! [-] minimum value of the shape parameter for the weibull law
+  !$OMP THREADPRIVATE(k0_pdf_lscp)
+  
+  REAL, SAVE, PROTECTED :: kappa_pdf_lscp=0.0147             ! [] 
+  !$OMP THREADPRIVATE(kappa_pdf_lscp)
+  
+  REAL, SAVE, PROTECTED :: rhl0_pdf_lscp=81.5                ! [%] 
+  !$OMP THREADPRIVATE(rhl0_pdf_lscp)
+  
+  REAL, SAVE, PROTECTED :: cond_thresh_pdf_lscp=1.E-7        ! [-] 
+  !$OMP THREADPRIVATE(cond_thresh_pdf_lscp)
+  
+  REAL, SAVE, PROTECTED :: a_homofreez=2.349                 ! [-] 
+  !$OMP THREADPRIVATE(a_homofreez)
+  
+  REAL, SAVE, PROTECTED :: b_homofreez=259.                  ! [K] 
+  !$OMP THREADPRIVATE(b_homofreez)
+
+  REAL, SAVE, PROTECTED :: delta_hetfreez=1.                 ! [-] value between 0 and 1 to simulate for heterogeneous freezing.
+  !$OMP THREADPRIVATE(delta_hetfreez)
+  
+  REAL, SAVE, PROTECTED :: coef_mixing_lscp=1e-7             ! [-] 
+  !$OMP THREADPRIVATE(coef_mixing_lscp)
+  
+  REAL, SAVE, PROTECTED :: coef_shear_lscp=0.1               ! [-] 
+  !$OMP THREADPRIVATE(coef_shear_lscp)
+  
+  REAL, SAVE, PROTECTED :: chi_mixing_lscp=1.1               ! [-] 
+  !$OMP THREADPRIVATE(chi_mixing_lscp)
+
+!  REAL, SAVE, PROTECTED :: contrail_cross_section=200000.
+!  !$OMP THREADPRIVATE(contrail_cross_section)
+  !--End of the parameters for condensation and ice supersaturation
 
   !--Parameters for poprecip
@@ -225,18 +287,18 @@
 CONTAINS
 
-SUBROUTINE lscp_ini(dtime,lunout_in,prt_level_in,ok_ice_sursat, iflag_ratqs, fl_cor_ebil_in, RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, &
-                    RVTMP2_in, RTT_in,RD_in,RG_in,RPI_in)
+SUBROUTINE lscp_ini(dtime,lunout_in,prt_level_in,ok_ice_supersat_in, iflag_ratqs, fl_cor_ebil_in, &
+                    RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, RVTMP2_in, &
+                    RTT_in, RD_in, RV_in, RG_in, RPI_in, EPS_W_in)
 
 
    USE ioipsl_getin_p_mod, ONLY : getin_p
-   USE ice_sursat_mod, ONLY: ice_sursat_init
    USE lmdz_cloudth_ini, ONLY : cloudth_ini
 
    REAL, INTENT(IN)      :: dtime
    INTEGER, INTENT(IN)   :: lunout_in,prt_level_in,iflag_ratqs,fl_cor_ebil_in
-   LOGICAL, INTENT(IN)   :: ok_ice_sursat
+   LOGICAL, INTENT(IN)   :: ok_ice_supersat_in
 
    REAL, INTENT(IN)      :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in
-   REAL, INTENT(IN)      ::  RVTMP2_in, RTT_in, RD_in, RG_in, RPI_in
+   REAL, INTENT(IN)      :: RVTMP2_in, RTT_in, RD_in, RV_in, RG_in, RPI_in, EPS_W_in
    character (len=20) :: modname='lscp_ini_mod'
    character (len=80) :: abort_message
@@ -247,6 +309,9 @@
     fl_cor_ebil=fl_cor_ebil_in
 
+    ok_ice_supersat=ok_ice_supersat_in
+
     RG=RG_in
     RD=RD_in
+    RV=RV_in
     RCPD=RCPD_in
     RLVTT=RLVTT_in
@@ -257,4 +322,5 @@
     RVTMP2=RVTMP2_in
     RPI=RPI_in
+    EPS_W=EPS_W_in
 
 
@@ -297,4 +363,5 @@
     CALL getin_p('tresh_cl',tresh_cl)
     CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
+    ! for poprecip
     CALL getin_p('ok_poprecip',ok_poprecip)
     CALL getin_p('ok_corr_vap_evasub',ok_corr_vap_evasub)
@@ -316,4 +383,24 @@
     CALL getin_p('snow_fallspeed_clr',snow_fallspeed_clr)
     CALL getin_p('snow_fallspeed_cld',snow_fallspeed_cld)
+    ! for condensation and ice supersaturation
+    CALL getin_p('ok_external_lognormal',ok_external_lognormal)
+    CALL getin_p('ok_unadjusted_clouds',ok_unadjusted_clouds)
+    CALL getin_p('ok_weibull_warm_clouds',ok_weibull_warm_clouds)
+    CALL getin_p('iflag_cloud_sublim_pdf',iflag_cloud_sublim_pdf)
+    CALL getin_p('mu_subl_pdf_lscp',mu_subl_pdf_lscp)
+    CALL getin_p('beta_pdf_lscp',beta_pdf_lscp)
+    CALL getin_p('temp_thresh_pdf_lscp',temp_thresh_pdf_lscp)
+    CALL getin_p('rhlmid_pdf_lscp',rhlmid_pdf_lscp)
+    CALL getin_p('k0_pdf_lscp',k0_pdf_lscp)
+    CALL getin_p('kappa_pdf_lscp',kappa_pdf_lscp)
+    CALL getin_p('rhl0_pdf_lscp',rhl0_pdf_lscp)
+    CALL getin_p('cond_thresh_pdf_lscp',cond_thresh_pdf_lscp)
+    CALL getin_p('a_homofreez',a_homofreez)
+    CALL getin_p('b_homofreez',b_homofreez)
+    CALL getin_p('delta_hetfreez',delta_hetfreez)
+    CALL getin_p('coef_mixing_lscp',coef_mixing_lscp)
+    CALL getin_p('coef_shear_lscp',coef_shear_lscp)
+    CALL getin_p('chi_mixing_lscp',chi_mixing_lscp)
+    !CALL getin_p('contrail_cross_section',contrail_cross_section)
 
 
@@ -354,4 +441,5 @@
     WRITE(lunout,*) 'lscp_ini, iflag_oldbug_fisrtilp', iflag_oldbug_fisrtilp
     WRITE(lunout,*) 'lscp_ini, fl_cor_ebil', fl_cor_ebil
+    ! for poprecip
     WRITE(lunout,*) 'lscp_ini, ok_poprecip', ok_poprecip
     WRITE(lunout,*) 'lscp_ini, ok_corr_vap_evasub', ok_corr_vap_evasub
@@ -367,4 +455,24 @@
     WRITE(lunout,*) 'lscp_ini, snow_fallspeed_clr:', snow_fallspeed_clr
     WRITE(lunout,*) 'lscp_ini, snow_fallspeed_cld:', snow_fallspeed_cld
+    ! for condensation and ice supersaturation
+    WRITE(lunout,*) 'lscp_ini, ok_external_lognormal:', ok_external_lognormal
+    WRITE(lunout,*) 'lscp_ini, ok_ice_supersat:', ok_ice_supersat
+    WRITE(lunout,*) 'lscp_ini, ok_unadjusted_clouds:', ok_unadjusted_clouds
+    WRITE(lunout,*) 'lscp_ini, ok_weibull_warm_clouds:', ok_weibull_warm_clouds
+    WRITE(lunout,*) 'lscp_ini, iflag_cloud_sublim_pdf:', iflag_cloud_sublim_pdf
+    WRITE(lunout,*) 'lscp_ini, mu_subl_pdf_lscp:', mu_subl_pdf_lscp
+    WRITE(lunout,*) 'lscp_ini, beta_pdf_lscp:', beta_pdf_lscp
+    WRITE(lunout,*) 'lscp_ini, temp_thresh_pdf_lscp:', temp_thresh_pdf_lscp
+    WRITE(lunout,*) 'lscp_ini, rhlmid_pdf_lscp:', rhlmid_pdf_lscp
+    WRITE(lunout,*) 'lscp_ini, k0_pdf_lscp:', k0_pdf_lscp
+    WRITE(lunout,*) 'lscp_ini, kappa_pdf_lscp:', kappa_pdf_lscp
+    WRITE(lunout,*) 'lscp_ini, rhl0_pdf_lscp:', rhl0_pdf_lscp
+    WRITE(lunout,*) 'lscp_ini, a_homofreez:', a_homofreez
+    WRITE(lunout,*) 'lscp_ini, b_homofreez:', b_homofreez
+    WRITE(lunout,*) 'lscp_ini, delta_hetfreez', delta_hetfreez
+    WRITE(lunout,*) 'lscp_ini, coef_mixing_lscp:', coef_mixing_lscp
+    WRITE(lunout,*) 'lscp_ini, coef_shear_lscp:', coef_shear_lscp
+    WRITE(lunout,*) 'lscp_ini, chi_mixing_lscp:', chi_mixing_lscp
+!    WRITE(lunout,*) 'lscp_ini, contrail_cross_section:', contrail_cross_section
 
 
@@ -389,4 +497,21 @@
 
 
+    !--Check flags for condensation and ice supersaturation
+    IF ( ok_external_lognormal .AND. ok_ice_supersat ) THEN
+      abort_message = 'in lscp, ok_external_lognormal=y is incompatible with ok_ice_supersat=y'
+      CALL abort_physic (modname,abort_message,1)
+    ENDIF
+
+    IF ( ok_weibull_warm_clouds .AND. .NOT. ok_ice_supersat ) THEN
+      abort_message = 'in lscp, ok_weibull_warm_clouds=y needs ok_ice_supersat=y'
+      CALL abort_physic (modname,abort_message,1)
+    ENDIF
+
+    IF ( ok_unadjusted_clouds .AND. .NOT. ok_ice_supersat ) THEN
+      abort_message = 'in lscp, ok_unadjusted_clouds=y needs ok_ice_supersat=y'
+      CALL abort_physic (modname,abort_message,1)
+    ENDIF
+
+
     !AA Temporary initialisation
     a_tr_sca(1) = -0.5
@@ -395,6 +520,4 @@
     a_tr_sca(4) = -0.5
     
-    IF (ok_ice_sursat) CALL ice_sursat_init()
-
     CALL cloudth_ini(iflag_cloudth_vert,iflag_ratqs)
 
Index: /LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_tools.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_tools.F90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_tools.F90	(revision 4951)
@@ -311,5 +311,6 @@
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
-    use lmdz_lscp_ini, only: iflag_gammasat, t_glace_min, RTT
+    use lmdz_lscp_ini, only: iflag_gammasat, temp_nowater, RTT
+    use lmdz_lscp_ini, only: a_homofreez, b_homofreez, delta_hetfreez
 
     IMPLICIT NONE
@@ -326,7 +327,5 @@
 
     REAL, DIMENSION(klon) ::  qsi,qsl,dqsl,dqsi
-    REAL  fcirrus, fac
-    REAL, PARAMETER :: acirrus=2.349
-    REAL, PARAMETER :: bcirrus=259.0
+    REAL  f_homofreez, fac
 
     INTEGER i
@@ -342,5 +341,5 @@
             dgammasatdt(i)=0.
 
-        ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN
+        ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. temp_nowater)) THEN
             
             IF (iflag_gammasat .GE. 2) THEN          
@@ -355,14 +354,14 @@
 
             IF (iflag_gammasat .GE.1) THEN
-            ! homogeneous freezing of aerosols, according to
-            ! Koop, 2000 and Karcher 2008, QJRMS
-            ! 'Cirrus regime'
-               fcirrus=acirrus-temp(i)/bcirrus
-               IF (fcirrus .GT. qsl(i)/qsi(i)) 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)=fcirrus
-                  dgammasatdt(i)=-1.0/bcirrus
+                  gammasat(i)= (1.-delta_hetfreez) + delta_hetfreez * f_homofreez
+                  dgammasatdt(i)=-delta_hetfreez/b_homofreez
                ENDIF
            
@@ -452,4 +451,180 @@
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+FUNCTION GAMMAINC ( p, x )
+
+!*****************************************************************************80
+!
+!! GAMMAINC computes the regularized lower incomplete Gamma Integral
+!
+!  Modified:
+!
+!    20 January 2008
+!
+!  Author:
+!
+!    Original FORTRAN77 version by B Shea.
+!    FORTRAN90 version by John Burkardt.
+!
+!  Reference:
+!
+!    B Shea,
+!    Algorithm AS 239:
+!    Chi-squared and Incomplete Gamma Integral,
+!    Applied Statistics,
+!    Volume 37, Number 3, 1988, pages 466-473.
+!
+!  Parameters:
+!
+!    Input, real X, P, the parameters of the incomplete
+!    gamma ratio.  0 <= X, and 0 < P.
+!
+!    Output, real GAMMAINC, the value of the incomplete
+!    Gamma integral.
+!
+  IMPLICIT NONE
+
+  REAL A
+  REAL AN
+  REAL ARG
+  REAL B
+  REAL C
+  REAL, PARAMETER :: ELIMIT = - 88.0E+00
+  REAL GAMMAINC
+  REAL, PARAMETER :: OFLO = 1.0E+37
+  REAL P
+  REAL, PARAMETER :: PLIMIT = 1000.0E+00
+  REAL PN1
+  REAL PN2
+  REAL PN3
+  REAL PN4
+  REAL PN5
+  REAL PN6
+  REAL RN
+  REAL, PARAMETER :: TOL = 1.0E-14
+  REAL X
+  REAL, PARAMETER :: XBIG = 1.0E+08
+
+  GAMMAINC = 0.0E+00
+
+  IF ( X == 0.0E+00 ) THEN
+    GAMMAINC = 0.0E+00
+    RETURN
+  END IF
+!
+!  IF P IS LARGE, USE A NORMAL APPROXIMATION.
+!
+  IF ( PLIMIT < P ) THEN
+
+    PN1 = 3.0E+00 * SQRT ( P ) * ( ( X / P )**( 1.0E+00 / 3.0E+00 ) &
+    + 1.0E+00 / ( 9.0E+00 * P ) - 1.0E+00 )
+
+    GAMMAINC = 0.5E+00 * ( 1. + ERF ( PN1 ) )
+    RETURN
+
+  END IF
+!
+!  IF X IS LARGE SET GAMMAD = 1.
+!
+  IF ( XBIG < X ) THEN
+    GAMMAINC = 1.0E+00
+    RETURN
+  END IF
+!
+!  USE PEARSON'S SERIES EXPANSION.
+!  (NOTE THAT P IS NOT LARGE ENOUGH TO FORCE OVERFLOW IN ALOGAM).
+!
+  IF ( X <= 1.0E+00 .OR. X < P ) THEN
+
+    ARG = P * LOG ( X ) - X - LOG_GAMMA ( P + 1.0E+00 )
+    C = 1.0E+00
+    GAMMAINC = 1.0E+00
+    A = P
+
+    DO
+
+      A = A + 1.0E+00
+      C = C * X / A
+      GAMMAINC = GAMMAINC + C
+
+      IF ( C <= TOL ) THEN
+        EXIT
+      END IF
+
+    END DO
+
+    ARG = ARG + LOG ( GAMMAINC )
+
+    IF ( ELIMIT <= ARG ) THEN
+      GAMMAINC = EXP ( ARG )
+    ELSE
+      GAMMAINC = 0.0E+00
+    END IF
+!
+!  USE A CONTINUED FRACTION EXPANSION.
+!
+  ELSE
+
+    ARG = P * LOG ( X ) - X - LOG_GAMMA ( P )
+    A = 1.0E+00 - P
+    B = A + X + 1.0E+00
+    C = 0.0E+00
+    PN1 = 1.0E+00
+    PN2 = X
+    PN3 = X + 1.0E+00
+    PN4 = X * B
+    GAMMAINC = PN3 / PN4
+
+    DO
+
+      A = A + 1.0E+00
+      B = B + 2.0E+00
+      C = C + 1.0E+00
+      AN = A * C
+      PN5 = B * PN3 - AN * PN1
+      PN6 = B * PN4 - AN * PN2
+
+      IF ( PN6 /= 0.0E+00 ) THEN
+
+        RN = PN5 / PN6
+
+        IF ( ABS ( GAMMAINC - RN ) <= MIN ( TOL, TOL * RN ) ) THEN
+          EXIT
+        END IF
+
+        GAMMAINC = RN
+
+      END IF
+
+      PN1 = PN3
+      PN2 = PN4
+      PN3 = PN5
+      PN4 = PN6
+!
+!  RE-SCALE TERMS IN CONTINUED FRACTION IF TERMS ARE LARGE.
+!
+      IF ( OFLO <= ABS ( PN5 ) ) THEN
+        PN1 = PN1 / OFLO
+        PN2 = PN2 / OFLO
+        PN3 = PN3 / OFLO
+        PN4 = PN4 / OFLO
+      END IF
+
+    END DO
+
+    ARG = ARG + LOG ( GAMMAINC )
+
+    IF ( ELIMIT <= ARG ) THEN
+      GAMMAINC = 1.0E+00 - EXP ( ARG )
+    ELSE
+      GAMMAINC = 1.0E+00
+    END IF
+
+  END IF
+
+  RETURN
+END FUNCTION GAMMAINC
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
 END MODULE lmdz_lscp_tools
 
Index: /LMDZ6/branches/cirrus/libf/phylmd/phyetat0_mod.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/phyetat0_mod.F90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/phylmd/phyetat0_mod.F90	(revision 4951)
@@ -21,5 +21,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, &
@@ -397,5 +398,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.)
@@ -412,4 +412,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
@@ -419,5 +428,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. &
@@ -432,4 +440,11 @@
        ancien_ok=.false.
     ENDIF
+  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
 
Index: /LMDZ6/branches/cirrus/libf/phylmd/phyredem.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/phyredem.F90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/phylmd/phyredem.F90	(revision 4951)
@@ -19,6 +19,7 @@
                                 zval, rugoro, t_ancien, q_ancien,            &
                                 prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien,      &
-                                ql_ancien, qs_ancien, qbs_ancien, rneb_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, awake_s,   &
@@ -252,5 +253,8 @@
     ENDIF
 
-    CALL put_field(pass,"RNEBANCIEN", "RNEBANCIEN", rneb_ancien)
+    IF ( ok_ice_supersat ) THEN
+      CALL put_field(pass,"CFANCIEN", "CFANCIEN", cf_ancien)
+      CALL put_field(pass,"RVCANCIEN", "RVCANCIEN", rvc_ancien)
+    ENDIF
 
     CALL put_field(pass,"PRWANCIEN", "PRWANCIEN", prw_ancien)
Index: /LMDZ6/branches/cirrus/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/phys_local_var_mod.F90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/phylmd/phys_local_var_mod.F90	(revision 4951)
@@ -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)
@@ -494,44 +494,33 @@
 !$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)
+!-- 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)
@@ -539,6 +528,5 @@
       !$OMP THREADPRIVATE(temp_cltop)
 
-
-!--POPRECIP variables
+!-- LSCP - POPRECIP variables
       REAL, SAVE, ALLOCATABLE :: qraindiag(:,:)
       !$OMP THREADPRIVATE(qraindiag)
@@ -668,4 +656,5 @@
       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))
@@ -678,4 +667,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))
@@ -954,15 +944,20 @@
       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
+!-- 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))
@@ -1022,4 +1017,5 @@
       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)
@@ -1030,4 +1026,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)
@@ -1270,15 +1267,20 @@
       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
+!-- 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)
Index: /LMDZ6/branches/cirrus/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 4951)
@@ -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) /))
@@ -2069,51 +2067,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/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/phys_output_write_mod.F90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/phylmd/phys_output_write_mod.F90	(revision 4951)
@@ -54,4 +54,5 @@
          o_SWdnTOAclr, o_nettop, o_SWup200, &
          o_SWup200clr, o_SWdn200, o_SWdn200clr, &
+         o_LWupTOA, o_LWupTOAclr, &
          o_LWup200, o_LWup200clr, o_LWdn200, &
          o_LWdn200clr, o_sols, o_sols0, &
@@ -216,9 +217,12 @@
          o_p_tropopause, o_z_tropopause, o_t_tropopause,  &
          o_col_O3_strato, o_col_O3_tropo,                 & 
-!--aviation & supersaturation
-         o_oclr, o_ocld, o_oss, o_ovc, o_rnebss, o_rnebclr, o_rnebseri, o_gammass, &
-         o_N1_ss, o_N2_ss, o_qsatl, o_qsats, &
-         o_drnebsub, o_drnebcon, o_drnebtur, o_drnebavi, o_flight_m, o_flight_h2o, &
+!-- LSCP - condensation and ice supersaturation variables
+         o_cfseri, o_dcfdyn, o_rvcseri, o_drvcdyn, &
+         o_qsub, o_qissr, o_qcld, o_subfra, o_issrfra, o_gammacond, &
+         o_dcfsub, o_dcfcon, o_dcfmix, o_dqiadj, o_dqisub, o_dqicon, o_dqimix, &
+         o_dqvcadj, o_dqvcsub, o_dqvccon, o_dqvcmix, o_qsatl, o_qsati, &
+!-- LSCP - aviation variables
          o_Tcontr, o_qcontr, o_qcontr2, o_fcontrN, o_fcontrP, &
+         o_dcfavi, o_dqiavi, o_dqvcavi, o_flight_dist, o_flight_h2o, &
 !--interactive CO2
          o_flx_co2_ocean, o_flx_co2_ocean_cor, &
@@ -258,5 +262,4 @@
 #endif
 
-    USE ice_sursat_mod, ONLY: flight_m, flight_h2o
     USE lmdz_lscp_ini, ONLY: ok_poprecip
 
@@ -323,5 +326,4 @@
          cldq, flwp, fiwp, ue, ve, uq, vq, &
          uwat, vwat, &
-         rneb_seri, d_rneb_dyn, &
          plcl, plfc, wbeff, convoccur, upwd, dnwd, dnwd0, prw, prlw, prsw, prbsw, water_budget, &
          s_pblh, s_pblt, s_lcl, s_therm, uwriteSTD, &
@@ -337,8 +339,12 @@
          wdtrainA, wdtrainS, wdtrainM, n2, s2, strig, zcong, zlcl_th, proba_notrig, &
          random_notrig, &
-         qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
-         N1_ss, N2_ss, zqsatl, zqsats, &
+         cf_seri, d_cf_dyn, rvc_seri, d_rvc_dyn, &
+         qsub, qissr, qcld, subfra, issrfra, gamma_cond, &
+         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, &
-         drneb_sub, drneb_con, drneb_tur, drneb_avi, &
+         dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
          alp_bl_det, alp_bl_fluct_m, alp_bl_conv, &
          alp_bl_stat, alp_bl_fluct_tke, slab_wfbils, &
@@ -1066,4 +1072,14 @@
        CALL histwrite_phy(o_LWupSFCclr, zx_tmp_fi2d)
        CALL histwrite_phy(o_LWdnSFCclr, sollwdownclr)
+       
+       IF (vars_defined) THEN
+          zx_tmp_fi2d(:) = lwup(:,klevp1)
+       ENDIF
+       CALL histwrite_phy(o_LWupTOA, zx_tmp_fi2d)
+       
+       IF (vars_defined) THEN
+          zx_tmp_fi2d(:) = lwup0(:,klevp1)
+       ENDIF
+       CALL histwrite_phy(o_LWupTOAclr, zx_tmp_fi2d)
 
        IF (type_trac/='inca' .OR. config_inca=='aeNP') THEN 
@@ -2006,22 +2022,32 @@
        ENDIF
 
-!--aviation & supersaturation
-       IF (ok_ice_sursat) THEN
-         CALL histwrite_phy(o_oclr, qclr)
-         CALL histwrite_phy(o_ocld, qcld)
-         CALL histwrite_phy(o_oss, qss)
-         CALL histwrite_phy(o_ovc, qvc)
-         CALL histwrite_phy(o_rnebclr, rnebclr)
-         CALL histwrite_phy(o_rnebss, rnebss)
-         CALL histwrite_phy(o_rnebseri, rneb_seri)
-         CALL histwrite_phy(o_gammass, gamma_ss)
-         CALL histwrite_phy(o_N1_ss, N1_ss)
-         CALL histwrite_phy(o_N2_ss, N2_ss)
-         CALL histwrite_phy(o_drnebsub, drneb_sub)
-         CALL histwrite_phy(o_drnebcon, drneb_con)
-         CALL histwrite_phy(o_drnebtur, drneb_tur)
-         CALL histwrite_phy(o_drnebavi, drneb_avi)
-         CALL histwrite_phy(o_qsatl, zqsatl)
-         CALL histwrite_phy(o_qsats, zqsats)
+!-- LSCP - condensation and supersaturation variables
+       IF (ok_ice_supersat) THEN
+         CALL histwrite_phy(o_cfseri, cf_seri)
+         CALL histwrite_phy(o_dcfdyn, d_cf_dyn)
+         CALL histwrite_phy(o_rvcseri, rvc_seri)
+         CALL histwrite_phy(o_drvcdyn, d_rvc_dyn)
+         CALL histwrite_phy(o_qsub, qsub)
+         CALL histwrite_phy(o_qissr, qissr)
+         CALL histwrite_phy(o_qcld, qcld)
+         CALL histwrite_phy(o_subfra, subfra)
+         CALL histwrite_phy(o_issrfra, issrfra)
+         CALL histwrite_phy(o_gammacond, gamma_cond)
+         CALL histwrite_phy(o_dcfsub, dcf_sub)
+         CALL histwrite_phy(o_dcfcon, dcf_con)
+         CALL histwrite_phy(o_dcfmix, dcf_mix)
+         CALL histwrite_phy(o_dqiadj, dqi_adj)
+         CALL histwrite_phy(o_dqisub, dqi_sub)
+         CALL histwrite_phy(o_dqicon, dqi_con)
+         CALL histwrite_phy(o_dqimix, dqi_mix)
+         CALL histwrite_phy(o_dqvcadj, dqvc_adj)
+         CALL histwrite_phy(o_dqvcsub, dqvc_sub)
+         CALL histwrite_phy(o_dqvccon, dqvc_con)
+         CALL histwrite_phy(o_dqvcmix, dqvc_mix)
+         CALL histwrite_phy(o_qsatl, qsatliq)
+         CALL histwrite_phy(o_qsati, qsatice)
+       ENDIF
+!-- LSCP - aviation variables
+       IF (ok_plane_contrail) THEN
          CALL histwrite_phy(o_Tcontr, Tcontr)
          CALL histwrite_phy(o_qcontr, qcontr)
@@ -2029,7 +2055,8 @@
          CALL histwrite_phy(o_fcontrN, fcontrN)
          CALL histwrite_phy(o_fcontrP, fcontrP)
-       ENDIF
-       IF (ok_plane_contrail) THEN
-         CALL histwrite_phy(o_flight_m, flight_m)
+         CALL histwrite_phy(o_dcfavi, dcf_avi)
+         CALL histwrite_phy(o_dqiavi, dqi_avi)
+         CALL histwrite_phy(o_dqvcavi, dqvc_avi)
+         CALL histwrite_phy(o_flight_dist, flight_dist)
        ENDIF
        IF (ok_plane_h2o) THEN
Index: /LMDZ6/branches/cirrus/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/phys_state_var_mod.F90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/phylmd/phys_state_var_mod.F90	(revision 4951)
@@ -92,4 +92,6 @@
       REAL, ALLOCATABLE, SAVE :: u_ancien(:,:), v_ancien(:,:)
 !$OMP THREADPRIVATE(u_ancien, v_ancien)
+      REAL, ALLOCATABLE, SAVE :: cf_ancien(:,:), rvc_ancien(:,:)
+!$OMP THREADPRIVATE(cf_ancien, rvc_ancien)
 !!! RomP >>>
       REAL, ALLOCATABLE, SAVE :: tr_ancien(:,:,:)
@@ -100,6 +102,4 @@
       REAL, ALLOCATABLE, SAVE :: clwcon(:,:),rnebcon(:,:)
 !$OMP THREADPRIVATE(clwcon,rnebcon)
-      REAL, ALLOCATABLE, SAVE :: rneb_ancien(:,:)
-!$OMP THREADPRIVATE(rneb_ancien)
       REAL, ALLOCATABLE, SAVE :: qtc_cv(:,:),sigt_cv(:,:),detrain_cv(:,:),fm_cv(:,:)
 !$OMP THREADPRIVATE(qtc_cv,sigt_cv,detrain_cv,fm_cv)
@@ -586,9 +586,9 @@
       ALLOCATE(prw_ancien(klon), prlw_ancien(klon), prsw_ancien(klon), prbsw_ancien(klon))
       ALLOCATE(u_ancien(klon,klev), v_ancien(klon,klev))
+      ALLOCATE(cf_ancien(klon,klev), rvc_ancien(klon,klev))
 !!! Rom P >>>
       ALLOCATE(tr_ancien(klon,klev,nbtr))
 !!! Rom P <<<
       ALLOCATE(clwcon(klon,klev),rnebcon(klon,klev))
-      ALLOCATE(rneb_ancien(klon,klev))
       ALLOCATE(qtc_cv(klon,klev),sigt_cv(klon,klev),detrain_cv(klon,klev),fm_cv(klon,klev))
       ALLOCATE(ratqs(klon,klev))
@@ -809,8 +809,9 @@
       DEALLOCATE(zthe, zpic, zval)
       DEALLOCATE(rugoro, t_ancien, q_ancien, clwcon, rnebcon)
-      DEALLOCATE(qs_ancien, ql_ancien, qbs_ancien, rneb_ancien)
+      DEALLOCATE(qs_ancien, ql_ancien, qbs_ancien)
       DEALLOCATE(prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien)
       DEALLOCATE(qtc_cv,sigt_cv,detrain_cv,fm_cv)
       DEALLOCATE(u_ancien, v_ancien)
+      DEALLOCATE(cf_ancien, rvc_ancien)
       DEALLOCATE(tr_ancien)                           !RomP
       DEALLOCATE(ratqs, pbl_tke,coefh,coefm)
Index: /LMDZ6/branches/cirrus/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/branches/cirrus/libf/phylmd/physiq_mod.F90	(revision 4950)
+++ /LMDZ6/branches/cirrus/libf/phylmd/physiq_mod.F90	(revision 4951)
@@ -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
@@ -157,8 +156,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
@@ -334,7 +335,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,  &
@@ -511,6 +517,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)
     !
     !
@@ -1362,6 +1368,7 @@
        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'))
+       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)
@@ -1416,25 +1423,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)
@@ -1442,5 +1449,5 @@
 
         IF (ok_bs) THEN
-         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
@@ -1870,5 +1877,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)
@@ -2342,4 +2350,6 @@
                     sollwdown(:))
 
+      !--Init for LSCP - condensation
+      ratio_qi_qtot(:,:) = 0.
 
 
@@ -2439,16 +2449,17 @@
           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
@@ -2526,5 +2537,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
@@ -2538,6 +2551,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
     ELSE
        d_u_dyn(:,:)  = 0.0
@@ -2547,4 +2558,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
@@ -2554,6 +2568,4 @@
        IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
        ! !! RomP <<<
-       d_rneb_dyn(:,:)=0.0
-       d_qbs_dyn(:,:)=0.0
        ancien_ok = .TRUE.
     ENDIF
@@ -2664,4 +2676,20 @@
        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
+
     !
     ! Re-evaporer l'eau liquide nuageuse
@@ -3888,11 +3916,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, &
@@ -3900,7 +3928,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, &
@@ -5623,7 +5655,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
 
@@ -5659,6 +5692,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(:,:)
     CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
     CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
