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

Last change on this file since 5625 was 5623, checked in by aborella, 7 months ago

Small bug corrections + temporary patch for children tracers in ICO (back to full parents tracers) + now truly no aviation data needed to run
File size: 26.8 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, eps ) / pdf_scale(i) ) * 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, eps ) / pdf_scale(i) ) * 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 read_aviation_emissions(klon, klev)
251    ! This subroutine allows to read the traffic density data read in the file aviation.nc
252    ! This file is defined in ./COMP/lmdz.card
253    ! Field names in context_input_lmdz.xml should be the same as in the file.
254    USE netcdf
255    USE mod_phys_lmdz_para
256    USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, nbp_lev_=>nbp_lev, &
257                                 klon_glo, grid2Dto1D_glo, grid_type, unstructured
258    USE iophy, ONLY : io_lon, io_lat
259    USE lmdz_xios
260    USE print_control_mod, ONLY: lunout
261    USE readaerosol_mod, ONLY: check_err
262    USE lmdz_lscp_ini, ONLY: EI_H2O_aviation
263    IMPLICIT NONE
264
265    INTEGER, INTENT(IN) :: klon, klev  ! number of horizontal grid points and vertical levels
266
267    !----------------------------------------------------
268    ! Local variable
269    !----------------------------------------------------
270    REAL, ALLOCATABLE :: flight_dist_mpi(:,:,:), flight_h2o_mpi(:,:,:)
271    INTEGER :: ierr
272
273    ! FOR REGULAR LON LAT
274    CHARACTER(len=30)     :: fname
275    INTEGER               :: nid, dimid, varid
276    INTEGER               :: imth, i, j, k
277    REAL                  :: npole, spole
278    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_dist_glo2D
279    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_h2o_glo2D
280    REAL, ALLOCATABLE, DIMENSION(:)       :: vartmp_lev
281    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: vartmp
282    REAL, ALLOCATABLE, DIMENSION(:)       :: lon_src              ! longitudes in file
283    REAL, ALLOCATABLE, DIMENSION(:)       :: lat_src, lat_src_inv ! latitudes in file
284    LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
285    INTEGER                               :: nbp_lon, nbp_lat
286
287
288IF (grid_type==unstructured) THEN
289
290    IF (is_omp_master) THEN
291      ! Activate aviation file
292      CALL xios_set_file_attr("aviation", enabled=.TRUE.)
293      ! Activate aviation fields
294      CALL xios_set_field_attr("KMFLOWN_read", enabled=.TRUE.)
295      CALL xios_set_field_attr("KMFLOWN_interp", enabled=.TRUE.)
296
297      ! Get number of vertical levels and level values
298      CALL xios_get_axis_attr( "aviation_lev", n_glo=nleva )
299    ENDIF
300    CALL bcast_omp(nleva)
301
302    ! Allocation of arrays
303    ALLOCATE(flight_dist_read(klon, nleva,1), STAT=ierr)
304    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_dist',1)
305    ALLOCATE(flight_h2o_read(klon, nleva,1), STAT=ierr)
306    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_h2o',1)
307    ALLOCATE(aviation_lev(nleva), STAT=ierr)
308    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate aviation_lev',1)
309
310    ! Read the data from the file
311    ! is_omp_master is necessary to make XIOS works
312    IF (is_omp_master) THEN
313        ALLOCATE(flight_dist_mpi(klon_mpi, nleva,1), STAT=ierr)
314        IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_dist_mpi',1)
315        ALLOCATE(flight_h2o_mpi(klon_mpi, nleva,1), STAT=ierr)
316        IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_h2o_mpi',1)
317        CALL xios_recv_field("KMFLOWN_interp", flight_dist_mpi(:,:,1))
318        !CALL xios_recv_field("KGH2O_interp", flight_h2o_mpi(:,:,1))
319        flight_h2o_mpi(:,:,:) = 0.
320        ! Get number of vertical levels and level values
321        CALL xios_get_axis_attr( "aviation_lev", value=aviation_lev(:))
322    ENDIF
323
324    ! Propagate to other OMP threads: flight_dist_mpi(klon_mpi,klev) to flight_dist(klon,klev)
325    CALL scatter_omp(flight_dist_mpi, flight_dist_read)
326    CALL scatter_omp(flight_h2o_mpi, flight_h2o_read)
327    CALL bcast_omp(aviation_lev)
328
329ELSE
330
331    nbp_lon=nbp_lon_
332    nbp_lat=nbp_lat_
333   
334    IF (is_mpi_root) THEN
335      ALLOCATE(lon_src(nbp_lon))
336      ALLOCATE(lat_src(nbp_lat))
337      ALLOCATE(lat_src_inv(nbp_lat))
338    ELSE
339      ALLOCATE(flight_dist_glo2D(0,0,0,0))
340      ALLOCATE(flight_h2o_glo2D(0,0,0,0))
341    ENDIF
342
343    IF (is_mpi_root .AND. is_omp_root) THEN
344
345! 1) Open file
346!****************************************************************************************
347      fname = 'aviation.nc'
348      CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, nid), "pb open "//TRIM(fname) )
349
350
351! Test for equal longitudes and latitudes in file and model
352!****************************************************************************************
353      ! Read and test longitudes
354      CALL check_err( nf90_inq_varid(nid, 'longitude', varid), "pb inq lon" )
355      CALL check_err( nf90_get_var(nid, varid, lon_src(:)), "pb get lon" )
356     
357      IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
358         WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
359         WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
360         WRITE(lunout,*) 'longitudes in model :', io_lon
361       
362         CALL abort_physic('lmdz_aviation', 'longitudes are not the same in file and model',1)
363      END IF
364
365      ! Read and test latitudes
366      CALL check_err( nf90_inq_varid(nid, 'latitude', varid),"pb inq lat" )
367      CALL check_err( nf90_get_var(nid, varid, lat_src(:)),"pb get lat" )
368
369      ! Invert source latitudes
370      DO j = 1, nbp_lat
371         lat_src_inv(j) = lat_src(nbp_lat+1-j)
372      END DO
373
374      IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
375         ! Latitudes are the same
376         invert_lat=.FALSE.
377      ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
378         ! Inverted source latitudes correspond to model latitudes
379         WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
380         invert_lat=.TRUE.
381      ELSE
382         WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
383         WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src     
384         WRITE(lunout,*) 'latitudes in model :', io_lat
385         CALL abort_physic('lmdz_aviation', 'latitudes do not correspond between file and model',1)
386      END IF
387
388       
389! 2) Find vertical dimension nleva
390!****************************************************************************************
391       ierr = nf90_inq_dimid(nid, 'pressure_Pa', dimid)
392       CALL check_err( nf90_inquire_dimension(nid, dimid, len = nleva), "pb inq dim for PRESNIVS or lev" )
393       
394     ! Allocate variables depending on the number of vertical levels
395       ALLOCATE(flight_dist_glo2D(nbp_lon, nbp_lat, nleva, 1), flight_h2o_glo2D(nbp_lon, nbp_lat, nleva, 1), stat=ierr)
396       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 1',1)
397
398       ALLOCATE(aviation_lev(nleva), stat=ierr)
399       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 2',1)
400
401! 3) Read all variables from file
402!****************************************************************************************
403
404       ! Get variable id
405       CALL check_err( nf90_inq_varid(nid, 'seg_length_m', varid),"pb inq var seg_length_m" )
406       ! Get the variable
407       CALL check_err( nf90_get_var(nid, varid, flight_dist_glo2D),"pb get var seg_length_m"  )
408
409!       ! Get variable id
410!       CALL check_err( nf90_inq_varid(nid, 'fuel_burn', varid),"pb inq var fuel_burn" )
411!       ! Get the variable
412!       CALL check_err( nf90_get_var(nid, varid, flight_h2o_glo2D),"pb get var fuel_burn"  )
413       flight_h2o_glo2D(:,:,:,:) = 0.
414
415       ! Get variable id
416       CALL check_err( nf90_inq_varid(nid, "pressure_Pa", varid),"pb inq var pressure_Pa" )
417       ! Get the variable
418       CALL check_err( nf90_get_var(nid, varid, aviation_lev),"pb get var pressure_Pa" )
419         
420
421! 4) Close file 
422!****************************************************************************************
423       CALL check_err( nf90_close(nid), "pb in close" )
424     
425
426! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
427!****************************************************************************************
428! Test if vertical levels have to be inversed
429
430       IF (aviation_lev(1) < aviation_lev(nleva)) THEN
431         
432          ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva), vartmp_lev(nleva))
433
434          ! Inverse vertical levels for flight_dist_glo2D
435          vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1)
436          DO k=1, nleva
437             DO j=1, nbp_lat
438                DO i=1, nbp_lon
439                   flight_dist_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
440                END DO
441             END DO
442          END DO
443
444          ! Inverse vertical levels for flight_h2o_glo2D
445          vartmp(:,:,:) = flight_h2o_glo2D(:,:,:,1)
446          DO k=1, nleva
447             DO j=1, nbp_lat
448                DO i=1, nbp_lon
449                   flight_h2o_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
450                END DO
451             END DO
452          END DO
453           
454          ALLOCATE(vartmp_lev(nleva))
455          ! Inverte vertical axes for aviation_lev
456          vartmp_lev(:) = aviation_lev(:)
457          DO k=1, nleva
458             aviation_lev(k) = vartmp_lev(nleva+1-k)
459          END DO
460
461          DEALLOCATE(vartmp, vartmp_lev)
462
463       ELSE
464          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
465       END IF
466
467
468!     - Invert latitudes if necessary
469       IF (invert_lat) THEN
470
471          ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva))
472
473          ! Invert latitudes for the variable
474          vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1) ! use varmth temporarly
475          DO k=1,nleva
476             DO j=1,nbp_lat
477                DO i=1,nbp_lon
478                   flight_dist_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
479                END DO
480             END DO
481          END DO
482
483          ! Invert latitudes for the variable
484          vartmp(:,:,:) = flight_h2o_glo2D(:,:,:,1) ! use varmth temporarly
485          DO k=1,nleva
486             DO j=1,nbp_lat
487                DO i=1,nbp_lon
488                   flight_h2o_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
489                END DO
490             END DO
491          END DO
492
493          DEALLOCATE(vartmp)
494       END IF ! invert_lat
495       
496       ! Do zonal mead at poles and distribut at whole first and last latitude
497       DO k=1, nleva
498          npole=0.  ! North pole, j=1
499          spole=0.  ! South pole, j=nbp_lat       
500          DO i=1,nbp_lon
501             npole = npole + flight_dist_glo2D(i,1,k,1)
502             spole = spole + flight_dist_glo2D(i,nbp_lat,k,1)
503          END DO
504          npole = npole/REAL(nbp_lon)
505          spole = spole/REAL(nbp_lon)
506          flight_dist_glo2D(:,1,    k,1) = npole
507          flight_dist_glo2D(:,nbp_lat,k,1) = spole
508       END DO
509
510       ! Do zonal mead at poles and distribut at whole first and last latitude
511       DO k=1, nleva
512          npole=0.  ! North pole, j=1
513          spole=0.  ! South pole, j=nbp_lat       
514          DO i=1,nbp_lon
515             npole = npole + flight_h2o_glo2D(i,1,k,1)
516             spole = spole + flight_h2o_glo2D(i,nbp_lat,k,1)
517          END DO
518          npole = npole/REAL(nbp_lon)
519          spole = spole/REAL(nbp_lon)
520          flight_h2o_glo2D(:,1,    k,1) = npole
521          flight_h2o_glo2D(:,nbp_lat,k,1) = spole
522       END DO
523     
524       ALLOCATE(flight_dist_mpi(klon_glo, nleva, 1), flight_h2o_mpi(klon_glo, nleva, 1), stat=ierr)
525       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 3',1)
526     
527       ! Transform from 2D to 1D field
528       CALL grid2Dto1D_glo(flight_dist_glo2D, flight_dist_mpi)
529       CALL grid2Dto1D_glo(flight_h2o_glo2D, flight_h2o_mpi)
530   
531    ELSE
532        ALLOCATE(flight_dist_mpi(0,0,0), flight_h2o_mpi(0,0,0))
533    END IF ! is_mpi_root .AND. is_omp_root
534
535!$OMP BARRIER
536 
537! 6) Distribute to all processes
538!    Scatter global field(klon_glo) to local process domain(klon)
539!    and distribute nleva to all processes
540!****************************************************************************************
541
542    ! Distribute nleva
543    CALL bcast(nleva)
544
545    ! Allocate and distribute pt_ap and pt_b
546    IF (.NOT. ALLOCATED(aviation_lev)) THEN  ! if pt_ap is allocated also pt_b is allocated
547       ALLOCATE(aviation_lev(nleva), stat=ierr)
548       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 4',1)
549    END IF
550    CALL bcast(aviation_lev)
551
552    ! Allocate space for output pointer variable at local process
553    ALLOCATE(flight_dist_read(klon, nleva, 1), flight_h2o_read(klon, nleva, 1), stat=ierr)
554    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 5',1)
555
556    ! Scatter global field to local domain at local process
557    CALL scatter(flight_dist_mpi, flight_dist_read)
558    CALL scatter(flight_h2o_mpi, flight_h2o_read)
559
560ENDIF
561
562END SUBROUTINE read_aviation_emissions
563
564SUBROUTINE vertical_interpolation_aviation(klon, klev, paprs, pplay, temp, flight_dist, flight_h2o)
565    ! This subroutine performs the vertical interpolation from the read data in aviation.nc
566    ! where there are nleva vertical levels described in aviation_lev to the klev levels or
567    ! the model.
568    ! flight_dist_read(klon,nleva) -> flight_dist(klon, klev)
569    ! flight_h2o_read(klon,nleva) -> flight_h2o(klon, klev)
570
571    USE lmdz_lscp_ini, ONLY: RD, RG
572    USE lmdz_lscp_ini, ONLY: aviation_coef
573
574    IMPLICIT NONE
575
576    INTEGER,                    INTENT(IN)  :: klon, klev  ! number of horizontal grid points and vertical levels
577    REAL, INTENT(IN)    :: paprs(klon, klev+1) ! inter-layer pressure [Pa]
578    REAL, INTENT(IN)    :: pplay(klon, klev) ! mid-layer pressure [Pa]
579    REAL, INTENT(IN)    :: temp(klon, klev) ! temperature [K]
580    REAL, INTENT(OUT)   :: flight_dist(klon,klev) ! Aviation distance flown within the mesh [m/s/mesh]
581    REAL, INTENT(OUT)   :: flight_h2o(klon,klev)  ! Aviation H2O emitted within the mesh [kgH2O/s/mesh]
582
583    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
584    !  Local variable
585    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
586    REAL :: aviation_interface(1:nleva+1) ! Pressure of aviation file interfaces [ Pa ]
587    INTEGER :: k, kori  ! Loop index for vertical layers
588    INTEGER :: i  ! Loop index for horizontal grid
589    REAL :: zfrac ! Fraction of layer kori in layer k
590    REAL :: width_read_layer(1:nleva) ! width of a given layer [ Pa ]
591    REAL :: rho, rhodz, dz
592
593    ! Initialisation
594    flight_dist(:,:) = 0.
595    flight_h2o(:,:) = 0.
596
597    ! Compute the array with the vertical interface
598    ! It starts at 1 and has length nleva + 1
599    ! Note that aviation_lev has nleva and gives the altitude in the middle of the layers
600    ! Surface pressure in standard atmosphere model [ Pa ]
601    aviation_interface(1) = 101325.
602    DO kori=2, nleva
603        aviation_interface(kori) = (aviation_lev(kori-1)+aviation_lev(kori))/2.0  ! [ Pa ]
604    ENDDO
605    ! Last interface - we assume the same spacing as the very last one
606    aviation_interface(nleva+1) = aviation_interface(nleva) - (aviation_lev(nleva-1) - aviation_lev(nleva))
607
608    ! Vertical width of each layer of the read file
609    ! It is positive
610    DO kori=1, nleva
611        width_read_layer(kori) = aviation_interface(kori) - aviation_interface(kori+1)
612    ENDDO
613
614    ! Vertical reprojection
615    ! The loop over klon is induced since it is done by MPI threads
616    ! zfrac is the fraction of layer kori (read file) included in layer k (model)
617    DO i=1,klon
618        DO k=1, klev
619            DO kori=1,nleva
620                 ! Which of the lower interfaces is the highest (<=> the lowest pressure) ?
621                 zfrac = min(paprs(i,k), aviation_interface(kori))
622                 ! Which of the upper interfaces is the lowest (<=> the greatest pressure) ?
623                 zfrac = zfrac - max(paprs(i,k+1), aviation_interface(kori+1))
624                 ! If zfrac is negative, the layers are not overlapping
625                 ! Otherwise, we get the fraction of layer kori that overlap with layer k
626                 ! after normalisation to the total kori layer width
627                 zfrac = max(0.0, zfrac) / width_read_layer(kori)
628                 
629                 ! Vertical reprojection for each desired array
630                 flight_dist(i,k) = flight_dist(i,k) + zfrac * flight_dist_read(i,kori,1)
631                 flight_h2o(i,k)  = flight_h2o(i,k) + zfrac * flight_h2o_read(i,kori,1)
632            ENDDO
633
634            !--Dry density [kg/m3]
635            rho = pplay(i,k) / temp(i,k) / RD
636            !--Dry air mass [kg/m2]
637            rhodz = ( paprs(i,k) - paprs(i,k+1) ) / RG
638            !--Cell thickness [m]
639            dz = rhodz / rho
640
641            !--Normalisation with the cell thickness
642            flight_dist(i,k) = flight_dist(i,k) / dz
643            flight_h2o(i,k) = flight_h2o(i,k) / dz
644           
645            !--Enhancement factor
646            flight_dist(i,k) = flight_dist(i,k) * aviation_coef
647            flight_h2o(i,k) = flight_h2o(i,k) * aviation_coef
648        ENDDO
649    ENDDO
650 
651END SUBROUTINE vertical_interpolation_aviation
652
653END MODULE lmdz_aviation
Note: See TracBrowser for help on using the repository browser.