Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |