Home

Context Navigation

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

Last change on this file since 4468 was 4403, checked in by banzhafs, 5 years ago
chemistry model: implemented on-demand emission read mode for LOD 2
Property svn:keywords set to `Id`
File size: 96.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |