Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |