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, pplay, paprsdn, paprsup, cfcon, & cldfra_in, qvc_in, qliq_in, qice_in, shear, pbl_eps, cell_area, & temp, qtot_in, qsat, gamma_cond, ratqs, keepgoing, pt_pron_clds, & cldfra_above, icesed_flux, & cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, dcf_sed, & dqi_adj, dqi_sub, dqi_con, dqi_mix, dqi_sed, & dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, dqvc_sed, & contfra_in, perscontfra_in, qva_in, qia_in, flight_dist, flight_h2o, & contfra, perscontfra, qcont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, & dcfa_ini, dqia_ini, dqta_ini, dcfa_sub, dqia_sub, dqta_sub, & dcfa_cir, dqta_cir, dcfa_mix, dqia_mix, dqta_mix) !---------------------------------------------------------------------- ! 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: RLSTT, RTT, RD, RG, RV, RPI, EPS_W USE lmdz_lscp_ini, ONLY: eps, temp_nowater, ok_unadjusted_clouds USE lmdz_lscp_ini, ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice USE lmdz_lscp_ini, ONLY: N_ice_volume, corr_incld_depsub, nu_iwc_pdf_lscp USE lmdz_lscp_ini, ONLY: beta_pdf_lscp, temp_thresh_pdf_lscp USE lmdz_lscp_ini, ONLY: std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp USE lmdz_lscp_ini, ONLY: coef_mixing_lscp, coef_shear_lscp, incld_ice_thresh USE lmdz_lscp_ini, ONLY: ok_plane_contrail, aspect_ratio_contrails USE lmdz_lscp_ini, ONLY: coef_mixing_contrails, coef_shear_contrails USE lmdz_lscp_ini, ONLY: chi_mixing_contrails, linear_contrails_lifetime USE lmdz_aviation, ONLY: contrails_formation IMPLICIT NONE ! ! Input ! INTEGER, INTENT(IN) :: klon ! number of horizontal grid points REAL, INTENT(IN) :: dtime ! time step [s] ! 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) :: cfcon ! cloud fraction from deep convection [-] REAL, INTENT(IN) , DIMENSION(klon) :: cldfra_in ! cloud fraction [-] REAL, INTENT(IN) , DIMENSION(klon) :: qvc_in ! gridbox-mean water vapor in cloud [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: qliq_in ! specific liquid water content [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: qice_in ! specific ice water content [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: shear ! vertical shear [s-1] REAL, INTENT(IN) , DIMENSION(klon) :: pbl_eps ! TKE dissipation [m2/s3] REAL, INTENT(IN) , DIMENSION(klon) :: cell_area ! cell area [m2] REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K] REAL, INTENT(IN) , DIMENSION(klon) :: qtot_in ! 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 LOGICAL, INTENT(IN) , DIMENSION(klon) :: pt_pron_clds ! .TRUE. if clouds are prognostic in this mesh REAL, INTENT(IN) , DIMENSION(klon) :: cldfra_above ! cloud fraction IN THE LAYER ABOVE [-] REAL, INTENT(IN) , DIMENSION(klon) :: icesed_flux ! sedimentated ice flux [kg/s/m2] ! ! Input for aviation ! REAL, INTENT(IN) , DIMENSION(klon) :: contfra_in ! input linear contrails fraction [-] REAL, INTENT(IN) , DIMENSION(klon) :: perscontfra_in! input contrail induced cirrus fraction [-] REAL, INTENT(IN) , DIMENSION(klon) :: qva_in ! input linear contrails total specific humidity [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: qia_in ! input linear contrails total specific humidity [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: flight_dist ! aviation distance flown concentration [m/s/m3] REAL, INTENT(IN) , DIMENSION(klon) :: flight_h2o ! aviation emitted H2O concentration [kgH2O/s/m3] ! ! 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) :: dcf_sed ! cloud fraction tendency because of sedimentation [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_adj ! specific ice content tendency because of temperature adjustment [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_sub ! specific ice content tendency because of sublimation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_con ! specific ice content tendency because of condensation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_mix ! specific ice content tendency because of cloud mixing [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_sed ! specific ice content tendency because of sedimentation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_adj ! specific cloud water vapor tendency because of temperature adjustment [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_sub ! specific cloud water vapor tendency because of sublimation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_con ! specific cloud water vapor tendency because of condensation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_sed ! specific cloud water vapor tendency because of sedimentation [kg/kg/s] ! ! Diagnostics for aviation ! NB. idem for the INOUT ! REAL, INTENT(INOUT), DIMENSION(klon) :: contfra ! linear contrail fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: perscontfra ! linear contrail induced cirrus fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: qcont ! linear contrail specific humidity [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: Tcritcont ! critical temperature for contrail formation [K] REAL, INTENT(INOUT), DIMENSION(klon) :: qcritcont ! critical specific humidity for contrail formation [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: potcontfraP ! potential persistent contrail fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: potcontfraNP ! potential non-persistent contrail fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: dcfa_ini ! contrails cloud fraction tendency because of initial formation [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqia_ini ! contrails ice specific humidity tendency because of initial formation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqta_ini ! contrails total specific humidity tendency because of initial formation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dcfa_sub ! contrails cloud fraction tendency because of sublimation [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqia_sub ! contrails ice specific humidity tendency because of sublimation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqta_sub ! contrails total specific humidity tendency because of sublimation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dcfa_cir ! contrails cloud fraction tendency because of conversion in cirrus [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqta_cir ! contrails total specific humidity tendency because of conversion in cirrus [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dcfa_mix ! contrails cloud fraction tendency because of mixing [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqia_mix ! contrails ice specific humidity tendency because of mixing [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqta_mix ! contrails total specific humidity tendency because of mixing [kg/kg/s] ! ! Local ! INTEGER :: i LOGICAL :: ok_warm_cloud REAL, DIMENSION(klon) :: qcld, qzero REAL, DIMENSION(klon) :: clrfra, qclr ! ! for lognormal REAL :: pdf_std, pdf_k, pdf_delta REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2 ! ! for unadjusted clouds REAL :: qiceincld, qvapincld, qvapincld_new ! ! for deposition / sublimation REAL :: pres_sat, kappa_depsub, tauinv_depsub REAL :: air_thermal_conduct, water_vapor_diff ! ! for dissipation REAL :: pdf_shape ! ! for condensation REAL, DIMENSION(klon) :: qsatl, dqsatl REAL, DIMENSION(klon) :: pdf_alpha, pdf_scale, pdf_gamma REAL :: rhl_clr, pdf_loc REAL :: pdf_e3, pdf_x, pdf_y REAL :: dqt_con ! ! for sedimentation REAL, DIMENSION(klon) :: qice_sedim REAL :: clrfra_sed ! ! for mixing REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix REAL :: cldfra_mix, clrfra_mix, sigma_mix REAL :: L_shear, shear_fra REAL :: qiceinmix, qvapinmix_lim, qvapinclr_lim REAL :: pdf_fra_above_nuc, pdf_q_above_nuc REAL :: pdf_fra_above_lim, pdf_q_above_lim REAL :: pdf_fra_below_lim REAL :: mixed_fraction ! ! for cell properties REAL :: rho, rhodz, dz ! ! for contrails REAL :: perscontfra_ratio REAL :: contrails_conversion_factor 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 ( .NOT. pt_pron_clds(i) ) THEN pdf_std = ratqs(i) * qtot_in(i) pdf_k = -SQRT( LOG( 1. + (pdf_std / qtot_in(i))**2 ) ) pdf_delta = LOG( qtot_in(i) / ( gamma_cond(i) * qsat(i) ) ) pdf_a = pdf_delta / ( pdf_k * SQRT(2.) ) pdf_b = pdf_k / (2. * SQRT(2.)) 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_in(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 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. ok_warm_cloud = ( temp(i) .GT. temp_nowater ) !--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. - cfcon(i), cldfra_in(i))) qcld(i) = MAX(0., MIN(qtot_in(i), qliq_in(i) + qice_in(i) + qvc_in(i))) qvc(i) = MAX(0., MIN(qcld(i), qvc_in(i))) !--Initialise clear fraction properties clrfra(i) = MAX(0., MIN(1., ( 1. - cfcon(i) ) - cldfra(i))) qclr(i) = qtot_in(i) - qcld(i) 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. dcf_sed(i) = 0. dqi_sed(i) = 0. dqvc_sed(i) = 0. IF ( icesed_flux(i) .GT. 0. ) THEN !--If ice sedimentation is activated, the quantity of sedimentated ice was added !--to the total water vapor in the precipitation routine. Here we remove it !--(it will be reincluded later) qice_sedim(i) = icesed_flux(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime qclr(i) = qclr(i) - qice_sedim(i) ENDIF !--Initialisation of the cell properties !--Dry density [kg/m3] rho = pplay(i) / temp(i) / RD !--Dry air mass [kg/m2] rhodz = ( paprsdn(i) - paprsup(i) ) / RG !--Cell thickness [m] dz = rhodz / rho !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment !--hypothesis is lost, and the vapor in the cloud is purely prognostic. ! !--The deposition equation is !-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) !--from Lohmann et al. (2016), where !--alpha is the deposition coefficient [-] !--mi is the mass of one ice crystal [kg] !--C is the capacitance of an ice crystal [m] !--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-] !--R_v is the specific gas constant for humid air [J/kg/K] !--T is the temperature [K] !--esi is the saturation pressure w.r.t. ice [Pa] !--Dv is the diffusivity of water vapor [m2/s] !--Ls is the specific latent heat of sublimation [J/kg/K] !--ka is the thermal conductivity of dry air [J/m/s/K] ! !--alpha is a coefficient to take into account the fact that during deposition, a water !--molecule cannot join the crystal from everywhere, it must do so that the crystal stays !--coherent (with the same structure). It has no impact for sublimation. !--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016)) !--during deposition, and alpha = 1. during sublimation. !--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus !-- C = capa_cond_cirrus * rm_ice ! !--We have qice = Nice * mi, where Nice is the ice crystal !--number concentration per kg of moist air !--HYPOTHESIS 1: the ice crystals are spherical, therefore !-- mi = 4/3 * pi * rm_ice**3 * rho_ice !--HYPOTHESIS 2: the ice crystals concentration is constant in the cloud ! !--The equation in terms of q_ice is valide locally, and the local ice water content !--follows a Gamma distribution with a factor nu_iwc_pdf_lscp. Therefore, by !--integrating the local equation over the PDF (entire cloud), a correcting factor !--must be included, equal to !-- corr_incld_depsub = GAMMA(nu + 1/3) / GAMMA(nu) / nu**(1/3) !--NB. this is equal to about 0.9, hence the correction is not big !--NB. to lighten the calculated, corr_incld_depsub is calculated in lmdz_lscp_ini ! !--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 !--hence dmi = dqi ! !--The deposition equation then reads for qi the in-cloud ice water content: !-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat * corr_incld_depsub & !-- / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice !-- dqi/dt = alpha*4pi*capa_cond_cirrus*Nice*corr_incld_depsub & !-- / ( 4pi/3 N_ice rho_ice )**(1/3) & !-- / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) & !-- qi**(1/3) * (qvc - qsat) / qsat !--and we have !-- dqvc/dt = - alpha * kappa(T) * qi**(1/3) * (qvc - qsat) !-- dqi/dt = alpha * kappa(T) * qi**(1/3) * (qvc - qsat) ! !--This system of equations can be resolved with an exact !--explicit numerical integration, having one variable resolved !--explicitly, the other exactly. qvc is always the variable solved exactly. ! !--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 !--NB. the greater kappa_depsub, the more efficient is the !--deposition/sublimation process kappa_depsub = 4. * RPI * capa_cond_cirrus * N_ice_volume / rho * corr_incld_depsub & / qsat(i) / ( 4. / 3. * RPI * N_ice_volume / rho * rho_ice )**(1./3.) & / ( RV * temp(i) / water_vapor_diff / pres_sat & + RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) ) !--If contrails are activated IF ( ok_plane_contrail ) THEN contfra(i) = MAX(0., MIN(cldfra(i), contfra_in(i))) perscontfra(i) = MAX(0., MIN(cldfra(i), perscontfra_in(i))) qcont(i) = MAX(0., MIN(qcld(i), qva_in(i) + qia_in(i))) dcfa_ini(i) = 0. dqia_ini(i) = 0. dqta_ini(i) = 0. dcfa_sub(i) = 0. dqia_sub(i) = 0. dqta_sub(i) = 0. dcfa_cir(i) = 0. dqta_cir(i) = 0. dcfa_mix(i) = 0. dqia_mix(i) = 0. dqta_mix(i) = 0. ELSE contfra(i) = 0. perscontfra(i) = 0. qcont(i) = 0. ENDIF !---------------------------------------------------------------------- !-- SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CONTRAIL -- !---------------------------------------------------------------------- !--If there is a contrail IF ( contfra(i) .GT. eps ) THEN !--We remove contrails from the main class cldfra(i) = cldfra(i) - contfra(i) qcld(i) = qcld(i) - qcont(i) qvc(i) = qvc(i) - qsat(i) * contfra(i) !--The contrail is always adjusted to saturation qiceincld = ( qcont(i) / contfra(i) - qsat(i) ) !--If the ice water content is too low, the cloud is purely sublimated IF ( qiceincld .LT. eps ) THEN dcfa_sub(i) = - contfra(i) dqia_sub(i) = - qiceincld * contfra(i) dqta_sub(i) = - qcont(i) contfra(i) = 0. qcont(i) = 0. clrfra(i) = MIN(1., clrfra(i) - dcfa_sub(i)) qclr(i) = qclr(i) - dqta_sub(i) ENDIF ! qiceincld .LT. eps ENDIF ! contfra(i) .GT. eps !--If there is a contrail induced cirrus, we save the ratio perscontfra_ratio = perscontfra(i) / MAX(eps, cldfra(i)) !------------------------------------------------------------------- !-- 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. clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i)) qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i) !--Else, the cloud is adjusted and sublimated ELSE !--Cirrus clouds cannot have an in-cloud ice water content lower than !--incld_ice_thresh times the saturation IF ( qiceincld .LT. ( incld_ice_thresh * qsat(i) ) ) THEN dcf_sub(i) = ( qiceincld / ( incld_ice_thresh * qsat(i) ) - 1. ) * cldfra(i) dqvc_sub(i) = qvapincld * dcf_sub(i) cldfra(i) = cldfra(i) + dcf_sub(i) qcld(i) = qcld(i) + dqvc_sub(i) qvc(i) = qvc(i) + dqvc_sub(i) clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i)) qclr(i) = qclr(i) - dqvc_sub(i) ENDIF IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN IF ( qvapincld .GE. qsat(i) ) THEN !--If the cloud is initially supersaturated !--Exact explicit integration (qvc exact, qice explicit) tauinv_depsub = depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub ELSE !--If the cloud is initially subsaturated !--Exact explicit integration (qvc exact, qice explicit) !--Same but depo_coef_cirrus = 1 tauinv_depsub = qiceincld**(1./3.) * kappa_depsub ENDIF ! qvapincld .GT. qsat qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime * tauinv_depsub ) !--If all the ice is sublimated IF ( qvapincld_new .GE. ( qvapincld + qiceincld ) ) qvapincld_new = 0. ELSE !--We keep the saturation adjustment hypothesis, and the vapor in the !--cloud is set equal to the saturation vapor IF ( ( qvapincld + qiceincld ) .GT. qsat(i) ) THEN qvapincld_new = qsat(i) ELSE qvapincld_new = 0. ENDIF ENDIF ! ok_unadjusted_clouds !------------------------------------ !-- DISSIPATION OF THE CLOUD -- !------------------------------------ !--If the dissipation process must be activated IF ( qvapincld_new .GT. qvapincld ) THEN !--Gamma distribution starting at qvapincld pdf_shape = nu_iwc_pdf_lscp / qiceincld pdf_y = pdf_shape * ( qvapincld_new - qvapincld ) pdf_e1 = GAMMAINC ( nu_iwc_pdf_lscp , pdf_y ) pdf_e2 = GAMMAINC ( nu_iwc_pdf_lscp + 1. , pdf_y ) !--Tendencies and diagnostics dcf_sub(i) = - cldfra(i) * pdf_e1 dqi_sub(i) = - cldfra(i) * pdf_e2 / pdf_shape dqvc_sub(i) = dcf_sub(i) * qvapincld !--Add tendencies cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i)) qcld(i) = qcld(i) + dqvc_sub(i) + dqi_sub(i) qvc(i) = qvc(i) + dqvc_sub(i) clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i)) qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i) ENDIF ! qvapincld_new .GT. 0. !------------------------------------ !-- PHASE ADJUSTMENT -- !------------------------------------ IF ( qvapincld_new .EQ. 0. ) THEN !--If all the ice has been sublimated, we sublimate !--completely the cloud and do not activate the dissipation !--process !--Tendencies and diagnostics dcf_sub(i) = - cldfra(i) dqvc_sub(i) = - qvc(i) dqi_sub(i) = - ( qcld(i) - qvc(i) ) cldfra(i) = 0. qcld(i) = 0. qvc(i) = 0. clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i)) qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i) ELSE !--Adjustment of the IWC to the new vapor in cloud !--(this can be either positive or negative) dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) ) dqi_adj(i) = - dqvc_adj(i) !--Add tendencies !--The vapor in the cloud is updated, but not qcld as it is constant !--through this process, as well as cldfra which is unmodified qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i))) ENDIF ENDIF ! qiceincld .LT. eps ENDIF ! cldfra(i) .GT. eps !--If there is a contrail induced cirrus, we restore it perscontfra(i) = perscontfra_ratio * cldfra(i) !-------------------------------------------------------------------------- !-- 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 ( clrfra(i) .GT. eps ) THEN !--New PDF rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. rhl_clr = MAX(0., MIN(150., rhl_clr)) !--Calculation of the properties of the PDF !--Parameterization from IAGOS observations !--pdf_alpha, pdf_scale and pdf_gamma will be reused below !--Coefficient for standard deviation: !-- tuning coef * (clear sky area**0.25) * (function of temperature) pdf_e1 = beta_pdf_lscp * ( clrfra(i) * cell_area(i) )**0.25 & * MAX( temp(i) - temp_thresh_pdf_lscp, 0. ) IF ( rhl_clr .GT. 50. ) THEN pdf_std = ( pdf_e1 - std100_pdf_lscp ) * ( 100. - rhl_clr ) / 50. + std100_pdf_lscp ELSE pdf_std = pdf_e1 * rhl_clr / 50. ENDIF pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. ) pdf_alpha(i) = EXP( rhl_clr / 100. ) * pdf_e3 pdf_alpha(i) = MIN(10., pdf_alpha(i)) !--Avoid overflows !IF ( ok_warm_cloud ) THEN ! !--If the statistical scheme is activated, the standard deviation is adapted ! !--to depend on the pressure level. It is multiplied by ratqs, so that near the ! !--surface std is almost 0, and upper than about 450 hPa the std is left untouched ! pdf_std = pdf_std * ratqs(i) !ENDIF pdf_gamma(i) = GAMMA(1. + 1. / pdf_alpha(i)) !--Barrier to avoid overflows pdf_scale(i) = MAX(0.01, pdf_std / SQRT( & GAMMA(1. + 2. / pdf_alpha(i)) - pdf_gamma(i)**2 )) pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) !--Calculation of the newly condensed water and fraction (pronostic) !--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 = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows pdf_fra_above_nuc = 0. pdf_q_above_nuc = 0. ELSEIF ( pdf_y .LT. -10. ) THEN pdf_fra_above_nuc = 1. pdf_q_above_nuc = qclr(i) / clrfra(i) ELSE pdf_y = EXP( pdf_y ) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) pdf_fra_above_nuc = EXP( - pdf_y ) pdf_q_above_nuc = ( pdf_e3 + pdf_loc * pdf_fra_above_nuc ) * qsatl(i) / 100. ENDIF IF ( pdf_fra_above_nuc .GT. eps ) THEN dcf_con(i) = clrfra(i) * pdf_fra_above_nuc dqt_con = clrfra(i) * pdf_q_above_nuc !--Barriers (should be useless dcf_con(i) = MIN(dcf_con(i), clrfra(i)) dqt_con = MIN(dqt_con, qclr(i)) IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute !--the new vapor qvapincld. The timestep is divided by two because we do not !--know when the condensation occurs qvapincld = gamma_cond(i) * qsat(i) qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i) tauinv_depsub = depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) & * EXP( - dtime / 2. * tauinv_depsub ) ELSE !--We keep the saturation adjustment hypothesis, and the vapor in the !--newly formed cloud is set equal to the saturation vapor. qvapincld_new = qsat(i) ENDIF !--Tendency on cloud vapor and diagnostic dqvc_con(i) = qvapincld_new * dcf_con(i) dqi_con(i) = dqt_con - dqvc_con(i) !--Add tendencies cldfra(i) = MIN(1., cldfra(i) + dcf_con(i)) qcld(i) = qcld(i) + dqt_con qvc(i) = qvc(i) + dqvc_con(i) clrfra(i) = MAX(0., clrfra(i) - dcf_con(i)) qclr(i) = qclr(i) - dqt_con ENDIF ! pdf_fra_above_nuc .GT. eps ELSE !--Default value for the clear sky distribution: homogeneous distribution pdf_alpha(i) = 1. pdf_gamma(i) = 1. pdf_scale(i) = 0.01 ENDIF ! clrfra(i) .GT. eps !--If there is a contrail induced cirrus, we save the ratio perscontfra_ratio = perscontfra(i) / MAX(eps, cldfra(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) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN !-- 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 in contact with clear sky, and can be non-integer. !--In particular, it is 0 if cldfra = 1. !--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)) !--P / a is a function of b / a only, that we can calculate !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) ) Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) ) !--The clouds perimeter is imposed using the formula from Morcrette 2012, !--based on observations. !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra ) !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have: !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra ) !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) ) a_mix = Povera / coef_mixing_lscp / RPI / ( 1. - cldfra(i) ) / bovera !--and finally, !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a ) N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - cldfra(i) ) * cell_area(i) & / Povera / a_mix !--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 clrfra_mix = N_cld_mix * RPI / cell_area(i) & * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 ) !--The fraction of clear sky mixed is !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area cldfra_mix = N_cld_mix * RPI / cell_area(i) & * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 ) !-- 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 fraction of clear sky mixed is !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area !--and the fraction of cloud mixed is !-- N_cld_mix * ( (a * b) - (a - L_shear) * b ) * RPI / 2. / cell_area !--(note that they are equal) shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i) !--and the environment and cloud mixed fractions are the same, !--which we add to the previous calculated mixed fractions. !--We therefore assume that the sheared clouds and the turbulent !--mixed clouds are different. clrfra_mix = clrfra_mix + shear_fra cldfra_mix = cldfra_mix + shear_fra !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES clrfra_mix = MAX(eps, MIN(clrfra(i), clrfra_mix)) cldfra_mix = MAX(eps, MIN(cldfra(i), cldfra_mix)) !--We compute the limit vapor in clear sky where the mixed cloud could not !--survive if all the ice crystals were sublimated. Note that here we assume, !--for growth or reduction of the cloud, that the mixed cloud is adjusted !--to saturation, ie the vapor in the mixed cloud is qsat. This is only a !--diagnostic, and if the cloud size is increased, we add the new vapor to the !--cloud's vapor without condensing or sublimating ice crystals IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN qiceinmix = ( qcld(i) - qvc(i) ) / cldfra(i) / ( 1. + clrfra_mix / cldfra_mix ) tauinv_depsub = qiceinmix**(1./3.) * kappa_depsub qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime * tauinv_depsub ) ) qvapinclr_lim = qvapinmix_lim * ( 1. + cldfra_mix / clrfra_mix ) & - qvc(i) / cldfra(i) * cldfra_mix / clrfra_mix ELSE !--NB. if tau_depsub = 0 (ie tauinv_depsub = inf), we get the same result as above qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) & - qcld(i) / cldfra(i) * cldfra_mix / clrfra_mix ENDIF IF ( qvapinclr_lim .LT. 0. ) THEN !--Whatever we do, the cloud will increase in size dcf_mix(i) = clrfra_mix dqvc_mix(i) = clrfra_mix * qclr(i) / clrfra(i) ELSE !--We then calculate the clear sky part where the humidity is lower than !--qvapinclr_lim, and the part where it is higher than qvapinclr_lim !--This is the clear-sky PDF calculated in the condensation section. Note !--that if we are here, we necessarily went through the condensation part !--because the clear sky fraction can only be reduced by condensation. !--Thus the `pdf_xxx` variables are well defined. rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qvapinclr_lim / qsatl(i) * 100. pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows pdf_fra_above_lim = 0. pdf_q_above_lim = 0. ELSEIF ( pdf_y .LT. -10. ) THEN pdf_fra_above_lim = clrfra(i) pdf_q_above_lim = qclr(i) ELSE pdf_y = EXP( pdf_y ) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) pdf_q_above_lim = ( pdf_e3 * clrfra(i) & + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. ENDIF pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim !--sigma_mix is the ratio of the surroundings of the clouds where mixing !--increases the size of the cloud, to the total surroundings of the clouds. !--This implies that ( 1. - sigma_mix ) quantifies the ratio where mixing !--decreases the size of the clouds sigma_mix = pdf_fra_above_lim / ( pdf_fra_below_lim + pdf_fra_above_lim ) IF ( pdf_fra_above_lim .GT. eps ) THEN dcf_mix(i) = clrfra_mix * sigma_mix dqvc_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim ENDIF IF ( pdf_fra_below_lim .GT. eps ) THEN dcf_mix(i) = dcf_mix(i) - cldfra_mix * ( 1. - sigma_mix ) dqvc_mix(i) = dqvc_mix(i) - cldfra_mix * ( 1. - sigma_mix ) & * qvc(i) / cldfra(i) dqi_mix(i) = dqi_mix(i) - cldfra_mix * ( 1. - sigma_mix ) & * ( qcld(i) - qvc(i) ) / cldfra(i) ENDIF ENDIF ENDIF ! ( cldfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) !-------------------------------------- !-- CONTRAIL MIXING -- !-------------------------------------- IF ( ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN !-- 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 in contact with clear sky, and can be non-integer. !--In particular, it is 0 if cldfra = 1. !--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 = aspect_ratio_contrails !--P / a is a function of b / a only, that we can calculate !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) ) Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) ) !--The clouds perimeter is imposed using the formula from Morcrette 2012, !--based on observations. !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra ) !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have: !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra ) !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) ) a_mix = Povera / coef_mixing_contrails / RPI / ( 1. - contfra(i) ) / bovera !--and finally, !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a ) N_cld_mix = coef_mixing_contrails * contfra(i) * ( 1. - contfra(i) ) * cell_area(i) & / Povera / a_mix !--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 clrfra_mix = N_cld_mix * RPI / cell_area(i) & * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 ) !--The fraction of clear sky mixed is !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area cldfra_mix = N_cld_mix * RPI / cell_area(i) & * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 ) !-- PART 2 - SHEARING !--The clouds are then sheared. We keep the shape and number !--assumptions from before. The clouds are sheared with a random orientation !--of the wind, on average we assume that the wind and the semi-major axis !--make a 45 degrees angle. Moreover, the contrails only mix !--along their semi-minor axis (b), because it is easier to compute. !--With this, the clouds increase in size along b only, by a factor !--L_shear * SQRT(2.) / 2. (to account for the 45 degrees orientation of the wind) L_shear = coef_shear_contrails * shear(i) * dz * dtime !--therefore, the fraction of clear sky mixed is !-- N_cld_mix * ( a * (b + L_shear * SQRT(2.) / 2.) - a * b ) * RPI / 2. / cell_area !--and the fraction of cloud mixed is !-- N_cld_mix * ( a * b - a * (b - L_shear * SQRT(2.) / 2.) ) * RPI / 2. / cell_area !--(note that they are equal) shear_fra = RPI * L_shear * a_mix * SQRT(2.) / 2. / 2. * N_cld_mix / cell_area(i) !--and the environment and cloud mixed fractions are the same, !--which we add to the previous calculated mixed fractions. !--We therefore assume that the sheared clouds and the turbulent !--mixed clouds are different. clrfra_mix = clrfra_mix + shear_fra cldfra_mix = cldfra_mix + shear_fra !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES clrfra_mix = MAX(eps, MIN(clrfra(i), clrfra_mix)) cldfra_mix = MAX(eps, MIN(cldfra(i), cldfra_mix)) !--We compute the limit vapor in clear sky where the mixed cloud could not !--survive if all the ice crystals were sublimated. Note that here we assume, !--for growth or reduction of the cloud, that the mixed cloud is adjusted !--to saturation, ie the vapor in the mixed cloud is qsat. This is only a !--diagnostic, and if the cloud size is increased, we add the new vapor to the !--cloud's vapor without condensing or sublimating ice crystals qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) & - qcont(i) / contfra(i) * cldfra_mix / clrfra_mix IF ( qvapinclr_lim .LT. 0. ) THEN !--Whatever we do, the cloud will increase in size dcfa_mix(i) = clrfra_mix dqta_mix(i) = clrfra_mix * qclr(i) / clrfra(i) ELSE !--We then calculate the clear sky part where the humidity is lower than !--qvapinclr_lim, and the part where it is higher than qvapinclr_lim !--This is the clear-sky PDF calculated in the condensation section. Note !--that if we are here, we necessarily went through the condensation part !--because the clear sky fraction can only be reduced by condensation. !--Thus the `pdf_xxx` variables are well defined. rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qvapinclr_lim / qsatl(i) * 100. pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows pdf_fra_above_lim = 0. pdf_q_above_lim = 0. ELSEIF ( pdf_y .LT. -10. ) THEN pdf_fra_above_lim = clrfra(i) pdf_q_above_lim = qclr(i) ELSE pdf_y = EXP( pdf_y ) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) pdf_q_above_lim = ( pdf_e3 * clrfra(i) & + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. ENDIF pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim !--sigma_mix is the ratio of the surroundings of the clouds where mixing !--increases the size of the cloud, to the total surroundings of the clouds. !--This implies that ( 1. - sigma_mix ) quantifies the ratio where mixing !--decreases the size of the clouds !--For aviation, we increase the chance that the air surrounding contrails !--is moist. This is quantified with chi_mixing_contrails sigma_mix = chi_mixing_contrails * pdf_fra_above_lim & / ( pdf_fra_below_lim + chi_mixing_contrails * pdf_fra_above_lim ) IF ( pdf_fra_above_lim .GT. eps ) THEN dcfa_mix(i) = clrfra_mix * sigma_mix dqta_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim ENDIF IF ( pdf_fra_below_lim .GT. eps ) THEN qvapincld = qcont(i) / contfra(i) qiceincld = qvapincld - qsat(i) dcfa_mix(i) = dcfa_mix(i) - cldfra_mix * ( 1. - sigma_mix ) dqta_mix(i) = dqta_mix(i) - cldfra_mix * ( 1. - sigma_mix ) * qvapincld dqia_mix(i) = dqia_mix(i) - cldfra_mix * ( 1. - sigma_mix ) * qiceincld ENDIF ENDIF ENDIF ! ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) IF ( contfra(i) .GT. eps ) THEN !--We balance the mixing tendencies between the different cloud classes mixed_fraction = dcf_mix(i) + dcfa_mix(i) IF ( mixed_fraction .GT. clrfra(i) ) THEN mixed_fraction = clrfra(i) / mixed_fraction dcf_mix(i) = dcf_mix(i) * mixed_fraction dqvc_mix(i) = dqvc_mix(i) * mixed_fraction dqi_mix(i) = dqi_mix(i) * mixed_fraction dcfa_mix(i) = dcfa_mix(i) * mixed_fraction dqta_mix(i) = dqta_mix(i) * mixed_fraction ENDIF !--Add tendencies contfra(i) = contfra(i) + dcfa_mix(i) qcont(i) = qcont(i) + dqta_mix(i) clrfra(i) = clrfra(i) - dcfa_mix(i) qclr(i) = qclr(i) - dqta_mix(i) ENDIF !--Add tendencies cldfra(i) = cldfra(i) + dcf_mix(i) qcld(i) = qcld(i) + dqvc_mix(i) + dqi_mix(i) qvc(i) = qvc(i) + dqvc_mix(i) clrfra(i) = clrfra(i) - dcf_mix(i) qclr(i) = qclr(i) - dqvc_mix(i) - dqi_mix(i) !--If there is a contrail induced cirrus, we restore it perscontfra(i) = perscontfra_ratio * cldfra(i) !--------------------------------------- !-- ICE SEDIMENTATION -- !--------------------------------------- ! !--If ice supersaturation is activated, the cloud properties are prognostic. !--The falling ice is then partly considered a new cloud in the gridbox. !--BEWARE with this parameterisation, we can create a new cloud with a much !--different ice water content and water vapor content than the existing cloud !--(if it exists). This implies than unphysical fluxes of ice and vapor !--occur between the existing cloud and the newly formed cloud (from sedimentation). !--Note also that currently, the precipitation scheme assume a maximum !--random overlap, meaning that the new formed clouds will not be affected !--by vertical wind shear. ! IF ( icesed_flux(i) .GT. 0. ) THEN clrfra_sed = MIN(clrfra(i), cldfra_above(i) - cldfra(i)) IF ( ( clrfra_sed .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN qiceincld = qice_sedim(i) / cldfra_above(i) IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN tauinv_depsub = qiceincld**(1./3.) * kappa_depsub qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) ) ELSE qvapinclr_lim = qsat(i) - qiceincld ENDIF rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qvapinclr_lim / qsatl(i) * 100. pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows pdf_fra_above_lim = 0. pdf_q_above_lim = 0. ELSEIF ( pdf_y .LT. -10. ) THEN pdf_fra_above_lim = clrfra(i) pdf_q_above_lim = qclr(i) ELSE pdf_y = EXP( pdf_y ) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) pdf_q_above_lim = ( pdf_e3 * clrfra(i) & + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. ENDIF IF ( pdf_fra_above_lim .GT. eps ) THEN dcf_sed(i) = clrfra_sed * pdf_fra_above_lim / clrfra(i) dqvc_sed(i) = dcf_sed(i) * pdf_q_above_lim / pdf_fra_above_lim !--We include the sedimentated ice: dqi_sed(i) = qiceincld & !--We include the sedimentated ice: * ( dcf_sed(i) & !--the part that falls in the newly formed cloud and + cldfra(i) ) !--the part that falls in the existing cloud ENDIF ELSE dqi_sed(i) = qice_sedim(i) ENDIF ! ( clrfra_sed .GT. eps .AND. ( clrfra(i) .GT. eps ) !--Add tendencies cldfra(i) = MIN(1., cldfra(i) + dcf_sed(i)) qcld(i) = qcld(i) + dqvc_sed(i) + dqi_sed(i) qvc(i) = qvc(i) + dqvc_sed(i) clrfra(i) = MAX(0., clrfra(i) - dcf_sed(i)) !--We re-include sublimated sedimentated ice into clear sky water vapor qclr(i) = qclr(i) - dqvc_sed(i) + qice_sedim(i) - dqi_sed(i) ENDIF ! icesed_flux(i) .GT. 0. !--We put back contrails in the clouds class cldfra(i) = cldfra(i) + contfra(i) qcld(i) = qcld(i) + qcont(i) qvc(i) = qvc(i) + qsat(i) * contfra(i) !--Diagnose ISSRs IF ( clrfra(i) .GT. eps ) THEN rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qsat(i) / qsatl(i) * 100. pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows issrfra(i) = 0. qissr(i) = 0. ELSEIF ( pdf_y .LT. -10. ) THEN issrfra(i) = clrfra(i) qissr(i) = qclr(i) ELSE pdf_y = EXP( pdf_y ) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) issrfra(i) = EXP( - pdf_y ) * clrfra(i) qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc * issrfra(i) ) * qsatl(i) / 100. ENDIF ELSE issrfra(i) = 0. qissr(i) = 0. 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 dcf_sed(i) = dcf_sed(i) / dtime dqi_adj(i) = dqi_adj(i) / dtime dqi_sub(i) = dqi_sub(i) / dtime dqi_con(i) = dqi_con(i) / dtime dqi_mix(i) = dqi_mix(i) / dtime dqi_sed(i) = dqi_sed(i) / dtime dqvc_adj(i) = dqvc_adj(i) / dtime dqvc_sub(i) = dqvc_sub(i) / dtime dqvc_con(i) = dqvc_con(i) / dtime dqvc_mix(i) = dqvc_mix(i) / dtime dqvc_sed(i) = dqvc_sed(i) / dtime ENDIF ! pt_pron_clds(i) ENDIF ! end keepgoing ENDDO ! end loop on i !---------------------------------------- !-- CONTRAILS AND AVIATION -- !---------------------------------------- IF ( ok_plane_contrail ) THEN CALL contrails_formation( & klon, dtime, pplay, temp, qsat, qsatl, gamma_cond, & flight_dist, clrfra, qclr, pdf_scale, pdf_alpha, pdf_gamma, & keepgoing, pt_pron_clds, & Tcritcont, qcritcont, potcontfraP, potcontfraNP, & dcfa_ini, dqia_ini, dqta_ini) DO i = 1, klon IF ( keepgoing(i) .AND. pt_pron_clds(i) ) THEN !--Convert existing contrail fraction into "natural" cirrus cloud fraction IF ( cldfra(i) .GE. ( 1. - cfcon(i) - eps ) ) THEN contrails_conversion_factor = 1. ELSE contrails_conversion_factor = ( 1. - & !--First, the linear contrails are continuously degraded in induced cirrus EXP( - dtime / linear_contrails_lifetime ) & !--Second, if the cloud fills the entire gridbox, the linear contrails !--cannot exist. The exponent is set so that this only happens for !--very cloudy gridboxes * ( 1. - cldfra(i) / ( 1. - cfcon(i) ) )**0.1 ) ENDIF dcfa_cir(i) = - contrails_conversion_factor * contfra(i) dqta_cir(i) = - contrails_conversion_factor * qcont(i) !--Add tendencies issrfra(i) = MAX(0., issrfra(i) - dcfa_ini(i)) qissr(i) = MAX(0., qissr(i) - dqta_ini(i)) clrfra(i) = MAX(0., clrfra(i) - dcfa_ini(i)) qclr(i) = MAX(0., qclr(i) - dqta_ini(i)) cldfra(i) = MAX(0., MIN(1. - cfcon(i), cldfra(i) + dcfa_ini(i))) qcld(i) = MAX(0., MIN(qtot_in(i), qcld(i) + dqta_ini(i))) qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dcfa_ini(i) * qsat(i))) contfra(i) = MAX(0., MIN(cldfra(i), contfra(i) + dcfa_cir(i) + dcfa_ini(i))) qcont(i) = MAX(0., MIN(qcld(i), qcont(i) + dqta_cir(i) + dqta_ini(i))) perscontfra(i) = perscontfra(i) - dcfa_cir(i) !--Diagnostics dcfa_ini(i) = dcfa_ini(i) / dtime dqia_ini(i) = dqia_ini(i) / dtime dqta_ini(i) = dqta_ini(i) / dtime dcfa_sub(i) = dcfa_sub(i) / dtime dqia_sub(i) = dqta_sub(i) / dtime dqta_sub(i) = dqta_sub(i) / dtime dcfa_cir(i) = dcfa_cir(i) / dtime dqta_cir(i) = dqta_cir(i) / dtime dcfa_mix(i) = dcfa_mix(i) / dtime dqia_mix(i) = dqia_mix(i) / dtime dqta_mix(i) = dqta_mix(i) / dtime !------------------------------------------- !-- FINAL BARRIERS AND OUTPUTS -- !------------------------------------------- IF ( cldfra(i) .LT. eps ) THEN !--If the cloud is too small, it is sublimated. cldfra(i) = 0. contfra(i)= 0. perscontfra(i) = 0. qcld(i) = 0. qvc(i) = 0. qcont(i) = 0. qincld(i) = qsat(i) ELSE qincld(i) = qcld(i) / cldfra(i) ENDIF ! cldfra .LT. eps IF ( contfra(i) .LT. eps ) THEN contfra(i) = 0. qcont(i) = 0. ENDIF IF ( perscontfra(i) .LT. eps ) perscontfra(i) = 0. ENDIF ! keepgoing ENDDO ! loop on klon ENDIF ! ok_plane_contrail END SUBROUTINE condensation_ice_supersat !********************************************************************************** END MODULE lmdz_lscp_condensation