!> @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:
! ------------------
! - Removed salsa dependency.
! - Enabled PARAMETRIZED mode for default surfaces when LSM is not applied but
! salsa is (M. Kurppa)
!
! Former revisions:
! -----------------
! $Id: chem_emissions_mod.f90 3582 2018-11-29 19:16:36Z sward $
! resler:
! Break lines at 132 characters
!
! 3483 2018-11-02 14:19:26Z raasch
! 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 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
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