Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |