Index: LMDZ6/branches/contrails/DefLists/field_def_lmdz.xml
===================================================================
--- LMDZ6/branches/contrails/DefLists/field_def_lmdz.xml	(revision 5602)
+++ LMDZ6/branches/contrails/DefLists/field_def_lmdz.xml	(revision 5609)
@@ -678,18 +678,4 @@
         <field id="tke_max_sic"    long_name="Max Turb. Kinetic Energy sic"    unit="m2/s2" operation="maximum"/>
 	
-	<field id="qrain_lsc"    long_name="LS specific rain content"    unit="kg/kg" />
-	<field id="qsnow_lsc"    long_name="LS specific snow content"    unit="kg/kg" />
-	<field id="dqreva"    long_name="LS rain tendency due to evaporation"    unit="kg/kg/s" />
-	<field id="dqssub"    long_name="LS snow tendency due to sublimation"    unit="kg/kg/s" />
-	<field id="dqrauto"    long_name="LS rain tendency due to autoconversion"    unit="kg/kg/s" />
-	<field id="dqrcol"    long_name="LS rain tendency due to collection"    unit="kg/kg/s" />
-	<field id="dqrmelt"    long_name="LS rain tendency due to melting"    unit="kg/kg/s" />
-	<field id="dqrfreez"    long_name="LS rain tendency due to freezing"    unit="kg/kg/s" />
-	<field id="dqsauto"    long_name="LS snow tendency due to autoconversion"    unit="kg/kg/s" />
-	<field id="dqsagg"    long_name="LS snow tendency due to aggregation"    unit="kg/kg/s" />
-	<field id="dqsrim"    long_name="LS snow tendency due to riming"    unit="kg/kg/s" />
-	<field id="dqsmelt"    long_name="LS snow tendency due to melting"    unit="kg/kg/s" />
-	<field id="dqsfreez"    long_name="LS snow tendency due to freezing"    unit="kg/kg/s" />
-
 	<field id="l_mix_ter"    long_name="PBL mixing length ter"    unit="m" />
         <field id="l_mix_lic"    long_name="PBL mixing length lic"    unit="m" />
@@ -909,24 +895,31 @@
         <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="dcfsed"     long_name="Sedimentation 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="dqised"     long_name="Sedimentation 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="dqvcsed"    long_name="Sedimentation 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="cfaseri"    long_name="Linear contrail fraction" unit="-" />
-	    <field id="dcfadyn"    long_name="Dynamics linear contrail fraction tendency" unit="s-1" />
-	    <field id="qtaseri"    long_name="Linear contrail total specific humidity" unit="s-1" />
-        <field id="dqtadyn"    long_name="Dynamics linear contrail total specific humidity tendency" unit="kg/kg/s" />
-	    <field id="Tcritcont"  long_name="Temperature threshold for contrail formation"    unit="K" />
-	    <field id="qcritcont"  long_name="Specific humidity threshold for contrail formation"    unit="kg/kg" />
+        <field id="cfaseri"    long_name="Linear contrail fraction" unit="-" />
+        <field id="dcfadyn"    long_name="Dynamics linear contrail fraction tendency" unit="s-1" />
+        <field id="pcfseri"    long_name="Contrail induced cirrus fraction" unit="-" />
+        <field id="dpcfdyn"    long_name="Dynamics contrail induced cirrus fraction tendency" unit="s-1" />
+        <field id="qvaseri"    long_name="Linear contrail vapor specific humidity" unit="kg/kg" />
+        <field id="dqvadyn"    long_name="Dynamics linear contrail vapor specific humidity tendency" unit="kg/kg/s" />
+        <field id="qiaseri"    long_name="Linear contrail ice specific humidity" unit="kg/kg" />
+        <field id="dqiadyn"    long_name="Dynamics linear contrail ice specific humidity tendency" unit="kg/kg/s" />
+        <field id="Tcritcont"  long_name="Temperature threshold for contrail formation"    unit="K" />
+        <field id="qcritcont"  long_name="Specific humidity threshold for contrail formation"    unit="kg/kg" />
         <field id="potcontfraP"  long_name="Fraction with pontential persistent contrail"    unit="-" />
         <field id="potcontfraNP" long_name="Fraction with potential non-persistent contrail"    unit="-" />
-        <field id="qice_cont"  long_name="Contrails ice specific humidity"    unit="kg/kg" />
+        <field id="qice_perscont" long_name="Contrail induced cirrus ice specific humidity"    unit="kg/kg" />
         <field id="dcfaini"    long_name="Initial formation linear contrail fraction tendency"    unit="s-1" />
         <field id="dqiaini"    long_name="Initial formation linear contrail ice specific humidity tendency"    unit="kg/kg/s" />
Index: LMDZ6/branches/contrails/libf/dyn3d_common/infotrac.f90
===================================================================
--- LMDZ6/branches/contrails/libf/dyn3d_common/infotrac.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/dyn3d_common/infotrac.f90	(revision 5609)
@@ -266,6 +266,12 @@
    !---------------------------------------------------------------------------------------------------------------------------
    nqtrue = SIZE(tracers)                                                                               !--- "true" tracers
-   nqo    =      COUNT(tracers(:)%component == 'lmdz' .AND. delPhase(tracers(:)%name)     == 'H2O')     !--- Water phases
-   nbtr = nqtrue-COUNT(tracers(:)%component == 'lmdz' .AND. delPhase(tracers(:)%gen0Name) == 'H2O')     !--- Passed to phytrac
+!   nqo    =      COUNT(tracers(:)%component == 'lmdz' .AND. delPhase(tracers(:)%name)     == 'H2O')     !--- Water phases
+   nqo    =      COUNT(tracers(:)%component == 'lmdz' .AND. &
+       ((delPhase(tracers(:)%name)     == 'H2O') .OR. &    !--- Water phases
+        (delPhase(tracers(:)%name)     == 'CLDFRA')))
+!   nbtr = nqtrue-COUNT(tracers(:)%component == 'lmdz' .AND. delPhase(tracers(:)%gen0Name) == 'H2O')     !--- Passed to phytrac
+   nbtr = nqtrue-COUNT(tracers(:)%component == 'lmdz' .AND. &
+       ((delPhase(tracers(:)%gen0Name)     == 'H2O') .OR. &    !--- Passed to phytrac
+        (delPhase(tracers(:)%gen0Name)     == 'CLDFRA')))
    nqCO2  =      COUNT( [type_trac == 'inco', type_trac == 'co2i'] )
 IF (CPPKEY_INCA) THEN
@@ -335,5 +341,7 @@
       t1%iadv       = iad
       t1%isAdvected = iad >= 0
-      t1%isInPhysics= delPhase(t1%gen0Name) /= 'H2O' .OR. t1%component /= 'lmdz' !=== MORE EXCEPTIONS ? CO2i, SURSAT CLOUD H2O
+!      t1%isInPhysics= delPhase(t1%gen0Name) /= 'H2O' .OR. t1%component /= 'lmdz' !=== MORE EXCEPTIONS ? CO2i, SURSAT CLOUD H2O
+      t1%isInPhysics=((delPhase(t1%gen0Name) /= 'H2O') .AND. &
+                      (delPhase(t1%gen0Name) /= 'CLDFRA')) .OR. t1%component /= 'lmdz'
       ttr(iq)       = t1
 
@@ -388,7 +396,13 @@
 
    !--- Note: nqtottr can differ from nbtr when nmom/=0
-   nqtottr = nqtot - COUNT(delPhase(tracers%gen0Name) == 'H2O' .AND. tracers%component == 'lmdz')
-   IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers%name) == 'H2O' .AND. tracers%component == 'lmdz') /= nqtottr) &
-      CALL abort_gcm(modname, 'problem with the computation of nqtottr', 1)
+!   nqtottr = nqtot - COUNT(delPhase(tracers%gen0Name) == 'H2O' .AND. tracers%component == 'lmdz')
+   nqtottr = nqtot - COUNT(tracers(:)%component == 'lmdz' .AND. &
+       ((delPhase(tracers(:)%gen0Name)     == 'H2O') .OR. & 
+        (delPhase(tracers(:)%gen0Name)     == 'CLDFRA')))
+!   IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers%name) == 'H2O' .AND. tracers%component == 'lmdz') /= nqtottr) &
+!   IF(COUNT(tracers%iso_iName == 0) - COUNT(tracers(:)%component == 'lmdz' .AND. &
+!       ((delPhase(tracers(:)%name)     == 'H2O') .OR. &
+!        (delPhase(tracers(:)%name)     == 'CLDFRA') /= nqtottr) &
+!      CALL abort_gcm(modname, 'problem with the computation of nqtottr', 1)
 
    !=== DISPLAY THE RESULTS
Index: LMDZ6/branches/contrails/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.f90
===================================================================
--- LMDZ6/branches/contrails/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.f90	(revision 5609)
@@ -44,5 +44,5 @@
           solsw, solswfdiff, radsol, t_ancien, wake_deltat, wake_s,  rain_fall, qsol, z0h, &
           sollw,sollwdown, rugoro, q_ancien, wake_deltaq, wake_pe, snow_fall, ratqs,w01, &
-    sig1, ftsol, clwcon, fm_therm, wake_Cstar,  pctsrf,  entr_therm,radpas, f0,&
+    sig1, ftsol, cwcon, clwcon, fm_therm, wake_Cstar,  pctsrf,  entr_therm,radpas, f0,&
     zmax0,fevap, rnebcon,falb_dir, falb_dif, wake_fip,    agesno,  detr_therm, pbl_tke,  &
     phys_state_var_init, ql_ancien, qs_ancien, prlw_ancien, prsw_ancien, &
@@ -50,5 +50,5 @@
     ale_bl, ale_bl_trig, alp_bl, &
     ale_wake, ale_bl_stat, AWAKE_S, &
-    cf_ancien, qvc_ancien, cfa_ancien, qta_ancien
+    cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien
 
   USE comconst_mod, ONLY: pi, dtvr
@@ -243,6 +243,9 @@
   cf_ancien = 0.
   qvc_ancien = 0.
+  cwcon = 0.
   cfa_ancien = 0.
-  qta_ancien = 0.
+  pcf_ancien = 0.
+  qva_ancien = 0.
+  qia_ancien = 0.
   
   z0m(:,:)=0 ! ym missing 5th subsurface initialization
Index: LMDZ6/branches/contrails/libf/misc/readTracFiles_mod.f90
===================================================================
--- LMDZ6/branches/contrails/libf/misc/readTracFiles_mod.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/misc/readTracFiles_mod.f90	(revision 5609)
@@ -63,5 +63,5 @@
     LOGICAL               :: isAdvected  = .FALSE.              !--- "true" tracers: iadv > 0.   COUNT(isAdvected )=nqtrue
     LOGICAL               :: isInPhysics = .TRUE.               !--- "true" tracers: in tr_seri. COUNT(isInPhysics)=nqtottr
-    INTEGER               :: iso_iGroup  = 0                    !--- Isotopes group index in isotopes(:)
+    INTEGER               :: iso_iGroup  = -1                   !--- Isotopes group index in isotopes(:)
     INTEGER               :: iso_iName   = 0                    !--- Isotope  name  index in isotopes(iso_iGroup)%trac(:)
     INTEGER               :: iso_iZone   = 0                    !--- Isotope  zone  index in isotopes(iso_iGroup)%zone(:)
@@ -122,9 +122,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   = 'vlibfcaq'    !--- Old phases for water (no separator)
-  CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsbfcaq'    !--- Known phases initials
+  CHARACTER(LEN=maxlen), PARAMETER :: old_phases   = 'vlib'     !--- Old phases for water (no separator)
+  CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsb'     !--- 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    ','blownSnow', 'fracCloud', 'cldVapor ', 'aviContra', 'qtContra ']
+                                = ['gaseous  ', 'liquid   ', 'solid    ','blownSnow']
   CHARACTER(LEN=1),      SAVE :: phases_sep  =  '_'             !--- Phase separator
   CHARACTER(LEN=maxlen), SAVE :: isoFile = 'isotopes_params.def'!--- Name of the isotopes parameters file
Index: LMDZ6/branches/contrails/libf/phylmd/create_etat0_unstruct_mod.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/create_etat0_unstruct_mod.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/create_etat0_unstruct_mod.f90	(revision 5609)
@@ -260,4 +260,5 @@
     cf_ancien = 0.
     qvc_ancien = 0.
+    cwcon = 0.
 
     wake_delta_pbl_TKE(:,:,:)=0
Index: LMDZ6/branches/contrails/libf/phylmd/dyn1d/old_lmdz1d.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/dyn1d/old_lmdz1d.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/dyn1d/old_lmdz1d.f90	(revision 5609)
@@ -10,5 +10,5 @@
       USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar,getin
    USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, &
-       clwcon, detr_therm, &
+       cwcon, clwcon, detr_therm, &
        qsol, fevap, z0m, z0h, agesno, &
        du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
@@ -23,5 +23,5 @@
        zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, &
        ql_ancien, qs_ancien, qbs_ancien, &
-       cf_ancien, qvc_ancien, cfa_ancien, qta_ancien, &
+       cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien, &
        prlw_ancien, prsw_ancien, prbsw_ancien, prw_ancien, &
        u10m,v10m,ale_wake,ale_bl_stat
@@ -874,8 +874,11 @@
           cf_ancien = 0.
           qvc_ancien = 0.
+          cwcon = 0.
         ENDIF
         IF ( ok_plane_contrail ) THEN
           cfa_ancien = 0.
-          qta_ancien = 0.
+          pcf_ancien = 0.
+          qva_ancien = 0.
+          qia_ancien = 0.
         ENDIF
 !jyg<
Index: LMDZ6/branches/contrails/libf/phylmd/dyn1d/scm.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/dyn1d/scm.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/dyn1d/scm.f90	(revision 5609)
@@ -6,5 +6,5 @@
    USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar,getin
    USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, &
-       clwcon, detr_therm, &
+       cwcon, clwcon, detr_therm, &
        qsol, fevap, z0m, z0h, agesno, &
        du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
@@ -19,5 +19,5 @@
        zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, &
        ql_ancien, qs_ancien, qbs_ancien, &
-       cf_ancien, qvc_ancien, cfa_ancien, qta_ancien, &
+       cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien, &
        prlw_ancien, prsw_ancien, prbsw_ancien, prw_ancien, &
        u10m,v10m,ale_wake,ale_bl_stat, ratqs_inter_
@@ -616,8 +616,11 @@
           cf_ancien = 0.
           qvc_ancien = 0.
+          cwcon = 0.
         ENDIF
         IF ( ok_plane_contrail ) THEN
           cfa_ancien = 0.
-          qta_ancien = 0.
+          pcf_ancien = 0.
+          qva_ancien = 0.
+          qia_ancien = 0.
         ENDIF
         rain_fall=0.
@@ -820,4 +823,7 @@
       ENDIF
       
+      ! calculation of potential temperature for the advection
+      teta=temp*(pzero/play)**rkappa
+
       ! vertical tendencies computed as d X / d t = -W  d X / d z
 
@@ -829,5 +835,5 @@
         d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1))
         ! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
-        d_t_vert_adv(l)=-w_adv(l)*(temp(l+1)-temp(l-1))/(z_adv(l+1)-z_adv(l-1)) 
+        d_t_vert_adv(l)=-w_adv(l)*(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)) / (pzero/play(l))**rkappa
         d_q_vert_adv(l,:)=-w_adv(l)*(q(l+1,:)-q(l-1,:))/(z_adv(l+1)-z_adv(l-1))
       enddo
@@ -842,10 +848,12 @@
           d_u_vert_adv(l)=-w_adv(l)*(u(l)-u(l-1))/(z_adv(l)-z_adv(l-1))
           d_v_vert_adv(l)=-w_adv(l)*(v(l)-v(l-1))/(z_adv(l)-z_adv(l-1))
-          d_t_vert_adv(l)=-w_adv(l)*(temp(l)-temp(l-1))/(z_adv(l)-z_adv(l-1))
+          ! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
+          d_t_vert_adv(l)=-w_adv(l)*(teta(l)-teta(l-1))/(z_adv(l)-z_adv(l-1)) / (pzero/play(l))**rkappa
           d_q_vert_adv(l,:)=-w_adv(l)*(q(l,:)-q(l-1,:))/(z_adv(l)-z_adv(l-1))
         ELSE
           d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l))/(z_adv(l+1)-z_adv(l))
           d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l))/(z_adv(l+1)-z_adv(l))
-          d_t_vert_adv(l)=-w_adv(l)*(temp(l+1)-temp(l))/(z_adv(l+1)-z_adv(l))
+          ! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
+          d_t_vert_adv(l)=-w_adv(l)*(teta(l+1)-teta(l))/(z_adv(l+1)-z_adv(l)) / (pzero/play(l))**rkappa
           d_q_vert_adv(l,:)=-w_adv(l)*(q(l+1,:)-q(l,:))/(z_adv(l+1)-z_adv(l))
         ENDIF
Index: LMDZ6/branches/contrails/libf/phylmd/infotrac_phy.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/infotrac_phy.F90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/infotrac_phy.F90	(revision 5609)
@@ -293,6 +293,10 @@
    !---------------------------------------------------------------------------------------------------------------------------
    nqtrue = SIZE(tracers)                                                                               !--- "true" tracers
-   nqo    =      COUNT(tracers(:)%component == 'lmdz' .AND. delPhase(tracers(:)%name)     == 'H2O')     !--- Water phases
-   nbtr = nqtrue-COUNT(tracers(:)%component == 'lmdz' .AND. delPhase(tracers(:)%gen0Name) == 'H2O')     !--- Passed to phytrac
+   nqo    =      COUNT(tracers(:)%component == 'lmdz' .AND. &
+       ((delPhase(tracers(:)%name)     == 'H2O') .OR. &    !--- Water phases
+        (delPhase(tracers(:)%name)     == 'CLDFRA')))
+   nbtr = nqtrue-COUNT(tracers(:)%component == 'lmdz' .AND. &
+       ((delPhase(tracers(:)%gen0Name)     == 'H2O') .OR. &    !--- Passed to phytrac
+        (delPhase(tracers(:)%gen0Name)     == 'CLDFRA')))
    nqCO2  =      COUNT( [type_trac == 'inco', type_trac == 'co2i'] )
 IF (CPPKEY_INCA) THEN
@@ -350,5 +354,7 @@
       t1%longName   = t1%name; IF(iad > 0) t1%longName=TRIM(t1%name)//descrq(iad)
       t1%isAdvected = iad >= 0
-      t1%isInPhysics= delPhase(t1%gen0Name) /= 'H2O' .OR. t1%component /= 'lmdz' !=== MORE EXCEPTIONS ? CO2i, SURSAT CLOUD H2O
+      !t1%isInPhysics= delPhase(t1%gen0Name) /= 'H2O' .OR. t1%component /= 'lmdz' !=== MORE EXCEPTIONS ? CO2i, SURSAT CLOUD H2O
+      t1%isInPhysics=((delPhase(t1%gen0Name) /= 'H2O') .AND. &
+                      (delPhase(t1%gen0Name) /= 'CLDFRA')) .OR. t1%component /= 'lmdz'
       ttr(iq)       = t1
 
@@ -397,7 +403,13 @@
 
    !--- Note: nqtottr can differ from nbtr when nmom/=0
-   nqtottr = nqtot - COUNT(delPhase(tracers%gen0Name) == 'H2O' .AND. tracers%component == 'lmdz')
-   IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers%name) == 'H2O' .AND. tracers%component == 'lmdz') /= nqtottr) &
-      CALL abort_physic(modname, 'problem with the computation of nqtottr', 1)
+!   nqtottr = nqtot - COUNT(delPhase(tracers%gen0Name) == 'H2O' .AND. tracers%component == 'lmdz')
+   nqtottr = nqtot - COUNT(tracers(:)%component == 'lmdz' .AND. &
+       ((delPhase(tracers(:)%gen0Name)     == 'H2O') .OR. &
+        (delPhase(tracers(:)%gen0Name)     == 'CLDFRA')))
+!   IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers%name) == 'H2O' .AND. tracers%component == 'lmdz') /= nqtottr) &
+!   IF(COUNT(tracers%iso_iName == 0) - COUNT(tracers(:)%component == 'lmdz' .AND. &
+!       ((delPhase(tracers(:)%name)     == 'H2O') .OR. &
+!        (delPhase(tracers(:)%name)     == 'CLDFRA'))) /= nqtottr) &
+!      CALL abort_physic(modname, 'problem with the computation of nqtottr', 1)
 
    !=== DISPLAY THE RESULTS
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90	(revision 5609)
@@ -73,5 +73,5 @@
 SUBROUTINE contrails_formation( &
       klon, dtime, pplay, temp, qsat, qsatl, gamma_cond, flight_dist, &
-      clrfra, pdf_loc, pdf_scale, pdf_alpha, keepgoing, &
+      clrfra, qclr, pdf_scale, pdf_alpha, pdf_gamma, keepgoing, pt_pron_clds, &
       Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
       dcfa_ini, dqia_ini, dqta_ini)
@@ -97,8 +97,10 @@
 REAL,     INTENT(IN) , DIMENSION(klon) :: flight_dist  ! aviation distance flown concentration [m/s/m3]
 REAL,     INTENT(IN) , DIMENSION(klon) :: clrfra       ! clear sky fraction [-]
-REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_loc      ! location parameter of the clear sky PDF [%]
+REAL,     INTENT(IN) , DIMENSION(klon) :: qclr         ! clear sky specific humidity [kg/kg]
 REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_scale    ! scale parameter of the clear sky PDF [%]
 REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_alpha    ! shape parameter of the clear sky PDF [-]
+REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_gamma    ! tmp parameter of the clear sky PDF [-]
 LOGICAL,  INTENT(IN) , DIMENSION(klon) :: keepgoing    ! .TRUE. if a new condensation loop should be computed
+LOGICAL,  INTENT(IN) , DIMENSION(klon) :: pt_pron_clds  ! .TRUE. if clouds are prognostic in this mesh
 !
 ! Output
@@ -119,7 +121,8 @@
 REAL, DIMENSION(klon) :: qsatl_crit, dqsatl_crit
 REAL :: Gcont, psatl_crit, pcrit
-REAL :: pdf_x, pdf_y, pdf_e2, pdf_e3
-REAL :: pdf_fra_above_qcritcont, pdf_fra_above_qsat, pdf_fra_above_qnuc
-REAL :: pdf_q_above_qcritcont, pdf_q_above_qsat, pdf_q_above_qnuc
+REAL :: rhl_clr, pdf_loc
+REAL :: pdf_x, pdf_y, pdf_e3
+REAL :: pdf_fra_above_qcritcont, pdf_fra_above_qsat
+REAL :: pdf_q_above_qcritcont, pdf_q_above_qsat
 REAL :: qpotcontP
 !
@@ -159,5 +162,5 @@
 
 DO i = 1, klon
-  IF ( keepgoing(i) .AND. ( temp(i) .LE. temp_nowater ) ) THEN
+  IF ( keepgoing(i) .AND. pt_pron_clds(i) .AND. ( temp(i) .LE. temp_nowater ) ) THEN
     
     psatl_crit = qsatl_crit(i) / ( EPS_W + ( 1. - EPS_W ) * qsatl_crit(i) ) * pplay(i)
@@ -168,36 +171,25 @@
     IF ( ( clrfra(i) .GT. eps ) .AND. ( temp(i) .LT. Tcritcont(i) ) ) THEN
     
+      rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
+      pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
       pdf_x = qcritcont(i) / qsatl(i) * 100.
-      pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
-      pdf_e2 = GAMMA(1. + 1. / pdf_alpha(i))
+      pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
       pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
-      pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
+      pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
       pdf_fra_above_qcritcont = EXP( - pdf_y ) * clrfra(i)
       pdf_q_above_qcritcont = ( pdf_e3 * clrfra(i) &
-          + pdf_loc(i) * pdf_fra_above_qcritcont ) * qsatl(i) / 100.
-    
-      pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
-      pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
-      pdf_e2 = GAMMA(1. + 1. / pdf_alpha(i))
+          + pdf_loc * pdf_fra_above_qcritcont ) * qsatl(i) / 100.
+    
+      pdf_x = qsat(i) / qsatl(i) * 100.
+      pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
       pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
-      pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
-      pdf_fra_above_qnuc = EXP( - pdf_y ) * clrfra(i)
-      pdf_q_above_qnuc = ( pdf_e3 * clrfra(i) &
-          + pdf_loc(i) * pdf_fra_above_qnuc ) * qsatl(i) / 100.
-    
-      pdf_x = qsat(i) / qsatl(i) * 100.
-      pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
-      pdf_e2 = GAMMA(1. + 1. / pdf_alpha(i))
-      pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
-      pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
+      pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
       pdf_fra_above_qsat = EXP( - pdf_y ) * clrfra(i)
       pdf_q_above_qsat = ( pdf_e3 * clrfra(i) &
-          + pdf_loc(i) * pdf_fra_above_qsat ) * qsatl(i) / 100.
+          + pdf_loc * pdf_fra_above_qsat ) * qsatl(i) / 100.
     
       potcontfraNP(i) = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat)
-      potcontfraP(i) = MAX(0., MIN(pdf_fra_above_qsat - pdf_fra_above_qnuc, &
-                        pdf_fra_above_qcritcont - pdf_fra_above_qnuc))
-      qpotcontP = MAX(0., MIN(pdf_q_above_qsat - pdf_q_above_qnuc, &
-                      pdf_q_above_qcritcont - pdf_q_above_qnuc))
+      potcontfraP(i) = MIN(pdf_fra_above_qsat, pdf_fra_above_qcritcont)
+      qpotcontP = MIN(pdf_q_above_qsat, pdf_q_above_qcritcont)
 
     ENDIF ! temp .LT. Tcritcont
@@ -219,129 +211,10 @@
 
 !**********************************************************************************
-FUNCTION QVAPINCLD_DEPSUB_CONTRAILS( &
-    qvapincld, qiceincld, temp, qsat, pplay, dtime)
-
-USE lmdz_lscp_ini,        ONLY: RV, RLSTT, RTT, EPS_W
-USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
-USE lmdz_lscp_ini,        ONLY: rm_ice_crystals_contrails
+FUNCTION contrail_cross_section_onera()
+
+USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails
 
 IMPLICIT NONE
 
-! inputs
-REAL :: qvapincld     ! 
-REAL :: qiceincld     ! 
-REAL :: temp          ! 
-REAL :: qsat          ! 
-REAL :: pplay         ! 
-REAL :: dtime         ! time step [s]
-! output
-REAL :: qvapincld_depsub_contrails
-! local
-REAL :: pres_sat, rho, kappa
-REAL :: air_thermal_conduct, water_vapor_diff
-REAL :: rm_ice
-
-!--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_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
-!--from Lohmann et al. (2016), where
-!--alpha is the deposition coefficient [-]
-!--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_v is the specific gas constant for humid air [J/kg/K]
-!--T is the temperature [K]
-!--esi is the saturation pressure w.r.t. ice [Pa]
-!--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]
-!
-!--alpha is a coefficient to take into account the fact that during deposition, a water
-!--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
-!--coherent (with the same structure). It has no impact for sublimation.
-!--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
-!--during deposition, and alpha = 1. during sublimation.
-!--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
-!-- C = capa_cond_cirrus * rm_ice
-!
-!--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
-!-- mi = 4/3 * pi * rm_ice**3 * rho_ice
-!--HYPOTHESIS 2: the ice crystals are monodisperse with the
-!--initial radius rm_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 rm_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 rm_ice_0**3 rho_ice )
-!
-!--The deposition equation then reads:
-!-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice
-!-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
-!--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
-!--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
-!-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
-!--  *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 rm_ice_0**2 rho_ice )
-!--and we have
-!-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
-!-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
-!--where kappa = 1/3*rho_ice/capa_cond_cirrus*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 / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay
-!--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
-air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184
-!--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
-water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4
-kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat &
-      * ( RV * temp / water_vapor_diff / pres_sat &
-        + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) )
-!--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
-
-!--Here, the initial vapor in the cloud is qvapincld, and we compute
-!--the new vapor qvapincld_depsub_contrails
-
-rm_ice = rm_ice_crystals_contrails
-
-IF ( qvapincld .GE. qsat ) THEN
-  !--If the cloud is initially supersaturated
-  !--Exact explicit integration (qvc exact, qice explicit)
-  qvapincld_depsub_contrails = qsat + ( qvapincld - qsat ) &
-                * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
-ELSE
-  !--If the cloud is initially subsaturated
-  !--Exact explicit integration (qvc exact, qice explicit)
-  !--Same but depo_coef_cirrus = 1
-  qvapincld_depsub_contrails = qsat + ( qvapincld - qsat ) &
-                * EXP( - dtime * qiceincld / kappa / rm_ice**2 )
-  IF ( qvapincld_depsub_contrails .GE. ( qvapincld + qiceincld ) ) &
-                qvapincld_depsub_contrails = 0.
-ENDIF ! qvapincld .GT. qsat
-
-END FUNCTION QVAPINCLD_DEPSUB_CONTRAILS
-
-
-!**********************************************************************************
-FUNCTION contrail_cross_section_onera()
-
-USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails
-
-IMPLICIT NONE
-
 !
 ! Output
@@ -355,4 +228,5 @@
 
 END FUNCTION contrail_cross_section_onera
+
 
 SUBROUTINE read_aviation_emissions(klon, klev)
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90	(revision 5609)
@@ -12,6 +12,6 @@
     icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, &
     !--AB contrails
-    contfra, qice_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
-    pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
+    contfra, perscontfra, qice_cont, qice_perscont, pclc_nocont, &
+    pcltau_nocont, pclemi_nocont, pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
     xfiwp_nocont, xfiwc_nocont, reice_nocont)
 
@@ -95,6 +95,8 @@
 
   !--AB for contrails. All these are used / outputed only if ok_plane_contrail=y
-  REAL, INTENT(IN)  :: contfra(klon, klev)       ! contrails fraction [-]
-  REAL, INTENT(IN)  :: qice_cont(klon, klev)     ! contrails condensed water [kg/kg]
+  REAL, INTENT(IN)  :: contfra(klon, klev)       ! linear contrails fraction [-]
+  REAL, INTENT(IN)  :: perscontfra(klon, klev)   ! contrail induced cirrus fraction [-]
+  REAL, INTENT(IN)  :: qice_cont(klon, klev)     ! linear contrails condensed water [kg/kg]
+  REAL, INTENT(IN)  :: qice_perscont(klon, klev) ! contrail induced cirrus condensed water [kg/kg]
   REAL, INTENT(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
   REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
@@ -137,6 +139,6 @@
     icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, &
     !--AB for contrails
-    contfra, qice_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
-    pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
+    contfra, perscontfra, qice_cont, qice_perscont, pclc_nocont, pcltau_nocont, &
+    pclemi_nocont, pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
     xfiwp_nocont, xfiwc_nocont, reice_nocont)
   ELSE
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop.f90	(revision 5609)
@@ -11,6 +11,6 @@
     icefrac_optics, dNovrN, ptconv, rnebcon, ccwcon, &
     !--AB contrails
-    contfra, qice_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
-    pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
+    contfra, perscontfra, qice_cont, qice_perscont, pclc_nocont, pcltau_nocont, &
+    pclemi_nocont, pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
     xfiwp_nocont, xfiwc_nocont, reice_nocont)
 
@@ -119,6 +119,8 @@
 
   !--AB for contrails. All these are used / outputed only if ok_plane_contrail=y
-  REAL, INTENT(IN)  :: contfra(klon, klev)       ! contrails fraction [-]
-  REAL, INTENT(IN)  :: qice_cont(klon, klev)     ! contrails condensed water [kg/kg]
+  REAL, INTENT(IN)  :: contfra(klon, klev)       ! linear contrails fraction [-]
+  REAL, INTENT(IN)  :: perscontfra(klon, klev)   ! contrail induced cirrus fraction [-]
+  REAL, INTENT(IN)  :: qice_cont(klon, klev)     ! linear contrails condensed water [kg/kg]
+  REAL, INTENT(IN)  :: qice_perscont(klon, klev) ! contrail induced cirrus condensed water [kg/kg]
   REAL, INTENT(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
   REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
@@ -171,5 +173,6 @@
   REAL zflwp_var, zfiwp_var
   REAL d_rei_dt
-  REAL pclc_var
+  REAL pclc_cont, pclc_perscont
+  REAL rei_cont, rei_perscont
 
 
@@ -252,6 +255,6 @@
       DO k = 1, klev
         DO i = 1, klon
-          pclc_nocont(i,k) = pclc(i, k) - contfra(i, k)
-          xfiwc_nocont(i, k) = xfiwc(i, k) - qice_cont(i, k)
+          pclc_nocont(i,k) = pclc(i, k) - contfra(i, k) - perscontfra(i, k)
+          xfiwc_nocont(i, k) = xfiwc(i, k) - qice_cont(i, k) - qice_perscont(i, k)
         ENDDO
       ENDDO
@@ -393,12 +396,5 @@
 
           IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
-          IF ( ok_plane_contrail ) THEN
-            !--If contrails are activated, rei is a weighted average between the natural
-            !--rei and the contrails rei, with the weights being the fraction of natural
-            !--vs contrail cirrus in the gridbox
-            !--Beware, re_ice_crystals_contrails is in m, while rei is in microns
-            rei = ( rei * pclc_nocont(i,k) &
-                + re_ice_crystals_contrails * 1.E6 * contfra(i,k) ) / pclc(i,k)
-          ENDIF
+
           pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + &
             zfiwp_var*(3.448E-03+2.431/rei)
@@ -687,5 +683,5 @@
           reice_nocont(i,k) = 1.
           pclc_nocont(i,k) = 0.
-          pclc(i,k) = contfra(i,k)
+          pclc(i,k) = contfra(i,k) + perscontfra(i,k)
           pcltau_nocont(i, k) = 0.
           pclemi_nocont(i, k) = 0.
@@ -695,24 +691,49 @@
 
         IF ( contfra(i,k) .GT. 0.01 * seuil_neb ) THEN
+          pclc_cont = contfra(i,k)
+          rei_cont = re_ice_crystals_contrails * 1.E6
+        ELSE
+          pclc_cont = 0.
+          rei_cont = 1.
+        ENDIF
+
+
+        IF ( perscontfra(i,k) .GT. 0.01 * seuil_neb ) THEN
           !--Everything is the same but with contrails
-          zfiwp_var = 1000.*qice_cont(i, k)/contfra(i, k)*rhodz(i, k)
-          pclc_var = contfra(i,k)
-
-          pcltau_cont(i, k) = zfiwp_var*(3.448E-03+2.431/(re_ice_crystals_contrails*1.E6))
+          pclc_perscont = perscontfra(i,k)
+
+          tc = temp(i, k) - 273.15
+          rei_perscont = d_rei_dt*tc + rei_max
+          IF (tc<=-81.4) rei_perscont = rei_min
+
+        ELSE
+          pclc_perscont = 0.
+          rei_perscont = 1.
+        ENDIF
+
+        IF ( MAX(contfra(i,k), perscontfra(i,k)) .GT. 0.01 * seuil_neb ) THEN
+
+          rei = ( rei_cont * pclc_cont + rei_perscont * pclc_perscont ) &
+              / ( pclc_cont + pclc_perscont )
+          zfiwp_var = 1000.*(xfiwc(i, k)-xfiwc_nocont(i, k))&
+                    / (pclc(i, k)-pclc_nocont(i, k))*rhodz(i, k)
+
+          pcltau_cont(i, k) = zfiwp_var*(3.448E-03+2.431/rei)
           ! -- cloud infrared emissivity:
           ! [the broadband infrared absorption coefficient is PARAMETERized
           ! as a function of the effective cld droplet radius]
           ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
-          k_ice = k_ice0 + 1.0/(re_ice_crystals_contrails*1.E6)
+          k_ice = k_ice0 + 1.0/rei
           pclemi_cont(i, k) = 1.0 - exp(-df*k_ice*zfiwp_var)
 
         ELSE
-          pclc_var = 0.
           pcltau_cont(i, k) = 0.
           pclemi_cont(i, k) = 0.
         ENDIF
 
-        rei = ( re_ice_crystals_contrails*1.E6 * pclc_var &
-            + reice_nocont(i, k) * pclc_nocont(i,k) ) / ( pclc_var + pclc_nocont(i,k) )
+
+        rei = ( rei_cont * pclc_cont + rei_perscont * pclc_perscont &
+            + reice_nocont(i, k) * pclc_nocont(i, k) ) &
+            / ( pclc_cont + pclc_perscont + pclc_nocont(i,k) )
 
         zflwp_var = 1000.*xflwc(i, k)/pclc(i, k)*rhodz(i, k)
@@ -861,8 +882,8 @@
     DO k = klev, 1, -1
       DO i = 1, klon
-        zclear(i) = zclear(i)*(1.-max(contfra(i,k),zcloud(i)))/(1.-min(real( &
-          zcloud(i),kind=8),1.-zepsec))
+        zclear(i) = zclear(i)*(1.-max(contfra(i,k)+perscontfra(i,k),zcloud(i)))/&
+          (1.-min(real(zcloud(i),kind=8),1.-zepsec))
         pct_cont(i) = 1. - zclear(i)
-        zcloud(i) = contfra(i,k)
+        zcloud(i) = contfra(i,k) + perscontfra(i,k)
         IF (paprs(i,k)<prmhc) THEN
           pch_nocont(i) = pch_nocont(i)*(1.-max(pclc_nocont(i,k),zcloudh(i)))/(1.-min(real( &
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90	(revision 5609)
@@ -8,5 +8,6 @@
 SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
      paprs, pplay, omega, temp, qt, ql_seri, qi_seri,   &
-     ptconv, cldfracv, qcondcv, ratqs, sigma_qtherm,    &
+     ratqs, sigma_qtherm, ptconv, cfcon_old, qvcon_old, &
+     qccon_old, cfcon, qvcon, qccon,                    &
      d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
      pfraclr, pfracld,                                  &
@@ -24,6 +25,6 @@
      dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
      dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
-     cfa_seri, qta_seri, flight_dist, flight_h2o,       &
-     qice_cont, Tcritcont,                              &
+     cfa_seri, pcf_seri, qva_seri, qia_seri,flight_dist,&
+     flight_h2o, qice_perscont, Tcritcont,              &
      qcritcont, potcontfraP, potcontfraNP,              &
      cloudth_sth,                                       &
@@ -117,9 +118,10 @@
 USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca
 USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat
-USE lmdz_lscp_ini, ONLY : t_glace_min, min_frac_th_cld
+USE lmdz_lscp_ini, ONLY : t_glace_min, temp_nowater, min_frac_th_cld
 USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG
 USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp
 USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
-USE lmdz_lscp_ini, ONLY : ok_plane_contrail
+USE lmdz_lscp_ini, ONLY : ok_weibull_warm_clouds, ok_no_issr_strato
+USE lmdz_lscp_ini, ONLY : ok_plane_contrail, ok_ice_sedim
 
 ! Temporary call for Lamquin et al (2012) diagnostics
@@ -153,6 +155,10 @@
                                                                    ! CR: if iflag_ice_thermo=2, only convection is active
   LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cldfracv        ! cloud fraction from deep convection [-]
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qcondcv         ! in-cloud condensed specific humidity from deep convection [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cfcon_old       ! cloud fraction from deep convection from previous timestep [-]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qvcon_old       ! in-cloud vapor specific humidity from deep convection from previous timestep [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qccon_old       ! in-cloud condensed specific humidity from deep convection from previous timestep [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cfcon           ! cloud fraction from deep convection [-]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qvcon           ! in-cloud vapor specific humidity from deep convection [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qccon           ! in-cloud condensed specific humidity from deep convection [kg/kg]
 
   !Inputs associated with thermal plumes
@@ -185,5 +191,7 @@
   !--------------------------------------------------
   REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cfa_seri         ! linear contrails fraction [-]
-  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qta_seri         ! linear contrails total specific humidity [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: pcf_seri         ! contrails induced cirrus fraction [-]
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qva_seri         ! linear contrails total specific humidity [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qia_seri         ! linear contrails total specific humidity [kg/kg]
   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 [kgH2O/s/mesh]
@@ -243,5 +251,5 @@
   ! for contrails and aviation
 
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_cont      !--condensed water in linear contrails used in the radiation scheme [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_perscont  !--condensed water in contrail induced cirrus used in the radiation scheme [kg/kg]
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcritcont      !--critical temperature for contrail formation [K]
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcritcont      !--critical specific humidity for contrail formation [kg/kg]
@@ -279,5 +287,5 @@
   ! LOCAL VARIABLES:
   !----------------
-  REAL, DIMENSION(klon) :: qliq_in, qice_in
+  REAL, DIMENSION(klon) :: qliq_in, qice_in, qvc_in, cldfra_in
   REAL, DIMENSION(klon,klev) :: ctot
   REAL, DIMENSION(klon,klev) :: ctot_vol
@@ -307,5 +315,5 @@
   REAL, DIMENSION(klon) :: ztupnew
   REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
-  REAL, DIMENSION(klon) :: zqice_sedim ! for sedimentation of ice crystals
+  REAL, DIMENSION(klon) :: cldfra_above, icesed_flux ! for sedimentation of ice crystals
   REAL, DIMENSION(klon) :: zrflclr, zrflcld
   REAL, DIMENSION(klon) :: ziflclr, ziflcld
@@ -324,5 +332,6 @@
   REAL :: delta_z
   ! for contrails
-  REAL, DIMENSION(klon) :: contfra, qcont
+  REAL, DIMENSION(klon) :: contfra, perscontfra, qcont
+  LOGICAL, DIMENSION(klon) :: pt_pron_clds
   !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrail)
   ! Constants used for calculating ratios that are advected (using a parent-child
@@ -424,12 +433,15 @@
 dcf_con(:,:)    = 0.
 dcf_mix(:,:)    = 0.
+dcfsed(:,:)     = 0.
 dqi_adj(:,:)    = 0.
 dqi_sub(:,:)    = 0.
 dqi_con(:,:)    = 0.
 dqi_mix(:,:)    = 0.
+dqised(:,:)     = 0.
 dqvc_adj(:,:)   = 0.
 dqvc_sub(:,:)   = 0.
 dqvc_con(:,:)   = 0.
 dqvc_mix(:,:)   = 0.
+dqvcsed(:,:)    = 0.
 qvc(:)          = 0.
 shear(:)        = 0.
@@ -467,5 +479,6 @@
 zqupnew(:)    = 0.
 zqvapclr(:)   = 0.
-zqice_sedim(:)= 0.
+cldfra_above(:)= 0.
+icesed_flux(:)= 0.
 
 
@@ -506,4 +519,8 @@
         !c_iso init of iso
     ENDDO
+    IF ( ok_ice_supersat ) THEN
+      cldfra_in(:) = cf_seri(:,k)
+      qvc_in(:) = qvc_seri(:,k)
+    ENDIF
 
     ! --------------------------------------------------------------------
@@ -519,10 +536,9 @@
       CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
                         zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
-                        zqvapclr, zqupnew, zqice_sedim, &
-                        cf_seri(:,k), qvc_seri(:,k), qliq_in, qice_in, &
+                        zqvapclr, zqupnew, icesed_flux, &
+                        cldfra_in, qvc_in, qliq_in, qice_in, &
                         zrfl, zrflclr, zrflcld, &
                         zifl, ziflclr, ziflcld, &
-                        dqreva(:,k), dqssub(:,k), &
-                        dcfsed(:,k), dqvcsed(:,k) &
+                        dqreva(:,k), dqssub(:,k) &
                         )
 
@@ -531,5 +547,5 @@
 
       CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
-                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, &
+                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, icesed_flux, &
                         zrfl, zrflclr, zrflcld, &
                         zifl, ziflclr, ziflcld, &
@@ -646,4 +662,74 @@
             lognormale(:) = .FALSE.
         
+        ENDIF
+
+
+        IF ( ok_ice_supersat ) THEN
+
+          !--Initialisation
+          IF ( ok_plane_contrail ) THEN
+            contfra(:)     = 0.
+            qcont(:)       = 0.
+            perscontfra(:) = 0.
+          ENDIF
+
+          DO i = 1, klon
+
+            pt_pron_clds(i) = ( ( ( ( zt(i) .LE. temp_nowater ) .OR. ok_weibull_warm_clouds ) &
+                .AND. ( .NOT. ok_no_issr_strato .OR. ( stratomask(i,k) .EQ. 0. ) ) ) &
+                .AND. ( cfcon(i,k) .LT. ( 1. - eps ) ) )
+
+            IF ( pt_pron_clds(i) ) THEN
+
+              !--If deep convection is activated, the condensation scheme activates
+              !--only in the environment. NB. the clear sky fraction will the be
+              !--maximised by 1. - cfcon(i,k)
+              IF ( ptconv(i,k) ) zq(i) = zq(i) - ( qvcon(i,k) + qccon(i,k) ) * cfcon(i,k)
+
+              IF ( ( cfcon(i,k) * qccon(i,k) ) .LT. ( cfcon_old(i,k) * qccon_old(i,k) ) ) THEN
+                !--If deep convection is weakening, we add the clouds that are not anymore
+                !--'in' deep convection to the advected clouds
+                cldfra_in(i) = cldfra_in(i) + MAX(0., cfcon_old(i,k) - cfcon(i,k))
+                qvc_in(i) = qvc_in(i) + qvcon_old(i,k) * MAX(0., cfcon_old(i,k) - cfcon(i,k))
+                qice_in(i) = qice_in(i) + ( qccon_old(i,k) * cfcon_old(i,k) &
+                                          - qccon(i,k) * cfcon(i,k) )
+              ELSEIF ( cfcon(i,k) .GT. cfcon_old(i,k) ) THEN
+                !--Else if deep convection is strengthening, it consumes the existing cloud
+                !--fraction (which does not at this moment represent deep convection)
+                !--NB. if deep convection is strengthening while the fraction decreases,
+                !--clear sky water vapor will be transfered in priority
+                cldfra_in(i) = cldfra_in(i) * ( 1. &
+                             - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
+                qvc_in(i)    = qvc_in(i)    * ( 1. &
+                             - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
+                qice_in(i)   = qice_in(i)   * ( 1. &
+                             - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
+              ENDIF
+
+              !--Barriers
+              cldfra_in(i) = MAX(0., MIN(1. - cfcon(i,k), cldfra_in(i)))
+              qvc_in(i)    = MAX(0., MIN(zq(i), qvc_in(i)))
+              qice_in(i)   = MAX(0., MIN(zq(i) - qvc_in(i), qice_in(i)))
+
+              !--Calculate the shear value (input for condensation and ice supersat)
+              !--Cell thickness [m]
+              delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * zt(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
+            ENDIF
+          ENDDO
         ENDIF
 
@@ -717,24 +803,9 @@
                   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
+                    IF ( iftop ) THEN
+                      cldfra_above(:) = 0.
+                    ELSE
+                      cldfra_above(:) = rneb(:,k+1)
+                    ENDIF
 
                     !---------------------------------------------
@@ -744,15 +815,15 @@
                     CALL condensation_ice_supersat( &
                         klon, dtime, pplay(:,k), paprs(:,k), paprs(:,k+1), &
-                        cldfracv(:,k), qcondcv(:,k), &
-                        cf_seri(:,k), qvc_seri(:,k), qliq_in, qice_in, &
-                        shear, tke_dissip(:,k), cell_area, stratomask(:,k), &
-                        Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, &
+                        cfcon(:,k), cldfra_in, qvc_in, qliq_in, qice_in, &
+                        shear, tke_dissip(:,k), cell_area, Tbef, zq, zqs, &
+                        gammasat, ratqs(:,k), keepgoing, pt_pron_clds, &
+                        cldfra_above, icesed_flux,&
                         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), &
-                        cfa_seri(:,k), qta_seri(:,k), flight_dist(:,k), flight_h2o(:,k), &
-                        contfra, qcont, Tcritcont(:,k), qcritcont(:,k), &
-                        potcontfraP(:,k), potcontfraNP(:,k), &
+                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), dcfsed(:,k), &
+                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), dqised(:,k), &
+                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), dqvcsed(:,k), &
+                        cfa_seri(:,k), pcf_seri(:,k), qva_seri(:,k), qia_seri(:,k), &
+                        flight_dist(:,k), flight_h2o(:,k), contfra, perscontfra, qcont, &
+                        Tcritcont(:,k), qcritcont(:,k), potcontfraP(:,k), potcontfraNP(:,k), &
                         dcfa_ini(:,k), dqia_ini(:,k), dqta_ini(:,k), &
                         dcfa_sub(:,k), dqia_sub(:,k), dqta_sub(:,k), &
@@ -951,9 +1022,4 @@
 
     IF (ok_plane_contrail) THEN
-      !!--Contrails do not precipitate. We remove then from the variables temporarily
-      !DO i = 1, klon
-      !  rneb(i,k) = rneb(i,k) - contfra(i)
-      !  zoliqi(i) = zoliqi(i) - ( qcont(i) - zqs(i) * contfra(i) )
-      !ENDDO
       !--Contrails precipitate as natural clouds. We save the partition of ice
       !--between natural clouds and contrails
@@ -975,5 +1041,5 @@
                             ctot_vol(:,k), ptconv(:,k), &
                             zt, zq, zoliql, zoliqi, zfice, &
-                            rneb(:,k), zqice_sedim, znebprecipclr, znebprecipcld, &
+                            rneb(:,k), icesed_flux, znebprecipclr, znebprecipcld, &
                             zrfl, zrflclr, zrflcld, &
                             zifl, ziflclr, ziflcld, &
@@ -991,10 +1057,10 @@
 
       CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
-                            ctot_vol(:,k), ptconv(:,k), zdqsdT_raw, &
-                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
+                            ctot_vol(:,k), ptconv(:,k), pt_pron_clds, zdqsdT_raw, &
+                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, icesed_flux, &
                             rneb(:,k), znebprecipclr, znebprecipcld, &
                             zneb, tot_zneb, zrho_up, zvelo_up, &
                             zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
-                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) &
+                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k), dqised(:,k) &
                             )
 
@@ -1002,9 +1068,4 @@
     
     IF (ok_plane_contrail) THEN
-      !!--Contrails are reintroduced in the variables
-      !DO i = 1, klon
-      !  rneb(i,k) = rneb(i,k) + contfra(i)
-      !  zoliqi(i) = zoliqi(i) + ( qcont(i) - zqs(i) * contfra(i) )
-      !ENDDO
       !--Contrails fraction is left unchanged, but contrails water has changed
       DO i = 1, klon
@@ -1106,19 +1167,32 @@
     ! P6 > write diagnostics and outputs
     !------------------------------------------------------------
+
+    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)
    
     !--AB Write diagnostics and tracers for ice supersaturation
     IF ( ok_plane_contrail ) THEN
       cfa_seri(:,k) = contfra(:)
-      qta_seri(:,k) = qcont(:)
-      qice_cont(:,k) = qcont(:) - zqs(:) * contfra(:)
+      pcf_seri(:,k) = perscontfra(:)
+      qva_seri(:,k) = zqs(:) * contfra(:)
+      qia_seri(:,k) = qcont(:) - zqs(:) * contfra(:)
+      DO i = 1, klon
+        IF ( ( rneb(i,k) - cfa_seri(i,k) ) .GT. eps ) THEN
+          qice_perscont(i,k) = ( radocond(i,k) - qia_seri(i,k) ) &
+              * perscontfra(i) / ( rneb(i,k) - cfa_seri(i,k) )
+        ELSE
+          qice_perscont(i,k) = 0.
+        ENDIF
+      ENDDO
     ENDIF
 
-    IF ( ok_ice_supersat .AND. .NOT. ok_unadjusted_clouds) qvc(:) = zqs(:) * rneb(:,k)
-
     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
+
+        !--If prognostic clouds are activated, deep convection vapor is
+        !--re-added to the total water vapor
+        IF ( ptconv(i,k) .AND. pt_pron_clds(i) ) &
+            zq(i) = zq(i) + ( qvcon(i,k) + qccon(i,k) ) * cfcon(i,k)
 
         cf_seri(i,k) = rneb(i,k)
@@ -1236,4 +1310,10 @@
   ENDDO
 
+  IF ( ok_ice_sedim ) THEN
+    DO i = 1, klon
+      snow(i) = snow(i) + icesed_flux(i)
+    ENDDO
+  ENDIF
+
   IF (ncoreczq>0) THEN
       WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90	(revision 5609)
@@ -94,11 +94,13 @@
 !**********************************************************************************
 SUBROUTINE condensation_ice_supersat( &
-      klon, dtime, pplay, paprsdn, paprsup, cldfracv, qcondcv, &
-      cldfra_in, qvc_in, qliq_in, qice_in, shear, pbl_eps, cell_area, stratomask, &
-      temp, qtot_in, 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, &
-      contfra_in, qcont_in, flight_dist, flight_h2o, contfra, qcont, &
-      Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
+      klon, dtime, pplay, paprsdn, paprsup, cfcon, &
+      cldfra_in, qvc_in, qliq_in, qice_in, shear, pbl_eps, cell_area, &
+      temp, qtot_in, qsat, gamma_cond, ratqs, keepgoing, pt_pron_clds, &
+      cldfra_above, icesed_flux, &
+      cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, dcf_sed, &
+      dqi_adj, dqi_sub, dqi_con, dqi_mix, dqi_sed, &
+      dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, dqvc_sed, &
+      contfra_in, perscontfra_in, qva_in, qia_in, flight_dist, flight_h2o, &
+      contfra, perscontfra, qcont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
       dcfa_ini, dqia_ini, dqta_ini, dcfa_sub, dqia_sub, dqta_sub, &
       dcfa_cir, dqta_cir, dcfa_mix, dqia_mix, dqta_mix)
@@ -118,10 +120,9 @@
 
 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, ok_no_issr_strato
-
-USE lmdz_lscp_ini,   ONLY: depo_coef_cirrus, capa_cond_cirrus, std_subl_pdf_lscp
-USE lmdz_lscp_ini,   ONLY: beta_pdf_lscp, temp_thresh_pdf_lscp
+USE lmdz_lscp_ini,   ONLY: RLSTT, RTT, RD, RG, RV, RPI, EPS_W
+USE lmdz_lscp_ini,   ONLY: eps, temp_nowater, ok_unadjusted_clouds
+
+USE lmdz_lscp_ini,   ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
+USE lmdz_lscp_ini,   ONLY: mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp
 USE lmdz_lscp_ini,   ONLY: std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp
 USE lmdz_lscp_ini,   ONLY: coef_mixing_lscp, coef_shear_lscp
@@ -143,6 +144,5 @@
 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) :: cldfracv      ! cloud fraction from deep convection [-]
-REAL,     INTENT(IN)   , DIMENSION(klon) :: qcondcv       ! in-cloud condensed specific humidity from deep convection [kg/kg]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: cfcon         ! cloud fraction from deep convection [-]
 REAL,     INTENT(IN)   , DIMENSION(klon) :: cldfra_in     ! cloud fraction [-]
 REAL,     INTENT(IN)   , DIMENSION(klon) :: qvc_in        ! gridbox-mean water vapor in cloud [kg/kg]
@@ -152,5 +152,4 @@
 REAL,     INTENT(IN)   , DIMENSION(klon) :: pbl_eps       ! TKE dissipation [m2/s3]
 REAL,     INTENT(IN)   , DIMENSION(klon) :: cell_area     ! cell area [m2]
-REAL,     INTENT(IN)   , DIMENSION(klon) :: stratomask    ! fraction of stratosphere in the mesh (1 or 0)
 REAL,     INTENT(IN)   , DIMENSION(klon) :: temp          ! temperature [K]
 REAL,     INTENT(IN)   , DIMENSION(klon) :: qtot_in       ! total specific humidity (without precip) [kg/kg]
@@ -159,9 +158,14 @@
 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
+LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: pt_pron_clds  ! .TRUE. if clouds are prognostic in this mesh
+REAL,     INTENT(IN)   , DIMENSION(klon) :: cldfra_above  ! cloud fraction IN THE LAYER ABOVE [-]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: icesed_flux   ! sedimentated ice flux [kg/s/m2]
 !
 ! Input for aviation
 !
 REAL,     INTENT(IN)   , DIMENSION(klon) :: contfra_in    ! input linear contrails fraction [-]
-REAL,     INTENT(IN)   , DIMENSION(klon) :: qcont_in      ! input linear contrails total specific humidity [kg/kg]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: perscontfra_in! input contrail induced cirrus fraction [-]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: qva_in        ! input linear contrails total specific humidity [kg/kg]
+REAL,     INTENT(IN)   , DIMENSION(klon) :: qia_in        ! input linear contrails total specific humidity [kg/kg]
 REAL,     INTENT(IN)   , DIMENSION(klon) :: flight_dist   ! aviation distance flown concentration [m/s/m3]
 REAL,     INTENT(IN)   , DIMENSION(klon) :: flight_h2o    ! aviation emitted H2O concentration [kgH2O/s/m3]
@@ -186,12 +190,15 @@
 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) :: dcf_sed  ! cloud fraction tendency because of sedimentation [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) :: dqi_sed  ! specific ice content tendency because of sedimentation [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]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_sed ! specific cloud water vapor tendency because of sedimentation [kg/kg/s]
 !
 !  Diagnostics for aviation
@@ -199,4 +206,5 @@
 !
 REAL,     INTENT(INOUT), DIMENSION(klon) :: contfra      ! linear contrail fraction [-]
+REAL,     INTENT(INOUT), DIMENSION(klon) :: perscontfra  ! linear contrail induced cirrus fraction [-]
 REAL,     INTENT(INOUT), DIMENSION(klon) :: qcont        ! linear contrail specific humidity [kg/kg]
 REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcritcont    ! critical temperature for contrail formation [K]
@@ -220,5 +228,5 @@
 INTEGER :: i
 LOGICAL :: ok_warm_cloud
-REAL, DIMENSION(klon) :: qtot, qcld, qzero
+REAL, DIMENSION(klon) :: qcld, qzero
 REAL, DIMENSION(klon) :: clrfra, qclr
 !
@@ -228,32 +236,43 @@
 !
 ! for unadjusted clouds
-REAL :: qvapincld, qiceincld, qvapincld_new
-!
-! for sublimation
-REAL :: dqt_sub
+REAL :: qiceincld, qvapincld, qvapincld_new
+!
+! for deposition / sublimation
+REAL :: pres_sat, kappa_depsub, tau_depsub
+REAL :: air_thermal_conduct, water_vapor_diff
+REAL :: N_ice
+!
+! for dissipation
+REAL :: pdf_shape
 !
 ! for condensation
 REAL, DIMENSION(klon) :: qsatl, dqsatl
-REAL, DIMENSION(klon) :: pdf_alpha, pdf_scale, pdf_loc
-REAL :: rhl_clr
-REAL :: pdf_ratqs, pdf_skew
-REAL :: pdf_x, pdf_y, pdf_T
-REAL :: pdf_e3, pdf_e4
-REAL :: cf_cond, qt_cond, dqt_con
+REAL, DIMENSION(klon) :: pdf_alpha, pdf_scale, pdf_gamma
+REAL :: rhl_clr, pdf_loc
+REAL :: pdf_e3, pdf_x, pdf_y
+REAL :: dqt_con
+!
+! for sedimentation
+REAL, DIMENSION(klon) :: qice_sedim
+REAL :: clrfra_sed
 !
 ! for mixing
 REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix
-REAL :: clrfra_mix, sigma_mix
+REAL :: cldfra_mix, clrfra_mix, sigma_mix
 REAL :: L_shear, shear_fra
-REAL :: qvapinclr_lim
+REAL :: qiceinmix, qvapinmix_lim, qvapinclr_lim
 REAL :: pdf_fra_above_nuc, pdf_q_above_nuc
 REAL :: pdf_fra_above_lim, pdf_q_above_lim
-REAL :: pdf_fra_below_lim
+REAL :: pdf_fra_below_lim, pdf_q_below_lim
+REAL :: mixed_fraction
 !
 ! for cell properties
-REAL, DIMENSION(klon) :: dz
+REAL :: rho, rhodz, dz
+!
+! for contrails
+REAL :: perscontfra_ratio
+REAL :: contrails_conversion_factor
 
 qzero(:) = 0.
-qtot = qtot_in
 
 !--Calculation of qsat w.r.t. liquid
@@ -269,10 +288,4 @@
   !--formed elsewhere)
   IF (keepgoing(i)) THEN
-
-    !--Initialisation
-    issrfra(i)  = 0.
-    qissr(i)    = 0.
-    contfra(i)  = 0.
-    qcont(i)    = 0.
 
     !--If the temperature is higher than the threshold below which
@@ -281,10 +294,9 @@
     !--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 ) .OR. &
-         ( ok_no_issr_strato .AND. ( stratomask(i) .EQ. 1. ) ) ) 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) ) )
+    IF ( .NOT. pt_pron_clds(i) ) THEN
+
+      pdf_std   = ratqs(i) * qtot_in(i)
+      pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot_in(i))**2 ) )
+      pdf_delta = LOG( qtot_in(i) / ( gamma_cond(i) * qsat(i) ) )
       pdf_a     = pdf_delta / ( pdf_k * SQRT(2.) )
       pdf_b     = pdf_k / (2. * SQRT(2.))
@@ -303,5 +315,5 @@
       ELSE
         cldfra(i) = 0.5 * pdf_e1
-        qincld(i) = qtot(i) * pdf_e2 / pdf_e1
+        qincld(i) = qtot_in(i) * pdf_e2 / pdf_e1
         qvc(i) = qsat(i) * cldfra(i)
       ENDIF
@@ -319,17 +331,14 @@
       ok_warm_cloud = ( temp(i) .GT. temp_nowater )
 
-      !--We remove the convection water from the total available water
-      qtot(i) = qtot(i) - ( qcondcv(i) + qsat(i) ) * cldfracv(i)
-
       !--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. - cldfracv(i), cldfra_in(i)))
-      qcld(i)   = MAX(0., MIN(qtot(i), qliq_in(i) + qice_in(i) + qvc_in(i)))
+      cldfra(i) = MAX(0., MIN(1. - cfcon(i), cldfra_in(i)))
+      qcld(i)   = MAX(0., MIN(qtot_in(i), qliq_in(i) + qice_in(i) + qvc_in(i)))
       qvc(i)    = MAX(0., MIN(qcld(i), qvc_in(i)))
 
       !--Initialise clear fraction properties
-      clrfra(i) = ( 1. - cldfracv(i) ) - cldfra(i)
-      qclr(i) = qtot(i) - qcld(i)
+      clrfra(i) = MAX(0., MIN(1., ( 1. - cfcon(i) ) - cldfra(i)))
+      qclr(i) = qtot_in(i) - qcld(i)
 
       dcf_sub(i)  = 0.
@@ -344,18 +353,105 @@
       dqi_mix(i)  = 0.
       dqvc_mix(i) = 0.
+      dcf_sed(i)  = 0.
+      dqi_sed(i)  = 0.
+      dqvc_sed(i) = 0.
+
+      IF ( icesed_flux(i) .GT. 0. ) THEN
+        !--If ice sedimentation is activated, the quantity of sedimentated ice was added
+        !--to the total water vapor in the precipitation routine. Here we remove it
+        !--(it will be reincluded later)
+        qice_sedim(i) = icesed_flux(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
+        qclr(i) = qclr(i) - qice_sedim(i)
+      ENDIF
 
       !--Initialisation of the cell properties
       !--Dry density [kg/m3]
-      !rho = pplay(i) / temp(i) / RD
+      rho = pplay(i) / temp(i) / RD
       !--Dry air mass [kg/m2]
-      !rhodz = ( paprsdn(i) - paprsup(i) ) / RG
+      rhodz = ( paprsdn(i) - paprsup(i) ) / RG
       !--Cell thickness [m]
-      !dz = rhodz / rho
-      dz(i) = ( paprsdn(i) - paprsup(i) ) / pplay(i) * temp(i) * RD / RG
+      dz = rhodz / rho
+
+      !--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_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
+      !--from Lohmann et al. (2016), where
+      !--alpha is the deposition coefficient [-]
+      !--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_v is the specific gas constant for humid air [J/kg/K]
+      !--T is the temperature [K]
+      !--esi is the saturation pressure w.r.t. ice [Pa]
+      !--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]
+      !
+      !--alpha is a coefficient to take into account the fact that during deposition, a water
+      !--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
+      !--coherent (with the same structure). It has no impact for sublimation.
+      !--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
+      !--during deposition, and alpha = 1. during sublimation.
+      !--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
+      !-- C = capa_cond_cirrus * rm_ice
+      !
+      !--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
+      !-- mi = 4/3 * pi * rm_ice**3 * rho_ice
+      !--HYPOTHESIS 2: the ice crystals are monodisperse with the
+      !--initial radius rm_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 rm_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 rm_ice_0**3 rho_ice )
+      !
+      !--The deposition equation then reads:
+      !-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice
+      !-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
+      !--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
+      !--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
+      !-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
+      !--  *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 rm_ice_0**2 rho_ice )
+      !--and we have
+      !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
+      !-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
+      !--where kappa = 1/3*rho_ice/capa_cond_cirrus*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
+      !--Median ice crystal concentration [#/m3] in cirrus clouds from Kramer et al (2020)
+      N_ice = 0.003 * 1e6
+      !--NB. the greater kappa_depsub, the lower the efficiency of the
+      !--deposition/sublimation process
+      kappa_depsub = ( 4. / 3. * RPI * N_ice / rho * rho_ice )**(1./3.) &
+                   * qsat(i) / ( 4. * RPI * capa_cond_cirrus * N_ice / rho ) &
+                   * ( RV * temp(i) / water_vapor_diff / pres_sat &
+                     + RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) )
 
       !--If contrails are activated
       IF ( ok_plane_contrail ) THEN
-        contfra(i) = contfra_in(i)
-        qcont(i) = qcont_in(i)
+        contfra(i) = MAX(0., MIN(cldfra(i), contfra_in(i)))
+        perscontfra(i) = MAX(0., MIN(cldfra(i), perscontfra_in(i)))
+        qcont(i) = MAX(0., MIN(qcld(i), qva_in(i) + qia_in(i)))
 
         dcfa_ini(i) = 0.
@@ -370,4 +466,8 @@
         dqia_mix(i) = 0.
         dqta_mix(i) = 0.
+      ELSE
+        contfra(i) = 0.
+        perscontfra(i) = 0.
+        qcont(i) = 0.
       ENDIF
 
@@ -379,9 +479,14 @@
       !--If there is a contrail
       IF ( contfra(i) .GT. eps ) THEN
+        !--We remove contrails from the main class
+        cldfra(i) = cldfra(i) - contfra(i)
+        qcld(i) = qcld(i) - qcont(i)
+        qvc(i) = qvc(i) - qsat(i) * contfra(i)
+
         !--The contrail is always adjusted to saturation
         qiceincld = ( qcont(i) / contfra(i) - qsat(i) )
         
         !--If the ice water content is too low, the cloud is purely sublimated
-        IF ( qiceincld .LT. eps ) THEN
+        IF ( qiceincld .LT. ( 0.005 * qsat(i) ) ) THEN
           dcfa_sub(i) = - contfra(i)
           dqia_sub(i) = - qiceincld * contfra(i)
@@ -389,13 +494,11 @@
           contfra(i) = 0.
           qcont(i) = 0.
-          clrfra(i) = clrfra(i) + dcfa_sub(i)
-          qclr(i) = qclr(i) + dqta_sub(i)
-        ELSE
-          !--We remove contrails from the main class
-          cldfra(i) = cldfra(i) - contfra(i)
-          qcld(i) = qcld(i) - qcont(i)
-          qvc(i) = qvc(i) - qsat(i) * contfra(i)
+          clrfra(i) = MIN(1., clrfra(i) - dcfa_sub(i))
+          qclr(i) = qclr(i) - dqta_sub(i)
         ENDIF   ! qiceincld .LT. eps
       ENDIF ! contfra(i) .GT. eps
+
+      !--If there is a contrail induced cirrus, we save the ratio
+      perscontfra_ratio = perscontfra(i) / MAX(eps, cldfra(i))
 
 
@@ -421,53 +524,44 @@
           qcld(i) = 0.
           qvc(i) = 0.
-          clrfra(i) = clrfra(i) - dcf_sub(i)
+          clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i))
           qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i)
-
-          !--If all the ice has been sublimated, we sublimate
-          !--completely the cloud and do not activate the sublimation
-          !--process
-          qvapincld_new = 0.
 
         !--Else, the cloud is adjusted and sublimated
         ELSE
 
-          !--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 ( qiceincld .LT. ( 0.005 * qsat(i) ) ) THEN
+          !  dcf_sub(i) = ( qiceincld / ( 0.005 * qsat(i) ) - 1. ) * cldfra(i)
+          !  dqvc_sub(i) = qvapincld * dcf_sub(i)
+
+          !  cldfra(i) = cldfra(i) + dcf_sub(i)
+          !  qcld(i) = qcld(i) + dqvc_sub(i)
+          !  qvc(i) = qvc(i) + dqvc_sub(i)
+          !  clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i))
+          !  qclr(i) = qclr(i) - dqvc_sub(i)
+          !ENDIF
 
           IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
-            qvapincld_new = QVAPINCLD_DEPSUB( &
-                qvapincld, qiceincld, temp(i), qsat(i), pplay(i), dtime)
-            IF ( qvapincld_new .EQ. 0. ) THEN
-              !--If all the ice has been sublimated, we sublimate
-              !--completely the cloud and do not activate the dissipation
-              !--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.
-              clrfra(i) = clrfra(i) - dcf_sub(i)
-              qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i)
-            ENDIF
+            IF ( qvapincld .GE. qsat(i) ) THEN
+              !--If the cloud is initially supersaturated
+              !--Exact explicit integration (qvc exact, qice explicit)
+              tau_depsub = 1. / ( depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub )
+            ELSE
+              !--If the cloud is initially subsaturated
+              !--Exact explicit integration (qvc exact, qice explicit)
+              !--Same but depo_coef_cirrus = 1
+              tau_depsub = 1. / ( qiceincld**(1./3.) * kappa_depsub )
+            ENDIF ! qvapincld .GT. qsat
+            qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime / tau_depsub )
+            !--If all the ice is sublimated
+            IF ( qvapincld_new .GE. ( qvapincld + qiceincld ) ) qvapincld_new = 0.
           ELSE
             !--We keep the saturation adjustment hypothesis, and the vapor in the
             !--cloud is set equal to the saturation vapor
-            qvapincld_new = qsat(i)
+            IF ( ( qvapincld + qiceincld ) .GT. qsat(i) ) THEN
+              qvapincld_new = qsat(i)
+            ELSE
+              qvapincld_new = 0.
+            ENDIF
           ENDIF ! ok_unadjusted_clouds
-          
-          !--Adjustment of the IWC 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)))
           
 
@@ -477,29 +571,61 @@
 
           !--If the dissipation process must be activated
-          IF ( qvapincld_new .GT. 0. ) THEN
-            !--Normal distribution around qvapincld + qiceincld with width
-            !--constant in the RHi space
-            pdf_y = ( qvapincld_new - qvapincld - qiceincld ) &
-                  / ( std_subl_pdf_lscp / 100. * qsat(i))
-            pdf_e1 = 0.5 * ( 1. + ERF( pdf_y / SQRT(2.) ) )
-            pdf_e2 = EXP( - pdf_y**2 / 2. ) / SQRT( 2. * RPI )
-
-            dcf_sub(i) = - cldfra(i) * pdf_e1
-            dqt_sub    = - cldfra(i) * ( ( qvapincld + qiceincld ) * pdf_e1 &
-                                       - std_subl_pdf_lscp / 100. * qsat(i) * pdf_e2 )
+          !IF ( cldfra(i) .GT. eps ) THEN
+          IF ( qvapincld_new .GT. qvapincld ) THEN
+            !--Gamma distribution starting at qvapincld
+            pdf_shape = ( mu_subl_pdf_lscp + 1. ) / qiceincld
+            pdf_y = pdf_shape * ( qvapincld_new - qvapincld )
+            pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y )
+            pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y )
 
             !--Tendencies and diagnostics
-            dqvc_sub(i) = dqt_sub
+            dcf_sub(i)  = - cldfra(i) * pdf_e1
+            dqi_sub(i)  = - cldfra(i) * pdf_e2 / pdf_shape
+            dqvc_sub(i) = dcf_sub(i) * qvapincld
 
             !--Add tendencies
-            cldfra(i) = cldfra(i) + dcf_sub(i)
-            qcld(i)   = qcld(i) + dqt_sub
+            cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i))
+            qcld(i)   = qcld(i) + dqvc_sub(i) + dqi_sub(i)
             qvc(i)    = qvc(i) + dqvc_sub(i)
-            clrfra(i) = clrfra(i) - dcf_sub(i)
-            qclr(i)   = qclr(i) - dqt_sub
+            clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i))
+            qclr(i)   = qclr(i) - dqvc_sub(i) - dqi_sub(i)
           ENDIF ! qvapincld_new .GT. 0.
+
+
+          !------------------------------------
+          !--       PHASE ADJUSTMENT        --
+          !------------------------------------
+
+          IF ( qvapincld_new .EQ. 0. ) THEN
+            !--If all the ice has been sublimated, we sublimate
+            !--completely the cloud and do not activate the dissipation
+            !--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.
+            clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i))
+            qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i)
+          ELSE
+            !--Adjustment of the IWC 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
 
         ENDIF   ! qiceincld .LT. eps
       ENDIF ! cldfra(i) .GT. eps
+
+      !--If there is a contrail induced cirrus, we restore it
+      perscontfra(i) = perscontfra_ratio * cldfra(i)
 
 
@@ -519,5 +645,5 @@
         !--Calculation of the properties of the PDF
         !--Parameterization from IAGOS observations
-        !--pdf_e1 and pdf_e2 will be reused below
+        !--pdf_alpha, pdf_scale and pdf_gamma will be reused below
 
         !--Coefficient for standard deviation:
@@ -537,14 +663,15 @@
         pdf_alpha(i) = EXP( rhl_clr / 100. ) * pdf_e3
         
-        IF ( ok_warm_cloud ) THEN
-          !--If the statistical scheme is activated, the standard deviation is adapted
-          !--to depend on the pressure level. It is multiplied by ratqs, so that near the
-          !--surface std is almost 0, and upper than about 450 hPa the std is left untouched
-          pdf_std = pdf_std * ratqs(i)
-        ENDIF
-
-        pdf_e2 = GAMMA(1. + 1. / pdf_alpha(i))
-        pdf_scale(i) = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha(i)) - pdf_e2**2 ))
-        pdf_loc(i) = rhl_clr - pdf_scale(i) * pdf_e2
+        !IF ( ok_warm_cloud ) THEN
+        !  !--If the statistical scheme is activated, the standard deviation is adapted
+        !  !--to depend on the pressure level. It is multiplied by ratqs, so that near the
+        !  !--surface std is almost 0, and upper than about 450 hPa the std is left untouched
+        !  pdf_std = pdf_std * ratqs(i)
+        !ENDIF
+
+        pdf_gamma(i) = GAMMA(1. + 1. / pdf_alpha(i))
+        pdf_scale(i) = MAX(eps, pdf_std / SQRT( &
+                                    GAMMA(1. + 2. / pdf_alpha(i)) - pdf_gamma(i)**2 ))
+        pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
 
         !--Calculation of the newly condensed water and fraction (pronostic)
@@ -553,14 +680,14 @@
 
         pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
-        pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
+        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
         pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
-        pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
-        cf_cond = EXP( - pdf_y )
-        qt_cond = ( pdf_e3 + pdf_loc(i) * cf_cond ) * qsatl(i) / 100.
-
-        IF ( cf_cond .GT. eps ) THEN
-
-          dcf_con(i) = clrfra(i) * cf_cond
-          dqt_con = clrfra(i) * qt_cond
+        pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
+        pdf_fra_above_nuc = EXP( - pdf_y )
+        pdf_q_above_nuc = ( pdf_e3 + pdf_loc * pdf_fra_above_nuc ) * qsatl(i) / 100.
+
+        IF ( pdf_fra_above_nuc .GT. eps ) THEN
+
+          dcf_con(i) = clrfra(i) * pdf_fra_above_nuc
+          dqt_con = clrfra(i) * pdf_q_above_nuc
 
           !--Barriers (should be useless
@@ -574,22 +701,21 @@
             qvapincld = gamma_cond(i) * qsat(i)
             qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
-            qvapincld_new = QVAPINCLD_DEPSUB( &
-                qvapincld, qiceincld, temp(i), qsat(i), pplay(i), dtime / 2.)
-            qvapincld = qvapincld_new
+            tau_depsub = 1. / ( depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub )
+            qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime / 2. / tau_depsub )
           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)
+            qvapincld_new = qsat(i)
           ENDIF
 
           !--Tendency on cloud vapor and diagnostic
-          dqvc_con(i) = qvapincld * dcf_con(i)
+          dqvc_con(i) = qvapincld_new * dcf_con(i)
           dqi_con(i) = dqt_con - dqvc_con(i)
 
           !--Add tendencies
-          cldfra(i)  = cldfra(i) + dcf_con(i)
+          cldfra(i)  = MIN(1., cldfra(i) + dcf_con(i))
           qcld(i)    = qcld(i) + dqt_con
           qvc(i)     = qvc(i) + dqvc_con(i)
-          clrfra(i)  = clrfra(i) - dcf_con(i)
+          clrfra(i)  = MAX(0., clrfra(i) - dcf_con(i))
           qclr(i)    = qclr(i) - dqt_con
 
@@ -597,4 +723,7 @@
       ENDIF ! ( 1. - cldfra(i) ) .GT. eps
 
+
+      !--If there is a contrail induced cirrus, we save the ratio
+      perscontfra_ratio = perscontfra(i) / MAX(eps, cldfra(i))
 
       !--------------------------------------
@@ -658,4 +787,8 @@
         clrfra_mix = N_cld_mix * RPI / cell_area(i) &
                    * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
+        !--The fraction of clear 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) &
+                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
 
 
@@ -666,5 +799,5 @@
         !--semi-major axis (a_mix), on the entire cell heigh dz.
         !--The increase in size is
-        L_shear = coef_shear_lscp * shear(i) * dz(i) * dtime
+        L_shear = coef_shear_lscp * shear(i) * dz * dtime
         !--therefore, the fraction of clear sky mixed is
         !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
@@ -678,9 +811,10 @@
         !--mixed clouds are different.
         clrfra_mix = clrfra_mix + shear_fra
-
+        cldfra_mix = cldfra_mix + shear_fra
 
         !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
 
-        clrfra_mix = MIN( clrfra_mix, clrfra(i) )
+        clrfra_mix = MAX(eps, MIN(clrfra(i), clrfra_mix))
+        cldfra_mix = MAX(eps, MIN(cldfra(i), cldfra_mix))
 
         !--We compute the limit vapor in clear sky where the mixed cloud could not
@@ -690,5 +824,15 @@
         !--diagnostic, and if the cloud size is increased, we add the new vapor to the
         !--cloud's vapor without condensing or sublimating ice crystals
-        qvapinclr_lim = ( ( cldfra(i) + clrfra_mix ) * qsat(i) - qcld(i) ) / clrfra_mix
+        IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
+          qiceinmix = ( qcld(i) - qvc(i) ) / cldfra(i) / ( 1. + clrfra_mix / cldfra_mix )
+          tau_depsub = 1. / ( qiceinmix**(1./3.) * kappa_depsub )
+          qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime / tau_depsub ) )
+          qvapinclr_lim = qvapinmix_lim * ( 1. + cldfra_mix / clrfra_mix ) &
+                        - qvc(i) / cldfra(i) * cldfra_mix / clrfra_mix
+        ELSE
+          !--NB. if tau_depsub = 0, we get the same result as above
+          qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) &
+                        - qcld(i) / cldfra(i) * cldfra_mix / clrfra_mix
+        ENDIF
 
         IF ( qvapinclr_lim .LT. 0. ) THEN
@@ -696,5 +840,4 @@
           dcf_mix(i)  = clrfra_mix
           dqvc_mix(i) = clrfra_mix * qclr(i) / clrfra(i)
-          dqi_mix(i)  = 0.
         ELSE
           !--We then calculate the clear sky part where the humidity is lower than
@@ -704,23 +847,17 @@
           !--because the clear sky fraction can only be reduced by condensation.
           !--Thus the `pdf_xxx` variables are well defined.
-          pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
-          pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
+
+          rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
+          pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
+          pdf_x = qvapinclr_lim / qsatl(i) * 100.
+          pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
           pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
-          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
-          pdf_fra_above_nuc = EXP( - pdf_y ) * clrfra(i)
-          pdf_q_above_nuc = ( pdf_e3 * clrfra(i) &
-                          + pdf_loc(i) * pdf_fra_above_nuc ) * qsatl(i) / 100.
-
-          pdf_x = qvapinclr_lim / qsatl(i) * 100.
-          pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
-          pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
-          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
+          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
           pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
           pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
-                          + pdf_loc(i) * pdf_fra_above_lim ) * qsatl(i) / 100.
+                          + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
 
           pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim
-          pdf_fra_above_lim = pdf_fra_above_lim - pdf_fra_above_nuc
-          pdf_q_above_lim = pdf_q_above_lim - pdf_q_above_nuc
+          pdf_q_below_lim = qclr(i) - pdf_q_above_lim
 
           !--sigma_mix is the ratio of the surroundings of the clouds where mixing
@@ -733,23 +870,15 @@
             dcf_mix(i)  = clrfra_mix * sigma_mix
             dqvc_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim
-            dqi_mix(i)  = 0.
           ENDIF
 
           IF ( pdf_fra_below_lim .GT. eps ) THEN
-            dcf_mix(i)  = dcf_mix(i)  - cldfra(i) * ( 1. - sigma_mix )
-            dqvc_mix(i) = dqvc_mix(i) - cldfra(i) * ( 1. - sigma_mix ) &
-                                                  * qvc(i) / cldfra(i)
-            dqi_mix(i)  = dqi_mix(i)  - cldfra(i) * ( 1. - sigma_mix ) &
-                                                  * ( qcld(i) - qvc(i) ) / cldfra(i)
+            dcf_mix(i)  = dcf_mix(i)  - cldfra_mix * ( 1. - sigma_mix )
+            dqvc_mix(i) = dqvc_mix(i) - cldfra_mix * ( 1. - sigma_mix ) &
+                                      * qvc(i) / cldfra(i)
+            dqi_mix(i)  = dqi_mix(i)  - cldfra_mix * ( 1. - sigma_mix ) &
+                                      * ( qcld(i) - qvc(i) ) / cldfra(i)
           ENDIF
+
         ENDIF
-
-        !--Add tendencies
-        cldfra(i)  = cldfra(i) + dcf_mix(i)
-        qcld(i)    = qcld(i) + dqvc_mix(i) + dqi_mix(i)
-        qvc(i)     = qvc(i) + dqvc_mix(i)
-        clrfra(i)  = clrfra(i) - dcf_mix(i)
-        qclr(i)    = qclr(i) - dqvc_mix(i) - dqi_mix(i)
-      
       ENDIF ! ( cldfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps )
 
@@ -812,4 +941,8 @@
         clrfra_mix = N_cld_mix * RPI / cell_area(i) &
                    * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
+        !--The fraction of clear 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) &
+                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
         
         
@@ -823,5 +956,5 @@
         !--With this, the clouds increase in size along b only, by a factor
         !--L_shear * SQRT(2.) / 2. (to account for the 45 degrees orientation of the wind)
-        L_shear = coef_shear_contrails * shear(i) * dz(i) * dtime
+        L_shear = coef_shear_contrails * shear(i) * dz * dtime
         !--therefore, the fraction of clear sky mixed is
         !-- N_cld_mix * ( a * (b + L_shear * SQRT(2.) / 2.) - a * b ) * RPI / 2. / cell_area
@@ -835,9 +968,11 @@
         !--mixed clouds are different.
         clrfra_mix = clrfra_mix + shear_fra
+        cldfra_mix = cldfra_mix + shear_fra
         
         
         !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
         
-        clrfra_mix = MIN( clrfra_mix, clrfra(i) )
+        clrfra_mix = MAX(eps, MIN(clrfra(i), clrfra_mix))
+        cldfra_mix = MAX(eps, MIN(cldfra(i), cldfra_mix))
         
         !--We compute the limit vapor in clear sky where the mixed cloud could not
@@ -847,5 +982,6 @@
         !--diagnostic, and if the cloud size is increased, we add the new vapor to the
         !--cloud's vapor without condensing or sublimating ice crystals
-        qvapinclr_lim = ( ( contfra(i) + clrfra_mix ) * qsat(i) - qcont(i) ) / clrfra_mix
+        qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) &
+                      - qcont(i) / contfra(i) * cldfra_mix / clrfra_mix
 
         IF ( qvapinclr_lim .LT. 0. ) THEN
@@ -853,5 +989,4 @@
           dcfa_mix(i) = clrfra_mix
           dqta_mix(i) = clrfra_mix * qclr(i) / clrfra(i)
-          dqia_mix(i) = 0.
         ELSE
           !--We then calculate the clear sky part where the humidity is lower than
@@ -861,23 +996,16 @@
           !--because the clear sky fraction can only be reduced by condensation.
           !--Thus the `pdf_xxx` variables are well defined.
-          pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
-          pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
+
+          rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
+          pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
+          pdf_x = qvapinclr_lim / qsatl(i) * 100.
+          pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
           pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
-          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
-          pdf_fra_above_nuc = EXP( - pdf_y ) * clrfra(i)
-          pdf_q_above_nuc = ( pdf_e3 * clrfra(i) &
-                          + pdf_loc(i) * pdf_fra_above_nuc ) * qsatl(i) / 100.
-
-          pdf_x = qvapinclr_lim / qsatl(i) * 100.
-          pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
-          pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
-          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
+          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
           pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
           pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
-                          + pdf_loc(i) * pdf_fra_above_lim ) * qsatl(i) / 100.
+                          + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
 
           pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim
-          pdf_fra_above_lim = pdf_fra_above_lim - pdf_fra_above_nuc
-          pdf_q_above_lim = pdf_q_above_lim - pdf_q_above_nuc
 
           !--sigma_mix is the ratio of the surroundings of the clouds where mixing
@@ -893,14 +1021,27 @@
             dcfa_mix(i) = clrfra_mix * sigma_mix
             dqta_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim
-            dqia_mix(i) = 0.
           ENDIF
 
           IF ( pdf_fra_below_lim .GT. eps ) THEN
-            dcfa_mix(i) = dcfa_mix(i) - contfra(i) * ( 1. - sigma_mix )
-            dqta_mix(i) = dqta_mix(i) - contfra(i) * ( 1. - sigma_mix ) &
-                                      * qcont(i) / contfra(i)
-            dqia_mix(i) = dqia_mix(i) - contfra(i) * ( 1. - sigma_mix ) &
-                                      * ( qcont(i) / contfra(i) - qsat(i) )
+            qvapincld = qcont(i) / contfra(i)
+            qiceincld = qvapincld - qsat(i)
+            dcfa_mix(i) = dcfa_mix(i) - cldfra_mix * ( 1. - sigma_mix )
+            dqta_mix(i) = dqta_mix(i) - cldfra_mix * ( 1. - sigma_mix ) * qvapincld
+            dqia_mix(i) = dqia_mix(i) - cldfra_mix * ( 1. - sigma_mix ) * qiceincld
           ENDIF
+
+        ENDIF
+      ENDIF ! ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps )
+
+      IF ( contfra(i) .GT. eps ) THEN
+        !--We balance the mixing tendencies between the different cloud classes
+        mixed_fraction = dcf_mix(i) + dcfa_mix(i)
+        IF ( mixed_fraction .GT. clrfra(i) ) THEN
+          mixed_fraction = clrfra(i) / mixed_fraction
+          dcf_mix(i)  = dcf_mix(i)  * mixed_fraction
+          dqvc_mix(i) = dqvc_mix(i) * mixed_fraction
+          dqi_mix(i)  = dqi_mix(i)  * mixed_fraction
+          dcfa_mix(i) = dcfa_mix(i) * mixed_fraction
+          dqta_mix(i) = dqta_mix(i) * mixed_fraction
         ENDIF
 
@@ -910,6 +1051,81 @@
         clrfra(i)  = clrfra(i) - dcfa_mix(i)
         qclr(i)    = qclr(i) - dqta_mix(i)
-
-      ENDIF ! ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps )
+      ENDIF
+
+      !--Add tendencies
+      cldfra(i)  = cldfra(i) + dcf_mix(i)
+      qcld(i)    = qcld(i) + dqvc_mix(i) + dqi_mix(i)
+      qvc(i)     = qvc(i) + dqvc_mix(i)
+      clrfra(i)  = clrfra(i) - dcf_mix(i)
+      qclr(i)    = qclr(i) - dqvc_mix(i) - dqi_mix(i)
+
+      !--If there is a contrail induced cirrus, we restore it
+      perscontfra(i) = perscontfra_ratio * cldfra(i)
+
+
+      !---------------------------------------
+      !--         ICE SEDIMENTATION         --
+      !---------------------------------------
+      !
+      !--If ice supersaturation is activated, the cloud properties are prognostic.
+      !--The falling ice is then considered a new cloud in the gridbox.
+      !--BEWARE with this parameterisation, we can create a new cloud with a much
+      !--different ice water content and water vapor content than the existing cloud
+      !--(if it exists). This implies than unphysical fluxes of ice and vapor
+      !--occur between the existing cloud and the newly formed cloud (from sedimentation).
+      !--Note also that currently, the precipitation scheme assume a maximum
+      !--random overlap, meaning that the new formed clouds will not be affected
+      !--by vertical wind shear.
+      !
+      IF ( icesed_flux(i) .GT. 0. ) THEN
+
+        clrfra_sed = MIN(clrfra(i), cldfra_above(i) - cldfra(i))
+
+        IF ( ( clrfra_sed .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN
+
+          qiceincld = qice_sedim(i) / cldfra_above(i)
+          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
+            tau_depsub = 1. / ( qiceincld**(1./3.) * kappa_depsub )
+            qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime / tau_depsub ) )
+          ELSE
+            qvapinclr_lim = qsat(i) - qiceincld
+          ENDIF
+
+          rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
+          pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
+          pdf_x = qvapinclr_lim / qsatl(i) * 100.
+          pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
+          pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
+          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
+          pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
+          pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
+                          + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
+
+          IF ( pdf_fra_above_lim .GT. eps ) THEN
+            dcf_sed(i)  = clrfra_sed * pdf_fra_above_lim / clrfra(i)
+            dqvc_sed(i) = dcf_sed(i) * pdf_q_above_lim / pdf_fra_above_lim
+            !--We include the sedimentated ice:
+            dqi_sed(i)  = qiceincld &    !--We include the sedimentated ice:
+                        * ( dcf_sed(i) & !--the part that falls in the newly formed cloud and
+                          + cldfra(i) )  !--the part that falls in the existing cloud
+                        
+          ENDIF
+
+        ELSE
+
+          dqi_sed(i) = qice_sedim(i)
+
+        ENDIF ! clrfra(i) .GT. eps
+
+        !--Add tendencies
+        cldfra(i)  = MIN(1., cldfra(i) + dcf_sed(i))
+        qcld(i)    = qcld(i) + dqvc_sed(i) + dqi_sed(i)
+        qvc(i)     = qvc(i) + dqvc_sed(i)
+        clrfra(i)  = MAX(0., clrfra(i) - dcf_sed(i))
+        !--We re-include sublimated sedimentated ice into clear sky water vapor
+        qclr(i)    = qclr(i) - dqvc_sed(i) + qice_sedim(i) - dqi_sed(i)
+
+      ENDIF ! qice_sedim(i) .GT. 0.
+
 
       !--We put back contrails in the clouds class
@@ -919,38 +1135,17 @@
 
 
-      !--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
-        
       !--Diagnose ISSRs
       IF ( clrfra(i) .GT. eps ) THEN
-        !--We then calculate the part that is greater than qsat
-        !--and lower than gamma_cond * qsat, which is the ice supersaturated region
-        pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
-        pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
+        rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
+        pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
+        pdf_x = qsat(i) / qsatl(i) * 100.
+        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
         pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
-        pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
-        pdf_fra_above_nuc = EXP( - pdf_y ) * clrfra(i)
-        pdf_q_above_nuc = ( pdf_e3 * clrfra(i) &
-                        + pdf_loc(i) * pdf_fra_above_nuc ) * qsatl(i) / 100.
-
-        pdf_x = qsat(i) / qsatl(i) * 100.
-        pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
-        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
-        pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
+        pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
         issrfra(i) = EXP( - pdf_y ) * clrfra(i)
-        qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc(i) * issrfra(i) ) * qsatl(i) / 100.
-
-        issrfra(i) = issrfra(i) - pdf_fra_above_nuc
-        qissr(i) = qissr(i) - pdf_q_above_nuc
+        qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
+      ELSE
+        issrfra(i) = 0.
+        qissr(i) = 0.
       ENDIF
 
@@ -969,4 +1164,20 @@
       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
+      dcf_sed(i)  = dcf_sed(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
+      dqi_sed(i)  = dqi_sed(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
+      dqvc_sed(i) = dqvc_sed(i) / dtime
+
     ENDIF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds
 
@@ -983,23 +1194,39 @@
   CALL contrails_formation( &
       klon, dtime, pplay, temp, qsat, qsatl, gamma_cond, &
-      flight_dist, clrfra, pdf_loc, pdf_scale, pdf_alpha, keepgoing, &
+      flight_dist, clrfra, qclr, pdf_scale, pdf_alpha, pdf_gamma, &
+      keepgoing, pt_pron_clds, &
       Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
       dcfa_ini, dqia_ini, dqta_ini)
 
   DO i = 1, klon
-    IF ( keepgoing(i) .AND. ( temp(i) .LE. temp_nowater ) ) THEN
+    IF ( keepgoing(i) .AND. pt_pron_clds(i) ) THEN
 
       !--Convert existing contrail fraction into "natural" cirrus cloud fraction
-      dcfa_cir(i) = contfra(i) * ( EXP( - dtime / linear_contrails_lifetime ) - 1. )
-      dqta_cir(i) = qcont(i)   * ( EXP( - dtime / linear_contrails_lifetime ) - 1. )
+      IF ( cldfra(i) .GE. ( 1. - cfcon(i) - eps ) ) THEN
+        contrails_conversion_factor = 1.
+      ELSE
+        contrails_conversion_factor = ( 1. - &
+            !--First, the linear contrails are continuously degraded in induced cirrus
+            EXP( - dtime / linear_contrails_lifetime ) &
+            !--Second, if the cloud fills the entire gridbox, the linear contrails
+            !--cannot exist. The exponent is set so that this only happens for
+            !--very cloudy gridboxes
+            * ( 1. - cldfra(i) / ( 1. - cfcon(i) ) )**0.1 )
+      ENDIF
+      dcfa_cir(i) = - contrails_conversion_factor * contfra(i)
+      dqta_cir(i) = - contrails_conversion_factor * qcont(i)
 
       !--Add tendencies
       issrfra(i) = MAX(0., issrfra(i) - dcfa_ini(i))
       qissr(i)   = MAX(0., qissr(i) - dqta_ini(i))
-      cldfra(i)  = MAX(0., MIN(1. - cldfracv(i), cldfra(i) + dcfa_ini(i)))
-      qcld(i)    = MAX(0., MIN(qtot(i), qcld(i) + dqta_ini(i)))
+      clrfra(i)  = MAX(0., clrfra(i) - dcfa_ini(i))
+      qclr(i)    = MAX(0., qclr(i) - dqta_ini(i))
+
+      cldfra(i)  = MAX(0., MIN(1. - cfcon(i), cldfra(i) + dcfa_ini(i)))
+      qcld(i)    = MAX(0., MIN(qtot_in(i), qcld(i) + dqta_ini(i)))
       qvc(i)     = MAX(0., MIN(qcld(i), qvc(i) + dcfa_ini(i) * qsat(i)))
       contfra(i) = MAX(0., MIN(cldfra(i), contfra(i) + dcfa_cir(i) + dcfa_ini(i)))
       qcont(i)   = MAX(0., MIN(qcld(i), qcont(i) + dqta_cir(i) + dqta_ini(i)))
+      perscontfra(i) = perscontfra(i) - dcfa_cir(i)
 
       !--Diagnostics
@@ -1024,6 +1251,8 @@
         cldfra(i) = 0.
         contfra(i)= 0.
+        perscontfra(i) = 0.
         qcld(i)   = 0.
         qvc(i)    = 0.
+        qcont(i)  = 0.
         qincld(i) = qsat(i)
       ELSE
@@ -1036,4 +1265,6 @@
       ENDIF
 
+      IF ( perscontfra(i) .LT. eps ) perscontfra(i) = 0.
+
     ENDIF ! keepgoing
   ENDDO
@@ -1044,146 +1275,3 @@
 !**********************************************************************************
 
-FUNCTION QVAPINCLD_DEPSUB( &
-    qvapincld, qiceincld, temp, qsat, pplay, dtime)
-
-USE lmdz_lscp_ini, ONLY: RV, RLSTT, RTT, RD, EPS_W
-USE lmdz_lscp_ini, ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
-USE lmdz_lscp_ini, ONLY: eps
-
-IMPLICIT NONE
-
-! inputs
-REAL :: qvapincld     ! 
-REAL :: qiceincld     ! 
-REAL :: temp          ! 
-REAL :: qsat          ! 
-REAL :: pplay         ! 
-REAL :: dtime         ! time step [s]
-! outpu
-REAL :: qvapincld_depsub
-
-!
-! local
-REAL :: pres_sat, rho, kappa
-REAL :: air_thermal_conduct, water_vapor_diff
-REAL :: iwc
-REAL :: iwc_log_inf100, iwc_log_sup100
-REAL :: iwc_inf100, alpha_inf100, coef_inf100
-REAL :: mu_sup100, sigma_sup100, coef_sup100
-REAL :: Dm_ice, rm_ice
-
-!--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_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
-!--from Lohmann et al. (2016), where
-!--alpha is the deposition coefficient [-]
-!--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_v is the specific gas constant for humid air [J/kg/K]
-!--T is the temperature [K]
-!--esi is the saturation pressure w.r.t. ice [Pa]
-!--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]
-!
-!--alpha is a coefficient to take into account the fact that during deposition, a water
-!--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
-!--coherent (with the same structure). It has no impact for sublimation.
-!--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
-!--during deposition, and alpha = 1. during sublimation.
-!--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
-!-- C = capa_cond_cirrus * rm_ice
-!
-!--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
-!-- mi = 4/3 * pi * rm_ice**3 * rho_ice
-!--HYPOTHESIS 2: the ice crystals are monodisperse with the
-!--initial radius rm_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 rm_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 rm_ice_0**3 rho_ice )
-!
-!--The deposition equation then reads:
-!-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice
-!-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
-!--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
-!--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
-!-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
-!--  *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 rm_ice_0**2 rho_ice )
-!--and we have
-!-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
-!-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
-!--where kappa = 1/3*rho_ice/capa_cond_cirrus*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 / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay
-!--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
-air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184
-!--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
-water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4
-kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat &
-      * ( RV * temp / water_vapor_diff / pres_sat &
-        + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) )
-!--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
-
-!--Dry density [kg/m3]
-rho = pplay / temp / RD
-
-!--Here, the initial vapor in the cloud is qvapincld, and we compute
-!--the new vapor qvapincld_new
-
-!--rm_ice formula from McFarquhar and Heymsfield (1997)
-iwc = qiceincld * rho * 1e3
-iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837)
-iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) )
-iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) )
-
-alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100
-coef_inf100 = iwc_inf100 * alpha_inf100**3 / 120.
-
-mu_sup100 = ( 5.2 + 0.0013 * ( temp - RTT ) ) &
-          + ( 0.026 - 1.2E-3 * ( temp - RTT ) ) * iwc_log_sup100
-sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp - RTT ) ) &
-             + ( 0.018 - 2.1E-4 * ( temp - RTT ) ) * iwc_log_sup100
-coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3. * mu_sup100 + 4.5 * sigma_sup100**2 )
-
-Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2 ) &
-       * coef_sup100 ) / ( coef_inf100 + coef_sup100 )
-rm_ice = Dm_ice / 2. * 1.E-6
-
-IF ( qvapincld .GE. qsat ) THEN
-  !--If the cloud is initially supersaturated
-  !--Exact explicit integration (qvc exact, qice explicit)
-  qvapincld_depsub = qsat + ( qvapincld - qsat ) &
-                   * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
-ELSE
-  !--If the cloud is initially subsaturated
-  !--Exact explicit integration (qvc exact, qice explicit)
-  !--Same but depo_coef_cirrus = 1
-  qvapincld_depsub = qsat + ( qvapincld - qsat ) &
-                   * EXP( - dtime * qiceincld / kappa / rm_ice**2 )
-  IF ( qvapincld_depsub .GE. ( qvapincld + qiceincld ) ) &
-                   qvapincld_depsub = 0.
-ENDIF ! qvapincld .GT. qsat
-
-END FUNCTION QVAPINCLD_DEPSUB
-
 END MODULE lmdz_lscp_condensation
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90	(revision 5609)
@@ -162,4 +162,7 @@
   !$OMP THREADPRIVATE(ok_weibull_warm_clouds)
 
+  REAL, SAVE, PROTECTED :: ffallv_issr                       ! tuning coefficient crystal fall velocity, cirrus clouds (with ISSR)
+  !$OMP THREADPRIVATE(ffallv_issr)
+
   REAL, SAVE, PROTECTED :: depo_coef_cirrus=.7               ! [-] deposition coefficient for growth of ice crystals in cirrus clouds
   !$OMP THREADPRIVATE(depo_coef_cirrus)
@@ -168,8 +171,8 @@
   !$OMP THREADPRIVATE(capa_cond_cirrus)
 
-  REAL, SAVE, PROTECTED :: std_subl_pdf_lscp=2.              ! [%] standard deviation of the gaussian distribution of water inside ice clouds
-  !$OMP THREADPRIVATE(std_subl_pdf_lscp)
-  
-  REAL, SAVE, PROTECTED :: beta_pdf_lscp=1.E-3               ! [SI] tuning coefficient for the standard deviation of the PDF of water vapor in the clear sky region
+  REAL, SAVE, PROTECTED :: mu_subl_pdf_lscp=1./3.            ! [-] factor for the ice distribution inside cirrus clouds
+  !$OMP THREADPRIVATE(mu_subl_pdf_lscp)
+  
+  REAL, SAVE, PROTECTED :: beta_pdf_lscp=8.75E-4             ! [SI] tuning coefficient for the standard deviation of the PDF of water vapor in the clear sky region
   !$OMP THREADPRIVATE(beta_pdf_lscp)
   
@@ -192,5 +195,5 @@
   !$OMP THREADPRIVATE(b_homofreez)
 
-  REAL, SAVE, PROTECTED :: delta_hetfreez=1.                 ! [-] value between 0 and 1 to simulate for heterogeneous freezing.
+  REAL, SAVE, PROTECTED :: delta_hetfreez=0.85               ! [-] value between 0 and 1 to simulate for heterogeneous freezing.
   !$OMP THREADPRIVATE(delta_hetfreez)
   
@@ -344,6 +347,6 @@
   !$OMP THREADPRIVATE(ok_ice_sedim)
 
-  REAL, SAVE, PROTECTED :: ice_fallspeed_sedim=0.1          ! Ice fallspeed velocity for sedimentation [m/s]
-  !$OMP THREADPRIVATE(ice_fallspeed_sedim)
+  REAL, SAVE, PROTECTED :: fallice_sedim=1.                 ! Tuning factor for ice fallspeed velocity for sedimentation [-]
+  !$OMP THREADPRIVATE(fallice_sedim)
   !--End of the parameters for poprecip
 
@@ -466,11 +469,13 @@
     CALL getin_p('snow_fallspeed_cld',snow_fallspeed_cld)
     CALL getin_p('ok_ice_sedim',ok_ice_sedim)
-    CALL getin_p('ice_fallspeed_sedim',ice_fallspeed_sedim)
+    CALL getin_p('fallice_sedim',fallice_sedim)
     ! for condensation and ice supersaturation
     CALL getin_p('ok_unadjusted_clouds',ok_unadjusted_clouds)
     CALL getin_p('ok_weibull_warm_clouds',ok_weibull_warm_clouds)
+    ffallv_issr=ffallv_lsc
+    CALL getin_p('ffallv_issr',ffallv_issr)
     CALL getin_p('depo_coef_cirrus',depo_coef_cirrus)
     CALL getin_p('capa_cond_cirrus',capa_cond_cirrus)
-    CALL getin_p('std_subl_pdf_lscp',std_subl_pdf_lscp)
+    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)
@@ -562,12 +567,13 @@
     WRITE(lunout,*) 'lscp_ini, snow_fallspeed_cld:', snow_fallspeed_cld
     WRITE(lunout,*) 'lscp_ini, ok_ice_sedim:', ok_ice_sedim
-    WRITE(lunout,*) 'lscp_ini, ice_fallspeed_sedim:', ice_fallspeed_sedim
+    WRITE(lunout,*) 'lscp_ini, fallice_sedim:', fallice_sedim
     ! for condensation and ice supersaturation
     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, ffallv_issr', ffallv_issr
     WRITE(lunout,*) 'lscp_ini, depo_coef_cirrus:', depo_coef_cirrus
     WRITE(lunout,*) 'lscp_ini, capa_cond_cirrus:', capa_cond_cirrus
-    WRITE(lunout,*) 'lscp_ini, std_subl_pdf_lscp:', std_subl_pdf_lscp
+    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
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_precip.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_precip.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_precip.f90	(revision 5609)
@@ -18,9 +18,9 @@
 SUBROUTINE histprecip_precld( &
            klon, dtime, iftop, paprsdn, paprsup, pplay, zt, ztupnew, zq, &
-           zmqc, zneb, znebprecip, znebprecipclr, &
+           zmqc, zneb, znebprecip, znebprecipclr, icesed_flux, &
            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, dqreva, dqssub &
            )
 
-USE lmdz_lscp_ini, ONLY : iflag_evap_prec
+USE lmdz_lscp_ini, ONLY : iflag_evap_prec, ok_ice_sedim
 USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, ztfondue
 USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
@@ -44,4 +44,5 @@
 REAL,    INTENT(INOUT), DIMENSION(klon) :: zq             !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
 REAL,    INTENT(INOUT), DIMENSION(klon) :: zmqc           !--specific humidity in the precipitation falling from the upper layer [kg/kg]
+REAL,    INTENT(IN),    DIMENSION(klon) :: icesed_flux    !--sedimentated ice flux [kg/s/m2]
 
 REAL,    INTENT(IN),    DIMENSION(klon) :: zneb           !--cloud fraction IN THE LAYER ABOVE [-]
@@ -74,4 +75,5 @@
 
 REAL :: zmelt
+REAL :: qice_sedim
 INTEGER :: i
 
@@ -118,4 +120,24 @@
   ENDDO
 
+ENDIF
+
+!--------------------
+!-- ICE SEDIMENTATION
+!--------------------
+IF ( ok_ice_sedim ) THEN
+  DO i =1, klon
+    !--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
+    !--added to the total water content of the gridbox
+    IF ( icesed_flux(i) .GT. 0. ) THEN
+      qice_sedim = icesed_flux(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
+      !--Vapor is updated after evaporation/sublimation (it is increased)
+      zq(i) = zq(i) + qice_sedim
+      !--Air and precip temperature (i.e., gridbox temperature)
+      !--is updated due to latent heat cooling
+      zt(i) = zt(i) &
+              - qice_sedim * RLSTT / RCPD &
+              / ( 1. + RVTMP2 * ( zq(i) + zmqc(i) ) )
+    ENDIF ! ok_ice_sedim
+  ENDDO
 ENDIF
 
@@ -307,9 +329,9 @@
 
 SUBROUTINE histprecip_postcld( &
-        klon, dtime, iftop, paprsdn, paprsup, pplay, ctot_vol, ptconv, zdqsdT_raw, &
-        zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
+        klon, dtime, iftop, paprsdn, paprsup, pplay, ctot_vol, ptconv, pt_pron_clds, &
+        zdqsdT_raw, zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, icesed_flux, &
         rneb, znebprecipclr, znebprecipcld, zneb, tot_zneb, zrho_up, zvelo_up, &
         zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
-        zradocond, zradoice, dqrauto, dqsauto &
+        zradocond, zradoice, dqrauto, dqsauto, dqised &
         )
 
@@ -321,4 +343,5 @@
                           iflag_autoconversion, ok_radocond_snow, ok_bug_phase_lscp, &
                           niter_lscp
+USE lmdz_lscp_ini, ONLY : ok_ice_sedim, fallice_sedim, eps
 USE lmdz_lscp_tools, ONLY : fallice_velocity
 
@@ -335,4 +358,5 @@
 REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--volumic cloud fraction [-]
 LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--true if we are in a convective point
+LOGICAL, INTENT(IN),    DIMENSION(klon) :: pt_pron_clds   !--true if we are in a prognostic clouds point
 REAL,    INTENT(IN),    DIMENSION(klon) :: zdqsdT_raw     !--derivative of qsat w.r.t. temperature times L/Cp [SI]
 
@@ -360,4 +384,5 @@
 REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflclr        !--flux of snow gridbox-mean in clear sky [kg/s/m2]
 REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflcld        !--flux of snow gridbox-mean in cloudy air [kg/s/m2]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: icesed_flux    !--flux of sedimentated ice [kg/s/m2]
 
 REAL,    INTENT(OUT),   DIMENSION(klon) :: zradocond      !--condensed water used in the radiation scheme [kg/kg]
@@ -365,4 +390,5 @@
 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: dqised         !--ice water content tendency due to sedimentation of ice crystals [kg/kg/s]
 
 
@@ -390,4 +416,8 @@
 REAL, DIMENSION(klon) :: ziflprev
 
+! Local variable for sedimentation
+REAL :: iwcg, tempc, icevelo
+REAL :: qice_sedim
+
 ! Misc
 INTEGER :: i, n
@@ -398,4 +428,6 @@
 zqpreci(:) = 0.
 ziflprev(:) = 0.
+
+icesed_flux(:) = 0.
 
 
@@ -474,5 +506,5 @@
   ENDDO
 
-  CALL fallice_velocity(klon, iwc, zt, zrho, pplay, ptconv, zvelo)
+  CALL fallice_velocity(klon, iwc, zt, zrho, pplay, ptconv, pt_pron_clds, zvelo)
 
   DO i = 1, klon
@@ -656,4 +688,39 @@
 ENDDO
 
+!---------------------------------------------------------
+!--                  ICE SEDIMENTATION
+!---------------------------------------------------------
+IF ( ok_ice_sedim ) THEN
+
+  DO i = 1, klon
+    IF ( ( zoliqi(i) .GT. 0. ) .AND. ( rneb(i) .GT. eps ) .AND. pt_pron_clds(i) ) THEN
+
+      iwcg = MAX( zrho(i) * zoliqi(i) / zneb(i) * 1000., eps ) ! iwc in g/m3
+      tempc = zt(i) - RTT ! celcius temp
+      ! so-called 'empirical parameterization' in Stubenrauch et al. 2019
+      IF ( tempc .GE. -60. ) THEN
+        icevelo = EXP( 1.9714 + 0.00216078 * tempc &
+                - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) )
+      ELSE
+        icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15
+      ENDIF
+      icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s
+
+      qice_sedim = zoliqi(i) * MIN(1., icevelo * dtime / zdz(i))
+      !--We remove the ice that was sedimentated in the gridbox (it cannot sedimentate two
+      !--times in the same timestep)
+      !qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime)
+      !--Add tendencies
+      zoliqi(i) = zoliqi(i) - qice_sedim
+      zoliq(i)  = zoliq(i)  - qice_sedim
+      !--Convert to flux
+      icesed_flux(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime
+
+      !--Diagnostic tendencies
+      dqised(i) = dqised(i) - qice_sedim / dtime
+    ENDIF
+  ENDDO
+ENDIF
+
 ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
 ! if iflag_evap_prec>=4
@@ -713,9 +780,8 @@
 SUBROUTINE poprecip_precld( &
            klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
-           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, qice_sedim, &
+           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, icesed_flux, &
            cldfra, qvc, qliq, qice, &
            rain, rainclr, raincld, snow, snowclr, snowcld, &
-           dqreva, dqssub, dcfsed, dqvcsed &
-           )
+           dqreva, dqssub)
 
 USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac
@@ -749,5 +815,5 @@
 REAL,    INTENT(IN),    DIMENSION(klon) :: qvapclrup      !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
 REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [kg/kg]
-REAL,    INTENT(IN),    DIMENSION(klon) :: qice_sedim     !--ice water content that sedimentated from the layer above [kg/kg]
+REAL,    INTENT(IN),    DIMENSION(klon) :: icesed_flux    !--sedimentated ice water flux [kg/s/m2]
 
 REAL,    INTENT(INOUT), DIMENSION(klon) :: cldfra         !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
@@ -766,6 +832,4 @@
 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqreva         !--rain tendency due to evaporation [kg/kg/s]
 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dcfsed         !--cloud fraction tendency due to sedimentation of ice crystals [s-1]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqvcsed        !--cloud water vapor tendency due to sedimentation of ice crystals [kg/kg/s]
 
 
@@ -791,4 +855,6 @@
 !--Specific heat constant
 REAL :: cpair, cpw
+!--Specific ice water content sedimentated
+REAL :: qice_sedim
 
 !--Initialisation
@@ -796,8 +862,4 @@
 dqreva(:) = 0.
 dqssub(:) = 0.
-IF ( ok_ice_sedim .AND. ok_ice_supersat ) THEN
-  dcfsed(:) = 0.
-  dqvcsed(:) = 0.
-ENDIF
 
 !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
@@ -943,34 +1005,13 @@
     !--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
     !--added to the total water content of the gridbox
-    IF ( ok_ice_sedim .AND. ( qice_sedim(i) .GT. 0. ) ) THEN
+    IF ( icesed_flux(i) .GT. 0. ) THEN
+      qice_sedim = icesed_flux(i) / dhum_to_dflux(i)
       !--Vapor is updated after evaporation/sublimation (it is increased)
-      qvap(i) = qvap(i) + qice_sedim(i)
+      qvap(i) = qvap(i) + qice_sedim
       !--Air and precip temperature (i.e., gridbox temperature)
       !--is updated due to latent heat cooling
       temp(i) = temp(i) &
-              - qice_sedim(i) * RLSTT / RCPD &
+              - qice_sedim * RLSTT / RCPD &
               / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) )
-    
-      IF ( ok_ice_supersat ) THEN
-        !--If ice supersaturation is activated, the cloud properties are prognostic.
-        !--The falling ice is then considered a new cloud in the gridbox.
-        !--BEWARE with this parameterisation, we can create a new cloud with a much
-        !--different ice water content and water vapor content than the existing cloud
-        !--(if it exists). This implies than unphysical fluxes of ice and vapor
-        !--occur between the existing cloud and the newly formed cloud (from sedimentation).
-        !--Note also that currently, the precipitation scheme assume a maximum
-        !--random overlap, meaning that the new formed clouds will not be affected
-        !--by vertical wind shear.
-        ! Maybe we should reduce dcfsed?
-        dcfsed(i) = precipfracclr_tmp(i)
-        dqvcsed(i) = qvapclr * dcfsed(i)
-        !--Add tendencies
-        cldfra(i) = cldfra(i) + dcfsed(i)
-        qvc(i) = qvc(i) + dqvcsed(i)
-        qice(i) = qice(i) + qice_sedim(i)
-        !--Normalisation
-        dcfsed(i) = dcfsed(i) / dtime
-        dqvcsed(i) = dqvcsed(i) / dtime
-      ENDIF
     ENDIF ! ok_ice_sedim
 
@@ -1326,5 +1367,5 @@
 SUBROUTINE poprecip_postcld( &
            klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
-           temp, qvap, qliq, qice, icefrac, cldfra, qice_sedim, &
+           temp, qvap, qliq, qice, icefrac, cldfra, icesed_flux, &
            precipfracclr, precipfraccld, &
            rain, rainclr, raincld, snow, snowclr, snowcld, &
@@ -1347,5 +1388,5 @@
                           rain_fallspeed_clr, rain_fallspeed_cld,              &
                           snow_fallspeed_clr, snow_fallspeed_cld,              &
-                          ok_ice_sedim, ice_fallspeed_sedim
+                          ok_ice_sedim, fallice_sedim
 
 
@@ -1368,5 +1409,4 @@
 REAL,    INTENT(IN),    DIMENSION(klon) :: icefrac        !--ice fraction [-]
 REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--cloud fraction [-]
-REAL,    INTENT(INOUT), DIMENSION(klon) :: qice_sedim     !--quantity of sedimentated ice [kg/kg]
 
 REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
@@ -1382,4 +1422,5 @@
 REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
 
+REAL,    INTENT(OUT),   DIMENSION(klon) :: icesed_flux    !--flux of sedimentated ice [kg/s/m2]
 REAL,    INTENT(OUT),   DIMENSION(klon) :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
 REAL,    INTENT(OUT),   DIMENSION(klon) :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
@@ -1393,5 +1434,5 @@
 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqised         !--ice water content tendency due to sedimentation of ice crystals [kg/kg/s]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: dqised         !--ice water content tendency due to sedimentation of ice crystals [kg/kg/s]
 
 
@@ -1436,5 +1477,6 @@
 
 !--Ice sedimentation
-REAL :: dz, qice_sedim_tmp
+REAL :: iwcg, tempc, icevelo
+REAL :: dz, qice_sedim
 
 
@@ -1443,4 +1485,5 @@
 
 qzero(:)    = 0.
+icesed_flux(:) = 0.
 
 dqrcol(:)   = 0.
@@ -1453,5 +1496,4 @@
 dqrfreez(:) = 0.
 dqsfreez(:) = 0.
-IF ( ok_ice_sedim ) dqised(:) = 0.
 
 
@@ -1917,18 +1959,29 @@
   !---------------------------------------------------------
   IF ( ok_ice_sedim ) THEN
-    qice_sedim_tmp = 0.
-
-    IF ( temp(i) .LE. temp_nowater ) THEN
+
+    IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. eps ) ) THEN
+
       rho = pplay(i) / temp(i) / RD
+      iwcg = MAX( rho * qice(i) / cldfra(i) * 1000., eps ) ! iwc in g/m3
+      tempc = temp(i) - RTT ! celcius temp
+      ! so-called 'empirical parameterization' in Stubenrauch et al. 2019
+      IF ( tempc .GE. -60. ) THEN
+        icevelo = EXP( 1.9714 + 0.00216078 * tempc &
+                - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) )
+      ELSE
+        icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15
+      ENDIF
+      icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s
+
       dz = ( paprsdn(i) - paprsup(i) ) / RG / rho
-      qice_sedim_tmp = qice(i) * MIN(1., ice_fallspeed_sedim * dtime / dz)
+      qice_sedim = qice(i) * MIN(1., icevelo * dtime / dz)
       !--Add tendencies
-      qice(i) = qice(i) - qice_sedim_tmp
-    ENDIF
-
-    !--Diagnostic tendencies
-    dqised(i) = ( qice_sedim(i) - qice_sedim_tmp ) / dtime
-    !--Save quantity that will be sedimentated in the layer below
-    qice_sedim(i) = qice_sedim_tmp
+      qice(i) = qice(i) - qice_sedim
+      !--Convert to flux
+      icesed_flux(i) = qice_sedim * dhum_to_dflux(i)
+
+      !--Diagnostic tendencies
+      dqised(i) = dqised(i) - qice_sedim / dtime
+    ENDIF
   ENDIF
 
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_tools.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_tools.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_tools.f90	(revision 5609)
@@ -6,5 +6,5 @@
 
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,velo)
+SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,ptpronclds,velo)
 
     ! Ref:
@@ -16,7 +16,7 @@
     ! 3212–3234. https://doi.org/10.1029/2019MS001642
     
-    use lmdz_lscp_ini, only: iflag_vice, ffallv_con, ffallv_lsc
+    use lmdz_lscp_ini, only: iflag_vice, ffallv_con, ffallv_lsc, ffallv_issr
     use lmdz_lscp_ini, only: cice_velo, dice_velo
-    use lmdz_lscp_ini, only: ok_bug_ice_fallspeed
+    use lmdz_lscp_ini, only: ok_bug_ice_fallspeed, eps
 
     IMPLICIT NONE
@@ -28,4 +28,5 @@
     REAL, INTENT(IN), DIMENSION(klon) :: pres      ! air pressure [Pa]
     LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv    ! convective point  [-]
+    LOGICAL, INTENT(IN), DIMENSION(klon) :: ptpronclds! prognostic clouds point  [-]
 
     REAL, INTENT(OUT), DIMENSION(klon) :: velo    ! fallspeed velocity of crystals [m/s]
@@ -42,4 +43,6 @@
         IF (ptconv(i)) THEN
             fallv_tun=ffallv_con
+        ELSEIF (ptpronclds(i)) THEN
+            fallv_tun=ffallv_issr
         ELSE
             fallv_tun=ffallv_lsc
@@ -51,5 +54,5 @@
         ELSE
             ! AB the threshold is way too high, we reduce it
-            iwcg=MAX(iwc(i)*1000.,1E-10) ! iwc in g/m3. We set a minimum value to prevent from division by 0
+            iwcg=MAX(iwc(i)*1000.,eps) ! iwc in g/m3. We set a minimum value to prevent from division by 0
         ENDIF
         phpa=pres(i)/100.    ! pressure in hPa
Index: LMDZ6/branches/contrails/libf/phylmd/phyetat0_mod.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phyetat0_mod.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/phyetat0_mod.f90	(revision 5609)
@@ -18,11 +18,12 @@
   USE surface_data,     ONLY : type_ocean, version_ocean
   USE phyetat0_get_mod, ONLY : phyetat0_get, phyetat0_srf
-  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, phys_tstep, &
+  USE phys_state_var_mod, ONLY : ancien_ok, cwcon, clwcon, detr_therm, phys_tstep, &
        qsol, fevap, z0m, z0h, agesno, &
        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, &
-       cf_ancien, qvc_ancien, cfa_ancien, qta_ancien, radpas, radsol, rain_fall, ratqs, &
-       rnebcon, rugoro, sig1, snow_fall, bs_fall, solaire_etat0, sollw, sollwdown, &
+       cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien, qva_ancien, qia_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, &
        wake_deltat, wake_delta_pbl_TKE, delta_tsurf, beta_aridity, wake_fip, wake_pe, &
@@ -415,4 +416,5 @@
     ancien_ok=ancien_ok.AND.phyetat0_get(cf_ancien,"CFANCIEN","CFANCIEN",0.)
     ancien_ok=ancien_ok.AND.phyetat0_get(qvc_ancien,"QVCANCIEN","QVCANCIEN",0.)
+    found=phyetat0_get(cwcon,"CWCON","CWCON",0.)
   ELSE
     cf_ancien(:,:)=0.
@@ -423,8 +425,12 @@
   IF ( ok_plane_contrail ) THEN
     ancien_ok=ancien_ok.AND.phyetat0_get(cfa_ancien,"CFAANCIEN","CFAANCIEN",0.)
-    ancien_ok=ancien_ok.AND.phyetat0_get(cfa_ancien,"QTAANCIEN","QTAANCIEN",0.)
+    ancien_ok=ancien_ok.AND.phyetat0_get(pcf_ancien,"PCFANCIEN","PCFANCIEN",0.)
+    ancien_ok=ancien_ok.AND.phyetat0_get(qva_ancien,"QVAANCIEN","QVAANCIEN",0.)
+    ancien_ok=ancien_ok.AND.phyetat0_get(qia_ancien,"QIAANCIEN","QIAANCIEN",0.)
   ELSE
     cfa_ancien(:,:)=0.
-    qta_ancien(:,:)=0.
+    pcf_ancien(:,:)=0.
+    qva_ancien(:,:)=0.
+    qia_ancien(:,:)=0.
   ENDIF
 
@@ -458,5 +464,7 @@
   IF ( ok_plane_contrail ) THEN
     IF ( ( maxval(cfa_ancien).EQ.minval(cfa_ancien) ) .OR. &
-         ( maxval(qta_ancien).EQ.minval(qta_ancien) ) ) THEN
+         ( maxval(pcf_ancien).EQ.minval(pcf_ancien) ) .OR. &
+         ( maxval(qva_ancien).EQ.minval(qva_ancien) ) .OR. &
+         ( maxval(qia_ancien).EQ.minval(qia_ancien) ) ) THEN
        ancien_ok=.false.
      ENDIF
Index: LMDZ6/branches/contrails/libf/phylmd/phyredem.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phyredem.f90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/phyredem.f90	(revision 5609)
@@ -23,7 +23,7 @@
                                 prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien,      &
                                 ql_ancien, qs_ancien, qbs_ancien, cf_ancien, &
-                                qvc_ancien, cfa_ancien, qta_ancien,          &
-                                u_ancien, v_ancien,                          &
-                                clwcon, rnebcon, ratqs, pbl_tke,             &
+                                qvc_ancien, cfa_ancien, pcf_ancien,          &
+                                qva_ancien, qia_ancien, u_ancien, v_ancien,  &
+                                cwcon, clwcon, rnebcon, ratqs, pbl_tke,      &
                                 wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
                                 wake_deltat, wake_deltaq, wake_s, awake_s,   &
@@ -255,9 +255,12 @@
       CALL put_field(pass,"CFANCIEN", "CFANCIEN", cf_ancien)
       CALL put_field(pass,"QVCANCIEN", "QVCANCIEN", qvc_ancien)
+      CALL put_field(pass,"CWCON", "CWCON", cwcon)
     ENDIF
 
     IF ( ok_plane_contrail ) THEN
       CALL put_field(pass,"CFAANCIEN", "CFAANCIEN", cfa_ancien)
-      CALL put_field(pass,"QTAANCIEN", "QTAANCIEN", qta_ancien)
+      CALL put_field(pass,"PCFANCIEN", "PCFANCIEN", pcf_ancien)
+      CALL put_field(pass,"QVAANCIEN", "QVAANCIEN", qva_ancien)
+      CALL put_field(pass,"QIAANCIEN", "QIAANCIEN", qia_ancien)
     ENDIF
 
Index: LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90	(revision 5609)
@@ -669,6 +669,10 @@
       REAL, SAVE, ALLOCATABLE :: cfa_seri(:,:), d_cfa_dyn(:,:)
       !$OMP THREADPRIVATE(cfa_seri, d_cfa_dyn)
-      REAL, SAVE, ALLOCATABLE :: qta_seri(:,:), d_qta_dyn(:,:)
-      !$OMP THREADPRIVATE(qta_seri, d_qta_dyn)
+      REAL, SAVE, ALLOCATABLE :: pcf_seri(:,:), d_pcf_dyn(:,:)
+      !$OMP THREADPRIVATE(pcf_seri, d_pcf_dyn)
+      REAL, SAVE, ALLOCATABLE :: qva_seri(:,:), d_qva_dyn(:,:)
+      !$OMP THREADPRIVATE(qva_seri, d_qva_dyn)
+      REAL, SAVE, ALLOCATABLE :: qia_seri(:,:), d_qia_dyn(:,:)
+      !$OMP THREADPRIVATE(qia_seri, d_qia_dyn)
       REAL, SAVE, ALLOCATABLE :: flight_dist(:,:), flight_h2o(:,:)
       !$OMP THREADPRIVATE(flight_dist, flight_h2o)
@@ -677,6 +681,8 @@
       REAL, SAVE, ALLOCATABLE :: potcontfraP(:,:), potcontfraNP(:,:)
       !$OMP THREADPRIVATE(potcontfraP, potcontfraNP)
-      REAL, SAVE, ALLOCATABLE :: qice_cont(:,:), dcfa_cir(:,:), dqta_cir(:,:)
-      !$OMP THREADPRIVATE(qice_cont, dcfa_cir, dqta_cir)
+      REAL, SAVE, ALLOCATABLE :: qice_perscont(:,:)
+      !$OMP THREADPRIVATE(qice_perscont)
+      REAL, SAVE, ALLOCATABLE :: dcfa_cir(:,:), dqta_cir(:,:)
+      !$OMP THREADPRIVATE(dcfa_cir, dqta_cir)
       REAL, SAVE, ALLOCATABLE :: dcfa_ini(:,:), dqia_ini(:,:), dqta_ini(:,:)
       !$OMP THREADPRIVATE(dcfa_ini, dqia_ini, dqta_ini)
@@ -1244,9 +1250,12 @@
       ALLOCATE(d_q_avi(klon,klev))
       ALLOCATE(cfa_seri(klon,klev), d_cfa_dyn(klon,klev))
-      ALLOCATE(qta_seri(klon,klev), d_qta_dyn(klon,klev))
+      ALLOCATE(pcf_seri(klon,klev), d_pcf_dyn(klon,klev))
+      ALLOCATE(qva_seri(klon,klev), d_qva_dyn(klon,klev))
+      ALLOCATE(qia_seri(klon,klev), d_qia_dyn(klon,klev))
       ALLOCATE(flight_dist(klon,klev), flight_h2o(klon,klev))
       ALLOCATE(Tcritcont(klon,klev), qcritcont(klon,klev))
       ALLOCATE(potcontfraP(klon,klev), potcontfraNP(klon,klev))
-      ALLOCATE(qice_cont(klon,klev), dcfa_cir(klon,klev), dqta_cir(klon,klev))
+      ALLOCATE(qice_perscont(klon,klev))
+      ALLOCATE(dcfa_cir(klon,klev), dqta_cir(klon,klev))
       ALLOCATE(dcfa_ini(klon,klev), dqia_ini(klon,klev), dqta_ini(klon,klev))
       ALLOCATE(dcfa_sub(klon,klev), dqia_sub(klon,klev), dqta_sub(klon,klev))
@@ -1661,7 +1670,9 @@
 
 !-- LSCP - aviation and contrails variables
-      DEALLOCATE(d_q_avi, cfa_seri, d_cfa_dyn, qta_seri, d_qta_dyn, flight_dist, flight_h2o)
+      DEALLOCATE(cfa_seri, d_cfa_dyn, pcf_seri, d_pcf_dyn)
+      DEALLOCATE(qva_seri, d_qva_dyn, qia_seri, d_qia_dyn)
+      DEALLOCATE(d_q_avi, flight_dist, flight_h2o)
       DEALLOCATE(Tcritcont, qcritcont, potcontfraP, potcontfraNP)
-      DEALLOCATE(qice_cont, dcfa_cir, dqta_cir, dcfa_ini, dqia_ini, dqta_ini)
+      DEALLOCATE(qice_perscont, dcfa_cir, dqta_cir, dcfa_ini, dqia_ini, dqta_ini)
       DEALLOCATE(dcfa_sub, dqia_sub, dqta_sub, dcfa_mix, dqia_mix, dqta_mix)
       DEALLOCATE(cldfra_nocont, cldtau_nocont, cldemi_nocont, conttau, contemi)
Index: LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5609)
@@ -2192,8 +2192,16 @@
   TYPE(ctrl_out), SAVE :: o_dcfadyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     'dcfadyn', 'Dynamics linear contrail fraction tendency', 's-1', (/ ('',i=1,10) /))
-  TYPE(ctrl_out), SAVE :: o_qtaseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
-    'qtaseri', 'Linear contrail total specific humidity', 'kg/kg', (/ ('',i=1,10) /))
-  TYPE(ctrl_out), SAVE :: o_dqtadyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
-    'dqtadyn', 'Dynamics linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('',i=1,10) /))
+  TYPE(ctrl_out), SAVE :: o_pcfseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
+    'pcfseri', 'Contrail induced cirrus fraction', '-', (/ ('',i=1,10) /))
+  TYPE(ctrl_out), SAVE :: o_dpcfdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
+    'dpcfdyn', 'Dynamics contrail induced cirrus fraction tendency', 's-1', (/ ('',i=1,10) /))
+  TYPE(ctrl_out), SAVE :: o_qvaseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
+    'qvaseri', 'Linear contrail total specific humidity', 'kg/kg', (/ ('',i=1,10) /))
+  TYPE(ctrl_out), SAVE :: o_dqvadyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
+    'dqvadyn', 'Dynamics linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('',i=1,10) /))
+  TYPE(ctrl_out), SAVE :: o_qiaseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
+    'qiaseri', 'Linear contrail total specific humidity', 'kg/kg', (/ ('',i=1,10) /))
+  TYPE(ctrl_out), SAVE :: o_dqiadyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
+    'dqiadyn', 'Dynamics linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('',i=1,10) /))
   TYPE(ctrl_out), SAVE :: o_Tcritcont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     'Tcritcont', 'Temperature threshold for contrail formation', 'K', (/ ('',i=1,10) /))
@@ -2204,6 +2212,6 @@
   TYPE(ctrl_out), SAVE :: o_potcontfraNP = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     'potcontfraNP', 'Potential non-persistent contrail fraction', '-', (/ ('', i=1,10)/))
-  TYPE(ctrl_out), SAVE :: o_qice_cont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
-    'qice_cont', 'Linear contrails ice specific humidity', 'kg/kg', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_qice_perscont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'qice_perscont', 'Contrail induced cirrus ice specific humidity', 'kg/kg', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_dcfaini = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     'dcfaini', 'Initial formation linear contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
Index: LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90	(revision 5609)
@@ -229,7 +229,8 @@
          o_issrfra250to300, o_issrfra300to400, o_issrfra400to500, &
 !-- LSCP - aviation variables
-         o_cfaseri, o_dcfadyn, o_qtaseri, o_dqtadyn, o_dqavi, &
+         o_cfaseri, o_dcfadyn, o_pcfseri, o_dpcfdyn, &
+         o_qvaseri, o_dqvadyn, o_qiaseri, o_dqiadyn, o_dqavi, &
          o_Tcritcont, o_qcritcont, o_potcontfraP, o_potcontfraNP, &
-         o_flight_dist, o_flight_h2o, o_qice_cont, o_dcfacir, o_dqtacir, &
+         o_flight_dist, o_flight_h2o, o_qice_perscont, o_dcfacir, o_dqtacir, &
          o_dcfaini, o_dqiaini, o_dqtaini, o_dcfasub, o_dqiasub, o_dqtasub, &
          o_dcfamix, o_dqiamix, o_dqtamix, &
@@ -360,7 +361,8 @@
          issrfra100to150, issrfra150to200, issrfra200to250, &
          issrfra250to300, issrfra300to400, issrfra400to500, &
-         cfa_seri, d_cfa_dyn, qta_seri, d_qta_dyn, d_q_avi, &
+         cfa_seri, d_cfa_dyn, pcf_seri, d_pcf_dyn, &
+         qva_seri, d_qva_dyn, qia_seri, d_qia_dyn, d_q_avi, &
          Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
-         flight_dist, flight_h2o, qice_cont, dcfa_cir, dqta_cir, &
+         flight_dist, flight_h2o, qice_perscont, dcfa_cir, dqta_cir, &
          dcfa_ini, dqia_ini, dqta_ini, dcfa_sub, dqia_sub, dqta_sub, &
          dcfa_mix, dqia_mix, dqta_mix, &
@@ -2105,9 +2107,4 @@
            CALL histwrite_phy(o_dqsfreez, dqsfreez)
            CALL histwrite_phy(o_dqsrim, dqsrim)
-           IF ( ok_ice_sedim ) THEN
-            CALL histwrite_phy(o_dqised, dqised)
-            CALL histwrite_phy(o_dcfsed, dcfsed)
-            CALL histwrite_phy(o_dqvcsed, dqvcsed)
-           ENDIF
            ELSE
             CALL histwrite_phy(o_dqreva, dqreva)
@@ -2116,4 +2113,6 @@
             CALL histwrite_phy(o_dqsauto, dqsauto)
            ENDIF
+           CALL histwrite_phy(o_qsatl, qsatliq)
+           CALL histwrite_phy(o_qsati, qsatice)
        ENDIF
 
@@ -2141,6 +2140,4 @@
          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)
          CALL histwrite_phy(o_issrfra100to150, issrfra100to150)
          CALL histwrite_phy(o_issrfra150to200, issrfra150to200)
@@ -2149,4 +2146,9 @@
          CALL histwrite_phy(o_issrfra300to400, issrfra300to400)
          CALL histwrite_phy(o_issrfra400to500, issrfra400to500)
+         IF ( ok_ice_sedim ) THEN
+          CALL histwrite_phy(o_dqised, dqised)
+          CALL histwrite_phy(o_dcfsed, dcfsed)
+          CALL histwrite_phy(o_dqvcsed, dqvcsed)
+         ENDIF
        ENDIF
 !-- LSCP - aviation variables
@@ -2154,6 +2156,10 @@
          CALL histwrite_phy(o_cfaseri, cfa_seri)
          CALL histwrite_phy(o_dcfadyn, d_cfa_dyn)
-         CALL histwrite_phy(o_qtaseri, qta_seri)
-         CALL histwrite_phy(o_dqtadyn, d_qta_dyn)
+         CALL histwrite_phy(o_pcfseri, pcf_seri)
+         CALL histwrite_phy(o_dpcfdyn, d_pcf_dyn)
+         CALL histwrite_phy(o_qvaseri, qva_seri)
+         CALL histwrite_phy(o_dqvadyn, d_qva_dyn)
+         CALL histwrite_phy(o_qiaseri, qia_seri)
+         CALL histwrite_phy(o_dqiadyn, d_qia_dyn)
          CALL histwrite_phy(o_flight_dist, flight_dist)
          CALL histwrite_phy(o_Tcritcont, Tcritcont)
@@ -2161,5 +2167,5 @@
          CALL histwrite_phy(o_potcontfraP, potcontfraP)
          CALL histwrite_phy(o_potcontfraNP, potcontfraNP)
-         CALL histwrite_phy(o_qice_cont, qice_cont)
+         CALL histwrite_phy(o_qice_perscont, qice_perscont)
          CALL histwrite_phy(o_dcfacir, dcfa_cir)
          CALL histwrite_phy(o_dqtacir, dqta_cir)
@@ -2219,5 +2225,5 @@
        CALL histwrite_phy(o_dqlphy2d,  zx_tmp_fi2d)
 
-       IF (nqo.EQ.3) THEN
+       IF (nqo .GE. 3) THEN
        CALL histwrite_phy(o_dqsphy,  d_qx(:,:,isol))
        IF (vars_defined) CALL water_int(klon,klev,d_qx(:,:,isol),zmasse,zx_tmp_fi2d)
Index: LMDZ6/branches/contrails/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phys_state_var_mod.F90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/phys_state_var_mod.F90	(revision 5609)
@@ -94,6 +94,10 @@
       REAL, ALLOCATABLE, SAVE :: cf_ancien(:,:), qvc_ancien(:,:)
 !$OMP THREADPRIVATE(cf_ancien, qvc_ancien)
-      REAL, ALLOCATABLE, SAVE :: cfa_ancien(:,:), qta_ancien(:,:)
-!$OMP THREADPRIVATE(cfa_ancien, qta_ancien)
+      REAL, ALLOCATABLE, SAVE :: cfa_ancien(:,:), pcf_ancien(:,:)
+!$OMP THREADPRIVATE(cfa_ancien, pcf_ancien)
+      REAL, ALLOCATABLE, SAVE :: qva_ancien(:,:), qia_ancien(:,:)
+!$OMP THREADPRIVATE(qva_ancien, qia_ancien)
+      REAL, ALLOCATABLE, SAVE :: cwcon(:,:), cwcon0(:,:)
+!$OMP THREADPRIVATE(cwcon, cwcon0)
 !!! RomP >>>
       REAL, ALLOCATABLE, SAVE :: tr_ancien(:,:,:)
@@ -590,5 +594,7 @@
       ALLOCATE(u_ancien(klon,klev), v_ancien(klon,klev))
       ALLOCATE(cf_ancien(klon,klev), qvc_ancien(klon,klev))
-      ALLOCATE(cfa_ancien(klon,klev), qta_ancien(klon,klev))
+      ALLOCATE(cfa_ancien(klon,klev), pcf_ancien(klon,klev))
+      ALLOCATE(qva_ancien(klon,klev), qia_ancien(klon,klev))
+      ALLOCATE(cwcon(klon,klev), cwcon0(klon,klev))
 !!! Rom P >>>
       ALLOCATE(tr_ancien(klon,klev,nbtr))
@@ -817,5 +823,7 @@
       DEALLOCATE(qtc_cv,sigt_cv,detrain_cv,fm_cv)
       DEALLOCATE(u_ancien, v_ancien)
-      DEALLOCATE(cf_ancien, qvc_ancien, cfa_ancien, qta_ancien)
+      DEALLOCATE(cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien)
+      DEALLOCATE(qva_ancien, qia_ancien)
+      DEALLOCATE(cwcon, cwcon0)
       DEALLOCATE(tr_ancien)                           !RomP
       DEALLOCATE(ratqs, pbl_tke,coefh,coefm)
Index: LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90	(revision 5602)
+++ LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90	(revision 5609)
@@ -328,7 +328,7 @@
        dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
        !-- LSCP - aviation and contrails variables
-       cfa_seri, qta_seri, d_cfa_dyn, d_qta_dyn, &
-       d_q_avi, flight_dist, flight_h2o, &
-       qice_cont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
+       cfa_seri, pcf_seri, qva_seri, qia_seri, d_cfa_dyn, d_pcf_dyn, d_qva_dyn, d_qia_dyn, &
+       d_q_avi, flight_dist, flight_h2o, qice_perscont, &
+       Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
        cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
        conttau, contemi, contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont, &
@@ -518,6 +518,6 @@
     !
     ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
-    INTEGER,SAVE :: ivap, iliq, isol, ibs, icf, iqvc, icfa, iqta
-!$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, iqvc, icfa, iqta)
+    INTEGER,SAVE :: ivap, iliq, isol, ibs, icf, iqvc, icfa, ipcf, iqia, iqva
+!$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, iqvc, icfa, ipcf, iqia, iqva)
     !
     !
@@ -1363,8 +1363,10 @@
        isol = strIdx(tracers(:)%name, addPhase('H2O', 's'))
        ibs  = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
-       icf  = strIdx(tracers(:)%name, addPhase('H2O', 'f'))
-       iqvc = strIdx(tracers(:)%name, addPhase('H2O', 'c'))
-       icfa = strIdx(tracers(:)%name, addPhase('H2O', 'a'))
-       iqta = strIdx(tracers(:)%name, addPhase('H2O', 'q'))
+       icf  = strIdx(tracers(:)%name, 'CLDFRA')
+       iqvc = strIdx(tracers(:)%name, 'CLDVAP_g')
+       icfa = strIdx(tracers(:)%name, 'CONTFRA')
+       ipcf = strIdx(tracers(:)%name, 'PERSCONTFRA')
+       iqva = strIdx(tracers(:)%name, 'CONTWATER_g')
+       iqia = strIdx(tracers(:)%name, 'CONTWATER_s')
 !       CALL init_etat0_limit_unstruct
 !       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
@@ -1414,7 +1416,8 @@
        ENDIF
 
-       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.'
+       IF (ok_ice_supersat.AND.((icf.EQ.0).OR.(iqvc.EQ.0))) THEN
+          WRITE (lunout, *) ' ok_ice_supersat=y requires 5 tracers ', &
+               '(H2O_g, H2O_l, H2O_s, CLDFRA, CLDVAP_g) but not all are ', &
+               'provided. Might as well stop here.'
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
@@ -1433,7 +1436,8 @@
        ENDIF
 
-       IF (ok_plane_contrail.AND.(nqo.LT.6)) THEN
-          WRITE (lunout, *) ' ok_plane_contrail=y requires 6 H2O tracers ', &
-              '(H2O_g, H2O_l, H2O_s, H2O_f, H2O_c, H2O_a) but nqo=', nqo, '. Might as well stop here.'
+       IF (ok_plane_contrail.AND.((icfa.EQ.0).OR.(iqva.EQ.0).OR.(iqia.EQ.0).OR.(ipcf.EQ.0))) THEN
+          WRITE (lunout, *) ' ok_plane_contrail=y requires 8 tracers ', &
+              '(H2O_g, H2O_l, H2O_s, CLDFRA, CLDVAP_g, CONTFRA, PERSCONTFRA, CONTWATER_s) ', &
+              'but not all are provided. Might as well stop here.'
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
@@ -1446,12 +1450,10 @@
        ENDIF
 
-        IF (ok_bs) 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
-             abort_message='see above'
-             CALL abort_physic(modname,abort_message, 1)
-         ENDIF
-        ENDIF
+       IF (ok_bs.AND.(nqo.LT.4)) THEN
+          WRITE (lunout, *) 'activation of blowing snow needs a specific H2O tracer', &
+                            'but nqo=', nqo
+          abort_message='see above'
+          CALL abort_physic(modname,abort_message, 1)
+       ENDIF
 
        Ncvpaseq1 = 0
@@ -1631,4 +1633,7 @@
        rnebcon(:,:) = 0.0
        clwcon(:,:) = 0.0
+       !--AB for prognostic clouds
+       cwcon0(:,:) = 0.0
+       cwcon(:,:) = 0.0
 
        !
@@ -2454,27 +2459,28 @@
           q_seri(i,k)  = qx(i,k,ivap)
           ql_seri(i,k) = qx(i,k,iliq)
+          qs_seri(i,k) = 0.
           qbs_seri(i,k)= 0.
           cf_seri(i,k) = 0.
           qvc_seri(i,k)= 0.
           cfa_seri(i,k)= 0.
-          qta_seri(i,k)= 0.
+          pcf_seri(i,k)= 0.
+          qva_seri(i,k)= 0.
+          qia_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.
-          ELSE IF (nqo.EQ.3) THEN        !--vapour, liquid and ice
+          IF (nqo .GE. 3) THEN        !--vapour, liquid and ice
              qs_seri(i,k) = qx(i,k,isol)
-          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_supersat) THEN
-               cf_seri(i,k) = qx(i,k,icf)
-               qvc_seri(i,k) = qx(i,k,iqvc)
-             ENDIF
-             IF (ok_plane_contrail) THEN
-               cfa_seri(i,k) = qx(i,k,icfa)
-               qta_seri(i,k) = qx(i,k,iqta)
-             ENDIF
-             IF (ok_bs) THEN
-               qbs_seri(i,k)= qx(i,k,ibs)
-             ENDIF
+          ENDIF
+          IF (ok_bs) THEN
+             qbs_seri(i,k) = qx(i,k,ibs)
+          ENDIF
+          IF (ok_ice_supersat) THEN
+             cf_seri(i,k)  = qx(i,k,icf)
+             qvc_seri(i,k) = qx(i,k,iqvc)
+          ENDIF
+          IF (ok_plane_contrail) THEN
+             cfa_seri(i,k) = qx(i,k,icfa)
+             pcf_seri(i,k) = qx(i,k,ipcf)
+             qva_seri(i,k) = qx(i,k,iqva)
+             qia_seri(i,k) = qx(i,k,iqia)
           ENDIF
        ENDDO
@@ -2552,5 +2558,7 @@
        d_qvc_dyn(:,:)= (qvc_seri(:,:)-qvc_ancien(:,:))/phys_tstep
        d_cfa_dyn(:,:)= (cfa_seri(:,:)-cfa_ancien(:,:))/phys_tstep
-       d_qta_dyn(:,:)= (qta_seri(:,:)-qta_ancien(:,:))/phys_tstep
+       d_pcf_dyn(:,:)= (pcf_seri(:,:)-pcf_ancien(:,:))/phys_tstep
+       d_qva_dyn(:,:)= (qva_seri(:,:)-qva_ancien(:,:))/phys_tstep
+       d_qia_dyn(:,:)= (qia_seri(:,:)-qia_ancien(:,:))/phys_tstep
        CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
        d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
@@ -2575,5 +2583,7 @@
        d_qvc_dyn(:,:)= 0.0
        d_cfa_dyn(:,:)= 0.0
-       d_qta_dyn(:,:)= 0.0
+       d_pcf_dyn(:,:)= 0.0
+       d_qva_dyn(:,:)= 0.0
+       d_qia_dyn(:,:)= 0.0
        d_q_dyn2d(:)  = 0.0
        d_ql_dyn2d(:) = 0.0
@@ -2700,10 +2710,12 @@
             qi_seri_lscp(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
             qvc_seri(i,k) = qvc_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
-            qta_seri(i,k) = qta_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
+            qva_seri(i,k) = qva_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
+            qia_seri(i,k) = qia_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
           ELSE
             ql_seri_lscp(i,k) = 0.
             qi_seri_lscp(i,k) = 0.
             qvc_seri(i,k) = 0.
-            qta_seri(i,k) = 0.
+            qva_seri(i,k) = 0.
+            qia_seri(i,k) = 0.
           ENDIF
         ENDDO
@@ -3927,4 +3939,5 @@
         DO i = 1, klon
           qvc_seri(i,k) = qvc_seri(i,k) * q_seri(i,k)
+          cwcon0(i,k) = zqsat(i,k)
         ENDDO
       ENDDO
@@ -3933,5 +3946,6 @@
       DO k = 1, klev
         DO i = 1, klon
-          qta_seri(i,k) = qta_seri(i,k) * q_seri(i,k)
+          qva_seri(i,k) = qva_seri(i,k) * q_seri(i,k)
+          qia_seri(i,k) = qia_seri(i,k) * q_seri(i,k)
         ENDDO
       ENDDO
@@ -3945,6 +3959,6 @@
 
     CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay,omega, &
-         t_seri, q_seri, ql_seri_lscp, qi_seri_lscp, &
-         ptconv, rnebcon0, clwcon0, ratqs, sigma_qtherm, &
+         t_seri, q_seri, ql_seri_lscp, qi_seri_lscp, ratqs, sigma_qtherm, &
+         ptconv, rnebcon, cwcon, clwcon, rnebcon0, cwcon0, clwcon0, &
          d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, &
          pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb,  &
@@ -3960,6 +3974,6 @@
          dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
          dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
-         cfa_seri, qta_seri, flight_dist, flight_h2o, &
-         qice_cont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
+         cfa_seri, pcf_seri, qva_seri, qia_seri, flight_dist, flight_h2o, &
+         qice_perscont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
          cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
          qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
@@ -3967,4 +3981,5 @@
          dqised, dcfsed, dqvcsed)
 
+    IF (ok_ice_supersat) cwcon(:,:) = cwcon0(:,:)
 
     ELSE
@@ -4513,6 +4528,6 @@
                zfice, dNovrN, ptconv, rnebcon, clwcon, &
                !--AB contrails
-               cfa_seri, qice_cont, cldfra_nocont, cldtau_nocont, cldemi_nocont, &
-               conttau, contemi, cldh_nocont, contcov, &
+               cfa_seri, pcf_seri, qia_seri, qice_perscont, cldfra_nocont, &
+               cldtau_nocont, cldemi_nocont, conttau, contemi, cldh_nocont, contcov, &
                fiwp_nocont, fiwc_nocont, ref_ice_nocont)
 
@@ -5709,18 +5724,17 @@
              d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
           ENDIF
-          !--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
+          IF (ok_bs) THEN
+             d_qx(i,k,ibs) = ( qbs_seri(i,k) - qx(i,k,ibs) ) / phys_tstep
+          ENDIF
+          IF (ok_ice_supersat) THEN
              d_qx(i,k,icf) = ( cf_seri(i,k) - qx(i,k,icf) ) / phys_tstep
              d_qx(i,k,iqvc) = ( qvc_seri(i,k) - qx(i,k,iqvc) ) / phys_tstep
-             IF (nqo.ge.6 .and. ok_plane_contrail) THEN
-               d_qx(i,k,icfa) = ( cfa_seri(i,k) - qx(i,k,icfa) ) / phys_tstep
-               d_qx(i,k,iqta) = ( qta_seri(i,k) - qx(i,k,iqta) ) / phys_tstep
-             ENDIF
           ENDIF
-
-           IF (nqo.ge.4 .and. ok_bs) THEN
-             d_qx(i,k,ibs) = ( qbs_seri(i,k) - qx(i,k,ibs) ) / phys_tstep
+          IF (ok_plane_contrail) THEN
+             d_qx(i,k,icfa) = ( cfa_seri(i,k) - qx(i,k,icfa) ) / phys_tstep
+             d_qx(i,k,ipcf) = ( pcf_seri(i,k) - qx(i,k,ipcf) ) / phys_tstep
+             d_qx(i,k,iqva) = ( qva_seri(i,k) - qx(i,k,iqva) ) / phys_tstep
+             d_qx(i,k,iqia) = ( qia_seri(i,k) - qx(i,k,iqia) ) / phys_tstep
           ENDIF
-
        ENDDO
     ENDDO
@@ -5754,5 +5768,7 @@
     qvc_ancien(:,:)= qvc_seri(:,:)
     cfa_ancien(:,:)= cfa_seri(:,:)
-    qta_ancien(:,:)= qta_seri(:,:)
+    pcf_ancien(:,:)= pcf_seri(:,:)
+    qva_ancien(:,:)= qva_seri(:,:)
+    qia_ancien(:,:)= qia_seri(:,:)
     CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
     CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
