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

Last change on this file since 5601 was 5601, checked in by aborella, 2 months ago

Multiple changes:

  • tracers which were ratios are now absolute quantities. This is needed because when the ratio

is not defined, some aberrations may occur

  • added a new tracer for total specific humidity in contrails
  • rework of the mixing process for cirrus clouds (and contrails)
  • changed the numerical integration of ice crystals' sublimation
  • subroutines do not take real inputs anymore (at least klon tables)
  • added more radiative diagnostics for contrails
File size: 31.7 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
12! flight_h2o_read  is the water content added to the air
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/mesh]
17!$OMP THREADPRIVATE(flight_dist_read)
18REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE, PRIVATE :: flight_h2o_read  ! Aviation H2O emitted within the mesh [kgH2O/s/mesh]
19!$OMP THREADPRIVATE(flight_h2o_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_h2o, d_q_avi &
28      )
29
30USE lmdz_lscp_ini, ONLY: RD
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_h2o ! aviation emitted H2O concentration [kgH2O/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
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
50    !--q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
51    ! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O ) 
52    !      = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) )
53    !--The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q)
54    !--flight_h2O is in kg H2O / s / m3
55   
56    !d_q_avi(i,k) = ( M_cell * qtot(i,k) + V_cell * flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
57    !             / ( M_cell             + V_cell * flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
58    !             - qtot(i,k)
59    !--NB., M_cell = V_cell * rho
60    !d_q_avi(i,k) = ( rho * qtot(i,k) + flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
61    !             / ( rho             + flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
62    !             - qtot(i,k)
63    !--Same formula, more computationally effective but less readable
64    d_q_avi(i,k) = flight_h2o(i,k) * ( 1. - qtot(i,k) ) &
65                 / ( rho / dtime / ( 1. - qtot(i,k) ) + flight_h2o(i,k) )
66  ENDDO
67ENDDO
68
69END SUBROUTINE aviation_water_emissions
70
71
72!**********************************************************************************
73SUBROUTINE contrails_formation( &
74      klon, dtime, pplay, temp, qsat, qsatl, gamma_cond, flight_dist, &
75      cldfra, pdf_loc, pdf_scale, pdf_alpha, keepgoing, &
76      Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
77      dcfa_ini, dqia_ini, dqta_ini)
78
79USE lmdz_lscp_ini, ONLY: RPI, RCPD, EPS_W, RTT
80USE lmdz_lscp_ini, ONLY: EI_H2O_aviation, qheat_fuel_aviation, prop_efficiency_aviation
81USE lmdz_lscp_ini, ONLY: eps, temp_nowater
82
83USE lmdz_lscp_tools, ONLY: GAMMAINC, calc_qsat_ecmwf
84
85IMPLICIT NONE
86
87!
88! Input
89!
90INTEGER,  INTENT(IN)                   :: klon         ! number of horizontal grid points
91REAL,     INTENT(IN)                   :: dtime        ! time step [s]
92REAL,     INTENT(IN) , DIMENSION(klon) :: pplay        ! layer pressure [Pa]
93REAL,     INTENT(IN) , DIMENSION(klon) :: temp         ! temperature [K]
94REAL,     INTENT(IN) , DIMENSION(klon) :: qsat         ! saturation specific humidity [kg/kg]
95REAL,     INTENT(IN) , DIMENSION(klon) :: qsatl        ! saturation specific humidity w.r.t. liquid [kg/kg]
96REAL,     INTENT(IN) , DIMENSION(klon) :: gamma_cond   ! condensation threshold w.r.t. qsat [-]
97REAL,     INTENT(IN) , DIMENSION(klon) :: flight_dist  ! aviation distance flown concentration [m/s/m3]
98REAL,     INTENT(IN) , DIMENSION(klon) :: cldfra       ! cloud fraction [-]
99REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_loc      ! location parameter of the clear sky PDF [%]
100REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_scale    ! scale parameter of the clear sky PDF [%]
101REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_alpha    ! shape parameter of the clear sky PDF [-]
102LOGICAL,  INTENT(IN) , DIMENSION(klon) :: keepgoing    ! .TRUE. if a new condensation loop should be computed
103!
104! Output
105!
106REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcritcont    ! critical temperature for contrail formation [K]
107REAL,     INTENT(INOUT), DIMENSION(klon) :: qcritcont    ! critical specific humidity for contrail formation [kg/kg]
108REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraP  ! potential persistent contrail fraction [-]
109REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraNP ! potential non-persistent contrail fraction [-]
110REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfa_ini     ! contrails cloud fraction tendency bec ause of initial formation [s-1]
111REAL,     INTENT(INOUT), DIMENSION(klon) :: dqia_ini     ! contrails ice specific humidity tende ncy because of initial formation [kg/kg/s]
112REAL,     INTENT(INOUT), DIMENSION(klon) :: dqta_ini     ! contrails total specific humidity ten dency because of initial formation [kg/kg/s]
113!
114! Local
115!
116INTEGER :: i
117! for Schmidt-Appleman criteria
118REAL, DIMENSION(klon) :: qzero
119REAL, DIMENSION(klon) :: qsatl_crit, dqsatl_crit
120REAL :: Gcont, psatl_crit, pcrit
121REAL :: pdf_x, pdf_y, pdf_e2, pdf_e3
122REAL :: pdf_fra_above_qcritcont, pdf_fra_above_qsat, pdf_fra_above_qnuc
123REAL :: pdf_q_above_qcritcont, pdf_q_above_qsat, pdf_q_above_qnuc
124REAL :: qpotcontP
125!
126! for new contrail formation
127REAL :: contrail_cross_section, contfra_new
128
129qzero(:) = 0.
130
131DO i = 1, klon
132  IF ( keepgoing(i) ) THEN
133    Tcritcont(i)    = 0.
134    qcritcont(i)    = 0.
135    potcontfraP(i)  = 0.
136    potcontfraNP(i) = 0.
137  ENDIF
138ENDDO
139
140!---------------------------------
141!--  SCHMIDT-APPLEMAN CRITERIA  --
142!---------------------------------
143!--Revised by Schumann (1996) and Rap et al. (2010)
144
145DO i = 1, klon
146  !--Gcont is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
147  !--in Pa / K. See Rap et al. (2010) in JGR.
148  !--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1)
149  Gcont = EI_H2O_aviation * RCPD * pplay(i) &
150         / ( EPS_W * qheat_fuel_aviation * ( 1. - prop_efficiency_aviation ) )
151  !--critical temperature below which no liquid contrail can form in exhaust
152  !--noted T_LM in Schumann (1996), their Eq. 31 in Appendix 2
153  !--in Kelvins
154  Tcritcont(i) = 226.69 + 9.43 * LOG( MAX(Gcont, 0.1) - 0.053 ) &
155         + 0.72 * LOG( MAX(Gcont, 0.1) - 0.053 )**2
156ENDDO
157
158CALL calc_qsat_ecmwf(klon, Tcritcont, qzero, pplay, RTT, 1, .FALSE., qsatl_crit, dqsatl_crit)
159
160DO i = 1, klon
161  IF ( keepgoing(i) .AND. ( temp(i) .LE. temp_nowater ) ) THEN
162   
163    psatl_crit = qsatl_crit(i) / ( EPS_W + ( 1. - EPS_W ) * qsatl_crit(i) ) * pplay(i)
164    pcrit = Gcont * ( temp(i) - Tcritcont(i) ) + psatl_crit
165    qcritcont(i) = EPS_W * pcrit / ( pplay(i) - ( 1. - EPS_W ) * pcrit )
166
167
168    IF ( ( ( 1. - cldfra(i) ) .GT. eps ) .AND. ( temp(i) .LT. Tcritcont(i) ) ) THEN
169   
170      pdf_x = qcritcont(i) / qsatl(i) * 100.
171      pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
172      pdf_e2 = GAMMA(1. + 1. / pdf_alpha(i))
173      pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
174      pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
175      pdf_fra_above_qcritcont = EXP( - pdf_y ) * ( 1. - cldfra(i) )
176      pdf_q_above_qcritcont = ( pdf_e3 * ( 1. - cldfra(i) ) &
177          + pdf_loc(i) * pdf_fra_above_qcritcont ) * qsatl(i) / 100.
178   
179      pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
180      pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
181      pdf_e2 = GAMMA(1. + 1. / pdf_alpha(i))
182      pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
183      pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
184      pdf_fra_above_qnuc = EXP( - pdf_y ) * ( 1. - cldfra(i) )
185      pdf_q_above_qnuc = ( pdf_e3 * ( 1. - cldfra(i) ) &
186          + pdf_loc(i) * pdf_fra_above_qnuc ) * qsatl(i) / 100.
187   
188      pdf_x = qsat(i) / qsatl(i) * 100.
189      pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
190      pdf_e2 = GAMMA(1. + 1. / pdf_alpha(i))
191      pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
192      pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
193      pdf_fra_above_qsat = EXP( - pdf_y ) * ( 1. - cldfra(i) )
194      pdf_q_above_qsat = ( pdf_e3 * ( 1. - cldfra(i) ) &
195          + pdf_loc(i) * pdf_fra_above_qsat ) * qsatl(i) / 100.
196   
197      potcontfraNP(i) = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat)
198      potcontfraP(i) = MAX(0., MIN(pdf_fra_above_qsat - pdf_fra_above_qnuc, &
199                        pdf_fra_above_qcritcont - pdf_fra_above_qnuc))
200      qpotcontP = MAX(0., MIN(pdf_q_above_qsat - pdf_q_above_qnuc, &
201                      pdf_q_above_qcritcont - pdf_q_above_qnuc))
202
203    ENDIF ! temp .LT. Tcritcont
204   
205    !--Add a source of contrails from aviation
206    IF ( potcontfraP(i) .GT. eps ) THEN
207      contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA()
208      contfra_new = MIN(1., flight_dist(i) * dtime * contrail_cross_section)
209      dcfa_ini(i) = potcontfraP(i) * contfra_new
210      dqta_ini(i) = qpotcontP * contfra_new
211      dqia_ini(i) = dqta_ini(i) - qsat(i) * dcfa_ini(i)
212    ENDIF
213
214  ENDIF ! keepgoing
215ENDDO
216
217END SUBROUTINE contrails_formation
218
219
220!**********************************************************************************
221FUNCTION QVAPINCLD_DEPSUB_CONTRAILS( &
222    qvapincld, qiceincld, temp, qsat, pplay, dtime)
223
224USE lmdz_lscp_ini,        ONLY: RV, RLSTT, RTT, EPS_W
225USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
226USE lmdz_lscp_ini,        ONLY: rm_ice_crystals_contrails
227
228IMPLICIT NONE
229
230! inputs
231REAL :: qvapincld     !
232REAL :: qiceincld     !
233REAL :: temp          !
234REAL :: qsat          !
235REAL :: pplay         !
236REAL :: dtime         ! time step [s]
237! output
238REAL :: qvapincld_depsub_contrails
239! local
240REAL :: pres_sat, rho, kappa
241REAL :: air_thermal_conduct, water_vapor_diff
242REAL :: rm_ice
243
244!--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
245!--hypothesis is lost, and the vapor in the cloud is purely prognostic.
246!
247!--The deposition equation is
248!-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
249!--from Lohmann et al. (2016), where
250!--alpha is the deposition coefficient [-]
251!--mi is the mass of one ice crystal [kg]
252!--C is the capacitance of an ice crystal [m]
253!--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
254!--R_v is the specific gas constant for humid air [J/kg/K]
255!--T is the temperature [K]
256!--esi is the saturation pressure w.r.t. ice [Pa]
257!--Dv is the diffusivity of water vapor [m2/s]
258!--Ls is the specific latent heat of sublimation [J/kg/K]
259!--ka is the thermal conductivity of dry air [J/m/s/K]
260!
261!--alpha is a coefficient to take into account the fact that during deposition, a water
262!--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
263!--coherent (with the same structure). It has no impact for sublimation.
264!--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
265!--during deposition, and alpha = 1. during sublimation.
266!--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
267!-- C = capa_cond_cirrus * rm_ice
268!
269!--We have qice = Nice * mi, where Nice is the ice crystal
270!--number concentration per kg of moist air
271!--HYPOTHESIS 1: the ice crystals are spherical, therefore
272!-- mi = 4/3 * pi * rm_ice**3 * rho_ice
273!--HYPOTHESIS 2: the ice crystals are monodisperse with the
274!--initial radius rm_ice_0.
275!--NB. this is notably different than the assumption
276!--of a distributed qice in the cloud made in the sublimation process;
277!--should it be consistent?
278!
279!--As the deposition process does not create new ice crystals,
280!--and because we assume a same rm_ice value for all crystals
281!--therefore the sublimation process does not destroy ice crystals
282!--(or, in a limit case, it destroys all ice crystals), then
283!--Nice is a constant during the sublimation/deposition process.
284!-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
285!
286!--The deposition equation then reads:
287!-- 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
288!-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
289!--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
290!--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
291!-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
292!--  *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 )
293!--and we have
294!-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
295!-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
296!--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
297!
298!--This system of equations can be resolved with an exact
299!--explicit numerical integration, having one variable resolved
300!--explicitly, the other exactly. The exactly resolved variable is
301!--the one decreasing, so it is qvc if the process is
302!--condensation, qi if it is sublimation.
303!
304!--kappa is computed as an initialisation constant, as it depends only
305!--on temperature and other pre-computed values
306pres_sat = qsat / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay
307!--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
308air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184
309!--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
310water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4
311kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat &
312      * ( RV * temp / water_vapor_diff / pres_sat &
313        + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) )
314!--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
315
316!--Here, the initial vapor in the cloud is qvapincld, and we compute
317!--the new vapor qvapincld_depsub_contrails
318
319rm_ice = rm_ice_crystals_contrails
320
321IF ( qvapincld .GE. qsat ) THEN
322  !--If the cloud is initially supersaturated
323  !--Exact explicit integration (qvc exact, qice explicit)
324  qvapincld_depsub_contrails = qsat + ( qvapincld - qsat ) &
325                * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
326ELSE
327  !--If the cloud is initially subsaturated
328  !--Exact explicit integration (qvc exact, qice explicit)
329  !--Same but depo_coef_cirrus = 1
330  qvapincld_depsub_contrails = qsat + ( qvapincld - qsat ) &
331                * EXP( - dtime * qiceincld / kappa / rm_ice**2 )
332ENDIF ! qvapincld .GT. qsat
333
334END FUNCTION QVAPINCLD_DEPSUB_CONTRAILS
335
336
337!**********************************************************************************
338FUNCTION contrail_cross_section_onera()
339
340USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails
341
342IMPLICIT NONE
343
344!
345! Output
346!
347REAL :: contrail_cross_section_onera ! [m2]
348!
349! Local
350!
351
352contrail_cross_section_onera = initial_width_contrails * initial_height_contrails
353
354END FUNCTION contrail_cross_section_onera
355
356SUBROUTINE read_aviation_emissions(klon, klev)
357    ! This subroutine allows to read the traffic density data read in the file aviation.nc
358    ! This file is defined in ./COMP/lmdz.card
359    ! Field names in context_input_lmdz.xml should be the same as in the file.
360    USE netcdf
361    USE mod_phys_lmdz_para
362    USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, nbp_lev_=>nbp_lev, &
363                                 klon_glo, grid2Dto1D_glo, grid_type, unstructured
364    USE iophy, ONLY : io_lon, io_lat
365    USE lmdz_xios
366    USE print_control_mod, ONLY: lunout
367    USE readaerosol_mod, ONLY: check_err
368    USE lmdz_lscp_ini, ONLY: EI_H2O_aviation
369    IMPLICIT NONE
370
371    INTEGER, INTENT(IN) :: klon, klev  ! number of horizontal grid points and vertical levels
372
373    !----------------------------------------------------
374    ! Local variable
375    !----------------------------------------------------
376    REAL, ALLOCATABLE :: flight_dist_mpi(:,:,:), flight_h2o_mpi(:,:,:)
377    INTEGER :: ierr
378
379    ! FOR REGULAR LON LAT
380    CHARACTER(len=30)     :: fname
381    INTEGER               :: nid, dimid, varid
382    INTEGER               :: imth, i, j, k
383    REAL                  :: npole, spole
384    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_dist_glo2D
385    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_h2o_glo2D
386    REAL, ALLOCATABLE, DIMENSION(:)       :: vartmp_lev
387    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: vartmp
388    REAL, ALLOCATABLE, DIMENSION(:)       :: lon_src              ! longitudes in file
389    REAL, ALLOCATABLE, DIMENSION(:)       :: lat_src, lat_src_inv ! latitudes in file
390    LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
391    INTEGER                               :: nbp_lon, nbp_lat
392
393
394IF (grid_type==unstructured) THEN
395
396    ! Get number of vertical levels and level values
397    IF (is_omp_master) CALL xios_get_axis_attr( "aviation_lev", n_glo=nleva )
398    CALL bcast_omp(nleva)
399
400    ! Allocation of arrays
401    ALLOCATE(flight_dist_read(klon, nleva,1), STAT=ierr)
402    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_dist',1)
403    ALLOCATE(flight_h2o_read(klon, nleva,1), STAT=ierr)
404    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_h2o',1)
405    ALLOCATE(aviation_lev(nleva), STAT=ierr)
406    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate aviation_lev',1)
407
408    ! Read the data from the file
409    ! is_omp_master is necessary to make XIOS works
410    IF (is_omp_master) THEN
411        ALLOCATE(flight_dist_mpi(klon_mpi, nleva,1), STAT=ierr)
412        IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_dist_mpi',1)
413        ALLOCATE(flight_h2o_mpi(klon_mpi, nleva,1), STAT=ierr)
414        IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_h2o_mpi',1)
415        CALL xios_recv_field("KMFLOWN_interp", flight_dist_mpi(:,:,1))
416        !CALL xios_recv_field("KGH2O_interp", flight_h2o_mpi(:,:,1))
417        flight_h2o_mpi(:,:,:) = 0.
418        ! Get number of vertical levels and level values
419        CALL xios_get_axis_attr( "aviation_lev", value=aviation_lev(:))
420    ENDIF
421
422    ! Propagate to other OMP threads: flight_dist_mpi(klon_mpi,klev) to flight_dist(klon,klev)
423    CALL scatter_omp(flight_dist_mpi, flight_dist_read)
424    CALL scatter_omp(flight_h2o_mpi, flight_h2o_read)
425    CALL bcast_omp(aviation_lev)
426
427ELSE
428
429    nbp_lon=nbp_lon_
430    nbp_lat=nbp_lat_
431   
432    IF (is_mpi_root) THEN
433      ALLOCATE(lon_src(nbp_lon))
434      ALLOCATE(lat_src(nbp_lat))
435      ALLOCATE(lat_src_inv(nbp_lat))
436    ELSE
437      ALLOCATE(flight_dist_glo2D(0,0,0,0))
438      ALLOCATE(flight_h2o_glo2D(0,0,0,0))
439    ENDIF
440
441    IF (is_mpi_root .AND. is_omp_root) THEN
442
443! 1) Open file
444!****************************************************************************************
445      fname = 'aviation.nc'
446      CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, nid), "pb open "//TRIM(fname) )
447
448
449! Test for equal longitudes and latitudes in file and model
450!****************************************************************************************
451      ! Read and test longitudes
452      CALL check_err( nf90_inq_varid(nid, 'longitude', varid), "pb inq lon" )
453      CALL check_err( nf90_get_var(nid, varid, lon_src(:)), "pb get lon" )
454     
455      IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
456         WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
457         WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
458         WRITE(lunout,*) 'longitudes in model :', io_lon
459       
460         CALL abort_physic('lmdz_aviation', 'longitudes are not the same in file and model',1)
461      END IF
462
463      ! Read and test latitudes
464      CALL check_err( nf90_inq_varid(nid, 'latitude', varid),"pb inq lat" )
465      CALL check_err( nf90_get_var(nid, varid, lat_src(:)),"pb get lat" )
466
467      ! Invert source latitudes
468      DO j = 1, nbp_lat
469         lat_src_inv(j) = lat_src(nbp_lat+1-j)
470      END DO
471
472      IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
473         ! Latitudes are the same
474         invert_lat=.FALSE.
475      ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
476         ! Inverted source latitudes correspond to model latitudes
477         WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
478         invert_lat=.TRUE.
479      ELSE
480         WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
481         WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src     
482         WRITE(lunout,*) 'latitudes in model :', io_lat
483         CALL abort_physic('lmdz_aviation', 'latitudes do not correspond between file and model',1)
484      END IF
485
486       
487! 2) Find vertical dimension nleva
488!****************************************************************************************
489       ierr = nf90_inq_dimid(nid, 'pressure_Pa', dimid)
490       CALL check_err( nf90_inquire_dimension(nid, dimid, len = nleva), "pb inq dim for PRESNIVS or lev" )
491       
492     ! Allocate variables depending on the number of vertical levels
493       ALLOCATE(flight_dist_glo2D(nbp_lon, nbp_lat, nleva, 1), flight_h2o_glo2D(nbp_lon, nbp_lat, nleva, 1), stat=ierr)
494       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 1',1)
495
496       ALLOCATE(aviation_lev(nleva), stat=ierr)
497       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 2',1)
498
499! 3) Read all variables from file
500!****************************************************************************************
501
502       ! Get variable id
503       CALL check_err( nf90_inq_varid(nid, 'seg_length_m', varid),"pb inq var seg_length_m" )
504       ! Get the variable
505       CALL check_err( nf90_get_var(nid, varid, flight_dist_glo2D),"pb get var seg_length_m"  )
506
507!       ! Get variable id
508!       CALL check_err( nf90_inq_varid(nid, 'fuel_burn', varid),"pb inq var fuel_burn" )
509!       ! Get the variable
510!       CALL check_err( nf90_get_var(nid, varid, flight_h2o_glo2D),"pb get var fuel_burn"  )
511       flight_h2o_glo2D(:,:,:,:) = 0.
512
513       ! Get variable id
514       CALL check_err( nf90_inq_varid(nid, "pressure_Pa", varid),"pb inq var pressure_Pa" )
515       ! Get the variable
516       CALL check_err( nf90_get_var(nid, varid, aviation_lev),"pb get var pressure_Pa" )
517         
518
519! 4) Close file 
520!****************************************************************************************
521       CALL check_err( nf90_close(nid), "pb in close" )
522     
523
524! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
525!****************************************************************************************
526! Test if vertical levels have to be inversed
527
528       IF (aviation_lev(1) < aviation_lev(nleva)) THEN
529         
530          ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva), vartmp_lev(nleva))
531
532          ! Inverse vertical levels for flight_dist_glo2D
533          vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1)
534          DO k=1, nleva
535             DO j=1, nbp_lat
536                DO i=1, nbp_lon
537                   flight_dist_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
538                END DO
539             END DO
540          END DO
541
542          ! Inverse vertical levels for flight_h2o_glo2D
543          vartmp(:,:,:) = flight_h2o_glo2D(:,:,:,1)
544          DO k=1, nleva
545             DO j=1, nbp_lat
546                DO i=1, nbp_lon
547                   flight_h2o_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
548                END DO
549             END DO
550          END DO
551           
552          ALLOCATE(vartmp_lev(nleva))
553          ! Inverte vertical axes for aviation_lev
554          vartmp_lev(:) = aviation_lev(:)
555          DO k=1, nleva
556             aviation_lev(k) = vartmp_lev(nleva+1-k)
557          END DO
558
559          DEALLOCATE(vartmp, vartmp_lev)
560
561       ELSE
562          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
563       END IF
564
565
566!     - Invert latitudes if necessary
567       IF (invert_lat) THEN
568
569          ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva))
570
571          ! Invert latitudes for the variable
572          vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1) ! use varmth temporarly
573          DO k=1,nleva
574             DO j=1,nbp_lat
575                DO i=1,nbp_lon
576                   flight_dist_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
577                END DO
578             END DO
579          END DO
580
581          ! Invert latitudes for the variable
582          vartmp(:,:,:) = flight_h2o_glo2D(:,:,:,1) ! use varmth temporarly
583          DO k=1,nleva
584             DO j=1,nbp_lat
585                DO i=1,nbp_lon
586                   flight_h2o_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
587                END DO
588             END DO
589          END DO
590
591          DEALLOCATE(vartmp)
592       END IF ! invert_lat
593       
594       ! Do zonal mead at poles and distribut at whole first and last latitude
595       DO k=1, nleva
596          npole=0.  ! North pole, j=1
597          spole=0.  ! South pole, j=nbp_lat       
598          DO i=1,nbp_lon
599             npole = npole + flight_dist_glo2D(i,1,k,1)
600             spole = spole + flight_dist_glo2D(i,nbp_lat,k,1)
601          END DO
602          npole = npole/REAL(nbp_lon)
603          spole = spole/REAL(nbp_lon)
604          flight_dist_glo2D(:,1,    k,1) = npole
605          flight_dist_glo2D(:,nbp_lat,k,1) = spole
606       END DO
607
608       ! Do zonal mead at poles and distribut at whole first and last latitude
609       DO k=1, nleva
610          npole=0.  ! North pole, j=1
611          spole=0.  ! South pole, j=nbp_lat       
612          DO i=1,nbp_lon
613             npole = npole + flight_h2o_glo2D(i,1,k,1)
614             spole = spole + flight_h2o_glo2D(i,nbp_lat,k,1)
615          END DO
616          npole = npole/REAL(nbp_lon)
617          spole = spole/REAL(nbp_lon)
618          flight_h2o_glo2D(:,1,    k,1) = npole
619          flight_h2o_glo2D(:,nbp_lat,k,1) = spole
620       END DO
621     
622       ALLOCATE(flight_dist_mpi(klon_glo, nleva, 1), flight_h2o_mpi(klon_glo, nleva, 1), stat=ierr)
623       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 3',1)
624     
625       ! Transform from 2D to 1D field
626       CALL grid2Dto1D_glo(flight_dist_glo2D, flight_dist_mpi)
627       CALL grid2Dto1D_glo(flight_h2o_glo2D, flight_h2o_mpi)
628   
629    ELSE
630        ALLOCATE(flight_dist_mpi(0,0,0), flight_h2o_mpi(0,0,0))
631    END IF ! is_mpi_root .AND. is_omp_root
632
633!$OMP BARRIER
634 
635! 6) Distribute to all processes
636!    Scatter global field(klon_glo) to local process domain(klon)
637!    and distribute nleva to all processes
638!****************************************************************************************
639
640    ! Distribute nleva
641    CALL bcast(nleva)
642
643    ! Allocate and distribute pt_ap and pt_b
644    IF (.NOT. ALLOCATED(aviation_lev)) THEN  ! if pt_ap is allocated also pt_b is allocated
645       ALLOCATE(aviation_lev(nleva), stat=ierr)
646       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 4',1)
647    END IF
648    CALL bcast(aviation_lev)
649
650    ! Allocate space for output pointer variable at local process
651    ALLOCATE(flight_dist_read(klon, nleva, 1), flight_h2o_read(klon, nleva, 1), stat=ierr)
652    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 5',1)
653
654    ! Scatter global field to local domain at local process
655    CALL scatter(flight_dist_mpi, flight_dist_read)
656    CALL scatter(flight_h2o_mpi, flight_h2o_read)
657
658ENDIF
659
660END SUBROUTINE read_aviation_emissions
661
662SUBROUTINE vertical_interpolation_aviation(klon, klev, paprs, pplay, temp, flight_dist, flight_h2o)
663    ! This subroutine performs the vertical interpolation from the read data in aviation.nc
664    ! where there are nleva vertical levels described in aviation_lev to the klev levels or
665    ! the model.
666    ! flight_dist_read(klon,nleva) -> flight_dist(klon, klev)
667    ! flight_h2o_read(klon,nleva) -> flight_h2o(klon, klev)
668
669    USE lmdz_lscp_ini, ONLY: RD, RG
670    USE lmdz_lscp_ini, ONLY: aviation_coef
671
672    IMPLICIT NONE
673
674    INTEGER,                    INTENT(IN)  :: klon, klev  ! number of horizontal grid points and vertical levels
675    REAL, INTENT(IN)    :: paprs(klon, klev+1) ! inter-layer pressure [Pa]
676    REAL, INTENT(IN)    :: pplay(klon, klev) ! mid-layer pressure [Pa]
677    REAL, INTENT(IN)    :: temp(klon, klev) ! temperature [K]
678    REAL, INTENT(OUT)   :: flight_dist(klon,klev) ! Aviation distance flown within the mesh [m/s/mesh]
679    REAL, INTENT(OUT)   :: flight_h2o(klon,klev)  ! Aviation H2O emitted within the mesh [kgH2O/s/mesh]
680
681    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
682    !  Local variable
683    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
684    REAL :: aviation_interface(1:nleva+1) ! Pressure of aviation file interfaces [ Pa ]
685    INTEGER :: k, kori  ! Loop index for vertical layers
686    INTEGER :: i  ! Loop index for horizontal grid
687    REAL :: zfrac ! Fraction of layer kori in layer k
688    REAL :: width_read_layer(1:nleva) ! width of a given layer [ Pa ]
689    REAL :: rho, rhodz, dz
690
691    ! Initialisation
692    flight_dist(:,:) = 0.
693    flight_h2o(:,:) = 0.
694
695    ! Compute the array with the vertical interface
696    ! It starts at 1 and has length nleva + 1
697    ! Note that aviation_lev has nleva and gives the altitude in the middle of the layers
698    ! Surface pressure in standard atmosphere model [ Pa ]
699    aviation_interface(1) = 101325.
700    DO kori=2, nleva
701        aviation_interface(kori) = (aviation_lev(kori-1)+aviation_lev(kori))/2.0  ! [ Pa ]
702    ENDDO
703    ! Last interface - we assume the same spacing as the very last one
704    aviation_interface(nleva+1) = aviation_interface(nleva) - (aviation_lev(nleva-1) - aviation_lev(nleva))
705
706    ! Vertical width of each layer of the read file
707    ! It is positive
708    DO kori=1, nleva
709        width_read_layer(kori) = aviation_interface(kori) - aviation_interface(kori+1)
710    ENDDO
711
712    ! Vertical reprojection
713    ! The loop over klon is induced since it is done by MPI threads
714    ! zfrac is the fraction of layer kori (read file) included in layer k (model)
715    DO i=1,klon
716        DO k=1, klev
717            DO kori=1,nleva
718                 ! Which of the lower interfaces is the highest (<=> the lowest pressure) ?
719                 zfrac = min(paprs(i,k), aviation_interface(kori))
720                 ! Which of the upper interfaces is the lowest (<=> the greatest pressure) ?
721                 zfrac = zfrac - max(paprs(i,k+1), aviation_interface(kori+1))
722                 ! If zfrac is negative, the layers are not overlapping
723                 ! Otherwise, we get the fraction of layer kori that overlap with layer k
724                 ! after normalisation to the total kori layer width
725                 zfrac = max(0.0, zfrac) / width_read_layer(kori)
726                 
727                 ! Vertical reprojection for each desired array
728                 flight_dist(i,k) = flight_dist(i,k) + zfrac * flight_dist_read(i,kori,1)
729                 flight_h2o(i,k)  = flight_h2o(i,k) + zfrac * flight_h2o_read(i,kori,1)
730            ENDDO
731
732            !--Dry density [kg/m3]
733            rho = pplay(i,k) / temp(i,k) / RD
734            !--Dry air mass [kg/m2]
735            rhodz = ( paprs(i,k) - paprs(i,k+1) ) / RG
736            !--Cell thickness [m]
737            dz = rhodz / rho
738
739            !--Normalisation with the cell thickness
740            flight_dist(i,k) = flight_dist(i,k) / dz
741            flight_h2o(i,k) = flight_h2o(i,k) / dz
742           
743            !--Enhancement factor
744            flight_dist(i,k) = flight_dist(i,k) * aviation_coef
745            flight_h2o(i,k) = flight_h2o(i,k) * aviation_coef
746        ENDDO
747    ENDDO
748 
749END SUBROUTINE vertical_interpolation_aviation
750
751END MODULE lmdz_aviation
Note: See TracBrowser for help on using the repository browser.