Home

Context Navigation

source: palm/trunk/SOURCE/chem_emissions_mod.f90 @ 3558

Last change on this file since 3558 was 3483, checked in by raasch, 6 years ago
bugfix: misplaced positions of cpp-directives for netCDF and MPI fixed; output format limited to a maximum line length of 80
Property svn:keywords set to `Id`
File size: 80.9 KB

Rev	Line
[3190]	1	!> @file chem_emissions_mod.f90
	2	!--------------------------------------------------------------------------------!
	3	! This file is part of PALM model system.
	4	!
	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[3298]	17	! Copyright 2018-2018 Leibniz Universitaet Hannover
	18	! Copyright 2018-2018 Freie Universitaet Berlin
	19	! Copyright 2018-2018 Karlsruhe Institute of Technology
[3190]	20	!--------------------------------------------------------------------------------!
	21	!
	22	! Current revisions:
[3298]	23	! ------------------
[3190]	24	!
[3298]	25	!
[3190]	26	! Former revisions:
	27	! -----------------
	28	! $Id: chem_emissions_mod.f90 3483 2018-11-02 14:19:26Z kanani $
[3483]	29	! bugfix: wrong locations of netCDF directives fixed
	30	!
	31	! 3467 2018-10-30 19:05:21Z suehring
[3467]	32	! Enabled PARAMETRIZED mode for default surfaces when LSM is not applied but
	33	! salsa is used
	34	!
	35	! 3458 2018-10-30 14:51:23Z kanani
[3458]	36	! from chemistry branch r3443, banzhafs, Russo:
	37	! Additional correction for index of input file of pre-processed mode
	38	! Removed atomic and molecular weights as now available in chem_modules and
	39	! added par_emis_time_factor (formerly in netcdf_data_input_mod)
	40	! Added correction for index of input file of pre-processed mode
	41	! Added a condition for default mode necessary for information read from ncdf file
	42	! in pre-processed and default mode
	43	! Correction of multiplying factor necessary for scaling emission values in time
	44	! Introduced correction for scaling units in the case of DEFAULT emission mode
	45	!
	46	! 3373 2018-10-18 15:25:56Z kanani
[3373]	47	! Fix wrong location of __netcdf directive
	48	!
	49	! 3337 2018-10-12 15:17:09Z kanani
[3337]	50	! (from branch resler)
	51	! Formatting
	52	!
	53	! 3312 2018-10-06 14:15:46Z knoop
[3298]	54	! Initial revision
[3190]	55	!
[3298]	56	! 3286 2018-09-28 07:41:39Z forkel
	57	!
[3190]	58	! Authors:
	59	! --------
[3458]	60	! @author Emmanuele Russo (Fu-Berlin)
	61	! @author Sabine Banzhaf (FU-Berlin)
	62	! @author Martijn Schaap (FU-Berlin, TNO Utrecht)
[3190]	63	!
	64	! Description:
	65	! ------------
	66	!> MODULE for reading-in Chemistry Emissions data
	67	!>
	68	!> @todo Check_parameters routine should be developed: for now it includes just one condition
	69	!> @todo Use of Restart files not contempled at the moment
[3458]	70	!> @todo revise indices of files read from the netcdf: preproc_emission_data and expert_emission_data
	71	!> @todo for now emission data may be passed on a singular vertical level: need to be more flexible
[3190]	72	!> @note <Enter notes on the module>
	73	!> @bug <Enter known bugs here>
	74	!------------------------------------------------------------------------------!
	75
	76	MODULE chem_emissions_mod
	77
[3266]	78	USE arrays_3d, &
	79	ONLY: rho_air
	80
[3190]	81	USE control_parameters, &
[3458]	82	ONLY: initializing_actions, end_time, message_string, &
[3190]	83	intermediate_timestep_count, dt_3d
	84
	85	USE indices
	86
	87	USE kinds
	88
[3458]	89	#if defined ( __netcdf )
[3483]	90	USE NETCDF
	91	#endif
[3458]	92
[3190]	93	USE netcdf_data_input_mod, &
	94	ONLY: chem_emis_att_type, chem_emis_val_type
	95
	96	USE date_and_time_mod, &
[3458]	97	ONLY: time_default_indices, time_preprocessed_indices, &
	98	month_of_year, day_of_month, hour_of_day, &
[3190]	99	index_mm, index_dd, index_hh
	100
	101	USE chem_gasphase_mod, &
	102	ONLY: spc_names
	103
	104	USE chem_modules
	105
	106	USE statistics, &
	107	ONLY: weight_pres
	108
	109	IMPLICIT NONE
	110
	111	!-- Declare all global variables within the module
	112
	113	CHARACTER (LEN=80) :: filename_emis !> Variable for the name of the netcdf input file
	114
	115	INTEGER(iwp) :: i !> index 1st selected dimension (some dims are not spatial)
	116	INTEGER(iwp) :: j !> index 2nd selected dimension
	117	INTEGER(iwp) :: i_start !> Index to start read variable from netcdf along one dims
	118	INTEGER(iwp) :: i_end !> Index to end read variable from netcdf in one dims
	119	INTEGER(iwp) :: j_start !> Index to start read variable from netcdf in additional dims
	120	INTEGER(iwp) :: j_end !> Index to end read variable from netcdf in additional dims
	121	INTEGER(iwp) :: z_start !> Index to start read variable from netcdf in additional dims
	122	INTEGER(iwp) :: z_end !> Index to end read variable from netcdf in additional dims
	123	INTEGER(iwp) :: dt_emis !> Time Step Emissions
	124	INTEGER(iwp) :: len_index !> length of index (used for several indices)
	125	INTEGER(iwp) :: len_index_voc !> length of voc index
	126	INTEGER(iwp) :: len_index_pm !> length of PMs index
	127	REAL(wp) :: con_factor !> Units Conversion Factor
	128
	129
	130	! ---------------------------------------------------------------
	131	! Set Values of molecules, mols, etc
	132	! ---------------------------------------------------------------
	133
	134	!> Avogadro number
	135	REAL, PARAMETER :: Avog = 6.02205e23 ! mlc/mol
	136
	137	!> Dobson units:
	138	REAL, PARAMETER :: Dobs = 2.68668e16 ! (mlc/cm2) / DU
	139
	140	!> sesalt composition:
	141	! (Seinfeld and Pandis, "Atmospheric Chemistry and Physics",
	142	! table 7.8 "Composition of Sea-Salt", p. 444)
	143	REAL, PARAMETER :: massfrac_Cl_in_seasalt = 0.5504 ! (kg Cl )/(kg seasalt)
	144	REAL, PARAMETER :: massfrac_Na_in_seasalt = 0.3061 ! (kg Na )/(kg seasalt)
	145	REAL, PARAMETER :: massfrac_SO4_in_seasalt = 0.0768 ! (kg SO4)/(kg seasalt)
	146
	147	!> other numbers
	148	REAL, PARAMETER :: xm_seasalt = 74.947e-3 ! kg/mol : NaCl, SO4, ..
	149	REAL, PARAMETER :: rho_seasalt = 2.2e3 ! kg/m3
	150
	151	!> * amonium sulphate
	152
	153	REAL, PARAMETER :: xm_NH4HSO4 = xm_NH4 + xm_H + xm_SO4 ! kg/mol
	154	REAL, PARAMETER :: rho_NH4HSO4a = 1.8e3 ! kg/m3
	155
	156	! ---------------------------------------------------------------
	157	! gas
	158	! ---------------------------------------------------------------
	159
	160	!> gas constant
	161	REAL, PARAMETER :: Rgas = 8.3144 ! J/mol/K
	162
	163	!> gas constant for dry air
	164	REAL, PARAMETER :: Rgas_air = Rgas / xm_air ! J/kg/K
	165
	166	!> water vapour
	167	REAL, PARAMETER :: Rgas_h2o = Rgas / xm_h2o ! J/kg/K
	168
	169	!> standard pressure
	170	REAL, PARAMETER :: p0 = 1.0e5 ! Pa
	171
	172	!> sea-level pressure
	173	REAL, PARAMETER :: p_sealevel = 1.01325e5 ! Pa <-- suggestion Bram Bregman
	174
	175	!> global mean pressure
	176	REAL, PARAMETER :: p_global = 98500.0 ! Pa
	177
	178	!> specific heat of dry air at constant pressure
	179	REAL, PARAMETER :: cp_air = 1004.0 ! J/kg/K
	180
	181	!> Latent heat of evaporation
	182	REAL, PARAMETER :: lvap = 2.5e6 ! [J kg-1]
	183
	184	!> Latent heat of condensation at 0 deg Celcius
	185	! (heat (J) necesarry to evaporate 1 kg of water)
	186	REAL, PARAMETER :: Lc = 22.6e5 ! J/kg
	187
	188	!> kappa = R/cp = 2/7
	189	REAL, PARAMETER :: kappa = 2.0/7.0
	190
	191	!> Von Karman constant (dry_dep)
	192	REAL, PARAMETER :: vkarman = 0.4
	193
	194	!> Boltzmann constant:
[3483]	195	REAL(wp), PARAMETER :: kbolz = 1.38066e-23_wp ! J/K
[3190]	196
	197	!> Inverse Reference Pressure (1/Pa)
	198	REAL(wp), PARAMETER :: pref_i = 1.0_wp / 100000.0_wp
	199
	200	!>
	201	REAL(wp), PARAMETER :: r_cp = 0.286_wp !< R / cp (exponent for potential temperature)
	202
	203
	204	SAVE
	205
	206	!-- Interfaces Part
	207
	208	!-- Checks Input parameters
	209	INTERFACE chem_emissions_check_parameters
	210	MODULE PROCEDURE chem_emissions_check_parameters
	211	END INTERFACE chem_emissions_check_parameters
	212
	213	!-- Reading of NAMELIST parameters
	214	! INTERFACE chem_emissions_parin
	215	! MODULE PROCEDURE chem_emissions_parin
	216	! END INTERFACE chem_emissions_parin
	217
	218	!-- Output of information to the header file
	219	! INTERFACE chem_emissions_header
	220	! MODULE PROCEDURE chem_emissions_header
	221	! END INTERFACE chem_emissions_header
	222
	223	!-- Matching Emissions actions
	224	INTERFACE chem_emissions_match
	225	MODULE PROCEDURE chem_emissions_match
	226	END INTERFACE chem_emissions_match
	227
	228	!-- Initialization actions
	229	INTERFACE chem_emissions_init
	230	MODULE PROCEDURE chem_emissions_init
	231	END INTERFACE chem_emissions_init
	232
	233	!-- Setup of Emissions
	234	INTERFACE chem_emissions_setup
	235	MODULE PROCEDURE chem_emissions_setup
	236	END INTERFACE chem_emissions_setup
	237
	238
	239	!-- Public Interfaces
	240	PUBLIC chem_emissions_init,chem_emissions_match,chem_emissions_setup
	241
	242	!-- Public Variables
	243
[3458]	244	PUBLIC con_factor, len_index,len_index_voc,len_index_pm
	245
[3190]	246	CONTAINS
	247
	248	!------------------------------------------------------------------------------!
	249	! Description:
	250	! ------------
	251	!> Routine for checking input parameters
	252	!------------------------------------------------------------------------------!
	253	SUBROUTINE chem_emissions_check_parameters
	254
	255
[3458]	256	!TBD: Where Should we put the call to chem_emissions_check_parameters? In chem_init or in check_parameters?
[3190]	257
	258	IMPLICIT NONE
	259
	260	INTEGER(iwp) :: tmp
	261
	262	TYPE(chem_emis_att_type) :: emt
	263	!
	264	!-- Call routine for controlling chem_emissions namelist
	265	! CALL chem_emissions_parin
	266
[3228]	267	!TBD: In case the given emission values do not coincide with the passed names we should think of a solution:
	268	! I would personally do that if the name of the species differ from the number of emission values, I would
	269	! print a warning that says that in correspondance of that particular species the value is zero.
	270	! An alternative would be to initialize the cs_surface_flux values to a negative number.
[3190]	271
	272	!-- Check Emission Species Number equal to number of passed names for the chemistry species:
	273	IF ( SIZE(emt%species_name) /= emt%nspec ) THEN
	274
	275	message_string = 'Numbers of input emission species names and number of species' // &
	276	'for which emission values are given do not match'
[3281]	277	CALL message( 'chem_emissions_check_parameters', 'CM0437', 2, 2, 0, 6, 0 )
[3190]	278
	279
	280	ENDIF
	281
	282	!-- Check Emission Species Number equals to number of passed names for the species
	283	!IF ( SIZE(emt%species_name) /= SIZE(emt%species_index) ) THEN
	284	! message_string = 'Number of input emission species names and ' // &
	285	! ' number of passed species indices do not match'
	286	! CALL message( 'chem_emissions_check_parameters', 'CM0101', 2, 2, 0, 6, 0 )
	287
	288	!ENDIF
	289
	290	!-- Check Emission Categories
	291	! IF ( SIZE(chem_emis%cat_name) /= SIZE(chem_emis%cat_index) ) THEN
	292	! WRITE( message_string, * ) &
	293	! 'Number of input emission categories name and categories index does not match\\'
	294	! CALL message( 'chem_emissions_check_parameters', 'CM0101', 1, 2, 0, 6, 0 )
	295	! ENDIF
	296
	297
	298
	299	!TBD: Check which other consistency controls should be included
	300
[3228]	301	!TBD: Include check for spatial dimension of netcdf file. If they do not match with the ones
	302	! of the simulation, the model should print an error.
	303
[3190]	304	END SUBROUTINE chem_emissions_check_parameters
	305
	306	!------------------------------------------------------------------------------!
	307	! Description:
	308	! ------------
[3228]	309	!> Matching the chemical species indices. The routine checks what are the indices of the emission input species
	310	! and the corresponding ones of the model species. The routine gives as output a vector containing the number
	311	! of common species: it is important to note that while the model species are distinct, their values could be
	312	! 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.
[3190]	313	!------------------------------------------------------------------------------!
	314	SUBROUTINE chem_emissions_match(emt_att,len_index)
	315
	316
	317	INTEGER(iwp), INTENT(INOUT) :: len_index !< Variable where to store the number of common species between the input dataset and the model species
	318
	319	TYPE(chem_emis_att_type), INTENT(INOUT) :: emt_att !< Chemistry Emission Array containing information for all the input chemical emission species
	320
	321	INTEGER(iwp) :: ind_mod, ind_inp !< Parameters for cycling through chemical model and input species
[3221]	322	INTEGER(iwp) :: nspec_emis_inp !< Variable where to store the number of the emission species in input
[3190]	323
	324	INTEGER(iwp) :: ind_voc !< Indices to check whether a split for voc should be done
	325
[3266]	326	INTEGER(iwp) :: ispec !> index for cycle over effective number of emission species
[3190]	327
[3266]	328
[3190]	329	!> Tell the user where we are!!
	330	CALL location_message( 'Matching input emissions and model chemistry species', .FALSE. )
	331
	332	!> Number of input emission species.
	333
	334	nspec_emis_inp=emt_att%nspec
	335
[3266]	336	!> Check the emission mode: DEFAULT, PRE-PROCESSED or PARAMETERIZED !TBD: Add option for capital or small letters
[3190]	337	SELECT CASE(TRIM(mode_emis))
	338
[3266]	339	!> 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.
	340	CASE ("PRE-PROCESSED")
[3190]	341
[3266]	342	CALL location_message( 'chem_emissions PRE-PROCESSED_mode_matching', .FALSE. )
[3190]	343
	344	len_index=0
	345	len_index_voc=0
	346
	347	!> The first condition is that both the number of model and input emissions species are not null
	348	IF ( (SIZE(spc_names) .GT. 0) .AND. (nspec_emis_inp .GT. 0)) THEN
	349
	350	!> Cycle over model species
	351	DO ind_mod = 1, SIZE(spc_names)
	352	!> Cycle over Input species
	353	DO ind_inp = 1, nspec_emis_inp
	354
[3266]	355	!> 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
[3190]	356	! > Check for VOC Species: In input in this case we only have one species (VOC)
	357	IF (TRIM(emt_att%species_name(ind_inp))=="VOC") THEN
	358	!> 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
	359	DO ind_voc= 1,emt_att%nvoc
	360	!> Determine the indices of the coinciding species in the two cases and assign them to matching arrays
	361	IF (TRIM(emt_att%voc_name(ind_voc))==TRIM(spc_names(ind_mod))) THEN
	362	len_index=len_index+1
	363	len_index_voc=len_index_voc+1
	364	ENDIF
	365	END DO
	366	ENDIF
	367	!> Other Species
	368	IF (TRIM(emt_att%species_name(ind_inp))==TRIM(spc_names(ind_mod))) THEN
	369	len_index=len_index+1
	370	ENDIF
	371	ENDDO
	372	ENDDO
	373
	374	!> Allocate array for storing the indices of the matched species
	375	IF (len_index>0) THEN
	376
	377	ALLOCATE(match_spec_input(len_index))
	378
	379	ALLOCATE(match_spec_model(len_index))
	380
	381	IF (len_index_voc>0) THEN
	382
	383	ALLOCATE(match_spec_voc_model(len_index_voc)) !> Contains indices of the VOC model species
	384
	385	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
	386
	387	ENDIF
	388
	389	!> Repeat the same cycle of above but passing the species indices to the newly declared arrays
	390	len_index=0
	391
	392	!> Cycle over model species
	393	DO ind_mod = 1, SIZE(spc_names)
	394	!> Cycle over Input species
	395	DO ind_inp = 1, nspec_emis_inp
	396
	397	!> Determine the indices of the coinciding species in the two cases and assign them to matching arrays
	398
	399	! > VOCs
	400	IF ( TRIM(emt_att%species_name(ind_inp))=="VOC" .AND. ALLOCATED(match_spec_voc_input) ) THEN
	401	DO ind_voc= 1,emt_att%nvoc
	402	IF (TRIM(emt_att%voc_name(ind_voc))==TRIM(spc_names(ind_mod))) THEN
	403	len_index=len_index+1
	404	len_index_voc=len_index_voc+1
	405
	406	match_spec_input(len_index)=ind_inp
	407	match_spec_model(len_index)=ind_mod
	408
	409	match_spec_voc_input(len_index_voc)=ind_voc
	410	match_spec_voc_model(len_index_voc)=ind_mod
	411	ENDIF
	412	END DO
	413	ENDIF
	414
	415	!>Other Species
	416	IF (TRIM(emt_att%species_name(ind_inp))==TRIM(spc_names(ind_mod))) THEN
	417	len_index=len_index+1
	418	match_spec_input(len_index)=ind_inp
	419	match_spec_model(len_index)=ind_mod
	420	ENDIF
	421	END DO
	422	END DO
	423
	424	!> In the case there are no species matching, the emission module should not be called
	425	ELSE
	426
	427	message_string = 'Given Chemistry Emission Species' // &
[3266]	428	' DO NOT MATCH' // &
	429	' model chemical species:' // &
	430	' Chemistry Emissions routine is not called'
[3281]	431	CALL message( 'chem_emissions_matching', 'CM0438', 0, 0, 0, 6, 0 )
[3190]	432	ENDIF !> IF (len_index>0)
	433
	434	ELSE
	435
	436	!In this case either spc_names is null or nspec_emis_inp is not allocated
	437
	438	message_string = 'Array of Emission species not allocated:' // &
[3266]	439	' Either no emission species are provided as input or' // &
	440	' no chemical species are used by PALM:' // &
	441	' Chemistry Emissions routine is not called'
[3281]	442	CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 )
[3190]	443
	444	ENDIF !> IF ( (SIZE(spc_names) .GT. 0) .AND. ALLOCATED(nspec_emis_inp))
	445
	446	!> ------------------------------------------------------------------
	447	!> DEFAULT case
	448
	449	CASE ("DEFAULT")
	450
	451	CALL location_message( 'chem_emissions DEFAULT_mode_matching', .FALSE. )
	452
	453	len_index=0 !> index for TOTAL number of species
	454	len_index_voc=0 !> index for TOTAL number of VOCs
	455	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.
	456
	457	!> The first condition is that both the number of model and input emissions species are not null
	458	IF ( (SIZE(spc_names) .GT. 0) .AND. (nspec_emis_inp .GT. 0) ) THEN
	459
	460	!> Cycle over model species
	461	DO ind_mod = 1, SIZE(spc_names)
	462	!> Cycle over input species
	463	DO ind_inp = 1, nspec_emis_inp
	464
	465	! > Check for VOC Species: In input in this case we only have one species (VOC)
	466	IF (TRIM(emt_att%species_name(ind_inp))=="VOC") THEN
	467	!> 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
	468	DO ind_voc= 1,emt_att%nvoc
	469	!> Determine the indices of the coinciding species in the two cases and assign them to matching arrays
	470	IF (TRIM(emt_att%voc_name(ind_voc))==TRIM(spc_names(ind_mod))) THEN
	471	len_index=len_index+1
	472	len_index_voc=len_index_voc+1
	473	ENDIF
	474	END DO
	475	ENDIF
	476
	477	!> 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.
	478	IF (TRIM(emt_att%species_name(ind_inp))=="PM") THEN
	479	!>pm1
	480	IF (TRIM(spc_names(ind_mod))=="PM1") THEN
	481	len_index=len_index+1
	482	!>pm2.5
	483	ELSE IF (TRIM(spc_names(ind_mod))=="PM25") THEN
	484	len_index=len_index+1
	485	!>pm10
	486	ELSE IF (TRIM(spc_names(ind_mod))=="PM10") THEN
	487	len_index=len_index+1
	488	ENDIF
	489	ENDIF
	490
	491	!> 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
	492	IF (TRIM(emt_att%species_name(ind_inp))=="NOx") THEN
	493	!>no
	494	IF (TRIM(spc_names(ind_mod))=="NO") THEN
	495	len_index=len_index+1
	496	!>no2
	497	ELSE IF (TRIM(spc_names(ind_mod))=="NO2") THEN
	498	len_index=len_index+1
	499	ENDIF
	500	ENDIF
	501
	502	!>SOX: same as for NOx, but with SO2 and SO4
	503	IF (TRIM(emt_att%species_name(ind_inp))=="SOx") THEN
	504	!>no
	505	IF (TRIM(spc_names(ind_mod))=="SO2") THEN
	506	len_index=len_index+1
	507	!>no2
	508	ELSE IF (TRIM(spc_names(ind_mod))=="SO4") THEN
	509	len_index=len_index+1
	510	ENDIF
	511	ENDIF
	512
	513	!> Other Species
	514	IF (TRIM(emt_att%species_name(ind_inp))==TRIM(spc_names(ind_mod))) THEN
	515	len_index=len_index+1
	516	ENDIF
	517	END DO !> number of emission input species
	518	END DO !> number of model species
	519
	520
	521	!> Allocate Arrays
	522	IF (len_index>0) THEN
	523
	524	ALLOCATE(match_spec_input(len_index))
	525	ALLOCATE(match_spec_model(len_index))
	526
	527	IF (len_index_voc>0) THEN
	528	ALLOCATE(match_spec_voc_model(len_index_voc)) !> contains indices of the VOC model species
[3458]	529	ALLOCATE(match_spec_voc_input(len_index_voc)) !> In input there is only one value for VOCs in the DEFAULT mode.
	530	! This array contains the indices of the different values of VOC compositions of the input variable VOC_composition
[3190]	531	ENDIF
	532
	533	!> ASSIGN INDICES
	534	!> Repeat the same cycles of above, but passing the species indices to the newly declared arrays
	535	len_index=0
	536	len_index_voc=0
	537
	538	DO ind_mod = 1, SIZE(spc_names)
	539	DO ind_inp = 1, nspec_emis_inp
	540
	541	! > VOCs
	542	IF ( TRIM(emt_att%species_name(ind_inp))=="VOC" .AND. ALLOCATED(match_spec_voc_input) ) THEN
	543	DO ind_voc= 1,emt_att%nvoc
	544	IF (TRIM(emt_att%voc_name(ind_voc))==TRIM(spc_names(ind_mod))) THEN
	545	len_index=len_index+1
	546	len_index_voc=len_index_voc+1
	547
	548	match_spec_input(len_index)=ind_inp
	549	match_spec_model(len_index)=ind_mod
	550
	551	match_spec_voc_input(len_index_voc)=ind_voc
	552	match_spec_voc_model(len_index_voc)=ind_mod
	553	ENDIF
	554	END DO
	555	ENDIF
	556
	557	!> 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.
	558	IF (TRIM(emt_att%species_name(ind_inp))=="PM") THEN
	559	!>pm1
	560	IF (TRIM(spc_names(ind_mod))=="PM1") THEN
	561	len_index=len_index+1
	562
	563	match_spec_input(len_index)=ind_inp
	564	match_spec_model(len_index)=ind_mod
	565
	566	!match_spec_pm(1)=ind_mod
	567	!>pm2.5
	568	ELSE IF (TRIM(spc_names(ind_mod))=="PM25") THEN
	569	len_index=len_index+1
	570
	571	match_spec_input(len_index)=ind_inp
	572	match_spec_model(len_index)=ind_mod
	573
	574	!match_spec_pm(2)=ind_mod
	575	!>pm10
	576	ELSE IF (TRIM(spc_names(ind_mod))=="PM10") THEN
	577	len_index=len_index+1
	578
	579	match_spec_input(len_index)=ind_inp
	580	match_spec_model(len_index)=ind_mod
	581
	582	!match_spec_pm(3)=ind_mod
	583	ENDIF
	584	ENDIF
	585
	586	!> NOx - The same as for PMs but here the model species are only 2: NO and NO2
	587	IF (TRIM(emt_att%species_name(ind_inp))=="NOx") THEN
	588	!>no
	589	IF (TRIM(spc_names(ind_mod))=="NO") THEN
	590	len_index=len_index+1
	591
	592	match_spec_input(len_index)=ind_inp
	593	match_spec_model(len_index)=ind_mod
	594
	595	!match_spec_nox(1)=ind_mod
	596	!>no2
	597	ELSE IF (TRIM(spc_names(ind_mod))=="NO2") THEN
	598	len_index=len_index+1
	599
	600	match_spec_input(len_index)=ind_inp
	601	match_spec_model(len_index)=ind_mod
	602
	603	! match_spec_nox(2)=ind_mod
	604	ENDIF
	605	ENDIF
	606
	607	!> SOx
	608	IF (TRIM(emt_att%species_name(ind_inp))=="SOx") THEN
	609	!>so2
	610	IF (TRIM(spc_names(ind_mod))=="SO2") THEN
	611	len_index=len_index+1
	612
	613	match_spec_input(len_index)=ind_inp
	614	match_spec_model(len_index)=ind_mod
	615
	616	! match_spec_sox(1)=ind_mod
	617	!>so4
	618	ELSE IF (TRIM(spc_names(ind_mod))=="SO4") THEN
	619	len_index=len_index+1
	620
	621	match_spec_input(len_index)=ind_inp
	622	match_spec_model(len_index)=ind_mod
	623
	624	! match_spec_sox(2)=ind_mod
	625	ENDIF
	626	ENDIF
	627
	628	!> Other Species
	629	IF (TRIM(emt_att%species_name(ind_inp))==TRIM(spc_names(ind_mod))) THEN
	630	len_index=len_index+1
	631
	632	match_spec_input(len_index)=ind_inp
	633	match_spec_model(len_index)=ind_mod
	634	ENDIF
	635	END DO
	636	END DO
	637
	638	ELSE
	639
	640	message_string = 'Given Chemistry Emission Species' // &
[3266]	641	' DO NOT MATCH' // &
	642	' model chemical species' // &
[3458]	643	' Chemistry Emissions routine is not called'
[3281]	644	CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 )
[3190]	645
	646	ENDIF
	647
	648	ELSE
	649
	650	message_string = 'Array of Emission species not allocated: ' // &
[3266]	651	' Either no emission species are provided as input or' // &
	652	' no chemical species are used by PALM:' // &
	653	' Chemistry Emissions routine is not called'
[3281]	654	CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 )
[3190]	655
	656	ENDIF
	657
	658	!-- 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.
	659
	660	CASE ("PARAMETERIZED")
	661
	662	len_index=0
	663
	664	!spc_names have to be greater than zero for proceeding
	665	IF ( (SIZE(spc_names) .GT. 0) .AND. (nspec_emis_inp .GT. 0) ) THEN
	666
[3266]	667
[3190]	668	!cycle over model species
	669	DO ind_mod = 1, SIZE(spc_names)
	670	ind_inp=1
	671	DO WHILE ( TRIM( surface_csflux_name( ind_inp ) ) /= 'novalue' ) !<'novalue' is the default
	672	IF (TRIM(surface_csflux_name( ind_inp ))==TRIM(spc_names(ind_mod))) THEN
	673	len_index=len_index+1
	674	ENDIF
	675	ind_inp=ind_inp+1
	676
	677	ENDDO
	678	ENDDO
	679
	680	!Proceed only if there are matched species
	681	IF ( len_index .GT. 0 ) THEN
	682
	683	!Allocation of Arrays of the matched species
	684	ALLOCATE(match_spec_input(len_index))
	685
	686	ALLOCATE(match_spec_model(len_index))
	687
	688	!Assign corresponding indices of input and model matched species
	689	len_index=0
	690
	691	DO ind_mod = 1, SIZE(spc_names)
	692	ind_inp = 1
	693	DO WHILE ( TRIM( surface_csflux_name( ind_inp ) ) /= 'novalue' ) !<'novalue' is the default
	694	IF (TRIM( surface_csflux_name( ind_inp ) ) == TRIM(spc_names(ind_mod))) THEN
	695	len_index=len_index+1
	696	match_spec_input(len_index)=ind_inp
	697	match_spec_model(len_index)=ind_mod
	698	ENDIF
	699	ind_inp=ind_inp+1
	700	END DO
	701	END DO
	702
[3266]	703	! also, Add AN if to check that when the surface_csflux are passed to the namelist, also the street factor values are provided
	704	DO ispec = 1 , len_index
	705
[3458]	706	IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND. emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN
[3266]	707
[3458]	708	message_string = 'PARAMETERIZED emissions mode selected:' // &
[3266]	709	' EMISSIONS POSSIBLE ONLY ON STREET SURFACES' // &
	710	' but values of scaling factors for street types' // &
	711	' emiss_factor_main AND emiss_factor_side' // &
	712	' not provided for each of the emissions species' // &
	713	' or not provided at all: PLEASE set a finite value' // &
	714	' for these parameters in the chemistry namelist'
[3281]	715	CALL message( 'chem_emissions_matching', 'CM0442', 2, 2, 0, 6, 0 )
[3266]	716	ENDIF
	717	END DO
	718
	719
[3190]	720	ELSE
	721
	722	message_string = 'Given Chemistry Emission Species' // &
[3266]	723	' DO NOT MATCH' // &
	724	' model chemical species' // &
	725	' Chemistry Emissions routine is not called'
[3281]	726	CALL message( 'chem_emissions_matching', 'CM0443', 0, 0, 0, 6, 0 )
[3190]	727	ENDIF
	728
	729	ELSE
[3266]	730
[3190]	731	message_string = 'Array of Emission species not allocated: ' // &
[3266]	732	' Either no emission species are provided as input or' // &
	733	' no chemical species are used by PALM.' // &
	734	' Chemistry Emissions routine is not called'
[3281]	735	CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 )
[3190]	736
	737	ENDIF
	738
	739
	740	!-- CASE when emission module is switched on but mode_emis is not specified or it is given the wrong name
	741	CASE DEFAULT
	742
	743	message_string = 'Chemistry Emissions Module switched ON, but' // &
[3266]	744	' either no emission mode specified or incorrectly given :' // &
	745	' please, pass the correct value to the namelist parameter "mode_emis"'
[3281]	746	CALL message( 'chem_emissions_matching', 'CM0445', 2, 2, 0, 6, 0 )
[3190]	747
	748	END SELECT
	749
	750	END SUBROUTINE chem_emissions_match
	751
	752	!------------------------------------------------------------------------------!
	753	! Description:
	754	! ------------
	755	!> Initialization:
	756	!> Netcdf reading, arrays allocation and first calculation of cssws fluxes at timestep 0
	757	!>
	758	!------------------------------------------------------------------------------!
	759	SUBROUTINE chem_emissions_init(emt_att,emt,nspec_out)
	760
	761	USE surface_mod, &
	762	ONLY: surf_lsm_h,surf_def_h,surf_usm_h
	763
	764	IMPLICIT NONE
	765
	766	TYPE(chem_emis_att_type),INTENT(INOUT) :: emt_att !> variable where to store all the information of
	767	! emission inputs whose values do not depend
	768	! on the considered species
	769
	770	TYPE(chem_emis_val_type),INTENT(INOUT), ALLOCATABLE, DIMENSION(:) :: emt !> variable where to store emission inputs values,
	771	! depending on the considered species
	772
	773	INTEGER(iwp),INTENT(INOUT) :: nspec_out !> number of outputs of chemistry emission routines
	774
	775	CHARACTER (LEN=80) :: units !> units of chemistry inputs
	776
	777	INTEGER(iwp) :: ispec !> Index to go through the emission chemical species
	778
[3458]	779
[3190]	780	!-- Actions for initial runs : TBD: needs to be updated
	781	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	782	!-- ...
	783
	784	!
	785	!-- Actions for restart runs
	786	ELSE
	787	!-- ...
	788
	789	ENDIF
	790
	791
	792	CALL location_message( 'Starting initialization of chemistry emissions arrays', .FALSE. )
	793
	794
	795	!-- Match Input and Model emissions
	796	CALL chem_emissions_match(emt_att,nspec_out)
	797
	798	!> If nspec_out==0, then do not use emission module: The corresponding message is already printed in the matching routine
	799	IF ( nspec_out == 0 ) THEN
	800
	801	emission_output_required=.FALSE.
	802
	803	ELSE
	804
	805	emission_output_required=.TRUE.
	806
	807
	808	!-----------------------------------------------------------------
	809	!> Set molecule masses'
	810	ALLOCATE(emt_att%xm(nspec_out))
	811	! loop over emisisons:
	812
	813	DO ispec = 1, nspec_out
	814	! switch:
	815	SELECT CASE ( TRIM(spc_names(match_spec_model(ispec))) )
	816	CASE ( 'SO2','SO4' ) ; emt_att%xm(ispec) = xm_S + xm_O * 2 ! kg/mole
	817	CASE ( 'NO','NO2' ) ; emt_att%xm(ispec) = xm_N + xm_O * 2 ! kg/mole NO2 equivalent
	818	CASE ( 'NH3' ) ; emt_att%xm(ispec) = xm_N + xm_H * 3 ! kg/mole
	819	CASE ( 'CO' ) ; emt_att%xm(ispec) = xm_C + xm_O ! kg/mole
	820	CASE ( 'CO2' ) ; emt_att%xm(ispec) = xm_C + xm_O * 2 ! kg/mole
	821	CASE ( 'CH4' ) ; emt_att%xm(ispec) = xm_C + xm_H * 4 ! kg/mole
	822	CASE ( 'HNO3' ); emt_att%xm(ispec) = xm_H + xm_N + xm_O*3 !kg/mole
	823	CASE DEFAULT
	824	!-- TBD: check if this hase to be removed
	825	!emt_att%xm(ispec) = 1.0_wp
	826	END SELECT
	827	ENDDO
	828
	829
	830	!> ASSIGN EMISSION VALUES for the different emission modes
	831	SELECT CASE(TRIM(mode_emis)) !TBD: Add the option for CApital or small letters
	832
	833
[3266]	834	!> PRE-PROCESSED case
	835	CASE ("PRE-PROCESSED")
[3190]	836
	837	IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(nzb:nzt+1,0:ny,0:nx,nspec_out))
	838
[3266]	839	CALL location_message( 'emis_distribution array allocated in PRE-PROCESSED mode', .FALSE. )
[3190]	840
	841	!> Calculate the values of the emissions at the first time step
	842	CALL chem_emissions_setup(emt_att,emt,nspec_out)
	843
[3458]	844	!> Default case
	845	CASE ("DEFAULT")
	846
	847	IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution( 1, 0:ny, 0:nx, nspec_out ) )
	848
	849	CALL location_message( 'emis_distribution array allocated in DEFAULT mode', .FALSE. )
	850
	851	!> Calculate the values of the emissions at the first time step
	852	CALL chem_emissions_setup(emt_att,emt,nspec_out)
	853
	854	!> PARAMETERIZED case
[3190]	855	CASE ("PARAMETERIZED")
	856
	857	CALL location_message( 'emis_distribution array allocated in PARAMETERIZED mode', .FALSE. )
	858
[3458]	859	! For now for PAR and DEF values only, first vertical level of emis_distribution is allocated, while for PRE-PROCESSED all.
[3190]	860	IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(1,0:ny,0:nx,nspec_out))
	861
	862	!> Calculate the values of the emissions at the first time step
	863	CALL chem_emissions_setup(emt_att,emt,nspec_out)
	864
	865	END SELECT
	866
	867	ENDIF
	868
	869	END SUBROUTINE chem_emissions_init
	870
	871
	872
	873	!------------------------------------------------------------------------------!
	874	! Description:
	875	! ------------
	876	!> Routine for Update of Emission values at each timestep
	877	!-------------------------------------------------------------------------------!
	878
	879	SUBROUTINE chem_emissions_setup(emt_att,emt,nspec_out)
	880
[3266]	881	USE arrays_3d, &
	882	ONLY: dzw
[3190]	883	USE grid_variables, &
	884	ONLY: dx, dy
	885	USE indices, &
	886	ONLY: nnx,nny,nnz
[3467]	887	USE salsa_mod, &
	888	ONLY: salsa
[3190]	889	USE surface_mod, &
	890	ONLY: surf_lsm_h,surf_def_h,surf_usm_h
	891	USE netcdf_data_input_mod, &
	892	ONLY: street_type_f
	893	USE arrays_3d, &
	894	ONLY: hyp, pt
	895
	896	IMPLICIT NONE
	897
	898	!--- IN/OUT
	899
	900	TYPE(chem_emis_att_type),INTENT(INOUT) :: emt_att !> variable where to store all the information of
	901	! emission inputs whose values do not depend
	902	! on the considered species
	903
	904	TYPE(chem_emis_val_type),INTENT(INOUT), ALLOCATABLE, DIMENSION(:) :: emt !> variable where to store emission inputs values,
	905	! depending on the considered species
	906
[3266]	907	INTEGER,INTENT(IN) :: nspec_out !> Output of routine chem_emis_matching with number
	908	! of matched species
[3190]	909	!---
	910
	911	REAL(wp),ALLOCATABLE,DIMENSION(:,:) :: delta_emis !> Term to add to the emission_distribution array
	912	! for each of the categories at each time step
	913	! in the default mode
	914
	915	CHARACTER(LEN=80) :: units !> Units of the emissions
	916
	917	INTEGER(iwp) :: icat !> Index for number of categories
	918	INTEGER(iwp) :: ispec !> index for number of species
	919	INTEGER(iwp) :: i_pm_comp !> index for number of PM components
	920	INTEGER(iwp) :: ivoc !> Index for number of components of VOCs
	921	INTEGER(iwp) :: time_index !> Index for time
	922
[3266]	923	REAL(wp),ALLOCATABLE, DIMENSION(:) :: time_factor !> factor for time scaling of emissions
[3190]	924	REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: emis
	925
[3458]	926	REAL(wp), DIMENSION(24) :: par_emis_time_factor !< time factors
	927	! for the parameterized mode: these are fixed for each hour
	928	! of a single day.
[3266]	929	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: conv_to_ratio !> factor used for converting input
[3190]	930	! to adimensional concentration ratio
[3266]	931
[3190]	932	! ------------------------------------------
	933	! variables for the conversion of indices i and j according to surface_mod
	934
	935	INTEGER(iwp) :: i !> running index for grid in x-direction
	936	INTEGER(iwp) :: j !> running index for grid in y-direction
	937	INTEGER(iwp) :: k !> running index for grid in z-direction
	938	INTEGER(iwp) :: m !> running index for horizontal surfaces
	939
[3221]	940	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: tmp_temp
[3190]	941
	942	! --- const -------------------------------
	943	!-CONVERSION FACTORS: TIME
[3458]	944	! number of sec per hour = 3600
	945	REAL, PARAMETER :: s_per_hour = 3600.0 ! (s)/(hour)
	946	! number of sec per day = 86400
	947	REAL, PARAMETER :: s_per_day = 86400.0 ! (s)/(day)
[3190]	948	! number of hours in a year of 365 days:
[3458]	949	REAL, PARAMETER :: hour_per_year = 8760.0 !> TBD: What about leap years?
	950	! number of hours in a day:
	951	REAL, PARAMETER :: hour_per_day = 24.0
	952
[3190]	953	! conversion from hours to seconds (s/hour) = 1/3600.0 ~ 0.2777778
[3458]	954	REAL, PARAMETER :: hour_to_s = 1/s_per_hour ! (hour)/(s)
[3190]	955	! conversion from day to seconds (s/day) = 1/86400.0 ~ 1.157407e-05
[3458]	956	REAL, PARAMETER :: day_to_s = 1/s_per_day ! (day)/(s)
[3190]	957	! conversion from year to sec (s/year) = 1/31536000.0 ~ 3.170979e-08
[3458]	958	REAL, PARAMETER :: year_to_s = 1/(s_per_hour*hour_per_year) ! (year)/(s)
[3190]	959
	960	!-CONVERSION FACTORS: WEIGHT
	961	! Conversion from tons to kg (kg/tons) = 100.0/1 ~ 100
	962	REAL, PARAMETER :: tons_to_kg = 100 ! (tons)/(kg)
	963	! Conversion from g to kg (kg/g) = 1/1000.0 ~ 0.001
	964	REAL, PARAMETER :: g_to_kg = 0.001 ! (g)/(kg)
	965	! Conversion from g to kg (kg/g) = 1/1000.0 ~ 0.001
	966	REAL, PARAMETER :: miug_to_kg = 0.000000001 ! (microng)/(kg)
	967
	968	!< CONVERSION FACTORS: fraction to ppm
[3266]	969	REAL(wp), PARAMETER :: ratio2ppm = 1.0e06_wp
[3190]	970	!------------------------------------------------------
[3458]	971
[3190]	972	IF ( emission_output_required ) THEN
	973
[3458]	974	!> 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.
[3190]	975
	976	IF ( call_chem_at_all_substeps ) THEN
	977
	978	dt_emis = dt_3d * weight_pres(intermediate_timestep_count)
	979
	980	ELSE
	981
	982	dt_emis = dt_3d
	983
	984	ENDIF
	985
	986
	987	! --- CHECK UNITS
	988	!>-----------------------------------------------------
[3286]	989	!> Conversion of the units to the ones employed in PALM.
	990	!> Possible temporal units of input emissions data are: year, hour and second;
	991	!> the mass could be expressed in: tons, kilograms or grams.
	992	!> IN the PARAMETERIZED mode no conversion is possible: in this case INPUTS have to have fixed units.
[3190]	993	!>-----------------------------------------------------
	994
[3266]	995	IF (TRIM(mode_emis)=="DEFAULT" .OR. TRIM(mode_emis)=="PRE-PROCESSED" ) THEN
[3190]	996
	997	SELECT CASE(TRIM(emt_att%units))
	998	!> kilograms
	999	CASE ( 'kg/m2/s','KG/M2/S')
	1000	CALL location_message( 'Units of the emissions are consistent with' // &
	1001	' the ones employed in the PALM-4U model ', .FALSE. )
	1002	con_factor=1
	1003
	1004	CASE ('kg/m2/hour','KG/M2/HOUR')
	1005	CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
	1006
	1007	con_factor=hour_to_s
	1008
	1009	CASE ( 'kg/m2/day','KG/M2/DAY' )
	1010	CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
	1011
	1012	con_factor=day_to_s
	1013
	1014	CASE ( 'kg/m2/year','KG/M2/YEAR' )
	1015	CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
	1016
	1017	con_factor=year_to_s
	1018
	1019	!> Tons
	1020	CASE ( 'ton/m2/s','TON/M2/S' )
	1021	CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
	1022
	1023	con_factor=tons_to_kg
	1024
	1025	CASE ( 'ton/m2/hour','TON/M2/HOUR' )
	1026	CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
	1027
	1028	con_factor=tons_to_kg*hour_to_s
	1029
	1030	CASE ( 'ton/m2/year','TON/M2/YEAR' )
	1031	CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
	1032
	1033	con_factor=tons_to_kg*year_to_s
	1034
	1035	!> Grams
	1036	CASE ( 'g/m2/s','G/M2/S' )
	1037	CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
	1038
	1039	con_factor=g_to_kg
	1040
	1041	CASE ( 'g/m2/hour','G/M2/HOUR' )
	1042	CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
	1043
	1044	con_factor=g_to_kg*hour_to_s
	1045
	1046	CASE ( 'g/m2/year','G/M2/YEAR' )
	1047	CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
	1048
	1049	con_factor=g_to_kg*year_to_s
	1050
	1051	!> Micro Grams
	1052	CASE ( 'micrograms/m2/s','MICROGRAMS/M2/S' )
	1053	CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
	1054
	1055	con_factor=miug_to_kg
	1056
	1057	CASE ( 'micrograms/m2/hour','MICROGRAMS/M2/HOUR' )
	1058	CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
	1059
	1060	con_factor=miug_to_kg*hour_to_s
	1061
	1062	CASE ( 'micrograms/m2/year','MICROGRAMS/M2/YEAR' )
	1063	CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
	1064
	1065	con_factor=miug_to_kg*year_to_s
	1066
	1067	CASE DEFAULT
	1068	message_string = 'The Units of the provided input chemistry emission species' // &
	1069	' are not the ones required by PALM-4U: please check ' // &
[3266]	1070	' chemistry emission module documentation.'
[3312]	1071	CALL message( 'chem_emissions_setup', 'CM0446', 2, 2, 0, 6, 0 )
[3190]	1072
	1073	END SELECT
	1074
	1075
[3266]	1076	!> PRE-PROCESSED and parameterized mode
[3190]	1077	ELSE
	1078
	1079	message_string = 'No Units conversion required for units of chemistry emissions' // &
[3458]	1080	' of the PARAMETERIZED mode: units have to be provided in' // &
	1081	' micromole/m**2/day for GASES and' // &
	1082	' kg/m**2/day for PMs'
[3281]	1083	CALL message( 'chem_emissions_setup', 'CM0447', 0, 0, 0, 6, 0 )
[3190]	1084
	1085	ENDIF
	1086
	1087	!> Conversion factors (if the units are kg/m**2/s) we have to convert them to ppm/s
	1088	DO i=nxl,nxr
	1089	DO j=nys,nyn
	1090	!> Derive Temperature from Potential Temperature
	1091	tmp_temp(nzb:nzt+1,j,i) = pt(nzb:nzt+1,j,i) * ( hyp(nzb:nzt+1) * pref_i )**r_cp
	1092
	1093	!> We need to pass to cssws<- (ppm/s) * dz.
[3266]	1094	! Input is Nmole/(m^2*s)
[3190]	1095	! To go to ppmdz basically we need to multiply the input by (m2/N)dz
[3266]	1096	! (m*2/N)dz == V/N
[3190]	1097	! V/N=RT/P
	1098
[3458]	1099	!> m*3/Nmole (J/mol)K^-1 K Pa
[3221]	1100	conv_to_ratio(nzb:nzt+1,j,i) = ( (Rgas * tmp_temp(nzb:nzt+1,j,i)) / ((hyp(nzb:nzt+1))) )
[3190]	1101	ENDDO
	1102	ENDDO
	1103
	1104	!>------------------------------------------------
	1105	!> Start The Processing of the input data
	1106
	1107	emis_distribution(:,nys:nyn,nxl:nxr,:) = 0.0_wp
	1108
	1109	!>-----------------------------------------------
[3266]	1110	!> 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
[3190]	1111
[3266]	1112	!> PRE-PROCESSED MODE
	1113	IF (TRIM(mode_emis)=="PRE-PROCESSED") THEN
[3190]	1114
[3458]	1115	!> Update time indices
	1116	CALL time_preprocessed_indices(index_hh)
	1117
[3266]	1118	CALL location_message( 'PRE-PROCESSED MODE: No time-factor specification required', .FALSE. )
[3190]	1119
	1120	ELSEIF (TRIM(mode_emis)=="DEFAULT") THEN
	1121
	1122	CALL location_message( 'DEFAULT MODE: time-factor specification required', .FALSE. )
	1123
	1124	!> Allocate array where to store temporary emission values
	1125	IF(.NOT. ALLOCATED(emis)) ALLOCATE(emis(nys:nyn,nxl:nxr))
	1126
	1127	!> Allocate time factor per emitted component category
	1128	ALLOCATE(time_factor(emt_att%ncat))
	1129
	1130	!> Read-in HOURLY emission time factor
	1131	IF (TRIM(time_fac_type)=="HOUR") THEN
	1132
	1133	!> Update time indices
	1134	CALL time_default_indices(month_of_year,day_of_month,hour_of_day,index_hh)
	1135
	1136	!> Check if the index is less or equal to the temporal dimension of HOURLY emission files
	1137	IF (index_hh .LE. SIZE(emt_att%hourly_emis_time_factor(1,:))) THEN
	1138
	1139	!> Read-in the correspondant time factor
	1140	time_factor(:)= emt_att%hourly_emis_time_factor(:,index_hh)
	1141
	1142	ELSE
	1143
[3266]	1144	message_string = 'The "HOUR" time-factors (DEFAULT mode) of the chemistry emission species' // &
[3190]	1145	' are not provided for each hour of the total simulation time'
[3281]	1146	CALL message( 'chem_emissions_setup', 'CM0448', 2, 2, 0, 6, 0 )
[3190]	1147
	1148	ENDIF
	1149
	1150	!> Read-in MDH emissions time factors
	1151	ELSEIF (TRIM(time_fac_type)=="MDH") THEN
	1152
	1153	!> Update time indices
	1154	CALL time_default_indices(daytype_mdh,month_of_year,day_of_month,hour_of_day,index_mm,index_dd,index_hh)
	1155
	1156
	1157	!> Check if the index is less or equal to the temporal dimension of MDH emission files
	1158	IF ((index_hh+index_dd+index_mm) .LE. SIZE(emt_att%mdh_emis_time_factor(1,:))) THEN
	1159
	1160	!> Read-in the correspondant time factor
[3266]	1161	time_factor(:)=emt_att%mdh_emis_time_factor(:,index_mm)emt_att%mdh_emis_time_factor(:,index_dd) &
	1162	emt_att%mdh_emis_time_factor(:,index_hh)
[3190]	1163
	1164	ELSE
	1165
	1166	message_string = 'The "MDH" time-factors (DEFAULT mode) of the chemistry emission species' // &
	1167	' are not provided for each hour/day/month of the total simulation time'
[3281]	1168	CALL message( 'chem_emissions_setup', 'CM0449', 2, 2, 0, 6, 0 )
[3190]	1169
	1170	ENDIF
	1171
	1172	ELSE
	1173	!> condition when someone used the DEFAULT mode but forgets to indicate the time-factor in the namelist
	1174
	1175	message_string = 'The time-factor (DEFAULT mode) of the chemistry emission species' // &
	1176	' is not provided in the NAMELIST'
[3281]	1177	CALL message( 'chem_emissions_setup', 'CM0450', 2, 2, 0, 6, 0 )
[3190]	1178
	1179	ENDIF
	1180
	1181	!> PARAMETERIZED MODE
	1182	ELSEIF (TRIM(mode_emis)=="PARAMETERIZED") THEN
	1183	CALL location_message( 'PARAMETERIZED MODE: time-factor specification is fixed: ' // &
	1184	' 24 values for every day of the year ', .FALSE. )
	1185
[3458]	1186	!Assign Constant Values of time factors, diurnal time profile for traffic sector:
	1187	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, &
	1188	0.059, 0.061, 0.064, 0.067, 0.069, 0.069, 0.049, 0.039, 0.039, 0.029, 0.024, 0.019 /)
	1189
[3190]	1190	!> in this case allocate time factor each hour in a day
	1191	IF (.NOT. ALLOCATED(time_factor)) ALLOCATE(time_factor(1))
	1192
	1193	!>Pass the values of time factors in the parameterized mode to the time_factor variable. in this case index_hh==hour_of_day
	1194	index_hh=hour_of_day
	1195
[3458]	1196	time_factor(1) = par_emis_time_factor(index_hh)
[3190]	1197
	1198	ENDIF
	1199
	1200	!--------------------------------
	1201	!-- EMISSION DISTRIBUTION Calculation
	1202
	1203	!> PARAMETERIZED CASE
	1204	IF ( mode_emis == "PARAMETERIZED" ) THEN
	1205
	1206	DO ispec=1,nspec_out
	1207
[3458]	1208	!> Values are still micromoles/(m*2s). Units in this case are always micromoles/m*2day (or kilograms/m*2day for PMs)
	1209	emis_distribution(1,nys:nyn,nxl:nxr,ispec)=surface_csflux(match_spec_input(ispec))time_factor(1)hour_to_s
[3190]	1210
	1211	ENDDO
	1212
[3266]	1213	!> PRE-PROCESSED CASE
	1214	ELSEIF (TRIM(mode_emis)=="PRE-PROCESSED") THEN
[3190]	1215
	1216	!> Start Cycle over Species
	1217	DO ispec=1,nspec_out !> nspec_out represents the number of species in common between
	1218	! the emission input data and the chemistry mechanism used
[3458]	1219
	1220	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emt(match_spec_input(ispec))% &
	1221	preproc_emission_data(index_hh,1,nys+1:nyn+1,nxl+1:nxr+1)* &
	1222	con_factor
	1223
[3190]	1224	ENDDO
	1225
	1226	!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.
	1227
	1228	ELSE IF (TRIM(mode_emis)=="DEFAULT") THEN
	1229
	1230	!> Allocate array for the emission value corresponding to a specific category and time factor
	1231	ALLOCATE(delta_emis(nys:nyn,nxl:nxr))
	1232
	1233
	1234	!> Start Cycle over categories
	1235	DO icat = 1, emt_att%ncat
	1236
	1237	!> Start Cycle over Species
	1238	DO ispec=1,nspec_out !> nspec_out represents the number of species in common between
	1239	! the emission input data and the chemistry mechanism used
	1240
	1241	emis(nys:nyn,nxl:nxr) = emt(match_spec_input(ispec))%default_emission_data(icat,nys+1:nyn+1,nxl+1:nxr+1)
	1242
	1243	!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.
	1244
[3458]	1245	!> 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.
	1246	! So given any units, we convert them to seconds and finally multiply them by 24 ((value/sec)(243600)=value/day ---- (value/day)*time_factor=value/hour ---(value/hour)/(3600)=value/sec )
	1247	! ((value/sec)(243600)time_factor)/3600=24(value/sec)*time_factor
	1248
[3190]	1249	!> NOX Compositions
	1250	IF (TRIM(spc_names(match_spec_model(ispec)))=="NO") THEN
[3458]	1251	!> Kg/m2s kg/m2s
	1252	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)time_factor(icat)emt_att%nox_comp(icat,1)con_factorhour_per_day
[3190]	1253
[3458]	1254	emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
[3190]	1255
	1256	ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="NO2") THEN
	1257
[3458]	1258	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)time_factor(icat)emt_att%nox_comp(icat,2)con_factorhour_per_day
[3190]	1259
[3458]	1260	emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
[3190]	1261
	1262	!> SOX Compositions
	1263
	1264	ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO2") THEN
[3458]	1265	!> Kg/m2s kg/m2s
	1266	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)time_factor(icat)emt_att%sox_comp(icat,1)con_factorhour_per_day
[3190]	1267
[3458]	1268	emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
[3190]	1269
	1270	ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO4") THEN
[3458]	1271	!> Kg/m2s kg/m2s
	1272	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)time_factor(icat)emt_att%sox_comp(icat,2)con_factorhour_per_day
[3190]	1273
[3458]	1274	emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
[3190]	1275
	1276
	1277	!> PMs should be in kg/m**2/s, so conversion factor is here still required
	1278	!> PM1 Compositions
	1279	ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1") THEN
	1280
	1281	!> Cycle over the different pm components for PM1 type
	1282	DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,1))
	1283
[3458]	1284	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)time_factor(icat)emt_att%pm_comp(icat,i_pm_comp,1)con_factorhour_per_day
[3190]	1285
	1286
[3458]	1287	emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
[3190]	1288
	1289	ENDDO
	1290
	1291	!> PM2.5 Compositions
	1292	ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="PM25") THEN
	1293
	1294	!> Cycle over the different pm components for PM2.5 type
	1295	DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,2))
	1296
[3458]	1297	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)time_factor(icat)emt_att%pm_comp(icat,i_pm_comp,2)con_factorhour_per_day
[3190]	1298
	1299
[3458]	1300	emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
[3190]	1301
	1302	ENDDO
	1303
	1304	!> PM10 Compositions
	1305	ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
	1306
	1307	!> Cycle over the different pm components for PM10 type
	1308	DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,3))
	1309
[3458]	1310	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)time_factor(icat)emt_att%pm_comp(icat,i_pm_comp,3)con_factorhour_per_day
[3190]	1311
	1312
[3458]	1313	emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
[3190]	1314
	1315	ENDDO
	1316
	1317	!> VOCs Compositions: for VOCs, the input value is provided in kg/m*2s but the composition is provided in mole/kg: a conversion factor for the input that could be eventually provided in, for example, tons/(m*2s) is still required
	1318	ELSE IF (SIZE(match_spec_voc_input) .GT. 0) THEN
	1319
	1320	!TBD: maybe we can avoid the cycle
	1321	DO ivoc= 1,SIZE(match_spec_voc_input)
	1322
	1323	IF (TRIM(spc_names(match_spec_model(ispec)))==TRIM(emt_att%voc_name(ivoc))) THEN
	1324
[3458]	1325	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)time_factor(icat) &
	1326	emt_att%voc_comp(icat,match_spec_voc_input(ivoc))con_factorhour_per_day
[3190]	1327
[3458]	1328	emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
[3190]	1329
	1330	ENDIF
	1331
	1332	ENDDO
	1333
	1334	!> General case (other species)
	1335	ELSE
	1336
[3458]	1337	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)time_factor(icat)con_factor*hour_per_day
[3190]	1338
[3458]	1339	emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
[3190]	1340
	1341	ENDIF ! IF (spc_names==)
	1342
	1343	!> for every species and category set emis to 0 so to avoid overwriting. TBD: discuss whether necessary
	1344
	1345	emis(:,:)= 0
	1346
	1347	ENDDO
	1348
	1349	!> for every category set delta_emis to 0 so to avoid overwriting. TBD: discuss whether necessary
	1350
	1351	delta_emis(:,:)=0
	1352
	1353	ENDDO
	1354
	1355	ENDIF !> mode_emis
	1356
[3266]	1357	!-------------------------------------------------------------------------------------------------------
[3190]	1358	!--- Cycle to transform x,y coordinates to the one of surface_mod and to assign emission values to cssws
[3266]	1359	!-------------------------------------------------------------------------------------------------------
[3190]	1360
	1361	!-- PARAMETERIZED mode
	1362	!> For the PARAMETERIZED mode units of inputs are always micromoles/(m*2s). In this case we do not need the molar mass for conversion into ppms
	1363	IF (TRIM(mode_emis)=="PARAMETERIZED") THEN
	1364
	1365	IF ( street_type_f%from_file ) THEN
	1366
	1367	!> Streets are lsm surfaces, hence, no usm surface treatment required
[3467]	1368	IF (surf_lsm_h%ns .GT. 0 ) THEN
[3190]	1369	DO m = 1, surf_lsm_h%ns
	1370	i = surf_lsm_h%i(m)
	1371	j = surf_lsm_h%j(m)
	1372	k = surf_lsm_h%k(m)
	1373
	1374
	1375	IF ( street_type_f%var(j,i) >= main_street_id .AND. street_type_f%var(j,i) < max_street_id ) THEN
	1376
	1377	!> Cycle over already matched species
	1378	DO ispec=1,nspec_out
	1379
[3458]	1380	!> PMs are already in mass units: kilograms
	1381	IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" &
[3190]	1382	.OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
	1383
[3266]	1384	! kg/(m^2s) kg/m^3
[3190]	1385	surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec)) * &
[3276]	1386	! kg/(m^2*s)
[3266]	1387	emis_distribution(1,j,i,ispec)* &
	1388	! kg/m^3
	1389	rho_air(k)
[3190]	1390
	1391	ELSE
	1392
[3266]	1393	!> Other Species: inputs are micromoles: they have to be converted in moles
	1394	! ppm/s m kg/m^3
[3458]	1395	surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec))* &
[3266]	1396	! micromoles/(m^2*s)
[3458]	1397	emis_distribution(1,j,i,ispec) * &
	1398	! m^3/Nmole
	1399	conv_to_ratio(k,j,i)* &
[3266]	1400	! kg/m^3
	1401	rho_air(k)
[3190]	1402
[3266]	1403
[3190]	1404	ENDIF
	1405
	1406	ENDDO
	1407
	1408
[3458]	1409	ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. street_type_f%var(j,i) < main_street_id ) THEN
[3190]	1410
	1411	!> Cycle over already matched species
	1412	DO ispec=1,nspec_out
	1413
	1414	!> PMs are already in mass units: micrograms
[3458]	1415	IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" &
[3190]	1416	.OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
	1417
[3266]	1418	! kg/(m^2s) kg/m^3
[3458]	1419	surf_lsm_h%cssws(match_spec_model(ispec),m)= emiss_factor_side(match_spec_input(ispec)) * &
[3276]	1420	! kg/(m^2*s)
[3266]	1421	emis_distribution(1,j,i,ispec)* &
	1422	! kg/m^3
	1423	rho_air(k)
[3190]	1424	ELSE
	1425
	1426
	1427	!>Other Species: inputs are micromoles
[3266]	1428	! ppm/s m kg/m^3
[3458]	1429	surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * &
[3266]	1430	! micromoles/(m^2*s)
[3458]	1431	emis_distribution(1,j,i,ispec) * &
	1432	! m^3/Nmole
	1433	conv_to_ratio(k,j,i)* &
[3266]	1434	! kg/m^3
	1435	rho_air(k)
[3190]	1436	ENDIF
	1437
	1438	ENDDO
	1439
	1440	ELSE
	1441
	1442	!> If no street type is defined, then assign null emissions to all the species
	1443	surf_lsm_h%cssws(:,m) = 0.0_wp
	1444
	1445	ENDIF
	1446
	1447	ENDDO
[3467]	1448	ELSEIF ( salsa ) THEN
	1449	DO m = 1, surf_def_h(0)%ns
	1450	i = surf_def_h(0)%i(m)
	1451	j = surf_def_h(0)%j(m)
	1452	k = surf_def_h(0)%k(m)
[3190]	1453
[3467]	1454
	1455	IF ( street_type_f%var(j,i) >= main_street_id .AND. &
	1456	street_type_f%var(j,i) < max_street_id ) &
	1457	THEN
	1458
	1459	!> Cycle over already matched species
	1460	DO ispec=1,nspec_out
	1461
	1462	!> PMs are already in mass units:micrograms: have to be converted to kilograms
	1463	IF ( TRIM(spc_names(match_spec_model(ispec)))=="PM1" &
	1464	.OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" &
	1465	.OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10")&
	1466	THEN
	1467
	1468	surf_def_h(0)%cssws(match_spec_model(ispec),m) = &
	1469	emiss_factor_main(match_spec_input(ispec)) * &
	1470	emis_distribution(1,j,i,ispec) * rho_air(k) /&
	1471	time_factor(1)
	1472	ELSE
	1473
	1474	!> Other Species: inputs are micromoles: have to be converted
	1475	surf_def_h(0)%cssws(match_spec_model(ispec),m) = &
	1476	emiss_factor_main(match_spec_input(ispec)) * &
	1477	emis_distribution(1,j,i,ispec) * &
	1478	conv_to_ratio(k,j,i) * rho_air(k) / time_factor(1)
	1479	ENDIF
	1480	ENDDO
	1481
	1482	ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. &
	1483	street_type_f%var(j,i) < main_street_id ) &
	1484	THEN
	1485
	1486	!> Cycle over already matched species
	1487	DO ispec=1,nspec_out
	1488
	1489	!> PMs are already in mass units: micrograms
	1490	IF ( TRIM(spc_names(match_spec_model(ispec)))=="PM1" &
	1491	.OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" &
	1492	.OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10")&
	1493	THEN
	1494
	1495	surf_def_h(0)%cssws(match_spec_model(ispec),m) = &
	1496	emiss_factor_side(match_spec_input(ispec)) * &
	1497	emis_distribution(1,j,i,ispec) * rho_air(k) / &
	1498	time_factor(1)
	1499	ELSE
	1500
	1501	surf_def_h(0)%cssws(match_spec_model(ispec),m) = &
	1502	emiss_factor_side(match_spec_input(ispec)) * &
	1503	emis_distribution(1,j,i,ispec) * &
	1504	conv_to_ratio(k,j,i) * rho_air(k) / time_factor(1)
	1505	ENDIF
	1506
	1507	ENDDO
	1508
	1509	ELSE
	1510
	1511	!> If no street type is defined, then assign null emissions to all the species
	1512	surf_def_h(0)%cssws(:,m) = 0.0_wp
	1513
	1514	ENDIF
	1515
	1516	ENDDO
	1517
[3190]	1518	ENDIF
	1519
	1520	ENDIF
	1521
[3266]	1522	!> For both DEFAULT and PRE-PROCESSED
[3190]	1523	ELSE
	1524
	1525
[3266]	1526	DO ispec=1,nspec_out
	1527	!TBD: for the PRE-PROCESSED mode in the future, values at different heights should be included!
[3190]	1528	!> Default surface type
	1529	IF (surf_def_h(0)%ns .GT. 0) THEN
	1530
	1531	DO m = 1, surf_def_h(0)%ns
	1532
	1533	i = surf_def_h(0)%i(m)
	1534	j = surf_def_h(0)%j(m)
	1535
[3458]	1536	IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN
[3190]	1537
[3458]	1538
	1539	!> Distinguish between PMs (no needing conversion in ppms),
	1540	! VOC (already converted to moles/(m*2s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and
	1541	! other Species (still expressed in Kg/(m*2s) at this point)
	1542
	1543	!> PMs
	1544	IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" &
	1545	.OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
[3190]	1546
[3458]	1547	! kg/(m^2s) kg/m^3 kg/(m^2*s)
	1548	surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)* &
	1549	! kg/m^3
	1550	rho_air(nzb)
[3190]	1551
[3266]	1552
[3458]	1553	ELSE
[3190]	1554
[3458]	1555	!> VOCs
	1556	IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN
	1557	! ( ppm/s) * m * kg/m^3 mole/(m^2/s)
	1558	surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * &
	1559	! m^3/mole ppm
	1560	conv_to_ratio(nzb,j,i)ratio2ppm &
	1561	! kg/m^3
	1562	rho_air(nzb)
[3190]	1563
[3266]	1564
[3458]	1565	!> OTHER SPECIES
	1566	ELSE
[3190]	1567
[3458]	1568	! ( ppm/s )m kg/m^3 kg/(m^2/s)
	1569	surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * &
	1570	! mole/Kg
	1571	(1/emt_att%xm(ispec))* &
	1572	! m^3/mole ppm
	1573	conv_to_ratio(nzb,j,i)ratio2ppm &
	1574	! kg/m^3
	1575	rho_air(nzb)
[3266]	1576
[3190]	1577
[3458]	1578	ENDIF
	1579
[3190]	1580	ENDIF
	1581
	1582	ENDIF
	1583
	1584	ENDDO
	1585
	1586	END IF
	1587
	1588	!-- Land Surface Mode surface type
	1589	IF (surf_lsm_h%ns .GT. 0) THEN
	1590
	1591	DO m = 1, surf_lsm_h%ns
	1592
	1593	i = surf_lsm_h%i(m)
	1594	j = surf_lsm_h%j(m)
	1595	k = surf_lsm_h%k(m)
	1596
[3458]	1597	IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN
[3190]	1598
[3458]	1599	!> Distinguish between PMs (no needing conversion in ppms),
	1600	! VOC (already converted to moles/(m*2s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and
	1601	! other Species (still expressed in Kg/(m*2s) at this point)
[3190]	1602
[3458]	1603	!> PMs
	1604	IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" &
	1605	.OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
	1606
	1607	! kg/(m^2s) kg/m^3 kg/(m^2*s)
	1608	surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * &
	1609	! kg/m^3
	1610	rho_air(k)
[3190]	1611
[3458]	1612	ELSE
[3190]	1613
[3458]	1614	!> VOCs
	1615	IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN
	1616	! ( ppm/s) * m * kg/m^3 mole/(m^2/s)
	1617	surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * &
	1618	! m^3/mole ppm
	1619	conv_to_ratio(k,j,i)ratio2ppm &
	1620	! kg/m^3
	1621	rho_air(k)
[3190]	1622
[3458]	1623	!> OTHER SPECIES
	1624	ELSE
	1625
	1626	! ( ppm/s) * m * kg/m^3 kg/(m^2/s)
	1627	surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * &
	1628	! mole/Kg
	1629	(1/emt_att%xm(ispec))* &
	1630	! m^3/mole ppm
	1631	conv_to_ratio(k,j,i)ratio2ppm &
	1632	! kg/m^3
	1633	rho_air(k)
	1634
	1635	ENDIF
[3266]	1636
[3190]	1637	ENDIF
	1638
	1639	ENDIF
[3458]	1640
[3190]	1641	ENDDO
	1642
	1643	END IF
	1644
	1645	!-- Urban Surface Mode surface type
	1646	IF (surf_usm_h%ns .GT. 0) THEN
	1647
	1648
	1649	DO m = 1, surf_usm_h%ns
	1650
	1651	i = surf_usm_h%i(m)
	1652	j = surf_usm_h%j(m)
	1653	k = surf_usm_h%k(m)
	1654
[3458]	1655	IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN
[3190]	1656
[3458]	1657	!> PMs
	1658	IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" &
	1659	.OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
	1660
	1661	! kg/(m^2s) kg/m^3 kg/(m^2*s)
	1662	surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)* &
	1663	! kg/m^3
	1664	rho_air(k)
[3266]	1665
[3190]	1666
[3458]	1667	ELSE
[3190]	1668
[3458]	1669	!> VOCs
	1670	IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN
	1671	! ( ppm/s) * m * kg/m^3 mole/(m^2/s)
	1672	surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * &
	1673	! m^3/mole ppm
	1674	conv_to_ratio(k,j,i)ratio2ppm &
	1675	! kg/m^3
	1676	rho_air(k)
[3190]	1677
[3458]	1678	!> OTHER SPECIES
	1679	ELSE
[3190]	1680
	1681
[3458]	1682	! ( ppm/s) * m * kg/m^3 kg/(m^2/s)
	1683	surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * &
	1684	! mole/Kg
	1685	(1/emt_att%xm(ispec))* &
	1686	! m^3/mole ppm
	1687	conv_to_ratio(k,j,i)ratio2ppm &
	1688	! kg/m^3
	1689	rho_air(k)
[3266]	1690
	1691
[3458]	1692	ENDIF
	1693
[3190]	1694	ENDIF
	1695
	1696	ENDIF
	1697
	1698	ENDDO
	1699	END IF
	1700	ENDDO
	1701
	1702	ENDIF
	1703
	1704
	1705	!> At the end of each call to chem_emissions setup, the time_factor is deallocated (ALLOCATED ONLY in the DEFAULT mode)
	1706
	1707	IF ( ALLOCATED ( time_factor ) ) DEALLOCATE( time_factor )
	1708
	1709	ENDIF !> emis_output_required
	1710
	1711
	1712	END SUBROUTINE chem_emissions_setup
	1713
	1714	END MODULE chem_emissions_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |