Changeset 5609 for LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90

Timestamp:

Apr 13, 2025, 7:10:19 PM (13 months ago)

Author:

aborella

Message:

• changed treatment of prognostic variables for prognostic clouds
• adapted sedimentation and autoconversion for prognostic cirrus clouds
• cloud mixing, ice sedimentation and ISSR diagnosis are now consistent with the water vapor PDF
• simplified assumptions for ice crystals deposition / sublimation
• first version of the coupling between prognostic cirrus clouds and deep convection
• added persistent contrail cirrus clouds in radiative diagnostics

File:

• LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90 (modified) (45 diffs)

LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90

@@ -94 +94 @@
 !**********************************************************************************
 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 +120 @@
 
 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 +144 @@
 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 +152 @@
 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 +158 @@
 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 +190 @@
 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 +206 @@
 !
 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 +228 @@
 INTEGER :: i
 LOGICAL :: ok_warm_cloud
-REAL, DIMENSION(klon) :: qtot, qcld, qzero
+REAL, DIMENSION(klon) :: qcld, qzero
 REAL, DIMENSION(klon) :: clrfra, qclr
 !
@@ -228 +236 @@
 !
 ! 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 +288 @@
   !--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 +294 @@
     !--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 +315 @@
       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 +331 @@
       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 +353 @@
       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 +466 @@
         dqia_mix(i) = 0.
         dqta_mix(i) = 0.
+      ELSE
+        contfra(i) = 0.
+        perscontfra(i) = 0.
+        qcont(i) = 0.
       ENDIF
 
@@ -379 +479 @@
       !--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 +494 @@
           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 +524 @@
           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 +571 @@
 
           !--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 +645 @@
         !--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 +663 @@
         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 +680 @@
 
         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 +701 @@
             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 +723 @@
       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 +787 @@
         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 +799 @@
         !--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 +811 @@
         !--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 +824 @@
         !--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 +840 @@
           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 +847 @@
           !--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 +870 @@
             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 +941 @@
         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 +956 @@
         !--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 +968 @@
         !--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 +982 @@
         !--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 +989 @@
           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 +996 @@
           !--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 +1021 @@
             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 +1051 @@
         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 +1135 @@
 
 
-      !--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 +1164 @@
       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 +1194 @@
   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 +1251 @@
         cldfra(i) = 0.
         contfra(i)= 0.
+        perscontfra(i) = 0.
         qcld(i)   = 0.
         qvc(i)    = 0.
+        qcont(i)  = 0.
         qincld(i) = qsat(i)
       ELSE
@@ -1036 +1265 @@
       ENDIF
 
+      IF ( perscontfra(i) .LT. eps ) perscontfra(i) = 0.
+
     ENDIF ! keepgoing
   ENDDO
@@ -1044 +1275 @@
 !**********************************************************************************
 
-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