MODULE lmdz_aviation !---------------------------------------------------------------- ! Module for aviation and contrails IMPLICIT NONE CONTAINS !********************************************************************************** SUBROUTINE contrails_formation_evolution( & dtime, pplay, temp, qsat, qsatl, gamma_cond, rcont_seri, flight_dist, & cldfra, qvc, V_cell, M_cell, pdf_loc, pdf_scale, pdf_alpha, & Tcritcont, qcritcont, potcontfraP, potcontfraNP, contfra, & dcf_avi, dqvc_avi, dqi_avi & ) USE lmdz_lscp_ini, ONLY: RCPD, EPS_W, RTT USE lmdz_lscp_ini, ONLY: EI_H2O_aviation, qheat_fuel_aviation, prop_efficiency_aviation USE lmdz_lscp_ini, ONLY: linear_contrails_lifetime USE lmdz_lscp_ini, ONLY: eps USE lmdz_lscp_tools, ONLY: GAMMAINC, calc_qsat_ecmwf IMPLICIT NONE ! ! Input ! REAL, INTENT(IN) :: dtime ! time step [s] REAL, INTENT(IN) :: pplay ! layer pressure [Pa] REAL, INTENT(IN) :: temp ! temperature [K] REAL, INTENT(IN) :: qsat ! saturation specific humidity [kg/kg] REAL, INTENT(IN) :: qsatl ! saturation specific humidity w.r.t. liquid [kg/kg] REAL, INTENT(IN) :: gamma_cond ! condensation threshold w.r.t. qsat [-] REAL, INTENT(IN) :: rcont_seri ! ratio of contrails fraction to total cloud fraction [-] REAL, INTENT(IN) :: flight_dist ! REAL, INTENT(IN) :: cldfra ! cloud fraction [-] REAL, INTENT(IN) :: qvc ! gridbox-mean vapor in the cloud [kg/kg] REAL, INTENT(IN) :: V_cell ! cell volume [m3] REAL, INTENT(IN) :: M_cell ! cell mass [kg] REAL, INTENT(IN) :: pdf_loc ! location parameter of the clear sky PDF [%] REAL, INTENT(IN) :: pdf_scale ! scale parameter of the clear sky PDF [%] REAL, INTENT(IN) :: pdf_alpha ! shape parameter of the clear sky PDF [-] ! ! Output ! REAL, INTENT(OUT) :: Tcritcont ! REAL, INTENT(OUT) :: qcritcont ! REAL, INTENT(OUT) :: potcontfraP ! REAL, INTENT(OUT) :: potcontfraNP ! REAL, INTENT(OUT) :: contfra ! REAL, INTENT(OUT) :: dcf_avi ! REAL, INTENT(OUT) :: dqvc_avi ! REAL, INTENT(OUT) :: dqi_avi ! ! ! Local ! ! for Schmidt-Applemant criteria REAL, DIMENSION(1) :: ztemp, zpplay, qzero, zqsatl, zdqsatl REAL :: Gcont, qsatl_crit, psatl_crit, pcrit REAL :: pdf_x, pdf_y, pdf_e2, pdf_e3 REAL :: pdf_fra_above_qcritcont, pdf_fra_above_qsat, pdf_fra_above_qnuc REAL :: pdf_q_above_qcritcont, pdf_q_above_qsat, pdf_q_above_qnuc REAL :: qpotcontP ! ! REAL :: contrail_cross_section, contfra_new qzero(:) = 0. !--more local variables for diagnostics !--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. !--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) !--------------------------------- !-- SCHIMDT-APPLEMAN CRITERIA -- !--------------------------------- !--Revised by Schumann (1996) and Rap et al. (2010) !--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1) !--in Pa / K Gcont = EI_H2O_aviation * RCPD * pplay & / ( EPS_W * qheat_fuel_aviation * ( 1. - prop_efficiency_aviation ) ) !--critical T_LM below which no liquid contrail can form in exhaust !--in Kelvins Tcritcont = 226.69 + 9.43 * LOG( MAX(Gcont, 0.1) - 0.053 ) & + 0.72 * LOG( MAX(Gcont, 0.1) - 0.053 )**2 ztemp(1) = Tcritcont zpplay(1) = pplay CALL calc_qsat_ecmwf(1, ztemp, qzero, zpplay, RTT, 1, .FALSE., zqsatl, zdqsatl) qsatl_crit = zqsatl(1) psatl_crit = qsatl_crit / ( EPS_W + ( 1. - EPS_W ) * qsatl_crit ) * pplay pcrit = Gcont * ( temp - Tcritcont ) + psatl_crit qcritcont = EPS_W * pcrit / ( pplay - ( 1. - EPS_W ) * pcrit ) IF ( temp .LT. Tcritcont ) THEN !--contrail formation is possible pdf_x = qcritcont / qsatl * 100. pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha pdf_e2 = GAMMA(1. + 1. / pdf_alpha) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y ) pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2 pdf_fra_above_qcritcont = EXP( - pdf_y ) * ( 1. - cldfra ) pdf_q_above_qcritcont = ( pdf_e3 * ( 1. - cldfra ) + pdf_loc * pdf_fra_above_qcritcont ) * qsatl / 100. pdf_x = gamma_cond * qsat / qsatl * 100. pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha pdf_e2 = GAMMA(1. + 1. / pdf_alpha) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y ) pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2 pdf_fra_above_qnuc = EXP( - pdf_y ) * ( 1. - cldfra ) pdf_q_above_qnuc = ( pdf_e3 * ( 1. - cldfra ) + pdf_loc * pdf_fra_above_qnuc ) * qsatl / 100. pdf_x = qsat / qsatl * 100. pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha pdf_e2 = GAMMA(1. + 1. / pdf_alpha) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y ) pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2 pdf_fra_above_qsat = EXP( - pdf_y ) * ( 1. - cldfra ) pdf_q_above_qsat = ( pdf_e3 * ( 1. - cldfra ) + pdf_loc * pdf_fra_above_qsat ) * qsatl / 100. potcontfraNP = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat) potcontfraP = MIN(pdf_fra_above_qsat - pdf_fra_above_qnuc, & pdf_fra_above_qcritcont - pdf_fra_above_qnuc) qpotcontP = MIN(pdf_q_above_qsat - pdf_q_above_qnuc, & pdf_q_above_qcritcont - pdf_q_above_qnuc) ELSE potcontfraNP = 0. potcontfraP = 0. ENDIF ! temp .LT. Tcritcont ! 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 !--Convert existing contrail fraction into "natural" cirrus cloud fraction contfra = rcont_seri * cldfra contfra = contfra * EXP( - dtime / linear_contrails_lifetime ) !--Add a source of contrails from aviation dcf_avi = 0. dqi_avi = 0. dqvc_avi = 0. IF ( potcontfraP .GT. eps ) THEN contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA() contfra_new = MIN(1., flight_dist * dtime * contrail_cross_section / V_cell) dcf_avi = potcontfraP * contfra_new IF ( cldfra .GT. eps ) THEN dqvc_avi = dcf_avi * qvc / cldfra ELSE dqvc_avi = dcf_avi * qsat ENDIF dqi_avi = dcf_avi * qpotcontP / potcontfraP - dqvc_avi !--Add tendencies contfra = contfra + contfra_new ENDIF END SUBROUTINE contrails_formation_evolution !********************************************************************************** SUBROUTINE contrails_mixing( & dtime, pplay, shear, pbl_eps, cell_area, dz, temp, qtot, qsat, & subfra, qsub, issrfra, qissr, cldfra, qcld, qvc, rcont_seri, & dcf_mix, dqvc_mix, dqi_mix, dqt_mix, dcf_mix_issr, dqt_mix_issr & ) !---------------------------------------------------------------------- ! This subroutine calculates the mixing of contrails in their environment. ! Authors: Audran Borella, Etienne Vignon, Olivier Boucher ! December 2024 !---------------------------------------------------------------------- USE lmdz_lscp_ini, ONLY: RPI, eps, ok_unadjusted_clouds USE lmdz_lscp_ini, ONLY: aspect_ratio_contrails, coef_mixing_contrails, & coef_shear_contrails, chi_mixing_contrails IMPLICIT NONE ! ! Input ! REAL, INTENT(IN) :: dtime ! time step [s] ! REAL, INTENT(IN) :: pplay ! layer pressure [Pa] REAL, INTENT(IN) :: shear ! vertical shear [s-1] REAL, INTENT(IN) :: pbl_eps ! TKE dissipation[m2/s3] REAL, INTENT(IN) :: cell_area ! cell area [m2] REAL, INTENT(IN) :: dz ! cell width [m] REAL, INTENT(IN) :: temp ! temperature [K] REAL, INTENT(IN) :: qtot ! total specific humidity (without precip) [kg/kg] REAL, INTENT(IN) :: qsat ! saturation specific humidity [kg/kg] REAL, INTENT(IN) :: subfra ! subsaturated clear sky fraction [-] REAL, INTENT(IN) :: qsub ! gridbox-mean subsaturated clear sky specific water [kg/kg] REAL, INTENT(IN) :: issrfra ! ISSR fraction [-] REAL, INTENT(IN) :: qissr ! gridbox-mean ISSR specific water [kg/kg] REAL, INTENT(IN) :: cldfra ! cloud fraction [-] REAL, INTENT(IN) :: qcld ! gridbox-mean cloudy specific total water [kg/kg] REAL, INTENT(IN) :: qvc ! gridbox-mean cloudy specific water vapor [kg/kg] REAL, INTENT(IN) :: rcont_seri ! ratio of contrails fraction to total cloud fraction [-] ! ! Input/Output ! REAL, INTENT(INOUT) :: dcf_mix ! cloud fraction tendency because of cloud mixing [s-1] REAL, INTENT(INOUT) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s] REAL, INTENT(INOUT) :: dqi_mix ! specific ice content tendency because of cloud mixing [kg/kg/s] REAL, INTENT(INOUT) :: dqt_mix ! total water tendency because of cloud mixing [kg/kg/s] REAL, INTENT(INOUT) :: dcf_mix_issr ! cloud fraction tendency because of cloud mixing in ISSR [s-1] REAL, INTENT(INOUT) :: dqt_mix_issr ! total water tendency because of cloud mixing in ISSR [kg/kg/s] ! ! Local ! REAL :: dqt_mix_sub_cont, dqt_mix_issr_cont REAL :: dcf_mix_sub_cont, dcf_mix_issr_cont REAL :: dqvc_mix_sub_cont, dqvc_mix_issr_cont REAL :: dcf_mix_cont, dqvc_mix_cont, dqi_mix_cont, dqt_mix_cont 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 REAL :: qvapincld, qvapinmix, qvapincld_new, qiceincld !--Initialisation dcf_mix_sub_cont = 0. dqt_mix_sub_cont = 0. dqvc_mix_sub_cont = 0. dcf_mix_issr_cont = 0. dqt_mix_issr_cont = 0. dqvc_mix_issr_cont = 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 = aspect_ratio_contrails !--P / a is a function of b / a only, that we can calculate !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) ) Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) ) !--The clouds perimeter is imposed using the formula from Morcrette 2012, !--based on observations. !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra ) !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have: !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra ) !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) ) a_mix = Povera / coef_mixing_contrails / RPI / ( 1. - cldfra ) / bovera !--and finally, !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a ) N_cld_mix = coef_mixing_contrails * cldfra * ( 1. - cldfra ) * cell_area & / 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 ) !--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 & * ( 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 & * ( 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_contrails * shear * 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 !--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 / ( subfra + chi_mixing_contrails * issrfra ) subfra_mix = MIN( sigma_mix * envfra_mix, subfra ) issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra ) cldfra_mix = MIN( cldfra_mix, cldfra ) !--First, we mix the subsaturated sky (subfra_mix) and the cloud close !--to this fraction (sigma_mix * cldfra_mix). IF ( subfra .GT. eps ) THEN !--We compute the total humidity in the mixed air, which !--can be either sub- or supersaturated. qvapinmix = ( qsub * subfra_mix / subfra & + qcld * cldfra_mix * sigma_mix / cldfra ) & / ( subfra_mix + cldfra_mix * sigma_mix ) IF ( ok_unadjusted_clouds ) THEN qiceincld = ( qcld - qvc ) * cldfra_mix * sigma_mix / cldfra & / ( subfra_mix + cldfra_mix * sigma_mix ) qvapincld_new = QVAPINCLD_DEPSUB_CONTRAILS( & qvapinmix, qiceincld, temp, qsat, pplay, dtime) 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_cont = - sigma_mix * cldfra_mix dqt_mix_sub_cont = dcf_mix_sub_cont * qcld / cldfra dqvc_mix_sub_cont = dcf_mix_sub_cont * qvc / cldfra ELSE dcf_mix_sub_cont = subfra_mix dqt_mix_sub_cont = dcf_mix_sub_cont * qsub / subfra dqvc_mix_sub_cont = dcf_mix_sub_cont * qvapincld_new ENDIF ELSE IF ( qvapinmix .GT. qsat ) THEN !--If the mixed air is supersaturated, we condense the subsaturated !--region which was mixed. dcf_mix_sub_cont = subfra_mix dqt_mix_sub_cont = dcf_mix_sub_cont * qsub / subfra dqvc_mix_sub_cont = dcf_mix_sub_cont * qsat ELSE !--Else, we sublimate the cloud which was mixed. dcf_mix_sub_cont = - sigma_mix * cldfra_mix dqt_mix_sub_cont = dcf_mix_sub_cont * qcld / cldfra dqvc_mix_sub_cont = dcf_mix_sub_cont * qsat 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 .GT. eps ) THEN IF ( ok_unadjusted_clouds ) THEN qvapinmix = ( qissr * issrfra_mix / issrfra & + qcld * cldfra_mix * (1. - sigma_mix) / cldfra ) & / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) ) qiceincld = ( qcld - qvc ) * cldfra_mix * (1. - sigma_mix) / cldfra & / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) ) qvapincld_new = QVAPINCLD_DEPSUB_CONTRAILS( & qvapinmix, qiceincld, temp, qsat, pplay, dtime) dcf_mix_issr_cont = issrfra_mix dqt_mix_issr_cont = dcf_mix_issr_cont * qissr / issrfra dqvc_mix_issr_cont = dcf_mix_issr_cont * qvapincld_new ELSE !--In this case, the additionnal vapor condenses dcf_mix_issr_cont = issrfra_mix dqt_mix_issr_cont = dcf_mix_issr_cont * qissr / issrfra dqvc_mix_issr_cont = dcf_mix_issr_cont * qsat ENDIF ! ok_unadjusted_clouds ENDIF ! issrfra .GT. eps !--Sum up the tendencies from subsaturated sky and supersaturated sky dcf_mix_cont = dcf_mix_sub_cont + dcf_mix_issr_cont dqt_mix_cont = dqt_mix_sub_cont + dqt_mix_issr_cont dqvc_mix_cont = dqvc_mix_sub_cont + dqvc_mix_issr_cont dqi_mix_cont = dqt_mix_cont - dqvc_mix_cont !--Outputs !--The mixing tendencies are a linear combination of the calculation done for "classical" cirrus !--and contrails dcf_mix = ( 1. - rcont_seri ) * dcf_mix + rcont_seri * dcf_mix_cont dqvc_mix = ( 1. - rcont_seri ) * dqvc_mix + rcont_seri * dqvc_mix_cont dqi_mix = ( 1. - rcont_seri ) * dqi_mix + rcont_seri * dqi_mix_cont dqt_mix = ( 1. - rcont_seri ) * dqt_mix + rcont_seri * dqt_mix_cont dcf_mix_issr = ( 1. - rcont_seri ) * dcf_mix_issr + rcont_seri * dcf_mix_issr_cont dqt_mix_issr = ( 1. - rcont_seri ) * dqt_mix_issr + rcont_seri * dqt_mix_issr_cont END SUBROUTINE contrails_mixing !********************************************************************************** FUNCTION qvapincld_depsub_contrails( & qvapincld, qiceincld, temp, qsat, pplay, dtime) USE lmdz_lscp_ini, ONLY: RV, RLSTT, RTT, EPS_W USE lmdz_lscp_ini, ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice USE lmdz_lscp_ini, ONLY: rm_ice_crystals_contrails IMPLICIT NONE ! inputs REAL :: qvapincld ! REAL :: qiceincld ! REAL :: temp ! REAL :: qsat ! REAL :: pplay ! REAL :: dtime ! time step [s] ! output REAL :: qvapincld_depsub_contrails ! local REAL :: qice_ratio REAL :: pres_sat, rho, kappa REAL :: air_thermal_conduct, water_vapor_diff REAL :: 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 !--Here, the initial vapor in the cloud is qvapincld, and we compute !--the new vapor qvapincld_depsub_contrails rm_ice = rm_ice_crystals_contrails IF ( qvapincld .GE. qsat ) THEN !--If the cloud is initially supersaturated !--Exact explicit integration (qvc exact, qice explicit) qvapincld_depsub_contrails = 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_depsub_contrails = 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_depsub_contrails = qvapincld + qiceincld * ( 1. - qice_ratio**1.5 ) !--The new vapor in the cloud cannot be greater than qsat qvapincld_depsub_contrails = MIN(qvapincld_depsub_contrails, qsat) ENDIF ! qice_ratio .LT. 0. ENDIF ! qvapincld .GT. qsat END FUNCTION qvapincld_depsub_contrails !********************************************************************************** FUNCTION contrail_cross_section_onera() IMPLICIT NONE ! input ! output REAL :: contrail_cross_section_onera ! [m2] ! local contrail_cross_section_onera = 200. * 200. END FUNCTION contrail_cross_section_onera END MODULE lmdz_aviation !******************************************************************* ! !SUBROUTINE airplane(debut,pphis,pplay,paprs,t_seri) ! ! USE dimphy ! USE mod_grid_phy_lmdz, ONLY: klon_glo ! USE geometry_mod, ONLY: cell_area ! USE phys_cal_mod, ONLY : mth_cur ! USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root ! USE mod_phys_lmdz_omp_data, ONLY: is_omp_root ! USE mod_phys_lmdz_para, ONLY: scatter, bcast ! USE print_control_mod, ONLY: lunout ! ! IMPLICIT NONE ! ! INCLUDE "YOMCST.h" ! INCLUDE 'netcdf.inc' ! ! !-------------------------------------------------------- ! !--input variables ! !-------------------------------------------------------- ! LOGICAL, INTENT(IN) :: debut ! REAL, INTENT(IN) :: pphis(klon), pplay(klon,klev), paprs(klon,klev+1), t_seri(klon,klev) ! ! !-------------------------------------------------------- ! ! ... Local variables ! !-------------------------------------------------------- ! ! CHARACTER (LEN=20) :: modname='airplane_mod' ! INTEGER :: i, k, kori, iret, varid, error, ncida, klona ! INTEGER,SAVE :: nleva, ntimea !!$OMP THREADPRIVATE(nleva,ntimea) ! REAL, ALLOCATABLE :: pkm_airpl_glo(:,:,:) !--km/s ! REAL, ALLOCATABLE :: ph2o_airpl_glo(:,:,:) !--molec H2O/cm3/s ! REAL, ALLOCATABLE, SAVE :: zmida(:), zinta(:) ! REAL, ALLOCATABLE, SAVE :: pkm_airpl(:,:,:) ! REAL, ALLOCATABLE, SAVE :: ph2o_airpl(:,:,:) !!$OMP THREADPRIVATE(pkm_airpl,ph2o_airpl,zmida,zinta) ! REAL :: zalt(klon,klev+1) ! REAL :: zrho, zdz(klon,klev), zfrac ! ! ! ! IF (debut) THEN ! !-------------------------------------------------------------------------------- ! ! ... Open the file and read airplane emissions ! !-------------------------------------------------------------------------------- ! ! ! IF (is_mpi_root .AND. is_omp_root) THEN ! ! ! iret = nf_open('aircraft_phy.nc', 0, ncida) ! IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to open aircraft_phy.nc file',1) ! ! ... Get lengths ! iret = nf_inq_dimid(ncida, 'time', varid) ! IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimid in aircraft_phy.nc file',1) ! iret = nf_inq_dimlen(ncida, varid, ntimea) ! IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimlen aircraft_phy.nc file',1) ! iret = nf_inq_dimid(ncida, 'vector', varid) ! IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimid aircraft_phy.nc file',1) ! iret = nf_inq_dimlen(ncida, varid, klona) ! IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimlen aircraft_phy.nc file',1) ! iret = nf_inq_dimid(ncida, 'lev', varid) ! IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1) ! iret = nf_inq_dimlen(ncida, varid, nleva) ! IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimlen aircraft_phy.nc file',1) ! ! ! IF ( klona /= klon_glo ) THEN ! WRITE(lunout,*) 'klona & klon_glo =', klona, klon_glo ! CALL abort_physic(modname,'problem klon in aircraft_phy.nc file',1) ! ENDIF ! ! ! IF ( ntimea /= 12 ) THEN ! WRITE(lunout,*) 'ntimea=', ntimea ! CALL abort_physic(modname,'problem ntime<>12 in aircraft_phy.nc file',1) ! ENDIF ! ! ! ALLOCATE(zmida(nleva), STAT=error) ! IF (error /= 0) CALL abort_physic(modname,'problem to allocate zmida',1) ! ALLOCATE(pkm_airpl_glo(klona,nleva,ntimea), STAT=error) ! IF (error /= 0) CALL abort_physic(modname,'problem to allocate pkm_airpl_glo',1) ! ALLOCATE(ph2o_airpl_glo(klona,nleva,ntimea), STAT=error) ! IF (error /= 0) CALL abort_physic(modname,'problem to allocate ph2o_airpl_glo',1) ! ! ! iret = nf_inq_varid(ncida, 'lev', varid) ! IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1) ! iret = nf_get_var_double(ncida, varid, zmida) ! IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read zmida file',1) ! ! ! iret = nf_inq_varid(ncida, 'emi_co2_aircraft', varid) !--CO2 as a proxy for m flown - ! IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_distance dimid aircraft_phy.nc file',1) ! iret = nf_get_var_double(ncida, varid, pkm_airpl_glo) ! IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read pkm_airpl file',1) ! ! ! iret = nf_inq_varid(ncida, 'emi_h2o_aircraft', varid) ! IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_h2o_aircraft dimid aircraft_phy.nc file',1) ! iret = nf_get_var_double(ncida, varid, ph2o_airpl_glo) ! IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read ph2o_airpl file',1) ! ! ! ENDIF !--is_mpi_root and is_omp_root ! ! ! CALL bcast(nleva) ! CALL bcast(ntimea) ! ! ! IF (.NOT.ALLOCATED(zmida)) ALLOCATE(zmida(nleva), STAT=error) ! IF (.NOT.ALLOCATED(zinta)) ALLOCATE(zinta(nleva+1), STAT=error) ! ! ! ALLOCATE(pkm_airpl(klon,nleva,ntimea)) ! ALLOCATE(ph2o_airpl(klon,nleva,ntimea)) ! ! ! ALLOCATE(flight_m(klon,klev)) ! ALLOCATE(flight_h2o(klon,klev)) ! ! ! CALL bcast(zmida) ! zinta(1)=0.0 !--surface ! DO k=2, nleva ! zinta(k) = (zmida(k-1)+zmida(k))/2.0*1000.0 !--conversion from km to m ! ENDDO ! zinta(nleva+1)=zinta(nleva)+(zmida(nleva)-zmida(nleva-1))*1000.0 !--extrapolation for last interface ! !print *,'zinta=', zinta ! ! ! CALL scatter(pkm_airpl_glo,pkm_airpl) ! CALL scatter(ph2o_airpl_glo,ph2o_airpl) ! ! !!$OMP MASTER ! IF (is_mpi_root .AND. is_omp_root) THEN ! DEALLOCATE(pkm_airpl_glo) ! DEALLOCATE(ph2o_airpl_glo) ! ENDIF !--is_mpi_root !!$OMP END MASTER ! ! ENDIF !--debut !! !!--compute altitude of model level interfaces !! ! DO i = 1, klon ! zalt(i,1)=pphis(i)/RG !--in m ! ENDDO !! ! DO k=1, klev ! DO i = 1, klon ! zrho=pplay(i,k)/t_seri(i,k)/RD ! zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho/RG ! zalt(i,k+1)=zalt(i,k)+zdz(i,k) !--in m ! ENDDO ! ENDDO !! !!--vertical reprojection !! ! flight_m(:,:)=0.0 ! flight_h2o(:,:)=0.0 !! ! DO k=1, klev ! DO kori=1, nleva ! DO i=1, klon ! !--fraction of layer kori included in layer k ! zfrac=max(0.0,min(zalt(i,k+1),zinta(kori+1))-max(zalt(i,k),zinta(kori)))/(zinta(kori+1)-zinta(kori)) ! !--reproject ! flight_m(i,k)=flight_m(i,k) + pkm_airpl(i,kori,mth_cur)*zfrac ! !--reproject ! flight_h2o(i,k)=flight_h2o(i,k) + ph2o_airpl(i,kori,mth_cur)*zfrac ! ENDDO ! ENDDO ! ENDDO !! ! DO k=1, klev ! DO i=1, klon ! !--molec.cm-3.s-1 / (molec/mol) * kg CO2/mol * m2 * m * cm3/m3 / (kg CO2/m) => m s-1 per cell ! flight_m(i,k)=flight_m(i,k)/RNAVO*44.e-3*cell_area(i)*zdz(i,k)*1.e6/16.37e-3 ! flight_m(i,k)=flight_m(i,k)*100.0 !--x100 to augment signal to noise ! !--molec.cm-3.s-1 / (molec/mol) * kg H2O/mol * m2 * m * cm3/m3 => kg H2O s-1 per cell ! flight_h2o(i,k)=flight_h2o(i,k)/RNAVO*18.e-3*cell_area(i)*zdz(i,k)*1.e6 ! ENDDO ! ENDDO !! !END SUBROUTINE airplane ! !!******************************************************************** !! simple routine to initialise flight_m and test a flight corridor !!--Olivier Boucher - 2021 !! !SUBROUTINE flight_init() ! USE dimphy ! USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg ! IMPLICIT NONE ! INTEGER :: i ! ! ALLOCATE(flight_m(klon,klev)) ! ALLOCATE(flight_h2o(klon,klev)) ! ! ! flight_m(:,:) = 0.0 !--initialisation ! flight_h2o(:,:) = 0.0 !--initialisation ! ! ! DO i=1, klon ! IF (latitude_deg(i).GE.42.0.AND.latitude_deg(i).LE.48.0) THEN ! flight_m(i,38) = 50000.0 !--5000 m of flight/second in grid cell x 10 scaling ! ENDIF ! ENDDO ! ! RETURN !END SUBROUTINE flight_init ! ! !--add a source of cirrus from aviation contrails ! IF (ok_plane_contrail) THEN ! drneb_avi(i,k) = rnebss*flight_m(i,k)*contrail_cross_section/V_cell !--tendency rneb due to aviation [s-1] ! drneb_avi(i,k) = MIN(drneb_avi(i,k), rnebss/dtime) !--majoration ! dqss_avi = qss*drneb_avi(i,k)/MAX(eps,rnebss) !--tendency q aviation [kg kg-1 s-1] ! rneb = rneb + drneb_avi(i,k)*dtime !--add tendency to rneb ! qcld = qcld + dqss_avi*dtime !--add tendency to qcld ! rnebss = rnebss - drneb_avi(i,k)*dtime !--add tendency to rnebss ! qss = qss - dqss_avi*dtime !--add tendency to qss ! ELSE ! drneb_avi(i,k)=0.0 ! ENDIF ! ! RETURN !END SUBROUTINE ice_sursat