Home

Context Navigation

chem_emissions_mod.f90 @ 4896

Last change on this file since 4896 was 4887, checked in by banzhafs, 4 years ago
Unnecessary comments removed from chemistry_model_mod and chem_emissions_mod
Property svn:keywords set to `Id`
File size: 94.4 KB

Rev	Line
[4055]	1	!> @file chem_emissions_mod.f90
[4559]	2	!--------------------------------------------------------------------------------------------------!
[4055]	3	! This file is part of PALM model system.
	4	!
[4559]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[4055]	8	!
[4559]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[4055]	12	!
[4559]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
[4055]	15	!
[4828]	16	! Copyright 2018-2021 Leibniz Universitaet Hannover
	17	! Copyright 2018-2021 Freie Universitaet Berlin
	18	! Copyright 2018-2021 Karlsruhe Institute of Technology
[4559]	19	!--------------------------------------------------------------------------------------------------!
[4055]	20	!
	21	! Current revisions:
	22	! ------------------
[4671]	23	!
	24	!
[4055]	25	! Former revisions:
	26	! -----------------
	27	! $Id: chem_emissions_mod.f90 4887 2021-02-26 16:22:32Z raasch $
[4887]	28	! Unnecessary comments removed
	29	!
	30	! 4828 2021-01-05 11:21:41Z Giersch
[4671]	31	! Implementation of downward facing USM and LSM surfaces
[4887]	32	!
[4671]	33	! 4559 2020-06-11 08:51:48Z raasch
[4559]	34	! file re-formatted to follow the PALM coding standard
	35	!
	36	! 4481 2020-03-31 18:55:54Z maronga
[4887]	37	! Implemented on-demand read mode for LOD 2 NB
[4403]	38	! - added following module global variables
	39	! - input_file_chem (namesake in netcdf_data_input_mod is local)
	40	! - timestamps
	41	! - added following public subroutines / interfaces
	42	! - chem_emisisons_header_init
	43	! - chem_emisisons_update_on_demand
	44	! - added following local subroutines
	45	! - chem_emisisons_header_init_lod2
	46	! - chem_emisisons_update_on_demand_lod2
	47	! - added following local auxiliary subroutines
	48	! - chem_emissions_init_species ( )
	49	! - chem_emissions_init_timestamps ( )
	50	! - chem_emissions_assign_surface_flux ( )
	51	! - added following local functions
	52	! - chem_emisisons_convert_base_units ( )
	53	! - chem_emissions_mass_2_molar_flux ( )
	54	! - chem_emissions_locate_species ( )
	55	! - chem_emissions_locate_timestep ( )
	56	! - added following error messages
	57	! - CM0468 - LOD mismatch (namelist / chemistry file)
	58	! - CM0469 - Timestamps no in choronological order
	59	! - depreciated unused module variable filename_emis
	60	!
	61	! 4356 2019-12-20 17:09:33Z suehring
[4356]	62	! Minor formatting adjustment
[4403]	63	!
[4356]	64	! 4242 2019-09-27 12:59:10Z suehring
[4559]	65	! Adjust index_hh access to new definition accompanied with new
[4242]	66	! palm_date_time_mod. Note, this is just a preleminary fix. (E Chan)
[4559]	67	!
[4242]	68	! 4230 2019-09-11 13:58:14Z suehring
[4559]	69	! Bugfix, consider that time_since_reference_point can be also negative when
[4230]	70	! time indices are determined.
[4559]	71	!
[4230]	72	! 4227 2019-09-10 18:04:34Z gronemeier
[4227]	73	! implement new palm_date_time_mod
[4559]	74	!
[4227]	75	! 4223 2019-09-10 09:20:47Z gronemeier
[4218]	76	! Unused routine chem_emissions_check_parameters commented out due to uninitialized content
[4559]	77	!
[4218]	78	! 4182 2019-08-22 15:20:23Z scharf
[4182]	79	! Corrected "Former revisions" section
[4559]	80	!
[4182]	81	! 4157 2019-08-14 09:19:12Z suehring
[4157]	82	! Replace global arrays also in mode_emis branch
[4559]	83	!
[4157]	84	! 4154 2019-08-13 13:35:59Z suehring
[4154]	85	! Replace global arrays for emissions by local ones.
[4559]	86	!
[4154]	87	! 4144 2019-08-06 09:11:47Z raasch
[4144]	88	! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
[4559]	89	!
[4144]	90	! 4055 2019-06-27 09:47:29Z suehring
[4559]	91	! - replaced local thermo. constants w/ module definitions in basic_constants_and_equations_mod
	92	! (rgas_univ, p_0, r_d_cp)
[4055]	93	! - initialize array emis_distribution immediately following allocation
[4559]	94	! - lots of minor formatting changes based on review sesson in 20190325 (E.C. Chan)
	95	!
[4055]	96	! 3968 2019-05-13 11:04:01Z suehring
[4559]	97	! - in subroutine chem_emissions_match replace all decision structures relating to mode_emis to
	98	! emiss_lod
[4055]	99	! - in subroutine chem_check_parameters replace emt%nspec with emt%n_emiss_species
	100	! - spring cleaning (E.C. Chan)
[4559]	101	!
[3968]	102	! 3885 2019-04-11 11:29:34Z kanani
[4559]	103	! Changes related to global restructuring of location messages and introduction of additional debug
	104	! messages
	105	!
[4055]	106	! 3831 2019-03-28 09:11:22Z forkel
[4559]	107	! added nvar to USE chem_gasphase_mod (chem_modules will not include nvar anymore)
	108	!
[4055]	109	! 3788 2019-03-07 11:40:09Z banzhafs
	110	! Removed unused variables from chem_emissions_mod
[4559]	111	!
	112	! 3772 2019-02-28 15:51:57Z suehring
	113	! - In case of parametrized emissions, assure that emissions are only on natural surfaces
	114	! (i.e. streets) and not on urban surfaces.
[4055]	115	! - some unnecessary if clauses removed
	116	!
	117	! 3685 2019 -01-21 01:02:11Z knoop
	118	! Some interface calls moved to module_interface + cleanup
[4182]	119	! 3286 2018-09-28 07:41:39Z forkel
	120	!
	121	! Authors:
	122	! --------
	123	! @author Emmanuele Russo (FU-Berlin)
	124	! @author Sabine Banzhaf (FU-Berlin)
	125	! @author Martijn Schaap (FU-Berlin, TNO Utrecht)
[4559]	126	!
[4055]	127	! Description:
	128	! ------------
	129	!> MODULE for reading-in Chemistry Emissions data
	130	!>
	131	!> @note <Enter notes on the module>
	132	!> @bug <Enter known bugs here>
[4559]	133	!--------------------------------------------------------------------------------------------------!
[4055]	134
	135	MODULE chem_emissions_mod
	136
[4559]	137	USE arrays_3d, &
[4055]	138	ONLY: rho_air
	139
[4559]	140	USE basic_constants_and_equations_mod, &
	141	ONLY: p_0, rd_d_cp, rgas_univ
[4055]	142
[4559]	143	USE control_parameters, &
	144	ONLY: debug_output, end_time, initializing_actions, intermediate_timestep_count, &
	145	message_string, dt_3d
	146
[4055]	147	USE indices
	148
	149	USE kinds
	150
	151	#if defined( __netcdf )
	152	USE netcdf
	153	#endif
	154
[4559]	155	USE netcdf_data_input_mod, &
[4055]	156	ONLY: chem_emis_att_type, chem_emis_val_type
	157
[4559]	158	USE chem_gasphase_mod, &
[4055]	159	ONLY: nvar, spc_names
[4559]	160
[4055]	161	USE chem_modules
	162
[4559]	163	USE statistics, &
[4055]	164	ONLY: weight_pres
	165
[4403]	166	!
[4887]	167	!-- 20200203 NB
[4403]	168	!-- Added new palm_date_time_mod for on-demand emission reading
	169
[4559]	170	USE palm_date_time_mod, &
[4403]	171	ONLY: get_date_time
	172
[4055]	173	IMPLICIT NONE
	174
	175	!
[4559]	176	!-- Declare all global variables within the module
[4403]	177
	178	!
[4887]	179	!-- 20200203 NB new variables for on-demand read mode
[4403]	180
	181	CHARACTER(LEN=*), PARAMETER :: input_file_chem = 'PIDS_CHEM' !< chemistry file
[4055]	182
[4559]	183	CHARACTER(LEN=512), ALLOCATABLE, DIMENSION(:) :: timestamps !< timestamps in chemistry file
	184
	185
	186	INTEGER(iwp) :: dt_emis !< Time Step Emissions
[4403]	187	INTEGER(iwp) :: i !< index 1st selected dimension (some dims are not spatial)
[4559]	188	INTEGER(iwp) :: j !< index 2nd selected dimension
	189	INTEGER(iwp) :: i_end !< Index to end read variable from netcdf in one dims
[4403]	190	INTEGER(iwp) :: i_start !< Index to start read variable from netcdf along one dims
[4559]	191	INTEGER(iwp) :: j_end !< Index to end read variable from netcdf in additional dims
[4403]	192	INTEGER(iwp) :: j_start !< Index to start read variable from netcdf in additional dims
	193	INTEGER(iwp) :: len_index !< length of index (used for several indices)
	194	INTEGER(iwp) :: len_index_pm !< length of PMs index
	195	INTEGER(iwp) :: len_index_voc !< length of voc index
	196	INTEGER(iwp) :: previous_timestamp_index !< index for current timestamp (20200203 ECC)
[4559]	197	INTEGER(iwp) :: z_end !< Index to end read variable from netcdf in additional dims
[4403]	198	INTEGER(iwp) :: z_start !< Index to start read variable from netcdf in additional dims
[4055]	199
	200	REAL(wp) :: conversion_factor !< Units Conversion Factor
	201
	202	SAVE
	203
	204	!
[4559]	205	!-- Matching Emissions actions
[4055]	206	INTERFACE chem_emissions_match
	207	MODULE PROCEDURE chem_emissions_match
	208	END INTERFACE chem_emissions_match
	209	!
[4559]	210	!-- Initialization actions
[4055]	211	INTERFACE chem_emissions_init
	212	MODULE PROCEDURE chem_emissions_init
	213	END INTERFACE chem_emissions_init
	214	!
	215	!-- Setup of Emissions
	216	INTERFACE chem_emissions_setup
	217	MODULE PROCEDURE chem_emissions_setup
	218	END INTERFACE chem_emissions_setup
[4403]	219
[4055]	220	!
[4887]	221	!-- 20200203 NB new interfaces for on-demand mode
[4403]	222
	223	!
[4559]	224	!-- initialization actions for on-demand mode
[4403]	225	INTERFACE chem_emissions_header_init
	226	MODULE PROCEDURE chem_emissions_header_init
	227	END INTERFACE chem_emissions_header_init
	228	!
	229	!-- load emission data for on-demand mode
	230	INTERFACE chem_emissions_update_on_demand
	231	MODULE PROCEDURE chem_emissions_update_on_demand
	232	END INTERFACE chem_emissions_update_on_demand
	233
	234	!
[4887]	235	!-- 20200203 NB update public routines
[4403]	236
[4559]	237	PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup, &
[4403]	238	chem_emissions_header_init, chem_emissions_update_on_demand
	239	!
[4055]	240	!-- Public Variables
	241	PUBLIC conversion_factor, len_index, len_index_pm, len_index_voc
	242
	243	CONTAINS
	244
	245
[4559]	246	!--------------------------------------------------------------------------------------------------!
[4055]	247	! Description:
	248	! ------------
[4559]	249	!> Matching the chemical species indices. The routine checks what are the indices of the emission
	250	!> input species and the corresponding ones of the model species. The routine gives as output a
	251	!> vector containing the number of common species: it is important to note that while the model
	252	!> species are distinct, their values could be given to a single species in input.
	253	!> For example, in the case of NO2 and NO, values may be passed in input as NOX values.
	254	!--------------------------------------------------------------------------------------------------!
[4055]	255
[4559]	256	SUBROUTINE chem_emissions_match( emt_att,len_index )
[4055]	257
[4559]	258	INTEGER(iwp) :: ind_inp !< Parameters for cycling through chemical input species
	259	INTEGER(iwp) :: ind_mod !< Parameters for cycling through chemical model species
[4055]	260	INTEGER(iwp) :: ind_voc !< Indices to check whether a split for voc should be done
	261	INTEGER(iwp) :: ispec !< index for cycle over effective number of emission species
	262	INTEGER(iwp) :: nspec_emis_inp !< Variable where to store # of emission species in input
	263
	264	INTEGER(iwp), INTENT(INOUT) :: len_index !< number of common species between input dataset & model species
	265
	266	TYPE(chem_emis_att_type), INTENT(INOUT) :: emt_att !< Chemistry Emission Array (decl. netcdf_data_input.f90)
	267
	268
	269	IF ( debug_output ) CALL debug_message( 'chem_emissions_match', 'start' )
	270
	271	!
	272	!-- Number of input emission species
	273
	274	nspec_emis_inp = emt_att%n_emiss_species
[4559]	275	! nspec_emis_inp=emt_att%nspec
[4055]	276
	277	!
	278	!-- Check the emission LOD: 0 (PARAMETERIZED), 1 (DEFAULT), 2 (PRE-PROCESSED)
	279	SELECT CASE (emiss_lod)
	280
	281	!
[4559]	282	!-- LOD 0 (PARAMETERIZED mode)
[4055]	283	CASE (0)
	284
	285	len_index = 0
[4559]	286	!
	287	!-- Number of species and number of matched species can be different but call is only made if
	288	!-- both are greater than zero.
[4055]	289	IF ( nvar > 0 .AND. nspec_emis_inp > 0 ) THEN
	290
	291	!
[4559]	292	!-- Cycle over model species
[4055]	293	DO ind_mod = 1, nvar
	294	ind_inp = 1
[4559]	295	DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' ) !< 'novalue' is the default
[4055]	296	IF ( TRIM( surface_csflux_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
	297	len_index = len_index + 1
	298	ENDIF
	299	ind_inp = ind_inp + 1
	300	ENDDO
	301	ENDDO
	302
	303	IF ( len_index > 0 ) THEN
	304
	305	!
[4559]	306	!-- Allocation of Arrays of the matched species
	307	ALLOCATE ( match_spec_input(len_index) )
[4055]	308	ALLOCATE ( match_spec_model(len_index) )
	309
	310	!
[4559]	311	!-- Pass species indices to declared arrays
[4055]	312	len_index = 0
	313
[4559]	314	DO ind_mod = 1, nvar
[4055]	315	ind_inp = 1
	316	DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )
[4559]	317	IF ( TRIM( surface_csflux_name(ind_inp) ) == TRIM(spc_names(ind_mod) ) ) &
	318	THEN
[4055]	319	len_index = len_index + 1
	320	match_spec_input(len_index) = ind_inp
	321	match_spec_model(len_index) = ind_mod
	322	ENDIF
	323	ind_inp = ind_inp + 1
	324	END DO
	325	END DO
	326
	327	!
[4559]	328	!-- Check
[4055]	329	DO ispec = 1, len_index
	330
[4559]	331	IF ( emiss_factor_main(match_spec_input(ispec) ) < 0 .AND. &
[4055]	332	emiss_factor_side(match_spec_input(ispec) ) < 0 ) THEN
	333
[4559]	334	message_string = 'PARAMETERIZED emissions mode selected:' // &
	335	' EMISSIONS POSSIBLE ONLY ON STREET SURFACES' // &
	336	' but values of scaling factors for street types' // &
	337	' emiss_factor_main AND emiss_factor_side' // &
	338	' not provided for each of the emissions species' // &
	339	' or not provided at all: PLEASE set a finite value' // &
	340	' for these parameters in the chemistry namelist'
[4055]	341	CALL message( 'chem_emissions_matching', 'CM0442', 2, 2, 0, 6, 0 )
[4559]	342
[4055]	343	ENDIF
	344
	345	END DO
	346
	347
	348	ELSE
[4559]	349
	350	message_string = 'Non of given Emission Species' // &
	351	' matches' // &
	352	' model chemical species' // &
	353	' Emission routine is not called'
[4055]	354	CALL message( 'chem_emissions_matching', 'CM0443', 0, 0, 0, 6, 0 )
[4559]	355
[4055]	356	ENDIF
	357
	358	ELSE
[4559]	359
	360	message_string = 'Array of Emission species not allocated: ' // &
	361	' Either no emission species are provided as input or' // &
	362	' no chemical species are used by PALM.' // &
	363	' Emission routine is not called'
	364	CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 )
	365
[4055]	366	ENDIF
	367
	368	!
[4559]	369	!-- LOD 1 (DEFAULT mode)
[4055]	370	CASE (1)
	371
[4559]	372	len_index = 0 ! total number of species (to be accumulated)
[4055]	373	len_index_voc = 0 ! total number of VOCs (to be accumulated)
	374	len_index_pm = 3 ! total number of PMs: PM1, PM2.5, PM10.
	375
	376	!
[4559]	377	!-- Number of model species and input species could be different but process this only when both are
	378	!-- non-zero
[4055]	379	IF ( nvar > 0 .AND. nspec_emis_inp > 0 ) THEN
	380
	381	!
[4559]	382	!-- Cycle over model species
[4055]	383	DO ind_mod = 1, nvar
	384
	385	!
[4559]	386	!-- Cycle over input species
[4055]	387	DO ind_inp = 1, nspec_emis_inp
	388
	389	!
[4559]	390	!-- Check for VOC Species
[4055]	391	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" ) THEN
	392	DO ind_voc= 1, emt_att%nvoc
[4559]	393
	394	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) ) &
	395	THEN
[4055]	396	len_index = len_index + 1
	397	len_index_voc = len_index_voc + 1
	398	ENDIF
[4559]	399
[4055]	400	END DO
	401	ENDIF
	402
	403	!
[4559]	404	!-- PMs: There is one input species name for all PM. This variable has 3 dimensions, one for PM1,
	405	!-- PM2.5 and PM10
[4055]	406	IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" ) THEN
	407
	408	IF ( TRIM( spc_names(ind_mod) ) == "PM1" ) THEN
	409	len_index = len_index + 1
	410	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" ) THEN
	411	len_index = len_index + 1
	412	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" ) THEN
	413	len_index = len_index + 1
	414	ENDIF
	415
	416	ENDIF
	417
	418	!
[4559]	419	!-- NOX: NO2 and NO
	420	IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOX" ) THEN
[4055]	421
	422	IF ( TRIM( spc_names(ind_mod) ) == "NO" ) THEN
	423	len_index = len_index + 1
	424	ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" ) THEN
	425	len_index = len_index + 1
	426	ENDIF
	427
	428	ENDIF
	429
	430	!
[4559]	431	!-- SOX: SO2 and SO4
[4055]	432	IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOX" ) THEN
	433
	434	IF ( TRIM( spc_names(ind_mod) ) == "SO2" ) THEN
	435	len_index = len_index + 1
	436	ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" ) THEN
	437	len_index = len_index + 1
	438	ENDIF
	439
	440	ENDIF
	441
	442	!
[4559]	443	!-- Other Species
	444	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
[4055]	445	len_index = len_index + 1
	446	ENDIF
	447
	448	END DO ! ind_inp ...
	449
	450	END DO ! ind_mod ...
	451
	452
	453	!
[4559]	454	!-- Allocate arrays
[4055]	455	IF ( len_index > 0 ) THEN
	456
	457	ALLOCATE ( match_spec_input(len_index) )
	458	ALLOCATE ( match_spec_model(len_index) )
	459
	460	IF ( len_index_voc > 0 ) THEN
	461	!
[4559]	462	!-- Contains indices of the VOC model species
[4055]	463	ALLOCATE( match_spec_voc_model(len_index_voc) )
	464	!
[4559]	465	!-- Contains the indices of different values of VOC composition of input variable
	466	!-- VOC_composition
[4055]	467	ALLOCATE( match_spec_voc_input(len_index_voc) )
	468
	469	ENDIF
	470
	471	!
[4559]	472	!-- Pass the species indices to declared arrays
[4055]	473	len_index = 0
	474	len_index_voc = 0
[4559]	475
[4055]	476	DO ind_mod = 1, nvar
[4559]	477	DO ind_inp = 1, nspec_emis_inp
[4055]	478
	479	!
[4559]	480	!-- VOCs
	481	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" .AND. &
	482	ALLOCATED( match_spec_voc_input ) ) THEN
[4055]	483
	484	DO ind_voc = 1, emt_att%nvoc
	485
[4559]	486	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )&
	487	THEN
[4055]	488
	489	len_index = len_index + 1
	490	len_index_voc = len_index_voc + 1
[4559]	491
[4055]	492	match_spec_input(len_index) = ind_inp
	493	match_spec_model(len_index) = ind_mod
	494
	495	match_spec_voc_input(len_index_voc) = ind_voc
	496	match_spec_voc_model(len_index_voc) = ind_mod
	497
	498	ENDIF
	499
	500	END DO
	501
	502	ENDIF
	503
	504	!
[4559]	505	!-- PMs
[4055]	506	IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" ) THEN
	507
	508	IF ( TRIM( spc_names(ind_mod) ) == "PM1" ) THEN
	509	len_index = len_index + 1
	510	match_spec_input(len_index) = ind_inp
	511	match_spec_model(len_index) = ind_mod
	512	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" ) THEN
	513	len_index = len_index + 1
	514	match_spec_input(len_index) = ind_inp
	515	match_spec_model(len_index) = ind_mod
	516	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" ) THEN
	517	len_index = len_index + 1
	518	match_spec_input(len_index) = ind_inp
	519	match_spec_model(len_index) = ind_mod
	520	ENDIF
	521
	522	ENDIF
	523
	524	!
[4559]	525	!-- NOX
[4055]	526	IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOX" ) THEN
	527
	528	IF ( TRIM( spc_names(ind_mod) ) == "NO" ) THEN
	529	len_index = len_index + 1
	530
	531	match_spec_input(len_index) = ind_inp
	532	match_spec_model(len_index) = ind_mod
	533
	534	ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" ) THEN
	535	len_index = len_index + 1
	536
	537	match_spec_input(len_index) = ind_inp
	538	match_spec_model(len_index) = ind_mod
[4559]	539
[4055]	540	ENDIF
	541
	542	ENDIF
	543
	544
	545	!
[4559]	546	!-- SOX
[4055]	547	IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOX" ) THEN
	548
	549	IF ( TRIM( spc_names(ind_mod) ) == "SO2" ) THEN
	550	len_index = len_index + 1
	551	match_spec_input(len_index) = ind_inp
	552	match_spec_model(len_index) = ind_mod
	553	ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" ) THEN
	554	len_index = len_index + 1
	555	match_spec_input(len_index) = ind_inp
	556	match_spec_model(len_index) = ind_mod
	557	ENDIF
	558
	559	ENDIF
	560
	561	!
[4559]	562	!-- Other Species
	563	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) &
	564	THEN
[4055]	565	len_index = len_index + 1
	566	match_spec_input(len_index) = ind_inp
	567	match_spec_model(len_index) = ind_mod
	568	ENDIF
	569
	570	END DO ! inp_ind
	571
	572	END DO ! inp_mod
	573
	574	!
[4559]	575	!-- Error reporting (no matching)
[4055]	576	ELSE
	577
[4559]	578	message_string = 'None of given Emission Species matches' // &
	579	' model chemical species' // &
	580	' Emission routine is not called'
	581	CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 )
[4055]	582
	583	ENDIF
	584
	585	!
[4559]	586	!-- Error reporting (no species)
[4055]	587	ELSE
	588
[4559]	589	message_string = 'Array of Emission species not allocated: ' // &
	590	' Either no emission species are provided as input or' // &
	591	' no chemical species are used by PALM:' // &
	592	' Emission routine is not called'
	593	CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 )
	594
[4055]	595	ENDIF
	596
	597	!
[4559]	598	!-- LOD 2 (PRE-PROCESSED mode)
[4055]	599	CASE (2)
	600
	601	len_index = 0
	602	len_index_voc = 0
	603
	604	IF ( nvar > 0 .AND. nspec_emis_inp > 0 ) THEN
	605	!
[4559]	606	!-- Cycle over model species
	607	DO ind_mod = 1, nvar
[4055]	608
	609	!
[4559]	610	!-- Cycle over input species
[4055]	611	DO ind_inp = 1, nspec_emis_inp
	612
	613	!
[4559]	614	!-- Check for VOC Species
	615	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" ) THEN
[4055]	616	DO ind_voc = 1, emt_att%nvoc
[4559]	617	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) ) &
	618	THEN
[4055]	619	len_index = len_index + 1
	620	len_index_voc = len_index_voc + 1
	621	ENDIF
	622	END DO
	623	ENDIF
	624
	625	!
[4559]	626	!-- Other Species
	627	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
[4055]	628	len_index = len_index + 1
	629	ENDIF
	630	ENDDO
	631	ENDDO
	632
	633	!
[4559]	634	!-- Allocate array for storing the indices of the matched species
	635	IF ( len_index > 0 ) THEN
[4055]	636
[4559]	637	ALLOCATE ( match_spec_input(len_index) )
	638
[4055]	639	ALLOCATE ( match_spec_model(len_index) )
	640
	641	IF ( len_index_voc > 0 ) THEN
	642	!
[4559]	643	!-- Contains indices of the VOC model species
[4055]	644	ALLOCATE( match_spec_voc_model(len_index_voc) )
	645	!
[4559]	646	!-- Contains the indices of different values of VOC composition of input variable
	647	!-- VOC_composition
[4055]	648	ALLOCATE( match_spec_voc_input(len_index_voc) )
	649
	650	ENDIF
	651
	652	!
[4559]	653	!-- Pass the species indices to declared arrays
[4055]	654	len_index = 0
	655
	656	!
[4559]	657	!-- Cycle over model species
	658	DO ind_mod = 1, nvar
[4055]	659
	660	!
[4559]	661	!-- Cycle over Input species
[4055]	662	DO ind_inp = 1, nspec_emis_inp
	663
	664	!
[4559]	665	!-- VOCs
	666	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" .AND. &
	667	ALLOCATED( match_spec_voc_input ) ) THEN
[4055]	668
	669	DO ind_voc= 1, emt_att%nvoc
[4559]	670	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )&
	671	THEN
[4055]	672	len_index = len_index + 1
	673	len_index_voc = len_index_voc + 1
[4559]	674
[4055]	675	match_spec_input(len_index) = ind_inp
	676	match_spec_model(len_index) = ind_mod
	677
	678	match_spec_voc_input(len_index_voc) = ind_voc
[4559]	679	match_spec_voc_model(len_index_voc) = ind_mod
[4055]	680	ENDIF
	681	END DO
	682	ENDIF
	683
	684	!
[4559]	685	!-- Other Species
	686	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) &
	687	THEN
[4055]	688	len_index = len_index + 1
	689	match_spec_input(len_index) = ind_inp
	690	match_spec_model(len_index) = ind_mod
	691	ENDIF
	692
	693	END DO ! ind_inp
	694	END DO ! ind_mod
	695
[4144]	696	ELSE ! if len_index_voc <= 0
[4055]	697
	698	!
[4559]	699	!-- In case there are no species matching (just informational message)
	700	message_string = 'Non of given emission species' // &
	701	' matches' // &
	702	' model chemical species:' // &
	703	' Emission routine is not called'
[4055]	704	CALL message( 'chem_emissions_matching', 'CM0438', 0, 0, 0, 6, 0 )
	705	ENDIF
	706
	707	!
[4559]	708	!-- Error check (no matching)
[4055]	709	ELSE
	710
	711	!
[4559]	712	!-- Either spc_names is zero or nspec_emis_inp is not allocated
	713	message_string = 'Array of Emission species not allocated:' // &
	714	' Either no emission species are provided as input or' // &
	715	' no chemical species are used by PALM:' // &
	716	' Emission routine is not called'
	717	CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 )
[4055]	718
[4559]	719	ENDIF
[4055]	720
	721	!
	722	!-- If emission module is switched on but mode_emis is not specified or it is given the wrong name
	723
	724	!
[4559]	725	!-- Error check (no species)
[4055]	726	CASE DEFAULT
	727
	728	message_string = 'Emission Module switched ON, but' // &
	729	' either no emission mode specified or incorrectly given :' // &
[4559]	730	' please, pass the correct value to the namelist parameter "mode_emis"'
	731	CALL message( 'chem_emissions_matching', 'CM0445', 2, 2, 0, 6, 0 )
[4055]	732
	733	END SELECT
	734
	735	IF ( debug_output ) CALL debug_message( 'chem_emissions_match', 'end' )
	736
	737	END SUBROUTINE chem_emissions_match
	738
[4559]	739
	740	!--------------------------------------------------------------------------------------------------!
[4055]	741	! Description:
	742	! ------------
	743	!> Initialization:
[4559]	744	!> Netcdf reading, arrays allocation and first calculation of cssws fluxes at timestep 0
	745	!--------------------------------------------------------------------------------------------------!
[4055]	746
	747	SUBROUTINE chem_emissions_init
	748
[4559]	749	USE netcdf_data_input_mod, &
[4055]	750	ONLY: chem_emis, chem_emis_att
[4559]	751
[4055]	752	IMPLICIT NONE
[4559]	753
[4055]	754	INTEGER(iwp) :: ispec !< running index
	755
	756
	757	IF ( debug_output ) CALL debug_message( 'chem_emissions_init', 'start' )
	758
	759	!
	760	!-- Matching
[4559]	761	CALL chem_emissions_match( chem_emis_att, n_matched_vars )
[4055]	762
	763	IF ( n_matched_vars == 0 ) THEN
[4559]	764
[4055]	765	emission_output_required = .FALSE.
	766
	767	ELSE
	768
	769	emission_output_required = .TRUE.
	770
	771
	772	!
[4559]	773	!-- Set molecule masses (in kg/mol)
[4055]	774	ALLOCATE( chem_emis_att%xm(n_matched_vars) )
	775
	776	DO ispec = 1, n_matched_vars
	777	SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
	778	CASE ( 'SO2' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 2
	779	CASE ( 'SO4' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 4
	780	CASE ( 'NO' ); chem_emis_att%xm(ispec) = xm_N + xm_O
	781	CASE ( 'NO2' ); chem_emis_att%xm(ispec) = xm_N + xm_O * 2
	782	CASE ( 'NH3' ); chem_emis_att%xm(ispec) = xm_N + xm_H * 3
	783	CASE ( 'CO' ); chem_emis_att%xm(ispec) = xm_C + xm_O
	784	CASE ( 'CO2' ); chem_emis_att%xm(ispec) = xm_C + xm_O * 2
	785	CASE ( 'CH4' ); chem_emis_att%xm(ispec) = xm_C + xm_H * 4
	786	CASE ( 'HNO3' ); chem_emis_att%xm(ispec) = xm_H + xm_N + xm_O*3
	787	CASE DEFAULT
	788	chem_emis_att%xm(ispec) = 1.0_wp
	789	END SELECT
	790	ENDDO
	791
[4559]	792
[4055]	793	!
[4559]	794	!-- Get emissions for the first time step base on LOD (if defined) or emission mode
	795	!-- (if no LOD defined)
[4055]	796
	797	!
[4559]	798	!-- NOTE - I could use a combined if ( lod = xxx .or. mode = 'XXX' ) type of decision structure but
	799	! I think it is much better to implement it this way (i.e., conditional on lod if it is
	800	! defined, and mode if not) as we can easily take out the case structure for mode_emis
	801	! later on.
[4055]	802
	803	IF ( emiss_lod < 0 ) THEN !-- no LOD defined (not likely)
	804
[4559]	805	SELECT CASE ( TRIM( mode_emis ) )
[4055]	806
	807	CASE ( 'PARAMETERIZED' ) ! LOD 0
	808
[4559]	809	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
[4154]	810	ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
[4055]	811	ENDIF
	812
	813	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars)
	814
	815	CASE ( 'DEFAULT' ) ! LOD 1
	816
[4559]	817	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
[4154]	818	ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
[4055]	819	ENDIF
	820
	821	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
	822
	823	CASE ( 'PRE-PROCESSED' ) ! LOD 2
	824
[4559]	825	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
[4154]	826	!
[4559]	827	!-- Note, at the moment emissions are considered only by surface fluxes rather than
	828	!-- by volume sources. Therefore, no vertical dimension is required and is thus
	829	!-- allocated with 1. Later when volume sources are considered, the vertical
	830	!-- dimension will increase.
[4154]	831	!ALLOCATE( emis_distribution(nzb:nzt+1,nys:nyn,nxl:nxr,n_matched_vars) )
	832	ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
[4055]	833	ENDIF
[4559]	834
[4055]	835	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
	836
	837	END SELECT
	838
	839	ELSE ! if LOD is defined
	840
	841	SELECT CASE ( emiss_lod )
	842
	843	CASE ( 0 ) ! parameterized mode
	844
[4559]	845	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
[4157]	846	ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
[4055]	847	ENDIF
	848
	849	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars)
	850
	851	CASE ( 1 ) ! default mode
	852
[4559]	853	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
[4157]	854	ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
[4055]	855	ENDIF
	856
	857	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
	858
	859	CASE ( 2 ) ! pre-processed mode
	860
[4559]	861	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
[4157]	862	ALLOCATE( emis_distribution(nzb:nzt+1,nys:nyn,nxl:nxr,n_matched_vars) )
[4055]	863	ENDIF
[4559]	864
[4055]	865	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
	866
	867	END SELECT
	868
	869	ENDIF
	870
	871	!
[4559]	872	! -- Initialize
[4055]	873	emis_distribution = 0.0_wp
	874
	875	ENDIF
	876
	877	IF ( debug_output ) CALL debug_message( 'chem_emissions_init', 'end' )
	878
	879	END SUBROUTINE chem_emissions_init
	880
	881
	882
[4559]	883	!--------------------------------------------------------------------------------------------------!
[4055]	884	! Description:
	885	! ------------
[4227]	886	!> Routine for Update of Emission values at each timestep.
	887	!>
[4559]	888	!> @todo Clarify the correct usage of index_dd, index_hh and index_mm. Consider renaming of these
	889	!> variables.
[4227]	890	!> @todo Clarify time used in emis_lod=2 mode. ATM, the used time seems strange.
[4559]	891	!--------------------------------------------------------------------------------------------------!
[4055]	892
	893	SUBROUTINE chem_emissions_setup( emt_att, emt, n_matched_vars )
[4559]	894
	895	USE surface_mod, &
[4055]	896	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
	897
[4559]	898	USE netcdf_data_input_mod, &
[4055]	899	ONLY: street_type_f
	900
[4559]	901	USE arrays_3d, &
	902	ONLY: hyp, pt
[4055]	903
[4559]	904	USE control_parameters, &
[4227]	905	ONLY: time_since_reference_point
	906
[4559]	907	USE palm_date_time_mod, &
[4227]	908	ONLY: days_per_week, get_date_time, hours_per_day, months_per_year, seconds_per_day
[4559]	909
[4356]	910	IMPLICIT NONE
[4055]	911
[4559]	912	INTEGER(iwp) :: day_of_month !< day of the month
	913	INTEGER(iwp) :: day_of_week !< day of the week
	914	INTEGER(iwp) :: day_of_year !< day of the year
	915	INTEGER(iwp) :: days_since_reference_point !< days since reference point
	916	INTEGER(iwp) :: i !< running index for grid in x-direction
	917	INTEGER(iwp) :: i_pm_comp !< index for number of PM components
	918	INTEGER(iwp) :: icat !< Index for number of categories
	919	INTEGER(iwp) :: index_dd !< index day
	920	INTEGER(iwp) :: index_hh !< index hour
	921	INTEGER(iwp) :: index_mm !< index month
	922	INTEGER(iwp) :: ispec !< index for number of species
	923	INTEGER(iwp) :: ivoc !< Index for number of VOCs
	924	INTEGER(iwp) :: hour_of_day !< hour of the day
	925	INTEGER(iwp) :: j !< running index for grid in y-direction
	926	INTEGER(iwp) :: k !< running index for grid in z-direction
	927	INTEGER(iwp) :: m !< running index for horizontal surfaces
	928	INTEGER(iwp) :: month_of_year !< month of the year
[4055]	929
[4559]	930	INTEGER,INTENT(IN) :: n_matched_vars !< Output of matching routine with number
	931	!< of matched species
[4055]	932
[4559]	933	REAL(wp) :: time_utc_init !< second of day of initial date
[4055]	934
[4559]	935	TYPE(chem_emis_att_type), INTENT(INOUT) :: emt_att !< variable to store emission information
[4055]	936
[4559]	937	TYPE(chem_emis_val_type), INTENT(INOUT), ALLOCATABLE, DIMENSION(:) :: emt !< variable to store emission input values,
	938	!< depending on the considered species
	939	!
	940	!-- CONVERSION FACTORS: TIME
	941	REAL(wp), PARAMETER :: hour_per_year = 8760.0_wp !< number of hours in a year of 365 days
	942	REAL(wp), PARAMETER :: s_per_hour = 3600.0_wp !< number of sec per hour (s)/(hour)
	943	REAL(wp), PARAMETER :: s_per_day = 86400.0_wp !< number of sec per day (s)/(day)
[4227]	944
[4055]	945	REAL(wp), PARAMETER :: day_to_s = 1.0_wp/s_per_day !< conversion day -> sec
	946	REAL(wp), PARAMETER :: hour_to_s = 1.0_wp/s_per_hour !< conversion hours -> sec
	947	REAL(wp), PARAMETER :: year_to_s = 1.0_wp/(s_per_hour*hour_per_year) !< conversion year -> sec
[4559]	948
	949
	950	!
	951	!-- CONVERSION FACTORS: MASS
	952	REAL(wp), PARAMETER :: g_to_kg = 1.0E-03_wp !< Conversion from g to kg (kg/g)
[4055]	953	REAL(wp), PARAMETER :: miug_to_kg = 1.0E-09_wp !< Conversion from g to kg (kg/g)
[4559]	954	REAL(wp), PARAMETER :: tons_to_kg = 100.0_wp !< Conversion from tons to kg (kg/tons)
	955	!
	956	!-- CONVERSION FACTORS: PPM
	957	REAL(wp), PARAMETER :: ratio2ppm = 1.0E+06_wp
	958
[4055]	959	REAL(wp), DIMENSION(24) :: par_emis_time_factor !< time factors for the parameterized mode:
	960	!< fixed houlry profile for example day
[4559]	961	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: conv_to_ratio !< factor used for converting input
[4055]	962	!< to concentration ratio
	963	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: tmp_temp !< temporary variable for abs. temperature
	964
	965	REAL(wp), DIMENSION(:), ALLOCATABLE :: time_factor !< factor for time scaling of emissions
[4559]	966
	967	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: delta_emis !< incremental emission factor
[4055]	968	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: emis !< emission factor
	969
[4559]	970
[4055]	971	IF ( emission_output_required ) THEN
	972
	973	!
[4559]	974	!-- Set emis_dt to be used - since chemistry ODEs can be stiff, the option to solve them at every
	975	!-- RK substep is present to help improve stability should the need arises
	976
[4055]	977	IF ( call_chem_at_all_substeps ) THEN
	978
	979	dt_emis = dt_3d * weight_pres(intermediate_timestep_count)
	980
	981	ELSE
	982
	983	dt_emis = dt_3d
	984
	985	ENDIF
	986
	987	!
[4559]	988	!-- Conversion of units to the ones employed in PALM
	989	!-- In PARAMETERIZED mode no conversion is performed: in this case input units are fixed
[4055]	990	IF ( TRIM( mode_emis ) == "DEFAULT" .OR. TRIM( mode_emis ) == "PRE-PROCESSED" ) THEN
	991
	992	SELECT CASE ( TRIM( emt_att%units ) )
	993
	994	CASE ( 'kg/m2/s', 'KG/M2/S' ); conversion_factor = 1.0_wp ! kg
	995	CASE ( 'kg/m2/hour', 'KG/M2/HOUR' ); conversion_factor = hour_to_s
	996	CASE ( 'kg/m2/day', 'KG/M2/DAY' ); conversion_factor = day_to_s
	997	CASE ( 'kg/m2/year', 'KG/M2/YEAR' ); conversion_factor = year_to_s
	998
	999	CASE ( 'ton/m2/s', 'TON/M2/S' ); conversion_factor = tons_to_kg ! tonnes
	1000	CASE ( 'ton/m2/hour', 'TON/M2/HOUR' ); conversion_factor = tons_to_kg*hour_to_s
	1001	CASE ( 'ton/m2/year', 'TON/M2/YEAR' ); conversion_factor = tons_to_kg*year_to_s
	1002
	1003	CASE ( 'g/m2/s', 'G/M2/S' ); conversion_factor = g_to_kg ! grams
	1004	CASE ( 'g/m2/hour', 'G/M2/HOUR' ); conversion_factor = g_to_kg*hour_to_s
	1005	CASE ( 'g/m2/year', 'G/M2/YEAR' ); conversion_factor = g_to_kg*year_to_s
	1006
	1007	CASE ( 'micrograms/m2/s', 'MICROGRAMS/M2/S' ); conversion_factor = miug_to_kg ! ug
	1008	CASE ( 'micrograms/m2/hour', 'MICROGRAMS/M2/HOUR' ); conversion_factor = miug_to_kg*hour_to_s
	1009	CASE ( 'micrograms/m2/year', 'MICROGRAMS/M2/YEAR' ); conversion_factor = miug_to_kg*year_to_s
	1010
	1011	!
[4559]	1012	!-- Error check (need units)
	1013	CASE DEFAULT
	1014	message_string = 'The Units of the provided emission input' // &
	1015	' are not the ones required by PALM-4U: please check ' // &
	1016	' emission module documentation.'
[4055]	1017	CALL message( 'chem_emissions_setup', 'CM0446', 2, 2, 0, 6, 0 )
	1018
	1019	END SELECT
	1020
	1021	ENDIF
	1022
	1023	!
[4559]	1024	!-- Conversion factor to convert kg/m**2/s to ppm/s
[4055]	1025	DO i = nxl, nxr
	1026	DO j = nys, nyn
	1027
	1028	!
[4559]	1029	!-- Derive Temperature from Potential Temperature
	1030	tmp_temp(nzb:nzt+1,j,i) = pt(nzb:nzt+1,j,i) * ( hyp(nzb:nzt+1) / p_0 )**rd_d_cp
[4055]	1031
[4559]	1032	!
	1033	!-- We need to pass to cssws <- (ppm/s) * dz
	1034	!-- Input is Nmole/(m^2*s)
	1035	!-- To go to ppmdz multiply the input by (m2/N)dz
	1036	!-- (m*2/N)dz == V/N
	1037	!-- V/N = RT/P
	1038	conv_to_ratio(nzb:nzt+1,j,i) = rgas_univ * & ! J K-1 mol-1
	1039	tmp_temp(nzb:nzt+1,j,i) / & ! K
	1040	hyp(nzb:nzt+1) ! Pa
[4055]	1041
	1042	ENDDO
	1043	ENDDO
	1044
	1045
	1046	!
	1047	!-- LOD 2 (PRE-PROCESSED MODE)
	1048
	1049	IF ( emiss_lod == 2 ) THEN
	1050
	1051	!
[4559]	1052	!-- Update time indices
[4227]	1053	CALL get_date_time( 0.0_wp, second_of_day=time_utc_init )
[4559]	1054	CALL get_date_time( MAX( 0.0_wp, time_since_reference_point ), hour=hour_of_day )
[4055]	1055
[4559]	1056	days_since_reference_point = INT( FLOOR( ( time_utc_init + &
	1057	MAX( 0.0_wp, time_since_reference_point ) ) &
	1058	/ seconds_per_day ) )
[4055]	1059
[4227]	1060	index_hh = days_since_reference_point * hours_per_day + hour_of_day
	1061
[4055]	1062	!
[4559]	1063	!-- LOD 1 (DEFAULT MODE)
[4055]	1064	ELSEIF ( emiss_lod == 1 ) THEN
	1065
	1066	!
[4559]	1067	!-- Allocate array where to store temporary emission values
	1068	IF ( .NOT. ALLOCATED(emis) ) ALLOCATE( emis(nys:nyn,nxl:nxr) )
[4055]	1069
	1070	!
[4559]	1071	!-- Allocate time factor per category
[4055]	1072	ALLOCATE( time_factor(emt_att%ncat) )
	1073
	1074	!
[4559]	1075	!-- Read-in hourly emission time factor
	1076	IF ( TRIM( time_fac_type ) == "HOUR" ) THEN
[4055]	1077
	1078	!
[4559]	1079	!-- Update time indices
	1080	CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ), &
[4227]	1081	day_of_year=day_of_year, hour=hour_of_day )
	1082	index_hh = ( day_of_year - 1_iwp ) * hour_of_day
[4055]	1083
	1084	!
[4559]	1085	!-- Check if the index is less or equal to the temporal dimension of HOURLY emission files
[4055]	1086	IF ( index_hh <= SIZE( emt_att%hourly_emis_time_factor(1,:) ) ) THEN
	1087
	1088	!
[4559]	1089	!-- Read-in the correspondant time factor
	1090	time_factor(:) = emt_att%hourly_emis_time_factor(:,index_hh+1)
[4055]	1091
	1092	!
[4559]	1093	!-- Error check (time out of range)
[4055]	1094	ELSE
	1095
[4559]	1096	message_string = 'The "HOUR" time-factors in the DEFAULT mode ' // &
	1097	' are not provided for each hour of the total simulation time'
	1098	CALL message( 'chem_emissions_setup', 'CM0448', 2, 2, 0, 6, 0 )
[4055]	1099
	1100	ENDIF
	1101
	1102	!
[4559]	1103	!-- Read-in MDH emissions time factors
[4055]	1104	ELSEIF ( TRIM( time_fac_type ) == "MDH" ) THEN
	1105
	1106	!
[4559]	1107	!-- Update time indices
	1108	CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ), &
	1109	month = month_of_year, &
	1110	day = day_of_month, &
	1111	hour = hour_of_day, &
	1112	day_of_week = day_of_week &
	1113	)
[4227]	1114	index_mm = month_of_year
	1115	index_dd = months_per_year + day_of_week
[4559]	1116	SELECT CASE( TRIM( daytype_mdh ) )
[4055]	1117
[4227]	1118	CASE ("workday")
	1119	index_hh = months_per_year + days_per_week + hour_of_day
	1120
	1121	CASE ("weekend")
	1122	index_hh = months_per_year + days_per_week + hours_per_day + hour_of_day
	1123
	1124	CASE ("holiday")
	1125	index_hh = months_per_year + days_per_week + 2*hours_per_day + hour_of_day
	1126
	1127	END SELECT
[4055]	1128
	1129	!
[4559]	1130	!-- Check if the index is less or equal to the temporal dimension of MDH emission files
	1131	IF ( ( index_hh + index_dd + index_mm) <= SIZE( emt_att%mdh_emis_time_factor(1,:) ) )&
	1132	THEN
	1133	!
	1134	!-- Read in corresponding time factor
	1135	time_factor(:) = emt_att%mdh_emis_time_factor(:,index_mm) * &
	1136	emt_att%mdh_emis_time_factor(:,index_dd) * &
[4242]	1137	emt_att%mdh_emis_time_factor(:,index_hh+1)
[4055]	1138
	1139	!
[4559]	1140	!-- Error check (MDH time factor not provided)
[4055]	1141	ELSE
	1142
[4559]	1143	message_string = 'The "MDH" time-factors in the DEFAULT mode ' // &
	1144	' are not provided for each hour/day/month of the total simulation time'
[4055]	1145	CALL message( 'chem_emissions_setup', 'CM0449', 2, 2, 0, 6, 0 )
	1146
[4559]	1147	ENDIF
[4055]	1148
	1149	!
[4559]	1150	!-- Error check (no time factor defined)
	1151	ELSE
[4055]	1152
[4559]	1153	message_string = 'In the DEFAULT mode the time factor' // &
	1154	' has to be defined in the NAMELIST'
[4055]	1155	CALL message( 'chem_emissions_setup', 'CM0450', 2, 2, 0, 6, 0 )
[4559]	1156
[4055]	1157	ENDIF
	1158
	1159	!
[4559]	1160	!-- PARAMETERIZED MODE
[4055]	1161	ELSEIF ( emiss_lod == 0 ) THEN
	1162
	1163
	1164	!
[4559]	1165	!-- Assign constant values of time factors, diurnal profile for traffic sector
	1166	par_emis_time_factor(:) = (/ 0.009, 0.004, 0.004, 0.009, 0.029, 0.039, &
	1167	0.056, 0.053, 0.051, 0.051, 0.052, 0.055, &
	1168	0.059, 0.061, 0.064, 0.067, 0.069, 0.069, &
[4055]	1169	0.049, 0.039, 0.039, 0.029, 0.024, 0.019 /)
	1170
[4559]	1171	IF ( .NOT. ALLOCATED( time_factor ) ) ALLOCATE( time_factor(1) )
	1172
[4055]	1173	!
[4227]	1174	!-- Get time-factor for specific hour
[4559]	1175	CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ), hour = hour_of_day )
[4055]	1176
	1177	index_hh = hour_of_day
[4242]	1178	time_factor(1) = par_emis_time_factor(index_hh+1)
[4055]	1179
	1180	ENDIF ! emiss_lod
	1181
[4227]	1182
[4055]	1183	!
	1184	!-- Emission distribution calculation
	1185
	1186	IF ( emiss_lod == 0 ) THEN
	1187
	1188	DO ispec = 1, n_matched_vars
	1189
	1190	!
	1191	!-- Units are micromoles/m*2day (or kilograms/m*2day for PMs)
[4559]	1192	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = surface_csflux( match_spec_input(ispec) )&
	1193	* time_factor(1) * hour_to_s
[4055]	1194
[4559]	1195	ENDDO
[4055]	1196
	1197
	1198	!
[4559]	1199	!-- LOD 1 (DEFAULT mode)
[4055]	1200	ELSEIF ( emiss_lod == 1 ) THEN
	1201
	1202	!
[4559]	1203	!-- Allocate array for the emission value corresponding to a specific category and time factor
[4055]	1204	ALLOCATE (delta_emis(nys:nyn,nxl:nxr))
	1205
	1206	!
[4559]	1207	!-- Cycle over categories
	1208	DO icat = 1, emt_att%ncat
[4055]	1209
	1210	!
[4559]	1211	!-- Cycle over Species: n_matched_vars represents the number of species in common between
	1212	!-- the emission input data and the chemistry mechanism used
[4055]	1213	DO ispec = 1, n_matched_vars
	1214
[4559]	1215	emis(nys:nyn,nxl:nxr) = emt( match_spec_input(ispec) )% &
	1216	default_emission_data(icat,nys+1:nyn+1,nxl+1:nxr+1)
[4055]	1217
	1218	!
[4559]	1219	!-- NO
	1220	IF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "NO" ) THEN
[4055]	1221
[4559]	1222	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
	1223	time_factor(icat) * &
	1224	emt_att%nox_comp(icat,1) * &
[4227]	1225	conversion_factor * hours_per_day
[4055]	1226
[4559]	1227	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
	1228	emis_distribution(1,nys:nyn,nxl:nxr,ispec)&
	1229	+ delta_emis(nys:nyn,nxl:nxr)
[4055]	1230	!
[4559]	1231	!-- NO2
	1232	ELSEIF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "NO2" ) THEN
[4055]	1233
[4559]	1234	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
	1235	time_factor(icat) * &
	1236	emt_att%nox_comp(icat,2) * &
[4227]	1237	conversion_factor * hours_per_day
[4055]	1238
[4559]	1239	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
	1240	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
	1241	+ delta_emis(nys:nyn,nxl:nxr)
	1242
[4055]	1243	!
[4559]	1244	!-- SO2
	1245	ELSEIF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "SO2" ) THEN
	1246
	1247	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
	1248	time_factor(icat) * &
	1249	emt_att%sox_comp(icat,1) * &
[4227]	1250	conversion_factor * hours_per_day
[4055]	1251
[4559]	1252	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
	1253	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
	1254	+ delta_emis(nys:nyn,nxl:nxr)
[4055]	1255
	1256	!
[4559]	1257	!-- SO4
	1258	ELSEIF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "SO4" ) THEN
[4055]	1259
	1260
[4559]	1261	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
	1262	time_factor(icat) * &
	1263	emt_att%sox_comp(icat,2) * &
[4227]	1264	conversion_factor * hours_per_day
[4055]	1265
[4559]	1266	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
	1267	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
	1268	+ delta_emis(nys:nyn,nxl:nxr)
[4055]	1269
[4559]	1270
[4055]	1271	!
[4559]	1272	!-- PM1
	1273	ELSEIF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM1" ) THEN
[4055]	1274
	1275	DO i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,1) ) ! cycle through components
	1276
[4559]	1277	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
	1278	time_factor(icat) * &
	1279	emt_att%pm_comp(icat,i_pm_comp,1) * &
	1280	conversion_factor * hours_per_day
[4055]	1281
[4559]	1282	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
	1283	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
	1284	+ delta_emis(nys:nyn,nxl:nxr)
[4055]	1285
	1286	ENDDO
	1287
	1288	!
[4559]	1289	!-- PM2.5
	1290	ELSEIF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM25" ) THEN
[4055]	1291
	1292	DO i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,2) ) ! cycle through components
	1293
[4559]	1294	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
	1295	time_factor(icat) * &
	1296	emt_att%pm_comp(icat,i_pm_comp,2) * &
	1297	conversion_factor * hours_per_day
[4055]	1298
[4559]	1299	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
	1300	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
	1301	+ delta_emis(nys:nyn,nxl:nxr)
	1302
[4055]	1303	ENDDO
	1304
	1305	!
[4559]	1306	!-- PM10
	1307	ELSEIF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM10" ) THEN
[4055]	1308
	1309	DO i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,3) ) ! cycle through components
	1310
[4559]	1311	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
	1312	time_factor(icat) * &
	1313	emt_att%pm_comp(icat,i_pm_comp,3) * &
	1314	conversion_factor * hours_per_day
[4055]	1315
[4559]	1316	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
	1317	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
	1318	+ delta_emis(nys:nyn,nxl:nxr)
[4055]	1319
	1320	ENDDO
	1321
	1322	!
[4559]	1323	!-- VOCs
[4055]	1324	ELSEIF ( SIZE( match_spec_voc_input ) > 0 ) THEN
	1325
	1326	DO ivoc = 1, SIZE( match_spec_voc_input ) ! cycle through components
	1327
[4559]	1328	IF ( TRIM( spc_names(match_spec_model(ispec) ) ) == &
	1329	TRIM( emt_att%voc_name(ivoc) ) ) THEN
[4055]	1330
[4559]	1331	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * &
	1332	time_factor(icat) * &
	1333	emt_att%voc_comp(icat,match_spec_voc_input(ivoc)) * &
[4227]	1334	conversion_factor * hours_per_day
[4055]	1335
[4559]	1336	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
	1337	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
	1338	+ delta_emis(nys:nyn,nxl:nxr)
[4055]	1339
[4559]	1340	ENDIF
[4055]	1341
	1342	ENDDO
[4559]	1343
[4055]	1344	!
[4559]	1345	!-- Any other species
[4055]	1346	ELSE
	1347
[4559]	1348	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &
[4227]	1349	conversion_factor * hours_per_day
[4055]	1350
[4559]	1351	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
	1352	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
	1353	+ delta_emis(nys:nyn,nxl:nxr)
	1354
[4055]	1355	ENDIF ! TRIM spc_names
[4559]	1356
	1357	emis = 0
	1358
[4055]	1359	ENDDO
[4559]	1360
	1361	delta_emis = 0
	1362
[4055]	1363	ENDDO
	1364
	1365	!
[4559]	1366	!-- LOD 2 (PRE-PROCESSED mode)
[4055]	1367	ELSEIF ( emiss_lod == 2 ) THEN
	1368
	1369	!
[4559]	1370	!-- Cycle over species: n_matched_vars represents the number of species in common between the
	1371	!-- emission input data and the chemistry mechanism used
	1372	DO ispec = 1, n_matched_vars
[4055]	1373
[4559]	1374	! (ecc)
	1375	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emt(match_spec_input(ispec))% &
	1376	preproc_emission_data(index_hh+1,1,nys+1:nyn+1,nxl+1:nxr+1) &
	1377	* conversion_factor
[4055]	1378
	1379	ENDDO
	1380
	1381	ENDIF ! emiss_lod
	1382
[4559]	1383
[4055]	1384	!
[4559]	1385	!-- Cycle to transform x,y coordinates to the one of surface_mod and to assign emission values to
	1386	!-- cssws
[4055]	1387
	1388	!
[4559]	1389	!-- LOD 0 (PARAMETERIZED mode)
	1390	!-- Units of inputs are micromoles/m2/s
[4055]	1391	IF ( emiss_lod == 0 ) THEN
	1392
	1393	IF (street_type_f%from_file) THEN
	1394
	1395	!
[4559]	1396	!-- Streets are lsm surfaces, hence, no usm surface treatment required.
	1397	!-- However, urban surface may be initialized via default initialization in surface_mod, e.g. at
	1398	!-- horizontal urban walls that are at k == 0 (building is lower than the first grid point). Hence,
	1399	!-- in order to have only emissions at streets, set the surfaces emissions to zero at urban walls.
[4671]	1400	IF ( surf_usm_h(0)%ns >=1 ) surf_usm_h(0)%cssws = 0.0_wp
[4055]	1401
	1402	!
[4559]	1403	!-- Treat land-surfaces.
[4671]	1404	DO m = 1, surf_lsm_h(0)%ns
[4055]	1405
[4671]	1406	i = surf_lsm_h(0)%i(m)
	1407	j = surf_lsm_h(0)%j(m)
	1408	k = surf_lsm_h(0)%k(m)
[4055]	1409
	1410	!
[4559]	1411	!-- Set everything to zero then reassign according to street type
[4671]	1412	surf_lsm_h(0)%cssws(:,m) = 0.0_wp
[4055]	1413
[4559]	1414	IF ( street_type_f%var(j,i) >= main_street_id .AND. &
[4055]	1415	street_type_f%var(j,i) < max_street_id ) THEN
	1416
	1417	!
[4559]	1418	!-- Cycle over matched species
[4055]	1419	DO ispec = 1, n_matched_vars
	1420
	1421	!
[4559]	1422	!-- PMs are already in kilograms
	1423	IF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM1" .OR. &
	1424	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM25" .OR. &
	1425	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM10" ) THEN
[4055]	1426
[4559]	1427	!
	1428	!-- kg/(m^2s) kg/m^3
[4671]	1429	surf_lsm_h(0)%cssws(match_spec_model(ispec),m) = &
[4559]	1430	emiss_factor_main(match_spec_input(ispec)) * &
	1431	emis_distribution(1,j,i,ispec) * & ! kg/(m^2*s)
	1432	rho_air(k) ! kg/m^3
[4055]	1433
	1434	!
[4559]	1435	!-- Other Species
	1436	!-- Inputs are micromoles
[4055]	1437	ELSE
	1438
[4559]	1439	!
	1440	!-- ppm/s m kg/m^3
[4671]	1441	surf_lsm_h(0)%cssws(match_spec_model(ispec),m) = &
[4559]	1442	emiss_factor_main( match_spec_input(ispec) ) * &
	1443	emis_distribution(1,j,i,ispec) * & ! micromoles/(m^2*s)
	1444	conv_to_ratio(k,j,i) * & ! m^3/Nmole
	1445	rho_air(k) ! kg/m^3
[4055]	1446
	1447	ENDIF
	1448
	1449	ENDDO ! ispec
	1450
	1451
[4559]	1452	ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. &
[4055]	1453	street_type_f%var(j,i) < main_street_id ) THEN
	1454
	1455	!
[4559]	1456	!-- Cycle over matched species
[4055]	1457	DO ispec = 1, n_matched_vars
	1458
	1459	!
[4559]	1460	!-- PMs are already in kilograms
	1461	IF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM1" .OR. &
	1462	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM25" .OR. &
	1463	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM10" ) THEN
[4055]	1464
	1465	!
[4559]	1466	!-- kg/(m^2s) kg/m^3
[4671]	1467	surf_lsm_h(0)%cssws(match_spec_model(ispec),m) = &
[4559]	1468	emiss_factor_side( match_spec_input(ispec) ) * &
	1469	emis_distribution(1,j,i,ispec) * & ! kg/(m^2*s)
	1470	rho_air(k) ! kg/m^3
[4055]	1471	!
[4559]	1472	!-- Other species
	1473	!-- Inputs are micromoles
[4055]	1474	ELSE
	1475
[4559]	1476	!
	1477	!-- ppm/s m kg/m^3
[4671]	1478	surf_lsm_h(0)%cssws(match_spec_model(ispec),m) = &
[4559]	1479	emiss_factor_side( match_spec_input(ispec) ) * &
	1480	emis_distribution(1,j,i,ispec) * & ! micromoles/(m^2*s)
	1481	conv_to_ratio(k,j,i) * & ! m^3/Nmole
	1482	rho_air(k) ! kg/m^3
[4055]	1483
	1484	ENDIF
	1485
	1486	ENDDO ! ispec
	1487
	1488	ENDIF ! street type
	1489
	1490	ENDDO ! m
	1491
	1492	ENDIF ! street_type_f%from_file
	1493
	1494
	1495	!
[4559]	1496	!-- LOD 1 (DEFAULT) and LOD 2 (PRE-PROCESSED)
	1497	ELSE
[4055]	1498
	1499
[4559]	1500	DO ispec = 1, n_matched_vars
[4055]	1501
	1502	!
	1503	!-- Default surfaces
	1504	DO m = 1, surf_def_h(0)%ns
	1505
	1506	i = surf_def_h(0)%i(m)
	1507	j = surf_def_h(0)%j(m)
	1508
	1509	IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN
	1510
	1511	!
[4559]	1512	!-- PMs
	1513	IF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM1" .OR. &
	1514	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM25" .OR. &
	1515	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM10" ) THEN
	1516
	1517	surf_def_h(0)%cssws(match_spec_model(ispec),m) = & ! kg/m2/s * kg/m3
	1518	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
	1519	rho_air(nzb) ! kg/m^3
	1520
[4055]	1521	ELSE
	1522
	1523	!
[4559]	1524	!-- VOCs
	1525	IF ( len_index_voc > 0 .AND. &
[4055]	1526	emt_att%species_name(match_spec_input(ispec)) == "VOC" ) THEN
	1527
[4559]	1528	surf_def_h(0)%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
	1529	emis_distribution(1,j,i,ispec) * & ! mole/m2/s
	1530	conv_to_ratio(nzb,j,i) * & ! m^3/mole
	1531	ratio2ppm * & ! ppm
	1532	rho_air(nzb) ! kg/m^3
[4055]	1533
	1534
	1535	!
[4559]	1536	!-- Other species
[4055]	1537	ELSE
	1538
[4559]	1539	surf_def_h(0)%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
	1540	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
	1541	( 1.0_wp / emt_att%xm(ispec) ) * & ! mole/kg
	1542	conv_to_ratio(nzb,j,i) * & ! m^3/mole
	1543	ratio2ppm * & ! ppm
	1544	rho_air(nzb) ! kg/m^3
	1545
[4055]	1546	ENDIF ! VOC
	1547
	1548	ENDIF ! PM
	1549
	1550	ENDIF ! emis_distribution > 0
	1551
	1552	ENDDO ! m
[4559]	1553
[4055]	1554	!
[4559]	1555	!-- LSM surfaces
[4671]	1556	DO m = 1, surf_lsm_h(0)%ns
[4055]	1557
[4671]	1558	i = surf_lsm_h(0)%i(m)
	1559	j = surf_lsm_h(0)%j(m)
	1560	k = surf_lsm_h(0)%k(m)
[4055]	1561
	1562	IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN
	1563
	1564	!
[4559]	1565	!-- PMs
	1566	IF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM1" .OR. &
	1567	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM25" .OR. &
	1568	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM10" ) THEN
[4055]	1569
[4671]	1570	surf_lsm_h(0)%cssws(match_spec_model(ispec),m) = & ! kg/m2/s * kg/m3
	1571	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
	1572	rho_air(k) ! kg/m^3
[4559]	1573
[4055]	1574	ELSE
	1575
	1576	!
[4559]	1577	!-- VOCs
[4055]	1578	IF ( len_index_voc > 0 .AND. &
	1579	emt_att%species_name(match_spec_input(ispec)) == "VOC" ) THEN
	1580
[4671]	1581	surf_lsm_h(0)%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
	1582	emis_distribution(1,j,i,ispec) * & ! mole/m2/s
	1583	conv_to_ratio(k,j,i) * & ! m^3/mole
	1584	ratio2ppm * & ! ppm
	1585	rho_air(k) ! kg/m^3
[4055]	1586
	1587	!
[4559]	1588	!-- Other species
[4055]	1589	ELSE
	1590
[4671]	1591	surf_lsm_h(0)%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
	1592	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
	1593	( 1.0_wp / emt_att%xm(ispec) ) * & ! mole/kg
	1594	conv_to_ratio(k,j,i) * & ! m^3/mole
	1595	ratio2ppm * & ! ppm
	1596	rho_air(k) ! kg/m^3
[4559]	1597
[4055]	1598	ENDIF ! VOC
	1599
	1600	ENDIF ! PM
	1601
	1602	ENDIF ! emis_distribution
	1603
	1604	ENDDO ! m
	1605
[4403]	1606
	1607
[4055]	1608	!
[4559]	1609	!-- USM surfaces
[4671]	1610	DO m = 1, surf_usm_h(0)%ns
[4055]	1611
[4671]	1612	i = surf_usm_h(0)%i(m)
	1613	j = surf_usm_h(0)%j(m)
	1614	k = surf_usm_h(0)%k(m)
[4055]	1615
	1616	IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN
	1617
	1618	!
[4559]	1619	!-- PMs
	1620	IF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM1" .OR. &
	1621	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM25" .OR. &
	1622	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM10" ) THEN
	1623
[4671]	1624	surf_usm_h(0)%cssws(match_spec_model(ispec),m) = & ! kg/m2/s * kg/m3
	1625	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
	1626	rho_air(k) ! kg/m^3
[4055]	1627
	1628	ELSE
[4559]	1629
[4055]	1630	!
	1631	!-- VOCs
[4559]	1632	IF ( len_index_voc > 0 .AND. &
[4055]	1633	emt_att%species_name(match_spec_input(ispec)) == "VOC" ) THEN
	1634
[4671]	1635	surf_usm_h(0)%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
	1636	emis_distribution(1,j,i,ispec) * & ! m2/s
	1637	conv_to_ratio(k,j,i) * & ! m^3/mole
	1638	ratio2ppm * & ! ppm
	1639	rho_air(k) ! kg/m^3
[4055]	1640
	1641	!
	1642	!-- Other species
	1643	ELSE
	1644
[4671]	1645	surf_usm_h(0)%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
	1646	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
	1647	( 1.0_wp / emt_att%xm(ispec) ) * & ! mole/kg
	1648	conv_to_ratio(k,j,i) * & ! m^3/mole
	1649	ratio2ppm * & ! ppm
	1650	rho_air(k) ! kg/m^3
[4055]	1651
	1652
	1653	ENDIF ! VOC
	1654
	1655	ENDIF ! PM
	1656
	1657	ENDIF ! emis_distribution
[4403]	1658
[4055]	1659	ENDDO ! m
	1660
	1661	ENDDO
	1662
	1663	ENDIF
	1664
	1665	!
[4559]	1666	!-- Deallocate time_factor in case of DEFAULT mode)
	1667	IF ( ALLOCATED( time_factor ) ) DEALLOCATE( time_factor )
[4055]	1668
	1669	ENDIF
	1670
	1671	END SUBROUTINE chem_emissions_setup
	1672
[4403]	1673
	1674	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1675	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1676	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1677	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1678	!!
[4887]	1679	!! 20200203 NB - ON DEMAND EMISSION UPDATE MODE
[4403]	1680	!!
	1681	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1682	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1683	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1684	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1685
	1686
	1687	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1688	!!
	1689	!! WRAPPER / INTERFACE FUNCTIONS
	1690	!!
	1691	!! NOTE - I find using an explicity wrapper provides much better flow control
	1692	!! over an interface block
	1693	!!
	1694	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1695
	1696	!
[4887]	1697	!-- 20200203 NB
[4403]	1698	!
[4559]	1699	!--------------------------------------------------------------------------------------------------!
[4403]	1700	! Description:
	1701	! ------------
	1702	!> interface for initiation of emission arrays based on emission LOD
	1703	!
[4559]	1704	!--------------------------------------------------------------------------------------------------!
[4403]	1705
	1706	SUBROUTINE chem_emissions_header_init
	1707
	1708	IMPLICIT NONE
	1709
	1710	SELECT CASE ( emiss_lod )
	1711	CASE ( 0 )
	1712	! do nothing at the moment
	1713	CASE ( 1 )
	1714	! do nothing at the moment
	1715	CASE ( 2 )
	1716	CALL chem_emissions_header_init_lod2
	1717	END SELECT
	1718
	1719	END SUBROUTINE chem_emissions_header_init
	1720
	1721
	1722	!
[4887]	1723	!-- 20200203 NB
[4403]	1724	!
[4559]	1725	!--------------------------------------------------------------------------------------------------!
[4403]	1726	! Description:
	1727	! ------------
	1728	!> interface for initiation of emission arrays based on emission LOD
[4559]	1729	!--------------------------------------------------------------------------------------------------!
[4403]	1730
	1731	SUBROUTINE chem_emissions_update_on_demand
	1732
	1733	IMPLICIT NONE
	1734
	1735	SELECT CASE ( emiss_lod )
	1736	CASE ( 0 )
	1737	! do nothing at the moment
	1738	CASE ( 1 )
	1739	! do nothing at the moment
	1740	CASE ( 2 )
	1741	CALL chem_emissions_update_on_demand_lod2
	1742	END SELECT
	1743
	1744	END SUBROUTINE ! chem_emisisons_update_on_demand
	1745
[4559]	1746
[4403]	1747	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1748	!!
	1749	!! SUBROUTINES SPECIFIC FOR LOD 2
	1750	!!
	1751	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1752
	1753	!
[4887]	1754	!-- 20200203 NB
[4403]	1755	!
[4559]	1756	!--------------------------------------------------------------------------------------------------!
[4403]	1757	! Description:
	1758	! ------------
	1759	!> Initiates header for emissions data attributes for LOD 2
[4559]	1760	!--------------------------------------------------------------------------------------------------!
[4403]	1761
	1762	SUBROUTINE chem_emissions_header_init_lod2
	1763
[4559]	1764	USE control_parameters, &
[4403]	1765	ONLY: coupling_char, message_string
	1766
[4559]	1767	USE netcdf_data_input_mod, &
	1768	ONLY: chem_emis_att, close_input_file, get_attribute, get_dimension_length, get_variable, &
	1769	open_read_file
[4403]	1770
[4559]	1771
[4403]	1772	IMPLICIT NONE
	1773
[4559]	1774
	1775	INTEGER(iwp) :: att_lod !< lod attribute in chemistry file
[4403]	1776	INTEGER(iwp) :: ncid !< chemistry file netCDF handle
	1777
[4559]	1778	IF ( debug_output ) CALL debug_message( 'chem_emissions_header_init_lod2', 'start' )
[4403]	1779
	1780	!
[4559]	1781	!-- Opens _chemistry input file and obtain header information
	1782	CALL open_read_file ( TRIM( input_file_chem ) // TRIM( coupling_char ), ncid )
[4403]	1783	!
[4559]	1784	!-- Check if LOD in chemistry file matches LOD in namelist
	1785	CALL get_attribute ( ncid, 'lod', att_lod, .TRUE. )
[4403]	1786
	1787	IF ( att_lod /= emiss_lod ) THEN
	1788	message_string = '' ! to get around unused variable warning / error
[4559]	1789	WRITE ( message_string, * ) 'LOD mismatch between namelist (emiss_lod) and', &
	1790	CHAR( 10 ), ' ', 'chemistry input file (global attributes>lod)'
[4403]	1791	CALL message( 'chem_emissions_header_init_lod2', 'CM0468', 1, 2, 0, 6, 0 )
	1792	ENDIF
	1793	!
[4559]	1794	!-- Obtain unit conversion factor
[4403]	1795	CALL get_attribute ( ncid, 'units', chem_emis_att%units, .FALSE., "emission_values" )
	1796	conversion_factor = chem_emissions_convert_base_units ( chem_emis_att%units )
	1797	!
[4559]	1798	!-- Obtain header attributes
[4403]	1799	CALL chem_emissions_init_species ( ncid )
	1800	CALL chem_emissions_init_timestamps ( ncid )
	1801	!
[4559]	1802	!-- Done reading file
[4403]	1803	CALL close_input_file (ncid)
	1804
	1805	!
[4559]	1806	!-- Set previous timestamp index to something different to trigger a read event later on
[4403]	1807	previous_timestamp_index = -1
	1808
[4559]	1809	IF ( debug_output ) CALL debug_message( 'chem_emissions_header_init_lod2', 'end' )
	1810
[4403]	1811	END SUBROUTINE chem_emissions_header_init_lod2
	1812
	1813	!
[4887]	1814	!-- 20200203 NB
[4403]	1815	!
[4559]	1816	!--------------------------------------------------------------------------------------------------!
[4403]	1817	! Description:
	1818	! ------------
	1819	!> Reads emission data on demand for LOD2
[4559]	1820	!--------------------------------------------------------------------------------------------------!
[4403]	1821
	1822	SUBROUTINE chem_emissions_update_on_demand_lod2
	1823
[4559]	1824	USE control_parameters, &
	1825	ONLY: coupling_char, time_since_reference_point
[4403]	1826
[4559]	1827	USE netcdf_data_input_mod, &
	1828	ONLY: chem_emis_att, close_input_file, get_variable, open_read_file
[4403]	1829
[4559]	1830	USE arrays_3d, &
	1831	ONLY: hyp, pt
[4403]	1832
[4559]	1833	USE surface_mod, &
[4403]	1834	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
	1835
	1836
	1837	IMPLICIT NONE
	1838
[4559]	1839	CHARACTER(LEN=80) :: this_timestamp !< writes out timestamp
[4403]	1840
	1841	INTEGER(iwp) :: i,j,k,m !< generic counters
	1842	INTEGER(iwp) :: kmatch !< index of matched species
	1843	INTEGER(iwp) :: ncid !< netCDF file handle (chemistry file)
	1844	INTEGER(iwp) :: time_index_location !< location of most recent timestamp
	1845
	1846	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: cssws_def_h !< dummy default surface array
	1847	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: cssws_lsm_h !< dummy LSM surface array
	1848	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: cssws_usm_h !< dummy USM surface array
[4559]	1849	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: mass2mole !< conversion factor mass 2 molar (ppm) flux
[4403]	1850
[4559]	1851	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: emis_distrib !< surface emissions
	1852
	1853	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:) :: emissions_raw !< raw emissions data
	1854
	1855	IF ( debug_output ) CALL debug_message ( 'chem_emissions_update_on_demand_lod2', 'start' )
[4403]	1856	!
[4559]	1857	!-- Obtain current timestamp and locate index for most recent timestamp element
	1858	!-- end subroutine (RETURN) if it is still the same index as the existing time index
[4403]	1859
	1860	this_timestamp = '' ! string must be initiated before using
	1861	CALL get_date_time( time_since_reference_point, date_time_str=this_timestamp )
	1862
[4559]	1863	time_index_location = chem_emissions_locate_timestep &
	1864	( this_timestamp, timestamps, 1, chem_emis_att%dt_emission )
[4403]	1865
	1866	IF ( time_index_location == previous_timestamp_index ) RETURN
	1867
	1868	!
[4559]	1869	!-- Begin extract emission data for matched species from netCDF file
[4403]	1870	previous_timestamp_index = time_index_location
	1871
	1872	ALLOCATE ( emis_distrib(n_matched_vars,nys:nyn,nxl:nxr) )
	1873	emis_distrib = 0.0_wp
	1874
	1875	!
[4559]	1876	!-- Open netCDF file and allocate temp memory
	1877	CALL open_read_file( TRIM( input_file_chem ) // TRIM( coupling_char ), ncid )
[4403]	1878	ALLOCATE( emissions_raw(1,1,nys:nyn,nxl:nxr,1) )
	1879
	1880	DO k = 1, n_matched_vars
	1881	!
[4559]	1882	!-- Get index for matching species
	1883	kmatch = chem_emissions_locate_species( spc_names(match_spec_model(k)), &
	1884	chem_emis_att%species_name )
[4403]	1885	!
[4559]	1886	!-- Extract variable as-is
	1887	!-- Note C index notations for nx and ny due to MPI and reversed index dimension order for netCDF
	1888	!-- Fortran API)
[4403]	1889	emissions_raw = 0.0_wp
	1890
[4559]	1891	CALL get_variable ( ncid, 'emission_values', emissions_raw, &
	1892	kmatch, nxl+1, nys+1, 1, time_index_location, &
[4403]	1893	1, nxr-nxl+1, nyn-nys+1, 1, 1, .FALSE. )
	1894	!
[4559]	1895	!-- Transfer emission data
[4403]	1896	DO j = nys,nyn
	1897	DO i = nxl,nxr
	1898	emis_distrib(k,j,i) = emissions_raw(1,1,j,i,1) * conversion_factor
	1899	ENDDO
	1900	ENDDO
	1901
	1902	ENDDO ! k = n_matched_vars
	1903	!
	1904	!-- netCDF handle and temp memory no longer needed
	1905	DEALLOCATE( emissions_raw )
	1906	CALL close_input_file( ncid )
	1907	!
[4559]	1908	!-- Set emis_dt since chemistry ODEs can be stiff, the option to solve them at every RK substep is
	1909	!-- present to help improve stability should the need arise
[4403]	1910	dt_emis = dt_3d
	1911
[4559]	1912	IF ( call_chem_at_all_substeps ) dt_emis = dt_emis * weight_pres(intermediate_timestep_count)
[4403]	1913	!
[4559]	1914	!-- Calculate conversion factor from mass flux to molar flux (mixing ratio)
[4403]	1915	ALLOCATE ( mass2mole(nys:nyn,nxl:nxr) )
	1916	mass2mole = 0.0_wp
	1917
	1918	DO i = nxl, nxr
	1919	DO j = nys, nyn
	1920	mass2mole(j,i) = mass_2_molar_flux ( hyp(nzb), pt(nzb,j,i) )
	1921	ENDDO
	1922	ENDDO
	1923
	1924	!
[4559]	1925	!-- Calculate surface fluxes
	1926	!-- NOTE - For some reason I can not pass surf_xxx%cssws as output argument into subroutine
	1927	!-- assign_surface_flux ( ). The contents got mixed up once the subroutine is finished. I
	1928	!-- don't know why and I don't have time to investigate. As workaround I declared dummy
	1929	!-- variables and reassign them one by one (i.e., in a loop)
[4403]	1930	!-- ECC 20200206
	1931
	1932	!
[4559]	1933	!-- Allocate and initialize dummy surface fluxes
[4403]	1934	ALLOCATE( cssws_def_h(n_matched_vars,surf_def_h(0)%ns) )
	1935	cssws_def_h = 0.0_wp
	1936
[4671]	1937	ALLOCATE( cssws_lsm_h(n_matched_vars,surf_lsm_h(0)%ns) )
[4403]	1938	cssws_lsm_h = 0.0_wp
	1939
[4671]	1940	ALLOCATE( cssws_usm_h(n_matched_vars,surf_usm_h(0)%ns) )
[4403]	1941	cssws_usm_h = 0.0_wp
	1942
	1943	!
[4559]	1944	!-- Assign and transfer emissions as surface fluxes
	1945	CALL assign_surface_flux ( cssws_def_h, surf_def_h(0)%ns, &
	1946	surf_def_h(0)%j, surf_def_h(0)%i, &
[4403]	1947	emis_distrib, mass2mole )
	1948
	1949
[4671]	1950	CALL assign_surface_flux ( cssws_lsm_h, surf_lsm_h(0)%ns, &
	1951	surf_lsm_h(0)%j, surf_lsm_h(0)%i, &
[4403]	1952	emis_distrib, mass2mole )
	1953
	1954
[4671]	1955	CALL assign_surface_flux ( cssws_usm_h, surf_usm_h(0)%ns, &
	1956	surf_usm_h(0)%j, surf_usm_h(0)%i, &
[4403]	1957	emis_distrib, mass2mole )
	1958
	1959	DO k = 1, n_matched_vars
	1960
	1961	DO m = 1, surf_def_h(0)%ns
	1962	surf_def_h(0)%cssws(k,m) = cssws_def_h(k,m)
	1963	ENDDO
	1964
[4671]	1965	DO m = 1, surf_lsm_h(0)%ns
	1966	surf_lsm_h(0)%cssws(k,m) = cssws_lsm_h(k,m)
[4403]	1967	ENDDO
	1968
[4671]	1969	DO m = 1, surf_usm_h(0)%ns
	1970	surf_usm_h(0)%cssws(k,m) = cssws_usm_h(k,m)
[4403]	1971	ENDDO
	1972
	1973	ENDDO
	1974
	1975	!
[4559]	1976	!-- Cleaning up
[4403]	1977	DEALLOCATE( cssws_def_h )
	1978	DEALLOCATE( cssws_lsm_h )
	1979	DEALLOCATE( cssws_usm_h )
	1980
	1981	DEALLOCATE ( emis_distrib )
	1982	DEALLOCATE ( mass2mole )
	1983
[4559]	1984	IF ( debug_output ) CALL debug_message ( 'chem_emissions_update_on_demand_lod2', 'end' )
[4403]	1985
	1986	END SUBROUTINE ! chem_emissions_update_on_demand_lod2
	1987
[4559]	1988
[4403]	1989	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1990	!!
	1991	!! AUXILIARY SUBROUTINES
	1992	!!
	1993	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1994
	1995	!
[4887]	1996	!-- 20200203 NB
[4403]	1997	!
[4559]	1998	!--------------------------------------------------------------------------------------------------!
[4403]	1999	! Description:
	2000	! ------------
[4559]	2001	!> Look for matched species between emissions attributes and selected chemical mechanisms and
	2002	!> determine corresponding molecular weights
	2003	!--------------------------------------------------------------------------------------------------!
[4403]	2004
	2005	SUBROUTINE chem_emissions_init_species ( ncid )
	2006
[4559]	2007	USE netcdf_data_input_mod, &
	2008	ONLY: chem_emis_att, close_input_file, get_dimension_length, get_variable, open_read_file
[4403]	2009
	2010	IMPLICIT NONE
	2011
	2012	INTEGER(iwp) :: ispec !< generic counter 4 species
	2013
	2014	INTEGER(iwp), INTENT(IN) :: ncid !< netcdf file ID
	2015
[4559]	2016	IF ( debug_output ) CALL debug_message( 'chem_emissions_init_species', 'start' )
[4403]	2017	!
[4559]	2018	!- Assign species
[4403]	2019	CALL get_dimension_length ( ncid, chem_emis_att%n_emiss_species, 'nspecies' )
[4559]	2020	CALL get_variable ( ncid, 'emission_name', chem_emis_att%species_name, &
	2021	chem_emis_att%n_emiss_species )
[4403]	2022	!
[4887]	2023	!- Backward compatibility for salsa_mod NB
[4403]	2024	chem_emis_att%nspec = chem_emis_att%n_emiss_species
	2025	!
[4559]	2026	!-- Get a list of matched species between emission_attributes and selected chemical mechanism
[4403]	2027	emission_output_required = .FALSE.
	2028	CALL chem_emissions_match( chem_emis_att, n_matched_vars )
	2029
	2030	!
[4559]	2031	!-- If matched species found (at least 1),
	2032	!-- allocate memory for emission attributes,
	2033	!-- assign molecular masses [kg/mol],
	2034	!-- see chemistry_model_mod.f90 for reference.
[4403]	2035	IF ( n_matched_vars > 0 ) THEN
	2036
	2037	emission_output_required = .TRUE.
	2038
	2039	ALLOCATE( chem_emis_att%xm(n_matched_vars) )
	2040
	2041	DO ispec = 1, n_matched_vars
	2042	chem_emis_att%xm(ispec) = 1.0_wp
	2043	SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
	2044	CASE ( 'SO2' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 2
	2045	CASE ( 'SO4' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 4
	2046	CASE ( 'NO' ); chem_emis_att%xm(ispec) = xm_N + xm_O
	2047	CASE ( 'NO2' ); chem_emis_att%xm(ispec) = xm_N + xm_O * 2
	2048	CASE ( 'NH3' ); chem_emis_att%xm(ispec) = xm_N + xm_H * 3
	2049	CASE ( 'CO' ); chem_emis_att%xm(ispec) = xm_C + xm_O
	2050	CASE ( 'CO2' ); chem_emis_att%xm(ispec) = xm_C + xm_O * 2
	2051	CASE ( 'CH4' ); chem_emis_att%xm(ispec) = xm_C + xm_H * 4
	2052	CASE ( 'HNO3' ); chem_emis_att%xm(ispec) = xm_H + xm_N + xm_O*3
	2053	END SELECT
	2054	ENDDO
[4559]	2055
[4403]	2056	ENDIF ! IF ( n_matched_vars > 0 )
	2057
[4559]	2058	IF ( debug_output ) CALL debug_message( 'chem_emissions_init_species', 'end' )
[4403]	2059
	2060	END SUBROUTINE chem_emissions_init_species
	2061
[4559]	2062
[4403]	2063	!
[4887]	2064	!-- 20200203 NB
[4403]	2065	!
[4559]	2066	!--------------------------------------------------------------------------------------------------!
[4403]	2067	! Description:
	2068	! ------------
[4559]	2069	!> Extract timestamps from netCDF input
	2070	!--------------------------------------------------------------------------------------------------!
[4403]	2071
	2072	SUBROUTINE chem_emissions_init_timestamps ( ncid )
	2073
[4559]	2074	USE control_parameters, &
[4403]	2075	ONLY: message_string
	2076
[4559]	2077	USE netcdf_data_input_mod, &
	2078	ONLY: chem_emis_att, close_input_file, get_dimension_length, get_variable, open_read_file
[4403]	2079
	2080	IMPLICIT NONE
	2081
	2082	INTEGER(iwp) :: fld_len !< string field length
	2083	INTEGER(iwp) :: itime !< generic counter (4 species)
	2084
	2085	INTEGER(iwp), INTENT(IN) :: ncid !< netcdf file handle
	2086
[4559]	2087	IF ( debug_output ) CALL debug_message( 'chem_emissions_read_timestamps', 'start' )
[4403]	2088	!
[4559]	2089	!-- Import timestamps from netCDF input
[4403]	2090	CALL get_dimension_length ( ncid, chem_emis_att%dt_emission, 'time' )
	2091	CALL get_dimension_length ( ncid, fld_len, 'field_length' )
	2092	CALL get_variable ( ncid, 'timestamp', timestamps, chem_emis_att%dt_emission, fld_len )
	2093	!
[4559]	2094	!-- Throw error at first instance of timestamps not in listed in chronological order.
[4403]	2095	DO itime = 2,chem_emis_att%dt_emission
	2096
	2097	IF ( timestamps(itime-1) > timestamps(itime) ) THEN
	2098
[4559]	2099	WRITE( message_string, * ) &
	2100	'input timestamps not in chronological order for', &
	2101	CHAR( 10 ), ' ', &
	2102	'index ', (itime-1), ' : ', TRIM( timestamps(itime-1) ), ' and', &
	2103	CHAR( 10 ), ' ', &
	2104	'index ', (itime), ' : ', TRIM( timestamps(itime) )
[4403]	2105
	2106	CALL message( 'chem_emissions_read_timestamps', 'CM0469', 1, 2, 0, 6, 0 )
	2107
	2108	ENDIF
	2109
	2110	ENDDO
	2111
[4559]	2112	IF ( debug_output ) CALL debug_message( 'chem_emissions_read_timestamps', 'end' )
[4403]	2113
	2114	END SUBROUTINE chem_emissions_init_timestamps
	2115
	2116
	2117	!
[4887]	2118	!-- 20200203 NB
[4403]	2119	!
[4559]	2120	!--------------------------------------------------------------------------------------------------!
[4403]	2121	! Description:
	2122	! ------------
[4559]	2123	!> Assign emissions as surface fluxes
[4403]	2124	!
[4559]	2125	!> NOTE: For arguments, I originally wanted to use unspecified dimensions, but I could not get
	2126	!> this to work properly, hence the dimensioned array arguments.
	2127	!--------------------------------------------------------------------------------------------------!
[4403]	2128
[4559]	2129	SUBROUTINE assign_surface_flux ( surf_array, nsurfs, surf_j, surf_i, emis_dist, conv_mole )
[4403]	2130
[4559]	2131	USE arrays_3d, &
[4403]	2132	ONLY: rho_air
	2133
[4559]	2134	USE netcdf_data_input_mod, &
[4403]	2135	ONLY: chem_emis_att
	2136
	2137	USE surface_mod !< for surf_type
	2138
	2139	IMPLICIT NONE
	2140	!
[4559]	2141	!-- Input arguments
	2142	INTEGER(iwp), INTENT(IN) :: nsurfs !< # surfaces in surf_array
[4403]	2143
	2144	INTEGER(iwp), DIMENSION(nsurfs), INTENT(IN) :: surf_i !< i indices 4 surf. elements
	2145	INTEGER(iwp), DIMENSION(nsurfs), INTENT(IN) :: surf_j !< j indices 4 surf. elements
	2146
	2147	REAL(wp), DIMENSION(nys:nyn,nxl:nxr), INTENT(IN) :: conv_mole !< conv. 2 molar flux
	2148	REAL(wp), DIMENSION(n_matched_vars,nys:nyn,nxl:nxr), INTENT(IN) :: emis_dist !< surf. emissions
	2149
	2150	REAL(wp), DIMENSION(n_matched_vars,nsurfs), INTENT(INOUT) :: surf_array !< surface listing
	2151
	2152	!
[4559]	2153	!-- Parameters (magic numbers)
[4403]	2154	CHARACTER(LEN=2), PARAMETER :: sp_PM = 'PM' !< id string for all PMs
	2155	CHARACTER(LEN=3), PARAMETER :: sp_VOC = 'VOC' !< id string for VOC
	2156
	2157	REAL(wp), PARAMETER :: mol2ppm = 1.0E+06_wp !< conversion from mole 2 ppm
	2158	!
[4559]	2159	!-- Local variables
[4403]	2160	CHARACTER(LEN=80) :: this_species_name !< matched species name
	2161
	2162	INTEGER(iwp) :: i,j,k,m !< generic counters
	2163
	2164	REAL(wp) :: flux_conv_factor !< conversion factor
	2165
[4559]	2166	IF ( debug_output ) CALL debug_message( 'chem_emissions_header_init_lod2', 'start' )
[4403]	2167
	2168	DO k = 1, n_matched_vars
	2169
	2170	this_species_name = spc_names(k) !< species already matched
	2171
	2172	DO m = 1, nsurfs
	2173
	2174	j = surf_j(m) ! get surface coordinates
	2175	i = surf_i(m)
	2176
	2177	!
[4559]	2178	!-- Calculate conversion factor depending on emission species type
[4403]	2179	flux_conv_factor = rho_air(nzb)
	2180	!
[4559]	2181	!-- Account for conversion to different types of emisison species
	2182	IF ( TRIM( this_species_name( 1:LEN( sp_PM ) ) ) == sp_PM ) THEN
[4403]	2183
	2184	! do nothing (use mass flux directly)
	2185
[4559]	2186	ELSE IF ( TRIM( this_species_name( 1:LEN( sp_VOC ) ) ) == sp_VOC ) THEN
[4403]	2187
[4559]	2188	flux_conv_factor = flux_conv_factor * conv_mole(j,i) * mol2ppm
[4403]	2189
	2190	ELSE
	2191
[4559]	2192	flux_conv_factor = flux_conv_factor * conv_mole(j,i) * mol2ppm / chem_emis_att%xm(k)
[4403]	2193
	2194	ENDIF
	2195	!
[4559]	2196	!-- Finally assign surface flux
[4403]	2197	surf_array(k,m) = emis_dist(k,j,i) * flux_conv_factor
	2198
	2199	ENDDO ! m = 1, nsurfs
	2200
	2201	ENDDO ! k = 1, n_matched_vars
	2202
	2203
[4559]	2204	IF ( debug_output ) CALL debug_message( 'chem_emissions_header_init_lod2', 'end' )
	2205
[4403]	2206	END SUBROUTINE assign_surface_flux
	2207
[4559]	2208
[4403]	2209	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2210	!!
	2211	!! AUXILIARY FUNCTIONS
	2212	!!
	2213	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2214
	2215	!
[4887]	2216	!-- 20200203 NB
[4403]	2217	!
[4559]	2218	!--------------------------------------------------------------------------------------------------!
[4403]	2219	! Description:
	2220	! ------------
[4559]	2221	!> Given incoming flux units ( mass / area / time ) provide single-valued onversion factor to
	2222	!> ( kg / m2 / s )
	2223	!--------------------------------------------------------------------------------------------------!
[4403]	2224
	2225	FUNCTION chem_emissions_convert_base_units ( units_in ) RESULT ( conv_factor )
	2226
	2227	IMPLICIT NONE
	2228	!
[4559]	2229	!-- Function arguments
	2230	CHARACTER(LEN=*), INTENT(IN) :: units_in !< incoming units (ie emt_att%units)
[4403]	2231
	2232	REAL(wp) :: conv_factor !< convertion factor
	2233
	2234	!
[4559]	2235	!-- Parameters (magic numbers)
[4403]	2236	INTEGER(iwp), PARAMETER :: up2lo = 32 !< convert letter to lower case
	2237	!
[4559]	2238	!-- Base unit conversion factors (should be self-explanatory)
	2239	REAL(wp), PARAMETER :: hour_per_year = 8760.0_wp
[4403]	2240	REAL(wp), PARAMETER :: g_to_kg = 1.0E-03_wp
	2241	REAL(wp), PARAMETER :: miug_to_kg = 1.0E-09_wp
	2242	REAL(wp), PARAMETER :: s_per_hour = 3600.0_wp
	2243	REAL(wp), PARAMETER :: s_per_day = 86400.0_wp
[4559]	2244	REAL(wp), PARAMETER :: tons_to_kg = 100.0_wp
[4403]	2245
	2246	REAL(wp), PARAMETER :: day_to_s = 1.0_wp / s_per_day
	2247	REAL(wp), PARAMETER :: hour_to_s = 1.0_wp / s_per_hour
[4559]	2248	REAL(wp), PARAMETER :: year_to_s = 1.0_wp / s_per_hour / hour_per_year
	2249
	2250
[4403]	2251	!
[4559]	2252	!-- Local variables
[4403]	2253	CHARACTER(LEN=LEN(units_in)) :: units_in_lo !< units in lower case
	2254
	2255	INTEGER(iwp) :: j,k !< generic counters
	2256	INTEGER(iwp) :: str_len !< length of unit string
	2257	!
[4559]	2258	!-- Turn units string to lower case
[4403]	2259	units_in_lo = ''
[4559]	2260	str_len = LEN( TRIM( units_in ) )
[4403]	2261
	2262	DO k = 1,str_len
	2263	j = IACHAR( units_in(k:k) )
[4559]	2264	units_in_lo(k:k) = ACHAR( j )
	2265	IF ( ( j >= IACHAR( "A" ) ) .AND. ( j <= IACHAR( "Z" ) ) ) &
	2266	units_in_lo(k:k) = ACHAR( j + up2lo )
[4403]	2267	ENDDO
	2268
	2269	conv_factor = 1.0_wp !< default value (kg/m2/s)
	2270
	2271	SELECT CASE ( TRIM( units_in_lo ) )
	2272	CASE ( 'kg/m2/s' ); conv_factor = 1.0_wp
	2273	CASE ( 'kg/m2/hour' ); conv_factor = hour_to_s
	2274	CASE ( 'kg/m2/day' ); conv_factor = day_to_s
	2275	CASE ( 'kg/m2/year' ); conv_factor = year_to_s
	2276	CASE ( 'ton/m2/s' ); conv_factor = tons_to_kg
	2277	CASE ( 'ton/m2/hour' ); conv_factor = tons_to_kg * hour_to_s
	2278	CASE ( 'ton/m2/day' ); conv_factor = tons_to_kg * day_to_s
	2279	CASE ( 'ton/m2/year' ); conv_factor = tons_to_kg * year_to_s
	2280	CASE ( 'g/m2/s' ); conv_factor = g_to_kg
	2281	CASE ( 'g/m2/hour' ); conv_factor = g_to_kg * hour_to_s
	2282	CASE ( 'g/m2/day' ); conv_factor = g_to_kg * day_to_s
	2283	CASE ( 'g/m2/year' ); conv_factor = g_to_kg * year_to_s
	2284	CASE ( 'micrograms/m2/s' ); conv_factor = miug_to_kg
	2285	CASE ( 'micrograms/m2/hour' ); conv_factor = miug_to_kg * hour_to_s
	2286	CASE ( 'micrograms/m2/day' ); conv_factor = miug_to_kg * day_to_s
	2287	CASE ( 'micrograms/m2/year' ); conv_factor = miug_to_kg * year_to_s
	2288	CASE DEFAULT
	2289	message_string = '' ! to get around unused variable warning / error
[4559]	2290	WRITE ( message_string, * ) 'Specified emission units (', TRIM( units_in ), &
[4403]	2291	') not recognized in PALM-4U'
	2292	CALL message ( 'chem_emission_convert_units', 'CM0446', 2, 2, 0, 6, 0 )
	2293	END SELECT
	2294
	2295	END FUNCTION chem_emissions_convert_base_units
	2296
	2297
	2298	!
[4887]	2299	!-- 20200203 NB
[4403]	2300	!
[4559]	2301	!--------------------------------------------------------------------------------------------------!
[4403]	2302	! Description:
	2303	! ------------
[4559]	2304	!> Calculates conversion factor from mass flux to ppm (molar flux)
	2305	!--------------------------------------------------------------------------------------------------!
[4403]	2306
	2307	FUNCTION mass_2_molar_flux ( rhogh, theta ) RESULT ( conv_factor )
	2308
[4559]	2309	USE basic_constants_and_equations_mod, &
	2310	ONLY: p_0, rd_d_cp, rgas_univ
[4403]	2311
	2312	IMPLICIT NONE
	2313	!
[4559]	2314	!-- Function arguments
[4403]	2315	REAL(wp) :: conv_factor !< conversion factor
	2316	REAL(wp), INTENT(IN) :: rhogh !< hydrostatic pressure
	2317	REAL(wp), INTENT(IN) :: theta !< potential temperature
	2318
	2319	conv_factor = ( rgas_univ / rhogh ) * theta * ( ( rhogh / p_0 ) ** rd_d_cp )
	2320
	2321	END FUNCTION mass_2_molar_flux
	2322
	2323
	2324	!
[4887]	2325	!-- 20200203 NB
[4403]	2326	!
[4559]	2327	!--------------------------------------------------------------------------------------------------!
[4403]	2328	! Description:
	2329	! ------------
[4559]	2330	!> Given target sepecies locate index in species array
[4403]	2331	!> returns 0 if none is found
[4559]	2332	!--------------------------------------------------------------------------------------------------!
[4403]	2333
[4559]	2334	FUNCTION chem_emissions_locate_species ( this_species, species_array ) RESULT ( species_index )
[4403]	2335
	2336	IMPLICIT NONE
	2337	!
[4559]	2338	!-- Function arguments
[4403]	2339	INTEGER(iwp) :: species_index !> index matching species
	2340
[4559]	2341	CHARACTER(LEN=25), INTENT(IN) :: species_array(:) !> array of species
[4403]	2342	CHARACTER(LEN=*), INTENT(IN) :: this_species !> target species
	2343	!
[4559]	2344	!-- Local variables
[4403]	2345	INTEGER(iwp) :: k !> generic counter
	2346	INTEGER(iwp) :: n_species !> number of species in species_array
	2347
	2348	n_species = SIZE( species_array, 1 )
	2349
	2350	DO k = 1, n_species
[4559]	2351	IF ( TRIM( species_array(k) ) == TRIM( this_species ) ) EXIT
[4403]	2352	ENDDO
	2353
	2354	species_index = 0 !> assume no matching index is found
	2355
[4559]	2356	IF ( TRIM( species_array(k) ) == TRIM( this_species ) ) specieS_index = k
[4403]	2357
	2358	END FUNCTION chem_emissions_locate_species
	2359
	2360
	2361	!
[4887]	2362	!-- 20200203 NB
[4403]	2363	!
[4559]	2364	!--------------------------------------------------------------------------------------------------!
[4403]	2365	! Description:
	2366	! ------------
	2367	!> given target timestamp locate most recent timestep in timestamp array
	2368	!> using bisection search (since array is sorted)
[4559]	2369	!--------------------------------------------------------------------------------------------------!
[4403]	2370
[4559]	2371	RECURSIVE FUNCTION chem_emissions_locate_timestep &
	2372	( this_timestamp, timestamp_array, lower_bound, upper_bound ) &
[4403]	2373	RESULT ( timestamp_index )
	2374
	2375	!
[4559]	2376	!-- Function arguments
[4403]	2377	CHARACTER(LEN=*), INTENT(IN) :: this_timestamp !> target timestamp
	2378	CHARACTER(LEN=512), INTENT(IN) :: timestamp_array(:) !> array of timestamps
	2379
	2380	INTEGER(iwp), INTENT(IN) :: lower_bound, upper_bound !> timestamp_array index bounds
	2381
[4559]	2382	INTEGER(iwp) :: timestamp_index !> index for most recent timestamp in timestamp_array
	2383
[4403]	2384	!
[4559]	2385	!-- Local variables
[4403]	2386	INTEGER(iwp) :: k0,km,k1 !> lower, central, and upper index bounds
	2387	!
[4559]	2388	!-- Assign bounds
[4403]	2389	k0 = lower_bound
	2390	k1 = upper_bound
	2391	!
[4559]	2392	!-- Make sure k1 is always not smaller than k0
[4403]	2393	IF ( k0 > k1 ) THEN
	2394	k0 = upper_bound
	2395	k1 = lower_bound
	2396	ENDIF
	2397	!
[4559]	2398	!-- Make sure k0 and k1 stay within array bounds by timestamp_array
[4403]	2399	IF ( k0 < 1 ) k0 = 1
[4559]	2400	IF ( k1 > SIZE( timestamp_array, 1 ) ) k1 = SIZE( timestamp_array, 1 )
[4403]	2401	!
[4559]	2402	!-- Terminate if target is contained within 2 consecutive indices otherwise calculate central bound
	2403	!-- (km) and determine new index bounds for the next iteration
[4403]	2404
	2405	IF ( ( k1 - k0 ) > 1 ) THEN
	2406	km = ( k0 + k1 ) / 2
[4559]	2407	IF ( TRIM( this_timestamp ) > TRIM( timestamp_array(km) ) ) THEN
[4403]	2408	k0 = km
	2409	ELSE
	2410	k1 = km
	2411	ENDIF
[4559]	2412	timestamp_index = chem_emissions_locate_timestep( this_timestamp, timestamp_array, k0, k1 )
[4403]	2413	ELSE
	2414	timestamp_index = k0
	2415	ENDIF
	2416
	2417	END FUNCTION chem_emissions_locate_timestep
	2418
	2419
	2420	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2421	!!
	2422	!! END OF MODULE
	2423	!!
	2424	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2425
[4887]	2426	END MODULE chem_emissions_mod

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

Context Navigation

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

Download in other formats: