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

Last change on this file since 5637 was 5637, checked in by aborella, 8 months ago

Fix for activating aviation emissions data file
File size: 27.3 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      clrfra, qclr, pdf_scale, pdf_alpha, pdf_gamma, keepgoing, pt_pron_clds, &
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) :: clrfra       ! clear sky fraction [-]
99REAL,     INTENT(IN) , DIMENSION(klon) :: qclr         ! clear sky specific humidity [kg/kg]
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 [-]
102REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_gamma    ! tmp parameter of the clear sky PDF [-]
103LOGICAL,  INTENT(IN) , DIMENSION(klon) :: keepgoing    ! .TRUE. if a new condensation loop should be computed
104LOGICAL,  INTENT(IN) , DIMENSION(klon) :: pt_pron_clds  ! .TRUE. if clouds are prognostic in this mesh
105!
106! Output
107!
108REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcritcont    ! critical temperature for contrail formation [K]
109REAL,     INTENT(INOUT), DIMENSION(klon) :: qcritcont    ! critical specific humidity for contrail formation [kg/kg]
110REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraP  ! potential persistent contrail fraction [-]
111REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraNP ! potential non-persistent contrail fraction [-]
112REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfa_ini     ! contrails cloud fraction tendency bec ause of initial formation [s-1]
113REAL,     INTENT(INOUT), DIMENSION(klon) :: dqia_ini     ! contrails ice specific humidity tende ncy because of initial formation [kg/kg/s]
114REAL,     INTENT(INOUT), DIMENSION(klon) :: dqta_ini     ! contrails total specific humidity ten dency because of initial formation [kg/kg/s]
115!
116! Local
117!
118INTEGER :: i
119! for Schmidt-Appleman criteria
120REAL, DIMENSION(klon) :: qzero
121REAL, DIMENSION(klon) :: qsatl_crit, dqsatl_crit
122REAL :: Gcont, psatl_crit, pcrit
123REAL :: rhl_clr, pdf_loc
124REAL :: pdf_x, pdf_y, pdf_e3
125REAL :: pdf_fra_above_qcritcont, pdf_fra_above_qsat
126REAL :: pdf_q_above_qcritcont, pdf_q_above_qsat
127REAL :: qpotcontP
128!
129! for new contrail formation
130REAL :: contrail_cross_section, contfra_new
131
132qzero(:) = 0.
133
134DO i = 1, klon
135  IF ( keepgoing(i) ) THEN
136    Tcritcont(i)    = 0.
137    qcritcont(i)    = 0.
138    potcontfraP(i)  = 0.
139    potcontfraNP(i) = 0.
140  ENDIF
141ENDDO
142
143!---------------------------------
144!--  SCHMIDT-APPLEMAN CRITERIA  --
145!---------------------------------
146!--Revised by Schumann (1996) and Rap et al. (2010)
147
148DO i = 1, klon
149  !--Gcont is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
150  !--in Pa / K. See Rap et al. (2010) in JGR.
151  !--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1)
152  Gcont = EI_H2O_aviation * RCPD * pplay(i) &
153         / ( EPS_W * qheat_fuel_aviation * ( 1. - prop_efficiency_aviation ) )
154  !--critical temperature below which no liquid contrail can form in exhaust
155  !--noted T_LM in Schumann (1996), their Eq. 31 in Appendix 2
156  !--in Kelvins
157  Tcritcont(i) = 226.69 + 9.43 * LOG( MAX(Gcont, 0.1) - 0.053 ) &
158         + 0.72 * LOG( MAX(Gcont, 0.1) - 0.053 )**2
159ENDDO
160
161CALL calc_qsat_ecmwf(klon, Tcritcont, qzero, pplay, RTT, 1, .FALSE., qsatl_crit, dqsatl_crit)
162
163DO i = 1, klon
164  IF ( keepgoing(i) .AND. pt_pron_clds(i) .AND. ( temp(i) .LE. temp_nowater ) ) THEN
165   
166    psatl_crit = qsatl_crit(i) / ( EPS_W + ( 1. - EPS_W ) * qsatl_crit(i) ) * pplay(i)
167    pcrit = Gcont * ( temp(i) - Tcritcont(i) ) + psatl_crit
168    qcritcont(i) = EPS_W * pcrit / ( pplay(i) - ( 1. - EPS_W ) * pcrit )
169
170
171    IF ( ( clrfra(i) .GT. eps ) .AND. ( temp(i) .LT. Tcritcont(i) ) ) THEN
172   
173      rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
174      pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
175      pdf_x = qcritcont(i) / qsatl(i) * 100.
176      pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
177      IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
178        pdf_fra_above_qcritcont = 0.
179        pdf_q_above_qcritcont = 0.
180      ELSEIF ( pdf_y .LT. -10. ) THEN
181        pdf_fra_above_qcritcont = clrfra(i)
182        pdf_q_above_qcritcont = qclr(i)
183      ELSE
184        pdf_y = EXP( pdf_y )
185        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
186        pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
187        pdf_fra_above_qcritcont = EXP( - pdf_y ) * clrfra(i)
188        pdf_q_above_qcritcont = ( pdf_e3 * clrfra(i) &
189                              + pdf_loc * pdf_fra_above_qcritcont ) * qsatl(i) / 100.
190      ENDIF
191   
192      pdf_x = qsat(i) / qsatl(i) * 100.
193      pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
194      IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
195        pdf_fra_above_qsat = 0.
196        pdf_q_above_qsat = 0.
197      ELSEIF ( pdf_y .LT. -10. ) THEN
198        pdf_fra_above_qsat = clrfra(i)
199        pdf_q_above_qsat = qclr(i)
200      ELSE
201        pdf_y = EXP( pdf_y )
202        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
203        pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
204        pdf_fra_above_qsat = EXP( - pdf_y ) * clrfra(i)
205        pdf_q_above_qsat = ( pdf_e3 * clrfra(i) &
206                         + pdf_loc * pdf_fra_above_qsat ) * qsatl(i) / 100.
207      ENDIF
208   
209      potcontfraNP(i) = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat)
210      potcontfraP(i) = MIN(pdf_fra_above_qsat, pdf_fra_above_qcritcont)
211      qpotcontP = MIN(pdf_q_above_qsat, pdf_q_above_qcritcont)
212
213    ENDIF ! temp .LT. Tcritcont
214   
215    !--Add a source of contrails from aviation
216    IF ( potcontfraP(i) .GT. eps ) THEN
217      contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA()
218      contfra_new = MIN(1., flight_dist(i) * dtime * contrail_cross_section)
219      dcfa_ini(i) = potcontfraP(i) * contfra_new
220      dqta_ini(i) = qpotcontP * contfra_new
221      dqia_ini(i) = dqta_ini(i) - qsat(i) * dcfa_ini(i)
222    ENDIF
223
224  ENDIF ! keepgoing
225ENDDO
226
227END SUBROUTINE contrails_formation
228
229
230!**********************************************************************************
231FUNCTION contrail_cross_section_onera()
232
233USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails
234
235IMPLICIT NONE
236
237!
238! Output
239!
240REAL :: contrail_cross_section_onera ! [m2]
241!
242! Local
243!
244
245contrail_cross_section_onera = initial_width_contrails * initial_height_contrails
246
247END FUNCTION contrail_cross_section_onera
248
249
250SUBROUTINE init_read_aviation_emissions
251
252  USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured
253  USE mod_phys_lmdz_para, ONLY: is_omp_master
254  USE lmdz_xios, ONLY: xios_set_file_attr, xios_set_field_attr
255
256  IF (grid_type==unstructured) THEN
257    IF (is_omp_master) THEN
258      ! Activate aviation file
259      CALL xios_set_file_attr("aviation_file", enabled=.TRUE.)
260      ! Activate aviation fields
261      CALL xios_set_field_attr("KMFLOWN_read", enabled=.TRUE.)
262      CALL xios_set_field_attr("KMFLOWN_interp", enabled=.TRUE.)
263    ENDIF
264  ENDIF
265
266END SUBROUTINE init_read_aviation_emissions
267
268
269SUBROUTINE read_aviation_emissions(klon, klev)
270    ! This subroutine allows to read the traffic density data read in the file aviation.nc
271    ! This file is defined in ./COMP/lmdz.card
272    ! Field names in context_input_lmdz.xml should be the same as in the file.
273    USE netcdf
274    USE mod_phys_lmdz_para
275    USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, nbp_lev_=>nbp_lev, &
276                                 klon_glo, grid2Dto1D_glo, grid_type, unstructured
277    USE iophy, ONLY : io_lon, io_lat
278    USE lmdz_xios
279    USE print_control_mod, ONLY: lunout
280    USE readaerosol_mod, ONLY: check_err
281    USE lmdz_lscp_ini, ONLY: EI_H2O_aviation
282    IMPLICIT NONE
283
284    INTEGER, INTENT(IN) :: klon, klev  ! number of horizontal grid points and vertical levels
285
286    !----------------------------------------------------
287    ! Local variable
288    !----------------------------------------------------
289    REAL, ALLOCATABLE :: flight_dist_mpi(:,:,:), flight_h2o_mpi(:,:,:)
290    INTEGER :: ierr
291    LOGICAL, SAVE :: first = .TRUE.
292    !$OMP THREADPRIVATE(first)
293
294    ! FOR REGULAR LON LAT
295    CHARACTER(len=30)     :: fname
296    INTEGER               :: nid, dimid, varid
297    INTEGER               :: imth, i, j, k
298    REAL                  :: npole, spole
299    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_dist_glo2D
300    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_h2o_glo2D
301    REAL, ALLOCATABLE, DIMENSION(:)       :: vartmp_lev
302    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: vartmp
303    REAL, ALLOCATABLE, DIMENSION(:)       :: lon_src              ! longitudes in file
304    REAL, ALLOCATABLE, DIMENSION(:)       :: lat_src, lat_src_inv ! latitudes in file
305    LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
306    INTEGER                               :: nbp_lon, nbp_lat
307
308
309IF (grid_type==unstructured) THEN
310
311    IF (first) THEN
312      IF (is_omp_master) THEN
313        ! Get number of vertical levels and level values
314        CALL xios_get_axis_attr("aviation_lev", n_glo=nleva)
315      ENDIF
316      CALL bcast_omp(nleva)
317
318      ! Allocation of arrays
319      ALLOCATE(flight_dist_read(klon, nleva, 1), STAT=ierr)
320      IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_dist',1)
321      ALLOCATE(flight_h2o_read(klon, nleva, 1), STAT=ierr)
322      IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_h2o',1)
323      ALLOCATE(aviation_lev(nleva), STAT=ierr)
324      IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate aviation_lev',1)
325
326      IF (is_omp_master) THEN
327        ! Get number of vertical levels and level values
328        CALL xios_get_axis_attr("aviation_lev", value=aviation_lev(:))
329      ENDIF
330      CALL bcast_omp(aviation_lev)
331      first = .FALSE.
332    ENDIF ! first
333
334    ! Read the data from the file
335    ! is_omp_master is necessary to make XIOS works
336    IF (is_omp_master) THEN
337        ALLOCATE(flight_dist_mpi(klon_mpi, nleva,1), STAT=ierr)
338        IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_dist_mpi',1)
339        ALLOCATE(flight_h2o_mpi(klon_mpi, nleva,1), STAT=ierr)
340        IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_h2o_mpi',1)
341        CALL xios_recv_field("KMFLOWN_interp", flight_dist_mpi(:,:,1))
342        !CALL xios_recv_field("KGH2O_interp", flight_h2o_mpi(:,:,1))
343        flight_h2o_mpi(:,:,:) = 0.
344    ENDIF
345
346    ! Propagate to other OMP threads: flight_dist_mpi(klon_mpi,klev) to flight_dist(klon,klev)
347    CALL scatter_omp(flight_dist_mpi, flight_dist_read)
348    CALL scatter_omp(flight_h2o_mpi, flight_h2o_read)
349
350ELSE
351
352    nbp_lon=nbp_lon_
353    nbp_lat=nbp_lat_
354   
355    IF (is_mpi_root) THEN
356      ALLOCATE(lon_src(nbp_lon))
357      ALLOCATE(lat_src(nbp_lat))
358      ALLOCATE(lat_src_inv(nbp_lat))
359    ELSE
360      ALLOCATE(flight_dist_glo2D(0,0,0,0))
361      ALLOCATE(flight_h2o_glo2D(0,0,0,0))
362    ENDIF
363
364    IF (is_mpi_root .AND. is_omp_root) THEN
365
366! 1) Open file
367!****************************************************************************************
368      fname = 'aviation.nc'
369      CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, nid), "pb open "//TRIM(fname) )
370
371
372! Test for equal longitudes and latitudes in file and model
373!****************************************************************************************
374      ! Read and test longitudes
375      CALL check_err( nf90_inq_varid(nid, 'longitude', varid), "pb inq lon" )
376      CALL check_err( nf90_get_var(nid, varid, lon_src(:)), "pb get lon" )
377     
378      IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
379         WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
380         WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
381         WRITE(lunout,*) 'longitudes in model :', io_lon
382       
383         CALL abort_physic('lmdz_aviation', 'longitudes are not the same in file and model',1)
384      END IF
385
386      ! Read and test latitudes
387      CALL check_err( nf90_inq_varid(nid, 'latitude', varid),"pb inq lat" )
388      CALL check_err( nf90_get_var(nid, varid, lat_src(:)),"pb get lat" )
389
390      ! Invert source latitudes
391      DO j = 1, nbp_lat
392         lat_src_inv(j) = lat_src(nbp_lat+1-j)
393      END DO
394
395      IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
396         ! Latitudes are the same
397         invert_lat=.FALSE.
398      ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
399         ! Inverted source latitudes correspond to model latitudes
400         WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
401         invert_lat=.TRUE.
402      ELSE
403         WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
404         WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src     
405         WRITE(lunout,*) 'latitudes in model :', io_lat
406         CALL abort_physic('lmdz_aviation', 'latitudes do not correspond between file and model',1)
407      END IF
408
409       
410! 2) Find vertical dimension nleva
411!****************************************************************************************
412       ierr = nf90_inq_dimid(nid, 'pressure_Pa', dimid)
413       CALL check_err( nf90_inquire_dimension(nid, dimid, len = nleva), "pb inq dim for PRESNIVS or lev" )
414       
415     ! Allocate variables depending on the number of vertical levels
416       ALLOCATE(flight_dist_glo2D(nbp_lon, nbp_lat, nleva, 1), flight_h2o_glo2D(nbp_lon, nbp_lat, nleva, 1), stat=ierr)
417       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 1',1)
418
419       ALLOCATE(aviation_lev(nleva), stat=ierr)
420       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 2',1)
421
422! 3) Read all variables from file
423!****************************************************************************************
424
425       ! Get variable id
426       CALL check_err( nf90_inq_varid(nid, 'seg_length_m', varid),"pb inq var seg_length_m" )
427       ! Get the variable
428       CALL check_err( nf90_get_var(nid, varid, flight_dist_glo2D),"pb get var seg_length_m"  )
429
430!       ! Get variable id
431!       CALL check_err( nf90_inq_varid(nid, 'fuel_burn', varid),"pb inq var fuel_burn" )
432!       ! Get the variable
433!       CALL check_err( nf90_get_var(nid, varid, flight_h2o_glo2D),"pb get var fuel_burn"  )
434       flight_h2o_glo2D(:,:,:,:) = 0.
435
436       ! Get variable id
437       CALL check_err( nf90_inq_varid(nid, "pressure_Pa", varid),"pb inq var pressure_Pa" )
438       ! Get the variable
439       CALL check_err( nf90_get_var(nid, varid, aviation_lev),"pb get var pressure_Pa" )
440         
441
442! 4) Close file 
443!****************************************************************************************
444       CALL check_err( nf90_close(nid), "pb in close" )
445     
446
447! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
448!****************************************************************************************
449! Test if vertical levels have to be inversed
450
451       IF (aviation_lev(1) < aviation_lev(nleva)) THEN
452         
453          ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva), vartmp_lev(nleva))
454
455          ! Inverse vertical levels for flight_dist_glo2D
456          vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1)
457          DO k=1, nleva
458             DO j=1, nbp_lat
459                DO i=1, nbp_lon
460                   flight_dist_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
461                END DO
462             END DO
463          END DO
464
465          ! Inverse vertical levels for flight_h2o_glo2D
466          vartmp(:,:,:) = flight_h2o_glo2D(:,:,:,1)
467          DO k=1, nleva
468             DO j=1, nbp_lat
469                DO i=1, nbp_lon
470                   flight_h2o_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
471                END DO
472             END DO
473          END DO
474           
475          ALLOCATE(vartmp_lev(nleva))
476          ! Inverte vertical axes for aviation_lev
477          vartmp_lev(:) = aviation_lev(:)
478          DO k=1, nleva
479             aviation_lev(k) = vartmp_lev(nleva+1-k)
480          END DO
481
482          DEALLOCATE(vartmp, vartmp_lev)
483
484       ELSE
485          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
486       END IF
487
488
489!     - Invert latitudes if necessary
490       IF (invert_lat) THEN
491
492          ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva))
493
494          ! Invert latitudes for the variable
495          vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1) ! use varmth temporarly
496          DO k=1,nleva
497             DO j=1,nbp_lat
498                DO i=1,nbp_lon
499                   flight_dist_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
500                END DO
501             END DO
502          END DO
503
504          ! Invert latitudes for the variable
505          vartmp(:,:,:) = flight_h2o_glo2D(:,:,:,1) ! use varmth temporarly
506          DO k=1,nleva
507             DO j=1,nbp_lat
508                DO i=1,nbp_lon
509                   flight_h2o_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
510                END DO
511             END DO
512          END DO
513
514          DEALLOCATE(vartmp)
515       END IF ! invert_lat
516       
517       ! Do zonal mead at poles and distribut at whole first and last latitude
518       DO k=1, nleva
519          npole=0.  ! North pole, j=1
520          spole=0.  ! South pole, j=nbp_lat       
521          DO i=1,nbp_lon
522             npole = npole + flight_dist_glo2D(i,1,k,1)
523             spole = spole + flight_dist_glo2D(i,nbp_lat,k,1)
524          END DO
525          npole = npole/REAL(nbp_lon)
526          spole = spole/REAL(nbp_lon)
527          flight_dist_glo2D(:,1,    k,1) = npole
528          flight_dist_glo2D(:,nbp_lat,k,1) = spole
529       END DO
530
531       ! Do zonal mead at poles and distribut at whole first and last latitude
532       DO k=1, nleva
533          npole=0.  ! North pole, j=1
534          spole=0.  ! South pole, j=nbp_lat       
535          DO i=1,nbp_lon
536             npole = npole + flight_h2o_glo2D(i,1,k,1)
537             spole = spole + flight_h2o_glo2D(i,nbp_lat,k,1)
538          END DO
539          npole = npole/REAL(nbp_lon)
540          spole = spole/REAL(nbp_lon)
541          flight_h2o_glo2D(:,1,    k,1) = npole
542          flight_h2o_glo2D(:,nbp_lat,k,1) = spole
543       END DO
544     
545       ALLOCATE(flight_dist_mpi(klon_glo, nleva, 1), flight_h2o_mpi(klon_glo, nleva, 1), stat=ierr)
546       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 3',1)
547     
548       ! Transform from 2D to 1D field
549       CALL grid2Dto1D_glo(flight_dist_glo2D, flight_dist_mpi)
550       CALL grid2Dto1D_glo(flight_h2o_glo2D, flight_h2o_mpi)
551   
552    ELSE
553        ALLOCATE(flight_dist_mpi(0,0,0), flight_h2o_mpi(0,0,0))
554    END IF ! is_mpi_root .AND. is_omp_root
555
556!$OMP BARRIER
557 
558! 6) Distribute to all processes
559!    Scatter global field(klon_glo) to local process domain(klon)
560!    and distribute nleva to all processes
561!****************************************************************************************
562
563    ! Distribute nleva
564    CALL bcast(nleva)
565
566    ! Allocate and distribute pt_ap and pt_b
567    IF (.NOT. ALLOCATED(aviation_lev)) THEN  ! if pt_ap is allocated also pt_b is allocated
568       ALLOCATE(aviation_lev(nleva), stat=ierr)
569       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 4',1)
570    END IF
571    CALL bcast(aviation_lev)
572
573    ! Allocate space for output pointer variable at local process
574    ALLOCATE(flight_dist_read(klon, nleva, 1), flight_h2o_read(klon, nleva, 1), stat=ierr)
575    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 5',1)
576
577    ! Scatter global field to local domain at local process
578    CALL scatter(flight_dist_mpi, flight_dist_read)
579    CALL scatter(flight_h2o_mpi, flight_h2o_read)
580
581ENDIF
582
583END SUBROUTINE read_aviation_emissions
584
585SUBROUTINE vertical_interpolation_aviation(klon, klev, paprs, pplay, temp, flight_dist, flight_h2o)
586    ! This subroutine performs the vertical interpolation from the read data in aviation.nc
587    ! where there are nleva vertical levels described in aviation_lev to the klev levels or
588    ! the model.
589    ! flight_dist_read(klon,nleva) -> flight_dist(klon, klev)
590    ! flight_h2o_read(klon,nleva) -> flight_h2o(klon, klev)
591
592    USE lmdz_lscp_ini, ONLY: RD, RG
593    USE lmdz_lscp_ini, ONLY: aviation_coef
594
595    IMPLICIT NONE
596
597    INTEGER,                    INTENT(IN)  :: klon, klev  ! number of horizontal grid points and vertical levels
598    REAL, INTENT(IN)    :: paprs(klon, klev+1) ! inter-layer pressure [Pa]
599    REAL, INTENT(IN)    :: pplay(klon, klev) ! mid-layer pressure [Pa]
600    REAL, INTENT(IN)    :: temp(klon, klev) ! temperature [K]
601    REAL, INTENT(OUT)   :: flight_dist(klon,klev) ! Aviation distance flown within the mesh [m/s/mesh]
602    REAL, INTENT(OUT)   :: flight_h2o(klon,klev)  ! Aviation H2O emitted within the mesh [kgH2O/s/mesh]
603
604    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
605    !  Local variable
606    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
607    REAL :: aviation_interface(1:nleva+1) ! Pressure of aviation file interfaces [ Pa ]
608    INTEGER :: k, kori  ! Loop index for vertical layers
609    INTEGER :: i  ! Loop index for horizontal grid
610    REAL :: zfrac ! Fraction of layer kori in layer k
611    REAL :: width_read_layer(1:nleva) ! width of a given layer [ Pa ]
612    REAL :: rho, rhodz, dz
613
614    ! Initialisation
615    flight_dist(:,:) = 0.
616    flight_h2o(:,:) = 0.
617
618    ! Compute the array with the vertical interface
619    ! It starts at 1 and has length nleva + 1
620    ! Note that aviation_lev has nleva and gives the altitude in the middle of the layers
621    ! Surface pressure in standard atmosphere model [ Pa ]
622    aviation_interface(1) = 101325.
623    DO kori=2, nleva
624        aviation_interface(kori) = (aviation_lev(kori-1)+aviation_lev(kori))/2.0  ! [ Pa ]
625    ENDDO
626    ! Last interface - we assume the same spacing as the very last one
627    aviation_interface(nleva+1) = aviation_interface(nleva) - (aviation_lev(nleva-1) - aviation_lev(nleva))
628
629    ! Vertical width of each layer of the read file
630    ! It is positive
631    DO kori=1, nleva
632        width_read_layer(kori) = aviation_interface(kori) - aviation_interface(kori+1)
633    ENDDO
634
635    ! Vertical reprojection
636    ! The loop over klon is induced since it is done by MPI threads
637    ! zfrac is the fraction of layer kori (read file) included in layer k (model)
638    DO i=1,klon
639        DO k=1, klev
640            DO kori=1,nleva
641                 ! Which of the lower interfaces is the highest (<=> the lowest pressure) ?
642                 zfrac = min(paprs(i,k), aviation_interface(kori))
643                 ! Which of the upper interfaces is the lowest (<=> the greatest pressure) ?
644                 zfrac = zfrac - max(paprs(i,k+1), aviation_interface(kori+1))
645                 ! If zfrac is negative, the layers are not overlapping
646                 ! Otherwise, we get the fraction of layer kori that overlap with layer k
647                 ! after normalisation to the total kori layer width
648                 zfrac = max(0.0, zfrac) / width_read_layer(kori)
649                 
650                 ! Vertical reprojection for each desired array
651                 flight_dist(i,k) = flight_dist(i,k) + zfrac * flight_dist_read(i,kori,1)
652                 flight_h2o(i,k)  = flight_h2o(i,k) + zfrac * flight_h2o_read(i,kori,1)
653            ENDDO
654
655            !--Dry density [kg/m3]
656            rho = pplay(i,k) / temp(i,k) / RD
657            !--Dry air mass [kg/m2]
658            rhodz = ( paprs(i,k) - paprs(i,k+1) ) / RG
659            !--Cell thickness [m]
660            dz = rhodz / rho
661
662            !--Normalisation with the cell thickness
663            flight_dist(i,k) = flight_dist(i,k) / dz
664            flight_h2o(i,k) = flight_h2o(i,k) / dz
665           
666            !--Enhancement factor
667            flight_dist(i,k) = flight_dist(i,k) * aviation_coef
668            flight_h2o(i,k) = flight_h2o(i,k) * aviation_coef
669        ENDDO
670    ENDDO
671 
672END SUBROUTINE vertical_interpolation_aviation
673
674END MODULE lmdz_aviation
Note: See TracBrowser for help on using the repository browser.