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

Last change on this file was 5790, checked in by aborella, 19 hours ago

Major modifs to treatment of contrails (from 2 classes to 2 moments) + diagnostics. Increased numerical efficiency
File size: 34.0 KB
Line 
1MODULE lmdz_aviation
2!----------------------------------------------------------------
3! Module for aviation and contrails
4
5IMPLICIT NONE
6
7! Arrays for the lecture of aviation files
8! The allocation is done in the read_aviation module
9! The size is (klon, nleva, 1) where
10! nleva            is the size of the vertical axis (read from file)
11! flight_dist_read is the number of km per second per m2
12! flight_fuel_read is the fuel consumed per second per m2
13! aviation_lev     is the value of the levels
14INTEGER, SAVE :: nleva  ! Size of the vertical axis in the file
15!$OMP THREADPRIVATE(nleva)
16REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE, PRIVATE :: flight_dist_read ! Aviation distance flown within the mesh [m/s/m2/vertmesh]
17!$OMP THREADPRIVATE(flight_dist_read)
18REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE, PRIVATE :: flight_fuel_read ! Aviation fuel consumed within the mesh [kg/s/m2/vertmesh]
19!$OMP THREADPRIVATE(flight_fuel_read)
20REAL, ALLOCATABLE, DIMENSION(:),     SAVE, PRIVATE :: aviation_lev     ! Pressure in the middle of the layers [Pa]
21!$OMP THREADPRIVATE(aviation_lev)
22
23CONTAINS
24
25SUBROUTINE aviation_water_emissions( &
26      klon, klev, dtime, pplay, temp, qtot, &
27      flight_fuel, d_q_avi &
28      )
29
30USE lmdz_lscp_ini, ONLY: RD, EI_H2O_aviation
31
32IMPLICIT NONE
33
34INTEGER,                      INTENT(IN)  :: klon, klev  ! number of horizontal grid points and vertical levels
35REAL,                         INTENT(IN)  :: dtime       ! time step [s]
36REAL, DIMENSION(klon,klev),   INTENT(IN)  :: pplay       ! mid-layer pressure [Pa]
37REAL, DIMENSION(klon,klev),   INTENT(IN)  :: temp        ! temperature (K)
38REAL, DIMENSION(klon,klev),   INTENT(IN)  :: qtot        ! total specific humidity (in vapor phase) [kg/kg]
39REAL, DIMENSION(klon,klev),   INTENT(IN)  :: flight_fuel ! aviation fuel consumption concentration [kg/s/m3]
40REAL, DIMENSION(klon,klev),   INTENT(OUT) :: d_q_avi     ! water vapor tendency from aviation [kg/kg]
41! Local
42INTEGER :: i, k
43REAL :: rho, flight_h2o
44
45DO i=1, klon
46  DO k=1, klev
47    !--Dry density [kg/m3]
48    rho = pplay(i,k) / temp(i,k) / RD
49    flight_h2o = flight_fuel(i,k) * EI_H2O_aviation
50
51    !--q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
52    ! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O ) 
53    !      = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) )
54    !--The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q)
55    !--flight_h2O is in kg H2O / s / m3
56   
57    !d_q_avi(i,k) = ( M_cell * qtot(i,k) + V_cell * flight_h2o * dtime * ( 1. - qtot(i,k) ) ) &
58    !             / ( M_cell             + V_cell * flight_h2o * dtime * ( 1. - qtot(i,k) ) ) &
59    !             - qtot(i,k)
60    !--NB., M_cell = V_cell * rho
61    !d_q_avi(i,k) = ( rho * qtot(i,k) + flight_h2o * dtime * ( 1. - qtot(i,k) ) ) &
62    !             / ( rho             + flight_h2o * dtime * ( 1. - qtot(i,k) ) ) &
63    !             - qtot(i,k)
64    !--Same formula, more computationally effective but less readable
65    d_q_avi(i,k) = flight_h2o * ( 1. - qtot(i,k) ) &
66                 / ( rho / dtime / ( 1. - qtot(i,k) ) + flight_h2o )
67  ENDDO
68ENDDO
69
70END SUBROUTINE aviation_water_emissions
71
72
73!**********************************************************************************
74SUBROUTINE contrails_formation( &
75      klon, dtime, pplay, temp, rho, qsat, qsatl, gamma_cond, flight_dist, flight_fuel, &
76      clrfra, qclr, pdf_scale, pdf_alpha, pdf_gamma, keepgoing, pt_pron_clds, &
77      Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
78      AEI_contrails, AEI_surv_contrails, fsurv_contrails, section_contrails, &
79      dcfc_ini, dqic_ini, dqtc_ini, dnic_ini)
80
81USE lmdz_lscp_ini, ONLY: RPI, RCPD, EPS_W, RTT
82USE lmdz_lscp_ini, ONLY: EI_H2O_aviation, qheat_fuel_aviation, prop_efficiency_aviation
83USE lmdz_lscp_ini, ONLY: eps, temp_nowater
84USE lmdz_lscp_ini, ONLY: Nice_init_contrails
85
86USE lmdz_lscp_tools, ONLY: GAMMAINC, calc_qsat_ecmwf
87
88IMPLICIT NONE
89
90!
91! Input
92!
93INTEGER,  INTENT(IN)                   :: klon         ! number of horizontal grid points
94REAL,     INTENT(IN)                   :: dtime        ! time step [s]
95REAL,     INTENT(IN) , DIMENSION(klon) :: pplay        ! layer pressure [Pa]
96REAL,     INTENT(IN) , DIMENSION(klon) :: temp         ! temperature [K]
97REAL,     INTENT(IN) , DIMENSION(klon) :: rho          ! air density [kg/m3]
98REAL,     INTENT(IN) , DIMENSION(klon) :: qsat         ! saturation specific humidity [kg/kg]
99REAL,     INTENT(IN) , DIMENSION(klon) :: qsatl        ! saturation specific humidity w.r.t. liquid [kg/kg]
100REAL,     INTENT(IN) , DIMENSION(klon) :: gamma_cond   ! condensation threshold w.r.t. qsat [-]
101REAL,     INTENT(IN) , DIMENSION(klon) :: flight_dist  ! aviation distance flown concentration [m/s/m3]
102REAL,     INTENT(IN) , DIMENSION(klon) :: flight_fuel  ! aviation fuel consumption concentration [kg/s/m3]
103REAL,     INTENT(IN) , DIMENSION(klon) :: clrfra       ! clear sky fraction [-]
104REAL,     INTENT(IN) , DIMENSION(klon) :: qclr         ! clear sky specific humidity [kg/kg]
105REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_scale    ! scale parameter of the clear sky PDF [%]
106REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_alpha    ! shape parameter of the clear sky PDF [-]
107REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_gamma    ! tmp parameter of the clear sky PDF [-]
108LOGICAL,  INTENT(IN) , DIMENSION(klon) :: keepgoing    ! .TRUE. if a new condensation loop should be computed
109LOGICAL,  INTENT(IN) , DIMENSION(klon) :: pt_pron_clds  ! .TRUE. if clouds are prognostic in this mesh
110!
111! Output
112!
113REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcritcont    ! critical temperature for contrail formation [K]
114REAL,     INTENT(INOUT), DIMENSION(klon) :: qcritcont    ! critical specific humidity for contrail formation [kg/kg]
115REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraP  ! potential persistent contrail fraction [-]
116REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraNP ! potential non-persistent contrail fraction [-]
117REAL,     INTENT(INOUT), DIMENSION(klon) :: AEI_contrails ! Apparent emission index contrails [#/kg]
118REAL,     INTENT(INOUT), DIMENSION(klon) :: AEI_surv_contrails ! Apparent emission index contrails after vortex loss [#/kg]
119REAL,     INTENT(INOUT), DIMENSION(klon) :: fsurv_contrails ! Survival fraction after vortex loss [-]
120REAL,     INTENT(INOUT), DIMENSION(klon) :: section_contrails ! Cross section of newly formed contrails [m2]
121REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfc_ini     ! contrails cloud fraction tendency bec ause of initial formation [s-1]
122REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_ini     ! contrails ice specific humidity tende ncy because of initial formation [kg/kg/s]
123REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtc_ini     ! contrails total specific humidity ten dency because of initial formation [kg/kg/s]
124REAL,     INTENT(INOUT), DIMENSION(klon) :: dnic_ini     ! contrails ice crystals concentration ten dency because of initial formation [#/kg/s]
125!
126! Local
127!
128INTEGER :: i
129! for Schmidt-Appleman criteria
130REAL, DIMENSION(klon) :: qzero
131REAL, DIMENSION(klon) :: qsatl_crit, dqsatl_crit
132REAL :: Gcont, psatl_crit, pcrit
133REAL :: rhl_clr, pdf_loc
134REAL :: pdf_x, pdf_y, pdf_e3
135REAL :: pdf_fra_above_qcritcont, pdf_fra_above_qsat
136REAL :: pdf_q_above_qcritcont, pdf_q_above_qsat
137REAL :: qpotcontP
138!
139! for new contrail formation
140REAL :: dist_contrails, Nice_per_m_init_contrails
141REAL :: icesat_ratio
142
143qzero(:) = 0.
144
145DO i = 1, klon
146  IF ( keepgoing(i) ) THEN
147    qcritcont(i)    = 0.
148    potcontfraP(i)  = 0.
149    potcontfraNP(i) = 0.
150  ENDIF
151ENDDO
152
153!---------------------------------
154!--  SCHMIDT-APPLEMAN CRITERIA  --
155!---------------------------------
156!--Revised by Schumann (1996) and Rap et al. (2010)
157
158DO i = 1, klon
159  !--Gcont is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
160  !--in Pa / K. See Rap et al. (2010) in JGR.
161  !--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1)
162  Gcont = EI_H2O_aviation * RCPD * pplay(i) &
163         / ( EPS_W * qheat_fuel_aviation * ( 1. - prop_efficiency_aviation ) )
164  !--critical temperature below which no liquid contrail can form in exhaust
165  !--noted T_LM in Schumann (1996), their Eq. 31 in Appendix 2
166  !--in Kelvins
167  Tcritcont(i) = 226.69 + 9.43 * LOG( MAX(Gcont, 0.1) - 0.053 ) &
168         + 0.72 * LOG( MAX(Gcont, 0.1) - 0.053 )**2
169ENDDO
170
171CALL calc_qsat_ecmwf(klon, Tcritcont, qzero, pplay, RTT, 1, .FALSE., qsatl_crit, dqsatl_crit)
172
173DO i = 1, klon
174  IF ( keepgoing(i) .AND. pt_pron_clds(i) ) THEN
175   
176    psatl_crit = qsatl_crit(i) / ( EPS_W + ( 1. - EPS_W ) * qsatl_crit(i) ) * pplay(i)
177    pcrit = Gcont * ( temp(i) - Tcritcont(i) ) + psatl_crit
178    qcritcont(i) = EPS_W * pcrit / ( pplay(i) - ( 1. - EPS_W ) * pcrit )
179
180
181    IF ( ( clrfra(i) .GT. eps ) .AND. ( temp(i) .LT. Tcritcont(i) ) ) THEN
182   
183      rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
184      pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
185      pdf_x = qcritcont(i) / qsatl(i) * 100.
186      pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
187      IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
188        pdf_fra_above_qcritcont = 0.
189        pdf_q_above_qcritcont = 0.
190      ELSEIF ( pdf_y .LT. -10. ) THEN
191        pdf_fra_above_qcritcont = clrfra(i)
192        pdf_q_above_qcritcont = qclr(i)
193      ELSE
194        pdf_y = EXP( pdf_y )
195        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
196        pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
197        pdf_fra_above_qcritcont = EXP( - pdf_y ) * clrfra(i)
198        pdf_q_above_qcritcont = ( pdf_e3 * clrfra(i) &
199                              + pdf_loc * pdf_fra_above_qcritcont ) * qsatl(i) / 100.
200      ENDIF
201   
202      pdf_x = qsat(i) / qsatl(i) * 100.
203      pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
204      IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
205        pdf_fra_above_qsat = 0.
206        pdf_q_above_qsat = 0.
207      ELSEIF ( pdf_y .LT. -10. ) THEN
208        pdf_fra_above_qsat = clrfra(i)
209        pdf_q_above_qsat = qclr(i)
210      ELSE
211        pdf_y = EXP( pdf_y )
212        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
213        pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
214        pdf_fra_above_qsat = EXP( - pdf_y ) * clrfra(i)
215        pdf_q_above_qsat = ( pdf_e3 * clrfra(i) &
216                         + pdf_loc * pdf_fra_above_qsat ) * qsatl(i) / 100.
217      ENDIF
218   
219      potcontfraNP(i) = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat)
220      potcontfraP(i) = MIN(pdf_fra_above_qsat, pdf_fra_above_qcritcont)
221      qpotcontP = MIN(pdf_q_above_qsat, pdf_q_above_qcritcont)
222
223    ENDIF ! temp .LT. Tcritcont
224   
225    !--Add a source of contrails from aviation
226    IF ( ( potcontfraP(i) .GT. eps ) .AND. ( flight_dist(i) .GT. 1e-20 ) ) THEN
227      !section_contrails(i) = CONTRAIL_CROSS_SECTION_ONERA()
228      section_contrails(i) = CONTRAIL_CROSS_SECTION_SCHUMANN( &
229              dtime, rho(i), flight_dist(i), flight_fuel(i))
230      icesat_ratio = qpotcontP / potcontfraP(i) / qsat(i)
231      !--If Nice init is fixed in the plume
232      !Nice_per_m_init_contrails = Nice_init_contrails * 1e6 * section_contrails(i)
233      !--Else if it is parameterized
234      CALL CONTRAIL_ICE_NUMBER_CONCENTRATION(temp(i), icesat_ratio, rho(i), &
235              flight_dist(i), flight_fuel(i), Nice_per_m_init_contrails, &
236              AEI_contrails(i), AEI_surv_contrails(i), fsurv_contrails(i))
237
238      !--If Nice_per_m_init_contrails .EQ. 0., all the crystals were lost in the vortex
239      IF ( Nice_per_m_init_contrails .GT. 0. ) THEN
240        !--dist_contrails is in meters of contrails per m3 (gridbox-average)
241        dist_contrails = flight_dist(i) * dtime * potcontfraP(i)
242        dcfc_ini(i) = dist_contrails * section_contrails(i)
243        dqtc_ini(i) = icesat_ratio * qsat(i) * dcfc_ini(i)
244        dqic_ini(i) = dqtc_ini(i) - qsat(i) * dcfc_ini(i)
245        dnic_ini(i) = Nice_per_m_init_contrails * dist_contrails / rho(i)
246      ENDIF
247    ENDIF
248
249  ENDIF ! keepgoing
250ENDDO
251
252END SUBROUTINE contrails_formation
253
254
255!**********************************************************************************
256FUNCTION CONTRAIL_CROSS_SECTION_ONERA()
257
258USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails
259
260IMPLICIT NONE
261
262!
263! Output
264!
265REAL :: contrail_cross_section_onera ! [m2]
266!
267! Local
268!
269
270contrail_cross_section_onera = initial_width_contrails * initial_height_contrails
271
272END FUNCTION CONTRAIL_CROSS_SECTION_ONERA
273
274
275!**********************************************************************************
276!--Based on Schumann (1998)
277FUNCTION CONTRAIL_CROSS_SECTION_SCHUMANN(dtime, rho_air, flight_dist, flight_fuel)
278
279USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails
280
281IMPLICIT NONE
282
283!
284! Input
285!
286REAL :: dtime       ! timestep [s]
287REAL :: rho_air     ! air density [kg/m3]
288REAL :: flight_dist ! flown distance [m/s/m3]
289REAL :: flight_fuel ! fuel consumed [kg/s/m3]
290!
291! Output
292!
293REAL :: contrail_cross_section_schumann ! [m2]
294!
295! Local
296!
297REAL :: dilution_factor
298
299!--The contrail is formed on average in the middle of the timestep
300dilution_factor = 7000. * (dtime / 2.)**0.8
301contrail_cross_section_schumann = flight_fuel / flight_dist * dilution_factor / rho_air
302
303END FUNCTION CONTRAIL_CROSS_SECTION_SCHUMANN
304
305
306!**********************************************************************************
307SUBROUTINE CONTRAIL_ICE_NUMBER_CONCENTRATION(temp, icesat_ratio, rho_air, &
308        flight_dist, flight_fuel, Nice_per_m_init_contrails, &
309        AEI_contrails, AEI_surv_contrails, fsurv_contrails)
310
311USE lmdz_lscp_ini, ONLY: RPI
312USE lmdz_lscp_ini, ONLY: EI_soot_aviation, air_to_fuel_ratio_engine, wingspan
313USE lmdz_lscp_ini, ONLY: Naer_amb, raer_amb_mean, raer_amb_std, r_soot_mean, r_soot_std
314USE lmdz_lscp_ini, ONLY: circ_0_loss, plume_area_loss, nice_init_ref_loss
315USE lmdz_lscp_ini, ONLY: N_Brunt_Vaisala_aviation, EI_H2O_aviation
316
317IMPLICIT NONE
318
319!
320! Input
321!
322REAL, INTENT(IN)  :: temp         ! air temperature [K]
323REAL, INTENT(IN)  :: icesat_ratio ! ice saturation ratio [-]
324REAL, INTENT(IN)  :: rho_air      ! air density [kg/m3]
325REAL, INTENT(IN)  :: flight_dist  ! flown distance [m/s/m3]
326REAL, INTENT(IN)  :: flight_fuel  ! fuel consumed [kg/s/m3]
327!
328! Output
329!
330REAL, INTENT(OUT) :: Nice_per_m_init_contrails ! [#/m]
331REAL, INTENT(OUT) :: AEI_contrails             ! [#/kg]
332REAL, INTENT(OUT) :: AEI_surv_contrails        ! [#/kg]
333REAL, INTENT(OUT) :: fsurv_contrails           ! [-]
334!
335! Local
336!
337! For initial ice nucleation
338REAL :: log_liqsat_ratio, coef_expo, dil_coef
339REAL :: rd_act_amb, rd_act_soot, phi_amb, phi_soot
340REAL :: AEI_soot, AEI_amb
341! For ice crystals loss
342REAL :: rho_emit, temp_205, nice_init
343REAL :: z_atm, z_emit, z_desc, z_delta
344REAL :: Nice_per_m, fuel_per_m, frac_surv
345
346!------------------------------
347!--  INITIAL ICE NUCLEATION  --
348!------------------------------
349!
350!--Karcher et al. (2015), Bier and Burkhardt (2019, 2022)
351!log_liqsat_ratio = LOG(liq_satratio)
352
353!! dry core radius in nm
354!! HERE SHOULD IT BE liqsat_ratio OR liqsat_ratio - 1. ?
355!rd_act_amb = (4. / 27. / LOG(liqsat_ratio)**2 / 0.5)**(1./3.)
356!! Integrate lognormal distribution between rd_act_amb and +inf
357!coef_expo = 4. / SQRT(2. * RPI) / LOG(raer_amb_std)
358!phi_amb = 1. / (1. + (rd_act_amb / raer_amb_mean)**coef_expo)
359!
360!! BU22, Eq. A1, dry core radius in nm
361!rd_act_soot = 0.96453426 + 1.21152973 / log_liqsat_ratio - 0.00520358 / log_liqsat_ratio**2 &
362!        + 2.32286432e-5 / log_liqsat_ratio**3 - 4.36979055e-8 / log_liqsat_ratio**4
363!rd_act_soot = MIN(150., MAX(1., rd_act_soot))
364!! Integrate lognormal distribution between rd_act_amb and +inf
365!coef_expo = 4. / SQRT(2. * RPI) / LOG(r_soot_std)
366!phi_amb = 1. / (1. + (rd_act_soot / r_soot_mean)**coef_expo)
367!
368!dil_coef = (0.01 / t0)**0.9
369!
370!! BU22, Eq. 5b
371!AEI_soot = phi_soot * EI_soot_aviation
372!AEI_amb = phi_amb * air_to_fuel_ratio_engine * (1. - dil_coef) / dil_coef &
373!        / rho_air * Naer_amb * 1e6
374!AEI_contrails = AEI_soot + AEI_amb
375AEI_contrails = EI_soot_aviation
376
377
378!----------------------------
379!--    ICE CRYSTALS LOSS   --
380!----------------------------
381!
382!--Based on Lottermoser and Unterstrasser (2025, LU25 hereinafter)
383!--which is an update of Unterstrasser (2016, U16 hereinafter)
384
385! fuel consumption in kg/m flown
386fuel_per_m = flight_fuel / flight_dist
387
388! LU25, Eq. A2
389z_atm = 607.46 * (icesat_ratio - 1.)**0.897 * (temp / 205.)**2.225
390
391! water vapor emissions [kg/m flown], LU25, Eq. 2
392! U16, Eq. 6
393rho_emit = fuel_per_m * EI_H2O_aviation / plume_area_loss
394! LU25, Eq. A3
395temp_205 = temp - 205.
396z_emit = 1106.6 * (rho_emit * 1e5)**(0.678 + 0.0116 * temp_205) &
397        * EXP(- (0.0807 + 0.000428 * temp_205) * temp_205)
398
399! U16, Eq. 4
400z_desc = SQRT(8. * circ_0_loss / RPI / N_Brunt_Vaisala_aviation )
401
402! initial number of ice crystals [#/m flown], LU25, Eq. 1
403Nice_per_m = fuel_per_m * AEI_contrails
404! ice crystals number concentration [#/m3], LU25, Eq. A1
405nice_init = Nice_per_m / plume_area_loss
406! LU25, Eq. 9, 13d, 13e, 13f and 13g
407z_delta = (nice_init / nice_init_ref_loss)**(-0.16) * (1.27 * z_atm + 0.42 * z_emit) - 0.49 * z_desc
408
409! LU25, Eq. 11, 12, 13a, 13b and 13c
410fsurv_contrails = 0.42 + 1.31 / RPI * ATAN(-1. + z_delta / 100.)
411fsurv_contrails = MIN(1., MAX(0., fsurv_contrails))
412
413Nice_per_m_init_contrails = Nice_per_m * fsurv_contrails
414AEI_surv_contrails = AEI_contrails * fsurv_contrails
415
416END SUBROUTINE CONTRAIL_ICE_NUMBER_CONCENTRATION
417
418
419SUBROUTINE init_read_aviation_emissions
420
421  USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured
422  USE mod_phys_lmdz_para, ONLY: is_omp_master
423  USE lmdz_xios, ONLY: xios_set_file_attr, xios_set_field_attr
424
425  IF (grid_type==unstructured) THEN
426    IF (is_omp_master) THEN
427      ! Activate aviation file
428      CALL xios_set_file_attr("aviation_file", enabled=.TRUE.)
429      ! Activate aviation fields
430      CALL xios_set_field_attr("DISTFLOWN_read", enabled=.TRUE.)
431      CALL xios_set_field_attr("DISTFLOWN_interp", enabled=.TRUE.)
432      CALL xios_set_field_attr("FUELBURN_read", enabled=.TRUE.)
433      CALL xios_set_field_attr("FUELBURN_interp", enabled=.TRUE.)
434    ENDIF
435  ENDIF
436
437END SUBROUTINE init_read_aviation_emissions
438
439
440SUBROUTINE read_aviation_emissions(klon, klev)
441    ! This subroutine allows to read the traffic density data read in the file aviation.nc
442    ! This file is defined in ./COMP/lmdz.card
443    ! Field names in context_input_lmdz.xml should be the same as in the file.
444    USE netcdf
445    USE mod_phys_lmdz_para
446    USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, nbp_lev_=>nbp_lev, &
447                                 klon_glo, grid2Dto1D_glo, grid_type, unstructured
448    USE iophy, ONLY : io_lon, io_lat
449    USE lmdz_xios
450    USE print_control_mod, ONLY: lunout
451    USE readaerosol_mod, ONLY: check_err
452    USE lmdz_lscp_ini, ONLY: EI_H2O_aviation
453    IMPLICIT NONE
454
455    INTEGER, INTENT(IN) :: klon, klev  ! number of horizontal grid points and vertical levels
456
457    !----------------------------------------------------
458    ! Local variable
459    !----------------------------------------------------
460    REAL, ALLOCATABLE :: flight_dist_mpi(:,:,:), flight_fuel_mpi(:,:,:)
461    INTEGER :: ierr
462    LOGICAL, SAVE :: first = .TRUE.
463    !$OMP THREADPRIVATE(first)
464
465    ! FOR REGULAR LON LAT
466    CHARACTER(len=30)     :: fname
467    INTEGER               :: nid, dimid, varid
468    INTEGER               :: imth, i, j, k
469    REAL                  :: npole, spole
470    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_dist_glo2D
471    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_fuel_glo2D
472    REAL, ALLOCATABLE, DIMENSION(:)       :: vartmp_lev
473    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: vartmp
474    REAL, ALLOCATABLE, DIMENSION(:)       :: lon_src              ! longitudes in file
475    REAL, ALLOCATABLE, DIMENSION(:)       :: lat_src, lat_src_inv ! latitudes in file
476    LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
477    INTEGER                               :: nbp_lon, nbp_lat
478
479
480IF (grid_type==unstructured) THEN
481
482    IF (first) THEN
483      IF (is_omp_master) THEN
484        ! Get number of vertical levels and level values
485        CALL xios_get_axis_attr("aviation_lev", n_glo=nleva)
486      ENDIF
487      CALL bcast_omp(nleva)
488
489      ! Allocation of arrays
490      ALLOCATE(flight_dist_read(klon, nleva, 1), STAT=ierr)
491      IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_dist',1)
492      ALLOCATE(flight_fuel_read(klon, nleva, 1), STAT=ierr)
493      IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_fuel',1)
494      ALLOCATE(aviation_lev(nleva), STAT=ierr)
495      IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate aviation_lev',1)
496
497      IF (is_omp_master) THEN
498        ! Get number of vertical levels and level values
499        CALL xios_get_axis_attr("aviation_lev", value=aviation_lev(:))
500      ENDIF
501      CALL bcast_omp(aviation_lev)
502      first = .FALSE.
503    ENDIF ! first
504
505    ! Read the data from the file
506    ! is_omp_master is necessary to make XIOS works
507    IF (is_omp_master) THEN
508        ALLOCATE(flight_dist_mpi(klon_mpi, nleva,1), STAT=ierr)
509        IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_dist_mpi',1)
510        ALLOCATE(flight_fuel_mpi(klon_mpi, nleva,1), STAT=ierr)
511        IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_fuel_mpi',1)
512        CALL xios_recv_field("DISTFLOWN_interp", flight_dist_mpi(:,:,1))
513        CALL xios_recv_field("FUELBURN_interp", flight_fuel_mpi(:,:,1))
514    ENDIF
515
516    ! Propagate to other OMP threads: flight_dist_mpi(klon_mpi,klev) to flight_dist(klon,klev)
517    CALL scatter_omp(flight_dist_mpi, flight_dist_read)
518    CALL scatter_omp(flight_fuel_mpi, flight_fuel_read)
519
520ELSE
521
522    nbp_lon=nbp_lon_
523    nbp_lat=nbp_lat_
524   
525    IF (is_mpi_root) THEN
526      ALLOCATE(lon_src(nbp_lon))
527      ALLOCATE(lat_src(nbp_lat))
528      ALLOCATE(lat_src_inv(nbp_lat))
529    ELSE
530      ALLOCATE(flight_dist_glo2D(0,0,0,0))
531      ALLOCATE(flight_fuel_glo2D(0,0,0,0))
532    ENDIF
533
534    IF (is_mpi_root .AND. is_omp_root) THEN
535
536! 1) Open file
537!****************************************************************************************
538      fname = 'aviation.nc'
539      CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, nid), "pb open "//TRIM(fname) )
540
541
542! Test for equal longitudes and latitudes in file and model
543!****************************************************************************************
544      ! Read and test longitudes
545      CALL check_err( nf90_inq_varid(nid, 'longitude', varid), "pb inq lon" )
546      CALL check_err( nf90_get_var(nid, varid, lon_src(:)), "pb get lon" )
547     
548      IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
549         WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
550         WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
551         WRITE(lunout,*) 'longitudes in model :', io_lon
552       
553         CALL abort_physic('lmdz_aviation', 'longitudes are not the same in file and model',1)
554      END IF
555
556      ! Read and test latitudes
557      CALL check_err( nf90_inq_varid(nid, 'latitude', varid),"pb inq lat" )
558      CALL check_err( nf90_get_var(nid, varid, lat_src(:)),"pb get lat" )
559
560      ! Invert source latitudes
561      DO j = 1, nbp_lat
562         lat_src_inv(j) = lat_src(nbp_lat+1-j)
563      END DO
564
565      IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
566         ! Latitudes are the same
567         invert_lat=.FALSE.
568      ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
569         ! Inverted source latitudes correspond to model latitudes
570         WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
571         invert_lat=.TRUE.
572      ELSE
573         WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
574         WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src     
575         WRITE(lunout,*) 'latitudes in model :', io_lat
576         CALL abort_physic('lmdz_aviation', 'latitudes do not correspond between file and model',1)
577      END IF
578
579       
580! 2) Find vertical dimension nleva
581!****************************************************************************************
582       ierr = nf90_inq_dimid(nid, 'pressure_Pa', dimid)
583       CALL check_err( nf90_inquire_dimension(nid, dimid, len = nleva), "pb inq dim for PRESNIVS or lev" )
584       
585     ! Allocate variables depending on the number of vertical levels
586       ALLOCATE(flight_dist_glo2D(nbp_lon, nbp_lat, nleva, 1), flight_fuel_glo2D(nbp_lon, nbp_lat, nleva, 1), stat=ierr)
587       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 1',1)
588
589       ALLOCATE(aviation_lev(nleva), stat=ierr)
590       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 2',1)
591
592! 3) Read all variables from file
593!****************************************************************************************
594
595       ! Get variable id
596       CALL check_err( nf90_inq_varid(nid, 'seg_length', varid),"pb inq var seg_length_m" )
597       ! Get the variable
598       CALL check_err( nf90_get_var(nid, varid, flight_dist_glo2D),"pb get var seg_length_m"  )
599
600       ! Get variable id
601       CALL check_err( nf90_inq_varid(nid, 'fuel_burn', varid),"pb inq var fuel_burn" )
602       ! Get the variable
603       CALL check_err( nf90_get_var(nid, varid, flight_fuel_glo2D),"pb get var fuel_burn"  )
604
605       ! Get variable id
606       CALL check_err( nf90_inq_varid(nid, "pressure_Pa", varid),"pb inq var pressure_Pa" )
607       ! Get the variable
608       CALL check_err( nf90_get_var(nid, varid, aviation_lev),"pb get var pressure_Pa" )
609         
610
611! 4) Close file 
612!****************************************************************************************
613       CALL check_err( nf90_close(nid), "pb in close" )
614     
615
616! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
617!****************************************************************************************
618! Test if vertical levels have to be inversed
619
620       IF (aviation_lev(1) < aviation_lev(nleva)) THEN
621         
622          ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva), vartmp_lev(nleva))
623
624          ! Inverse vertical levels for flight_dist_glo2D
625          vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1)
626          DO k=1, nleva
627             DO j=1, nbp_lat
628                DO i=1, nbp_lon
629                   flight_dist_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
630                END DO
631             END DO
632          END DO
633
634          ! Inverse vertical levels for flight_fuel_glo2D
635          vartmp(:,:,:) = flight_fuel_glo2D(:,:,:,1)
636          DO k=1, nleva
637             DO j=1, nbp_lat
638                DO i=1, nbp_lon
639                   flight_fuel_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
640                END DO
641             END DO
642          END DO
643           
644          ALLOCATE(vartmp_lev(nleva))
645          ! Inverte vertical axes for aviation_lev
646          vartmp_lev(:) = aviation_lev(:)
647          DO k=1, nleva
648             aviation_lev(k) = vartmp_lev(nleva+1-k)
649          END DO
650
651          DEALLOCATE(vartmp, vartmp_lev)
652
653       ELSE
654          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
655       END IF
656
657
658!     - Invert latitudes if necessary
659       IF (invert_lat) THEN
660
661          ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva))
662
663          ! Invert latitudes for the variable
664          vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1) ! use varmth temporarly
665          DO k=1,nleva
666             DO j=1,nbp_lat
667                DO i=1,nbp_lon
668                   flight_dist_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
669                END DO
670             END DO
671          END DO
672
673          ! Invert latitudes for the variable
674          vartmp(:,:,:) = flight_fuel_glo2D(:,:,:,1) ! use varmth temporarly
675          DO k=1,nleva
676             DO j=1,nbp_lat
677                DO i=1,nbp_lon
678                   flight_fuel_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
679                END DO
680             END DO
681          END DO
682
683          DEALLOCATE(vartmp)
684       END IF ! invert_lat
685       
686       ! Do zonal mead at poles and distribut at whole first and last latitude
687       DO k=1, nleva
688          npole=0.  ! North pole, j=1
689          spole=0.  ! South pole, j=nbp_lat       
690          DO i=1,nbp_lon
691             npole = npole + flight_dist_glo2D(i,1,k,1)
692             spole = spole + flight_dist_glo2D(i,nbp_lat,k,1)
693          END DO
694          npole = npole/REAL(nbp_lon)
695          spole = spole/REAL(nbp_lon)
696          flight_dist_glo2D(:,1,    k,1) = npole
697          flight_dist_glo2D(:,nbp_lat,k,1) = spole
698       END DO
699
700       ! Do zonal mead at poles and distribut at whole first and last latitude
701       DO k=1, nleva
702          npole=0.  ! North pole, j=1
703          spole=0.  ! South pole, j=nbp_lat       
704          DO i=1,nbp_lon
705             npole = npole + flight_fuel_glo2D(i,1,k,1)
706             spole = spole + flight_fuel_glo2D(i,nbp_lat,k,1)
707          END DO
708          npole = npole/REAL(nbp_lon)
709          spole = spole/REAL(nbp_lon)
710          flight_fuel_glo2D(:,1,    k,1) = npole
711          flight_fuel_glo2D(:,nbp_lat,k,1) = spole
712       END DO
713     
714       ALLOCATE(flight_dist_mpi(klon_glo, nleva, 1), flight_fuel_mpi(klon_glo, nleva, 1), stat=ierr)
715       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 3',1)
716     
717       ! Transform from 2D to 1D field
718       CALL grid2Dto1D_glo(flight_dist_glo2D, flight_dist_mpi)
719       CALL grid2Dto1D_glo(flight_fuel_glo2D, flight_fuel_mpi)
720   
721    ELSE
722        ALLOCATE(flight_dist_mpi(0,0,0), flight_fuel_mpi(0,0,0))
723    END IF ! is_mpi_root .AND. is_omp_root
724
725!$OMP BARRIER
726 
727! 6) Distribute to all processes
728!    Scatter global field(klon_glo) to local process domain(klon)
729!    and distribute nleva to all processes
730!****************************************************************************************
731
732    ! Distribute nleva
733    CALL bcast(nleva)
734
735    ! Allocate and distribute pt_ap and pt_b
736    IF (.NOT. ALLOCATED(aviation_lev)) THEN  ! if pt_ap is allocated also pt_b is allocated
737       ALLOCATE(aviation_lev(nleva), stat=ierr)
738       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 4',1)
739    END IF
740    CALL bcast(aviation_lev)
741
742    ! Allocate space for output pointer variable at local process
743    ALLOCATE(flight_dist_read(klon, nleva, 1), flight_fuel_read(klon, nleva, 1), stat=ierr)
744    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 5',1)
745
746    ! Scatter global field to local domain at local process
747    CALL scatter(flight_dist_mpi, flight_dist_read)
748    CALL scatter(flight_fuel_mpi, flight_fuel_read)
749
750ENDIF
751
752END SUBROUTINE read_aviation_emissions
753
754SUBROUTINE vertical_interpolation_aviation(klon, klev, paprs, pplay, temp, flight_dist, flight_fuel)
755    ! This subroutine performs the vertical interpolation from the read data in aviation.nc
756    ! where there are nleva vertical levels described in aviation_lev to the klev levels or
757    ! the model.
758    ! flight_dist_read(klon,nleva) -> flight_dist(klon, klev)
759    ! flight_fuel_read(klon,nleva) -> flight_fuel(klon, klev)
760
761    USE lmdz_lscp_ini, ONLY: RD, RG
762    USE lmdz_lscp_ini, ONLY: aviation_coef
763
764    IMPLICIT NONE
765
766    INTEGER,                    INTENT(IN)  :: klon, klev  ! number of horizontal grid points and vertical levels
767    REAL, INTENT(IN)    :: paprs(klon, klev+1) ! inter-layer pressure [Pa]
768    REAL, INTENT(IN)    :: pplay(klon, klev) ! mid-layer pressure [Pa]
769    REAL, INTENT(IN)    :: temp(klon, klev) ! temperature [K]
770    REAL, INTENT(OUT)   :: flight_dist(klon,klev) ! Aviation distance flown concentration [m/s/m3]
771    REAL, INTENT(OUT)   :: flight_fuel(klon,klev) ! Aviation fuel burn concentration [kg/s/m3]
772
773    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
774    !  Local variable
775    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
776    REAL :: aviation_interface(1:nleva+1) ! Pressure of aviation file interfaces [ Pa ]
777    INTEGER :: k, kori  ! Loop index for vertical layers
778    INTEGER :: i  ! Loop index for horizontal grid
779    REAL :: zfrac ! Fraction of layer kori in layer k
780    REAL :: width_read_layer(1:nleva) ! width of a given layer [ Pa ]
781    REAL :: rho, rhodz, dz
782
783    ! Initialisation
784    flight_dist(:,:) = 0.
785    flight_fuel(:,:) = 0.
786
787    ! Compute the array with the vertical interface
788    ! It starts at 1 and has length nleva + 1
789    ! Note that aviation_lev has nleva and gives the altitude in the middle of the layers
790    ! Surface pressure in standard atmosphere model [ Pa ]
791    aviation_interface(1) = 101325.
792    DO kori=2, nleva
793        aviation_interface(kori) = (aviation_lev(kori-1)+aviation_lev(kori))/2.0  ! [ Pa ]
794    ENDDO
795    ! Last interface - we assume the same spacing as the very last one
796    aviation_interface(nleva+1) = aviation_interface(nleva) - (aviation_lev(nleva-1) - aviation_lev(nleva))
797
798    ! Vertical width of each layer of the read file
799    ! It is positive
800    DO kori=1, nleva
801        width_read_layer(kori) = aviation_interface(kori) - aviation_interface(kori+1)
802    ENDDO
803
804    ! Vertical reprojection
805    ! The loop over klon is induced since it is done by MPI threads
806    ! zfrac is the fraction of layer kori (read file) included in layer k (model)
807    DO i=1,klon
808        DO k=1, klev
809            DO kori=1,nleva
810                 ! Which of the lower interfaces is the highest (<=> the lowest pressure) ?
811                 zfrac = min(paprs(i,k), aviation_interface(kori))
812                 ! Which of the upper interfaces is the lowest (<=> the greatest pressure) ?
813                 zfrac = zfrac - max(paprs(i,k+1), aviation_interface(kori+1))
814                 ! If zfrac is negative, the layers are not overlapping
815                 ! Otherwise, we get the fraction of layer kori that overlap with layer k
816                 ! after normalisation to the total kori layer width
817                 zfrac = max(0.0, zfrac) / width_read_layer(kori)
818                 
819                 ! Vertical reprojection for each desired array
820                 flight_dist(i,k) = flight_dist(i,k) + zfrac * flight_dist_read(i,kori,1)
821                 flight_fuel(i,k) = flight_fuel(i,k) + zfrac * flight_fuel_read(i,kori,1)
822            ENDDO
823
824            !--Dry density [kg/m3]
825            rho = pplay(i,k) / temp(i,k) / RD
826            !--Dry air mass [kg/m2]
827            rhodz = ( paprs(i,k) - paprs(i,k+1) ) / RG
828            !--Cell thickness [m]
829            dz = rhodz / rho
830
831            !--Normalisation with the cell thickness
832            flight_dist(i,k) = flight_dist(i,k) / dz
833            flight_fuel(i,k) = flight_fuel(i,k) / dz
834           
835            !--Enhancement factor
836            flight_dist(i,k) = flight_dist(i,k) * aviation_coef
837            flight_fuel(i,k) = flight_fuel(i,k) * aviation_coef
838        ENDDO
839    ENDDO
840 
841END SUBROUTINE vertical_interpolation_aviation
842
843END MODULE lmdz_aviation
Note: See TracBrowser for help on using the repository browser.