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 !********************************************************************************** !********************************************************************************** SUBROUTINE condensation_cloudth(klon, & & temp,qt,qt_th,frac_th,zpspsk,play,thetal_th, & & ratqs,sigma_qtherm,qsth,qsenv,qcloud,ctot,ctotth,ctot_vol, & & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) ! This routine computes the condensation of clouds in convective boundary layers ! with thermals assuming two separate distribution of the saturation deficit in ! the thermal plumes and in the environment ! It is based on the work of Arnaud Jam (Jam et al. 2013, BLM) ! Author : Etienne Vignon (LMDZ/CNRS) ! Date: February 2025 ! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam ! IMPORTANT NOTE: we assume iflag_cloudth_vert=7 !----------------------------------------------------------------------------------- use lmdz_lscp_ini, only: iflag_cloudth_vert,iflag_ratqs,iflag_cloudth_vert_noratqs use lmdz_lscp_ini, only: vert_alpha, vert_alpha_th ,sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin use lmdz_lscp_ini, only: RTT, RG, RPI, RD, RV, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th IMPLICIT NONE !------------------------------------------------------------------------------ ! Declarations !------------------------------------------------------------------------------ ! INPUT/OUTPUT INTEGER, INTENT(IN) :: klon REAL, DIMENSION(klon), INTENT(IN) :: temp ! Temperature (liquid temperature) in the mesh [K] : has seen evap of precip REAL, DIMENSION(klon), INTENT(IN) :: qt ! total water specific humidity in the mesh [kg/kg]: has seen evap of precip REAL, DIMENSION(klon), INTENT(IN) :: qt_th ! total water specific humidity in thermals [kg/kg]: has not seen evap of precip REAL, DIMENSION(klon), INTENT(IN) :: thetal_th ! Liquid potential temperature in thermals [K]: has not seen the evap of precip REAL, DIMENSION(klon), INTENT(IN) :: frac_th ! Fraction of the mesh covered by thermals [0-1] REAL, DIMENSION(klon), INTENT(IN) :: zpspsk ! Exner potential REAL, DIMENSION(klon), INTENT(IN) :: play ! Pressure of layers [Pa] REAL, DIMENSION(klon), INTENT(IN) :: ratqs ! Parameter that determines the width of the water distrib [-] REAL, DIMENSION(klon), INTENT(IN) :: sigma_qtherm ! Parameter determining the width of the distrib in thermals [-] REAL, DIMENSION(klon), INTENT(IN) :: qsth ! Saturation specific humidity in thermals REAL, DIMENSION(klon), INTENT(IN) :: qsenv ! Saturation specific humidity in environment REAL, DIMENSION(klon), INTENT(INOUT) :: ctot ! Cloud fraction [0-1] REAL, DIMENSION(klon), INTENT(INOUT) :: ctotth ! Cloud fraction [0-1] in thermals REAL, DIMENSION(klon), INTENT(INOUT) :: ctot_vol ! Volume cloud fraction [0-1] REAL, DIMENSION(klon), INTENT(INOUT) :: qcloud ! In cloud total water content [kg/kg] REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment ! LOCAL VARIABLES INTEGER itap,ind1,l,ig,iter,k INTEGER iflag_topthermals, niter REAL qcth(klon) REAL qcenv(klon) REAL qctot(klon) REAL cth(klon) REAL cenv(klon) REAL cth_vol(klon) REAL cenv_vol(klon) REAL qt_env(klon), thetal_env(klon) REAL sqrtpi,sqrt2,sqrt2pi REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2 REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth REAL zdelta,qsatbef,zcor REAL Tbefth(klon), Tbefenv(klon) REAL qlbef REAL dqsatenv(klon), dqsatth(klon) REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) REAL qincloud(klon) REAL alenvl, aenvl REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc !------------------------------------------------------------------------------ ! Initialisation !------------------------------------------------------------------------------ sqrt2pi=sqrt(2.*rpi) sqrt2=sqrt(2.) sqrtpi=sqrt(rpi) !------------------------------------------------------------------------------- ! Thermal fraction calculation and standard deviation of the distribution !------------------------------------------------------------------------------- ! initialisations and calculation of temperature, humidity and saturation specific humidity cloudth_senv(:) = 0. cloudth_sth(:) = 0. cloudth_sigmaenv(:) = 0. cloudth_sigmath(:) = 0. DO ind1=1,klon Tbefenv(ind1) = temp(ind1) thetal_env(ind1) = Tbefenv(ind1)/zpspsk(ind1) Tbefth(ind1) = thetal_th(ind1)*zpspsk(ind1) qt_env(ind1) = (qt(ind1)-frac_th(ind1)*qt_th(ind1))/(1.-frac_th(ind1)) !qt = a*qtth + (1-a)*qtenv ENDDO DO ind1=1,klon IF (frac_th(ind1).GT.min_frac_th_cld) THEN !Thermal and environnement ! Environment: alenv=(RD/RV*RLVTT*qsenv(ind1))/(rd*thetal_env(ind1)**2) aenv=1./(1.+(alenv*RLVTT/rcpd)) senv=aenv*(qt_env(ind1)-qsenv(ind1)) ! Thermals: alth=(RD/RV*RLVTT*qsth(ind1))/(rd*thetal_th(ind1)**2) ath=1./(1.+(alth*RLVTT/rcpd)) sth=ath*(qt_th(ind1)-qsth(ind1)) ! Standard deviation of the distributions ! environment sigma1s_fraca = (sigma1s_factor**0.5)*(frac_th(ind1)**sigma1s_power) / & & (1-frac_th(ind1))*((sth-senv)**2)**0.5 IF (cloudth_ratqsmin>0.) THEN sigma1s_ratqs = cloudth_ratqsmin*qt(ind1) ELSE sigma1s_ratqs = ratqs(ind1)*qt(ind1) ENDIF sigma1s = sigma1s_fraca + sigma1s_ratqs IF (iflag_ratqs.eq.10.or.iflag_ratqs.eq.11) then sigma1s = ratqs(ind1)*qt(ind1)*aenv ENDIF ! thermals sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((frac_th(ind1)+0.02)**sigma2s_power))+0.002*qt_th(ind1) IF (iflag_ratqs.eq.10.and.sigma_qtherm(ind1).ne.0) then sigma2s = sigma_qtherm(ind1)*ath ENDIF ! surface cloud fraction deltasenv=aenv*vert_alpha*sigma1s deltasth=ath*vert_alpha_th*sigma2s xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s) xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s) exp_xenv1 = exp(-1.*xenv1**2) exp_xenv2 = exp(-1.*xenv2**2) xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s) xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s) exp_xth1 = exp(-1.*xth1**2) exp_xth2 = exp(-1.*xth2**2) cth(ind1)=0.5*(1.-1.*erf(xth1)) cenv(ind1)=0.5*(1.-1.*erf(xenv1)) ctot(ind1)=frac_th(ind1)*cth(ind1)+(1.-1.*frac_th(ind1))*cenv(ind1) ctotth(ind1)=frac_th(ind1)*cth(ind1) !volume cloud fraction and condensed water !environnement IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2 IntJ_CF=0.5*(1.-1.*erf(xenv2)) IF (deltasenv .LT. 1.e-10) THEN qcenv(ind1)=IntJ cenv_vol(ind1)=IntJ_CF ELSE IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2) IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2) IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv) IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv) qcenv(ind1)=IntJ+IntI1+IntI2+IntI3 cenv_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF ENDIF !thermals IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 IntJ_CF=0.5*(1.-1.*erf(xth2)) IF (deltasth .LT. 1.e-10) THEN qcth(ind1)=IntJ cth_vol(ind1)=IntJ_CF ELSE IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) qcth(ind1)=IntJ+IntI1+IntI2+IntI3 cth_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF ENDIF ! total qctot(ind1)=frac_th(ind1)*qcth(ind1)+(1.-1.*frac_th(ind1))*qcenv(ind1) ctot_vol(ind1)=frac_th(ind1)*cth_vol(ind1)+(1.-1.*frac_th(ind1))*cenv_vol(ind1) IF (cenv(ind1).LT.min_neb_th.and.cth(ind1).LT.min_neb_th) THEN ctot(ind1)=0. ctot_vol(ind1)=0. qcloud(ind1)=qsenv(ind1) qincloud(ind1)=0. ELSE qincloud(ind1)=qctot(ind1)/ctot(ind1) !to prevent situations with cloud condensed water greater than available total water qincloud(ind1)=min(qincloud(ind1),qt(ind1)/ctot(ind1)) ! we assume that water vapor in cloud is qsenv qcloud(ind1)=qincloud(ind1)+qsenv(ind1) ENDIF ! Outputs used to check the PDFs cloudth_senv(ind1) = senv cloudth_sth(ind1) = sth cloudth_sigmaenv(ind1) = sigma1s cloudth_sigmath(ind1) = sigma2s ENDIF ! selection of grid points concerned by thermals ENDDO !loop on klon RETURN END SUBROUTINE condensation_cloudth !***************************************************************************************** !***************************************************************************************** ! pre-cmip7 routines are below and are becoming obsolete !***************************************************************************************** !***************************************************************************************** SUBROUTINE cloudth(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,t, & & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) use lmdz_lscp_ini, only: iflag_cloudth_vert,iflag_ratqs USE yomcst_mod_h USE yoethf_mod_h IMPLICIT NONE !=========================================================================== ! Auteur : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== INCLUDE "FCTTRE.h" INTEGER itap,ind1,ind2 INTEGER ngrid,klev,klon,l,ig real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv REAL ztv(ngrid,klev) REAL po(ngrid) REAL zqenv(ngrid) REAL zqta(ngrid,klev) REAL fraca(ngrid,klev+1) REAL zpspsk(ngrid,klev) REAL paprs(ngrid,klev+1) REAL pplay(ngrid,klev) REAL ztla(ngrid,klev) REAL zthl(ngrid,klev) REAL zqsatth(ngrid,klev) REAL zqsatenv(ngrid,klev) REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL ctot(ngrid,klev) REAL rneb(ngrid,klev) REAL t(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi REAL rdd,cppd,Lv REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,xth,xenv REAL Tbef,zdelta,qsatbef,zcor REAL qlbef REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) REAL zqs(ngrid), qcloud(ngrid) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Gestion de deux versions de cloudth !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (iflag_cloudth_vert.GE.1) THEN CALL cloudth_vert(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,t) RETURN ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !------------------------------------------------------------------------------- ! Initialisation des variables r?elles !------------------------------------------------------------------------------- sigma1(:,ind2)=0. sigma2(:,ind2)=0. qlth(:,ind2)=0. qlenv(:,ind2)=0. qltot(:,ind2)=0. rneb(:,ind2)=0. qcloud(:)=0. cth(:,ind2)=0. cenv(:,ind2)=0. ctot(:,ind2)=0. qsatmmussig1=0. qsatmmussig2=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2pi=sqrt(2.*pi) !------------------------------------------------------------------------------- ! Calcul de la fraction du thermique et des ?cart-types des distributions !------------------------------------------------------------------------------- do ind1=1,ngrid if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) ! zqenv(ind1)=po(ind1) Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) ath=1./(1.+(alth*Lv/cppd)) sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !------------------------------------------------------------------------------ ! Calcul des ?cart-types pour s !------------------------------------------------------------------------------ ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2) ! if (paprs(ind1,ind2).gt.90000) then ! ratqs(ind1,ind2)=0.002 ! else ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 ! endif sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) ! sigma1s=ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 !------------------------------------------------------------------------------ ! Calcul de l'eau condens?e et de la couverture nuageuse !------------------------------------------------------------------------------ sqrt2pi=sqrt(2.*pi) xth=sth/(sqrt(2.)*sigma2s) xenv=senv/(sqrt(2.)*sigma1s) cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2)) qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (ctot(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) endif else ! gaussienne environnement seule zqenv(ind1)=po(ind1) Tbef=t(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) sigma1s=ratqs(ind1,ind2)*zqenv(ind1) sqrt2pi=sqrt(2.*pi) xenv=senv/(sqrt(2.)*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) if (ctot(ind1,ind2).lt.1.e-3) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) endif endif enddo return ! end END SUBROUTINE cloudth !=========================================================================== SUBROUTINE cloudth_vert(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,t) !=========================================================================== ! Auteur : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== USE yoethf_mod_h use lmdz_lscp_ini, only: iflag_cloudth_vert, vert_alpha USE yomcst_mod_h IMPLICIT NONE INCLUDE "FCTTRE.h" INTEGER itap,ind1,ind2 INTEGER ngrid,klev,klon,l,ig REAL ztv(ngrid,klev) REAL po(ngrid) REAL zqenv(ngrid) REAL zqta(ngrid,klev) REAL fraca(ngrid,klev+1) REAL zpspsk(ngrid,klev) REAL paprs(ngrid,klev+1) REAL pplay(ngrid,klev) REAL ztla(ngrid,klev) REAL zthl(ngrid,klev) REAL zqsatth(ngrid,klev) REAL zqsatenv(ngrid,klev) REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL ctot(ngrid,klev) REAL rneb(ngrid,klev) REAL t(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi REAL rdd,cppd,Lv,sqrt2,sqrtpi REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,xth,xenv REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth REAL Tbef,zdelta,qsatbef,zcor REAL qlbef REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur ! Change the width of the PDF used for vertical subgrid scale heterogeneity ! (J Jouhaud, JL Dufresne, JB Madeleine) REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) REAL zqs(ngrid), qcloud(ngrid) !------------------------------------------------------------------------------ ! Initialisation des variables r?elles !------------------------------------------------------------------------------ sigma1(:,ind2)=0. sigma2(:,ind2)=0. qlth(:,ind2)=0. qlenv(:,ind2)=0. qltot(:,ind2)=0. rneb(:,ind2)=0. qcloud(:)=0. cth(:,ind2)=0. cenv(:,ind2)=0. ctot(:,ind2)=0. qsatmmussig1=0. qsatmmussig2=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2pi=sqrt(2.*pi) sqrt2=sqrt(2.) sqrtpi=sqrt(pi) !------------------------------------------------------------------------------- ! Calcul de la fraction du thermique et des ?cart-types des distributions !------------------------------------------------------------------------------- do ind1=1,ngrid if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) ! zqenv(ind1)=po(ind1) Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) ath=1./(1.+(alth*Lv/cppd)) sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !------------------------------------------------------------------------------ ! Calcul des ?cart-types pour s !------------------------------------------------------------------------------ sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) ! if (paprs(ind1,ind2).gt.90000) then ! ratqs(ind1,ind2)=0.002 ! else ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 ! endif ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) ! sigma1s=ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 !------------------------------------------------------------------------------ ! Calcul de l'eau condens?e et de la couverture nuageuse !------------------------------------------------------------------------------ sqrt2pi=sqrt(2.*pi) xth=sth/(sqrt(2.)*sigma2s) xenv=senv/(sqrt(2.)*sigma1s) cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2)) qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) IF (iflag_cloudth_vert == 1) THEN !------------------------------------------------------------------------------- ! Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs !------------------------------------------------------------------------------- ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) ! deltasenv=aenv*0.01*po(ind1) ! deltasth=ath*0.01*zqta(ind1,ind2) xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s) xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s) xth1=(sth-deltasth)/(sqrt(2.)*sigma2s) xth2=(sth+deltasth)/(sqrt(2.)*sigma2s) coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv) coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth) cth(ind1,ind2)=0.5*(1.+1.*erf(xth2)) cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1)) IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! qlenv(ind1,ind2)=IntJ ! print*, qlenv(ind1,ind2),'VERIF EAU' IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1)) ! IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2)) ! IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1)) IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2)) IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1)) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! qlth(ind1,ind2)=IntJ ! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2' qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) ELSE IF (iflag_cloudth_vert == 2) THEN !------------------------------------------------------------------------------- ! Version 3: Modification Jean Jouhaud. On condense a partir de -delta s !------------------------------------------------------------------------------- ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) ! deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) ! deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) deltasenv=aenv*vert_alpha*sigma1s deltasth=ath*vert_alpha*sigma2s xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s) xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s) xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s) xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s) ! coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv) ! coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth) cth(ind1,ind2)=0.5*(1.-1.*erf(xth1)) cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2) IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2)) ! IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) ! IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) ! IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! qlenv(ind1,ind2)=IntJ ! print*, qlenv(ind1,ind2),'VERIF EAU' IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2) IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2)) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! qlth(ind1,ind2)=IntJ ! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2' qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) ENDIF ! of if (iflag_cloudth_vert==1 or 2) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) ! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) & ! & +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) endif ! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif' else ! gaussienne environnement seule zqenv(ind1)=po(ind1) Tbef=t(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) sigma1s=ratqs(ind1,ind2)*zqenv(ind1) sqrt2pi=sqrt(2.*pi) xenv=senv/(sqrt(2.)*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) if (ctot(ind1,ind2).lt.1.e-3) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) endif endif enddo return ! end END SUBROUTINE cloudth_vert SUBROUTINE cloudth_v3(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,sigma_qtherm,zqs,t, & & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) use lmdz_lscp_ini, only: iflag_cloudth_vert USE yomcst_mod_h USE yoethf_mod_h IMPLICIT NONE !=========================================================================== ! Author : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== INCLUDE "FCTTRE.h" integer, intent(in) :: ind2 integer, intent(in) :: ngrid,klev real, dimension(ngrid,klev), intent(in) :: ztv real, dimension(ngrid), intent(in) :: po real, dimension(ngrid,klev), intent(in) :: zqta real, dimension(ngrid,klev+1), intent(in) :: fraca real, dimension(ngrid), intent(out) :: qcloud real, dimension(ngrid,klev), intent(out) :: ctot real, dimension(ngrid,klev), intent(out) :: ctot_vol real, dimension(ngrid,klev), intent(in) :: zpspsk real, dimension(ngrid,klev+1), intent(in) :: paprs real, dimension(ngrid,klev), intent(in) :: pplay real, dimension(ngrid,klev), intent(in) :: ztla real, dimension(ngrid,klev), intent(inout) :: zthl real, dimension(ngrid,klev), intent(in) :: ratqs,sigma_qtherm real, dimension(ngrid), intent(in) :: zqs real, dimension(ngrid,klev), intent(in) :: t real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv REAL zqenv(ngrid) REAL zqsatth(ngrid,klev) REAL zqsatenv(ngrid,klev) REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL cth_vol(ngrid,klev) REAL cenv_vol(ngrid,klev) REAL rneb(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi REAL rdd,cppd,Lv REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2 REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks REAL Tbef,zdelta,qsatbef,zcor REAL qlbef REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) INTEGER :: ind1,l, ig IF (iflag_cloudth_vert.GE.1) THEN CALL cloudth_vert_v3(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,sigma_qtherm,zqs,t, & & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) RETURN ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !------------------------------------------------------------------------------- ! Initialisation des variables r?elles !------------------------------------------------------------------------------- sigma1(:,ind2)=0. sigma2(:,ind2)=0. qlth(:,ind2)=0. qlenv(:,ind2)=0. qltot(:,ind2)=0. rneb(:,ind2)=0. qcloud(:)=0. cth(:,ind2)=0. cenv(:,ind2)=0. ctot(:,ind2)=0. cth_vol(:,ind2)=0. cenv_vol(:,ind2)=0. ctot_vol(:,ind2)=0. qsatmmussig1=0. qsatmmussig2=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2pi=sqrt(2.*pi) sqrt2=sqrt(2.) sqrtpi=sqrt(pi) !------------------------------------------------------------------------------- ! Cloud fraction in the thermals and standard deviation of the PDFs !------------------------------------------------------------------------------- do ind1=1,ngrid if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 !po = qt de l'environnement ET des thermique !zqenv = qt environnement !zqsatenv = qsat environnement !zthl = Tl environnement Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 ath=1./(1.+(alth*Lv/cppd)) !al, p84 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 !zqta = qt thermals !zqsatth = qsat thermals !ztla = Tl thermals !------------------------------------------------------------------------------ ! s standard deviations !------------------------------------------------------------------------------ ! tests ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) ! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1) ! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2) ! final option sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) !------------------------------------------------------------------------------ ! Condensed water and cloud cover !------------------------------------------------------------------------------ xth=sth/(sqrt2*sigma2s) xenv=senv/(sqrt2*sigma1s) cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) ctot_vol(ind1,ind2)=ctot(ind1,ind2) qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2)) qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) if (ctot(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) endif else ! Environnement only, follow the if l.110 zqenv(ind1)=po(ind1) Tbef=t(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) sigma1s=ratqs(ind1,ind2)*zqenv(ind1) xenv=senv/(sqrt2*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot_vol(ind1,ind2)=ctot(ind1,ind2) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2)) if (ctot(ind1,ind2).lt.1.e-3) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) endif endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183 enddo ! from the loop on ngrid l.108 return ! end END SUBROUTINE cloudth_v3 !=========================================================================== SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,sigma_qtherm,zqs,t, & & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) !=========================================================================== ! Auteur : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== use yoethf_mod_h use lmdz_lscp_ini, only : iflag_cloudth_vert,iflag_ratqs use lmdz_lscp_ini, only : vert_alpha,vert_alpha_th, sigma1s_factor, sigma1s_power , sigma2s_factor , sigma2s_power , cloudth_ratqsmin , iflag_cloudth_vert_noratqs USE yomcst_mod_h IMPLICIT NONE INCLUDE "FCTTRE.h" INTEGER itap,ind1,ind2 INTEGER ngrid,klev,klon,l,ig real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv REAL ztv(ngrid,klev) REAL po(ngrid) REAL zqenv(ngrid) REAL zqta(ngrid,klev) REAL fraca(ngrid,klev+1) REAL zpspsk(ngrid,klev) REAL paprs(ngrid,klev+1) REAL pplay(ngrid,klev) REAL ztla(ngrid,klev) REAL zthl(ngrid,klev) REAL zqsatth(ngrid,klev) REAL zqsatenv(ngrid,klev) REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL ctot(ngrid,klev) REAL cth_vol(ngrid,klev) REAL cenv_vol(ngrid,klev) REAL ctot_vol(ngrid,klev) REAL rneb(ngrid,klev) REAL t(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi REAL rdd,cppd,Lv REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2 REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth REAL Tbef,zdelta,qsatbef,zcor REAL qlbef REAL ratqs(ngrid,klev),sigma_qtherm(ngrid,klev) ! determine la largeur de distribution de vapeur ! Change the width of the PDF used for vertical subgrid scale heterogeneity ! (J Jouhaud, JL Dufresne, JB Madeleine) REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) REAL zqs(ngrid), qcloud(ngrid) REAL rhodz(ngrid,klev) REAL zrho(ngrid,klev) REAL dz(ngrid,klev) DO ind1 = 1, ngrid !Layer calculation rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2 zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3 dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre END DO !------------------------------------------------------------------------------ ! Initialize !------------------------------------------------------------------------------ sigma1(:,ind2)=0. sigma2(:,ind2)=0. qlth(:,ind2)=0. qlenv(:,ind2)=0. qltot(:,ind2)=0. rneb(:,ind2)=0. qcloud(:)=0. cth(:,ind2)=0. cenv(:,ind2)=0. ctot(:,ind2)=0. cth_vol(:,ind2)=0. cenv_vol(:,ind2)=0. ctot_vol(:,ind2)=0. qsatmmussig1=0. qsatmmussig2=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2pi=sqrt(2.*pi) sqrt2=sqrt(2.) sqrtpi=sqrt(pi) !------------------------------------------------------------------------------- ! Calcul de la fraction du thermique et des ecart-types des distributions !------------------------------------------------------------------------------- do ind1=1,ngrid if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 !zqenv = qt environnement !zqsatenv = qsat environnement !zthl = Tl environnement Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 ath=1./(1.+(alth*Lv/cppd)) !al, p84 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 !zqta = qt thermals !zqsatth = qsat thermals !ztla = Tl thermals !------------------------------------------------------------------------------ ! s standard deviation !------------------------------------------------------------------------------ sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / & & (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5 ! sigma1s_fraca = (1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5 IF (cloudth_ratqsmin>0.) THEN sigma1s_ratqs = cloudth_ratqsmin*po(ind1) ELSE sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1) ENDIF sigma1s = sigma1s_fraca + sigma1s_ratqs sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2) IF (iflag_ratqs.eq.10.or.iflag_ratqs.eq.11) then sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv IF (iflag_ratqs.eq.10.and.sigma_qtherm(ind1,ind2).ne.0) then sigma2s = sigma_qtherm(ind1,ind2)*ath ENDIF ENDIF ! tests ! sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) ! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1) ! sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) ! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2) ! if (paprs(ind1,ind2).gt.90000) then ! ratqs(ind1,ind2)=0.002 ! else ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 ! endif ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) ! sigma1s=ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 IF (iflag_cloudth_vert == 1) THEN !------------------------------------------------------------------------------- ! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs !------------------------------------------------------------------------------- deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s) xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s) xth1=(sth-deltasth)/(sqrt(2.)*sigma2s) xth2=(sth+deltasth)/(sqrt(2.)*sigma2s) coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv) coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth) cth(ind1,ind2)=0.5*(1.+1.*erf(xth2)) cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) ! Environment IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1)) IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! Thermal IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1)) IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2)) IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1)) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) ELSE IF (iflag_cloudth_vert >= 3) THEN IF (iflag_cloudth_vert < 5) THEN !------------------------------------------------------------------------------- ! Version 3: Changes by J. Jouhaud; condensation for q > -delta s !------------------------------------------------------------------------------- ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) ! deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) ! deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) IF (iflag_cloudth_vert == 3) THEN deltasenv=aenv*vert_alpha*sigma1s deltasth=ath*vert_alpha_th*sigma2s ELSE IF (iflag_cloudth_vert == 4) THEN IF (iflag_cloudth_vert_noratqs == 1) THEN deltasenv=vert_alpha*max(sigma1s_fraca,1e-10) deltasth=vert_alpha_th*sigma2s ELSE deltasenv=vert_alpha*sigma1s deltasth=vert_alpha_th*sigma2s ENDIF ENDIF xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s) xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s) exp_xenv1 = exp(-1.*xenv1**2) exp_xenv2 = exp(-1.*xenv2**2) xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s) xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s) exp_xth1 = exp(-1.*xth1**2) exp_xth2 = exp(-1.*xth2**2) !CF_surfacique cth(ind1,ind2)=0.5*(1.-1.*erf(xth1)) cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) !CF_volumique & eau condense !environnement IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2 IntJ_CF=0.5*(1.-1.*erf(xenv2)) if (deltasenv .lt. 1.e-10) then qlenv(ind1,ind2)=IntJ cenv_vol(ind1,ind2)=IntJ_CF else IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2) IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2) IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv) IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv) qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF endif !thermique IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 IntJ_CF=0.5*(1.-1.*erf(xth2)) if (deltasth .lt. 1.e-10) then qlth(ind1,ind2)=IntJ cth_vol(ind1,ind2)=IntJ_CF else IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF endif qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) ELSE IF (iflag_cloudth_vert == 5) THEN sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) & /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) & +ratqs(ind1,ind2)*po(ind1) !Environment sigma2s=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.02)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2) !Thermals !sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) !sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) xth=sth/(sqrt(2.)*sigma2s) xenv=senv/(sqrt(2.)*sigma1s) !Volumique cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth)) cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) !print *,'jeanjean_CV=',ctot_vol(ind1,ind2) qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth_vol(ind1,ind2)) qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv_vol(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) !Surfacique !Neggers !beta=0.0044 !inverse_rho=1.+beta*dz(ind1,ind2) !print *,'jeanjean : beta=',beta !cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho !cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho !ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) !Brooks a_Brooks=0.6694 b_Brooks=0.1882 A_Maj_Brooks=0.1635 !-- sans shear !A_Maj_Brooks=0.17 !-- ARM LES !A_Maj_Brooks=0.18 !-- RICO LES !A_Maj_Brooks=0.19 !-- BOMEX LES Dx_Brooks=200000. f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks)) !print *,'jeanjean_f=',f_Brooks cth(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.)) cenv(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.)) ctot(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.)) !print *,'JJ_ctot_1',ctot(ind1,ind2) ENDIF ! of if (iflag_cloudth_vert<5) ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4) ! if (ctot(ind1,ind2).lt.1.e-10) then if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. ctot_vol(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) ! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) & ! & +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) endif else ! gaussienne environnement seule zqenv(ind1)=po(ind1) Tbef=t(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) sth=0. sigma1s=ratqs(ind1,ind2)*zqenv(ind1) sigma2s=0. sqrt2pi=sqrt(2.*pi) xenv=senv/(sqrt(2.)*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot_vol(ind1,ind2)=ctot(ind1,ind2) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) if (ctot(ind1,ind2).lt.1.e-3) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ! ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) endif endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492 ! Outputs used to check the PDFs cloudth_senv(ind1,ind2) = senv cloudth_sth(ind1,ind2) = sth cloudth_sigmaenv(ind1,ind2) = sigma1s cloudth_sigmath(ind1,ind2) = sigma2s enddo ! from the loop on ngrid l.333 return ! end END SUBROUTINE cloudth_vert_v3 ! SUBROUTINE cloudth_v6(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,T, & & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) USE yoethf_mod_h USE lmdz_lscp_ini, only: iflag_cloudth_vert USE yomcst_mod_h IMPLICIT NONE INCLUDE "FCTTRE.h" !Domain variables INTEGER ngrid !indice Max lat-lon INTEGER klev !indice Max alt real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv INTEGER ind1 !indice in [1:ngrid] INTEGER ind2 !indice in [1:klev] !thermal plume fraction REAL fraca(ngrid,klev+1) !thermal plumes fraction in the gridbox !temperatures REAL T(ngrid,klev) !temperature REAL zpspsk(ngrid,klev) !factor (p/p0)**kappa (used for potential variables) REAL ztv(ngrid,klev) !potential temperature (voir thermcell_env.F90) REAL ztla(ngrid,klev) !liquid temperature in the thermals (Tl_th) REAL zthl(ngrid,klev) !liquid temperature in the environment (Tl_env) !pressure REAL paprs(ngrid,klev+1) !pressure at the interface of levels REAL pplay(ngrid,klev) !pressure at the middle of the level !humidity REAL ratqs(ngrid,klev) !width of the total water subgrid-scale distribution REAL po(ngrid) !total water (qt) REAL zqenv(ngrid) !total water in the environment (qt_env) REAL zqta(ngrid,klev) !total water in the thermals (qt_th) REAL zqsatth(ngrid,klev) !water saturation level in the thermals (q_sat_th) REAL zqsatenv(ngrid,klev) !water saturation level in the environment (q_sat_env) REAL qlth(ngrid,klev) !condensed water in the thermals REAL qlenv(ngrid,klev) !condensed water in the environment REAL qltot(ngrid,klev) !condensed water in the gridbox !cloud fractions REAL cth_vol(ngrid,klev) !cloud fraction by volume in the thermals REAL cenv_vol(ngrid,klev) !cloud fraction by volume in the environment REAL ctot_vol(ngrid,klev) !cloud fraction by volume in the gridbox REAL cth_surf(ngrid,klev) !cloud fraction by surface in the thermals REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox !PDF of saturation deficit variables REAL rdd,cppd,Lv REAL Tbef,zdelta,qsatbef,zcor REAL alth,alenv,ath,aenv REAL sth,senv !saturation deficits in the thermals and environment REAL sigma_env,sigma_th !standard deviations of the biGaussian PDF !cloud fraction variables REAL xth,xenv REAL inverse_rho,beta !Neggers et al. (2011) method REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method !Incloud total water variables REAL zqs(ngrid) !q_sat REAL qcloud(ngrid) !eau totale dans le nuage !Some arithmetic variables REAL pi,sqrt2,sqrt2pi !Depth of the layer REAL dz(ngrid,klev) !epaisseur de la couche en metre REAL rhodz(ngrid,klev) REAL zrho(ngrid,klev) DO ind1 = 1, ngrid rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2] zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd ![kg/m3] dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) ![m] END DO !------------------------------------------------------------------------------ ! Initialization !------------------------------------------------------------------------------ qlth(:,ind2)=0. qlenv(:,ind2)=0. qltot(:,ind2)=0. cth_vol(:,ind2)=0. cenv_vol(:,ind2)=0. ctot_vol(:,ind2)=0. cth_surf(:,ind2)=0. cenv_surf(:,ind2)=0. ctot_surf(:,ind2)=0. qcloud(:)=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2=sqrt(2.) sqrt2pi=sqrt(2.*pi) DO ind1=1,ngrid !------------------------------------------------------------------------------- !Both thermal and environment in the gridbox !------------------------------------------------------------------------------- IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN !-------------------------------------------- !calcul de qsat_env !-------------------------------------------- Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef !-------------------------------------------- !calcul de s_env !-------------------------------------------- alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 these Arnaud Jam aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam !-------------------------------------------- !calcul de qsat_th !-------------------------------------------- Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef !-------------------------------------------- !calcul de s_th !-------------------------------------------- alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 these Arnaud Jam ath=1./(1.+(alth*Lv/cppd)) !al, p84 these Arnaud Jam sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 these Arnaud Jam !-------------------------------------------- !calcul standard deviations bi-Gaussian PDF !-------------------------------------------- sigma_th=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.01)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2) sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) & /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) & +ratqs(ind1,ind2)*po(ind1) xth=sth/(sqrt2*sigma_th) xenv=senv/(sqrt2*sigma_env) !-------------------------------------------- !Cloud fraction by volume CF_vol !-------------------------------------------- cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth)) cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) !-------------------------------------------- !Condensed water qc !-------------------------------------------- qlth(ind1,ind2)=sigma_th*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth_vol(ind1,ind2)) qlenv(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv_vol(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) !-------------------------------------------- !Cloud fraction by surface CF_surf !-------------------------------------------- !Method Neggers et al. (2011) : ok for cumulus clouds only !beta=0.0044 (Jouhaud et al.2018) !inverse_rho=1.+beta*dz(ind1,ind2) !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho !Method Brooks et al. (2005) : ok for all types of clouds a_Brooks=0.6694 b_Brooks=0.1882 A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent Dx_Brooks=200000. !-- si l'on considere des mailles de 200km de cote f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks)) ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.)) !-------------------------------------------- !Incloud Condensed water qcloud !-------------------------------------------- if (ctot_surf(ind1,ind2) .lt. 1.e-10) then ctot_vol(ind1,ind2)=0. ctot_surf(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1) endif !------------------------------------------------------------------------------- !Environment only in the gridbox !------------------------------------------------------------------------------- ELSE !-------------------------------------------- !calcul de qsat_env !-------------------------------------------- Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef !-------------------------------------------- !calcul de s_env !-------------------------------------------- alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 these Arnaud Jam aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam !-------------------------------------------- !calcul standard deviations Gaussian PDF !-------------------------------------------- zqenv(ind1)=po(ind1) sigma_env=ratqs(ind1,ind2)*zqenv(ind1) xenv=senv/(sqrt2*sigma_env) !-------------------------------------------- !Cloud fraction by volume CF_vol !-------------------------------------------- ctot_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv)) !-------------------------------------------- !Condensed water qc !-------------------------------------------- qltot(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*ctot_vol(ind1,ind2)) !-------------------------------------------- !Cloud fraction by surface CF_surf !-------------------------------------------- !Method Neggers et al. (2011) : ok for cumulus clouds only !beta=0.0044 (Jouhaud et al.2018) !inverse_rho=1.+beta*dz(ind1,ind2) !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho !Method Brooks et al. (2005) : ok for all types of clouds a_Brooks=0.6694 b_Brooks=0.1882 A_Maj_Brooks=0.1635 !-- sans dependence au shear Dx_Brooks=200000. f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks)) ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.)) !-------------------------------------------- !Incloud Condensed water qcloud !-------------------------------------------- if (ctot_surf(ind1,ind2) .lt. 1.e-8) then ctot_vol(ind1,ind2)=0. ctot_surf(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2) endif END IF ! From the separation (thermal/envrionnement) et (environnement only) ! Outputs used to check the PDFs cloudth_senv(ind1,ind2) = senv cloudth_sth(ind1,ind2) = sth cloudth_sigmaenv(ind1,ind2) = sigma_env cloudth_sigmath(ind1,ind2) = sigma_th END DO ! From the loop on ngrid return END SUBROUTINE cloudth_v6 END MODULE lmdz_lscp_condensation