Changeset 5406 for LMDZ6/trunk/libf

Timestamp:

Dec 13, 2024, 11:40:07 AM (2 days ago)

Author:

evignon

Message:

changement de la parametrisation de la dissipation des cirrus et ajouts de commentaires.

1. Borella

Location:

LMDZ6/trunk/libf/phylmd

Files:

• lmdz_lscp.f90 (modified) (8 diffs)
• lmdz_lscp_condensation.f90 (modified) (26 diffs)
• lmdz_lscp_ini.f90 (modified) (9 diffs)

LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90

@@ -782 +782 @@
     
     ! --------------------------------------------------------------------
-    ! P2> Cloud formation
+    ! P3> Cloud formation
     !---------------------------------------------------------------------
     !
@@ -793 +793 @@
     ! -------------------------------------------------------------------
 
-        ! P2.1> With the PDFs (log-normal, bigaussian) 
+        ! P3.1> With the PDFs (log-normal, bigaussian) 
         ! cloud properties calculation with the initial values of t and q
         ! ----------------------------------------------------------------
@@ -965 +965 @@
                         pplay(:,k), paprs(:,k), paprs(:,k+1), &
                         cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
-                        shear(:), tke_dissip(:,k), cell_area(:), &
-                        Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), keepgoing(:), &
-                        rneb(:,k), zqn(:), qvc(:), issrfra(:,k), qissr(:,k), &
+                        shear, tke_dissip(:,k), cell_area, &
+                        Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, &
+                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
                         dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
                         dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
@@ -1038 +1038 @@
         ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
 
-        ! P2.3> Final quantities calculation 
+        ! P3.2> Final quantities calculation 
         ! Calculated variables: 
         !   rneb : cloud fraction
@@ -1133 +1133 @@
           rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
         ENDIF
-    
-    !--AB Write diagnostics and tracers for ice supersaturation
-    IF ( ok_ice_supersat ) THEN
-      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
-      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
-
-      DO i = 1, klon
-
-        cf_seri(i,k) = rneb(i,k)
-
-        IF ( .NOT. ok_unadjusted_clouds ) THEN
-          qvc(i) = zqs(i) * rneb(i,k)
-        ENDIF
-        IF ( zq(i) .GT. min_qParent ) THEN
-          rvc_seri(i,k) = qvc(i) / zq(i)
-        ELSE
-          rvc_seri(i,k) = min_ratio
-        ENDIF
-        !--The MIN barrier is NEEDED because of:
-        !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
-        !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
-        !--The MAX barrier is a safeguard that should not be activated
-        rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
-
-        !--Diagnostics
-        gamma_cond(i,k) = gammasat(i)
-        IF ( issrfra(i,k) .LT. eps ) THEN
-          issrfra(i,k) = 0.
-          qissr(i,k) = 0.
-        ENDIF
-        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
-        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
-        IF ( subfra(i,k) .LT. eps ) THEN
-          subfra(i,k) = 0.
-          qsub(i,k) = 0.
-        ENDIF
-        qcld(i,k) = qvc(i) + zcond(i)
-        IF ( cf_seri(i,k) .LT. eps ) THEN
-          qcld(i,k) = 0.
-        ENDIF
-      ENDDO
-    ENDIF
 
 
     ! ----------------------------------------------------------------
-    ! P3> Precipitation formation 
+    ! P4> Precipitation formation 
     ! ----------------------------------------------------------------
 
@@ -1514 +1472 @@
 
 
-    ! Calculation of cloud condensates amount seen by radiative scheme
-    !-----------------------------------------------------------------
+    !----------------------------------------------------------------------
+    ! P5> Calculation of cloud condensates amount seen by radiative scheme
+    !----------------------------------------------------------------------
 
     ! Calculation of the concentration of condensates seen by the radiation scheme
@@ -1547 +1506 @@
 
     ! ----------------------------------------------------------------
-    ! P4> Wet scavenging
+    ! P6> Wet scavenging
     ! ----------------------------------------------------------------
 
@@ -1602 +1561 @@
 
     ENDDO
+    
+    !------------------------------------------------------------
+    ! P7 > write diagnostics and outputs
+    !------------------------------------------------------------
+   
+    !--AB Write diagnostics and tracers for ice supersaturation
+    IF ( ok_ice_supersat ) THEN
+      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
+      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
+
+      DO i = 1, klon
+
+        IF ( zoliq(i) .LE. 0. ) THEN
+          !--If everything was precipitated, the remaining empty cloud is dissipated
+          !--and everything is transfered to the subsaturated clear sky region
+          rneb(i,k) = 0.
+          qvc(i) = 0.
+        ENDIF
+
+        cf_seri(i,k) = rneb(i,k)
+
+        IF ( .NOT. ok_unadjusted_clouds ) THEN
+          qvc(i) = zqs(i) * rneb(i,k)
+        ENDIF
+        IF ( zq(i) .GT. min_qParent ) THEN
+          rvc_seri(i,k) = qvc(i) / zq(i)
+        ELSE
+          rvc_seri(i,k) = min_ratio
+        ENDIF
+        !--The MIN barrier is NEEDED because of:
+        !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
+        !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
+        !--The MAX barrier is a safeguard that should not be activated
+        rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
+
+        !--Diagnostics
+        gamma_cond(i,k) = gammasat(i)
+        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
+        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
+        qcld(i,k) = qvc(i) + zfice(i) * zoliq(i)
+      ENDDO
+    ENDIF
 
     ! Outputs:

LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.f90

@@ -8 +8 @@
 
 CONTAINS
-
 
 !**********************************************************************************
@@ -122 +121 @@
 USE lmdz_lscp_ini,        ONLY: lunout
 
-USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, &
+USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, std_subl_pdf_lscp, &
                                 mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, &
-                                rhlmid_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, rhl0_pdf_lscp, &
+                                std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, &
                                 coef_mixing_lscp, coef_shear_lscp, & 
                                 chi_mixing_lscp, rho_ice
@@ -151 +150 @@
 REAL,     INTENT(IN)   , DIMENSION(klon) :: qsat          ! saturation specific humidity [kg/kg]
 REAL,     INTENT(IN)   , DIMENSION(klon) :: gamma_cond    ! condensation threshold w.r.t. qsat [-]
-REAL,     INTENT(INOUT), DIMENSION(klon) :: ratqs         ! ratio between the variance of the total water distribution and its average [-]
+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
 !
@@ -211 +210 @@
 ! for unadjusted clouds
 REAL :: qvapincld, qvapincld_new
-REAL :: qiceincld, qice_ratio
-REAL :: pres_sat, 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
+REAL :: qiceincld
 !
 ! for sublimation
@@ -309 +301 @@
 
       pdf_std   = ratqs(i) * qtot(i)
-      pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2. ) )
+      pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2 ) )
       pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
       pdf_a     = pdf_delta / ( pdf_k * SQRT(2.) )
@@ -342 +334 @@
         !--and the condensation process is slightly adapted
         !--This can happen only if ok_weibull_warm_clouds = .TRUE.
-        cldfra(i) = 0.
-        qvc(i)    = 0.
-        qcld(i)   = 0.
+        ! AB WARM CLOUD
+        !cldfra(i) = 0.
+        !qcld(i)   = 0.
+        !qvc(i)    = 0.
+        cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
+        qcld(i)   = MAX(0., MIN(qtot(i), ql_seri(i) + qi_seri(i) + rvc_seri(i) * qtot(i)))
+        qvc(i)    = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
         ok_warm_cloud = .TRUE.
       ELSE
@@ -384 +380 @@
       !--Cell dry air mass [kg]
       !M_cell = rhodz * cell_area(i)
-
-
-      IF ( ok_unadjusted_clouds ) THEN
-        !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
-        !--hypothesis is lost, and the vapor in the cloud is purely prognostic.
-        !
-        !--The deposition equation is
-        !-- dmi/dt = alpha*4pi*C*Svi / ( R_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
-        kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat(i) &
-              * ( RV * temp(i) / water_vapor_diff / pres_sat &
-                + RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) )
-        !--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
-      ENDIF
 
 
@@ -491 +412 @@
           qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
 
-          IF ( ok_unadjusted_clouds ) THEN
-            !--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(i) - RTT ) ) &
-                      + ( 0.026 - 1.2E-3 * ( temp(i) - RTT ) ) * iwc_log_sup100
-            sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp(i) - RTT ) ) &
-                         + ( 0.018 - 2.1E-4 * ( temp(i) - 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(i) ) THEN
-              !--If the cloud is initially supersaturated
-              !--Exact explicit integration (qvc exact, qice explicit)
-              qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) &
-                            * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2. )
-            ELSE
-              !--If the cloud is initially subsaturated
-              !--Exact explicit integration (qice exact, qvc explicit)
-              !--The barrier is set so that the resulting vapor in cloud
-              !--cannot be greater than qsat
-              !--qice_ratio is the ratio between the new ice content and
-              !--the old one, it is comprised between 0 and 1
-              qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2. * dtime * ( qsat(i) - qvapincld ) )
-
-              IF ( qice_ratio .LT. 0. ) THEN
-                !--If all the ice has been sublimated, we sublimate
-                !--completely the cloud and do not activate the sublimation
-                !--process
-                !--Tendencies and diagnostics
-                dcf_sub(i) = - cldfra(i)
-                dqvc_sub(i) = - qvc(i)
-                dqi_sub(i) = - ( qcld(i) - qvc(i) )
-
-                cldfra(i) = 0.
-                qcld(i) = 0.
-                qvc(i) = 0.
-
-                !--The new vapor in cloud is set to 0 so that the
-                !--sublimation process does not activate
-                qvapincld_new = 0.
-              ELSE
-                !--Else, the sublimation process is activated with the
-                !--diagnosed new cloud water vapor
-                !--The new vapor in the cloud is increased with the
-                !--sublimated ice
-                qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**(3./2.) )
-                !--The new vapor in the cloud cannot be greater than qsat
-                qvapincld_new = MIN(qvapincld_new, qsat(i))
-              ENDIF ! qice_ratio .LT. 0.
-            ENDIF ! qvapincld .GT. qsat(i)
+          ! AB WARM CLOUD
+          !IF ( ok_unadjusted_clouds ) THEN
+          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
+            CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
+                                        pplay(i), dtime, qvapincld_new)
+            IF ( qvapincld_new .EQ. 0. ) THEN
+              !--If all the ice has been sublimated, we sublimate
+              !--completely the cloud and do not activate the sublimation
+              !--process
+              !--Tendencies and diagnostics
+              dcf_sub(i) = - cldfra(i)
+              dqvc_sub(i) = - qvc(i)
+              dqi_sub(i) = - ( qcld(i) - qvc(i) )
+
+              cldfra(i) = 0.
+              qcld(i) = 0.
+              qvc(i) = 0.
+            ENDIF
           ELSE
             !--We keep the saturation adjustment hypothesis, and the vapor in the
@@ -576 +452 @@
 
           !--If the vapor in cloud is below vapor needed for the cloud to survive
-          IF ( qvapincld .LT. qvapincld_new ) THEN
+          IF ( ( qvapincld .LT. qvapincld_new ) .OR. ( iflag_cloud_sublim_pdf .GE. 4 ) ) THEN
             !--Sublimation of the subsaturated cloud
             !--iflag_cloud_sublim_pdf selects the PDF of the ice water content
@@ -589 +465 @@
 
               dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1
-              dqt_sub    = - cldfra(i) * ( qvapincld_new**2. - qvapincld**2. ) &
+              dqt_sub    = - cldfra(i) * ( qvapincld_new**2 - qvapincld**2 ) &
                                        * pdf_e1 / 2.
 
@@ -601 +477 @@
                                        + qvapincld - qvapincld_new * pdf_e1 )
 
-            ELSEIF ( iflag_cloud_sublim_pdf .GE. 3 ) THEN
+            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 3 ) THEN
               !--Gamma distribution starting at qvapincld
               pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld
@@ -610 +486 @@
               dcf_sub(i) = - cldfra(i) * pdf_e1
               dqt_sub    = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 )
+
+            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 4 ) 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 )
+                                                                                       
             ENDIF
 
@@ -640 +529 @@
         !--New PDF
         rhl_clr = qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100.
-        rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp)
 
         !--Calculation of the properties of the PDF
@@ -648 +536 @@
         !--Coefficient for standard deviation:
         !--  tuning coef * (clear sky area**0.25) * (function of temperature)
-        pdf_e1 = beta_pdf_lscp &
-               * ( ( 1. - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) &
-               * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
-        IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN
-          pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp
+        !pdf_e1 = beta_pdf_lscp &
+        !       * ( ( 1. - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) &
+        !       * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
+        !--  tuning coef * (cell area**0.25) * (function of temperature)
+        pdf_e1 = beta_pdf_lscp * ( ( 1. - cldfra(i) ) * cell_area(i) )**0.25 &
+                * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
+        IF ( rhl_clr .GT. 50. ) THEN
+          pdf_std = ( pdf_e1 - std100_pdf_lscp ) * ( 100. - rhl_clr ) / 50. + std100_pdf_lscp
         ELSE
-          pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp
+          pdf_std = pdf_e1 * rhl_clr / 50.
         ENDIF
         pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
-        pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3
+        pdf_alpha = EXP( rhl_clr / 100. ) * pdf_e3
+        pdf_alpha = MIN(10., pdf_alpha)
+        
+        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)
-        pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. ))
+        pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2 ))
         pdf_loc = rhl_clr - pdf_scale * pdf_e2
-
-        !--Diagnostics of ratqs
-        ratqs(i) = pdf_std / ( qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100. )
 
         !--Calculation of the newly condensed water and fraction (pronostic)
@@ -677 +573 @@
         qt_cond = ( pdf_e3 + pdf_loc * cf_cond ) * qsatl(i) / 100.
 
-        IF ( ok_warm_cloud ) THEN
-          !--If the statistical scheme is activated, the calculated increase is equal
-          !--to the cloud fraction, we assume saturation adjustment, and the
-          !--condensation process stops
-          cldfra(i) = cf_cond
-          qcld(i) = qt_cond
-          qvc(i) = cldfra(i) * qsat(i)
-
-        ELSEIF ( cf_cond .GT. eps ) THEN
+        ! AB WARM CLOUD
+        !IF ( ok_warm_cloud ) THEN
+        !  !--If the statistical scheme is activated, the calculated increase is equal
+        !  !--to the cloud fraction, we assume saturation adjustment, and the
+        !  !--condensation process stops
+        !  cldfra(i) = cf_cond
+        !  qcld(i) = qt_cond
+        !  qvc(i) = cldfra(i) * qsat(i)
+
+        !ELSEIF ( cf_cond .GT. eps ) THEN
+        IF ( cf_cond .GT. eps ) THEN
 
           dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond
@@ -695 +593 @@
 
 
-          IF ( ok_unadjusted_clouds ) THEN
+          ! AB WARM CLOUD
+          !IF ( ok_unadjusted_clouds ) THEN
+          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
             !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute
             !--the new vapor qvapincld. The timestep is divided by two because we do not
@@ -701 +601 @@
             qvapincld = gamma_cond(i) * qsat(i)
             qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
-
-            !--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(i) - RTT ) ) &
-                      + ( 0.026 - 1.2E-3 * ( temp(i) - RTT ) ) * iwc_log_sup100
-            sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp(i) - RTT ) ) &
-                         + ( 0.018 - 2.1E-4 * ( temp(i) - 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
-            !--As qvapincld is necessarily greater than qsat, we only
-            !--use the exact explicit formulation
-            !--Exact explicit version
-            qvapincld = qsat(i) + ( qvapincld - qsat(i) ) &
-                      * EXP( - depo_coef_cirrus * dtime / 2. * qiceincld / kappa / rm_ice**2. )
+            CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
+                                        pplay(i), dtime / 2., qvapincld_new)
+            qvapincld = qvapincld_new
           ELSE
             !--We keep the saturation adjustment hypothesis, and the vapor in the
@@ -735 +614 @@
           dqi_con(i) = dqt_con - dqvc_con(i)
 
-          !--Note that the tendencies are NOT added because they are
-          !--added after the mixing process. In the following, the gridbox fraction is
-          !-- 1. - dcf_con(i), and the total water in the gridbox is
-          !-- qtot(i) - dqi_con(i) - dqvc_con(i)
+          !--Add tendencies
+          cldfra(i) = MIN(1., cldfra(i) + dcf_con(i))
+          qcld(i)   = MIN(qtot(i), qcld(i) + dqt_con)
+          qvc(i)    = MIN(qcld(i), qvc(i) + dqvc_con(i))
 
         ENDIF ! ok_warm_cloud, cf_cond .GT. eps
-      ENDIF ! ( 1. - cldfra(i) ) .GT. eps
-
-      !--If there is still clear sky, we diagnose the ISSR
-      !--We recalculte the PDF properties (after the condensation process)
-      IF ( ( ( 1. - dcf_con(i) - cldfra(i) ) .GT. eps ) .AND. .NOT. ok_warm_cloud ) THEN
-        !--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
-        qclr = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i)
-
-        !--New PDF
-        rhl_clr = qclr / ( 1. - dcf_con(i) - cldfra(i) ) / qsatl(i) * 100.
-        rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp)
-
-        !--Calculation of the properties of the PDF
-        !--Parameterization from IAGOS observations
-        !--pdf_e1 and pdf_e2 will be reused below
-
-        !--Coefficient for standard deviation:
-        !--  tuning coef * (clear sky area**0.25) * (function of temperature)
-        pdf_e1 = beta_pdf_lscp &
-               * ( ( 1. - dcf_con(i) - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) &
-               * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
-        IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN
-          pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp
-        ELSE
-          pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp
-        ENDIF
-        pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
-        pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3
-
-        pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
-        pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. ))
-        pdf_loc = rhl_clr - pdf_scale * pdf_e2
-
+        
         !--We then calculate the part that is greater than qsat
-        !--and consider it supersaturated
-
+        !--and lower than gamma_cond * qsat, which is the ice supersaturated region
         pdf_x = qsat(i) / qsatl(i) * 100.
         pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
         pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
         pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
-        issrfra(i) = EXP( - pdf_y ) * ( 1. - dcf_con(i) - cldfra(i) )
-        qissr(i) = ( pdf_e3 * ( 1. - dcf_con(i) - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
-      ENDIF
+        issrfra(i) = EXP( - pdf_y ) * ( 1. - cldfra(i) )
+        qissr(i) = ( pdf_e3 * ( 1. - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
+
+        issrfra(i) = issrfra(i) - dcf_con(i)
+        qissr(i) = qissr(i) - dqvc_con(i) - dqi_con(i)
+
+      ENDIF ! ( 1. - cldfra(i) ) .GT. eps
 
       !--Calculation of the subsaturated clear sky fraction and water
-      subfra(i) = 1. - dcf_con(i) - cldfra(i) - issrfra(i)
-      qsub(i) = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i) - qissr(i)
+      subfra(i) = 1. - cldfra(i) - issrfra(i)
+      qsub(i) = qtot(i) - qcld(i) - qissr(i)
 
 
@@ -797 +647 @@
       !--but does not cover the entire mesh.
       !
-      IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) &
-           .AND. .NOT. ok_warm_cloud ) THEN
+      ! AB WARM CLOUD
+      !IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) &
+      !        .AND. .NOT. ok_warm_cloud ) THEN
+      IF ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) ) THEN
 
         !--Initialisation
@@ -845 +697 @@
         !--             * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
         !--and finally
-        N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - dcf_con(i) - cldfra(i) )**2. &
-                  * cell_area(i) * ( 1. - dcf_con(i) ) * bovera / RPI &
-                  / ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )**2.
-        !--where coef_mix_lscp = ( alpha * norm_constant )**2.
+        N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - cldfra(i) )**2 &
+                  * cell_area(i) * bovera / RPI &
+                  / ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )**2
+        !--where coef_mix_lscp = ( alpha * norm_constant )**2
         !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer
         !--In particular, it is 0 if cldfra = 1
-        a_mix = SQRT( cell_area(i) * ( 1. - dcf_con(i) ) * cldfra(i) / bovera / N_cld_mix / RPI )
+        a_mix = SQRT( cell_area(i) * cldfra(i) / bovera / N_cld_mix / RPI )
 
         !--The time required for turbulent diffusion to homogenize a region of size
-        !--L_mix is defined as (L_mix**2./tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
+        !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
         !--We compute L_mix and assume that the cloud is mixed over this length
-        L_mix = SQRT( dtime**3. * pbl_eps(i) )
+        L_mix = SQRT( dtime**3 * pbl_eps(i) )
         !--The mixing length cannot be greater than the semi-minor axis. In this case,
         !--the entire cloud is mixed.
@@ -863 +715 @@
         !--The fraction of clear sky mixed is
         !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
-        envfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) &
-                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2. )
+        envfra_mix = N_cld_mix * RPI / cell_area(i) &
+                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
         !--The fraction of cloudy sky mixed is
         !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
-        cldfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) &
-                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2. )
+        cldfra_mix = N_cld_mix * RPI / cell_area(i) &
+                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
 
 
@@ -878 +730 @@
         !--The increase in size is
         L_shear = coef_shear_lscp * shear(i) * dz * dtime
-        !--therefore, the total increase in fraction is
+        !--therefore, the fraction of clear sky mixed is
         !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
-        shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix &
-                  / cell_area(i) / ( 1. - dcf_con(i) )
+        !--and the fraction of cloud mixed is
+        !-- N_cld_mix * ( (a * b) - (a - L_shear) * b ) * RPI / 2. / cell_area
+        !--(note that they are equal)
+        shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i)
         !--and the environment and cloud mixed fractions are the same,
         !--which we add to the previous calculated mixed fractions.
@@ -906 +760 @@
         !--to this fraction (sigma_mix * cldfra_mix).
         IF ( subfra(i) .GT. eps ) THEN
-
-          IF ( ok_unadjusted_clouds ) THEN
-            !--The subsaturated air is simply added to the cloud,
-            !--with the corresponding cloud fraction
-            !--If the cloud is too subsaturated, the sublimation process
-            !--activated in the following timestep will reduce the cloud
-            !--fraction
-            dcf_mix_sub = subfra_mix
-            dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
-            dqvc_mix_sub = dqt_mix_sub
-
+          !--We compute the total humidity in the mixed air, which
+          !--can be either sub- or supersaturated.
+          qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &
+                    + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &
+                    / ( subfra_mix + cldfra_mix * sigma_mix )
+
+          ! AB WARM CLOUD
+          !IF ( ok_unadjusted_clouds ) THEN
+          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
+            qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * sigma_mix / cldfra(i) &
+                      / ( subfra_mix + cldfra_mix * sigma_mix )
+            CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
+                                        pplay(i), dtime, qvapincld_new)
+            IF ( qvapincld_new .EQ. 0. ) THEN
+              !--If all the ice has been sublimated, we sublimate
+              !--completely the cloud and do not activate the sublimation
+              !--process
+              !--Tendencies and diagnostics
+              dcf_mix_sub = - sigma_mix * cldfra_mix
+              dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
+              dqvc_mix_sub = dcf_mix_sub * qvc(i) / cldfra(i)
+            ELSE
+              dcf_mix_sub = subfra_mix
+              dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
+              dqvc_mix_sub = dcf_mix_sub * qvapincld_new
+            ENDIF
           ELSE
-            !--We compute the total humidity in the mixed air, which
-            !--can be either sub- or supersaturated.
-            qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &
-                      + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &
-                      / ( subfra_mix + cldfra_mix * sigma_mix )
-
             IF ( qvapinmix .GT. qsat(i) ) THEN
               !--If the mixed air is supersaturated, we condense the subsaturated
@@ -944 +807 @@
         IF ( issrfra(i) .GT. eps ) THEN
 
-          IF ( ok_unadjusted_clouds ) THEN
-            !--The ice supersaturated air is simply added to the
-            !--cloud, and supersaturated vapor will be deposited on the
-            !--cloud ice crystals by the deposition process in the
-            !--following timestep
+          ! AB WARM CLOUD
+          !IF ( ok_unadjusted_clouds ) THEN
+          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
+            qvapinmix = ( qissr(i) * issrfra_mix / issrfra(i) &
+                      + qcld(i) * cldfra_mix * (1. - sigma_mix) / cldfra(i) ) &
+                      / ( issrfra_mix + cldfra_mix * (1. -  sigma_mix) )
+            qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * (1. - sigma_mix) / cldfra(i) &
+                      / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) )
+            CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
+                                        pplay(i), dtime, qvapincld_new)
             dcf_mix_issr = issrfra_mix
             dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
-            dqvc_mix_issr = dqt_mix_issr
+            dqvc_mix_issr = dcf_mix_issr * qvapincld_new
           ELSE
             !--In this case, the additionnal vapor condenses
@@ -971 +839 @@
         issrfra(i) = MAX(0., issrfra(i) - dcf_mix_issr)
         qissr(i)   = MAX(0., qissr(i) - dqt_mix_issr)
-        cldfra(i)  = MAX(0., MIN(1. - dcf_con(i), cldfra(i) + dcf_mix(i)))
-        qcld(i)    = MAX(0., MIN(qtot(i) - dqi_con(i) - dqvc_con(i), qcld(i) + dqt_mix))
+        cldfra(i)  = MAX(0., MIN(1., cldfra(i) + dcf_mix(i)))
+        qcld(i)    = MAX(0., MIN(qtot(i), qcld(i) + dqt_mix))
         qvc(i)     = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i)))
 
-      ENDIF ! ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) )
-
-      !--Finally, we add the tendencies of condensation
-      cldfra(i) = MIN(1., cldfra(i) + dcf_con(i))
-      qcld(i)   = MIN(qtot(i), qcld(i) + dqvc_con(i) + dqi_con(i))
-      qvc(i)    = MIN(qcld(i), qvc(i) + dqvc_con(i))
+      ENDIF ! ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) )
 
 
@@ -1052 +915 @@
 !**********************************************************************************
 
+SUBROUTINE deposition_sublimation( &
+    qvapincld, qiceincld, temp, qsat, pplay, dtime, &
+    qvapincld_new)
+
+USE lmdz_lscp_ini,        ONLY: RV, RLSTT, RTT, RD, EPS_W
+USE lmdz_lscp_ini,        ONLY: eps
+USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
+
+REAL,     INTENT(IN) :: qvapincld     ! 
+REAL,     INTENT(IN) :: qiceincld     ! 
+REAL,     INTENT(IN) :: temp          ! 
+REAL,     INTENT(IN) :: qsat          ! 
+REAL,     INTENT(IN) :: pplay         ! 
+REAL,     INTENT(IN) :: dtime         ! time step [s]
+REAL,     INTENT(OUT):: qvapincld_new ! 
+
+!
+! for unadjusted clouds
+REAL :: qice_ratio
+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_new = qsat + ( qvapincld - qsat ) &
+                * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
+ELSE
+  !--If the cloud is initially subsaturated
+  !--Exact explicit integration (qice exact, qvc explicit)
+  !--The barrier is set so that the resulting vapor in cloud
+  !--cannot be greater than qsat
+  !--qice_ratio is the ratio between the new ice content and
+  !--the old one, it is comprised between 0 and 1
+  qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2 * dtime * ( qsat - qvapincld ) )
+
+  IF ( qice_ratio .LT. 0. ) THEN
+    !--The new vapor in cloud is set to 0 so that the
+    !--sublimation process does not activate
+    qvapincld_new = 0.
+  ELSE
+    !--Else, the sublimation process is activated with the
+    !--diagnosed new cloud water vapor
+    !--The new vapor in the cloud is increased with the
+    !--sublimated ice
+    qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**1.5 )
+    !--The new vapor in the cloud cannot be greater than qsat
+    qvapincld_new = MIN(qvapincld_new, qsat)
+  ENDIF ! qice_ratio .LT. 0.
+ENDIF ! qvapincld .GT. qsat
+
+END SUBROUTINE deposition_sublimation
+
 END MODULE lmdz_lscp_condensation

LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.f90

@@ -154 +154 @@
   !$OMP THREADPRIVATE(ok_weibull_warm_clouds)
 
-  INTEGER, SAVE, PROTECTED :: iflag_cloud_sublim_pdf=3       ! iflag for the distribution of water inside ice clouds
+  INTEGER, SAVE, PROTECTED :: iflag_cloud_sublim_pdf=4       ! iflag for the distribution of water inside ice clouds
   !$OMP THREADPRIVATE(iflag_cloud_sublim_pdf)
 
-  REAL, SAVE, PROTECTED :: depo_coef_cirrus=.5               ! [-] deposition coefficient for growth of ice crystals in cirrus clouds
+  REAL, SAVE, PROTECTED :: depo_coef_cirrus=.7               ! [-] deposition coefficient for growth of ice crystals in cirrus clouds
   !$OMP THREADPRIVATE(depo_coef_cirrus)
 
@@ -163 +163 @@
   !$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 :: mu_subl_pdf_lscp=1./3.            ! [-] shape factor of the gamma distribution of water inside ice clouds
   !$OMP THREADPRIVATE(mu_subl_pdf_lscp)
   
-  REAL, SAVE, PROTECTED :: beta_pdf_lscp=8.75E-4             ! [] tuning coefficient for the standard deviation of the PDF of water vapor in the clear sky region
+  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
   !$OMP THREADPRIVATE(beta_pdf_lscp)
   
-  REAL, SAVE, PROTECTED :: temp_thresh_pdf_lscp=188.         ! [K] factor for the PDF fit of water vapor in UTLS - below this temperature, water vapor is homogeneously distributed in the clear sky region
+  REAL, SAVE, PROTECTED :: temp_thresh_pdf_lscp=189.         ! [K] factor for the PDF fit of water vapor in UTLS - below this temperature, water vapor is homogeneously distributed in the clear sky region
   !$OMP THREADPRIVATE(temp_thresh_pdf_lscp)
   
-  REAL, SAVE, PROTECTED :: rhlmid_pdf_lscp=52.8              ! [%] factor for the PDF fit of water vapor in UTLS - below this rel hum wrt liq, std increases with RHliq, above it decreases with RHliq
-  !$OMP THREADPRIVATE(rhlmid_pdf_lscp)
-  
-  REAL, SAVE, PROTECTED :: k0_pdf_lscp=2.80                  ! [-] factor for the PDF fit of water vapor in UTLS
+  REAL, SAVE, PROTECTED :: k0_pdf_lscp=3.01                  ! [-] factor for the PDF fit of water vapor in UTLS
   !$OMP THREADPRIVATE(k0_pdf_lscp)
   
-  REAL, SAVE, PROTECTED :: kappa_pdf_lscp=0.0236             ! [] factor for the PDF fit of water vapor in UTLS
+  REAL, SAVE, PROTECTED :: kappa_pdf_lscp=0.0192             ! [] factor for the PDF fit of water vapor in UTLS
   !$OMP THREADPRIVATE(kappa_pdf_lscp)
   
-  REAL, SAVE, PROTECTED :: rhl0_pdf_lscp=88.7                ! [%] factor for the PDF fit of water vapor in UTLS
-  !$OMP THREADPRIVATE(rhl0_pdf_lscp)
+  REAL, SAVE, PROTECTED :: std100_pdf_lscp=4.08              ! [%] standard deviation at RHliq=100% of the PDF fit of water vapor in UTLS
+  !$OMP THREADPRIVATE(std100_pdf_lscp)
   
   REAL, SAVE, PROTECTED :: a_homofreez=2.349                 ! [-] factor for the Koop homogeneous freezing fit
@@ -193 +193 @@
   !$OMP THREADPRIVATE(delta_hetfreez)
   
-  REAL, SAVE, PROTECTED :: coef_mixing_lscp=1e-7             ! [-] tuning coefficient for the mixing process
+  REAL, SAVE, PROTECTED :: coef_mixing_lscp=9.E-8            ! [-] tuning coefficient for the mixing process
   !$OMP THREADPRIVATE(coef_mixing_lscp)
   
-  REAL, SAVE, PROTECTED :: coef_shear_lscp=0.1               ! [-] additional coefficient for the shearing process (subprocess of the mixing process)
+  REAL, SAVE, PROTECTED :: coef_shear_lscp=0.72              ! [-] additional coefficient for the shearing process (subprocess of the mixing process)
   !$OMP THREADPRIVATE(coef_shear_lscp)
   
-  REAL, SAVE, PROTECTED :: chi_mixing_lscp=1.1               ! [-] factor for the macro distribution of ISSRs wrt clouds in a gridbox
+  REAL, SAVE, PROTECTED :: chi_mixing_lscp=1.                ! [-] factor for the macro distribution of ISSRs wrt clouds in a gridbox
   !$OMP THREADPRIVATE(chi_mixing_lscp)
-
-!  REAL, SAVE, PROTECTED :: contrail_cross_section=200000.
-!  !$OMP THREADPRIVATE(contrail_cross_section)
   !--End of the parameters for condensation and ice supersaturation
 
@@ -234 +231 @@
   !$OMP THREADPRIVATE(gamma_snwretro)
 
-  REAL, SAVE, PROTECTED :: gamma_taud = 1.                  ! Tuning coeff for Lagrangian decorrelation timescale in lscp_icefrac_turb [-] 
+  REAL, SAVE, PROTECTED :: gamma_taud = 1.                  ! Tuning coeff for Lagrangian decorrelation timescale in lscp_icefrac_turb [-]
   !$OMP THREADPRIVATE(gamma_taud)
 
@@ -245 +242 @@
   REAL, SAVE, PROTECTED :: gamma_rim=1.                     ! Tuning coefficient for riming efficiency (poprecip) [-]
   !$OMP THREADPRIVATE(gamma_rim)
-  
+
   REAL, SAVE, PROTECTED :: gamma_melt=1.                    ! Tuning coefficient for snow melting efficiency (poprecip) [-]
   !$OMP THREADPRIVATE(gamma_melt)
@@ -414 +411 @@
     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)
-    CALL getin_p('rhlmid_pdf_lscp',rhlmid_pdf_lscp)
     CALL getin_p('k0_pdf_lscp',k0_pdf_lscp)
     CALL getin_p('kappa_pdf_lscp',kappa_pdf_lscp)
-    CALL getin_p('rhl0_pdf_lscp',rhl0_pdf_lscp)
+    CALL getin_p('std100_pdf_lscp',std100_pdf_lscp)
     CALL getin_p('a_homofreez',a_homofreez)
     CALL getin_p('b_homofreez',b_homofreez)
@@ -427 +424 @@
     CALL getin_p('coef_shear_lscp',coef_shear_lscp)
     CALL getin_p('chi_mixing_lscp',chi_mixing_lscp)
-    !CALL getin_p('contrail_cross_section',contrail_cross_section)
 
 
@@ -495 +491 @@
     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
-    WRITE(lunout,*) 'lscp_ini, rhlmid_pdf_lscp:', rhlmid_pdf_lscp
     WRITE(lunout,*) 'lscp_ini, k0_pdf_lscp:', k0_pdf_lscp
     WRITE(lunout,*) 'lscp_ini, kappa_pdf_lscp:', kappa_pdf_lscp
-    WRITE(lunout,*) 'lscp_ini, rhl0_pdf_lscp:', rhl0_pdf_lscp
+    WRITE(lunout,*) 'lscp_ini, std100_pdf_lscp:', std100_pdf_lscp
     WRITE(lunout,*) 'lscp_ini, a_homofreez:', a_homofreez
     WRITE(lunout,*) 'lscp_ini, b_homofreez:', b_homofreez
@@ -508 +504 @@
     WRITE(lunout,*) 'lscp_ini, coef_shear_lscp:', coef_shear_lscp
     WRITE(lunout,*) 'lscp_ini, chi_mixing_lscp:', chi_mixing_lscp
-!    WRITE(lunout,*) 'lscp_ini, contrail_cross_section:', contrail_cross_section