Changeset 5406 for LMDZ6/trunk/libf


Ignore:
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:
3 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90

    r5396 r5406  
    782782   
    783783    ! --------------------------------------------------------------------
    784     ! P2> Cloud formation
     784    ! P3> Cloud formation
    785785    !---------------------------------------------------------------------
    786786    !
     
    793793    ! -------------------------------------------------------------------
    794794
    795         ! P2.1> With the PDFs (log-normal, bigaussian)
     795        ! P3.1> With the PDFs (log-normal, bigaussian)
    796796        ! cloud properties calculation with the initial values of t and q
    797797        ! ----------------------------------------------------------------
     
    965965                        pplay(:,k), paprs(:,k), paprs(:,k+1), &
    966966                        cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
    967                         shear(:), tke_dissip(:,k), cell_area(:), &
    968                         Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), keepgoing(:), &
    969                         rneb(:,k), zqn(:), qvc(:), issrfra(:,k), qissr(:,k), &
     967                        shear, tke_dissip(:,k), cell_area, &
     968                        Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, &
     969                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
    970970                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
    971971                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
     
    10381038        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
    10391039
    1040         ! P2.3> Final quantities calculation
     1040        ! P3.2> Final quantities calculation
    10411041        ! Calculated variables:
    10421042        !   rneb : cloud fraction
     
    11331133          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
    11341134        ENDIF
    1135    
    1136     !--AB Write diagnostics and tracers for ice supersaturation
    1137     IF ( ok_ice_supersat ) THEN
    1138       CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
    1139       CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
    1140 
    1141       DO i = 1, klon
    1142 
    1143         cf_seri(i,k) = rneb(i,k)
    1144 
    1145         IF ( .NOT. ok_unadjusted_clouds ) THEN
    1146           qvc(i) = zqs(i) * rneb(i,k)
    1147         ENDIF
    1148         IF ( zq(i) .GT. min_qParent ) THEN
    1149           rvc_seri(i,k) = qvc(i) / zq(i)
    1150         ELSE
    1151           rvc_seri(i,k) = min_ratio
    1152         ENDIF
    1153         !--The MIN barrier is NEEDED because of:
    1154         !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
    1155         !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
    1156         !--The MAX barrier is a safeguard that should not be activated
    1157         rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
    1158 
    1159         !--Diagnostics
    1160         gamma_cond(i,k) = gammasat(i)
    1161         IF ( issrfra(i,k) .LT. eps ) THEN
    1162           issrfra(i,k) = 0.
    1163           qissr(i,k) = 0.
    1164         ENDIF
    1165         subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
    1166         qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
    1167         IF ( subfra(i,k) .LT. eps ) THEN
    1168           subfra(i,k) = 0.
    1169           qsub(i,k) = 0.
    1170         ENDIF
    1171         qcld(i,k) = qvc(i) + zcond(i)
    1172         IF ( cf_seri(i,k) .LT. eps ) THEN
    1173           qcld(i,k) = 0.
    1174         ENDIF
    1175       ENDDO
    1176     ENDIF
    11771135
    11781136
    11791137    ! ----------------------------------------------------------------
    1180     ! P3> Precipitation formation
     1138    ! P4> Precipitation formation
    11811139    ! ----------------------------------------------------------------
    11821140
     
    15141472
    15151473
    1516     ! Calculation of cloud condensates amount seen by radiative scheme
    1517     !-----------------------------------------------------------------
     1474    !----------------------------------------------------------------------
     1475    ! P5> Calculation of cloud condensates amount seen by radiative scheme
     1476    !----------------------------------------------------------------------
    15181477
    15191478    ! Calculation of the concentration of condensates seen by the radiation scheme
     
    15471506
    15481507    ! ----------------------------------------------------------------
    1549     ! P4> Wet scavenging
     1508    ! P6> Wet scavenging
    15501509    ! ----------------------------------------------------------------
    15511510
     
    16021561
    16031562    ENDDO
     1563   
     1564    !------------------------------------------------------------
     1565    ! P7 > write diagnostics and outputs
     1566    !------------------------------------------------------------
     1567   
     1568    !--AB Write diagnostics and tracers for ice supersaturation
     1569    IF ( ok_ice_supersat ) THEN
     1570      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
     1571      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
     1572
     1573      DO i = 1, klon
     1574
     1575        IF ( zoliq(i) .LE. 0. ) THEN
     1576          !--If everything was precipitated, the remaining empty cloud is dissipated
     1577          !--and everything is transfered to the subsaturated clear sky region
     1578          rneb(i,k) = 0.
     1579          qvc(i) = 0.
     1580        ENDIF
     1581
     1582        cf_seri(i,k) = rneb(i,k)
     1583
     1584        IF ( .NOT. ok_unadjusted_clouds ) THEN
     1585          qvc(i) = zqs(i) * rneb(i,k)
     1586        ENDIF
     1587        IF ( zq(i) .GT. min_qParent ) THEN
     1588          rvc_seri(i,k) = qvc(i) / zq(i)
     1589        ELSE
     1590          rvc_seri(i,k) = min_ratio
     1591        ENDIF
     1592        !--The MIN barrier is NEEDED because of:
     1593        !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
     1594        !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
     1595        !--The MAX barrier is a safeguard that should not be activated
     1596        rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
     1597
     1598        !--Diagnostics
     1599        gamma_cond(i,k) = gammasat(i)
     1600        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
     1601        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
     1602        qcld(i,k) = qvc(i) + zfice(i) * zoliq(i)
     1603      ENDDO
     1604    ENDIF
    16041605
    16051606    ! Outputs:
  • LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.f90

    r5396 r5406  
    88
    99CONTAINS
    10 
    1110
    1211!**********************************************************************************
     
    122121USE lmdz_lscp_ini,        ONLY: lunout
    123122
    124 USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, &
     123USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, std_subl_pdf_lscp, &
    125124                                mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, &
    126                                 rhlmid_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, rhl0_pdf_lscp, &
     125                                std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, &
    127126                                coef_mixing_lscp, coef_shear_lscp, &
    128127                                chi_mixing_lscp, rho_ice
     
    151150REAL,     INTENT(IN)   , DIMENSION(klon) :: qsat          ! saturation specific humidity [kg/kg]
    152151REAL,     INTENT(IN)   , DIMENSION(klon) :: gamma_cond    ! condensation threshold w.r.t. qsat [-]
    153 REAL,     INTENT(INOUT), DIMENSION(klon) :: ratqs         ! ratio between the variance of the total water distribution and its average [-]
     152REAL,     INTENT(IN)   , DIMENSION(klon) :: ratqs         ! ratio between the variance of the total water distribution and its average [-]
    154153LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: keepgoing     ! .TRUE. if a new condensation loop should be computed
    155154!
     
    211210! for unadjusted clouds
    212211REAL :: qvapincld, qvapincld_new
    213 REAL :: qiceincld, qice_ratio
    214 REAL :: pres_sat, kappa
    215 REAL :: air_thermal_conduct, water_vapor_diff
    216 REAL :: iwc
    217 REAL :: iwc_log_inf100, iwc_log_sup100
    218 REAL :: iwc_inf100, alpha_inf100, coef_inf100
    219 REAL :: mu_sup100, sigma_sup100, coef_sup100
    220 REAL :: Dm_ice, rm_ice
     212REAL :: qiceincld
    221213!
    222214! for sublimation
     
    309301
    310302      pdf_std   = ratqs(i) * qtot(i)
    311       pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2. ) )
     303      pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2 ) )
    312304      pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
    313305      pdf_a     = pdf_delta / ( pdf_k * SQRT(2.) )
     
    342334        !--and the condensation process is slightly adapted
    343335        !--This can happen only if ok_weibull_warm_clouds = .TRUE.
    344         cldfra(i) = 0.
    345         qvc(i)    = 0.
    346         qcld(i)   = 0.
     336        ! AB WARM CLOUD
     337        !cldfra(i) = 0.
     338        !qcld(i)   = 0.
     339        !qvc(i)    = 0.
     340        cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
     341        qcld(i)   = MAX(0., MIN(qtot(i), ql_seri(i) + qi_seri(i) + rvc_seri(i) * qtot(i)))
     342        qvc(i)    = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
    347343        ok_warm_cloud = .TRUE.
    348344      ELSE
     
    384380      !--Cell dry air mass [kg]
    385381      !M_cell = rhodz * cell_area(i)
    386 
    387 
    388       IF ( ok_unadjusted_clouds ) THEN
    389         !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
    390         !--hypothesis is lost, and the vapor in the cloud is purely prognostic.
    391         !
    392         !--The deposition equation is
    393         !-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
    394         !--from Lohmann et al. (2016), where
    395         !--alpha is the deposition coefficient [-]
    396         !--mi is the mass of one ice crystal [kg]
    397         !--C is the capacitance of an ice crystal [m]
    398         !--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
    399         !--R_v is the specific gas constant for humid air [J/kg/K]
    400         !--T is the temperature [K]
    401         !--esi is the saturation pressure w.r.t. ice [Pa]
    402         !--Dv is the diffusivity of water vapor [m2/s]
    403         !--Ls is the specific latent heat of sublimation [J/kg/K]
    404         !--ka is the thermal conductivity of dry air [J/m/s/K]
    405         !
    406         !--alpha is a coefficient to take into account the fact that during deposition, a water
    407         !--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
    408         !--coherent (with the same structure). It has no impact for sublimation.
    409         !--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
    410         !--during deposition, and alpha = 1. during sublimation.
    411         !--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
    412         !-- C = capa_cond_cirrus * rm_ice
    413         !
    414         !--We have qice = Nice * mi, where Nice is the ice crystal
    415         !--number concentration per kg of moist air
    416         !--HYPOTHESIS 1: the ice crystals are spherical, therefore
    417         !-- mi = 4/3 * pi * rm_ice**3 * rho_ice
    418         !--HYPOTHESIS 2: the ice crystals are monodisperse with the
    419         !--initial radius rm_ice_0.
    420         !--NB. this is notably different than the assumption
    421         !--of a distributed qice in the cloud made in the sublimation process;
    422         !--should it be consistent?
    423         !
    424         !--As the deposition process does not create new ice crystals,
    425         !--and because we assume a same rm_ice value for all crystals
    426         !--therefore the sublimation process does not destroy ice crystals
    427         !--(or, in a limit case, it destroys all ice crystals), then
    428         !--Nice is a constant during the sublimation/deposition process.
    429         !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
    430         !
    431         !--The deposition equation then reads:
    432         !-- 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
    433         !-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
    434         !--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
    435         !--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
    436         !-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
    437         !--  *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 )
    438         !--and we have
    439         !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
    440         !-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
    441         !--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
    442         !
    443         !--This system of equations can be resolved with an exact
    444         !--explicit numerical integration, having one variable resolved
    445         !--explicitly, the other exactly. The exactly resolved variable is
    446         !--the one decreasing, so it is qvc if the process is
    447         !--condensation, qi if it is sublimation.
    448         !
    449         !--kappa is computed as an initialisation constant, as it depends only
    450         !--on temperature and other pre-computed values
    451         pres_sat = qsat(i) / ( EPS_W + ( 1. - EPS_W ) * qsat(i) ) * pplay(i)
    452         !--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
    453         air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
    454         !--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
    455         water_vapor_diff = 0.211 * ( temp(i) / RTT )**1.94 * ( 101325. / pplay(i) ) * 1.e-4
    456         kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat(i) &
    457               * ( RV * temp(i) / water_vapor_diff / pres_sat &
    458                 + RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) )
    459         !--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
    460       ENDIF
    461382
    462383
     
    491412          qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
    492413
    493           IF ( ok_unadjusted_clouds ) THEN
    494             !--Here, the initial vapor in the cloud is qvapincld, and we compute
    495             !--the new vapor qvapincld_new
    496 
    497             !--rm_ice formula from McFarquhar and Heymsfield (1997)
    498             iwc = qiceincld * rho * 1e3
    499             iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837)
    500             iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) )
    501             iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) )
    502 
    503             alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100
    504             coef_inf100 = iwc_inf100 * alpha_inf100**3. / 120.
    505 
    506             mu_sup100 = ( 5.2 + 0.0013 * ( temp(i) - RTT ) ) &
    507                       + ( 0.026 - 1.2E-3 * ( temp(i) - RTT ) ) * iwc_log_sup100
    508             sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp(i) - RTT ) ) &
    509                          + ( 0.018 - 2.1E-4 * ( temp(i) - RTT ) ) * iwc_log_sup100
    510             coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3 * mu_sup100 + 4.5 * sigma_sup100**2. )
    511 
    512             Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2. ) &
    513                    * coef_sup100 ) / ( coef_inf100 + coef_sup100 )
    514             rm_ice = Dm_ice / 2. * 1.E-6
    515 
    516             IF ( qvapincld .GE. qsat(i) ) THEN
    517               !--If the cloud is initially supersaturated
    518               !--Exact explicit integration (qvc exact, qice explicit)
    519               qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) &
    520                             * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2. )
    521             ELSE
    522               !--If the cloud is initially subsaturated
    523               !--Exact explicit integration (qice exact, qvc explicit)
    524               !--The barrier is set so that the resulting vapor in cloud
    525               !--cannot be greater than qsat
    526               !--qice_ratio is the ratio between the new ice content and
    527               !--the old one, it is comprised between 0 and 1
    528               qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2. * dtime * ( qsat(i) - qvapincld ) )
    529 
    530               IF ( qice_ratio .LT. 0. ) THEN
    531                 !--If all the ice has been sublimated, we sublimate
    532                 !--completely the cloud and do not activate the sublimation
    533                 !--process
    534                 !--Tendencies and diagnostics
    535                 dcf_sub(i) = - cldfra(i)
    536                 dqvc_sub(i) = - qvc(i)
    537                 dqi_sub(i) = - ( qcld(i) - qvc(i) )
    538 
    539                 cldfra(i) = 0.
    540                 qcld(i) = 0.
    541                 qvc(i) = 0.
    542 
    543                 !--The new vapor in cloud is set to 0 so that the
    544                 !--sublimation process does not activate
    545                 qvapincld_new = 0.
    546               ELSE
    547                 !--Else, the sublimation process is activated with the
    548                 !--diagnosed new cloud water vapor
    549                 !--The new vapor in the cloud is increased with the
    550                 !--sublimated ice
    551                 qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**(3./2.) )
    552                 !--The new vapor in the cloud cannot be greater than qsat
    553                 qvapincld_new = MIN(qvapincld_new, qsat(i))
    554               ENDIF ! qice_ratio .LT. 0.
    555             ENDIF ! qvapincld .GT. qsat(i)
     414          ! AB WARM CLOUD
     415          !IF ( ok_unadjusted_clouds ) THEN
     416          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
     417            CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
     418                                        pplay(i), dtime, qvapincld_new)
     419            IF ( qvapincld_new .EQ. 0. ) THEN
     420              !--If all the ice has been sublimated, we sublimate
     421              !--completely the cloud and do not activate the sublimation
     422              !--process
     423              !--Tendencies and diagnostics
     424              dcf_sub(i) = - cldfra(i)
     425              dqvc_sub(i) = - qvc(i)
     426              dqi_sub(i) = - ( qcld(i) - qvc(i) )
     427
     428              cldfra(i) = 0.
     429              qcld(i) = 0.
     430              qvc(i) = 0.
     431            ENDIF
    556432          ELSE
    557433            !--We keep the saturation adjustment hypothesis, and the vapor in the
     
    576452
    577453          !--If the vapor in cloud is below vapor needed for the cloud to survive
    578           IF ( qvapincld .LT. qvapincld_new ) THEN
     454          IF ( ( qvapincld .LT. qvapincld_new ) .OR. ( iflag_cloud_sublim_pdf .GE. 4 ) ) THEN
    579455            !--Sublimation of the subsaturated cloud
    580456            !--iflag_cloud_sublim_pdf selects the PDF of the ice water content
     
    589465
    590466              dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1
    591               dqt_sub    = - cldfra(i) * ( qvapincld_new**2. - qvapincld**2. ) &
     467              dqt_sub    = - cldfra(i) * ( qvapincld_new**2 - qvapincld**2 ) &
    592468                                       * pdf_e1 / 2.
    593469
     
    601477                                       + qvapincld - qvapincld_new * pdf_e1 )
    602478
    603             ELSEIF ( iflag_cloud_sublim_pdf .GE. 3 ) THEN
     479            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 3 ) THEN
    604480              !--Gamma distribution starting at qvapincld
    605481              pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld
     
    610486              dcf_sub(i) = - cldfra(i) * pdf_e1
    611487              dqt_sub    = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 )
     488
     489            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 4 ) THEN
     490              !--Normal distribution around qvapincld + qiceincld with width
     491              !--constant in the RHi space
     492              pdf_y = ( qvapincld_new - qvapincld - qiceincld ) &
     493                    / ( std_subl_pdf_lscp / 100. * qsat(i))
     494              pdf_e1 = 0.5 * ( 1. + ERF( pdf_y / SQRT(2.) ) )
     495              pdf_e2 = EXP( - pdf_y**2 / 2. ) / SQRT( 2. * RPI )
     496
     497              dcf_sub(i) = - cldfra(i) * pdf_e1
     498              dqt_sub    = - cldfra(i) * ( ( qvapincld + qiceincld ) * pdf_e1 &
     499                                         - std_subl_pdf_lscp / 100. * qsat(i) * pdf_e2 )
     500                                                                                       
    612501            ENDIF
    613502
     
    640529        !--New PDF
    641530        rhl_clr = qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100.
    642         rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp)
    643531
    644532        !--Calculation of the properties of the PDF
     
    648536        !--Coefficient for standard deviation:
    649537        !--  tuning coef * (clear sky area**0.25) * (function of temperature)
    650         pdf_e1 = beta_pdf_lscp &
    651                * ( ( 1. - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) &
    652                * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
    653         IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN
    654           pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp
     538        !pdf_e1 = beta_pdf_lscp &
     539        !       * ( ( 1. - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) &
     540        !       * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
     541        !--  tuning coef * (cell area**0.25) * (function of temperature)
     542        pdf_e1 = beta_pdf_lscp * ( ( 1. - cldfra(i) ) * cell_area(i) )**0.25 &
     543                * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
     544        IF ( rhl_clr .GT. 50. ) THEN
     545          pdf_std = ( pdf_e1 - std100_pdf_lscp ) * ( 100. - rhl_clr ) / 50. + std100_pdf_lscp
    655546        ELSE
    656           pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp
     547          pdf_std = pdf_e1 * rhl_clr / 50.
    657548        ENDIF
    658549        pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
    659         pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3
     550        pdf_alpha = EXP( rhl_clr / 100. ) * pdf_e3
     551        pdf_alpha = MIN(10., pdf_alpha)
     552       
     553        IF ( ok_warm_cloud ) THEN
     554          !--If the statistical scheme is activated, the standard deviation is adapted
     555          !--to depend on the pressure level. It is multiplied by ratqs, so that near the
     556          !--surface std is almost 0, and upper than about 450 hPa the std is left untouched
     557          pdf_std = pdf_std * ratqs(i)
     558        ENDIF
    660559
    661560        pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
    662         pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. ))
     561        pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2 ))
    663562        pdf_loc = rhl_clr - pdf_scale * pdf_e2
    664 
    665         !--Diagnostics of ratqs
    666         ratqs(i) = pdf_std / ( qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100. )
    667563
    668564        !--Calculation of the newly condensed water and fraction (pronostic)
     
    677573        qt_cond = ( pdf_e3 + pdf_loc * cf_cond ) * qsatl(i) / 100.
    678574
    679         IF ( ok_warm_cloud ) THEN
    680           !--If the statistical scheme is activated, the calculated increase is equal
    681           !--to the cloud fraction, we assume saturation adjustment, and the
    682           !--condensation process stops
    683           cldfra(i) = cf_cond
    684           qcld(i) = qt_cond
    685           qvc(i) = cldfra(i) * qsat(i)
    686 
    687         ELSEIF ( cf_cond .GT. eps ) THEN
     575        ! AB WARM CLOUD
     576        !IF ( ok_warm_cloud ) THEN
     577        !  !--If the statistical scheme is activated, the calculated increase is equal
     578        !  !--to the cloud fraction, we assume saturation adjustment, and the
     579        !  !--condensation process stops
     580        !  cldfra(i) = cf_cond
     581        !  qcld(i) = qt_cond
     582        !  qvc(i) = cldfra(i) * qsat(i)
     583
     584        !ELSEIF ( cf_cond .GT. eps ) THEN
     585        IF ( cf_cond .GT. eps ) THEN
    688586
    689587          dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond
     
    695593
    696594
    697           IF ( ok_unadjusted_clouds ) THEN
     595          ! AB WARM CLOUD
     596          !IF ( ok_unadjusted_clouds ) THEN
     597          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
    698598            !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute
    699599            !--the new vapor qvapincld. The timestep is divided by two because we do not
     
    701601            qvapincld = gamma_cond(i) * qsat(i)
    702602            qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
    703 
    704             !--rm_ice formula from McFarquhar and Heymsfield (1997)
    705             iwc = qiceincld * rho * 1e3
    706             iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837)
    707             iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) )
    708             iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) )
    709 
    710             alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100
    711             coef_inf100 = iwc_inf100 * alpha_inf100**3. / 120.
    712 
    713             mu_sup100 = ( 5.2 + 0.0013 * ( temp(i) - RTT ) ) &
    714                       + ( 0.026 - 1.2E-3 * ( temp(i) - RTT ) ) * iwc_log_sup100
    715             sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp(i) - RTT ) ) &
    716                          + ( 0.018 - 2.1E-4 * ( temp(i) - RTT ) ) * iwc_log_sup100
    717             coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3. * mu_sup100 + 4.5 * sigma_sup100**2. )
    718 
    719             Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2. ) &
    720                    * coef_sup100 ) / ( coef_inf100 + coef_sup100 )
    721             rm_ice = Dm_ice / 2. * 1.E-6
    722             !--As qvapincld is necessarily greater than qsat, we only
    723             !--use the exact explicit formulation
    724             !--Exact explicit version
    725             qvapincld = qsat(i) + ( qvapincld - qsat(i) ) &
    726                       * EXP( - depo_coef_cirrus * dtime / 2. * qiceincld / kappa / rm_ice**2. )
     603            CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
     604                                        pplay(i), dtime / 2., qvapincld_new)
     605            qvapincld = qvapincld_new
    727606          ELSE
    728607            !--We keep the saturation adjustment hypothesis, and the vapor in the
     
    735614          dqi_con(i) = dqt_con - dqvc_con(i)
    736615
    737           !--Note that the tendencies are NOT added because they are
    738           !--added after the mixing process. In the following, the gridbox fraction is
    739           !-- 1. - dcf_con(i), and the total water in the gridbox is
    740           !-- qtot(i) - dqi_con(i) - dqvc_con(i)
     616          !--Add tendencies
     617          cldfra(i) = MIN(1., cldfra(i) + dcf_con(i))
     618          qcld(i)   = MIN(qtot(i), qcld(i) + dqt_con)
     619          qvc(i)    = MIN(qcld(i), qvc(i) + dqvc_con(i))
    741620
    742621        ENDIF ! ok_warm_cloud, cf_cond .GT. eps
    743       ENDIF ! ( 1. - cldfra(i) ) .GT. eps
    744 
    745       !--If there is still clear sky, we diagnose the ISSR
    746       !--We recalculte the PDF properties (after the condensation process)
    747       IF ( ( ( 1. - dcf_con(i) - cldfra(i) ) .GT. eps ) .AND. .NOT. ok_warm_cloud ) THEN
    748         !--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
    749         qclr = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i)
    750 
    751         !--New PDF
    752         rhl_clr = qclr / ( 1. - dcf_con(i) - cldfra(i) ) / qsatl(i) * 100.
    753         rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp)
    754 
    755         !--Calculation of the properties of the PDF
    756         !--Parameterization from IAGOS observations
    757         !--pdf_e1 and pdf_e2 will be reused below
    758 
    759         !--Coefficient for standard deviation:
    760         !--  tuning coef * (clear sky area**0.25) * (function of temperature)
    761         pdf_e1 = beta_pdf_lscp &
    762                * ( ( 1. - dcf_con(i) - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) &
    763                * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
    764         IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN
    765           pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp
    766         ELSE
    767           pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp
    768         ENDIF
    769         pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
    770         pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3
    771 
    772         pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
    773         pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. ))
    774         pdf_loc = rhl_clr - pdf_scale * pdf_e2
    775 
     622       
    776623        !--We then calculate the part that is greater than qsat
    777         !--and consider it supersaturated
    778 
     624        !--and lower than gamma_cond * qsat, which is the ice supersaturated region
    779625        pdf_x = qsat(i) / qsatl(i) * 100.
    780626        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
    781627        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
    782628        pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
    783         issrfra(i) = EXP( - pdf_y ) * ( 1. - dcf_con(i) - cldfra(i) )
    784         qissr(i) = ( pdf_e3 * ( 1. - dcf_con(i) - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
    785       ENDIF
     629        issrfra(i) = EXP( - pdf_y ) * ( 1. - cldfra(i) )
     630        qissr(i) = ( pdf_e3 * ( 1. - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
     631
     632        issrfra(i) = issrfra(i) - dcf_con(i)
     633        qissr(i) = qissr(i) - dqvc_con(i) - dqi_con(i)
     634
     635      ENDIF ! ( 1. - cldfra(i) ) .GT. eps
    786636
    787637      !--Calculation of the subsaturated clear sky fraction and water
    788       subfra(i) = 1. - dcf_con(i) - cldfra(i) - issrfra(i)
    789       qsub(i) = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i) - qissr(i)
     638      subfra(i) = 1. - cldfra(i) - issrfra(i)
     639      qsub(i) = qtot(i) - qcld(i) - qissr(i)
    790640
    791641
     
    797647      !--but does not cover the entire mesh.
    798648      !
    799       IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) &
    800            .AND. .NOT. ok_warm_cloud ) THEN
     649      ! AB WARM CLOUD
     650      !IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) &
     651      !        .AND. .NOT. ok_warm_cloud ) THEN
     652      IF ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) ) THEN
    801653
    802654        !--Initialisation
     
    845697        !--             * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
    846698        !--and finally
    847         N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - dcf_con(i) - cldfra(i) )**2. &
    848                   * cell_area(i) * ( 1. - dcf_con(i) ) * bovera / RPI &
    849                   / ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )**2.
    850         !--where coef_mix_lscp = ( alpha * norm_constant )**2.
     699        N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - cldfra(i) )**2 &
     700                  * cell_area(i) * bovera / RPI &
     701                  / ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )**2
     702        !--where coef_mix_lscp = ( alpha * norm_constant )**2
    851703        !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer
    852704        !--In particular, it is 0 if cldfra = 1
    853         a_mix = SQRT( cell_area(i) * ( 1. - dcf_con(i) ) * cldfra(i) / bovera / N_cld_mix / RPI )
     705        a_mix = SQRT( cell_area(i) * cldfra(i) / bovera / N_cld_mix / RPI )
    854706
    855707        !--The time required for turbulent diffusion to homogenize a region of size
    856         !--L_mix is defined as (L_mix**2./tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
     708        !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
    857709        !--We compute L_mix and assume that the cloud is mixed over this length
    858         L_mix = SQRT( dtime**3. * pbl_eps(i) )
     710        L_mix = SQRT( dtime**3 * pbl_eps(i) )
    859711        !--The mixing length cannot be greater than the semi-minor axis. In this case,
    860712        !--the entire cloud is mixed.
     
    863715        !--The fraction of clear sky mixed is
    864716        !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
    865         envfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) &
    866                    * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2. )
     717        envfra_mix = N_cld_mix * RPI / cell_area(i) &
     718                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
    867719        !--The fraction of cloudy sky mixed is
    868720        !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
    869         cldfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) &
    870                    * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2. )
     721        cldfra_mix = N_cld_mix * RPI / cell_area(i) &
     722                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
    871723
    872724
     
    878730        !--The increase in size is
    879731        L_shear = coef_shear_lscp * shear(i) * dz * dtime
    880         !--therefore, the total increase in fraction is
     732        !--therefore, the fraction of clear sky mixed is
    881733        !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
    882         shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix &
    883                   / cell_area(i) / ( 1. - dcf_con(i) )
     734        !--and the fraction of cloud mixed is
     735        !-- N_cld_mix * ( (a * b) - (a - L_shear) * b ) * RPI / 2. / cell_area
     736        !--(note that they are equal)
     737        shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i)
    884738        !--and the environment and cloud mixed fractions are the same,
    885739        !--which we add to the previous calculated mixed fractions.
     
    906760        !--to this fraction (sigma_mix * cldfra_mix).
    907761        IF ( subfra(i) .GT. eps ) THEN
    908 
    909           IF ( ok_unadjusted_clouds ) THEN
    910             !--The subsaturated air is simply added to the cloud,
    911             !--with the corresponding cloud fraction
    912             !--If the cloud is too subsaturated, the sublimation process
    913             !--activated in the following timestep will reduce the cloud
    914             !--fraction
    915             dcf_mix_sub = subfra_mix
    916             dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
    917             dqvc_mix_sub = dqt_mix_sub
    918 
     762          !--We compute the total humidity in the mixed air, which
     763          !--can be either sub- or supersaturated.
     764          qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &
     765                    + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &
     766                    / ( subfra_mix + cldfra_mix * sigma_mix )
     767
     768          ! AB WARM CLOUD
     769          !IF ( ok_unadjusted_clouds ) THEN
     770          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
     771            qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * sigma_mix / cldfra(i) &
     772                      / ( subfra_mix + cldfra_mix * sigma_mix )
     773            CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
     774                                        pplay(i), dtime, qvapincld_new)
     775            IF ( qvapincld_new .EQ. 0. ) THEN
     776              !--If all the ice has been sublimated, we sublimate
     777              !--completely the cloud and do not activate the sublimation
     778              !--process
     779              !--Tendencies and diagnostics
     780              dcf_mix_sub = - sigma_mix * cldfra_mix
     781              dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
     782              dqvc_mix_sub = dcf_mix_sub * qvc(i) / cldfra(i)
     783            ELSE
     784              dcf_mix_sub = subfra_mix
     785              dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
     786              dqvc_mix_sub = dcf_mix_sub * qvapincld_new
     787            ENDIF
    919788          ELSE
    920             !--We compute the total humidity in the mixed air, which
    921             !--can be either sub- or supersaturated.
    922             qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &
    923                       + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &
    924                       / ( subfra_mix + cldfra_mix * sigma_mix )
    925 
    926789            IF ( qvapinmix .GT. qsat(i) ) THEN
    927790              !--If the mixed air is supersaturated, we condense the subsaturated
     
    944807        IF ( issrfra(i) .GT. eps ) THEN
    945808
    946           IF ( ok_unadjusted_clouds ) THEN
    947             !--The ice supersaturated air is simply added to the
    948             !--cloud, and supersaturated vapor will be deposited on the
    949             !--cloud ice crystals by the deposition process in the
    950             !--following timestep
     809          ! AB WARM CLOUD
     810          !IF ( ok_unadjusted_clouds ) THEN
     811          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
     812            qvapinmix = ( qissr(i) * issrfra_mix / issrfra(i) &
     813                      + qcld(i) * cldfra_mix * (1. - sigma_mix) / cldfra(i) ) &
     814                      / ( issrfra_mix + cldfra_mix * (1. -  sigma_mix) )
     815            qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * (1. - sigma_mix) / cldfra(i) &
     816                      / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) )
     817            CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
     818                                        pplay(i), dtime, qvapincld_new)
    951819            dcf_mix_issr = issrfra_mix
    952820            dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
    953             dqvc_mix_issr = dqt_mix_issr
     821            dqvc_mix_issr = dcf_mix_issr * qvapincld_new
    954822          ELSE
    955823            !--In this case, the additionnal vapor condenses
     
    971839        issrfra(i) = MAX(0., issrfra(i) - dcf_mix_issr)
    972840        qissr(i)   = MAX(0., qissr(i) - dqt_mix_issr)
    973         cldfra(i)  = MAX(0., MIN(1. - dcf_con(i), cldfra(i) + dcf_mix(i)))
    974         qcld(i)    = MAX(0., MIN(qtot(i) - dqi_con(i) - dqvc_con(i), qcld(i) + dqt_mix))
     841        cldfra(i)  = MAX(0., MIN(1., cldfra(i) + dcf_mix(i)))
     842        qcld(i)    = MAX(0., MIN(qtot(i), qcld(i) + dqt_mix))
    975843        qvc(i)     = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i)))
    976844
    977       ENDIF ! ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) )
    978 
    979       !--Finally, we add the tendencies of condensation
    980       cldfra(i) = MIN(1., cldfra(i) + dcf_con(i))
    981       qcld(i)   = MIN(qtot(i), qcld(i) + dqvc_con(i) + dqi_con(i))
    982       qvc(i)    = MIN(qcld(i), qvc(i) + dqvc_con(i))
     845      ENDIF ! ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) )
    983846
    984847
     
    1052915!**********************************************************************************
    1053916
     917SUBROUTINE deposition_sublimation( &
     918    qvapincld, qiceincld, temp, qsat, pplay, dtime, &
     919    qvapincld_new)
     920
     921USE lmdz_lscp_ini,        ONLY: RV, RLSTT, RTT, RD, EPS_W
     922USE lmdz_lscp_ini,        ONLY: eps
     923USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
     924
     925REAL,     INTENT(IN) :: qvapincld     !
     926REAL,     INTENT(IN) :: qiceincld     !
     927REAL,     INTENT(IN) :: temp          !
     928REAL,     INTENT(IN) :: qsat          !
     929REAL,     INTENT(IN) :: pplay         !
     930REAL,     INTENT(IN) :: dtime         ! time step [s]
     931REAL,     INTENT(OUT):: qvapincld_new !
     932
     933!
     934! for unadjusted clouds
     935REAL :: qice_ratio
     936REAL :: pres_sat, rho, kappa
     937REAL :: air_thermal_conduct, water_vapor_diff
     938REAL :: iwc
     939REAL :: iwc_log_inf100, iwc_log_sup100
     940REAL :: iwc_inf100, alpha_inf100, coef_inf100
     941REAL :: mu_sup100, sigma_sup100, coef_sup100
     942REAL :: Dm_ice, rm_ice
     943
     944!--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
     945!--hypothesis is lost, and the vapor in the cloud is purely prognostic.
     946!
     947!--The deposition equation is
     948!-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
     949!--from Lohmann et al. (2016), where
     950!--alpha is the deposition coefficient [-]
     951!--mi is the mass of one ice crystal [kg]
     952!--C is the capacitance of an ice crystal [m]
     953!--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
     954!--R_v is the specific gas constant for humid air [J/kg/K]
     955!--T is the temperature [K]
     956!--esi is the saturation pressure w.r.t. ice [Pa]
     957!--Dv is the diffusivity of water vapor [m2/s]
     958!--Ls is the specific latent heat of sublimation [J/kg/K]
     959!--ka is the thermal conductivity of dry air [J/m/s/K]
     960!
     961!--alpha is a coefficient to take into account the fact that during deposition, a water
     962!--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
     963!--coherent (with the same structure). It has no impact for sublimation.
     964!--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
     965!--during deposition, and alpha = 1. during sublimation.
     966!--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
     967!-- C = capa_cond_cirrus * rm_ice
     968!
     969!--We have qice = Nice * mi, where Nice is the ice crystal
     970!--number concentration per kg of moist air
     971!--HYPOTHESIS 1: the ice crystals are spherical, therefore
     972!-- mi = 4/3 * pi * rm_ice**3 * rho_ice
     973!--HYPOTHESIS 2: the ice crystals are monodisperse with the
     974!--initial radius rm_ice_0.
     975!--NB. this is notably different than the assumption
     976!--of a distributed qice in the cloud made in the sublimation process;
     977!--should it be consistent?
     978!
     979!--As the deposition process does not create new ice crystals,
     980!--and because we assume a same rm_ice value for all crystals
     981!--therefore the sublimation process does not destroy ice crystals
     982!--(or, in a limit case, it destroys all ice crystals), then
     983!--Nice is a constant during the sublimation/deposition process.
     984!-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
     985!
     986!--The deposition equation then reads:
     987!-- 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
     988!-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
     989!--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
     990!--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
     991!-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
     992!--  *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 )
     993!--and we have
     994!-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
     995!-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
     996!--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
     997!
     998!--This system of equations can be resolved with an exact
     999!--explicit numerical integration, having one variable resolved
     1000!--explicitly, the other exactly. The exactly resolved variable is
     1001!--the one decreasing, so it is qvc if the process is
     1002!--condensation, qi if it is sublimation.
     1003!
     1004!--kappa is computed as an initialisation constant, as it depends only
     1005!--on temperature and other pre-computed values
     1006pres_sat = qsat / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay
     1007!--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
     1008air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184
     1009!--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
     1010water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4
     1011kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat &
     1012      * ( RV * temp / water_vapor_diff / pres_sat &
     1013        + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) )
     1014!--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
     1015
     1016!--Dry density [kg/m3]
     1017rho = pplay / temp / RD
     1018
     1019!--Here, the initial vapor in the cloud is qvapincld, and we compute
     1020!--the new vapor qvapincld_new
     1021
     1022!--rm_ice formula from McFarquhar and Heymsfield (1997)
     1023iwc = qiceincld * rho * 1e3
     1024iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837)
     1025iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) )
     1026iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) )
     1027
     1028alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100
     1029coef_inf100 = iwc_inf100 * alpha_inf100**3 / 120.
     1030
     1031mu_sup100 = ( 5.2 + 0.0013 * ( temp - RTT ) ) &
     1032          + ( 0.026 - 1.2E-3 * ( temp - RTT ) ) * iwc_log_sup100
     1033sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp - RTT ) ) &
     1034             + ( 0.018 - 2.1E-4 * ( temp - RTT ) ) * iwc_log_sup100
     1035coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3. * mu_sup100 + 4.5 * sigma_sup100**2 )
     1036
     1037Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2 ) &
     1038       * coef_sup100 ) / ( coef_inf100 + coef_sup100 )
     1039rm_ice = Dm_ice / 2. * 1.E-6
     1040
     1041IF ( qvapincld .GE. qsat ) THEN
     1042  !--If the cloud is initially supersaturated
     1043  !--Exact explicit integration (qvc exact, qice explicit)
     1044  qvapincld_new = qsat + ( qvapincld - qsat ) &
     1045                * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
     1046ELSE
     1047  !--If the cloud is initially subsaturated
     1048  !--Exact explicit integration (qice exact, qvc explicit)
     1049  !--The barrier is set so that the resulting vapor in cloud
     1050  !--cannot be greater than qsat
     1051  !--qice_ratio is the ratio between the new ice content and
     1052  !--the old one, it is comprised between 0 and 1
     1053  qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2 * dtime * ( qsat - qvapincld ) )
     1054
     1055  IF ( qice_ratio .LT. 0. ) THEN
     1056    !--The new vapor in cloud is set to 0 so that the
     1057    !--sublimation process does not activate
     1058    qvapincld_new = 0.
     1059  ELSE
     1060    !--Else, the sublimation process is activated with the
     1061    !--diagnosed new cloud water vapor
     1062    !--The new vapor in the cloud is increased with the
     1063    !--sublimated ice
     1064    qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**1.5 )
     1065    !--The new vapor in the cloud cannot be greater than qsat
     1066    qvapincld_new = MIN(qvapincld_new, qsat)
     1067  ENDIF ! qice_ratio .LT. 0.
     1068ENDIF ! qvapincld .GT. qsat
     1069
     1070END SUBROUTINE deposition_sublimation
     1071
    10541072END MODULE lmdz_lscp_condensation
  • LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.f90

    r5383 r5406  
    154154  !$OMP THREADPRIVATE(ok_weibull_warm_clouds)
    155155
    156   INTEGER, SAVE, PROTECTED :: iflag_cloud_sublim_pdf=3       ! iflag for the distribution of water inside ice clouds
     156  INTEGER, SAVE, PROTECTED :: iflag_cloud_sublim_pdf=4       ! iflag for the distribution of water inside ice clouds
    157157  !$OMP THREADPRIVATE(iflag_cloud_sublim_pdf)
    158158
    159   REAL, SAVE, PROTECTED :: depo_coef_cirrus=.5               ! [-] deposition coefficient for growth of ice crystals in cirrus clouds
     159  REAL, SAVE, PROTECTED :: depo_coef_cirrus=.7               ! [-] deposition coefficient for growth of ice crystals in cirrus clouds
    160160  !$OMP THREADPRIVATE(depo_coef_cirrus)
    161161
     
    163163  !$OMP THREADPRIVATE(capa_cond_cirrus)
    164164
     165  REAL, SAVE, PROTECTED :: std_subl_pdf_lscp=2.              ! [%] standard deviation of the gaussian distribution of water inside ice clouds
     166  !$OMP THREADPRIVATE(std_subl_pdf_lscp)
     167
    165168  REAL, SAVE, PROTECTED :: mu_subl_pdf_lscp=1./3.            ! [-] shape factor of the gamma distribution of water inside ice clouds
    166169  !$OMP THREADPRIVATE(mu_subl_pdf_lscp)
    167170 
    168   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
     171  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
    169172  !$OMP THREADPRIVATE(beta_pdf_lscp)
    170173 
    171   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
     174  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
    172175  !$OMP THREADPRIVATE(temp_thresh_pdf_lscp)
    173176 
    174   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
    175   !$OMP THREADPRIVATE(rhlmid_pdf_lscp)
    176  
    177   REAL, SAVE, PROTECTED :: k0_pdf_lscp=2.80                  ! [-] factor for the PDF fit of water vapor in UTLS
     177  REAL, SAVE, PROTECTED :: k0_pdf_lscp=3.01                  ! [-] factor for the PDF fit of water vapor in UTLS
    178178  !$OMP THREADPRIVATE(k0_pdf_lscp)
    179179 
    180   REAL, SAVE, PROTECTED :: kappa_pdf_lscp=0.0236             ! [] factor for the PDF fit of water vapor in UTLS
     180  REAL, SAVE, PROTECTED :: kappa_pdf_lscp=0.0192             ! [] factor for the PDF fit of water vapor in UTLS
    181181  !$OMP THREADPRIVATE(kappa_pdf_lscp)
    182182 
    183   REAL, SAVE, PROTECTED :: rhl0_pdf_lscp=88.7                ! [%] factor for the PDF fit of water vapor in UTLS
    184   !$OMP THREADPRIVATE(rhl0_pdf_lscp)
     183  REAL, SAVE, PROTECTED :: std100_pdf_lscp=4.08              ! [%] standard deviation at RHliq=100% of the PDF fit of water vapor in UTLS
     184  !$OMP THREADPRIVATE(std100_pdf_lscp)
    185185 
    186186  REAL, SAVE, PROTECTED :: a_homofreez=2.349                 ! [-] factor for the Koop homogeneous freezing fit
     
    193193  !$OMP THREADPRIVATE(delta_hetfreez)
    194194 
    195   REAL, SAVE, PROTECTED :: coef_mixing_lscp=1e-7             ! [-] tuning coefficient for the mixing process
     195  REAL, SAVE, PROTECTED :: coef_mixing_lscp=9.E-8            ! [-] tuning coefficient for the mixing process
    196196  !$OMP THREADPRIVATE(coef_mixing_lscp)
    197197 
    198   REAL, SAVE, PROTECTED :: coef_shear_lscp=0.1               ! [-] additional coefficient for the shearing process (subprocess of the mixing process)
     198  REAL, SAVE, PROTECTED :: coef_shear_lscp=0.72              ! [-] additional coefficient for the shearing process (subprocess of the mixing process)
    199199  !$OMP THREADPRIVATE(coef_shear_lscp)
    200200 
    201   REAL, SAVE, PROTECTED :: chi_mixing_lscp=1.1               ! [-] factor for the macro distribution of ISSRs wrt clouds in a gridbox
     201  REAL, SAVE, PROTECTED :: chi_mixing_lscp=1.                ! [-] factor for the macro distribution of ISSRs wrt clouds in a gridbox
    202202  !$OMP THREADPRIVATE(chi_mixing_lscp)
    203 
    204 !  REAL, SAVE, PROTECTED :: contrail_cross_section=200000.
    205 !  !$OMP THREADPRIVATE(contrail_cross_section)
    206203  !--End of the parameters for condensation and ice supersaturation
    207204
     
    234231  !$OMP THREADPRIVATE(gamma_snwretro)
    235232
    236   REAL, SAVE, PROTECTED :: gamma_taud = 1.                  ! Tuning coeff for Lagrangian decorrelation timescale in lscp_icefrac_turb [-] 
     233  REAL, SAVE, PROTECTED :: gamma_taud = 1.                  ! Tuning coeff for Lagrangian decorrelation timescale in lscp_icefrac_turb [-]
    237234  !$OMP THREADPRIVATE(gamma_taud)
    238235
     
    245242  REAL, SAVE, PROTECTED :: gamma_rim=1.                     ! Tuning coefficient for riming efficiency (poprecip) [-]
    246243  !$OMP THREADPRIVATE(gamma_rim)
    247  
     244
    248245  REAL, SAVE, PROTECTED :: gamma_melt=1.                    ! Tuning coefficient for snow melting efficiency (poprecip) [-]
    249246  !$OMP THREADPRIVATE(gamma_melt)
     
    414411    CALL getin_p('depo_coef_cirrus',depo_coef_cirrus)
    415412    CALL getin_p('capa_cond_cirrus',capa_cond_cirrus)
     413    CALL getin_p('std_subl_pdf_lscp',std_subl_pdf_lscp)
    416414    CALL getin_p('mu_subl_pdf_lscp',mu_subl_pdf_lscp)
    417415    CALL getin_p('beta_pdf_lscp',beta_pdf_lscp)
    418416    CALL getin_p('temp_thresh_pdf_lscp',temp_thresh_pdf_lscp)
    419     CALL getin_p('rhlmid_pdf_lscp',rhlmid_pdf_lscp)
    420417    CALL getin_p('k0_pdf_lscp',k0_pdf_lscp)
    421418    CALL getin_p('kappa_pdf_lscp',kappa_pdf_lscp)
    422     CALL getin_p('rhl0_pdf_lscp',rhl0_pdf_lscp)
     419    CALL getin_p('std100_pdf_lscp',std100_pdf_lscp)
    423420    CALL getin_p('a_homofreez',a_homofreez)
    424421    CALL getin_p('b_homofreez',b_homofreez)
     
    427424    CALL getin_p('coef_shear_lscp',coef_shear_lscp)
    428425    CALL getin_p('chi_mixing_lscp',chi_mixing_lscp)
    429     !CALL getin_p('contrail_cross_section',contrail_cross_section)
    430426
    431427
     
    495491    WRITE(lunout,*) 'lscp_ini, depo_coef_cirrus:', depo_coef_cirrus
    496492    WRITE(lunout,*) 'lscp_ini, capa_cond_cirrus:', capa_cond_cirrus
     493    WRITE(lunout,*) 'lscp_ini, std_subl_pdf_lscp:', std_subl_pdf_lscp
    497494    WRITE(lunout,*) 'lscp_ini, mu_subl_pdf_lscp:', mu_subl_pdf_lscp
    498495    WRITE(lunout,*) 'lscp_ini, beta_pdf_lscp:', beta_pdf_lscp
    499496    WRITE(lunout,*) 'lscp_ini, temp_thresh_pdf_lscp:', temp_thresh_pdf_lscp
    500     WRITE(lunout,*) 'lscp_ini, rhlmid_pdf_lscp:', rhlmid_pdf_lscp
    501497    WRITE(lunout,*) 'lscp_ini, k0_pdf_lscp:', k0_pdf_lscp
    502498    WRITE(lunout,*) 'lscp_ini, kappa_pdf_lscp:', kappa_pdf_lscp
    503     WRITE(lunout,*) 'lscp_ini, rhl0_pdf_lscp:', rhl0_pdf_lscp
     499    WRITE(lunout,*) 'lscp_ini, std100_pdf_lscp:', std100_pdf_lscp
    504500    WRITE(lunout,*) 'lscp_ini, a_homofreez:', a_homofreez
    505501    WRITE(lunout,*) 'lscp_ini, b_homofreez:', b_homofreez
     
    508504    WRITE(lunout,*) 'lscp_ini, coef_shear_lscp:', coef_shear_lscp
    509505    WRITE(lunout,*) 'lscp_ini, chi_mixing_lscp:', chi_mixing_lscp
    510 !    WRITE(lunout,*) 'lscp_ini, contrail_cross_section:', contrail_cross_section
    511506
    512507
Note: See TracChangeset for help on using the changeset viewer.