Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |