MODULE lmdz_lscp_condensation !---------------------------------------------------------------- ! Module for condensation of clouds routines ! that are called in LSCP IMPLICIT NONE CONTAINS !********************************************************************************** SUBROUTINE condensation_lognormal( & klon, temp, qtot, qsat, gamma_cond, ratqs, & keepgoing, cldfra, qincld, qvc) !---------------------------------------------------------------------- ! This subroutine calculates the formation of clouds, using a ! statistical cloud scheme. It uses a generalised lognormal distribution ! of total water in the gridbox ! See Bony and Emanuel, 2001 !---------------------------------------------------------------------- USE lmdz_lscp_ini, ONLY: eps IMPLICIT NONE ! ! Input ! INTEGER, INTENT(IN) :: klon ! number of horizontal grid points ! REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K] REAL, INTENT(IN) , DIMENSION(klon) :: qtot ! total specific humidity (without precip) [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: qsat ! saturation specific humidity [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-] REAL, INTENT(IN) , DIMENSION(klon) :: ratqs ! ratio between the variance of the total water distribution and its average [-] LOGICAL, INTENT(IN) , DIMENSION(klon) :: keepgoing ! .TRUE. if a new condensation loop should be computed ! ! Output ! NB. those are in INOUT because of the convergence loop on temperature ! (in some cases, the values are not re-computed) but the values ! are never used explicitely ! REAL, INTENT(INOUT), DIMENSION(klon) :: cldfra ! cloud fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: qincld ! cloud-mean in-cloud total specific water [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: qvc ! gridbox-mean vapor in the cloud [kg/kg] ! ! Local ! INTEGER :: i REAL :: pdf_std, pdf_k, pdf_delta REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2 ! !--Loop on klon ! DO i = 1, klon !--If a new calculation of the condensation is needed, !--i.e., temperature has not yet converged (or the cloud is !--formed elsewhere) IF (keepgoing(i)) THEN pdf_std = ratqs(i) * qtot(i) pdf_k = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2 ) ) pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) ) pdf_a = pdf_delta / ( pdf_k * SQRT(2.) ) pdf_b = pdf_k / (2. * SQRT(2.)) pdf_e1 = pdf_a - pdf_b pdf_e1 = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 ) pdf_e1 = 1. - ERF(pdf_e1) pdf_e2 = pdf_a + pdf_b pdf_e2 = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 ) pdf_e2 = 1. - ERF(pdf_e2) IF ( pdf_e1 .LT. eps ) THEN cldfra(i) = 0. qincld(i) = qsat(i) !--AB grid-mean vapor in the cloud - we assume saturation adjustment qvc(i) = 0. ELSE cldfra(i) = 0.5 * pdf_e1 qincld(i) = qtot(i) * pdf_e2 / pdf_e1 !--AB grid-mean vapor in the cloud - we assume saturation adjustment qvc(i) = qsat(i) * cldfra(i) ENDIF ENDIF ! end keepgoing ENDDO ! end loop on i END SUBROUTINE condensation_lognormal !********************************************************************************** !********************************************************************************** SUBROUTINE condensation_ice_supersat( & klon, dtime, missing_val, pplay, paprsdn, paprsup, & cf_seri, rvc_seri, ql_seri, qi_seri, shear, pbl_eps, cell_area, & temp, qtot, qsat, gamma_cond, ratqs, keepgoing, & cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, & dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, & Tcontr, qcontr, qcontr2, fcontrN, fcontrP, flight_dist, flight_h2o, & dcf_avi, dqi_avi, dqvc_avi) !---------------------------------------------------------------------- ! This subroutine calculates the formation, evolution and dissipation ! of clouds, using a process-oriented treatment of the cloud properties ! (cloud fraction, vapor in the cloud, condensed water in the cloud). ! It allows for ice supersaturation in cold regions, in clear sky. ! If ok_unadjusted_clouds, it also allows for sub- and supersaturated ! cloud water vapors. ! It also allows for the formation and evolution of condensation trails ! (contrails) from aviation. ! Authors: Audran Borella, Etienne Vignon, Olivier Boucher ! April 2024 !---------------------------------------------------------------------- USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf, calc_gammasat, GAMMAINC USE lmdz_lscp_ini, ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI, EPS_W USE lmdz_lscp_ini, ONLY: eps, temp_nowater, ok_weibull_warm_clouds USE lmdz_lscp_ini, ONLY: ok_unadjusted_clouds, iflag_cloud_sublim_pdf USE lmdz_lscp_ini, ONLY: lunout USE lmdz_lscp_ini, ONLY: depo_coef_cirrus, capa_cond_cirrus, & mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, & rhlmid_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, rhl0_pdf_lscp, & coef_mixing_lscp, coef_shear_lscp, & chi_mixing_lscp, rho_ice IMPLICIT NONE ! ! Input ! INTEGER, INTENT(IN) :: klon ! number of horizontal grid points REAL, INTENT(IN) :: dtime ! time step [s] REAL, INTENT(IN) :: missing_val ! missing value for output ! REAL, INTENT(IN) , DIMENSION(klon) :: pplay ! layer pressure [Pa] REAL, INTENT(IN) , DIMENSION(klon) :: paprsdn ! pressure at the lower interface [Pa] REAL, INTENT(IN) , DIMENSION(klon) :: paprsup ! pressure at the upper interface [Pa] REAL, INTENT(IN) , DIMENSION(klon) :: cf_seri ! cloud fraction [-] REAL, INTENT(IN) , DIMENSION(klon) :: rvc_seri ! gridbox-mean water vapor in cloud [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: ql_seri ! specific liquid water content [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: qi_seri ! specific ice water content [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: shear ! vertical shear [s-1] REAL, INTENT(IN) , DIMENSION(klon) :: pbl_eps ! REAL, INTENT(IN) , DIMENSION(klon) :: cell_area ! REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K] REAL, INTENT(IN) , DIMENSION(klon) :: qtot ! total specific humidity (without precip) [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: qsat ! saturation specific humidity [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-] REAL, INTENT(INOUT), DIMENSION(klon) :: ratqs ! ratio between the variance of the total water distribution and its average [-] LOGICAL, INTENT(IN) , DIMENSION(klon) :: keepgoing ! .TRUE. if a new condensation loop should be computed ! ! Input for aviation ! REAL, INTENT(IN), DIMENSION(klon) :: flight_dist ! REAL, INTENT(IN), DIMENSION(klon) :: flight_h2o ! ! ! Output ! NB. cldfra and qincld should be outputed as cf_seri and qi_seri, ! or as tendencies (maybe in the future) ! NB. those are in INOUT because of the convergence loop on temperature ! (in some cases, the values are not re-computed) but the values ! are never used explicitely ! REAL, INTENT(INOUT), DIMENSION(klon) :: cldfra ! cloud fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: qincld ! cloud-mean in-cloud total specific water [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: qvc ! gridbox-mean vapor in the cloud [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: issrfra ! ISSR fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: qissr ! gridbox-mean ISSR specific water [kg/kg] ! ! Diagnostics for condensation and ice supersaturation ! NB. idem for the INOUT ! REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_sub ! cloud fraction tendency because of sublimation [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_con ! cloud fraction tendency because of condensation [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_mix ! cloud fraction tendency because of cloud mixing [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_adj ! specific ice content tendency because of temperature adjustment [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_sub ! specific ice content tendency because of sublimation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_con ! specific ice content tendency because of condensation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_mix ! specific ice content tendency because of cloud mixing [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_adj ! specific cloud water vapor tendency because of temperature adjustment [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_sub ! specific cloud water vapor tendency because of sublimation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_con ! specific cloud water vapor tendency because of condensation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s] ! ! Diagnostics for aviation ! NB. idem for the INOUT ! REAL, INTENT(INOUT), DIMENSION(klon) :: Tcontr ! critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) [K] REAL, INTENT(INOUT), DIMENSION(klon) :: qcontr ! REAL, INTENT(INOUT), DIMENSION(klon) :: qcontr2 ! REAL, INTENT(INOUT), DIMENSION(klon) :: fcontrN ! REAL, INTENT(INOUT), DIMENSION(klon) :: fcontrP ! REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_avi ! cloud fraction tendency because of aviation [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_avi ! specific ice content tendency because of aviation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_avi ! specific cloud water vapor tendency because of aviation [kg/kg/s] ! ! Local ! INTEGER :: i LOGICAL :: ok_warm_cloud REAL, DIMENSION(klon) :: qcld, qzero ! ! for lognormal REAL :: pdf_std, pdf_k, pdf_delta REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2 ! ! for unadjusted clouds REAL :: qvapincld, qvapincld_new REAL :: qiceincld, qice_ratio REAL :: pres_sat, kappa REAL :: air_thermal_conduct, water_vapor_diff REAL :: 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 ! ! for sublimation REAL :: pdf_alpha REAL :: dqt_sub ! ! for condensation REAL, DIMENSION(klon) :: qsatl, dqsatl REAL :: clrfra, qclr, sl_clr, rhl_clr REAL :: pdf_ratqs, pdf_skew, pdf_scale, pdf_loc REAL :: pdf_x, pdf_y, pdf_T REAL :: pdf_e3, pdf_e4 REAL :: cf_cond, qt_cond, dqt_con ! ! for mixing REAL, DIMENSION(klon) :: subfra, qsub REAL :: dqt_mix_sub, dqt_mix_issr REAL :: dcf_mix_sub, dcf_mix_issr REAL :: dqvc_mix_sub, dqvc_mix_issr REAL :: dqt_mix REAL :: a_mix, bovera, N_cld_mix, L_mix REAL :: envfra_mix, cldfra_mix REAL :: L_shear, shear_fra REAL :: sigma_mix, issrfra_mix, subfra_mix, qvapinmix ! ! for cell properties REAL :: rho, rhodz, dz !REAL :: V_cell, M_cell ! ! for aviation and cell properties !REAL :: dqt_avi !REAL :: contrail_fra ! ! !--more local variables for diagnostics !--imported from YOMCST.h !--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1) !--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air !--values from Schumann, Meteorol Zeitschrift, 1996 !--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen !--Qheat = 43. / 50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen !REAL, PARAMETER :: EiH2O=1.25 !--emission index of water vapour for kerosene (kg kg-1) !REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1) !REAL, PARAMETER :: eta=0.3 !--average propulsion efficiency of the aircraft !--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute !--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010. !REAL :: Gcontr !--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) !--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1) !REAL :: qsatliqcontr !----------------------------------------------- ! Initialisations !----------------------------------------------- ! Ajout des émissions de H2O dues à l'aviation ! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q ! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O ) ! = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) ) ! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q) ! flight_h2O is in kg H2O / s / cell ! !IF (ok_plane_h2o) THEN ! q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) ) !ENDIF qzero(:) = 0. !--Calculation of qsat w.r.t. liquid CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 1, .FALSE., qsatl, dqsatl) ! !--Loop on klon ! DO i = 1, klon !--If a new calculation of the condensation is needed, !--i.e., temperature has not yet converged (or the cloud is !--formed elsewhere) IF (keepgoing(i)) THEN !--If the temperature is higher than the threshold below which !--there is no liquid in the gridbox, we activate the usual scheme !--(generalised lognormal from Bony and Emanuel 2001) !--If ok_weibull_warm_clouds = .TRUE., the Weibull law is used for !--all clouds, and the lognormal scheme is not activated IF ( ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds ) THEN pdf_std = ratqs(i) * qtot(i) pdf_k = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2. ) ) pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) ) pdf_a = pdf_delta / ( pdf_k * SQRT(2.) ) pdf_b = pdf_k / (2. * SQRT(2.)) pdf_e1 = pdf_a - pdf_b pdf_e1 = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 ) pdf_e1 = 1. - ERF(pdf_e1) pdf_e2 = pdf_a + pdf_b pdf_e2 = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 ) pdf_e2 = 1. - ERF(pdf_e2) IF ( pdf_e1 .LT. eps ) THEN cldfra(i) = 0. qincld(i) = qsat(i) qvc(i) = 0. ELSE cldfra(i) = 0.5 * pdf_e1 qincld(i) = qtot(i) * pdf_e2 / pdf_e1 qvc(i) = qsat(i) * cldfra(i) ENDIF !--If the temperature is lower than temp_nowater, we use the new !--condensation scheme that allows for ice supersaturation ELSE !--Initialisation IF ( temp(i) .GT. temp_nowater ) THEN !--If the air mass is warm (liquid water can exist), !--all the memory is lost and the scheme becomes statistical, !--i.e., the sublimation and mixing processes are deactivated, !--and the condensation process is slightly adapted !--This can happen only if ok_weibull_warm_clouds = .TRUE. cldfra(i) = 0. qvc(i) = 0. qcld(i) = 0. ok_warm_cloud = .TRUE. ELSE !--The following barriers ensure that the traced cloud properties !--are consistent. In some rare cases, i.e. the cloud water vapor !--can be greater than the total water in the gridbox cldfra(i) = MAX(0., MIN(1., cf_seri(i))) qcld(i) = MAX(0., MIN(qtot(i), rvc_seri(i) * qtot(i) + qi_seri(i))) qvc(i) = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i))) ok_warm_cloud = .FALSE. ENDIF dcf_sub(i) = 0. dqi_sub(i) = 0. dqvc_sub(i) = 0. dqi_adj(i) = 0. dqvc_adj(i) = 0. dcf_con(i) = 0. dqi_con(i) = 0. dqvc_con(i) = 0. dcf_mix(i) = 0. dqi_mix(i) = 0. dqvc_mix(i) = 0. issrfra(i) = 0. qissr(i) = 0. subfra(i) = 0. qsub(i) = 0. !--Initialisation of the cell properties !--Dry density [kg/m3] rho = pplay(i) / temp(i) / RD !--Dry air mass [kg/m2] rhodz = ( paprsdn(i) - paprsup(i) ) / RG !--Cell thickness [m] dz = rhodz / rho !--Cell volume [m3] !V_cell = dz * cell_area(i) !--Cell dry air mass [kg] !M_cell = rhodz * cell_area(i) IF ( ok_unadjusted_clouds ) THEN !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment !--hypothesis is lost, and the vapor in the cloud is purely prognostic. ! !--The deposition equation is !-- dmi/dt = alpha*4pi*C*Svi / ( R_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 !------------------------------------------------------------------- !-- SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CLOUD -- !------------------------------------------------------------------- !--If there is a cloud IF ( cldfra(i) .GT. eps ) THEN qvapincld = qvc(i) / cldfra(i) qiceincld = ( qcld(i) / cldfra(i) - qvapincld ) !--If the ice water content is too low, the cloud is purely sublimated !--Most probably, we advected a cloud with no ice water content (possible !--if the entire cloud precipited for example) IF ( qiceincld .LT. eps ) THEN dcf_sub(i) = - cldfra(i) dqvc_sub(i) = - qvc(i) dqi_sub(i) = - ( qcld(i) - qvc(i) ) cldfra(i) = 0. qcld(i) = 0. qvc(i) = 0. !--Else, the cloud is adjusted and sublimated ELSE !--The vapor in cloud cannot be higher than the !--condensation threshold qvapincld = MIN(qvapincld, gamma_cond(i) * qsat(i)) qiceincld = ( qcld(i) / cldfra(i) - qvapincld ) IF ( 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) ELSE !--We keep the saturation adjustment hypothesis, and the vapor in the !--cloud is set equal to the saturation vapor qvapincld_new = qsat(i) ENDIF ! ok_unadjusted_clouds !--Adjustment of the IWC to the new vapor in cloud !--(this can be either positive or negative) dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) ) dqi_adj(i) = - dqvc_adj(i) !--Add tendencies !--The vapor in the cloud is updated, but not qcld as it is constant !--through this process, as well as cldfra which is unmodified qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i))) !------------------------------------ !-- DISSIPATION OF THE CLOUD -- !------------------------------------ !--If the vapor in cloud is below vapor needed for the cloud to survive IF ( qvapincld .LT. qvapincld_new ) THEN !--Sublimation of the subsaturated cloud !--iflag_cloud_sublim_pdf selects the PDF of the ice water content !--to use. !--iflag = 1 --> uniform distribution !--iflag = 2 --> exponential distribution !--iflag = 3 --> gamma distribution (Karcher et al 2018) IF ( iflag_cloud_sublim_pdf .EQ. 1 ) THEN !--Uniform distribution starting at qvapincld pdf_e1 = 1. / ( 2. * qiceincld ) dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1 dqt_sub = - cldfra(i) * ( qvapincld_new**2. - qvapincld**2. ) & * pdf_e1 / 2. ELSEIF ( iflag_cloud_sublim_pdf .EQ. 2 ) THEN !--Exponential distribution starting at qvapincld pdf_alpha = 1. / qiceincld pdf_e1 = EXP( - pdf_alpha * ( qvapincld_new - qvapincld ) ) dcf_sub(i) = - cldfra(i) * ( 1. - pdf_e1 ) dqt_sub = - cldfra(i) * ( ( 1. - pdf_e1 ) / pdf_alpha & + qvapincld - qvapincld_new * pdf_e1 ) ELSEIF ( iflag_cloud_sublim_pdf .GE. 3 ) THEN !--Gamma distribution starting at qvapincld pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld pdf_y = pdf_alpha * ( qvapincld_new - qvapincld ) pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y ) pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y ) dcf_sub(i) = - cldfra(i) * pdf_e1 dqt_sub = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 ) ENDIF !--Tendencies and diagnostics dqvc_sub(i) = dqt_sub !--Add tendencies cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i)) qcld(i) = MAX(0., qcld(i) + dqt_sub) qvc(i) = MAX(0., qvc(i) + dqvc_sub(i)) ENDIF ! qvapincld .LT. qvapincld_new ENDIF ! qiceincld .LT. eps ENDIF ! cldfra(i) .GT. eps !-------------------------------------------------------------------------- !-- CONDENSATION AND DIAGNOTICS OF SUB- AND SUPERSATURATED REGIONS -- !-------------------------------------------------------------------------- !--This section relies on a distribution of water in the clear-sky region of !--the mesh. !--If there is a clear-sky region IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN !--Water quantity in the clear-sky + potential liquid cloud (gridbox average) qclr = qtot(i) - qcld(i) !--New PDF rhl_clr = qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100. rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp) !--Calculation of the properties of the PDF !--Parameterization from IAGOS observations !--pdf_e1 and pdf_e2 will be reused below !--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 ELSE pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp ENDIF pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. ) pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3 pdf_e2 = GAMMA(1. + 1. / pdf_alpha) pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. )) pdf_loc = rhl_clr - pdf_scale * pdf_e2 !--Diagnostics of ratqs ratqs(i) = pdf_std / ( qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100. ) !--Calculation of the newly condensed water and fraction (pronostic) !--Integration of the clear sky PDF between gamma_cond*qsat and +inf !--NB. the calculated values are clear-sky averaged pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100. pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y ) pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2 cf_cond = EXP( - pdf_y ) qt_cond = ( pdf_e3 + pdf_loc * cf_cond ) * qsatl(i) / 100. IF ( ok_warm_cloud ) THEN !--If the statistical scheme is activated, the 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 dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond dqt_con = ( 1. - cldfra(i) ) * qt_cond !--Barriers dcf_con(i) = MIN(dcf_con(i), 1. - cldfra(i)) dqt_con = MIN(dqt_con, qclr) IF ( ok_unadjusted_clouds ) THEN !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute !--the new vapor qvapincld. The timestep is divided by two because we do not !--know when the condensation occurs qvapincld = gamma_cond(i) * qsat(i) qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i) !--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. ) ELSE !--We keep the saturation adjustment hypothesis, and the vapor in the !--newly formed cloud is set equal to the saturation vapor. qvapincld = qsat(i) ENDIF !--Tendency on cloud vapor and diagnostic dqvc_con(i) = qvapincld * dcf_con(i) dqi_con(i) = dqt_con - dqvc_con(i) !--Note that the tendencies are NOT added because they are !--added after the mixing process. In the following, the gridbox fraction is !-- 1. - dcf_con(i), and the total water in the gridbox is !-- qtot(i) - dqi_con(i) - dqvc_con(i) ENDIF ! ok_warm_cloud, cf_cond .GT. 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 pdf_x = qsat(i) / qsatl(i) * 100. pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y ) pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2 issrfra(i) = EXP( - pdf_y ) * ( 1. - dcf_con(i) - cldfra(i) ) qissr(i) = ( pdf_e3 * ( 1. - dcf_con(i) - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100. ENDIF !--Calculation of the subsaturated clear sky fraction and water subfra(i) = 1. - dcf_con(i) - cldfra(i) - issrfra(i) qsub(i) = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i) - qissr(i) !-------------------------------------- !-- CLOUD MIXING -- !-------------------------------------- !--This process mixes the cloud with its surroundings: the subsaturated clear sky, !--and the supersaturated clear sky. It is activated if the cloud is big enough, !--but does not cover the entire mesh. ! IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) & .AND. .NOT. ok_warm_cloud ) THEN !--Initialisation dcf_mix_sub = 0. dqt_mix_sub = 0. dqvc_mix_sub = 0. dcf_mix_issr = 0. dqt_mix_issr = 0. dqvc_mix_issr = 0. !-- PART 1 - TURBULENT DIFFUSION !--Clouds within the mesh are assumed to be ellipses. The length of the !--semi-major axis is a and the length of the semi-minor axis is b. !--N_cld_mix is the number of clouds within the mesh, and !--clouds_perim is the total perimeter of the clouds within the mesh, !--not considering interfaces with other meshes (only the interfaces with clear !--sky are taken into account). !-- !--The area of each cloud is A = a * b * RPI, !--and the perimeter of each cloud is !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) ) !-- !--With cell_area the area of the cell, we have: !-- cldfra = A * N_cld_mix / cell_area !-- clouds_perim = P * N_cld_mix !-- !--We assume that the ratio between b and a is a function of !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds !--are spherical. !-- b / a = bovera = MAX(0.1, cldfra) bovera = MAX(0.1, cldfra(i)) !--The clouds perimeter is imposed using the formula from Morcrette 2012, !--based on observations. !-- clouds_perim_normalized = alpha * cldfra * ( 1. - cldfra ) = clouds_perim / ( norm_constant * cell_area ) !-- !--With all this, we have !-- a = SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) ) !-- P = RPI * a * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) ) !--and therefore, using the perimeter !-- alpha * cldfra * ( 1. - cldfra ) * norm_constant * cell_area !-- = N_cld_mix * RPI & !-- * SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) ) & !-- * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) ) !--and finally N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - dcf_con(i) - cldfra(i) )**2. & * cell_area(i) * ( 1. - dcf_con(i) ) * bovera / RPI & / ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )**2. !--where coef_mix_lscp = ( alpha * norm_constant )**2. !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer !--In particular, it is 0 if cldfra = 1 a_mix = SQRT( cell_area(i) * ( 1. - dcf_con(i) ) * cldfra(i) / bovera / N_cld_mix / RPI ) !--The time required for turbulent diffusion to homogenize a region of size !--L_mix is defined as (L_mix**2./tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014) !--We compute L_mix and assume that the cloud is mixed over this length L_mix = SQRT( dtime**3. * pbl_eps(i) ) !--The mixing length cannot be greater than the semi-minor axis. In this case, !--the entire cloud is mixed. L_mix = MIN(L_mix, a_mix * bovera) !--The fraction of clear sky mixed is !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area envfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) & * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2. ) !--The fraction of cloudy sky mixed is !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area cldfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) & * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2. ) !-- PART 2 - SHEARING !--The clouds are then sheared. We keep the shape and number !--assumptions from before. The clouds are sheared along their !--semi-major axis (a_mix), on the entire cell heigh dz. !--The increase in size is L_shear = coef_shear_lscp * shear(i) * dz * dtime !--therefore, the total increase in fraction is !-- 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 environment and cloud mixed fractions are the same, !--which we add to the previous calculated mixed fractions. !--We therefore assume that the sheared clouds and the turbulent !--mixed clouds are different. envfra_mix = envfra_mix + shear_fra cldfra_mix = cldfra_mix + shear_fra !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES !--The environment fraction is allocated to subsaturated sky or supersaturated sky, !--according to the factor sigma_mix. This is computed as the ratio of the !--subsaturated sky fraction to the environment fraction, corrected by a factor !--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the !--supersaturated sky is favoured. Physically, this means that it is more likely !--to have supersaturated sky around the cloud than subsaturated sky. sigma_mix = subfra(i) / ( subfra(i) + chi_mixing_lscp * issrfra(i) ) subfra_mix = MIN( sigma_mix * envfra_mix, subfra(i) ) issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra(i) ) cldfra_mix = MIN( cldfra_mix, cldfra(i) ) !--First, we mix the subsaturated sky (subfra_mix) and the cloud close !--to this fraction (sigma_mix * cldfra_mix). IF ( subfra(i) .GT. eps ) THEN IF ( ok_unadjusted_clouds ) THEN !--The subsaturated air is simply added to the cloud, !--with the corresponding cloud fraction !--If the cloud is too subsaturated, the sublimation process !--activated in the following timestep will reduce the cloud !--fraction dcf_mix_sub = subfra_mix dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i) dqvc_mix_sub = dqt_mix_sub ELSE !--We compute the total humidity in the mixed air, which !--can be either sub- or supersaturated. qvapinmix = ( qsub(i) * subfra_mix / subfra(i) & + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) & / ( subfra_mix + cldfra_mix * sigma_mix ) IF ( qvapinmix .GT. qsat(i) ) THEN !--If the mixed air is supersaturated, we condense the subsaturated !--region which was mixed. dcf_mix_sub = subfra_mix dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i) dqvc_mix_sub = dcf_mix_sub * qsat(i) ELSE !--Else, we sublimate the cloud which was mixed. dcf_mix_sub = - sigma_mix * cldfra_mix dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i) dqvc_mix_sub = dcf_mix_sub * qsat(i) ENDIF ENDIF ! ok_unadjusted_clouds ENDIF ! subfra .GT. eps !--We then mix the supersaturated sky (issrfra_mix) and the cloud, !--for which the mixed air is always supersatured, therefore !--the cloud necessarily expands IF ( issrfra(i) .GT. eps ) THEN IF ( ok_unadjusted_clouds ) THEN !--The ice supersaturated air is simply added to the !--cloud, and supersaturated vapor will be deposited on the !--cloud ice crystals by the deposition process in the !--following timestep dcf_mix_issr = issrfra_mix dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i) dqvc_mix_issr = dqt_mix_issr ELSE !--In this case, the additionnal vapor condenses dcf_mix_issr = issrfra_mix dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i) dqvc_mix_issr = dcf_mix_issr * qsat(i) ENDIF ! ok_unadjusted_clouds ENDIF ! issrfra .GT. eps !--Sum up the tendencies from subsaturated sky and supersaturated sky dcf_mix(i) = dcf_mix_sub + dcf_mix_issr dqt_mix = dqt_mix_sub + dqt_mix_issr dqvc_mix(i) = dqvc_mix_sub + dqvc_mix_issr dqi_mix(i) = dqt_mix - dqvc_mix(i) !--Add tendencies issrfra(i) = MAX(0., issrfra(i) - dcf_mix_issr) qissr(i) = MAX(0., qissr(i) - dqt_mix_issr) cldfra(i) = MAX(0., MIN(1. - dcf_con(i), cldfra(i) + dcf_mix(i))) qcld(i) = MAX(0., MIN(qtot(i) - dqi_con(i) - dqvc_con(i), qcld(i) + dqt_mix)) qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i))) ENDIF ! ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) ) !--Finally, we add the tendencies of condensation cldfra(i) = MIN(1., cldfra(i) + dcf_con(i)) qcld(i) = MIN(qtot(i), qcld(i) + dqvc_con(i) + dqi_con(i)) qvc(i) = MIN(qcld(i), qvc(i) + dqvc_con(i)) !---------------------------------------- !-- CONTRAILS AND AVIATION -- !---------------------------------------- !--Add a source of cirrus from aviation contrails !IF ( ok_plane_contrail ) THEN ! dcf_avi(i) = 0. ! dqi_avi(i) = 0. ! dqvc_avi(i) = 0. ! ! TODO implement ok_unadjusted_clouds ! IF ( issrfra(i) .GT. eps ) THEN ! contrail_fra = MIN(1., flight_m(i,k) * dtime * contrail_cross_section / V_cell) ! dcf_avi(i) = issrfra(i) * contrail_fra ! dqt_avi = dcf_avi(i) * qissr(i) / issrfra(i) ! dqvc_avi(i) = qsat(i) * dcf_avi(i) ! ! !--Add tendencies ! cldfra(i) = cldfra(i) + dcf_avi(i) ! issrfra(i) = issrfra(i) - dcf_avi(i) ! qcld(i) = qcld(i) + dqt_avi ! qvc(i) = qvc(i) + dqvc_avi(i) ! qissr(i) = qissr(i) - dqt_avi ! !--Diagnostics ! dqi_avi(i) = dqt_avi - qsat(i) * dcf_avi(i) ! ENDIF ! dcf_avi(i) = dcf_avi(i) / dtime ! dqi_avi(i) = dqi_avi(i) / dtime ! dqvc_avi(i) = dqvc_avi(i) / dtime !ENDIF !------------------------------------------- !-- FINAL BARRIERS AND OUTPUTS -- !------------------------------------------- IF ( cldfra(i) .LT. eps ) THEN !--If the cloud is too small, it is sublimated. cldfra(i) = 0. qcld(i) = 0. qvc(i) = 0. qincld(i) = qsat(i) ELSE qincld(i) = qcld(i) / cldfra(i) ENDIF ! cldfra .LT. eps !--Diagnostics dcf_sub(i) = dcf_sub(i) / dtime dcf_con(i) = dcf_con(i) / dtime dcf_mix(i) = dcf_mix(i) / dtime dqi_adj(i) = dqi_adj(i) / dtime dqi_sub(i) = dqi_sub(i) / dtime dqi_con(i) = dqi_con(i) / dtime dqi_mix(i) = dqi_mix(i) / dtime dqvc_adj(i) = dqvc_adj(i) / dtime dqvc_sub(i) = dqvc_sub(i) / dtime dqvc_con(i) = dqvc_con(i) / dtime dqvc_mix(i) = dqvc_mix(i) / dtime ENDIF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds ENDIF ! end keepgoing ENDDO ! end loop on i END SUBROUTINE condensation_ice_supersat !********************************************************************************** END MODULE lmdz_lscp_condensation