source: LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90 @ 5452

Last change on this file since 5452 was 5452, checked in by aborella, 5 days ago

First implementation of the contrails parameterisation
Lacks the emission of H2O + the impact on radiative transfer
File size: 33.3 KB
Line 
1MODULE lmdz_aviation
2!----------------------------------------------------------------
3! Module for aviation and contrails
4
5IMPLICIT NONE
6
7CONTAINS
8
9
10!**********************************************************************************
11SUBROUTINE contrails_formation_evolution( &
12      dtime, pplay, temp, qsat, qsatl, gamma_cond, rcont_seri, flight_dist, &
13      cldfra, qvc, V_cell, M_cell, pdf_loc, pdf_scale, pdf_alpha, &
14      Tcritcont, qcritcont, potcontfraP, potcontfraNP, contfra, &
15      dcf_avi, dqvc_avi, dqi_avi &
16      )
17
18USE lmdz_lscp_ini, ONLY: RCPD, EPS_W, RTT
19USE lmdz_lscp_ini, ONLY: EI_H2O_aviation, qheat_fuel_aviation, prop_efficiency_aviation
20USE lmdz_lscp_ini, ONLY: linear_contrails_lifetime
21USE lmdz_lscp_ini, ONLY: eps
22
23USE lmdz_lscp_tools, ONLY: GAMMAINC, calc_qsat_ecmwf
24
25IMPLICIT NONE
26
27!
28! Input
29!
30REAL, INTENT(IN)  :: dtime        ! time step [s]
31REAL, INTENT(IN)  :: pplay        ! layer pressure [Pa]
32REAL, INTENT(IN)  :: temp         ! temperature [K]
33REAL, INTENT(IN)  :: qsat         ! saturation specific humidity [kg/kg]
34REAL, INTENT(IN)  :: qsatl        ! saturation specific humidity w.r.t. liquid [kg/kg]
35REAL, INTENT(IN)  :: gamma_cond   ! condensation threshold w.r.t. qsat [-]
36REAL, INTENT(IN)  :: rcont_seri   ! ratio of contrails fraction to total cloud fraction [-]
37REAL, INTENT(IN)  :: flight_dist  !
38REAL, INTENT(IN)  :: cldfra       ! cloud fraction [-]
39REAL, INTENT(IN)  :: qvc          ! gridbox-mean vapor in the cloud [kg/kg]
40REAL, INTENT(IN)  :: V_cell       ! cell volume [m3]
41REAL, INTENT(IN)  :: M_cell       ! cell mass [kg]
42REAL, INTENT(IN)  :: pdf_loc      ! location parameter of the clear sky PDF [%]
43REAL, INTENT(IN)  :: pdf_scale    ! scale parameter of the clear sky PDF [%]
44REAL, INTENT(IN)  :: pdf_alpha    ! shape parameter of the clear sky PDF [-]
45!
46! Output
47!
48REAL, INTENT(OUT) :: Tcritcont    !
49REAL, INTENT(OUT) :: qcritcont    !
50REAL, INTENT(OUT) :: potcontfraP  !
51REAL, INTENT(OUT) :: potcontfraNP !
52REAL, INTENT(OUT) :: contfra      !
53REAL, INTENT(OUT) :: dcf_avi      !
54REAL, INTENT(OUT) :: dqvc_avi     !
55REAL, INTENT(OUT) :: dqi_avi      !
56!
57! Local
58!
59! for Schmidt-Applemant criteria
60REAL, DIMENSION(1) :: ztemp, zpplay, qzero, zqsatl, zdqsatl
61REAL :: Gcont, qsatl_crit, psatl_crit, pcrit
62REAL :: pdf_x, pdf_y, pdf_e2, pdf_e3
63REAL :: pdf_fra_above_qcritcont, pdf_fra_above_qsat, pdf_fra_above_qnuc
64REAL :: pdf_q_above_qcritcont, pdf_q_above_qsat, pdf_q_above_qnuc
65REAL :: qpotcontP
66!
67!
68REAL :: contrail_cross_section, contfra_new
69
70qzero(:) = 0.
71
72!--more local variables for diagnostics
73!--values from Schumann, Meteorol Zeitschrift, 1996
74!--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen
75!--Qheat = 43.  /  50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen
76!REAL, PARAMETER :: EiH2O=1.25  !--emission index of water vapour for kerosene (kg kg-1)
77!REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1)
78!REAL, PARAMETER :: eta=0.3     !--average propulsion efficiency of the aircraft
79!--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
80!--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010.
81!--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2)
82
83!---------------------------------
84!--  SCHIMDT-APPLEMAN CRITERIA  --
85!---------------------------------
86!--Revised by Schumann (1996) and Rap et al. (2010)
87
88!--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1)
89!--in Pa / K
90Gcont = EI_H2O_aviation * RCPD * pplay &
91       / ( EPS_W * qheat_fuel_aviation * ( 1. - prop_efficiency_aviation ) )
92!--critical T_LM below which no liquid contrail can form in exhaust
93!--in Kelvins
94Tcritcont = 226.69 + 9.43 * LOG( MAX(Gcont, 0.1) - 0.053 ) &
95       + 0.72 * LOG( MAX(Gcont, 0.1) - 0.053 )**2
96
97ztemp(1) = Tcritcont
98zpplay(1) = pplay
99CALL calc_qsat_ecmwf(1, ztemp, qzero, zpplay, RTT, 1, .FALSE., zqsatl, zdqsatl)
100qsatl_crit = zqsatl(1)
101
102psatl_crit = qsatl_crit / ( EPS_W + ( 1. - EPS_W ) * qsatl_crit ) * pplay
103pcrit = Gcont * ( temp - Tcritcont ) + psatl_crit
104qcritcont = EPS_W * pcrit / ( pplay - ( 1. - EPS_W ) * pcrit )
105
106
107IF ( temp .LT. Tcritcont ) THEN !--contrail formation is possible
108
109  pdf_x = qcritcont / qsatl * 100.
110  pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
111  pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
112  pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
113  pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
114  pdf_fra_above_qcritcont = EXP( - pdf_y ) * ( 1. - cldfra )
115  pdf_q_above_qcritcont = ( pdf_e3 * ( 1. - cldfra ) + pdf_loc * pdf_fra_above_qcritcont ) * qsatl / 100.
116
117  pdf_x = gamma_cond * qsat / qsatl * 100.
118  pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
119  pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
120  pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
121  pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
122  pdf_fra_above_qnuc = EXP( - pdf_y ) * ( 1. - cldfra )
123  pdf_q_above_qnuc = ( pdf_e3 * ( 1. - cldfra ) + pdf_loc * pdf_fra_above_qnuc ) * qsatl / 100.
124
125  pdf_x = qsat / qsatl * 100.
126  pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
127  pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
128  pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
129  pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
130  pdf_fra_above_qsat = EXP( - pdf_y ) * ( 1. - cldfra )
131  pdf_q_above_qsat = ( pdf_e3 * ( 1. - cldfra ) + pdf_loc * pdf_fra_above_qsat ) * qsatl / 100.
132
133  potcontfraNP = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat)
134  potcontfraP = MIN(pdf_fra_above_qsat - pdf_fra_above_qnuc, &
135                    pdf_fra_above_qcritcont - pdf_fra_above_qnuc)
136  qpotcontP = MIN(pdf_q_above_qsat - pdf_q_above_qnuc, &
137                  pdf_q_above_qcritcont - pdf_q_above_qnuc)
138
139ELSE
140
141  potcontfraNP = 0.
142  potcontfraP = 0.
143
144ENDIF ! temp .LT. Tcritcont
145
146
147! Ajout des émissions de H2O dues à l'aviation
148! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
149! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O ) 
150!      = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) )
151! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q)
152! flight_h2O is in kg H2O / s / cell
153!
154!IF (ok_plane_h2o) THEN
155!  q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) )
156!ENDIF
157
158!--Convert existing contrail fraction into "natural" cirrus cloud fraction
159contfra = rcont_seri * cldfra
160contfra = contfra * EXP( - dtime / linear_contrails_lifetime )
161
162!--Add a source of contrails from aviation
163dcf_avi = 0.
164dqi_avi = 0.
165dqvc_avi = 0.
166IF ( potcontfraP .GT. eps ) THEN
167  contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA()
168  contfra_new = MIN(1., flight_dist * dtime * contrail_cross_section / V_cell)
169  dcf_avi = potcontfraP * contfra_new
170  IF ( cldfra .GT. eps ) THEN
171    dqvc_avi = dcf_avi * qvc / cldfra
172  ELSE
173    dqvc_avi = dcf_avi * qsat
174  ENDIF
175  dqi_avi = dcf_avi * qpotcontP / potcontfraP - dqvc_avi
176
177  !--Add tendencies
178  contfra = contfra + contfra_new
179ENDIF
180
181END SUBROUTINE contrails_formation_evolution
182
183
184!**********************************************************************************
185SUBROUTINE contrails_mixing( &
186      dtime, pplay, shear, pbl_eps, cell_area, dz, temp, qtot, qsat, &
187      subfra, qsub, issrfra, qissr, cldfra, qcld, qvc, rcont_seri, &
188      dcf_mix, dqvc_mix, dqi_mix, dqt_mix, dcf_mix_issr, dqt_mix_issr &
189      )
190
191!----------------------------------------------------------------------
192! This subroutine calculates the mixing of contrails in their environment.
193! Authors: Audran Borella, Etienne Vignon, Olivier Boucher
194! December 2024
195!----------------------------------------------------------------------
196
197USE lmdz_lscp_ini,        ONLY: RPI, eps, ok_unadjusted_clouds
198USE lmdz_lscp_ini,        ONLY: aspect_ratio_contrails, coef_mixing_contrails, &
199                                coef_shear_contrails, chi_mixing_contrails
200
201IMPLICIT NONE
202
203!
204! Input
205!
206REAL, INTENT(IN)    :: dtime        ! time step [s]
207!
208REAL, INTENT(IN)    :: pplay        ! layer pressure [Pa]
209REAL, INTENT(IN)    :: shear        ! vertical shear [s-1]
210REAL, INTENT(IN)    :: pbl_eps      ! TKE dissipation[m2/s3]
211REAL, INTENT(IN)    :: cell_area    ! cell area [m2]
212REAL, INTENT(IN)    :: dz           ! cell width [m]
213REAL, INTENT(IN)    :: temp         ! temperature [K]
214REAL, INTENT(IN)    :: qtot         ! total specific humidity (without precip) [kg/kg]
215REAL, INTENT(IN)    :: qsat         ! saturation specific humidity [kg/kg]
216REAL, INTENT(IN)    :: subfra       ! subsaturated clear sky fraction [-]
217REAL, INTENT(IN)    :: qsub         ! gridbox-mean subsaturated clear sky specific water [kg/kg]
218REAL, INTENT(IN)    :: issrfra      ! ISSR fraction [-]
219REAL, INTENT(IN)    :: qissr        ! gridbox-mean ISSR specific water [kg/kg]
220REAL, INTENT(IN)    :: cldfra       ! cloud fraction [-]
221REAL, INTENT(IN)    :: qcld         ! gridbox-mean cloudy specific total water [kg/kg]
222REAL, INTENT(IN)    :: qvc          ! gridbox-mean cloudy specific water vapor [kg/kg]
223REAL, INTENT(IN)    :: rcont_seri   ! ratio of contrails fraction to total cloud fraction [-]
224!
225!  Input/Output
226!
227REAL, INTENT(INOUT) :: dcf_mix      ! cloud fraction tendency because of cloud mixing [s-1]
228REAL, INTENT(INOUT) :: dqvc_mix     ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
229REAL, INTENT(INOUT) :: dqi_mix      ! specific ice content tendency because of cloud mixing [kg/kg/s]
230REAL, INTENT(INOUT) :: dqt_mix      ! total water tendency because of cloud mixing [kg/kg/s]
231REAL, INTENT(INOUT) :: dcf_mix_issr ! cloud fraction tendency because of cloud mixing in ISSR [s-1]
232REAL, INTENT(INOUT) :: dqt_mix_issr ! total water tendency because of cloud mixing in ISSR [kg/kg/s]
233!
234! Local
235!
236REAL :: dqt_mix_sub_cont, dqt_mix_issr_cont
237REAL :: dcf_mix_sub_cont, dcf_mix_issr_cont
238REAL :: dqvc_mix_sub_cont, dqvc_mix_issr_cont
239REAL :: dcf_mix_cont, dqvc_mix_cont, dqi_mix_cont, dqt_mix_cont
240REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix
241REAL :: envfra_mix, cldfra_mix
242REAL :: L_shear, shear_fra
243REAL :: sigma_mix, issrfra_mix, subfra_mix
244REAL :: qvapincld, qvapinmix, qvapincld_new, qiceincld
245
246!--Initialisation
247dcf_mix_sub_cont   = 0.
248dqt_mix_sub_cont   = 0.
249dqvc_mix_sub_cont  = 0.
250dcf_mix_issr_cont  = 0.
251dqt_mix_issr_cont  = 0.
252dqvc_mix_issr_cont = 0.
253
254!-- PART 1 - TURBULENT DIFFUSION
255
256!--Clouds within the mesh are assumed to be ellipses. The length of the
257!--semi-major axis is a and the length of the semi-minor axis is b.
258!--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer.
259!--In particular, it is 0 if cldfra = 1.
260!--clouds_perim is the total perimeter of the clouds within the mesh,
261!--not considering interfaces with other meshes (only the interfaces with clear
262!--sky are taken into account).
263!--
264!--The area of each cloud is A = a * b * RPI,
265!--and the perimeter of each cloud is
266!-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
267!--
268!--With cell_area the area of the cell, we have:
269!-- cldfra = A * N_cld_mix / cell_area
270!-- clouds_perim = P * N_cld_mix
271!--
272!--We assume that the ratio between b and a is a function of
273!--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
274!--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
275!--are spherical.
276!-- b / a = bovera = MAX(0.1, cldfra)
277bovera = aspect_ratio_contrails
278!--P / a is a function of b / a only, that we can calculate
279!-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
280Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )
281!--The clouds perimeter is imposed using the formula from Morcrette 2012,
282!--based on observations.
283!-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra )
284!--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have:
285!-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra )
286!-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) )
287a_mix = Povera / coef_mixing_contrails / RPI / ( 1. - cldfra ) / bovera
288!--and finally,
289!-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a )
290N_cld_mix = coef_mixing_contrails * cldfra * ( 1. - cldfra ) * cell_area &
291          / Povera / a_mix
292
293!--The time required for turbulent diffusion to homogenize a region of size
294!--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
295!--We compute L_mix and assume that the cloud is mixed over this length
296L_mix = SQRT( dtime**3 * pbl_eps )
297!--The mixing length cannot be greater than the semi-minor axis. In this case,
298!--the entire cloud is mixed.
299L_mix = MIN(L_mix, a_mix * bovera)
300
301!--The fraction of clear sky mixed is
302!-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
303envfra_mix = N_cld_mix * RPI / cell_area &
304           * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
305!--The fraction of cloudy sky mixed is
306!-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
307cldfra_mix = N_cld_mix * RPI / cell_area &
308           * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
309
310
311!-- PART 2 - SHEARING
312
313!--The clouds are then sheared. We keep the shape and number
314!--assumptions from before. The clouds are sheared along their
315!--semi-major axis (a_mix), on the entire cell heigh dz.
316!--The increase in size is
317L_shear = coef_shear_contrails * shear * dz * dtime
318!--therefore, the fraction of clear sky mixed is
319!-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
320!--and the fraction of cloud mixed is
321!-- N_cld_mix * ( (a * b) - (a - L_shear) * b ) * RPI / 2. / cell_area
322!--(note that they are equal)
323shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area
324!--and the environment and cloud mixed fractions are the same,
325!--which we add to the previous calculated mixed fractions.
326!--We therefore assume that the sheared clouds and the turbulent
327!--mixed clouds are different.
328envfra_mix = envfra_mix + shear_fra
329cldfra_mix = cldfra_mix + shear_fra
330
331
332!-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
333
334!--The environment fraction is allocated to subsaturated sky or supersaturated sky,
335!--according to the factor sigma_mix. This is computed as the ratio of the
336!--subsaturated sky fraction to the environment fraction, corrected by a factor
337!--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the
338!--supersaturated sky is favoured. Physically, this means that it is more likely
339!--to have supersaturated sky around the cloud than subsaturated sky.
340sigma_mix = subfra / ( subfra + chi_mixing_contrails * issrfra )
341subfra_mix = MIN( sigma_mix * envfra_mix, subfra )
342issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra )
343cldfra_mix = MIN( cldfra_mix, cldfra )
344
345!--First, we mix the subsaturated sky (subfra_mix) and the cloud close
346!--to this fraction (sigma_mix * cldfra_mix).
347IF ( subfra .GT. eps ) THEN
348  !--We compute the total humidity in the mixed air, which
349  !--can be either sub- or supersaturated.
350  qvapinmix = ( qsub * subfra_mix / subfra &
351            + qcld * cldfra_mix * sigma_mix / cldfra ) &
352            / ( subfra_mix + cldfra_mix * sigma_mix )
353
354  IF ( ok_unadjusted_clouds ) THEN
355    qiceincld = ( qcld - qvc ) * cldfra_mix * sigma_mix / cldfra &
356              / ( subfra_mix + cldfra_mix * sigma_mix )
357    qvapincld_new = QVAPINCLD_DEPSUB_CONTRAILS( &
358        qvapinmix, qiceincld, temp, qsat, pplay, dtime)
359    IF ( qvapincld_new .EQ. 0. ) THEN
360      !--If all the ice has been sublimated, we sublimate
361      !--completely the cloud and do not activate the sublimation
362      !--process
363      !--Tendencies and diagnostics
364      dcf_mix_sub_cont = - sigma_mix * cldfra_mix
365      dqt_mix_sub_cont = dcf_mix_sub_cont * qcld / cldfra
366      dqvc_mix_sub_cont = dcf_mix_sub_cont * qvc / cldfra
367    ELSE
368      dcf_mix_sub_cont = subfra_mix
369      dqt_mix_sub_cont = dcf_mix_sub_cont * qsub / subfra
370      dqvc_mix_sub_cont = dcf_mix_sub_cont * qvapincld_new
371    ENDIF
372  ELSE
373    IF ( qvapinmix .GT. qsat ) THEN
374      !--If the mixed air is supersaturated, we condense the subsaturated
375      !--region which was mixed.
376      dcf_mix_sub_cont = subfra_mix
377      dqt_mix_sub_cont = dcf_mix_sub_cont * qsub / subfra
378      dqvc_mix_sub_cont = dcf_mix_sub_cont * qsat
379    ELSE
380      !--Else, we sublimate the cloud which was mixed.
381      dcf_mix_sub_cont = - sigma_mix * cldfra_mix
382      dqt_mix_sub_cont = dcf_mix_sub_cont * qcld / cldfra
383      dqvc_mix_sub_cont = dcf_mix_sub_cont * qsat
384    ENDIF
385  ENDIF ! ok_unadjusted_clouds
386ENDIF ! subfra .GT. eps
387
388!--We then mix the supersaturated sky (issrfra_mix) and the cloud,
389!--for which the mixed air is always supersatured, therefore
390!--the cloud necessarily expands
391IF ( issrfra .GT. eps ) THEN
392
393  IF ( ok_unadjusted_clouds ) THEN
394    qvapinmix = ( qissr * issrfra_mix / issrfra &
395              + qcld * cldfra_mix * (1. - sigma_mix) / cldfra ) &
396              / ( issrfra_mix + cldfra_mix * (1. -  sigma_mix) )
397    qiceincld = ( qcld - qvc ) * cldfra_mix * (1. - sigma_mix) / cldfra &
398              / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) )
399    qvapincld_new = QVAPINCLD_DEPSUB_CONTRAILS( &
400        qvapinmix, qiceincld, temp, qsat, pplay, dtime)
401    dcf_mix_issr_cont = issrfra_mix
402    dqt_mix_issr_cont = dcf_mix_issr_cont * qissr / issrfra
403    dqvc_mix_issr_cont = dcf_mix_issr_cont * qvapincld_new
404  ELSE
405    !--In this case, the additionnal vapor condenses
406    dcf_mix_issr_cont = issrfra_mix
407    dqt_mix_issr_cont = dcf_mix_issr_cont * qissr / issrfra
408    dqvc_mix_issr_cont = dcf_mix_issr_cont * qsat
409  ENDIF ! ok_unadjusted_clouds
410
411ENDIF ! issrfra .GT. eps
412
413!--Sum up the tendencies from subsaturated sky and supersaturated sky
414dcf_mix_cont  = dcf_mix_sub_cont  + dcf_mix_issr_cont
415dqt_mix_cont  = dqt_mix_sub_cont  + dqt_mix_issr_cont
416dqvc_mix_cont = dqvc_mix_sub_cont + dqvc_mix_issr_cont
417dqi_mix_cont  = dqt_mix_cont - dqvc_mix_cont
418
419!--Outputs
420!--The mixing tendencies are a linear combination of the calculation done for "classical" cirrus
421!--and contrails
422dcf_mix = ( 1. - rcont_seri ) * dcf_mix + rcont_seri * dcf_mix_cont
423dqvc_mix = ( 1. - rcont_seri ) * dqvc_mix + rcont_seri * dqvc_mix_cont
424dqi_mix = ( 1. - rcont_seri ) * dqi_mix + rcont_seri * dqi_mix_cont
425dqt_mix = ( 1. - rcont_seri ) * dqt_mix + rcont_seri * dqt_mix_cont
426dcf_mix_issr = ( 1. - rcont_seri ) * dcf_mix_issr + rcont_seri * dcf_mix_issr_cont
427dqt_mix_issr = ( 1. - rcont_seri ) * dqt_mix_issr + rcont_seri * dqt_mix_issr_cont
428
429END SUBROUTINE contrails_mixing
430
431
432!**********************************************************************************
433FUNCTION qvapincld_depsub_contrails( &
434    qvapincld, qiceincld, temp, qsat, pplay, dtime)
435
436USE lmdz_lscp_ini,        ONLY: RV, RLSTT, RTT, EPS_W
437USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
438USE lmdz_lscp_ini,        ONLY: rm_ice_crystals_contrails
439
440IMPLICIT NONE
441
442! inputs
443REAL :: qvapincld     !
444REAL :: qiceincld     !
445REAL :: temp          !
446REAL :: qsat          !
447REAL :: pplay         !
448REAL :: dtime         ! time step [s]
449! output
450REAL :: qvapincld_depsub_contrails
451! local
452REAL :: qice_ratio
453REAL :: pres_sat, rho, kappa
454REAL :: air_thermal_conduct, water_vapor_diff
455REAL :: rm_ice
456
457!--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
458!--hypothesis is lost, and the vapor in the cloud is purely prognostic.
459!
460!--The deposition equation is
461!-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
462!--from Lohmann et al. (2016), where
463!--alpha is the deposition coefficient [-]
464!--mi is the mass of one ice crystal [kg]
465!--C is the capacitance of an ice crystal [m]
466!--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
467!--R_v is the specific gas constant for humid air [J/kg/K]
468!--T is the temperature [K]
469!--esi is the saturation pressure w.r.t. ice [Pa]
470!--Dv is the diffusivity of water vapor [m2/s]
471!--Ls is the specific latent heat of sublimation [J/kg/K]
472!--ka is the thermal conductivity of dry air [J/m/s/K]
473!
474!--alpha is a coefficient to take into account the fact that during deposition, a water
475!--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
476!--coherent (with the same structure). It has no impact for sublimation.
477!--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
478!--during deposition, and alpha = 1. during sublimation.
479!--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
480!-- C = capa_cond_cirrus * rm_ice
481!
482!--We have qice = Nice * mi, where Nice is the ice crystal
483!--number concentration per kg of moist air
484!--HYPOTHESIS 1: the ice crystals are spherical, therefore
485!-- mi = 4/3 * pi * rm_ice**3 * rho_ice
486!--HYPOTHESIS 2: the ice crystals are monodisperse with the
487!--initial radius rm_ice_0.
488!--NB. this is notably different than the assumption
489!--of a distributed qice in the cloud made in the sublimation process;
490!--should it be consistent?
491!
492!--As the deposition process does not create new ice crystals,
493!--and because we assume a same rm_ice value for all crystals
494!--therefore the sublimation process does not destroy ice crystals
495!--(or, in a limit case, it destroys all ice crystals), then
496!--Nice is a constant during the sublimation/deposition process.
497!-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
498!
499!--The deposition equation then reads:
500!-- 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
501!-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
502!--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
503!--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
504!-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
505!--  *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 )
506!--and we have
507!-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
508!-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
509!--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
510!
511!--This system of equations can be resolved with an exact
512!--explicit numerical integration, having one variable resolved
513!--explicitly, the other exactly. The exactly resolved variable is
514!--the one decreasing, so it is qvc if the process is
515!--condensation, qi if it is sublimation.
516!
517!--kappa is computed as an initialisation constant, as it depends only
518!--on temperature and other pre-computed values
519pres_sat = qsat / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay
520!--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
521air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184
522!--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
523water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4
524kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat &
525      * ( RV * temp / water_vapor_diff / pres_sat &
526        + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) )
527!--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
528
529!--Here, the initial vapor in the cloud is qvapincld, and we compute
530!--the new vapor qvapincld_depsub_contrails
531
532rm_ice = rm_ice_crystals_contrails
533
534IF ( qvapincld .GE. qsat ) THEN
535  !--If the cloud is initially supersaturated
536  !--Exact explicit integration (qvc exact, qice explicit)
537  qvapincld_depsub_contrails = qsat + ( qvapincld - qsat ) &
538                * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
539ELSE
540  !--If the cloud is initially subsaturated
541  !--Exact explicit integration (qice exact, qvc explicit)
542  !--The barrier is set so that the resulting vapor in cloud
543  !--cannot be greater than qsat
544  !--qice_ratio is the ratio between the new ice content and
545  !--the old one, it is comprised between 0 and 1
546  qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2 * dtime * ( qsat - qvapincld ) )
547
548  IF ( qice_ratio .LT. 0. ) THEN
549    !--The new vapor in cloud is set to 0 so that the
550    !--sublimation process does not activate
551    qvapincld_depsub_contrails = 0.
552  ELSE
553    !--Else, the sublimation process is activated with the
554    !--diagnosed new cloud water vapor
555    !--The new vapor in the cloud is increased with the
556    !--sublimated ice
557    qvapincld_depsub_contrails = qvapincld + qiceincld * ( 1. - qice_ratio**1.5 )
558    !--The new vapor in the cloud cannot be greater than qsat
559    qvapincld_depsub_contrails = MIN(qvapincld_depsub_contrails, qsat)
560  ENDIF ! qice_ratio .LT. 0.
561ENDIF ! qvapincld .GT. qsat
562
563END FUNCTION qvapincld_depsub_contrails
564
565
566!**********************************************************************************
567FUNCTION contrail_cross_section_onera()
568
569IMPLICIT NONE
570
571! input
572! output
573REAL :: contrail_cross_section_onera ! [m2]
574! local
575
576contrail_cross_section_onera = 200. * 200.
577
578END FUNCTION contrail_cross_section_onera
579
580END MODULE lmdz_aviation
581
582!*******************************************************************
583!
584!SUBROUTINE airplane(debut,pphis,pplay,paprs,t_seri)
585!
586!  USE dimphy
587!  USE mod_grid_phy_lmdz,  ONLY: klon_glo
588!  USE geometry_mod, ONLY: cell_area
589!  USE phys_cal_mod, ONLY : mth_cur
590!  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
591!  USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
592!  USE mod_phys_lmdz_para, ONLY: scatter, bcast
593!  USE print_control_mod, ONLY: lunout
594!
595!  IMPLICIT NONE
596!
597!  INCLUDE "YOMCST.h"
598!  INCLUDE 'netcdf.inc'
599!
600!  !--------------------------------------------------------
601!  !--input variables
602!  !--------------------------------------------------------
603!  LOGICAL, INTENT(IN) :: debut
604!  REAL, INTENT(IN)    :: pphis(klon), pplay(klon,klev), paprs(klon,klev+1), t_seri(klon,klev)
605!
606!  !--------------------------------------------------------
607!  !    ... Local variables
608!  !--------------------------------------------------------
609!
610!  CHARACTER (LEN=20) :: modname='airplane_mod'
611!  INTEGER :: i, k, kori, iret, varid, error, ncida, klona
612!  INTEGER,SAVE :: nleva, ntimea
613!!$OMP THREADPRIVATE(nleva,ntimea)
614!  REAL, ALLOCATABLE :: pkm_airpl_glo(:,:,:)    !--km/s
615!  REAL, ALLOCATABLE :: ph2o_airpl_glo(:,:,:)   !--molec H2O/cm3/s
616!  REAL, ALLOCATABLE, SAVE :: zmida(:), zinta(:)
617!  REAL, ALLOCATABLE, SAVE :: pkm_airpl(:,:,:)
618!  REAL, ALLOCATABLE, SAVE :: ph2o_airpl(:,:,:)
619!!$OMP THREADPRIVATE(pkm_airpl,ph2o_airpl,zmida,zinta)
620!  REAL :: zalt(klon,klev+1)
621!  REAL :: zrho, zdz(klon,klev), zfrac
622!
623!  !
624!  IF (debut) THEN
625!  !--------------------------------------------------------------------------------
626!  !       ... Open the file and read airplane emissions
627!  !--------------------------------------------------------------------------------
628!  !
629!  IF (is_mpi_root .AND. is_omp_root) THEN
630!      !
631!      iret = nf_open('aircraft_phy.nc', 0, ncida)
632!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to open aircraft_phy.nc file',1)
633!      ! ... Get lengths
634!      iret = nf_inq_dimid(ncida, 'time', varid)
635!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimid in aircraft_phy.nc file',1)
636!      iret = nf_inq_dimlen(ncida, varid, ntimea)
637!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimlen aircraft_phy.nc file',1)
638!      iret = nf_inq_dimid(ncida, 'vector', varid)
639!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimid aircraft_phy.nc file',1)
640!      iret = nf_inq_dimlen(ncida, varid, klona)
641!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimlen aircraft_phy.nc file',1)
642!      iret = nf_inq_dimid(ncida, 'lev', varid)
643!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1)
644!      iret = nf_inq_dimlen(ncida, varid, nleva)
645!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimlen aircraft_phy.nc file',1)
646!      !
647!      IF ( klona /= klon_glo ) THEN
648!        WRITE(lunout,*) 'klona & klon_glo =', klona, klon_glo
649!        CALL abort_physic(modname,'problem klon in aircraft_phy.nc file',1)
650!      ENDIF
651!      !
652!      IF ( ntimea /= 12 ) THEN
653!        WRITE(lunout,*) 'ntimea=', ntimea
654!        CALL abort_physic(modname,'problem ntime<>12 in aircraft_phy.nc file',1)
655!      ENDIF
656!      !
657!      ALLOCATE(zmida(nleva), STAT=error)
658!      IF (error /= 0) CALL abort_physic(modname,'problem to allocate zmida',1)
659!      ALLOCATE(pkm_airpl_glo(klona,nleva,ntimea), STAT=error)
660!      IF (error /= 0) CALL abort_physic(modname,'problem to allocate pkm_airpl_glo',1)
661!      ALLOCATE(ph2o_airpl_glo(klona,nleva,ntimea), STAT=error)
662!      IF (error /= 0) CALL abort_physic(modname,'problem to allocate ph2o_airpl_glo',1)
663!      !
664!      iret = nf_inq_varid(ncida, 'lev', varid)
665!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1)
666!      iret = nf_get_var_double(ncida, varid, zmida)
667!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read zmida file',1)
668!      !
669!      iret = nf_inq_varid(ncida, 'emi_co2_aircraft', varid)  !--CO2 as a proxy for m flown -
670!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_distance dimid aircraft_phy.nc file',1)
671!      iret = nf_get_var_double(ncida, varid, pkm_airpl_glo)
672!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read pkm_airpl file',1)
673!      !
674!      iret = nf_inq_varid(ncida, 'emi_h2o_aircraft', varid)
675!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_h2o_aircraft dimid aircraft_phy.nc file',1)
676!      iret = nf_get_var_double(ncida, varid, ph2o_airpl_glo)
677!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read ph2o_airpl file',1)
678!      !
679!     ENDIF    !--is_mpi_root and is_omp_root
680!     !
681!     CALL bcast(nleva)
682!     CALL bcast(ntimea)
683!     !
684!     IF (.NOT.ALLOCATED(zmida)) ALLOCATE(zmida(nleva), STAT=error)
685!     IF (.NOT.ALLOCATED(zinta)) ALLOCATE(zinta(nleva+1), STAT=error)
686!     !
687!     ALLOCATE(pkm_airpl(klon,nleva,ntimea))
688!     ALLOCATE(ph2o_airpl(klon,nleva,ntimea))
689!     !
690!     ALLOCATE(flight_m(klon,klev))
691!     ALLOCATE(flight_h2o(klon,klev))
692!     !
693!     CALL bcast(zmida)
694!     zinta(1)=0.0                                   !--surface
695!     DO k=2, nleva
696!       zinta(k) = (zmida(k-1)+zmida(k))/2.0*1000.0  !--conversion from km to m
697!     ENDDO
698!     zinta(nleva+1)=zinta(nleva)+(zmida(nleva)-zmida(nleva-1))*1000.0 !--extrapolation for last interface
699!     !print *,'zinta=', zinta
700!     !
701!     CALL scatter(pkm_airpl_glo,pkm_airpl)
702!     CALL scatter(ph2o_airpl_glo,ph2o_airpl)
703!     !
704!!$OMP MASTER
705!     IF (is_mpi_root .AND. is_omp_root) THEN
706!        DEALLOCATE(pkm_airpl_glo)
707!        DEALLOCATE(ph2o_airpl_glo)
708!     ENDIF   !--is_mpi_root
709!!$OMP END MASTER
710!
711!  ENDIF !--debut
712!!
713!!--compute altitude of model level interfaces
714!!
715!  DO i = 1, klon
716!    zalt(i,1)=pphis(i)/RG         !--in m
717!  ENDDO
718!!
719!  DO k=1, klev
720!    DO i = 1, klon
721!      zrho=pplay(i,k)/t_seri(i,k)/RD
722!      zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho/RG
723!      zalt(i,k+1)=zalt(i,k)+zdz(i,k)   !--in m
724!    ENDDO
725!  ENDDO
726!!
727!!--vertical reprojection
728!!
729!  flight_m(:,:)=0.0
730!  flight_h2o(:,:)=0.0
731!!
732!  DO k=1, klev
733!    DO kori=1, nleva
734!      DO i=1, klon
735!        !--fraction of layer kori included in layer k
736!        zfrac=max(0.0,min(zalt(i,k+1),zinta(kori+1))-max(zalt(i,k),zinta(kori)))/(zinta(kori+1)-zinta(kori))
737!        !--reproject
738!        flight_m(i,k)=flight_m(i,k) + pkm_airpl(i,kori,mth_cur)*zfrac
739!        !--reproject
740!        flight_h2o(i,k)=flight_h2o(i,k) + ph2o_airpl(i,kori,mth_cur)*zfrac   
741!      ENDDO
742!    ENDDO
743!  ENDDO
744!!
745!  DO k=1, klev
746!    DO i=1, klon
747!      !--molec.cm-3.s-1 / (molec/mol) * kg CO2/mol * m2 * m * cm3/m3 / (kg CO2/m) => m s-1 per cell
748!      flight_m(i,k)=flight_m(i,k)/RNAVO*44.e-3*cell_area(i)*zdz(i,k)*1.e6/16.37e-3
749!      flight_m(i,k)=flight_m(i,k)*100.0  !--x100 to augment signal to noise
750!      !--molec.cm-3.s-1 / (molec/mol) * kg H2O/mol * m2 * m * cm3/m3 => kg H2O s-1 per cell
751!      flight_h2o(i,k)=flight_h2o(i,k)/RNAVO*18.e-3*cell_area(i)*zdz(i,k)*1.e6
752!    ENDDO
753!  ENDDO
754!!
755!END SUBROUTINE airplane
756!
757!!********************************************************************
758!! simple routine to initialise flight_m and test a flight corridor
759!!--Olivier Boucher - 2021
760!!
761!SUBROUTINE flight_init()
762!  USE dimphy
763!  USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
764!  IMPLICIT NONE
765!  INTEGER :: i
766!
767!  ALLOCATE(flight_m(klon,klev))
768!  ALLOCATE(flight_h2o(klon,klev))
769!  !
770!  flight_m(:,:) = 0.0    !--initialisation
771!  flight_h2o(:,:) = 0.0  !--initialisation
772!  !
773!  DO i=1, klon
774!   IF (latitude_deg(i).GE.42.0.AND.latitude_deg(i).LE.48.0) THEN
775!     flight_m(i,38) = 50000.0  !--5000 m of flight/second in grid cell x 10 scaling
776!   ENDIF
777!  ENDDO
778! 
779!  RETURN
780!END SUBROUTINE flight_init
781!
782!     !--add a source of cirrus from aviation contrails
783!     IF (ok_plane_contrail) THEN
784!       drneb_avi(i,k) = rnebss*flight_m(i,k)*contrail_cross_section/V_cell    !--tendency rneb due to aviation [s-1]
785!       drneb_avi(i,k) = MIN(drneb_avi(i,k), rnebss/dtime)                     !--majoration
786!       dqss_avi = qss*drneb_avi(i,k)/MAX(eps,rnebss)                          !--tendency q aviation [kg kg-1 s-1]
787!       rneb = rneb + drneb_avi(i,k)*dtime                                     !--add tendency to rneb
788!       qcld = qcld + dqss_avi*dtime                                           !--add tendency to qcld
789!       rnebss = rnebss - drneb_avi(i,k)*dtime                                 !--add tendency to rnebss
790!       qss = qss - dqss_avi*dtime                                             !--add tendency to qss
791!     ELSE
792!       drneb_avi(i,k)=0.0
793!     ENDIF
794!
795!  RETURN
796!END SUBROUTINE ice_sursat
Note: See TracBrowser for help on using the repository browser.