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

Last change on this file since 5467 was 5466, checked in by aborella, 3 days ago

Modification to the shearing process
File size: 35.0 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-Appleman 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 with a random orientation
346!--of the wind, on average we assume that the wind and the semi-major axis
347!--make a 45 degrees angle. Moreover, the contrails only mix
348!--along their semi-minor axis (b), because it is easier to compute.
349!--With this, the clouds increase in size along b only, by a factor
350!--L_shear * SQRT(2.) / 2. (to account for the 45 degrees orientation of the wind)
351L_shear = coef_shear_contrails * shear * dz * dtime
352!--therefore, the fraction of clear sky mixed is
353!-- N_cld_mix * ( a * (b + L_shear * SQRT(2.) / 2.) - a * b ) * RPI / 2. / cell_area
354!--and the fraction of cloud mixed is
355!-- N_cld_mix * ( a * b - a * (b - L_shear * SQRT(2.) / 2.) ) * RPI / 2. / cell_area
356!--(note that they are equal)
357shear_fra = RPI * SQRT(2.) / 2. * L_shear * a_mix / 2. * N_cld_mix / cell_area
358!--and the environment and cloud mixed fractions are the same,
359!--which we add to the previous calculated mixed fractions.
360!--We therefore assume that the sheared clouds and the turbulent
361!--mixed clouds are different.
362envfra_mix = envfra_mix + shear_fra
363cldfra_mix = cldfra_mix + shear_fra
364
365
366!-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
367
368!--The environment fraction is allocated to subsaturated sky or supersaturated sky,
369!--according to the factor sigma_mix. This is computed as the ratio of the
370!--subsaturated sky fraction to the environment fraction, corrected by a factor
371!--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the
372!--supersaturated sky is favoured. Physically, this means that it is more likely
373!--to have supersaturated sky around the cloud than subsaturated sky.
374sigma_mix = subfra / ( subfra + chi_mixing_contrails * issrfra )
375subfra_mix = MIN( sigma_mix * envfra_mix, subfra )
376issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra )
377cldfra_mix = MIN( cldfra_mix, cldfra )
378
379!--First, we mix the subsaturated sky (subfra_mix) and the cloud close
380!--to this fraction (sigma_mix * cldfra_mix).
381IF ( subfra .GT. eps ) THEN
382  !--We compute the total humidity in the mixed air, which
383  !--can be either sub- or supersaturated.
384  qvapinmix = ( qsub * subfra_mix / subfra &
385            + qcld * cldfra_mix * sigma_mix / cldfra ) &
386            / ( subfra_mix + cldfra_mix * sigma_mix )
387
388  IF ( ok_unadjusted_clouds ) THEN
389    qiceincld = ( qcld - qvc ) * cldfra_mix * sigma_mix / cldfra &
390              / ( subfra_mix + cldfra_mix * sigma_mix )
391    qvapincld_new = QVAPINCLD_DEPSUB_CONTRAILS( &
392        qvapinmix, qiceincld, temp, qsat, pplay, dtime)
393    IF ( qvapincld_new .EQ. 0. ) THEN
394      !--If all the ice has been sublimated, we sublimate
395      !--completely the cloud and do not activate the sublimation
396      !--process
397      !--Tendencies and diagnostics
398      dcf_mix_sub_cont = - sigma_mix * cldfra_mix
399      dqt_mix_sub_cont = dcf_mix_sub_cont * qcld / cldfra
400      dqvc_mix_sub_cont = dcf_mix_sub_cont * qvc / cldfra
401    ELSE
402      dcf_mix_sub_cont = subfra_mix
403      dqt_mix_sub_cont = dcf_mix_sub_cont * qsub / subfra
404      dqvc_mix_sub_cont = dcf_mix_sub_cont * qvapincld_new
405    ENDIF
406  ELSE
407    IF ( qvapinmix .GT. qsat ) THEN
408      !--If the mixed air is supersaturated, we condense the subsaturated
409      !--region which was mixed.
410      dcf_mix_sub_cont = subfra_mix
411      dqt_mix_sub_cont = dcf_mix_sub_cont * qsub / subfra
412      dqvc_mix_sub_cont = dcf_mix_sub_cont * qsat
413    ELSE
414      !--Else, we sublimate the cloud which was mixed.
415      dcf_mix_sub_cont = - sigma_mix * cldfra_mix
416      dqt_mix_sub_cont = dcf_mix_sub_cont * qcld / cldfra
417      dqvc_mix_sub_cont = dcf_mix_sub_cont * qsat
418    ENDIF
419  ENDIF ! ok_unadjusted_clouds
420ENDIF ! subfra .GT. eps
421
422!--We then mix the supersaturated sky (issrfra_mix) and the cloud,
423!--for which the mixed air is always supersatured, therefore
424!--the cloud necessarily expands
425IF ( issrfra .GT. eps ) THEN
426
427  IF ( ok_unadjusted_clouds ) THEN
428    qvapinmix = ( qissr * issrfra_mix / issrfra &
429              + qcld * cldfra_mix * (1. - sigma_mix) / cldfra ) &
430              / ( issrfra_mix + cldfra_mix * (1. -  sigma_mix) )
431    qiceincld = ( qcld - qvc ) * cldfra_mix * (1. - sigma_mix) / cldfra &
432              / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) )
433    qvapincld_new = QVAPINCLD_DEPSUB_CONTRAILS( &
434        qvapinmix, qiceincld, temp, qsat, pplay, dtime)
435    dcf_mix_issr_cont = issrfra_mix
436    dqt_mix_issr_cont = dcf_mix_issr_cont * qissr / issrfra
437    dqvc_mix_issr_cont = dcf_mix_issr_cont * qvapincld_new
438  ELSE
439    !--In this case, the additionnal vapor condenses
440    dcf_mix_issr_cont = issrfra_mix
441    dqt_mix_issr_cont = dcf_mix_issr_cont * qissr / issrfra
442    dqvc_mix_issr_cont = dcf_mix_issr_cont * qsat
443  ENDIF ! ok_unadjusted_clouds
444
445ENDIF ! issrfra .GT. eps
446
447!--Sum up the tendencies from subsaturated sky and supersaturated sky
448dcf_mix_cont  = dcf_mix_sub_cont  + dcf_mix_issr_cont
449dqt_mix_cont  = dqt_mix_sub_cont  + dqt_mix_issr_cont
450dqvc_mix_cont = dqvc_mix_sub_cont + dqvc_mix_issr_cont
451dqi_mix_cont  = dqt_mix_cont - dqvc_mix_cont
452
453!--Outputs
454!--The mixing tendencies are a linear combination of the calculation done for "classical" cirrus
455!--and contrails
456dcf_mix = ( 1. - rcont_seri ) * dcf_mix + rcont_seri * dcf_mix_cont
457dqvc_mix = ( 1. - rcont_seri ) * dqvc_mix + rcont_seri * dqvc_mix_cont
458dqi_mix = ( 1. - rcont_seri ) * dqi_mix + rcont_seri * dqi_mix_cont
459dqt_mix = ( 1. - rcont_seri ) * dqt_mix + rcont_seri * dqt_mix_cont
460dcf_mix_issr = ( 1. - rcont_seri ) * dcf_mix_issr + rcont_seri * dcf_mix_issr_cont
461dqt_mix_issr = ( 1. - rcont_seri ) * dqt_mix_issr + rcont_seri * dqt_mix_issr_cont
462
463END SUBROUTINE contrails_mixing
464
465
466!**********************************************************************************
467FUNCTION qvapincld_depsub_contrails( &
468    qvapincld, qiceincld, temp, qsat, pplay, dtime)
469
470USE lmdz_lscp_ini,        ONLY: RV, RLSTT, RTT, EPS_W
471USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
472USE lmdz_lscp_ini,        ONLY: rm_ice_crystals_contrails
473
474IMPLICIT NONE
475
476! inputs
477REAL :: qvapincld     !
478REAL :: qiceincld     !
479REAL :: temp          !
480REAL :: qsat          !
481REAL :: pplay         !
482REAL :: dtime         ! time step [s]
483! output
484REAL :: qvapincld_depsub_contrails
485! local
486REAL :: qice_ratio
487REAL :: pres_sat, rho, kappa
488REAL :: air_thermal_conduct, water_vapor_diff
489REAL :: rm_ice
490
491!--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
492!--hypothesis is lost, and the vapor in the cloud is purely prognostic.
493!
494!--The deposition equation is
495!-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
496!--from Lohmann et al. (2016), where
497!--alpha is the deposition coefficient [-]
498!--mi is the mass of one ice crystal [kg]
499!--C is the capacitance of an ice crystal [m]
500!--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
501!--R_v is the specific gas constant for humid air [J/kg/K]
502!--T is the temperature [K]
503!--esi is the saturation pressure w.r.t. ice [Pa]
504!--Dv is the diffusivity of water vapor [m2/s]
505!--Ls is the specific latent heat of sublimation [J/kg/K]
506!--ka is the thermal conductivity of dry air [J/m/s/K]
507!
508!--alpha is a coefficient to take into account the fact that during deposition, a water
509!--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
510!--coherent (with the same structure). It has no impact for sublimation.
511!--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
512!--during deposition, and alpha = 1. during sublimation.
513!--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
514!-- C = capa_cond_cirrus * rm_ice
515!
516!--We have qice = Nice * mi, where Nice is the ice crystal
517!--number concentration per kg of moist air
518!--HYPOTHESIS 1: the ice crystals are spherical, therefore
519!-- mi = 4/3 * pi * rm_ice**3 * rho_ice
520!--HYPOTHESIS 2: the ice crystals are monodisperse with the
521!--initial radius rm_ice_0.
522!--NB. this is notably different than the assumption
523!--of a distributed qice in the cloud made in the sublimation process;
524!--should it be consistent?
525!
526!--As the deposition process does not create new ice crystals,
527!--and because we assume a same rm_ice value for all crystals
528!--therefore the sublimation process does not destroy ice crystals
529!--(or, in a limit case, it destroys all ice crystals), then
530!--Nice is a constant during the sublimation/deposition process.
531!-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
532!
533!--The deposition equation then reads:
534!-- 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
535!-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
536!--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
537!--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
538!-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
539!--  *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 )
540!--and we have
541!-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
542!-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
543!--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
544!
545!--This system of equations can be resolved with an exact
546!--explicit numerical integration, having one variable resolved
547!--explicitly, the other exactly. The exactly resolved variable is
548!--the one decreasing, so it is qvc if the process is
549!--condensation, qi if it is sublimation.
550!
551!--kappa is computed as an initialisation constant, as it depends only
552!--on temperature and other pre-computed values
553pres_sat = qsat / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay
554!--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
555air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184
556!--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
557water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4
558kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat &
559      * ( RV * temp / water_vapor_diff / pres_sat &
560        + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) )
561!--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
562
563!--Here, the initial vapor in the cloud is qvapincld, and we compute
564!--the new vapor qvapincld_depsub_contrails
565
566rm_ice = rm_ice_crystals_contrails
567
568IF ( qvapincld .GE. qsat ) THEN
569  !--If the cloud is initially supersaturated
570  !--Exact explicit integration (qvc exact, qice explicit)
571  qvapincld_depsub_contrails = qsat + ( qvapincld - qsat ) &
572                * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
573ELSE
574  !--If the cloud is initially subsaturated
575  !--Exact explicit integration (qice exact, qvc explicit)
576  !--The barrier is set so that the resulting vapor in cloud
577  !--cannot be greater than qsat
578  !--qice_ratio is the ratio between the new ice content and
579  !--the old one, it is comprised between 0 and 1
580  qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2 * dtime * ( qsat - qvapincld ) )
581
582  IF ( qice_ratio .LT. 0. ) THEN
583    !--The new vapor in cloud is set to 0 so that the
584    !--sublimation process does not activate
585    qvapincld_depsub_contrails = 0.
586  ELSE
587    !--Else, the sublimation process is activated with the
588    !--diagnosed new cloud water vapor
589    !--The new vapor in the cloud is increased with the
590    !--sublimated ice
591    qvapincld_depsub_contrails = qvapincld + qiceincld * ( 1. - qice_ratio**1.5 )
592    !--The new vapor in the cloud cannot be greater than qsat
593    qvapincld_depsub_contrails = MIN(qvapincld_depsub_contrails, qsat)
594  ENDIF ! qice_ratio .LT. 0.
595ENDIF ! qvapincld .GT. qsat
596
597END FUNCTION qvapincld_depsub_contrails
598
599
600!**********************************************************************************
601FUNCTION contrail_cross_section_onera(dz)
602
603USE lmdz_lscp_ini, ONLY: initial_width_contrails
604
605IMPLICIT NONE
606
607!
608! Input
609!
610REAL :: dz ! cell width [m]
611!
612! Output
613!
614REAL :: contrail_cross_section_onera ! [m2]
615!
616! Local
617!
618
619contrail_cross_section_onera = initial_width_contrails * dz
620
621END FUNCTION contrail_cross_section_onera
622
623
624SUBROUTINE read_aviation_emissions( &
625        klon, klev, latitude_deg, longitude_deg, pplay, &
626        flight_dist, flight_h2o &
627        )
628
629IMPLICIT NONE
630!
631! Input
632!
633INTEGER, INTENT(IN)                        :: klon, klev       ! number of horizontal grid points and vertical levels
634REAL,    INTENT(IN),  DIMENSION(klon)      :: latitude_deg     ! latitude of the grid points [deg]
635REAL,    INTENT(IN),  DIMENSION(klon)      :: longitude_deg    ! longitude of the grid points [deg]
636REAL,    INTENT(IN),  DIMENSION(klon,klev) :: pplay            ! layer pressure [Pa]
637!
638! Output
639!
640REAL,    INTENT(OUT), DIMENSION(klon,klev) :: flight_dist      ! aviation distance flown within the mesh [m/s/mesh]
641REAL,    INTENT(OUT), DIMENSION(klon,klev) :: flight_h2o       ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
642!
643! Local
644!
645INTEGER :: i, k
646
647!--Initialisation
648flight_dist(:,:) = 0.
649flight_h2o(:,:) = 0.
650
651DO i=1, klon
652 IF ( ( latitude_deg(i) .GE. 42. ) .AND. ( latitude_deg(i) .LE. 48. ) ) THEN
653   flight_dist(i,15) = 50000.  !--5000 m of flight/second in grid cell x 10 scaling
654   flight_h2o(i,15) = 100.  !--10 kgH2O/second in grid cell x 10 scaling
655 ENDIF
656ENDDO
657
658END SUBROUTINE read_aviation_emissions
659
660END MODULE lmdz_aviation
661
662!*******************************************************************
663!
664!SUBROUTINE airplane(debut,pphis,pplay,paprs,t_seri)
665!
666!  USE dimphy
667!  USE mod_grid_phy_lmdz,  ONLY: klon_glo
668!  USE geometry_mod, ONLY: cell_area
669!  USE phys_cal_mod, ONLY : mth_cur
670!  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
671!  USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
672!  USE mod_phys_lmdz_para, ONLY: scatter, bcast
673!  USE print_control_mod, ONLY: lunout
674!
675!  IMPLICIT NONE
676!
677!  INCLUDE "YOMCST.h"
678!  INCLUDE 'netcdf.inc'
679!
680!  !--------------------------------------------------------
681!  !--input variables
682!  !--------------------------------------------------------
683!  LOGICAL, INTENT(IN) :: debut
684!  REAL, INTENT(IN)    :: pphis(klon), pplay(klon,klev), paprs(klon,klev+1), t_seri(klon,klev)
685!
686!  !--------------------------------------------------------
687!  !    ... Local variables
688!  !--------------------------------------------------------
689!
690!  CHARACTER (LEN=20) :: modname='airplane_mod'
691!  INTEGER :: i, k, kori, iret, varid, error, ncida, klona
692!  INTEGER,SAVE :: nleva, ntimea
693!!$OMP THREADPRIVATE(nleva,ntimea)
694!  REAL, ALLOCATABLE :: pkm_airpl_glo(:,:,:)    !--km/s
695!  REAL, ALLOCATABLE :: ph2o_airpl_glo(:,:,:)   !--molec H2O/cm3/s
696!  REAL, ALLOCATABLE, SAVE :: zmida(:), zinta(:)
697!  REAL, ALLOCATABLE, SAVE :: pkm_airpl(:,:,:)
698!  REAL, ALLOCATABLE, SAVE :: ph2o_airpl(:,:,:)
699!!$OMP THREADPRIVATE(pkm_airpl,ph2o_airpl,zmida,zinta)
700!  REAL :: zalt(klon,klev+1)
701!  REAL :: zrho, zdz(klon,klev), zfrac
702!
703!  !
704!  IF (debut) THEN
705!  !--------------------------------------------------------------------------------
706!  !       ... Open the file and read airplane emissions
707!  !--------------------------------------------------------------------------------
708!  !
709!  IF (is_mpi_root .AND. is_omp_root) THEN
710!      !
711!      iret = nf_open('aircraft_phy.nc', 0, ncida)
712!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to open aircraft_phy.nc file',1)
713!      ! ... Get lengths
714!      iret = nf_inq_dimid(ncida, 'time', varid)
715!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimid in aircraft_phy.nc file',1)
716!      iret = nf_inq_dimlen(ncida, varid, ntimea)
717!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimlen aircraft_phy.nc file',1)
718!      iret = nf_inq_dimid(ncida, 'vector', varid)
719!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimid aircraft_phy.nc file',1)
720!      iret = nf_inq_dimlen(ncida, varid, klona)
721!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimlen aircraft_phy.nc file',1)
722!      iret = nf_inq_dimid(ncida, 'lev', varid)
723!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1)
724!      iret = nf_inq_dimlen(ncida, varid, nleva)
725!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimlen aircraft_phy.nc file',1)
726!      !
727!      IF ( klona /= klon_glo ) THEN
728!        WRITE(lunout,*) 'klona & klon_glo =', klona, klon_glo
729!        CALL abort_physic(modname,'problem klon in aircraft_phy.nc file',1)
730!      ENDIF
731!      !
732!      IF ( ntimea /= 12 ) THEN
733!        WRITE(lunout,*) 'ntimea=', ntimea
734!        CALL abort_physic(modname,'problem ntime<>12 in aircraft_phy.nc file',1)
735!      ENDIF
736!      !
737!      ALLOCATE(zmida(nleva), STAT=error)
738!      IF (error /= 0) CALL abort_physic(modname,'problem to allocate zmida',1)
739!      ALLOCATE(pkm_airpl_glo(klona,nleva,ntimea), STAT=error)
740!      IF (error /= 0) CALL abort_physic(modname,'problem to allocate pkm_airpl_glo',1)
741!      ALLOCATE(ph2o_airpl_glo(klona,nleva,ntimea), STAT=error)
742!      IF (error /= 0) CALL abort_physic(modname,'problem to allocate ph2o_airpl_glo',1)
743!      !
744!      iret = nf_inq_varid(ncida, 'lev', varid)
745!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1)
746!      iret = nf_get_var_double(ncida, varid, zmida)
747!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read zmida file',1)
748!      !
749!      iret = nf_inq_varid(ncida, 'emi_co2_aircraft', varid)  !--CO2 as a proxy for m flown -
750!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_distance dimid aircraft_phy.nc file',1)
751!      iret = nf_get_var_double(ncida, varid, pkm_airpl_glo)
752!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read pkm_airpl file',1)
753!      !
754!      iret = nf_inq_varid(ncida, 'emi_h2o_aircraft', varid)
755!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_h2o_aircraft dimid aircraft_phy.nc file',1)
756!      iret = nf_get_var_double(ncida, varid, ph2o_airpl_glo)
757!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read ph2o_airpl file',1)
758!      !
759!     ENDIF    !--is_mpi_root and is_omp_root
760!     !
761!     CALL bcast(nleva)
762!     CALL bcast(ntimea)
763!     !
764!     IF (.NOT.ALLOCATED(zmida)) ALLOCATE(zmida(nleva), STAT=error)
765!     IF (.NOT.ALLOCATED(zinta)) ALLOCATE(zinta(nleva+1), STAT=error)
766!     !
767!     ALLOCATE(pkm_airpl(klon,nleva,ntimea))
768!     ALLOCATE(ph2o_airpl(klon,nleva,ntimea))
769!     !
770!     ALLOCATE(flight_m(klon,klev))
771!     ALLOCATE(flight_h2o(klon,klev))
772!     !
773!     CALL bcast(zmida)
774!     zinta(1)=0.0                                   !--surface
775!     DO k=2, nleva
776!       zinta(k) = (zmida(k-1)+zmida(k))/2.0*1000.0  !--conversion from km to m
777!     ENDDO
778!     zinta(nleva+1)=zinta(nleva)+(zmida(nleva)-zmida(nleva-1))*1000.0 !--extrapolation for last interface
779!     !print *,'zinta=', zinta
780!     !
781!     CALL scatter(pkm_airpl_glo,pkm_airpl)
782!     CALL scatter(ph2o_airpl_glo,ph2o_airpl)
783!     !
784!!$OMP MASTER
785!     IF (is_mpi_root .AND. is_omp_root) THEN
786!        DEALLOCATE(pkm_airpl_glo)
787!        DEALLOCATE(ph2o_airpl_glo)
788!     ENDIF   !--is_mpi_root
789!!$OMP END MASTER
790!
791!  ENDIF !--debut
792!!
793!!--compute altitude of model level interfaces
794!!
795!  DO i = 1, klon
796!    zalt(i,1)=pphis(i)/RG         !--in m
797!  ENDDO
798!!
799!  DO k=1, klev
800!    DO i = 1, klon
801!      zrho=pplay(i,k)/t_seri(i,k)/RD
802!      zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho/RG
803!      zalt(i,k+1)=zalt(i,k)+zdz(i,k)   !--in m
804!    ENDDO
805!  ENDDO
806!!
807!!--vertical reprojection
808!!
809!  flight_m(:,:)=0.0
810!  flight_h2o(:,:)=0.0
811!!
812!  DO k=1, klev
813!    DO kori=1, nleva
814!      DO i=1, klon
815!        !--fraction of layer kori included in layer k
816!        zfrac=max(0.0,min(zalt(i,k+1),zinta(kori+1))-max(zalt(i,k),zinta(kori)))/(zinta(kori+1)-zinta(kori))
817!        !--reproject
818!        flight_m(i,k)=flight_m(i,k) + pkm_airpl(i,kori,mth_cur)*zfrac
819!        !--reproject
820!        flight_h2o(i,k)=flight_h2o(i,k) + ph2o_airpl(i,kori,mth_cur)*zfrac   
821!      ENDDO
822!    ENDDO
823!  ENDDO
824!!
825!  DO k=1, klev
826!    DO i=1, klon
827!      !--molec.cm-3.s-1 / (molec/mol) * kg CO2/mol * m2 * m * cm3/m3 / (kg CO2/m) => m s-1 per cell
828!      flight_m(i,k)=flight_m(i,k)/RNAVO*44.e-3*cell_area(i)*zdz(i,k)*1.e6/16.37e-3
829!      flight_m(i,k)=flight_m(i,k)*100.0  !--x100 to augment signal to noise
830!      !--molec.cm-3.s-1 / (molec/mol) * kg H2O/mol * m2 * m * cm3/m3 => kg H2O s-1 per cell
831!      flight_h2o(i,k)=flight_h2o(i,k)/RNAVO*18.e-3*cell_area(i)*zdz(i,k)*1.e6
832!    ENDDO
833!  ENDDO
834!!
835!END SUBROUTINE airplane
Note: See TracBrowser for help on using the repository browser.