!> @file chem_emissions_mod.f90 !--------------------------------------------------------------------------------! ! This file is part of PALM model system. ! ! PALM is free software: you can redistribute it and/or modify it under the ! terms of the GNU General Public License as published by the Free Software ! Foundation, either version 3 of the License, or (at your option) any later ! version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 2018-2018 Leibniz Universitaet Hannover ! Copyright 2018-2018 Freie Universitaet Berlin ! Copyright 2018-2018 Karlsruhe Institute of Technology !--------------------------------------------------------------------------------! ! ! Current revisions: ! ------------------ ! ! ! Former revisions: ! ----------------- ! $Id: chem_emissions_mod.f90 3483 2018-11-02 14:19:26Z gronemeier $ ! bugfix: wrong locations of netCDF directives fixed ! ! 3467 2018-10-30 19:05:21Z suehring ! Enabled PARAMETRIZED mode for default surfaces when LSM is not applied but ! salsa is used ! ! 3458 2018-10-30 14:51:23Z kanani ! from chemistry branch r3443, banzhafs, Russo: ! Additional correction for index of input file of pre-processed mode ! Removed atomic and molecular weights as now available in chem_modules and ! added par_emis_time_factor (formerly in netcdf_data_input_mod) ! Added correction for index of input file of pre-processed mode ! Added a condition for default mode necessary for information read from ncdf file ! in pre-processed and default mode ! Correction of multiplying factor necessary for scaling emission values in time ! Introduced correction for scaling units in the case of DEFAULT emission mode ! ! 3373 2018-10-18 15:25:56Z kanani ! Fix wrong location of __netcdf directive ! ! 3337 2018-10-12 15:17:09Z kanani ! (from branch resler) ! Formatting ! ! 3312 2018-10-06 14:15:46Z knoop ! Initial revision ! ! 3286 2018-09-28 07:41:39Z forkel ! ! Authors: ! -------- ! @author Emmanuele Russo (Fu-Berlin) ! @author Sabine Banzhaf (FU-Berlin) ! @author Martijn Schaap (FU-Berlin, TNO Utrecht) ! ! Description: ! ------------ !> MODULE for reading-in Chemistry Emissions data !> !> @todo Check_parameters routine should be developed: for now it includes just one condition !> @todo Use of Restart files not contempled at the moment !> @todo revise indices of files read from the netcdf: preproc_emission_data and expert_emission_data !> @todo for now emission data may be passed on a singular vertical level: need to be more flexible !> @note !> @bug !------------------------------------------------------------------------------! MODULE chem_emissions_mod USE arrays_3d, & ONLY: rho_air USE control_parameters, & ONLY: initializing_actions, end_time, message_string, & intermediate_timestep_count, dt_3d USE indices USE kinds #if defined ( __netcdf ) USE NETCDF #endif USE netcdf_data_input_mod, & ONLY: chem_emis_att_type, chem_emis_val_type USE date_and_time_mod, & ONLY: time_default_indices, time_preprocessed_indices, & month_of_year, day_of_month, hour_of_day, & index_mm, index_dd, index_hh USE chem_gasphase_mod, & ONLY: spc_names USE chem_modules USE statistics, & ONLY: weight_pres IMPLICIT NONE !-- Declare all global variables within the module CHARACTER (LEN=80) :: filename_emis !> Variable for the name of the netcdf input file INTEGER(iwp) :: i !> index 1st selected dimension (some dims are not spatial) INTEGER(iwp) :: j !> index 2nd selected dimension INTEGER(iwp) :: i_start !> Index to start read variable from netcdf along one dims INTEGER(iwp) :: i_end !> Index to end read variable from netcdf in one dims INTEGER(iwp) :: j_start !> Index to start read variable from netcdf in additional dims INTEGER(iwp) :: j_end !> Index to end read variable from netcdf in additional dims INTEGER(iwp) :: z_start !> Index to start read variable from netcdf in additional dims INTEGER(iwp) :: z_end !> Index to end read variable from netcdf in additional dims INTEGER(iwp) :: dt_emis !> Time Step Emissions INTEGER(iwp) :: len_index !> length of index (used for several indices) INTEGER(iwp) :: len_index_voc !> length of voc index INTEGER(iwp) :: len_index_pm !> length of PMs index REAL(wp) :: con_factor !> Units Conversion Factor ! --------------------------------------------------------------- ! Set Values of molecules, mols, etc ! --------------------------------------------------------------- !> Avogadro number REAL, PARAMETER :: Avog = 6.02205e23 ! mlc/mol !> Dobson units: REAL, PARAMETER :: Dobs = 2.68668e16 ! (mlc/cm2) / DU !> sesalt composition: ! (Seinfeld and Pandis, "Atmospheric Chemistry and Physics", ! table 7.8 "Composition of Sea-Salt", p. 444) REAL, PARAMETER :: massfrac_Cl_in_seasalt = 0.5504 ! (kg Cl )/(kg seasalt) REAL, PARAMETER :: massfrac_Na_in_seasalt = 0.3061 ! (kg Na )/(kg seasalt) REAL, PARAMETER :: massfrac_SO4_in_seasalt = 0.0768 ! (kg SO4)/(kg seasalt) !> other numbers REAL, PARAMETER :: xm_seasalt = 74.947e-3 ! kg/mol : NaCl, SO4, .. REAL, PARAMETER :: rho_seasalt = 2.2e3 ! kg/m3 !> * amonium sulphate REAL, PARAMETER :: xm_NH4HSO4 = xm_NH4 + xm_H + xm_SO4 ! kg/mol REAL, PARAMETER :: rho_NH4HSO4a = 1.8e3 ! kg/m3 ! --------------------------------------------------------------- ! gas ! --------------------------------------------------------------- !> gas constant REAL, PARAMETER :: Rgas = 8.3144 ! J/mol/K !> gas constant for dry air REAL, PARAMETER :: Rgas_air = Rgas / xm_air ! J/kg/K !> water vapour REAL, PARAMETER :: Rgas_h2o = Rgas / xm_h2o ! J/kg/K !> standard pressure REAL, PARAMETER :: p0 = 1.0e5 ! Pa !> sea-level pressure REAL, PARAMETER :: p_sealevel = 1.01325e5 ! Pa <-- suggestion Bram Bregman !> global mean pressure REAL, PARAMETER :: p_global = 98500.0 ! Pa !> specific heat of dry air at constant pressure REAL, PARAMETER :: cp_air = 1004.0 ! J/kg/K !> Latent heat of evaporation REAL, PARAMETER :: lvap = 2.5e6 ! [J kg-1] !> Latent heat of condensation at 0 deg Celcius ! (heat (J) necesarry to evaporate 1 kg of water) REAL, PARAMETER :: Lc = 22.6e5 ! J/kg !> kappa = R/cp = 2/7 REAL, PARAMETER :: kappa = 2.0/7.0 !> Von Karman constant (dry_dep) REAL, PARAMETER :: vkarman = 0.4 !> Boltzmann constant: REAL(wp), PARAMETER :: kbolz = 1.38066e-23_wp ! J/K !> Inverse Reference Pressure (1/Pa) REAL(wp), PARAMETER :: pref_i = 1.0_wp / 100000.0_wp !> REAL(wp), PARAMETER :: r_cp = 0.286_wp !< R / cp (exponent for potential temperature) SAVE !-- Interfaces Part !-- Checks Input parameters INTERFACE chem_emissions_check_parameters MODULE PROCEDURE chem_emissions_check_parameters END INTERFACE chem_emissions_check_parameters !-- Reading of NAMELIST parameters ! INTERFACE chem_emissions_parin ! MODULE PROCEDURE chem_emissions_parin ! END INTERFACE chem_emissions_parin !-- Output of information to the header file ! INTERFACE chem_emissions_header ! MODULE PROCEDURE chem_emissions_header ! END INTERFACE chem_emissions_header !-- Matching Emissions actions INTERFACE chem_emissions_match MODULE PROCEDURE chem_emissions_match END INTERFACE chem_emissions_match !-- Initialization actions INTERFACE chem_emissions_init MODULE PROCEDURE chem_emissions_init END INTERFACE chem_emissions_init !-- Setup of Emissions INTERFACE chem_emissions_setup MODULE PROCEDURE chem_emissions_setup END INTERFACE chem_emissions_setup !-- Public Interfaces PUBLIC chem_emissions_init,chem_emissions_match,chem_emissions_setup !-- Public Variables PUBLIC con_factor, len_index,len_index_voc,len_index_pm CONTAINS !------------------------------------------------------------------------------! ! Description: ! ------------ !> Routine for checking input parameters !------------------------------------------------------------------------------! SUBROUTINE chem_emissions_check_parameters !TBD: Where Should we put the call to chem_emissions_check_parameters? In chem_init or in check_parameters? IMPLICIT NONE INTEGER(iwp) :: tmp TYPE(chem_emis_att_type) :: emt ! !-- Call routine for controlling chem_emissions namelist ! CALL chem_emissions_parin !TBD: In case the given emission values do not coincide with the passed names we should think of a solution: ! I would personally do that if the name of the species differ from the number of emission values, I would ! print a warning that says that in correspondance of that particular species the value is zero. ! An alternative would be to initialize the cs_surface_flux values to a negative number. !-- Check Emission Species Number equal to number of passed names for the chemistry species: IF ( SIZE(emt%species_name) /= emt%nspec ) THEN message_string = 'Numbers of input emission species names and number of species' // & 'for which emission values are given do not match' CALL message( 'chem_emissions_check_parameters', 'CM0437', 2, 2, 0, 6, 0 ) ENDIF !-- Check Emission Species Number equals to number of passed names for the species !IF ( SIZE(emt%species_name) /= SIZE(emt%species_index) ) THEN ! message_string = 'Number of input emission species names and ' // & ! ' number of passed species indices do not match' ! CALL message( 'chem_emissions_check_parameters', 'CM0101', 2, 2, 0, 6, 0 ) !ENDIF !-- Check Emission Categories ! IF ( SIZE(chem_emis%cat_name) /= SIZE(chem_emis%cat_index) ) THEN ! WRITE( message_string, * ) & ! 'Number of input emission categories name and categories index does not match\\' ! CALL message( 'chem_emissions_check_parameters', 'CM0101', 1, 2, 0, 6, 0 ) ! ENDIF !TBD: Check which other consistency controls should be included !TBD: Include check for spatial dimension of netcdf file. If they do not match with the ones ! of the simulation, the model should print an error. END SUBROUTINE chem_emissions_check_parameters !------------------------------------------------------------------------------! ! Description: ! ------------ !> Matching the chemical species indices. The routine checks what are the indices of the emission input species ! and the corresponding ones of the model species. The routine gives as output a vector containing the number ! of common species: it is important to note that while the model species are distinct, their values could be ! given to a single species in input: for example, in the case of NO2 and NO, values may be passed in input as NOx values. !------------------------------------------------------------------------------! SUBROUTINE chem_emissions_match(emt_att,len_index) INTEGER(iwp), INTENT(INOUT) :: len_index !< Variable where to store the number of common species between the input dataset and the model species TYPE(chem_emis_att_type), INTENT(INOUT) :: emt_att !< Chemistry Emission Array containing information for all the input chemical emission species INTEGER(iwp) :: ind_mod, ind_inp !< Parameters for cycling through chemical model and input species INTEGER(iwp) :: nspec_emis_inp !< Variable where to store the number of the emission species in input INTEGER(iwp) :: ind_voc !< Indices to check whether a split for voc should be done INTEGER(iwp) :: ispec !> index for cycle over effective number of emission species !> Tell the user where we are!! CALL location_message( 'Matching input emissions and model chemistry species', .FALSE. ) !> Number of input emission species. nspec_emis_inp=emt_att%nspec !> Check the emission mode: DEFAULT, PRE-PROCESSED or PARAMETERIZED !TBD: Add option for capital or small letters SELECT CASE(TRIM(mode_emis)) !> PRE-PROCESSED case: in this case the input species have to coincide with the ones of the model (except VOCs for which we have the VOC split): NO and NO2 in input and not NOx. CASE ("PRE-PROCESSED") CALL location_message( 'chem_emissions PRE-PROCESSED_mode_matching', .FALSE. ) len_index=0 len_index_voc=0 !> The first condition is that both the number of model and input emissions species are not null IF ( (SIZE(spc_names) .GT. 0) .AND. (nspec_emis_inp .GT. 0)) THEN !> Cycle over model species DO ind_mod = 1, SIZE(spc_names) !> Cycle over Input species DO ind_inp = 1, nspec_emis_inp !> In the PRE-PROCESSED mode each emission species is given input values, made exception for the VOCs, having the total number of VOCs in input: the model then executes a split through the different VOC species ! > Check for VOC Species: In input in this case we only have one species (VOC) IF (TRIM(emt_att%species_name(ind_inp))=="VOC") THEN !> Cycle over the input voc species: they have to be one of the VOCs of the mechanism used. A list of VOC species for different mechanisms is provided in the module documentation DO ind_voc= 1,emt_att%nvoc !> Determine the indices of the coinciding species in the two cases and assign them to matching arrays IF (TRIM(emt_att%voc_name(ind_voc))==TRIM(spc_names(ind_mod))) THEN len_index=len_index+1 len_index_voc=len_index_voc+1 ENDIF END DO ENDIF !> Other Species IF (TRIM(emt_att%species_name(ind_inp))==TRIM(spc_names(ind_mod))) THEN len_index=len_index+1 ENDIF ENDDO ENDDO !> Allocate array for storing the indices of the matched species IF (len_index>0) THEN ALLOCATE(match_spec_input(len_index)) ALLOCATE(match_spec_model(len_index)) IF (len_index_voc>0) THEN ALLOCATE(match_spec_voc_model(len_index_voc)) !> Contains indices of the VOC model species ALLOCATE(match_spec_voc_input(len_index_voc)) !> In input there is only one value for VOCs in the DEFAULT mode. This array contains the indices of the different values of VOC compositions of the input variable VOC_composition ENDIF !> Repeat the same cycle of above but passing the species indices to the newly declared arrays len_index=0 !> Cycle over model species DO ind_mod = 1, SIZE(spc_names) !> Cycle over Input species DO ind_inp = 1, nspec_emis_inp !> Determine the indices of the coinciding species in the two cases and assign them to matching arrays ! > VOCs IF ( TRIM(emt_att%species_name(ind_inp))=="VOC" .AND. ALLOCATED(match_spec_voc_input) ) THEN DO ind_voc= 1,emt_att%nvoc IF (TRIM(emt_att%voc_name(ind_voc))==TRIM(spc_names(ind_mod))) THEN len_index=len_index+1 len_index_voc=len_index_voc+1 match_spec_input(len_index)=ind_inp match_spec_model(len_index)=ind_mod match_spec_voc_input(len_index_voc)=ind_voc match_spec_voc_model(len_index_voc)=ind_mod ENDIF END DO ENDIF !>Other Species IF (TRIM(emt_att%species_name(ind_inp))==TRIM(spc_names(ind_mod))) THEN len_index=len_index+1 match_spec_input(len_index)=ind_inp match_spec_model(len_index)=ind_mod ENDIF END DO END DO !> In the case there are no species matching, the emission module should not be called ELSE message_string = 'Given Chemistry Emission Species' // & ' DO NOT MATCH' // & ' model chemical species:' // & ' Chemistry Emissions routine is not called' CALL message( 'chem_emissions_matching', 'CM0438', 0, 0, 0, 6, 0 ) ENDIF !> IF (len_index>0) ELSE !In this case either spc_names is null or nspec_emis_inp is not allocated message_string = 'Array of Emission species not allocated:' // & ' Either no emission species are provided as input or' // & ' no chemical species are used by PALM:' // & ' Chemistry Emissions routine is not called' CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 ) ENDIF !> IF ( (SIZE(spc_names) .GT. 0) .AND. ALLOCATED(nspec_emis_inp)) !> ------------------------------------------------------------------ !> DEFAULT case CASE ("DEFAULT") CALL location_message( 'chem_emissions DEFAULT_mode_matching', .FALSE. ) len_index=0 !> index for TOTAL number of species len_index_voc=0 !> index for TOTAL number of VOCs len_index_pm=3 !> index for TOTAL number of PM Types: PM1, PM2.5, PM10. It is needed because the number of emission inputs and the one of PMs in the model may not be the same. !> The first condition is that both the number of model and input emissions species are not null IF ( (SIZE(spc_names) .GT. 0) .AND. (nspec_emis_inp .GT. 0) ) THEN !> Cycle over model species DO ind_mod = 1, SIZE(spc_names) !> Cycle over input species DO ind_inp = 1, nspec_emis_inp ! > Check for VOC Species: In input in this case we only have one species (VOC) IF (TRIM(emt_att%species_name(ind_inp))=="VOC") THEN !> Cycle over the voc species given in input: they have to be one of the VOCs of the mechanism used. A list of VOC species for different mechanisms is provided in the module documentation DO ind_voc= 1,emt_att%nvoc !> Determine the indices of the coinciding species in the two cases and assign them to matching arrays IF (TRIM(emt_att%voc_name(ind_voc))==TRIM(spc_names(ind_mod))) THEN len_index=len_index+1 len_index_voc=len_index_voc+1 ENDIF END DO ENDIF !> PMs: For PMs there is only one input species name for all the PM types. This variable has 3 dimensions, one for each PM type. IF (TRIM(emt_att%species_name(ind_inp))=="PM") THEN !>pm1 IF (TRIM(spc_names(ind_mod))=="PM1") THEN len_index=len_index+1 !>pm2.5 ELSE IF (TRIM(spc_names(ind_mod))=="PM25") THEN len_index=len_index+1 !>pm10 ELSE IF (TRIM(spc_names(ind_mod))=="PM10") THEN len_index=len_index+1 ENDIF ENDIF !> NOx: for NOx by definition we have 2 splits. The Emission Input Name in this case is only one: NOx, while in the model we can have NO2 and NO IF (TRIM(emt_att%species_name(ind_inp))=="NOx") THEN !>no IF (TRIM(spc_names(ind_mod))=="NO") THEN len_index=len_index+1 !>no2 ELSE IF (TRIM(spc_names(ind_mod))=="NO2") THEN len_index=len_index+1 ENDIF ENDIF !>SOX: same as for NOx, but with SO2 and SO4 IF (TRIM(emt_att%species_name(ind_inp))=="SOx") THEN !>no IF (TRIM(spc_names(ind_mod))=="SO2") THEN len_index=len_index+1 !>no2 ELSE IF (TRIM(spc_names(ind_mod))=="SO4") THEN len_index=len_index+1 ENDIF ENDIF !> Other Species IF (TRIM(emt_att%species_name(ind_inp))==TRIM(spc_names(ind_mod))) THEN len_index=len_index+1 ENDIF END DO !> number of emission input species END DO !> number of model species !> Allocate Arrays IF (len_index>0) THEN ALLOCATE(match_spec_input(len_index)) ALLOCATE(match_spec_model(len_index)) IF (len_index_voc>0) THEN ALLOCATE(match_spec_voc_model(len_index_voc)) !> contains indices of the VOC model species ALLOCATE(match_spec_voc_input(len_index_voc)) !> In input there is only one value for VOCs in the DEFAULT mode. ! This array contains the indices of the different values of VOC compositions of the input variable VOC_composition ENDIF !> ASSIGN INDICES !> Repeat the same cycles of above, but passing the species indices to the newly declared arrays len_index=0 len_index_voc=0 DO ind_mod = 1, SIZE(spc_names) DO ind_inp = 1, nspec_emis_inp ! > VOCs IF ( TRIM(emt_att%species_name(ind_inp))=="VOC" .AND. ALLOCATED(match_spec_voc_input) ) THEN DO ind_voc= 1,emt_att%nvoc IF (TRIM(emt_att%voc_name(ind_voc))==TRIM(spc_names(ind_mod))) THEN len_index=len_index+1 len_index_voc=len_index_voc+1 match_spec_input(len_index)=ind_inp match_spec_model(len_index)=ind_mod match_spec_voc_input(len_index_voc)=ind_voc match_spec_voc_model(len_index_voc)=ind_mod ENDIF END DO ENDIF !> PMs:In this case the Inputs have only PM while the model has three different possible types: PM1, PM2.5, PM10. We need an additional index for matching each PM index in the model. IF (TRIM(emt_att%species_name(ind_inp))=="PM") THEN !>pm1 IF (TRIM(spc_names(ind_mod))=="PM1") THEN len_index=len_index+1 match_spec_input(len_index)=ind_inp match_spec_model(len_index)=ind_mod !match_spec_pm(1)=ind_mod !>pm2.5 ELSE IF (TRIM(spc_names(ind_mod))=="PM25") THEN len_index=len_index+1 match_spec_input(len_index)=ind_inp match_spec_model(len_index)=ind_mod !match_spec_pm(2)=ind_mod !>pm10 ELSE IF (TRIM(spc_names(ind_mod))=="PM10") THEN len_index=len_index+1 match_spec_input(len_index)=ind_inp match_spec_model(len_index)=ind_mod !match_spec_pm(3)=ind_mod ENDIF ENDIF !> NOx - The same as for PMs but here the model species are only 2: NO and NO2 IF (TRIM(emt_att%species_name(ind_inp))=="NOx") THEN !>no IF (TRIM(spc_names(ind_mod))=="NO") THEN len_index=len_index+1 match_spec_input(len_index)=ind_inp match_spec_model(len_index)=ind_mod !match_spec_nox(1)=ind_mod !>no2 ELSE IF (TRIM(spc_names(ind_mod))=="NO2") THEN len_index=len_index+1 match_spec_input(len_index)=ind_inp match_spec_model(len_index)=ind_mod ! match_spec_nox(2)=ind_mod ENDIF ENDIF !> SOx IF (TRIM(emt_att%species_name(ind_inp))=="SOx") THEN !>so2 IF (TRIM(spc_names(ind_mod))=="SO2") THEN len_index=len_index+1 match_spec_input(len_index)=ind_inp match_spec_model(len_index)=ind_mod ! match_spec_sox(1)=ind_mod !>so4 ELSE IF (TRIM(spc_names(ind_mod))=="SO4") THEN len_index=len_index+1 match_spec_input(len_index)=ind_inp match_spec_model(len_index)=ind_mod ! match_spec_sox(2)=ind_mod ENDIF ENDIF !> Other Species IF (TRIM(emt_att%species_name(ind_inp))==TRIM(spc_names(ind_mod))) THEN len_index=len_index+1 match_spec_input(len_index)=ind_inp match_spec_model(len_index)=ind_mod ENDIF END DO END DO ELSE message_string = 'Given Chemistry Emission Species' // & ' DO NOT MATCH' // & ' model chemical species' // & ' Chemistry Emissions routine is not called' CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 ) ENDIF ELSE message_string = 'Array of Emission species not allocated: ' // & ' Either no emission species are provided as input or' // & ' no chemical species are used by PALM:' // & ' Chemistry Emissions routine is not called' CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 ) ENDIF !-- CASE parameterized: In the parameterized case the user gives in input the exact species names of the chemical mechanism: no split is required for NOx, SOx, PMs and VOCs, but units of emissions inputs must be in (mole/m**)/s, made exception for PMs. CASE ("PARAMETERIZED") len_index=0 !spc_names have to be greater than zero for proceeding IF ( (SIZE(spc_names) .GT. 0) .AND. (nspec_emis_inp .GT. 0) ) THEN !cycle over model species DO ind_mod = 1, SIZE(spc_names) ind_inp=1 DO WHILE ( TRIM( surface_csflux_name( ind_inp ) ) /= 'novalue' ) !<'novalue' is the default IF (TRIM(surface_csflux_name( ind_inp ))==TRIM(spc_names(ind_mod))) THEN len_index=len_index+1 ENDIF ind_inp=ind_inp+1 ENDDO ENDDO !Proceed only if there are matched species IF ( len_index .GT. 0 ) THEN !Allocation of Arrays of the matched species ALLOCATE(match_spec_input(len_index)) ALLOCATE(match_spec_model(len_index)) !Assign corresponding indices of input and model matched species len_index=0 DO ind_mod = 1, SIZE(spc_names) ind_inp = 1 DO WHILE ( TRIM( surface_csflux_name( ind_inp ) ) /= 'novalue' ) !<'novalue' is the default IF (TRIM( surface_csflux_name( ind_inp ) ) == TRIM(spc_names(ind_mod))) THEN len_index=len_index+1 match_spec_input(len_index)=ind_inp match_spec_model(len_index)=ind_mod ENDIF ind_inp=ind_inp+1 END DO END DO ! also, Add AN if to check that when the surface_csflux are passed to the namelist, also the street factor values are provided DO ispec = 1 , len_index IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND. emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN message_string = 'PARAMETERIZED emissions mode selected:' // & ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES' // & ' but values of scaling factors for street types' // & ' emiss_factor_main AND emiss_factor_side' // & ' not provided for each of the emissions species' // & ' or not provided at all: PLEASE set a finite value' // & ' for these parameters in the chemistry namelist' CALL message( 'chem_emissions_matching', 'CM0442', 2, 2, 0, 6, 0 ) ENDIF END DO ELSE message_string = 'Given Chemistry Emission Species' // & ' DO NOT MATCH' // & ' model chemical species' // & ' Chemistry Emissions routine is not called' CALL message( 'chem_emissions_matching', 'CM0443', 0, 0, 0, 6, 0 ) ENDIF ELSE message_string = 'Array of Emission species not allocated: ' // & ' Either no emission species are provided as input or' // & ' no chemical species are used by PALM.' // & ' Chemistry Emissions routine is not called' CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 ) ENDIF !-- CASE when emission module is switched on but mode_emis is not specified or it is given the wrong name CASE DEFAULT message_string = 'Chemistry Emissions Module switched ON, but' // & ' either no emission mode specified or incorrectly given :' // & ' please, pass the correct value to the namelist parameter "mode_emis"' CALL message( 'chem_emissions_matching', 'CM0445', 2, 2, 0, 6, 0 ) END SELECT END SUBROUTINE chem_emissions_match !------------------------------------------------------------------------------! ! Description: ! ------------ !> Initialization: !> Netcdf reading, arrays allocation and first calculation of cssws fluxes at timestep 0 !> !------------------------------------------------------------------------------! SUBROUTINE chem_emissions_init(emt_att,emt,nspec_out) USE surface_mod, & ONLY: surf_lsm_h,surf_def_h,surf_usm_h IMPLICIT NONE TYPE(chem_emis_att_type),INTENT(INOUT) :: emt_att !> variable where to store all the information of ! emission inputs whose values do not depend ! on the considered species TYPE(chem_emis_val_type),INTENT(INOUT), ALLOCATABLE, DIMENSION(:) :: emt !> variable where to store emission inputs values, ! depending on the considered species INTEGER(iwp),INTENT(INOUT) :: nspec_out !> number of outputs of chemistry emission routines CHARACTER (LEN=80) :: units !> units of chemistry inputs INTEGER(iwp) :: ispec !> Index to go through the emission chemical species !-- Actions for initial runs : TBD: needs to be updated IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN !-- ... ! !-- Actions for restart runs ELSE !-- ... ENDIF CALL location_message( 'Starting initialization of chemistry emissions arrays', .FALSE. ) !-- Match Input and Model emissions CALL chem_emissions_match(emt_att,nspec_out) !> If nspec_out==0, then do not use emission module: The corresponding message is already printed in the matching routine IF ( nspec_out == 0 ) THEN emission_output_required=.FALSE. ELSE emission_output_required=.TRUE. !----------------------------------------------------------------- !> Set molecule masses' ALLOCATE(emt_att%xm(nspec_out)) ! loop over emisisons: DO ispec = 1, nspec_out ! switch: SELECT CASE ( TRIM(spc_names(match_spec_model(ispec))) ) CASE ( 'SO2','SO4' ) ; emt_att%xm(ispec) = xm_S + xm_O * 2 ! kg/mole CASE ( 'NO','NO2' ) ; emt_att%xm(ispec) = xm_N + xm_O * 2 ! kg/mole NO2 equivalent CASE ( 'NH3' ) ; emt_att%xm(ispec) = xm_N + xm_H * 3 ! kg/mole CASE ( 'CO' ) ; emt_att%xm(ispec) = xm_C + xm_O ! kg/mole CASE ( 'CO2' ) ; emt_att%xm(ispec) = xm_C + xm_O * 2 ! kg/mole CASE ( 'CH4' ) ; emt_att%xm(ispec) = xm_C + xm_H * 4 ! kg/mole CASE ( 'HNO3' ); emt_att%xm(ispec) = xm_H + xm_N + xm_O*3 !kg/mole CASE DEFAULT !-- TBD: check if this hase to be removed !emt_att%xm(ispec) = 1.0_wp END SELECT ENDDO !> ASSIGN EMISSION VALUES for the different emission modes SELECT CASE(TRIM(mode_emis)) !TBD: Add the option for CApital or small letters !> PRE-PROCESSED case CASE ("PRE-PROCESSED") IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(nzb:nzt+1,0:ny,0:nx,nspec_out)) CALL location_message( 'emis_distribution array allocated in PRE-PROCESSED mode', .FALSE. ) !> Calculate the values of the emissions at the first time step CALL chem_emissions_setup(emt_att,emt,nspec_out) !> Default case CASE ("DEFAULT") IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution( 1, 0:ny, 0:nx, nspec_out ) ) CALL location_message( 'emis_distribution array allocated in DEFAULT mode', .FALSE. ) !> Calculate the values of the emissions at the first time step CALL chem_emissions_setup(emt_att,emt,nspec_out) !> PARAMETERIZED case CASE ("PARAMETERIZED") CALL location_message( 'emis_distribution array allocated in PARAMETERIZED mode', .FALSE. ) ! For now for PAR and DEF values only, first vertical level of emis_distribution is allocated, while for PRE-PROCESSED all. IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(1,0:ny,0:nx,nspec_out)) !> Calculate the values of the emissions at the first time step CALL chem_emissions_setup(emt_att,emt,nspec_out) END SELECT ENDIF END SUBROUTINE chem_emissions_init !------------------------------------------------------------------------------! ! Description: ! ------------ !> Routine for Update of Emission values at each timestep !-------------------------------------------------------------------------------! SUBROUTINE chem_emissions_setup(emt_att,emt,nspec_out) USE arrays_3d, & ONLY: dzw USE grid_variables, & ONLY: dx, dy USE indices, & ONLY: nnx,nny,nnz USE salsa_mod, & ONLY: salsa USE surface_mod, & ONLY: surf_lsm_h,surf_def_h,surf_usm_h USE netcdf_data_input_mod, & ONLY: street_type_f USE arrays_3d, & ONLY: hyp, pt IMPLICIT NONE !--- IN/OUT TYPE(chem_emis_att_type),INTENT(INOUT) :: emt_att !> variable where to store all the information of ! emission inputs whose values do not depend ! on the considered species TYPE(chem_emis_val_type),INTENT(INOUT), ALLOCATABLE, DIMENSION(:) :: emt !> variable where to store emission inputs values, ! depending on the considered species INTEGER,INTENT(IN) :: nspec_out !> Output of routine chem_emis_matching with number ! of matched species !--- REAL(wp),ALLOCATABLE,DIMENSION(:,:) :: delta_emis !> Term to add to the emission_distribution array ! for each of the categories at each time step ! in the default mode CHARACTER(LEN=80) :: units !> Units of the emissions INTEGER(iwp) :: icat !> Index for number of categories INTEGER(iwp) :: ispec !> index for number of species INTEGER(iwp) :: i_pm_comp !> index for number of PM components INTEGER(iwp) :: ivoc !> Index for number of components of VOCs INTEGER(iwp) :: time_index !> Index for time REAL(wp),ALLOCATABLE, DIMENSION(:) :: time_factor !> factor for time scaling of emissions REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: emis REAL(wp), DIMENSION(24) :: par_emis_time_factor !< time factors ! for the parameterized mode: these are fixed for each hour ! of a single day. REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: conv_to_ratio !> factor used for converting input ! to adimensional concentration ratio ! ------------------------------------------ ! variables for the conversion of indices i and j according to surface_mod INTEGER(iwp) :: i !> running index for grid in x-direction INTEGER(iwp) :: j !> running index for grid in y-direction INTEGER(iwp) :: k !> running index for grid in z-direction INTEGER(iwp) :: m !> running index for horizontal surfaces REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: tmp_temp ! --- const ------------------------------- !-CONVERSION FACTORS: TIME ! number of sec per hour = 3600 REAL, PARAMETER :: s_per_hour = 3600.0 ! (s)/(hour) ! number of sec per day = 86400 REAL, PARAMETER :: s_per_day = 86400.0 ! (s)/(day) ! number of hours in a year of 365 days: REAL, PARAMETER :: hour_per_year = 8760.0 !> TBD: What about leap years? ! number of hours in a day: REAL, PARAMETER :: hour_per_day = 24.0 ! conversion from hours to seconds (s/hour) = 1/3600.0 ~ 0.2777778 REAL, PARAMETER :: hour_to_s = 1/s_per_hour ! (hour)/(s) ! conversion from day to seconds (s/day) = 1/86400.0 ~ 1.157407e-05 REAL, PARAMETER :: day_to_s = 1/s_per_day ! (day)/(s) ! conversion from year to sec (s/year) = 1/31536000.0 ~ 3.170979e-08 REAL, PARAMETER :: year_to_s = 1/(s_per_hour*hour_per_year) ! (year)/(s) !-CONVERSION FACTORS: WEIGHT ! Conversion from tons to kg (kg/tons) = 100.0/1 ~ 100 REAL, PARAMETER :: tons_to_kg = 100 ! (tons)/(kg) ! Conversion from g to kg (kg/g) = 1/1000.0 ~ 0.001 REAL, PARAMETER :: g_to_kg = 0.001 ! (g)/(kg) ! Conversion from g to kg (kg/g) = 1/1000.0 ~ 0.001 REAL, PARAMETER :: miug_to_kg = 0.000000001 ! (microng)/(kg) !< CONVERSION FACTORS: fraction to ppm REAL(wp), PARAMETER :: ratio2ppm = 1.0e06_wp !------------------------------------------------------ IF ( emission_output_required ) THEN !> Set emis_dt !TBD: this is the same as dt_chem. We should consider the fact that dt_emis should be the timestep of input emissions or better defined, the timestep at which the emission routines are called: for now one hour. It should be made changeable. IF ( call_chem_at_all_substeps ) THEN dt_emis = dt_3d * weight_pres(intermediate_timestep_count) ELSE dt_emis = dt_3d ENDIF ! --- CHECK UNITS !>----------------------------------------------------- !> Conversion of the units to the ones employed in PALM. !> Possible temporal units of input emissions data are: year, hour and second; !> the mass could be expressed in: tons, kilograms or grams. !> IN the PARAMETERIZED mode no conversion is possible: in this case INPUTS have to have fixed units. !>----------------------------------------------------- IF (TRIM(mode_emis)=="DEFAULT" .OR. TRIM(mode_emis)=="PRE-PROCESSED" ) THEN SELECT CASE(TRIM(emt_att%units)) !> kilograms CASE ( 'kg/m2/s','KG/M2/S') CALL location_message( 'Units of the emissions are consistent with' // & ' the ones employed in the PALM-4U model ', .FALSE. ) con_factor=1 CASE ('kg/m2/hour','KG/M2/HOUR') CALL location_message( 'Units of emission inputs need conversion', .FALSE. ) con_factor=hour_to_s CASE ( 'kg/m2/day','KG/M2/DAY' ) CALL location_message( 'Units of emission inputs need conversion', .FALSE. ) con_factor=day_to_s CASE ( 'kg/m2/year','KG/M2/YEAR' ) CALL location_message( 'Units of emission inputs need conversion', .FALSE. ) con_factor=year_to_s !> Tons CASE ( 'ton/m2/s','TON/M2/S' ) CALL location_message( 'Units of emission inputs need conversion', .FALSE. ) con_factor=tons_to_kg CASE ( 'ton/m2/hour','TON/M2/HOUR' ) CALL location_message( 'Units of emission inputs need conversion', .FALSE. ) con_factor=tons_to_kg*hour_to_s CASE ( 'ton/m2/year','TON/M2/YEAR' ) CALL location_message( 'Units of emission inputs need conversion', .FALSE. ) con_factor=tons_to_kg*year_to_s !> Grams CASE ( 'g/m2/s','G/M2/S' ) CALL location_message( 'Units of emission inputs need conversion', .FALSE. ) con_factor=g_to_kg CASE ( 'g/m2/hour','G/M2/HOUR' ) CALL location_message( 'Units of emission inputs need conversion', .FALSE. ) con_factor=g_to_kg*hour_to_s CASE ( 'g/m2/year','G/M2/YEAR' ) CALL location_message( 'Units of emission inputs need conversion', .FALSE. ) con_factor=g_to_kg*year_to_s !> Micro Grams CASE ( 'micrograms/m2/s','MICROGRAMS/M2/S' ) CALL location_message( 'Units of emission inputs need conversion', .FALSE. ) con_factor=miug_to_kg CASE ( 'micrograms/m2/hour','MICROGRAMS/M2/HOUR' ) CALL location_message( 'Units of emission inputs need conversion', .FALSE. ) con_factor=miug_to_kg*hour_to_s CASE ( 'micrograms/m2/year','MICROGRAMS/M2/YEAR' ) CALL location_message( 'Units of emission inputs need conversion', .FALSE. ) con_factor=miug_to_kg*year_to_s CASE DEFAULT message_string = 'The Units of the provided input chemistry emission species' // & ' are not the ones required by PALM-4U: please check ' // & ' chemistry emission module documentation.' CALL message( 'chem_emissions_setup', 'CM0446', 2, 2, 0, 6, 0 ) END SELECT !> PRE-PROCESSED and parameterized mode ELSE message_string = 'No Units conversion required for units of chemistry emissions' // & ' of the PARAMETERIZED mode: units have to be provided in' // & ' micromole/m**2/day for GASES and' // & ' kg/m**2/day for PMs' CALL message( 'chem_emissions_setup', 'CM0447', 0, 0, 0, 6, 0 ) ENDIF !> Conversion factors (if the units are kg/m**2/s) we have to convert them to ppm/s DO i=nxl,nxr DO j=nys,nyn !> Derive Temperature from Potential Temperature tmp_temp(nzb:nzt+1,j,i) = pt(nzb:nzt+1,j,i) * ( hyp(nzb:nzt+1) * pref_i )**r_cp !> We need to pass to cssws<- (ppm/s) * dz. ! Input is Nmole/(m^2*s) ! To go to ppm*dz basically we need to multiply the input by (m**2/N)*dz ! (m**2/N)*dz == V/N ! V/N=RT/P !> m**3/Nmole (J/mol)*K^-1 K Pa conv_to_ratio(nzb:nzt+1,j,i) = ( (Rgas * tmp_temp(nzb:nzt+1,j,i)) / ((hyp(nzb:nzt+1))) ) ENDDO ENDDO !>------------------------------------------------ !> Start The Processing of the input data emis_distribution(:,nys:nyn,nxl:nxr,:) = 0.0_wp !>----------------------------------------------- !> Distinguish between DEFAULT, PRE-PROCESSED and PARAMETERIZED mode when calculating time_factors: only done for DEFAULT mode. For the PARAMETERIZED mode there is a time factor but it is fixed in the model !> PRE-PROCESSED MODE IF (TRIM(mode_emis)=="PRE-PROCESSED") THEN !> Update time indices CALL time_preprocessed_indices(index_hh) CALL location_message( 'PRE-PROCESSED MODE: No time-factor specification required', .FALSE. ) ELSEIF (TRIM(mode_emis)=="DEFAULT") THEN CALL location_message( 'DEFAULT MODE: time-factor specification required', .FALSE. ) !> Allocate array where to store temporary emission values IF(.NOT. ALLOCATED(emis)) ALLOCATE(emis(nys:nyn,nxl:nxr)) !> Allocate time factor per emitted component category ALLOCATE(time_factor(emt_att%ncat)) !> Read-in HOURLY emission time factor IF (TRIM(time_fac_type)=="HOUR") THEN !> Update time indices CALL time_default_indices(month_of_year,day_of_month,hour_of_day,index_hh) !> Check if the index is less or equal to the temporal dimension of HOURLY emission files IF (index_hh .LE. SIZE(emt_att%hourly_emis_time_factor(1,:))) THEN !> Read-in the correspondant time factor time_factor(:)= emt_att%hourly_emis_time_factor(:,index_hh) ELSE message_string = 'The "HOUR" time-factors (DEFAULT mode) of the chemistry emission species' // & ' are not provided for each hour of the total simulation time' CALL message( 'chem_emissions_setup', 'CM0448', 2, 2, 0, 6, 0 ) ENDIF !> Read-in MDH emissions time factors ELSEIF (TRIM(time_fac_type)=="MDH") THEN !> Update time indices CALL time_default_indices(daytype_mdh,month_of_year,day_of_month,hour_of_day,index_mm,index_dd,index_hh) !> Check if the index is less or equal to the temporal dimension of MDH emission files IF ((index_hh+index_dd+index_mm) .LE. SIZE(emt_att%mdh_emis_time_factor(1,:))) THEN !> Read-in the correspondant time factor time_factor(:)=emt_att%mdh_emis_time_factor(:,index_mm)*emt_att%mdh_emis_time_factor(:,index_dd)* & emt_att%mdh_emis_time_factor(:,index_hh) ELSE message_string = 'The "MDH" time-factors (DEFAULT mode) of the chemistry emission species' // & ' are not provided for each hour/day/month of the total simulation time' CALL message( 'chem_emissions_setup', 'CM0449', 2, 2, 0, 6, 0 ) ENDIF ELSE !> condition when someone used the DEFAULT mode but forgets to indicate the time-factor in the namelist message_string = 'The time-factor (DEFAULT mode) of the chemistry emission species' // & ' is not provided in the NAMELIST' CALL message( 'chem_emissions_setup', 'CM0450', 2, 2, 0, 6, 0 ) ENDIF !> PARAMETERIZED MODE ELSEIF (TRIM(mode_emis)=="PARAMETERIZED") THEN CALL location_message( 'PARAMETERIZED MODE: time-factor specification is fixed: ' // & ' 24 values for every day of the year ', .FALSE. ) !Assign Constant Values of time factors, diurnal time profile for traffic sector: par_emis_time_factor( : ) = (/ 0.009, 0.004, 0.004, 0.009, 0.029, 0.039, 0.056, 0.053, 0.051, 0.051, 0.052, 0.055, & 0.059, 0.061, 0.064, 0.067, 0.069, 0.069, 0.049, 0.039, 0.039, 0.029, 0.024, 0.019 /) !> in this case allocate time factor each hour in a day IF (.NOT. ALLOCATED(time_factor)) ALLOCATE(time_factor(1)) !>Pass the values of time factors in the parameterized mode to the time_factor variable. in this case index_hh==hour_of_day index_hh=hour_of_day time_factor(1) = par_emis_time_factor(index_hh) ENDIF !-------------------------------- !-- EMISSION DISTRIBUTION Calculation !> PARAMETERIZED CASE IF ( mode_emis == "PARAMETERIZED" ) THEN DO ispec=1,nspec_out !> Values are still micromoles/(m**2*s). Units in this case are always micromoles/m**2*day (or kilograms/m**2*day for PMs) emis_distribution(1,nys:nyn,nxl:nxr,ispec)=surface_csflux(match_spec_input(ispec))*time_factor(1)*hour_to_s ENDDO !> PRE-PROCESSED CASE ELSEIF (TRIM(mode_emis)=="PRE-PROCESSED") THEN !> Start Cycle over Species DO ispec=1,nspec_out !> nspec_out represents the number of species in common between ! the emission input data and the chemistry mechanism used emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emt(match_spec_input(ispec))% & preproc_emission_data(index_hh,1,nys+1:nyn+1,nxl+1:nxr+1)* & con_factor ENDDO !TBD: At the moment the default mode considers a single vertical level (the surface). So we need to change it accordingly or eventually include the variable vertical_profile to be used in the default mode i we want to consider additional levels. ELSE IF (TRIM(mode_emis)=="DEFAULT") THEN !> Allocate array for the emission value corresponding to a specific category and time factor ALLOCATE(delta_emis(nys:nyn,nxl:nxr)) !> Start Cycle over categories DO icat = 1, emt_att%ncat !> Start Cycle over Species DO ispec=1,nspec_out !> nspec_out represents the number of species in common between ! the emission input data and the chemistry mechanism used emis(nys:nyn,nxl:nxr) = emt(match_spec_input(ispec))%default_emission_data(icat,nys+1:nyn+1,nxl+1:nxr+1) !TBD: The consideration of dt_emis of the input data is still missing. Basically the emissions could be provided every 10, 30 minutes and not always at one hour. This should be eventually solved in the date_and_time mode routine. !> the time factors are 24 for each day. When multiplied by a daily value, they allow to have an hourly value. Then to convert it to seconds, we still have to divide this value by 3600. ! So given any units, we convert them to seconds and finally multiply them by 24 ((value/sec)*(24*3600)=value/day ---- (value/day)*time_factor=value/hour ---(value/hour)/(3600)=value/sec ) ! ((value/sec)*(24*3600)*time_factor)/3600=24*(value/sec)*time_factor !> NOX Compositions IF (TRIM(spc_names(match_spec_model(ispec)))=="NO") THEN !> Kg/m2*s kg/m2*s delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,1)*con_factor*hour_per_day emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="NO2") THEN delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,2)*con_factor*hour_per_day emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) !> SOX Compositions ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO2") THEN !> Kg/m2*s kg/m2*s delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,1)*con_factor*hour_per_day emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO4") THEN !> Kg/m2*s kg/m2*s delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,2)*con_factor*hour_per_day emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) !> PMs should be in kg/m**2/s, so conversion factor is here still required !> PM1 Compositions ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1") THEN !> Cycle over the different pm components for PM1 type DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,1)) delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,1)*con_factor*hour_per_day emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) ENDDO !> PM2.5 Compositions ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="PM25") THEN !> Cycle over the different pm components for PM2.5 type DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,2)) delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,2)*con_factor*hour_per_day emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) ENDDO !> PM10 Compositions ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN !> Cycle over the different pm components for PM10 type DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,3)) delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,3)*con_factor*hour_per_day emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) ENDDO !> VOCs Compositions: for VOCs, the input value is provided in kg/m**2*s but the composition is provided in mole/kg: a conversion factor for the input that could be eventually provided in, for example, tons/(m**2*s) is still required ELSE IF (SIZE(match_spec_voc_input) .GT. 0) THEN !TBD: maybe we can avoid the cycle DO ivoc= 1,SIZE(match_spec_voc_input) IF (TRIM(spc_names(match_spec_model(ispec)))==TRIM(emt_att%voc_name(ivoc))) THEN delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)* & emt_att%voc_comp(icat,match_spec_voc_input(ivoc))*con_factor*hour_per_day emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) ENDIF ENDDO !> General case (other species) ELSE delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*con_factor*hour_per_day emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) ENDIF ! IF (spc_names==) !> for every species and category set emis to 0 so to avoid overwriting. TBD: discuss whether necessary emis(:,:)= 0 ENDDO !> for every category set delta_emis to 0 so to avoid overwriting. TBD: discuss whether necessary delta_emis(:,:)=0 ENDDO ENDIF !> mode_emis !------------------------------------------------------------------------------------------------------- !--- Cycle to transform x,y coordinates to the one of surface_mod and to assign emission values to cssws !------------------------------------------------------------------------------------------------------- !-- PARAMETERIZED mode !> For the PARAMETERIZED mode units of inputs are always micromoles/(m**2*s). In this case we do not need the molar mass for conversion into ppms IF (TRIM(mode_emis)=="PARAMETERIZED") THEN IF ( street_type_f%from_file ) THEN !> Streets are lsm surfaces, hence, no usm surface treatment required IF (surf_lsm_h%ns .GT. 0 ) THEN DO m = 1, surf_lsm_h%ns i = surf_lsm_h%i(m) j = surf_lsm_h%j(m) k = surf_lsm_h%k(m) IF ( street_type_f%var(j,i) >= main_street_id .AND. street_type_f%var(j,i) < max_street_id ) THEN !> Cycle over already matched species DO ispec=1,nspec_out !> PMs are already in mass units: kilograms IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN ! kg/(m^2*s) *kg/m^3 surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec)) * & ! kg/(m^2*s) emis_distribution(1,j,i,ispec)* & ! kg/m^3 rho_air(k) ELSE !> Other Species: inputs are micromoles: they have to be converted in moles ! ppm/s *m *kg/m^3 surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec))* & ! micromoles/(m^2*s) emis_distribution(1,j,i,ispec) * & ! m^3/Nmole conv_to_ratio(k,j,i)* & ! kg/m^3 rho_air(k) ENDIF ENDDO ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. street_type_f%var(j,i) < main_street_id ) THEN !> Cycle over already matched species DO ispec=1,nspec_out !> PMs are already in mass units: micrograms IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN ! kg/(m^2*s) *kg/m^3 surf_lsm_h%cssws(match_spec_model(ispec),m)= emiss_factor_side(match_spec_input(ispec)) * & ! kg/(m^2*s) emis_distribution(1,j,i,ispec)* & ! kg/m^3 rho_air(k) ELSE !>Other Species: inputs are micromoles ! ppm/s *m *kg/m^3 surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * & ! micromoles/(m^2*s) emis_distribution(1,j,i,ispec) * & ! m^3/Nmole conv_to_ratio(k,j,i)* & ! kg/m^3 rho_air(k) ENDIF ENDDO ELSE !> If no street type is defined, then assign null emissions to all the species surf_lsm_h%cssws(:,m) = 0.0_wp ENDIF ENDDO ELSEIF ( salsa ) THEN DO m = 1, surf_def_h(0)%ns i = surf_def_h(0)%i(m) j = surf_def_h(0)%j(m) k = surf_def_h(0)%k(m) IF ( street_type_f%var(j,i) >= main_street_id .AND. & street_type_f%var(j,i) < max_street_id ) & THEN !> Cycle over already matched species DO ispec=1,nspec_out !> PMs are already in mass units:micrograms: have to be converted to kilograms IF ( TRIM(spc_names(match_spec_model(ispec)))=="PM1" & .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10")& THEN surf_def_h(0)%cssws(match_spec_model(ispec),m) = & emiss_factor_main(match_spec_input(ispec)) * & emis_distribution(1,j,i,ispec) * rho_air(k) /& time_factor(1) ELSE !> Other Species: inputs are micromoles: have to be converted surf_def_h(0)%cssws(match_spec_model(ispec),m) = & emiss_factor_main(match_spec_input(ispec)) * & emis_distribution(1,j,i,ispec) * & conv_to_ratio(k,j,i) * rho_air(k) / time_factor(1) ENDIF ENDDO ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. & street_type_f%var(j,i) < main_street_id ) & THEN !> Cycle over already matched species DO ispec=1,nspec_out !> PMs are already in mass units: micrograms IF ( TRIM(spc_names(match_spec_model(ispec)))=="PM1" & .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10")& THEN surf_def_h(0)%cssws(match_spec_model(ispec),m) = & emiss_factor_side(match_spec_input(ispec)) * & emis_distribution(1,j,i,ispec) * rho_air(k) / & time_factor(1) ELSE surf_def_h(0)%cssws(match_spec_model(ispec),m) = & emiss_factor_side(match_spec_input(ispec)) * & emis_distribution(1,j,i,ispec) * & conv_to_ratio(k,j,i) * rho_air(k) / time_factor(1) ENDIF ENDDO ELSE !> If no street type is defined, then assign null emissions to all the species surf_def_h(0)%cssws(:,m) = 0.0_wp ENDIF ENDDO ENDIF ENDIF !> For both DEFAULT and PRE-PROCESSED ELSE DO ispec=1,nspec_out !TBD: for the PRE-PROCESSED mode in the future, values at different heights should be included! !> Default surface type IF (surf_def_h(0)%ns .GT. 0) THEN DO m = 1, surf_def_h(0)%ns i = surf_def_h(0)%i(m) j = surf_def_h(0)%j(m) IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN !> Distinguish between PMs (no needing conversion in ppms), ! VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and ! other Species (still expressed in Kg/(m**2*s) at this point) !> PMs IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN ! kg/(m^2*s) *kg/m^3 kg/(m^2*s) surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)* & ! kg/m^3 rho_air(nzb) ELSE !> VOCs IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN ! ( ppm/s) * m * kg/m^3 mole/(m^2/s) surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & ! m^3/mole ppm conv_to_ratio(nzb,j,i)*ratio2ppm * & ! kg/m^3 rho_air(nzb) !> OTHER SPECIES ELSE ! ( ppm/s )*m * kg/m^3 kg/(m^2/s) surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & ! mole/Kg (1/emt_att%xm(ispec))* & ! m^3/mole ppm conv_to_ratio(nzb,j,i)*ratio2ppm* & ! kg/m^3 rho_air(nzb) ENDIF ENDIF ENDIF ENDDO END IF !-- Land Surface Mode surface type IF (surf_lsm_h%ns .GT. 0) THEN DO m = 1, surf_lsm_h%ns i = surf_lsm_h%i(m) j = surf_lsm_h%j(m) k = surf_lsm_h%k(m) IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN !> Distinguish between PMs (no needing conversion in ppms), ! VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and ! other Species (still expressed in Kg/(m**2*s) at this point) !> PMs IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN ! kg/(m^2*s) * kg/m^3 kg/(m^2*s) surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & ! kg/m^3 rho_air(k) ELSE !> VOCs IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN ! ( ppm/s) * m * kg/m^3 mole/(m^2/s) surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & ! m^3/mole ppm conv_to_ratio(k,j,i)*ratio2ppm* & ! kg/m^3 rho_air(k) !> OTHER SPECIES ELSE ! ( ppm/s) * m * kg/m^3 kg/(m^2/s) surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & ! mole/Kg (1/emt_att%xm(ispec))* & ! m^3/mole ppm conv_to_ratio(k,j,i)*ratio2ppm* & ! kg/m^3 rho_air(k) ENDIF ENDIF ENDIF ENDDO END IF !-- Urban Surface Mode surface type IF (surf_usm_h%ns .GT. 0) THEN DO m = 1, surf_usm_h%ns i = surf_usm_h%i(m) j = surf_usm_h%j(m) k = surf_usm_h%k(m) IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN !> PMs IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN ! kg/(m^2*s) *kg/m^3 kg/(m^2*s) surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)* & ! kg/m^3 rho_air(k) ELSE !> VOCs IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN ! ( ppm/s) * m * kg/m^3 mole/(m^2/s) surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & ! m^3/mole ppm conv_to_ratio(k,j,i)*ratio2ppm* & ! kg/m^3 rho_air(k) !> OTHER SPECIES ELSE ! ( ppm/s) * m * kg/m^3 kg/(m^2/s) surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & ! mole/Kg (1/emt_att%xm(ispec))* & ! m^3/mole ppm conv_to_ratio(k,j,i)*ratio2ppm* & ! kg/m^3 rho_air(k) ENDIF ENDIF ENDIF ENDDO END IF ENDDO ENDIF !> At the end of each call to chem_emissions setup, the time_factor is deallocated (ALLOCATED ONLY in the DEFAULT mode) IF ( ALLOCATED ( time_factor ) ) DEALLOCATE( time_factor ) ENDIF !> emis_output_required END SUBROUTINE chem_emissions_setup END MODULE chem_emissions_mod