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

Last change on this file since 5456 was 5456, checked in by aborella, 37 hours ago

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