Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |