1 | MODULE lmdz_aviation |
---|
2 | !---------------------------------------------------------------- |
---|
3 | ! Module for aviation and contrails |
---|
4 | |
---|
5 | IMPLICIT NONE |
---|
6 | |
---|
7 | CONTAINS |
---|
8 | |
---|
9 | SUBROUTINE aviation_water_emissions( & |
---|
10 | klon, klev, dtime, paprs, pplay, temp, qtot, cell_area, & |
---|
11 | flight_h2o, d_q_avi & |
---|
12 | ) |
---|
13 | |
---|
14 | USE lmdz_lscp_ini, ONLY: RD, RG |
---|
15 | |
---|
16 | IMPLICIT NONE |
---|
17 | |
---|
18 | INTEGER, INTENT(IN) :: klon, klev ! number of horizontal grid points and vertical levels |
---|
19 | REAL, INTENT(IN) :: dtime ! time step [s] |
---|
20 | REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa] |
---|
21 | REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa] |
---|
22 | REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature (K) |
---|
23 | REAL, DIMENSION(klon,klev), INTENT(IN) :: qtot ! total specific humidity (in vapor phase) [kg/kg] |
---|
24 | REAL, DIMENSION(klon), INTENT(IN) :: cell_area ! area of each cell [m2] |
---|
25 | REAL, DIMENSION(klon,klev), INTENT(IN) :: flight_h2o ! aviation H2O emitted within the mesh [kgH2O/s/mesh] |
---|
26 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q_avi ! water vapor tendency from aviation [kg/kg] |
---|
27 | ! Local |
---|
28 | INTEGER :: i, k |
---|
29 | REAL :: rho, rhodz, dz, M_cell |
---|
30 | |
---|
31 | DO 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 |
---|
55 | ENDDO |
---|
56 | |
---|
57 | END SUBROUTINE aviation_water_emissions |
---|
58 | |
---|
59 | |
---|
60 | !********************************************************************************** |
---|
61 | SUBROUTINE 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 | |
---|
68 | USE lmdz_lscp_ini, ONLY: RCPD, EPS_W, RTT |
---|
69 | USE lmdz_lscp_ini, ONLY: EI_H2O_aviation, qheat_fuel_aviation, prop_efficiency_aviation |
---|
70 | USE lmdz_lscp_ini, ONLY: linear_contrails_lifetime |
---|
71 | USE lmdz_lscp_ini, ONLY: eps |
---|
72 | |
---|
73 | USE lmdz_lscp_tools, ONLY: GAMMAINC, calc_qsat_ecmwf |
---|
74 | |
---|
75 | IMPLICIT NONE |
---|
76 | |
---|
77 | ! |
---|
78 | ! Input |
---|
79 | ! |
---|
80 | REAL, INTENT(IN) :: dtime ! time step [s] |
---|
81 | REAL, INTENT(IN) :: pplay ! layer pressure [Pa] |
---|
82 | REAL, INTENT(IN) :: temp ! temperature [K] |
---|
83 | REAL, INTENT(IN) :: qsat ! saturation specific humidity [kg/kg] |
---|
84 | REAL, INTENT(IN) :: qsatl ! saturation specific humidity w.r.t. liquid [kg/kg] |
---|
85 | REAL, INTENT(IN) :: gamma_cond ! condensation threshold w.r.t. qsat [-] |
---|
86 | REAL, INTENT(IN) :: rcont_seri ! ratio of contrails fraction to total cloud fraction [-] |
---|
87 | REAL, INTENT(IN) :: flight_dist ! aviation distance flown within the mesh [m/s/mesh] |
---|
88 | REAL, INTENT(IN) :: cldfra ! cloud fraction [-] |
---|
89 | REAL, INTENT(IN) :: qvc ! gridbox-mean vapor in the cloud [kg/kg] |
---|
90 | REAL, INTENT(IN) :: dz ! cell width [m] |
---|
91 | REAL, INTENT(IN) :: V_cell ! cell volume [m3] |
---|
92 | REAL, INTENT(IN) :: pdf_loc ! location parameter of the clear sky PDF [%] |
---|
93 | REAL, INTENT(IN) :: pdf_scale ! scale parameter of the clear sky PDF [%] |
---|
94 | REAL, INTENT(IN) :: pdf_alpha ! shape parameter of the clear sky PDF [-] |
---|
95 | ! |
---|
96 | ! Output |
---|
97 | ! |
---|
98 | REAL, INTENT(OUT) :: Tcritcont ! critical temperature for contrail formation [K] |
---|
99 | REAL, INTENT(OUT) :: qcritcont ! critical specific humidity for contrail formation [kg/kg] |
---|
100 | REAL, INTENT(OUT) :: potcontfraP ! potential persistent contrail fraction [-] |
---|
101 | REAL, INTENT(OUT) :: potcontfraNP ! potential non-persistent contrail fraction [-] |
---|
102 | REAL, INTENT(OUT) :: contfra ! linear contrail fraction [-] |
---|
103 | REAL, INTENT(OUT) :: dcontfra_cir ! linear contrail fraction to cirrus cloud fraction tendency [s-1] |
---|
104 | REAL, INTENT(OUT) :: dcf_avi ! cloud fraction tendency because of aviation [s-1] |
---|
105 | REAL, INTENT(OUT) :: dqvc_avi ! specific ice content tendency because of aviation [kg/kg/s] |
---|
106 | REAL, INTENT(OUT) :: dqi_avi ! specific cloud water vapor tendency because of aviation [kg/kg/s] |
---|
107 | ! |
---|
108 | ! Local |
---|
109 | ! |
---|
110 | ! for Schmidt-Applemant criteria |
---|
111 | REAL, DIMENSION(1) :: ztemp, zpplay, qzero, zqsatl, zdqsatl |
---|
112 | REAL :: Gcont, qsatl_crit, psatl_crit, pcrit |
---|
113 | REAL :: pdf_x, pdf_y, pdf_e2, pdf_e3 |
---|
114 | REAL :: pdf_fra_above_qcritcont, pdf_fra_above_qsat, pdf_fra_above_qnuc |
---|
115 | REAL :: pdf_q_above_qcritcont, pdf_q_above_qsat, pdf_q_above_qnuc |
---|
116 | REAL :: qpotcontP |
---|
117 | ! |
---|
118 | ! for new contrail formation |
---|
119 | REAL :: contrail_cross_section, contfra_new |
---|
120 | |
---|
121 | qzero(:) = 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) |
---|
131 | Gcont = 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 |
---|
136 | Tcritcont = 226.69 + 9.43 * LOG( MAX(Gcont, 0.1) - 0.053 ) & |
---|
137 | + 0.72 * LOG( MAX(Gcont, 0.1) - 0.053 )**2 |
---|
138 | |
---|
139 | ztemp(1) = Tcritcont |
---|
140 | zpplay(1) = pplay |
---|
141 | CALL calc_qsat_ecmwf(1, ztemp, qzero, zpplay, RTT, 1, .FALSE., zqsatl, zdqsatl) |
---|
142 | qsatl_crit = zqsatl(1) |
---|
143 | |
---|
144 | psatl_crit = qsatl_crit / ( EPS_W + ( 1. - EPS_W ) * qsatl_crit ) * pplay |
---|
145 | pcrit = Gcont * ( temp - Tcritcont ) + psatl_crit |
---|
146 | qcritcont = EPS_W * pcrit / ( pplay - ( 1. - EPS_W ) * pcrit ) |
---|
147 | |
---|
148 | |
---|
149 | IF ( ( ( 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 | |
---|
181 | ELSE |
---|
182 | |
---|
183 | potcontfraNP = 0. |
---|
184 | potcontfraP = 0. |
---|
185 | |
---|
186 | ENDIF ! temp .LT. Tcritcont |
---|
187 | |
---|
188 | !--Convert existing contrail fraction into "natural" cirrus cloud fraction |
---|
189 | contfra = rcont_seri * cldfra |
---|
190 | dcontfra_cir = contfra * ( EXP( - dtime / linear_contrails_lifetime ) - 1. ) |
---|
191 | contfra = contfra + dcontfra_cir |
---|
192 | |
---|
193 | !--Add a source of contrails from aviation |
---|
194 | dcf_avi = 0. |
---|
195 | dqi_avi = 0. |
---|
196 | dqvc_avi = 0. |
---|
197 | IF ( 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 |
---|
210 | ENDIF |
---|
211 | |
---|
212 | END SUBROUTINE contrails_formation_evolution |
---|
213 | |
---|
214 | |
---|
215 | !********************************************************************************** |
---|
216 | SUBROUTINE 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 | |
---|
228 | USE lmdz_lscp_ini, ONLY: RPI, eps, ok_unadjusted_clouds |
---|
229 | USE lmdz_lscp_ini, ONLY: aspect_ratio_contrails, coef_mixing_contrails, & |
---|
230 | coef_shear_contrails, chi_mixing_contrails |
---|
231 | |
---|
232 | IMPLICIT NONE |
---|
233 | |
---|
234 | ! |
---|
235 | ! Input |
---|
236 | ! |
---|
237 | REAL, INTENT(IN) :: dtime ! time step [s] |
---|
238 | ! |
---|
239 | REAL, INTENT(IN) :: pplay ! layer pressure [Pa] |
---|
240 | REAL, INTENT(IN) :: shear ! vertical shear [s-1] |
---|
241 | REAL, INTENT(IN) :: pbl_eps ! TKE dissipation[m2/s3] |
---|
242 | REAL, INTENT(IN) :: cell_area ! cell area [m2] |
---|
243 | REAL, INTENT(IN) :: dz ! cell width [m] |
---|
244 | REAL, INTENT(IN) :: temp ! temperature [K] |
---|
245 | REAL, INTENT(IN) :: qtot ! total specific humidity (without precip) [kg/kg] |
---|
246 | REAL, INTENT(IN) :: qsat ! saturation specific humidity [kg/kg] |
---|
247 | REAL, INTENT(IN) :: subfra ! subsaturated clear sky fraction [-] |
---|
248 | REAL, INTENT(IN) :: qsub ! gridbox-mean subsaturated clear sky specific water [kg/kg] |
---|
249 | REAL, INTENT(IN) :: issrfra ! ISSR fraction [-] |
---|
250 | REAL, INTENT(IN) :: qissr ! gridbox-mean ISSR specific water [kg/kg] |
---|
251 | REAL, INTENT(IN) :: cldfra ! cloud fraction [-] |
---|
252 | REAL, INTENT(IN) :: qcld ! gridbox-mean cloudy specific total water [kg/kg] |
---|
253 | REAL, INTENT(IN) :: qvc ! gridbox-mean cloudy specific water vapor [kg/kg] |
---|
254 | REAL, INTENT(IN) :: rcont_seri ! ratio of contrails fraction to total cloud fraction [-] |
---|
255 | ! |
---|
256 | ! Input/Output |
---|
257 | ! |
---|
258 | REAL, INTENT(INOUT) :: dcf_mix ! cloud fraction tendency because of cloud mixing [s-1] |
---|
259 | REAL, INTENT(INOUT) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s] |
---|
260 | REAL, INTENT(INOUT) :: dqi_mix ! specific ice content tendency because of cloud mixing [kg/kg/s] |
---|
261 | REAL, INTENT(INOUT) :: dqt_mix ! total water tendency because of cloud mixing [kg/kg/s] |
---|
262 | REAL, INTENT(INOUT) :: dcf_mix_issr ! cloud fraction tendency because of cloud mixing in ISSR [s-1] |
---|
263 | REAL, INTENT(INOUT) :: dqt_mix_issr ! total water tendency because of cloud mixing in ISSR [kg/kg/s] |
---|
264 | ! |
---|
265 | ! Local |
---|
266 | ! |
---|
267 | REAL :: dqt_mix_sub_cont, dqt_mix_issr_cont |
---|
268 | REAL :: dcf_mix_sub_cont, dcf_mix_issr_cont |
---|
269 | REAL :: dqvc_mix_sub_cont, dqvc_mix_issr_cont |
---|
270 | REAL :: dcf_mix_cont, dqvc_mix_cont, dqi_mix_cont, dqt_mix_cont |
---|
271 | REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix |
---|
272 | REAL :: envfra_mix, cldfra_mix |
---|
273 | REAL :: L_shear, shear_fra |
---|
274 | REAL :: sigma_mix, issrfra_mix, subfra_mix |
---|
275 | REAL :: qvapincld, qvapinmix, qvapincld_new, qiceincld |
---|
276 | |
---|
277 | !--Initialisation |
---|
278 | dcf_mix_sub_cont = 0. |
---|
279 | dqt_mix_sub_cont = 0. |
---|
280 | dqvc_mix_sub_cont = 0. |
---|
281 | dcf_mix_issr_cont = 0. |
---|
282 | dqt_mix_issr_cont = 0. |
---|
283 | dqvc_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) |
---|
308 | bovera = 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) ) ) |
---|
311 | Povera = 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) ) |
---|
318 | a_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 ) |
---|
321 | N_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 |
---|
327 | L_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. |
---|
330 | L_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 |
---|
334 | envfra_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 |
---|
338 | cldfra_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 |
---|
348 | L_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) |
---|
354 | shear_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. |
---|
359 | envfra_mix = envfra_mix + shear_fra |
---|
360 | cldfra_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. |
---|
371 | sigma_mix = subfra / ( subfra + chi_mixing_contrails * issrfra ) |
---|
372 | subfra_mix = MIN( sigma_mix * envfra_mix, subfra ) |
---|
373 | issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra ) |
---|
374 | cldfra_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). |
---|
378 | IF ( 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 |
---|
417 | ENDIF ! 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 |
---|
422 | IF ( 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 | |
---|
442 | ENDIF ! issrfra .GT. eps |
---|
443 | |
---|
444 | !--Sum up the tendencies from subsaturated sky and supersaturated sky |
---|
445 | dcf_mix_cont = dcf_mix_sub_cont + dcf_mix_issr_cont |
---|
446 | dqt_mix_cont = dqt_mix_sub_cont + dqt_mix_issr_cont |
---|
447 | dqvc_mix_cont = dqvc_mix_sub_cont + dqvc_mix_issr_cont |
---|
448 | dqi_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 |
---|
453 | dcf_mix = ( 1. - rcont_seri ) * dcf_mix + rcont_seri * dcf_mix_cont |
---|
454 | dqvc_mix = ( 1. - rcont_seri ) * dqvc_mix + rcont_seri * dqvc_mix_cont |
---|
455 | dqi_mix = ( 1. - rcont_seri ) * dqi_mix + rcont_seri * dqi_mix_cont |
---|
456 | dqt_mix = ( 1. - rcont_seri ) * dqt_mix + rcont_seri * dqt_mix_cont |
---|
457 | dcf_mix_issr = ( 1. - rcont_seri ) * dcf_mix_issr + rcont_seri * dcf_mix_issr_cont |
---|
458 | dqt_mix_issr = ( 1. - rcont_seri ) * dqt_mix_issr + rcont_seri * dqt_mix_issr_cont |
---|
459 | |
---|
460 | END SUBROUTINE contrails_mixing |
---|
461 | |
---|
462 | |
---|
463 | !********************************************************************************** |
---|
464 | FUNCTION qvapincld_depsub_contrails( & |
---|
465 | qvapincld, qiceincld, temp, qsat, pplay, dtime) |
---|
466 | |
---|
467 | USE lmdz_lscp_ini, ONLY: RV, RLSTT, RTT, EPS_W |
---|
468 | USE lmdz_lscp_ini, ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice |
---|
469 | USE lmdz_lscp_ini, ONLY: rm_ice_crystals_contrails |
---|
470 | |
---|
471 | IMPLICIT NONE |
---|
472 | |
---|
473 | ! inputs |
---|
474 | REAL :: qvapincld ! |
---|
475 | REAL :: qiceincld ! |
---|
476 | REAL :: temp ! |
---|
477 | REAL :: qsat ! |
---|
478 | REAL :: pplay ! |
---|
479 | REAL :: dtime ! time step [s] |
---|
480 | ! output |
---|
481 | REAL :: qvapincld_depsub_contrails |
---|
482 | ! local |
---|
483 | REAL :: qice_ratio |
---|
484 | REAL :: pres_sat, rho, kappa |
---|
485 | REAL :: air_thermal_conduct, water_vapor_diff |
---|
486 | REAL :: 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 |
---|
550 | pres_sat = qsat / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay |
---|
551 | !--This formula for air thermal conductivity comes from Beard and Pruppacher (1971) |
---|
552 | air_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) |
---|
554 | water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4 |
---|
555 | kappa = 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 | |
---|
563 | rm_ice = rm_ice_crystals_contrails |
---|
564 | |
---|
565 | IF ( 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 ) |
---|
570 | ELSE |
---|
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. |
---|
592 | ENDIF ! qvapincld .GT. qsat |
---|
593 | |
---|
594 | END FUNCTION qvapincld_depsub_contrails |
---|
595 | |
---|
596 | |
---|
597 | !********************************************************************************** |
---|
598 | FUNCTION contrail_cross_section_onera(dz) |
---|
599 | |
---|
600 | USE lmdz_lscp_ini, ONLY: initial_width_contrails |
---|
601 | |
---|
602 | IMPLICIT NONE |
---|
603 | |
---|
604 | ! |
---|
605 | ! Input |
---|
606 | ! |
---|
607 | REAL :: dz ! cell width [m] |
---|
608 | ! |
---|
609 | ! Output |
---|
610 | ! |
---|
611 | REAL :: contrail_cross_section_onera ! [m2] |
---|
612 | ! |
---|
613 | ! Local |
---|
614 | ! |
---|
615 | |
---|
616 | contrail_cross_section_onera = initial_width_contrails * dz |
---|
617 | |
---|
618 | END FUNCTION contrail_cross_section_onera |
---|
619 | |
---|
620 | |
---|
621 | SUBROUTINE read_aviation_emissions( & |
---|
622 | klon, klev, latitude_deg, longitude_deg, pplay, & |
---|
623 | flight_dist, flight_h2o & |
---|
624 | ) |
---|
625 | |
---|
626 | IMPLICIT NONE |
---|
627 | ! |
---|
628 | ! Input |
---|
629 | ! |
---|
630 | INTEGER, INTENT(IN) :: klon, klev ! number of horizontal grid points and vertical levels |
---|
631 | REAL, INTENT(IN), DIMENSION(klon) :: latitude_deg ! latitude of the grid points [deg] |
---|
632 | REAL, INTENT(IN), DIMENSION(klon) :: longitude_deg ! longitude of the grid points [deg] |
---|
633 | REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! layer pressure [Pa] |
---|
634 | ! |
---|
635 | ! Output |
---|
636 | ! |
---|
637 | REAL, INTENT(OUT), DIMENSION(klon,klev) :: flight_dist ! aviation distance flown within the mesh [m/s/mesh] |
---|
638 | REAL, INTENT(OUT), DIMENSION(klon,klev) :: flight_h2o ! aviation H2O emitted within the mesh [kgH2O/s/mesh] |
---|
639 | ! |
---|
640 | ! Local |
---|
641 | ! |
---|
642 | INTEGER :: i, k |
---|
643 | |
---|
644 | !--Initialisation |
---|
645 | flight_dist(:,:) = 0. |
---|
646 | flight_h2o(:,:) = 0. |
---|
647 | |
---|
648 | DO 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 |
---|
653 | ENDDO |
---|
654 | |
---|
655 | END SUBROUTINE read_aviation_emissions |
---|
656 | |
---|
657 | END 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 |
---|