MODULE lmdz_aviation !---------------------------------------------------------------- ! Module for aviation and contrails IMPLICIT NONE ! Arrays for the lecture of aviation files ! The allocation is done in the read_aviation module ! The size is (klon, nleva, 1) where ! nleva is the size of the vertical axis (read from file) ! flight_dist_read is the number of km per second per m2 ! flight_fuel_read is the fuel consumed per second per m2 ! aviation_lev is the value of the levels INTEGER, SAVE :: nleva ! Size of the vertical axis in the file !$OMP THREADPRIVATE(nleva) REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE, PRIVATE :: flight_dist_read ! Aviation distance flown within the mesh [m/s/m2/vertmesh] !$OMP THREADPRIVATE(flight_dist_read) REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE, PRIVATE :: flight_fuel_read ! Aviation fuel consumed within the mesh [kg/s/m2/vertmesh] !$OMP THREADPRIVATE(flight_fuel_read) REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE :: aviation_lev ! Pressure in the middle of the layers [Pa] !$OMP THREADPRIVATE(aviation_lev) CONTAINS SUBROUTINE aviation_water_emissions( & klon, klev, dtime, pplay, temp, qtot, & flight_fuel, d_q_avi & ) USE lmdz_lscp_ini, ONLY: RD, EI_H2O_aviation IMPLICIT NONE INTEGER, INTENT(IN) :: klon, klev ! number of horizontal grid points and vertical levels REAL, INTENT(IN) :: dtime ! time step [s] REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa] REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature (K) REAL, DIMENSION(klon,klev), INTENT(IN) :: qtot ! total specific humidity (in vapor phase) [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: flight_fuel ! aviation fuel consumption concentration [kg/s/m3] REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q_avi ! water vapor tendency from aviation [kg/kg] ! Local INTEGER :: i, k REAL :: rho, flight_h2o DO i=1, klon DO k=1, klev !--Dry density [kg/m3] rho = pplay(i,k) / temp(i,k) / RD flight_h2o = flight_fuel(i,k) * EI_H2O_aviation !--q is the specific humidity (kg/kg humid air) hence the complicated equation to update q ! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O ) ! = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) ) !--The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q) !--flight_h2O is in kg H2O / s / m3 !d_q_avi(i,k) = ( M_cell * qtot(i,k) + V_cell * flight_h2o * dtime * ( 1. - qtot(i,k) ) ) & ! / ( M_cell + V_cell * flight_h2o * dtime * ( 1. - qtot(i,k) ) ) & ! - qtot(i,k) !--NB., M_cell = V_cell * rho !d_q_avi(i,k) = ( rho * qtot(i,k) + flight_h2o * dtime * ( 1. - qtot(i,k) ) ) & ! / ( rho + flight_h2o * dtime * ( 1. - qtot(i,k) ) ) & ! - qtot(i,k) !--Same formula, more computationally effective but less readable d_q_avi(i,k) = flight_h2o * ( 1. - qtot(i,k) ) & / ( rho / dtime / ( 1. - qtot(i,k) ) + flight_h2o ) ENDDO ENDDO END SUBROUTINE aviation_water_emissions !********************************************************************************** SUBROUTINE contrails_formation( & klon, dtime, pplay, temp, rho, qsat, qsatl, gamma_cond, flight_dist, flight_fuel, & clrfra, qclr, pdf_scale, pdf_alpha, pdf_gamma, keepgoing, pt_pron_clds, & Tcritcont, qcritcont, potcontfraP, potcontfraNP, & AEI_contrails, AEI_surv_contrails, fsurv_contrails, section_contrails, & dcfc_ini, dqic_ini, dqtc_ini, dnic_ini) USE lmdz_lscp_ini, ONLY: RPI, RCPD, EPS_W, RTT USE lmdz_lscp_ini, ONLY: EI_H2O_aviation, qheat_fuel_aviation, prop_efficiency_aviation USE lmdz_lscp_ini, ONLY: eps, temp_nowater USE lmdz_lscp_ini, ONLY: Nice_init_contrails USE lmdz_lscp_tools, ONLY: GAMMAINC, calc_qsat_ecmwf IMPLICIT NONE ! ! Input ! INTEGER, INTENT(IN) :: klon ! number of horizontal grid points REAL, INTENT(IN) :: dtime ! time step [s] REAL, INTENT(IN) , DIMENSION(klon) :: pplay ! layer pressure [Pa] REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K] REAL, INTENT(IN) , DIMENSION(klon) :: rho ! air density [kg/m3] REAL, INTENT(IN) , DIMENSION(klon) :: qsat ! saturation specific humidity [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: qsatl ! saturation specific humidity w.r.t. liquid [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-] REAL, INTENT(IN) , DIMENSION(klon) :: flight_dist ! aviation distance flown concentration [m/s/m3] REAL, INTENT(IN) , DIMENSION(klon) :: flight_fuel ! aviation fuel consumption concentration [kg/s/m3] REAL, INTENT(IN) , DIMENSION(klon) :: clrfra ! clear sky fraction [-] REAL, INTENT(IN) , DIMENSION(klon) :: qclr ! clear sky specific humidity [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: pdf_scale ! scale parameter of the clear sky PDF [%] REAL, INTENT(IN) , DIMENSION(klon) :: pdf_alpha ! shape parameter of the clear sky PDF [-] REAL, INTENT(IN) , DIMENSION(klon) :: pdf_gamma ! tmp parameter of the clear sky PDF [-] LOGICAL, INTENT(IN) , DIMENSION(klon) :: keepgoing ! .TRUE. if a new condensation loop should be computed LOGICAL, INTENT(IN) , DIMENSION(klon) :: pt_pron_clds ! .TRUE. if clouds are prognostic in this mesh ! ! Output ! REAL, INTENT(INOUT), DIMENSION(klon) :: Tcritcont ! critical temperature for contrail formation [K] REAL, INTENT(INOUT), DIMENSION(klon) :: qcritcont ! critical specific humidity for contrail formation [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: potcontfraP ! potential persistent contrail fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: potcontfraNP ! potential non-persistent contrail fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: AEI_contrails ! Apparent emission index contrails [#/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: AEI_surv_contrails ! Apparent emission index contrails after vortex loss [#/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: fsurv_contrails ! Survival fraction after vortex loss [-] REAL, INTENT(INOUT), DIMENSION(klon) :: section_contrails ! Cross section of newly formed contrails [m2] REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_ini ! contrails cloud fraction tendency bec ause of initial formation [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_ini ! contrails ice specific humidity tende ncy because of initial formation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_ini ! contrails total specific humidity ten dency because of initial formation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dnic_ini ! contrails ice crystals concentration ten dency because of initial formation [#/kg/s] ! ! Local ! INTEGER :: i ! for Schmidt-Appleman criteria REAL, DIMENSION(klon) :: qzero REAL, DIMENSION(klon) :: qsatl_crit, dqsatl_crit REAL :: Gcont, psatl_crit, pcrit REAL :: rhl_clr, pdf_loc REAL :: pdf_x, pdf_y, pdf_e3 REAL :: pdf_fra_above_qcritcont, pdf_fra_above_qsat REAL :: pdf_q_above_qcritcont, pdf_q_above_qsat REAL :: qpotcontP ! ! for new contrail formation REAL :: dist_contrails, Nice_per_m_init_contrails REAL :: icesat_ratio qzero(:) = 0. DO i = 1, klon IF ( keepgoing(i) ) THEN qcritcont(i) = 0. potcontfraP(i) = 0. potcontfraNP(i) = 0. ENDIF ENDDO !--------------------------------- !-- SCHMIDT-APPLEMAN CRITERIA -- !--------------------------------- !--Revised by Schumann (1996) and Rap et al. (2010) DO i = 1, klon !--Gcont is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute !--in Pa / K. See Rap et al. (2010) in JGR. !--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1) Gcont = EI_H2O_aviation * RCPD * pplay(i) & / ( EPS_W * qheat_fuel_aviation * ( 1. - prop_efficiency_aviation ) ) !--critical temperature below which no liquid contrail can form in exhaust !--noted T_LM in Schumann (1996), their Eq. 31 in Appendix 2 !--in Kelvins Tcritcont(i) = 226.69 + 9.43 * LOG( MAX(Gcont, 0.1) - 0.053 ) & + 0.72 * LOG( MAX(Gcont, 0.1) - 0.053 )**2 ENDDO CALL calc_qsat_ecmwf(klon, Tcritcont, qzero, pplay, RTT, 1, .FALSE., qsatl_crit, dqsatl_crit) DO i = 1, klon IF ( keepgoing(i) .AND. pt_pron_clds(i) ) THEN psatl_crit = qsatl_crit(i) / ( EPS_W + ( 1. - EPS_W ) * qsatl_crit(i) ) * pplay(i) pcrit = Gcont * ( temp(i) - Tcritcont(i) ) + psatl_crit qcritcont(i) = EPS_W * pcrit / ( pplay(i) - ( 1. - EPS_W ) * pcrit ) IF ( ( clrfra(i) .GT. eps ) .AND. ( temp(i) .LT. Tcritcont(i) ) ) THEN rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qcritcont(i) / qsatl(i) * 100. pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows pdf_fra_above_qcritcont = 0. pdf_q_above_qcritcont = 0. ELSEIF ( pdf_y .LT. -10. ) THEN pdf_fra_above_qcritcont = clrfra(i) pdf_q_above_qcritcont = qclr(i) ELSE pdf_y = EXP( pdf_y ) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) pdf_fra_above_qcritcont = EXP( - pdf_y ) * clrfra(i) pdf_q_above_qcritcont = ( pdf_e3 * clrfra(i) & + pdf_loc * pdf_fra_above_qcritcont ) * qsatl(i) / 100. ENDIF pdf_x = qsat(i) / qsatl(i) * 100. pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows pdf_fra_above_qsat = 0. pdf_q_above_qsat = 0. ELSEIF ( pdf_y .LT. -10. ) THEN pdf_fra_above_qsat = clrfra(i) pdf_q_above_qsat = qclr(i) ELSE pdf_y = EXP( pdf_y ) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) pdf_fra_above_qsat = EXP( - pdf_y ) * clrfra(i) pdf_q_above_qsat = ( pdf_e3 * clrfra(i) & + pdf_loc * pdf_fra_above_qsat ) * qsatl(i) / 100. ENDIF potcontfraNP(i) = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat) potcontfraP(i) = MIN(pdf_fra_above_qsat, pdf_fra_above_qcritcont) qpotcontP = MIN(pdf_q_above_qsat, pdf_q_above_qcritcont) ENDIF ! temp .LT. Tcritcont !--Add a source of contrails from aviation IF ( ( potcontfraP(i) .GT. eps ) .AND. ( flight_dist(i) .GT. 1e-20 ) ) THEN !section_contrails(i) = CONTRAIL_CROSS_SECTION_ONERA() section_contrails(i) = CONTRAIL_CROSS_SECTION_SCHUMANN( & dtime, rho(i), flight_dist(i), flight_fuel(i)) icesat_ratio = qpotcontP / potcontfraP(i) / qsat(i) !--If Nice init is fixed in the plume !Nice_per_m_init_contrails = Nice_init_contrails * 1e6 * section_contrails(i) !--Else if it is parameterized CALL CONTRAIL_ICE_NUMBER_CONCENTRATION(temp(i), icesat_ratio, rho(i), & flight_dist(i), flight_fuel(i), Nice_per_m_init_contrails, & AEI_contrails(i), AEI_surv_contrails(i), fsurv_contrails(i)) !--If Nice_per_m_init_contrails .EQ. 0., all the crystals were lost in the vortex IF ( Nice_per_m_init_contrails .GT. 0. ) THEN !--dist_contrails is in meters of contrails per m3 (gridbox-average) dist_contrails = flight_dist(i) * dtime * potcontfraP(i) dcfc_ini(i) = dist_contrails * section_contrails(i) dqtc_ini(i) = icesat_ratio * qsat(i) * dcfc_ini(i) dqic_ini(i) = dqtc_ini(i) - qsat(i) * dcfc_ini(i) dnic_ini(i) = Nice_per_m_init_contrails * dist_contrails / rho(i) ENDIF ENDIF ENDIF ! keepgoing ENDDO END SUBROUTINE contrails_formation !********************************************************************************** FUNCTION CONTRAIL_CROSS_SECTION_ONERA() USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails IMPLICIT NONE ! ! Output ! REAL :: contrail_cross_section_onera ! [m2] ! ! Local ! contrail_cross_section_onera = initial_width_contrails * initial_height_contrails END FUNCTION CONTRAIL_CROSS_SECTION_ONERA !********************************************************************************** !--Based on Schumann (1998) FUNCTION CONTRAIL_CROSS_SECTION_SCHUMANN(dtime, rho_air, flight_dist, flight_fuel) USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails IMPLICIT NONE ! ! Input ! REAL :: dtime ! timestep [s] REAL :: rho_air ! air density [kg/m3] REAL :: flight_dist ! flown distance [m/s/m3] REAL :: flight_fuel ! fuel consumed [kg/s/m3] ! ! Output ! REAL :: contrail_cross_section_schumann ! [m2] ! ! Local ! REAL :: dilution_factor !--The contrail is formed on average in the middle of the timestep dilution_factor = 7000. * (dtime / 2.)**0.8 contrail_cross_section_schumann = flight_fuel / flight_dist * dilution_factor / rho_air END FUNCTION CONTRAIL_CROSS_SECTION_SCHUMANN !********************************************************************************** SUBROUTINE CONTRAIL_ICE_NUMBER_CONCENTRATION(temp, icesat_ratio, rho_air, & flight_dist, flight_fuel, Nice_per_m_init_contrails, & AEI_contrails, AEI_surv_contrails, fsurv_contrails) USE lmdz_lscp_ini, ONLY: RPI USE lmdz_lscp_ini, ONLY: EI_soot_aviation, air_to_fuel_ratio_engine, wingspan USE lmdz_lscp_ini, ONLY: Naer_amb, raer_amb_mean, raer_amb_std, r_soot_mean, r_soot_std USE lmdz_lscp_ini, ONLY: circ_0_loss, plume_area_loss, nice_init_ref_loss USE lmdz_lscp_ini, ONLY: N_Brunt_Vaisala_aviation, EI_H2O_aviation IMPLICIT NONE ! ! Input ! REAL, INTENT(IN) :: temp ! air temperature [K] REAL, INTENT(IN) :: icesat_ratio ! ice saturation ratio [-] REAL, INTENT(IN) :: rho_air ! air density [kg/m3] REAL, INTENT(IN) :: flight_dist ! flown distance [m/s/m3] REAL, INTENT(IN) :: flight_fuel ! fuel consumed [kg/s/m3] ! ! Output ! REAL, INTENT(OUT) :: Nice_per_m_init_contrails ! [#/m] REAL, INTENT(OUT) :: AEI_contrails ! [#/kg] REAL, INTENT(OUT) :: AEI_surv_contrails ! [#/kg] REAL, INTENT(OUT) :: fsurv_contrails ! [-] ! ! Local ! ! For initial ice nucleation REAL :: log_liqsat_ratio, coef_expo, dil_coef REAL :: rd_act_amb, rd_act_soot, phi_amb, phi_soot REAL :: AEI_soot, AEI_amb ! For ice crystals loss REAL :: rho_emit, temp_205, nice_init REAL :: z_atm, z_emit, z_desc, z_delta REAL :: Nice_per_m, fuel_per_m, frac_surv !------------------------------ !-- INITIAL ICE NUCLEATION -- !------------------------------ ! !--Karcher et al. (2015), Bier and Burkhardt (2019, 2022) !log_liqsat_ratio = LOG(liq_satratio) !! dry core radius in nm !! HERE SHOULD IT BE liqsat_ratio OR liqsat_ratio - 1. ? !rd_act_amb = (4. / 27. / LOG(liqsat_ratio)**2 / 0.5)**(1./3.) !! Integrate lognormal distribution between rd_act_amb and +inf !coef_expo = 4. / SQRT(2. * RPI) / LOG(raer_amb_std) !phi_amb = 1. / (1. + (rd_act_amb / raer_amb_mean)**coef_expo) ! !! BU22, Eq. A1, dry core radius in nm !rd_act_soot = 0.96453426 + 1.21152973 / log_liqsat_ratio - 0.00520358 / log_liqsat_ratio**2 & ! + 2.32286432e-5 / log_liqsat_ratio**3 - 4.36979055e-8 / log_liqsat_ratio**4 !rd_act_soot = MIN(150., MAX(1., rd_act_soot)) !! Integrate lognormal distribution between rd_act_amb and +inf !coef_expo = 4. / SQRT(2. * RPI) / LOG(r_soot_std) !phi_amb = 1. / (1. + (rd_act_soot / r_soot_mean)**coef_expo) ! !dil_coef = (0.01 / t0)**0.9 ! !! BU22, Eq. 5b !AEI_soot = phi_soot * EI_soot_aviation !AEI_amb = phi_amb * air_to_fuel_ratio_engine * (1. - dil_coef) / dil_coef & ! / rho_air * Naer_amb * 1e6 !AEI_contrails = AEI_soot + AEI_amb AEI_contrails = EI_soot_aviation !---------------------------- !-- ICE CRYSTALS LOSS -- !---------------------------- ! !--Based on Lottermoser and Unterstrasser (2025, LU25 hereinafter) !--which is an update of Unterstrasser (2016, U16 hereinafter) ! fuel consumption in kg/m flown fuel_per_m = flight_fuel / flight_dist ! LU25, Eq. A2 z_atm = 607.46 * (icesat_ratio - 1.)**0.897 * (temp / 205.)**2.225 ! water vapor emissions [kg/m flown], LU25, Eq. 2 ! U16, Eq. 6 rho_emit = fuel_per_m * EI_H2O_aviation / plume_area_loss ! LU25, Eq. A3 temp_205 = temp - 205. z_emit = 1106.6 * (rho_emit * 1e5)**(0.678 + 0.0116 * temp_205) & * EXP(- (0.0807 + 0.000428 * temp_205) * temp_205) ! U16, Eq. 4 z_desc = SQRT(8. * circ_0_loss / RPI / N_Brunt_Vaisala_aviation ) ! initial number of ice crystals [#/m flown], LU25, Eq. 1 Nice_per_m = fuel_per_m * AEI_contrails ! ice crystals number concentration [#/m3], LU25, Eq. A1 nice_init = Nice_per_m / plume_area_loss ! LU25, Eq. 9, 13d, 13e, 13f and 13g z_delta = (nice_init / nice_init_ref_loss)**(-0.16) * (1.27 * z_atm + 0.42 * z_emit) - 0.49 * z_desc ! LU25, Eq. 11, 12, 13a, 13b and 13c fsurv_contrails = 0.42 + 1.31 / RPI * ATAN(-1. + z_delta / 100.) fsurv_contrails = MIN(1., MAX(0., fsurv_contrails)) Nice_per_m_init_contrails = Nice_per_m * fsurv_contrails AEI_surv_contrails = AEI_contrails * fsurv_contrails END SUBROUTINE CONTRAIL_ICE_NUMBER_CONCENTRATION SUBROUTINE init_read_aviation_emissions USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured USE mod_phys_lmdz_para, ONLY: is_omp_master USE lmdz_xios, ONLY: xios_set_file_attr, xios_set_field_attr IF (grid_type==unstructured) THEN IF (is_omp_master) THEN ! Activate aviation file CALL xios_set_file_attr("aviation_file", enabled=.TRUE.) ! Activate aviation fields CALL xios_set_field_attr("DISTFLOWN_read", enabled=.TRUE.) CALL xios_set_field_attr("DISTFLOWN_interp", enabled=.TRUE.) CALL xios_set_field_attr("FUELBURN_read", enabled=.TRUE.) CALL xios_set_field_attr("FUELBURN_interp", enabled=.TRUE.) ENDIF ENDIF END SUBROUTINE init_read_aviation_emissions SUBROUTINE read_aviation_emissions(klon, klev) ! This subroutine allows to read the traffic density data read in the file aviation.nc ! This file is defined in ./COMP/lmdz.card ! Field names in context_input_lmdz.xml should be the same as in the file. USE netcdf USE mod_phys_lmdz_para USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, nbp_lev_=>nbp_lev, & klon_glo, grid2Dto1D_glo, grid_type, unstructured USE iophy, ONLY : io_lon, io_lat USE lmdz_xios USE print_control_mod, ONLY: lunout USE readaerosol_mod, ONLY: check_err USE lmdz_lscp_ini, ONLY: EI_H2O_aviation IMPLICIT NONE INTEGER, INTENT(IN) :: klon, klev ! number of horizontal grid points and vertical levels !---------------------------------------------------- ! Local variable !---------------------------------------------------- REAL, ALLOCATABLE :: flight_dist_mpi(:,:,:), flight_fuel_mpi(:,:,:) INTEGER :: ierr LOGICAL, SAVE :: first = .TRUE. !$OMP THREADPRIVATE(first) ! FOR REGULAR LON LAT CHARACTER(len=30) :: fname INTEGER :: nid, dimid, varid INTEGER :: imth, i, j, k REAL :: npole, spole REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_dist_glo2D REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_fuel_glo2D REAL, ALLOCATABLE, DIMENSION(:) :: vartmp_lev REAL, ALLOCATABLE, DIMENSION(:,:,:) :: vartmp REAL, ALLOCATABLE, DIMENSION(:) :: lon_src ! longitudes in file REAL, ALLOCATABLE, DIMENSION(:) :: lat_src, lat_src_inv ! latitudes in file LOGICAL :: invert_lat ! true if the field has to be inverted for latitudes INTEGER :: nbp_lon, nbp_lat IF (grid_type==unstructured) THEN IF (first) THEN IF (is_omp_master) THEN ! Get number of vertical levels and level values CALL xios_get_axis_attr("aviation_lev", n_glo=nleva) ENDIF CALL bcast_omp(nleva) ! Allocation of arrays ALLOCATE(flight_dist_read(klon, nleva, 1), STAT=ierr) IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_dist',1) ALLOCATE(flight_fuel_read(klon, nleva, 1), STAT=ierr) IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_fuel',1) ALLOCATE(aviation_lev(nleva), STAT=ierr) IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate aviation_lev',1) IF (is_omp_master) THEN ! Get number of vertical levels and level values CALL xios_get_axis_attr("aviation_lev", value=aviation_lev(:)) ENDIF CALL bcast_omp(aviation_lev) first = .FALSE. ENDIF ! first ! Read the data from the file ! is_omp_master is necessary to make XIOS works IF (is_omp_master) THEN ALLOCATE(flight_dist_mpi(klon_mpi, nleva,1), STAT=ierr) IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_dist_mpi',1) ALLOCATE(flight_fuel_mpi(klon_mpi, nleva,1), STAT=ierr) IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_fuel_mpi',1) CALL xios_recv_field("DISTFLOWN_interp", flight_dist_mpi(:,:,1)) CALL xios_recv_field("FUELBURN_interp", flight_fuel_mpi(:,:,1)) ENDIF ! Propagate to other OMP threads: flight_dist_mpi(klon_mpi,klev) to flight_dist(klon,klev) CALL scatter_omp(flight_dist_mpi, flight_dist_read) CALL scatter_omp(flight_fuel_mpi, flight_fuel_read) ELSE nbp_lon=nbp_lon_ nbp_lat=nbp_lat_ IF (is_mpi_root) THEN ALLOCATE(lon_src(nbp_lon)) ALLOCATE(lat_src(nbp_lat)) ALLOCATE(lat_src_inv(nbp_lat)) ELSE ALLOCATE(flight_dist_glo2D(0,0,0,0)) ALLOCATE(flight_fuel_glo2D(0,0,0,0)) ENDIF IF (is_mpi_root .AND. is_omp_root) THEN ! 1) Open file !**************************************************************************************** fname = 'aviation.nc' CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, nid), "pb open "//TRIM(fname) ) ! Test for equal longitudes and latitudes in file and model !**************************************************************************************** ! Read and test longitudes CALL check_err( nf90_inq_varid(nid, 'longitude', varid), "pb inq lon" ) CALL check_err( nf90_get_var(nid, varid, lon_src(:)), "pb get lon" ) IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname) WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src WRITE(lunout,*) 'longitudes in model :', io_lon CALL abort_physic('lmdz_aviation', 'longitudes are not the same in file and model',1) END IF ! Read and test latitudes CALL check_err( nf90_inq_varid(nid, 'latitude', varid),"pb inq lat" ) CALL check_err( nf90_get_var(nid, varid, lat_src(:)),"pb get lat" ) ! Invert source latitudes DO j = 1, nbp_lat lat_src_inv(j) = lat_src(nbp_lat+1-j) END DO IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN ! Latitudes are the same invert_lat=.FALSE. ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN ! Inverted source latitudes correspond to model latitudes WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname) invert_lat=.TRUE. ELSE WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname) WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src WRITE(lunout,*) 'latitudes in model :', io_lat CALL abort_physic('lmdz_aviation', 'latitudes do not correspond between file and model',1) END IF ! 2) Find vertical dimension nleva !**************************************************************************************** ierr = nf90_inq_dimid(nid, 'pressure_Pa', dimid) CALL check_err( nf90_inquire_dimension(nid, dimid, len = nleva), "pb inq dim for PRESNIVS or lev" ) ! Allocate variables depending on the number of vertical levels ALLOCATE(flight_dist_glo2D(nbp_lon, nbp_lat, nleva, 1), flight_fuel_glo2D(nbp_lon, nbp_lat, nleva, 1), stat=ierr) IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 1',1) ALLOCATE(aviation_lev(nleva), stat=ierr) IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 2',1) ! 3) Read all variables from file !**************************************************************************************** ! Get variable id CALL check_err( nf90_inq_varid(nid, 'seg_length', varid),"pb inq var seg_length_m" ) ! Get the variable CALL check_err( nf90_get_var(nid, varid, flight_dist_glo2D),"pb get var seg_length_m" ) ! Get variable id CALL check_err( nf90_inq_varid(nid, 'fuel_burn', varid),"pb inq var fuel_burn" ) ! Get the variable CALL check_err( nf90_get_var(nid, varid, flight_fuel_glo2D),"pb get var fuel_burn" ) ! Get variable id CALL check_err( nf90_inq_varid(nid, "pressure_Pa", varid),"pb inq var pressure_Pa" ) ! Get the variable CALL check_err( nf90_get_var(nid, varid, aviation_lev),"pb get var pressure_Pa" ) ! 4) Close file !**************************************************************************************** CALL check_err( nf90_close(nid), "pb in close" ) ! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo) !**************************************************************************************** ! Test if vertical levels have to be inversed IF (aviation_lev(1) < aviation_lev(nleva)) THEN ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva), vartmp_lev(nleva)) ! Inverse vertical levels for flight_dist_glo2D vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1) DO k=1, nleva DO j=1, nbp_lat DO i=1, nbp_lon flight_dist_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k) END DO END DO END DO ! Inverse vertical levels for flight_fuel_glo2D vartmp(:,:,:) = flight_fuel_glo2D(:,:,:,1) DO k=1, nleva DO j=1, nbp_lat DO i=1, nbp_lon flight_fuel_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k) END DO END DO END DO ALLOCATE(vartmp_lev(nleva)) ! Inverte vertical axes for aviation_lev vartmp_lev(:) = aviation_lev(:) DO k=1, nleva aviation_lev(k) = vartmp_lev(nleva+1-k) END DO DEALLOCATE(vartmp, vartmp_lev) ELSE WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done' END IF ! - Invert latitudes if necessary IF (invert_lat) THEN ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva)) ! Invert latitudes for the variable vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1) ! use varmth temporarly DO k=1,nleva DO j=1,nbp_lat DO i=1,nbp_lon flight_dist_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k) END DO END DO END DO ! Invert latitudes for the variable vartmp(:,:,:) = flight_fuel_glo2D(:,:,:,1) ! use varmth temporarly DO k=1,nleva DO j=1,nbp_lat DO i=1,nbp_lon flight_fuel_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k) END DO END DO END DO DEALLOCATE(vartmp) END IF ! invert_lat ! Do zonal mead at poles and distribut at whole first and last latitude DO k=1, nleva npole=0. ! North pole, j=1 spole=0. ! South pole, j=nbp_lat DO i=1,nbp_lon npole = npole + flight_dist_glo2D(i,1,k,1) spole = spole + flight_dist_glo2D(i,nbp_lat,k,1) END DO npole = npole/REAL(nbp_lon) spole = spole/REAL(nbp_lon) flight_dist_glo2D(:,1, k,1) = npole flight_dist_glo2D(:,nbp_lat,k,1) = spole END DO ! Do zonal mead at poles and distribut at whole first and last latitude DO k=1, nleva npole=0. ! North pole, j=1 spole=0. ! South pole, j=nbp_lat DO i=1,nbp_lon npole = npole + flight_fuel_glo2D(i,1,k,1) spole = spole + flight_fuel_glo2D(i,nbp_lat,k,1) END DO npole = npole/REAL(nbp_lon) spole = spole/REAL(nbp_lon) flight_fuel_glo2D(:,1, k,1) = npole flight_fuel_glo2D(:,nbp_lat,k,1) = spole END DO ALLOCATE(flight_dist_mpi(klon_glo, nleva, 1), flight_fuel_mpi(klon_glo, nleva, 1), stat=ierr) IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 3',1) ! Transform from 2D to 1D field CALL grid2Dto1D_glo(flight_dist_glo2D, flight_dist_mpi) CALL grid2Dto1D_glo(flight_fuel_glo2D, flight_fuel_mpi) ELSE ALLOCATE(flight_dist_mpi(0,0,0), flight_fuel_mpi(0,0,0)) END IF ! is_mpi_root .AND. is_omp_root !$OMP BARRIER ! 6) Distribute to all processes ! Scatter global field(klon_glo) to local process domain(klon) ! and distribute nleva to all processes !**************************************************************************************** ! Distribute nleva CALL bcast(nleva) ! Allocate and distribute pt_ap and pt_b IF (.NOT. ALLOCATED(aviation_lev)) THEN ! if pt_ap is allocated also pt_b is allocated ALLOCATE(aviation_lev(nleva), stat=ierr) IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 4',1) END IF CALL bcast(aviation_lev) ! Allocate space for output pointer variable at local process ALLOCATE(flight_dist_read(klon, nleva, 1), flight_fuel_read(klon, nleva, 1), stat=ierr) IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 5',1) ! Scatter global field to local domain at local process CALL scatter(flight_dist_mpi, flight_dist_read) CALL scatter(flight_fuel_mpi, flight_fuel_read) ENDIF END SUBROUTINE read_aviation_emissions SUBROUTINE vertical_interpolation_aviation(klon, klev, paprs, pplay, temp, flight_dist, flight_fuel) ! This subroutine performs the vertical interpolation from the read data in aviation.nc ! where there are nleva vertical levels described in aviation_lev to the klev levels or ! the model. ! flight_dist_read(klon,nleva) -> flight_dist(klon, klev) ! flight_fuel_read(klon,nleva) -> flight_fuel(klon, klev) USE lmdz_lscp_ini, ONLY: RD, RG USE lmdz_lscp_ini, ONLY: aviation_coef IMPLICIT NONE INTEGER, INTENT(IN) :: klon, klev ! number of horizontal grid points and vertical levels REAL, INTENT(IN) :: paprs(klon, klev+1) ! inter-layer pressure [Pa] REAL, INTENT(IN) :: pplay(klon, klev) ! mid-layer pressure [Pa] REAL, INTENT(IN) :: temp(klon, klev) ! temperature [K] REAL, INTENT(OUT) :: flight_dist(klon,klev) ! Aviation distance flown concentration [m/s/m3] REAL, INTENT(OUT) :: flight_fuel(klon,klev) ! Aviation fuel burn concentration [kg/s/m3] !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Local variable !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! REAL :: aviation_interface(1:nleva+1) ! Pressure of aviation file interfaces [ Pa ] INTEGER :: k, kori ! Loop index for vertical layers INTEGER :: i ! Loop index for horizontal grid REAL :: zfrac ! Fraction of layer kori in layer k REAL :: width_read_layer(1:nleva) ! width of a given layer [ Pa ] REAL :: rho, rhodz, dz ! Initialisation flight_dist(:,:) = 0. flight_fuel(:,:) = 0. ! Compute the array with the vertical interface ! It starts at 1 and has length nleva + 1 ! Note that aviation_lev has nleva and gives the altitude in the middle of the layers ! Surface pressure in standard atmosphere model [ Pa ] aviation_interface(1) = 101325. DO kori=2, nleva aviation_interface(kori) = (aviation_lev(kori-1)+aviation_lev(kori))/2.0 ! [ Pa ] ENDDO ! Last interface - we assume the same spacing as the very last one aviation_interface(nleva+1) = aviation_interface(nleva) - (aviation_lev(nleva-1) - aviation_lev(nleva)) ! Vertical width of each layer of the read file ! It is positive DO kori=1, nleva width_read_layer(kori) = aviation_interface(kori) - aviation_interface(kori+1) ENDDO ! Vertical reprojection ! The loop over klon is induced since it is done by MPI threads ! zfrac is the fraction of layer kori (read file) included in layer k (model) DO i=1,klon DO k=1, klev DO kori=1,nleva ! Which of the lower interfaces is the highest (<=> the lowest pressure) ? zfrac = min(paprs(i,k), aviation_interface(kori)) ! Which of the upper interfaces is the lowest (<=> the greatest pressure) ? zfrac = zfrac - max(paprs(i,k+1), aviation_interface(kori+1)) ! If zfrac is negative, the layers are not overlapping ! Otherwise, we get the fraction of layer kori that overlap with layer k ! after normalisation to the total kori layer width zfrac = max(0.0, zfrac) / width_read_layer(kori) ! Vertical reprojection for each desired array flight_dist(i,k) = flight_dist(i,k) + zfrac * flight_dist_read(i,kori,1) flight_fuel(i,k) = flight_fuel(i,k) + zfrac * flight_fuel_read(i,kori,1) ENDDO !--Dry density [kg/m3] rho = pplay(i,k) / temp(i,k) / RD !--Dry air mass [kg/m2] rhodz = ( paprs(i,k) - paprs(i,k+1) ) / RG !--Cell thickness [m] dz = rhodz / rho !--Normalisation with the cell thickness flight_dist(i,k) = flight_dist(i,k) / dz flight_fuel(i,k) = flight_fuel(i,k) / dz !--Enhancement factor flight_dist(i,k) = flight_dist(i,k) * aviation_coef flight_fuel(i,k) = flight_fuel(i,k) * aviation_coef ENDDO ENDDO END SUBROUTINE vertical_interpolation_aviation END MODULE lmdz_aviation