Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |