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, std_subl_pdf_lscp, & mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, & std100_pdf_lscp, k0_pdf_lscp, kappa_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(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 ! ! 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 ! ! 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, Povera, 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 !--Initialisation issrfra(i) = 0. qissr(i) = 0. !--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. ! AB WARM CLOUD !cldfra(i) = 0. !qcld(i) = 0. !qvc(i) = 0. cldfra(i) = MAX(0., MIN(1., cf_seri(i))) qcld(i) = MAX(0., MIN(qtot(i), ql_seri(i) + qi_seri(i) + rvc_seri(i) * qtot(i))) qvc(i) = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i))) ok_warm_cloud = .TRUE. ELSE !--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. !--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) !------------------------------------------------------------------- !-- 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 ) ! AB WARM CLOUD !IF ( ok_unadjusted_clouds ) THEN IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), & pplay(i), dtime, qvapincld_new) IF ( qvapincld_new .EQ. 0. ) THEN !--If all the ice has been sublimated, we sublimate !--completely the cloud and do not activate the sublimation !--process !--Tendencies and diagnostics dcf_sub(i) = - cldfra(i) dqvc_sub(i) = - qvc(i) dqi_sub(i) = - ( qcld(i) - qvc(i) ) cldfra(i) = 0. qcld(i) = 0. qvc(i) = 0. ENDIF ELSE !--We keep the saturation adjustment hypothesis, and the vapor in the !--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 ) .OR. ( iflag_cloud_sublim_pdf .GE. 4 ) ) THEN !--Sublimation of the subsaturated cloud !--iflag_cloud_sublim_pdf selects the PDF of the ice water content !--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 .EQ. 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 ) ELSEIF ( iflag_cloud_sublim_pdf .EQ. 4 ) THEN !--Normal distribution around qvapincld + qiceincld with width !--constant in the RHi space pdf_y = ( qvapincld_new - qvapincld - qiceincld ) & / ( std_subl_pdf_lscp / 100. * qsat(i)) pdf_e1 = 0.5 * ( 1. + ERF( pdf_y / SQRT(2.) ) ) pdf_e2 = EXP( - pdf_y**2 / 2. ) / SQRT( 2. * RPI ) dcf_sub(i) = - cldfra(i) * pdf_e1 dqt_sub = - cldfra(i) * ( ( qvapincld + qiceincld ) * pdf_e1 & - std_subl_pdf_lscp / 100. * qsat(i) * pdf_e2 ) ENDIF !--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. !--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. ) !-- tuning coef * (cell area**0.25) * (function of temperature) pdf_e1 = beta_pdf_lscp * ( ( 1. - cldfra(i) ) * cell_area(i) )**0.25 & * MAX( temp(i) - temp_thresh_pdf_lscp, 0. ) IF ( rhl_clr .GT. 50. ) THEN pdf_std = ( pdf_e1 - std100_pdf_lscp ) * ( 100. - rhl_clr ) / 50. + std100_pdf_lscp ELSE pdf_std = pdf_e1 * rhl_clr / 50. ENDIF pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. ) pdf_alpha = EXP( rhl_clr / 100. ) * pdf_e3 pdf_alpha = MIN(10., pdf_alpha) IF ( ok_warm_cloud ) THEN !--If the statistical scheme is activated, the standard deviation is adapted !--to depend on the pressure level. It is multiplied by ratqs, so that near the !--surface std is almost 0, and upper than about 450 hPa the std is left untouched pdf_std = pdf_std * ratqs(i) ENDIF pdf_e2 = GAMMA(1. + 1. / pdf_alpha) pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2 )) pdf_loc = rhl_clr - pdf_scale * pdf_e2 !--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. ! AB WARM CLOUD !IF ( ok_warm_cloud ) THEN ! !--If the statistical scheme is activated, the calculated increase is equal ! !--to the cloud fraction, we assume saturation adjustment, and the ! !--condensation process stops ! cldfra(i) = cf_cond ! qcld(i) = qt_cond ! qvc(i) = cldfra(i) * qsat(i) !ELSEIF ( cf_cond .GT. eps ) THEN IF ( cf_cond .GT. eps ) THEN dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond dqt_con = ( 1. - cldfra(i) ) * qt_cond !--Barriers dcf_con(i) = MIN(dcf_con(i), 1. - cldfra(i)) dqt_con = MIN(dqt_con, qclr) ! AB WARM CLOUD !IF ( ok_unadjusted_clouds ) THEN IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute !--the new vapor qvapincld. The timestep is divided by two because we do not !--know when the condensation occurs qvapincld = gamma_cond(i) * qsat(i) qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i) CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), & pplay(i), dtime / 2., qvapincld_new) qvapincld = qvapincld_new ELSE !--We keep the saturation adjustment hypothesis, and the vapor in the !--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) !--Add tendencies cldfra(i) = MIN(1., cldfra(i) + dcf_con(i)) qcld(i) = MIN(qtot(i), qcld(i) + dqt_con) qvc(i) = MIN(qcld(i), qvc(i) + dqvc_con(i)) ENDIF ! ok_warm_cloud, cf_cond .GT. eps !--We then calculate the part that is greater than qsat !--and lower than gamma_cond * qsat, which is the ice supersaturated region pdf_x = qsat(i) / qsatl(i) * 100. pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y ) pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2 issrfra(i) = EXP( - pdf_y ) * ( 1. - cldfra(i) ) qissr(i) = ( pdf_e3 * ( 1. - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100. issrfra(i) = issrfra(i) - dcf_con(i) qissr(i) = qissr(i) - dqvc_con(i) - dqi_con(i) ENDIF ! ( 1. - cldfra(i) ) .GT. eps !--Calculation of the subsaturated clear sky fraction and water subfra(i) = 1. - cldfra(i) - issrfra(i) qsub(i) = qtot(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. ! ! AB WARM CLOUD !IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) & ! .AND. .NOT. ok_warm_cloud ) THEN IF ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) ) THEN !--Initialisation 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 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 envfra_mix = N_cld_mix * RPI / cell_area(i) & * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 ) !--The fraction of cloudy sky mixed is !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area cldfra_mix = N_cld_mix * RPI / cell_area(i) & * ( 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. 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 !--We compute the total humidity in the mixed air, which !--can be either sub- or supersaturated. qvapinmix = ( qsub(i) * subfra_mix / subfra(i) & + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) & / ( subfra_mix + cldfra_mix * sigma_mix ) ! AB WARM CLOUD !IF ( ok_unadjusted_clouds ) THEN IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * sigma_mix / cldfra(i) & / ( subfra_mix + cldfra_mix * sigma_mix ) CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), & pplay(i), dtime, qvapincld_new) IF ( qvapincld_new .EQ. 0. ) THEN !--If all the ice has been sublimated, we sublimate !--completely the cloud and do not activate the sublimation !--process !--Tendencies and diagnostics dcf_mix_sub = - sigma_mix * cldfra_mix dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i) dqvc_mix_sub = dcf_mix_sub * qvc(i) / cldfra(i) ELSE dcf_mix_sub = subfra_mix dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i) dqvc_mix_sub = dcf_mix_sub * qvapincld_new ENDIF ELSE 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 ! AB WARM CLOUD !IF ( ok_unadjusted_clouds ) THEN IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN qvapinmix = ( qissr(i) * issrfra_mix / issrfra(i) & + qcld(i) * cldfra_mix * (1. - sigma_mix) / cldfra(i) ) & / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) ) qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * (1. - sigma_mix) / cldfra(i) & / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) ) CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), & pplay(i), dtime, qvapincld_new) dcf_mix_issr = issrfra_mix dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i) dqvc_mix_issr = dcf_mix_issr * qvapincld_new 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., cldfra(i) + dcf_mix(i))) qcld(i) = MAX(0., MIN(qtot(i), qcld(i) + dqt_mix)) qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i))) ENDIF ! ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) ) !---------------------------------------- !-- 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 !********************************************************************************** SUBROUTINE deposition_sublimation( & qvapincld, qiceincld, temp, qsat, pplay, dtime, & qvapincld_new) USE lmdz_lscp_ini, ONLY: RV, RLSTT, RTT, RD, EPS_W USE lmdz_lscp_ini, ONLY: eps USE lmdz_lscp_ini, ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice REAL, INTENT(IN) :: qvapincld ! REAL, INTENT(IN) :: qiceincld ! REAL, INTENT(IN) :: temp ! REAL, INTENT(IN) :: qsat ! REAL, INTENT(IN) :: pplay ! REAL, INTENT(IN) :: dtime ! time step [s] REAL, INTENT(OUT):: qvapincld_new ! ! ! for unadjusted clouds REAL :: qice_ratio REAL :: pres_sat, rho, kappa REAL :: air_thermal_conduct, water_vapor_diff REAL :: iwc REAL :: iwc_log_inf100, iwc_log_sup100 REAL :: iwc_inf100, alpha_inf100, coef_inf100 REAL :: mu_sup100, sigma_sup100, coef_sup100 REAL :: Dm_ice, rm_ice !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment !--hypothesis is lost, and the vapor in the cloud is purely prognostic. ! !--The deposition equation is !-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) !--from Lohmann et al. (2016), where !--alpha is the deposition coefficient [-] !--mi is the mass of one ice crystal [kg] !--C is the capacitance of an ice crystal [m] !--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-] !--R_v is the specific gas constant for humid air [J/kg/K] !--T is the temperature [K] !--esi is the saturation pressure w.r.t. ice [Pa] !--Dv is the diffusivity of water vapor [m2/s] !--Ls is the specific latent heat of sublimation [J/kg/K] !--ka is the thermal conductivity of dry air [J/m/s/K] ! !--alpha is a coefficient to take into account the fact that during deposition, a water !--molecule cannot join the crystal from everywhere, it must do so that the crystal stays !--coherent (with the same structure). It has no impact for sublimation. !--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016)) !--during deposition, and alpha = 1. during sublimation. !--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus !-- C = capa_cond_cirrus * rm_ice ! !--We have qice = Nice * mi, where Nice is the ice crystal !--number concentration per kg of moist air !--HYPOTHESIS 1: the ice crystals are spherical, therefore !-- mi = 4/3 * pi * rm_ice**3 * rho_ice !--HYPOTHESIS 2: the ice crystals are monodisperse with the !--initial radius rm_ice_0. !--NB. this is notably different than the assumption !--of a distributed qice in the cloud made in the sublimation process; !--should it be consistent? ! !--As the deposition process does not create new ice crystals, !--and because we assume a same rm_ice value for all crystals !--therefore the sublimation process does not destroy ice crystals !--(or, in a limit case, it destroys all ice crystals), then !--Nice is a constant during the sublimation/deposition process. !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice ) ! !--The deposition equation then reads: !-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice !-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat & !-- / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) & !-- * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice ) !-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) & !-- *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 rm_ice_0**2 rho_ice ) !--and we have !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2 !-- dqi/dt = qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2 !--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1)) ! !--This system of equations can be resolved with an exact !--explicit numerical integration, having one variable resolved !--explicitly, the other exactly. The exactly resolved variable is !--the one decreasing, so it is qvc if the process is !--condensation, qi if it is sublimation. ! !--kappa is computed as an initialisation constant, as it depends only !--on temperature and other pre-computed values pres_sat = qsat / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay !--This formula for air thermal conductivity comes from Beard and Pruppacher (1971) air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184 !--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976) water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4 kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat & * ( RV * temp / water_vapor_diff / pres_sat & + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) ) !--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process !--Dry density [kg/m3] rho = pplay / temp / RD !--Here, the initial vapor in the cloud is qvapincld, and we compute !--the new vapor qvapincld_new !--rm_ice formula from McFarquhar and Heymsfield (1997) iwc = qiceincld * rho * 1e3 iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837) iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) ) iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) ) alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100 coef_inf100 = iwc_inf100 * alpha_inf100**3 / 120. mu_sup100 = ( 5.2 + 0.0013 * ( temp - RTT ) ) & + ( 0.026 - 1.2E-3 * ( temp - RTT ) ) * iwc_log_sup100 sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp - RTT ) ) & + ( 0.018 - 2.1E-4 * ( temp - RTT ) ) * iwc_log_sup100 coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3. * mu_sup100 + 4.5 * sigma_sup100**2 ) Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2 ) & * coef_sup100 ) / ( coef_inf100 + coef_sup100 ) rm_ice = Dm_ice / 2. * 1.E-6 IF ( qvapincld .GE. qsat ) THEN !--If the cloud is initially supersaturated !--Exact explicit integration (qvc exact, qice explicit) qvapincld_new = qsat + ( qvapincld - qsat ) & * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 ) ELSE !--If the cloud is initially subsaturated !--Exact explicit integration (qice exact, qvc explicit) !--The barrier is set so that the resulting vapor in cloud !--cannot be greater than qsat !--qice_ratio is the ratio between the new ice content and !--the old one, it is comprised between 0 and 1 qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2 * dtime * ( qsat - qvapincld ) ) IF ( qice_ratio .LT. 0. ) THEN !--The new vapor in cloud is set to 0 so that the !--sublimation process does not activate qvapincld_new = 0. ELSE !--Else, the sublimation process is activated with the !--diagnosed new cloud water vapor !--The new vapor in the cloud is increased with the !--sublimated ice qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**1.5 ) !--The new vapor in the cloud cannot be greater than qsat qvapincld_new = MIN(qvapincld_new, qsat) ENDIF ! qice_ratio .LT. 0. ENDIF ! qvapincld .GT. qsat END SUBROUTINE deposition_sublimation END MODULE lmdz_lscp_condensation