Home

Context Navigation

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

Last change on this file since 4016 was 3968, checked in by suehring, 5 years ago
Updates from chemistriy branched merged into trunk: code cleaning and formatting, code structure optimizations
Property svn:keywords set to `Id`
File size: 65.0 KB

Rev	Line
[3968]	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	!
[3968]	26	! Former revisions:
	27	! -----------------
	28	! $Id: chem_emissions_mod.f90 3968 2019-05-13 11:04:01Z forkel $
	29	! - in subroutine chem_emissions_match replace all decision structures relating to
	30	! mode_emis to emiss_lod
	31	! - in subroutine chem_check_parameters replace emt%nspec with emt%n_emiss_species
	32	! - spring cleaning (E.C. Chan)
	33	!
	34	! 3885 2019-04-11 11:29:34Z kanani
	35	! Changes related to global restructuring of location messages and introduction
	36	! of additional debug messages
	37	!
	38	! 3831 2019-03-28 09:11:22Z forkel
	39	! added nvar to USE chem_gasphase_mod (chem_modules will not include nvar anymore)
	40	!
	41	! 3788 2019-03-07 11:40:09Z banzhafs
	42	! Removed unused variables from chem_emissions_mod
	43	!
	44	! 3772 2019-02-28 15:51:57Z suehring
	45	! - In case of parametrized emissions, assure that emissions are only on natural
	46	! surfaces (i.e. streets) and not on urban surfaces.
	47	! - some unnecessary if clauses removed
	48	!
	49	! 3685 2019 -01-21 01:02:11Z knoop
	50	! Some interface calls moved to module_interface + cleanup
	51	!
	52	! 3611 2018-12-07 14:14:11Z banzhafs
	53	! Code update to comply PALM coding rules
	54	! removed unnecessary informative messsages/location
	55	! messages and corrected comments on PM units from to kg
	56	! bug fix: spcs_name replaced by nvar in DO loops
	57	!
	58	! 3591 2018-11-30 17:37:32Z suehring
	59	! - Removed salsa dependency.
	60	! - Enabled PARAMETRIZED mode for default surfaces when LSM is not applied but
	61	! salsa is (M. Kurppa)
	62	!
	63	! 3582 2018-11-29 19:16:36Z suehring
	64	! resler:
	65	! Break lines at 132 characters
	66	!
	67	! 3483 2018-11-02 14:19:26Z raasch
	68	! bugfix: wrong locations of netCDF directives fixed
	69	!
	70	! 3467 2018-10-30 19:05:21Z suehring
	71	! Enabled PARAMETRIZED mode for default surfaces when LSM is not applied but
	72	! salsa is used
	73	!
	74	! 3458 2018-10-30 14:51:23Z kanani
	75	! from chemistry branch r3443, banzhafs, Russo:
	76	! Additional correction for index of input file of pre-processed mode
	77	! Removed atomic and molecular weights as now available in chem_modules and
	78	! added par_emis_time_factor (formerly in netcdf_data_input_mod)
	79	! Added correction for index of input file of pre-processed mode
	80	! Added a condition for default mode necessary for information read from ncdf file
	81	! in pre-processed and default mode
	82	! Correction of multiplying factor necessary for scaling emission values in time
	83	! Introduced correction for scaling units in the case of DEFAULT emission mode
	84	!
	85	! 3373 2018-10-18 15:25:56Z kanani
	86	! Fix wrong location of __netcdf directive
	87	!
	88	! 3337 2018-10-12 15:17:09Z kanani
	89	! (from branch resler)
	90	! Formatting
	91	!
	92	! 3312 2018-10-06 14:15:46Z knoop
	93	! Initial revision
	94	!
	95	! 3286 2018-09-28 07:41:39Z forkel
	96	!
	97	! Authors:
	98	! --------
	99	! @author Emmanuele Russo (FU-Berlin)
	100	! @author Sabine Banzhaf (FU-Berlin)
	101	! @author Martijn Schaap (FU-Berlin, TNO Utrecht)
	102	!
	103	! Description:
	104	! ------------
	105	!> MODULE for reading-in Chemistry Emissions data
	106	!>
	107	!> @todo Rename nspec to n_emis to avoid inteferece with nspec from chem_gasphase_mod
	108	!> @todo Check_parameters routine should be developed: for now it includes just one condition
	109	!> @todo Use of Restart files not contempled at the moment
	110	!> @todo revise indices of files read from the netcdf: preproc_emission_data and expert_emission_data
	111	!> @todo for now emission data may be passed on a singular vertical level: need to be more flexible
	112	!> @todo fill/activate restart option in chem_emissions_init
	113	!> @todo discuss dt_emis
	114	!> @note <Enter notes on the module>
	115	!> @bug <Enter known bugs here>
	116	!------------------------------------------------------------------------------!
	117
	118	MODULE chem_emissions_mod
	119
	120	USE arrays_3d, &
	121	ONLY: rho_air
	122
	123	USE control_parameters, &
	124	ONLY: debug_output, &
	125	end_time, message_string, initializing_actions, &
	126	intermediate_timestep_count, dt_3d
	127
	128	USE indices
	129
	130	USE kinds
	131
	132	#if defined( __netcdf )
	133	USE netcdf
	134	#endif
	135
	136	USE netcdf_data_input_mod, &
	137	ONLY: chem_emis_att_type, chem_emis_val_type
	138
	139	USE date_and_time_mod, &
	140	ONLY: day_of_month, hour_of_day, &
	141	index_mm, index_dd, index_hh, &
	142	month_of_year, hour_of_day, &
	143	time_default_indices, time_preprocessed_indices
	144
	145	USE chem_gasphase_mod, &
	146	ONLY: nvar, spc_names
	147
	148	USE chem_modules
	149
	150	USE statistics, &
	151	ONLY: weight_pres
	152
	153
	154	IMPLICIT NONE
	155
	156	!
	157	!-- Declare all global variables within the module
	158
	159	CHARACTER (LEN=80) :: filename_emis !< Variable for the name of the netcdf input file
	160
	161	INTEGER(iwp) :: i !< index 1st selected dimension (some dims are not spatial)
	162	INTEGER(iwp) :: j !< index 2nd selected dimension
	163	INTEGER(iwp) :: i_start !< Index to start read variable from netcdf along one dims
	164	INTEGER(iwp) :: i_end !< Index to end read variable from netcdf in one dims
	165	INTEGER(iwp) :: j_start !< Index to start read variable from netcdf in additional dims
	166	INTEGER(iwp) :: j_end !< Index to end read variable from netcdf in additional dims
	167	INTEGER(iwp) :: z_start !< Index to start read variable from netcdf in additional dims
	168	INTEGER(iwp) :: z_end !< Index to end read variable from netcdf in additional dims
	169	INTEGER(iwp) :: dt_emis !< Time Step Emissions
	170	INTEGER(iwp) :: len_index !< length of index (used for several indices)
	171	INTEGER(iwp) :: len_index_voc !< length of voc index
	172	INTEGER(iwp) :: len_index_pm !< length of PMs index
	173
	174	REAL(wp) :: con_factor !< Units Conversion Factor
	175
	176	REAL(wp), PARAMETER :: Rgas = 8.3144 !< gas constant in J/mol/K
	177	REAL(wp), PARAMETER :: pref_i = 1.0_wp / 100000.0_wp !< Inverse Reference Pressure (1/Pa)
	178	REAL(wp), PARAMETER :: r_cp = 0.286_wp !< R / cp (exponent for potential temperature)
	179
	180	SAVE
	181
	182	!
	183	!-- Checks Input parameters
	184	INTERFACE chem_emissions_check_parameters
	185	MODULE PROCEDURE chem_emissions_check_parameters
	186	END INTERFACE chem_emissions_check_parameters
	187	!
	188	!-- Matching Emissions actions
	189	INTERFACE chem_emissions_match
	190	MODULE PROCEDURE chem_emissions_match
	191	END INTERFACE chem_emissions_match
	192	!
	193	!-- Initialization actions
	194	INTERFACE chem_emissions_init
	195	MODULE PROCEDURE chem_emissions_init
	196	END INTERFACE chem_emissions_init
	197	!
	198	!-- Setup of Emissions
	199	INTERFACE chem_emissions_setup
	200	MODULE PROCEDURE chem_emissions_setup
	201	END INTERFACE chem_emissions_setup
	202
	203	PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup
	204	!
	205	!-- Public Variables
	206	PUBLIC con_factor, len_index, len_index_pm, len_index_voc
	207
	208	CONTAINS
	209
	210	!------------------------------------------------------------------------------!
	211	! Description:
	212	! ------------
	213	!> Routine for checking input parameters
	214	!------------------------------------------------------------------------------!
	215	SUBROUTINE chem_emissions_check_parameters
	216
	217	IMPLICIT NONE
	218
	219	TYPE(chem_emis_att_type) :: emt
	220
	221	!
	222	!-- Check if species count matches the number of names
	223	!-- passed for the chemiscal species
	224
	225	IF ( SIZE(emt%species_name) /= emt%n_emiss_species ) THEN
	226	! IF ( SIZE(emt%species_name) /= emt%nspec ) THEN
	227
	228	message_string = 'Numbers of input emission species names and number of species' // &
	229	'for which emission values are given do not match'
	230	CALL message( 'chem_emissions_check_parameters', 'CM0437', 2, 2, 0, 6, 0 )
	231
	232	ENDIF
	233
	234	END SUBROUTINE chem_emissions_check_parameters
	235
	236
	237	!------------------------------------------------------------------------------!
	238	! Description:
	239	! ------------
	240	!> Matching the chemical species indices. The routine checks what are the
	241	!> indices of the emission input species and the corresponding ones of the
	242	!> model species. The routine gives as output a vector containing the number
	243	!> of common species: it is important to note that while the model species
	244	!> are distinct, their values could be given to a single species in input.
	245	!> For example, in the case of NO2 and NO, values may be passed in input as
	246	!> NOX values.
	247	!------------------------------------------------------------------------------!
	248
	249	SUBROUTINE chem_emissions_match( emt_att,len_index )
	250
	251
	252
	253	INTEGER(iwp) :: ind_inp !< Parameters for cycling through chemical input species
	254	INTEGER(iwp) :: ind_mod !< Parameters for cycling through chemical model species
	255	INTEGER(iwp) :: ind_voc !< Indices to check whether a split for voc should be done
	256	INTEGER(iwp) :: ispec !< index for cycle over effective number of emission species
	257	INTEGER(iwp) :: nspec_emis_inp !< Variable where to store # of emission species in input
	258
	259	INTEGER(iwp), INTENT(INOUT) :: len_index !< number of common species between input dataset & model species
	260
	261	TYPE(chem_emis_att_type), INTENT(INOUT) :: emt_att !< Chemistry Emission Array (decl. netcdf_data_input.f90)
	262
	263
	264	IF ( debug_output ) CALL debug_message( 'chem_emissions_match', 'start' )
	265
	266	!
	267	!-- Number of input emission species
	268
	269	nspec_emis_inp = emt_att%n_emiss_species
	270	! nspec_emis_inp=emt_att%nspec
	271
	272	!
	273	!-- Check the emission LOD: 0 (PARAMETERIZED), 1 (DEFAULT), 2 (PRE-PROCESSED)
	274	!
	275	SELECT CASE (emiss_lod)
	276
	277	!
	278	!-- LOD 0 (PARAMETERIZED mode)
	279
	280	CASE (0)
	281
	282	len_index = 0
	283
	284	IF ( nvar > 0 .AND. nspec_emis_inp > 0 ) THEN
	285
	286	!
	287	!-- Cycle over model species
	288
	289	DO ind_mod = 1, nvar
	290	ind_inp = 1
	291	DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' ) !< 'novalue' is the default
	292	IF ( TRIM( surface_csflux_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
	293	len_index = len_index + 1
	294	ENDIF
	295	ind_inp = ind_inp + 1
	296	ENDDO
	297	ENDDO
	298
	299	IF ( len_index > 0 ) THEN
	300
	301	!
	302	!-- Allocation of Arrays of the matched species
	303
	304	ALLOCATE ( match_spec_input(len_index) )
	305	ALLOCATE ( match_spec_model(len_index) )
	306
	307	!
	308	!-- Pass species indices to declared arrays
	309
	310	len_index = 0
	311
	312	DO ind_mod = 1, nvar
	313	ind_inp = 1
	314	DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )
	315	IF ( TRIM(surface_csflux_name(ind_inp)) == &
	316	TRIM(spc_names(ind_mod)) ) THEN
	317	len_index = len_index + 1
	318	match_spec_input(len_index) = ind_inp
	319	match_spec_model(len_index) = ind_mod
	320	ENDIF
	321	ind_inp = ind_inp + 1
	322	END DO
	323	END DO
	324
	325	!
	326	!-- Check
	327
	328	DO ispec = 1, len_index
	329
	330	IF ( emiss_factor_main(match_spec_input(ispec) ) < 0 .AND. &
	331	emiss_factor_side(match_spec_input(ispec) ) < 0 ) THEN
	332
	333	message_string = 'PARAMETERIZED emissions mode selected:' // &
	334	' EMISSIONS POSSIBLE ONLY ON STREET SURFACES' // &
	335	' but values of scaling factors for street types' // &
	336	' emiss_factor_main AND emiss_factor_side' // &
	337	' not provided for each of the emissions species' // &
	338	' or not provided at all: PLEASE set a finite value' // &
	339	' for these parameters in the chemistry namelist'
	340	CALL message( 'chem_emissions_matching', 'CM0442', 2, 2, 0, 6, 0 )
	341
	342	ENDIF
	343
	344	END DO
	345
	346
	347	ELSE
	348
	349	message_string = 'Non of given Emission Species' // &
	350	' matches' // &
	351	' model chemical species' // &
	352	' Emission routine is not called'
	353	CALL message( 'chem_emissions_matching', 'CM0443', 0, 0, 0, 6, 0 )
	354
	355	ENDIF
	356
	357	ELSE
	358
	359	message_string = 'Array of Emission species not allocated: ' // &
	360	' Either no emission species are provided as input or' // &
	361	' no chemical species are used by PALM.' // &
	362	' Emission routine is not called'
	363	CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 )
	364
	365	ENDIF
	366
	367	!
	368	!-- LOD 1 (DEFAULT mode)
	369
	370	CASE (1)
	371
	372	len_index = 0 ! TOTAL number of species (to be accumulated)
	373	len_index_voc = 0 ! TOTAL number of VOCs (to be accumulated)
	374	len_index_pm = 3 ! TOTAL number of PMs: PM1, PM2.5, PM10.
	375
	376	IF ( nvar > 0 .AND. nspec_emis_inp > 0 ) THEN
	377
	378	!
	379	!-- Cycle over model species
	380	DO ind_mod = 1, nvar
	381
	382	!
	383	!-- Cycle over input species
	384
	385	DO ind_inp = 1, nspec_emis_inp
	386
	387	!
	388	!-- Check for VOC Species
	389
	390	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" ) THEN
	391	DO ind_voc= 1, emt_att%nvoc
	392
	393	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) ) THEN
	394	len_index = len_index + 1
	395	len_index_voc = len_index_voc + 1
	396	ENDIF
	397
	398	END DO
	399	ENDIF
	400
	401	!
	402	!-- PMs: There is one input species name for all PM
	403	!-- This variable has 3 dimensions, one for PM1, PM2.5 and PM10
	404
	405	IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" ) THEN
	406
	407	IF ( TRIM( spc_names(ind_mod) ) == "PM1" ) THEN
	408	len_index = len_index + 1
	409	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" ) THEN
	410	len_index = len_index + 1
	411	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" ) THEN
	412	len_index = len_index + 1
	413	ENDIF
	414
	415	ENDIF
	416
	417	!
	418	!-- NOX: NO2 and NO
	419
	420	IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOX" ) THEN
	421
	422	IF ( TRIM( spc_names(ind_mod) ) == "NO" ) THEN
	423	len_index = len_index + 1
	424	ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" ) THEN
	425	len_index = len_index + 1
	426	ENDIF
	427
	428	ENDIF
	429
	430	!
	431	!-- SOX: SO2 and SO4
	432
	433	IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOX" ) THEN
	434
	435	IF ( TRIM( spc_names(ind_mod) ) == "SO2" ) THEN
	436	len_index = len_index + 1
	437	ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" ) THEN
	438	len_index = len_index + 1
	439	ENDIF
	440
	441	ENDIF
	442
	443	!
	444	!-- Other Species
	445
	446	IF ( TRIM( emt_att%species_name(ind_inp) ) == &
	447	TRIM( spc_names(ind_mod) ) ) THEN
	448	len_index = len_index + 1
	449	ENDIF
	450
	451	END DO ! ind_inp ...
	452
	453	END DO ! ind_mod ...
	454
	455
	456	!
	457	!-- Allocate arrays
	458
	459	IF ( len_index > 0 ) THEN
	460
	461	ALLOCATE ( match_spec_input(len_index) )
	462	ALLOCATE ( match_spec_model(len_index) )
	463
	464	IF ( len_index_voc > 0 ) THEN
	465
	466	!
	467	!-- Contains indices of the VOC model species
	468
	469	ALLOCATE( match_spec_voc_model(len_index_voc) )
	470
	471	!
	472	!-- Contains the indices of different values of VOC composition
	473	!-- of input variable VOC_composition
	474
	475	ALLOCATE( match_spec_voc_input(len_index_voc) )
	476
	477	ENDIF
	478
	479	!
	480	!-- Pass the species indices to declared arrays
	481
	482	len_index = 0
	483	len_index_voc = 0
	484
	485	DO ind_mod = 1, nvar
	486	DO ind_inp = 1, nspec_emis_inp
	487
	488	!
	489	!-- VOCs
	490
	491	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" .AND. &
	492	ALLOCATED (match_spec_voc_input) ) THEN
	493
	494	DO ind_voc = 1, emt_att%nvoc
	495
	496	IF ( TRIM( emt_att%voc_name(ind_voc) ) == &
	497	TRIM( spc_names(ind_mod) ) ) THEN
	498
	499	len_index = len_index + 1
	500	len_index_voc = len_index_voc + 1
	501
	502	match_spec_input(len_index) = ind_inp
	503	match_spec_model(len_index) = ind_mod
	504
	505	match_spec_voc_input(len_index_voc) = ind_voc
	506	match_spec_voc_model(len_index_voc) = ind_mod
	507
	508	ENDIF
	509
	510	END DO
	511
	512	ENDIF
	513
	514	!
	515	!-- PMs
	516
	517	IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" ) THEN
	518
	519	IF ( TRIM( spc_names(ind_mod) ) == "PM1" ) THEN
	520	len_index = len_index + 1
	521	match_spec_input(len_index) = ind_inp
	522	match_spec_model(len_index) = ind_mod
	523	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" ) THEN
	524	len_index = len_index + 1
	525	match_spec_input(len_index) = ind_inp
	526	match_spec_model(len_index) = ind_mod
	527	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" ) THEN
	528	len_index = len_index + 1
	529	match_spec_input(len_index) = ind_inp
	530	match_spec_model(len_index) = ind_mod
	531	ENDIF
	532
	533	ENDIF
	534
	535	!
	536	!-- NOX
	537	IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOX" ) THEN
	538
	539	IF ( TRIM( spc_names(ind_mod) ) == "NO" ) THEN
	540	len_index = len_index + 1
	541
	542	match_spec_input(len_index) = ind_inp
	543	match_spec_model(len_index) = ind_mod
	544
	545	ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" ) THEN
	546	len_index = len_index + 1
	547
	548	match_spec_input(len_index) = ind_inp
	549	match_spec_model(len_index) = ind_mod
	550
	551	ENDIF
	552
	553	ENDIF
	554
	555
	556	!
	557	!-- SOX
	558
	559	IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOX" ) THEN
	560
	561	IF ( TRIM( spc_names(ind_mod) ) == "SO2" ) THEN
	562	len_index = len_index + 1
	563	match_spec_input(len_index) = ind_inp
	564	match_spec_model(len_index) = ind_mod
	565	ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" ) THEN
	566	len_index = len_index + 1
	567	match_spec_input(len_index) = ind_inp
	568	match_spec_model(len_index) = ind_mod
	569	ENDIF
	570
	571	ENDIF
	572
	573	!
	574	!-- Other Species
	575
	576	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) 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	END DO ! inp_ind
	583
	584	END DO ! inp_mod
	585
	586	!
	587	!-- Error reporting (no matching)
	588
	589	ELSE
	590
	591	message_string = 'None of given Emission Species matches' // &
	592	' model chemical species' // &
	593	' Emission routine is not called'
	594	CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 )
	595
	596	ENDIF
	597
	598	!
	599	!-- Error reporting (no species)
	600
	601	ELSE
	602
	603	message_string = 'Array of Emission species not allocated: ' // &
	604	' Either no emission species are provided as input or' // &
	605	' no chemical species are used by PALM:' // &
	606	' Emission routine is not called'
	607	CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 )
	608
	609	ENDIF
	610
	611	!
	612	!-- LOD 2 (PRE-PROCESSED mode)
	613
	614	CASE (2)
	615
	616	len_index = 0
	617	len_index_voc = 0
	618
	619	IF ( nvar > 0 .AND. (nspec_emis_inp > 0) ) THEN
	620	!
	621	!-- Cycle over model species
	622	DO ind_mod = 1, nvar
	623
	624	!
	625	!-- Cycle over input species
	626	DO ind_inp = 1, nspec_emis_inp
	627
	628	!
	629	!-- Check for VOC Species
	630
	631	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" ) THEN
	632	DO ind_voc = 1, emt_att%nvoc
	633	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) ) THEN
	634	len_index = len_index + 1
	635	len_index_voc = len_index_voc + 1
	636	ENDIF
	637	END DO
	638	ENDIF
	639
	640	!
	641	!-- Other Species
	642
	643	IF ( TRIM(emt_att%species_name(ind_inp)) == TRIM(spc_names(ind_mod)) ) THEN
	644	len_index = len_index + 1
	645	ENDIF
	646	ENDDO
	647	ENDDO
	648
	649	!
	650	!-- Allocate array for storing the indices of the matched species
	651
	652	IF ( len_index > 0 ) THEN
	653
	654	ALLOCATE ( match_spec_input(len_index) )
	655
	656	ALLOCATE ( match_spec_model(len_index) )
	657
	658	IF ( len_index_voc > 0 ) THEN
	659	!
	660	!-- contains indices of the VOC model species
	661	ALLOCATE( match_spec_voc_model(len_index_voc) )
	662	!
	663	!-- contains the indices of different values of VOC composition of input variable VOC_composition
	664	ALLOCATE( match_spec_voc_input(len_index_voc) )
	665
	666	ENDIF
	667
	668	!
	669	!-- pass the species indices to declared arrays
	670
	671	len_index = 0
	672
	673	!
	674	!-- Cycle over model species
	675
	676	DO ind_mod = 1, nvar
	677
	678	!
	679	!-- Cycle over Input species
	680
	681	DO ind_inp = 1, nspec_emis_inp
	682
	683	!
	684	!-- VOCs
	685
	686	IF ( TRIM(emt_att%species_name(ind_inp) ) == "VOC" .AND. &
	687	ALLOCATED(match_spec_voc_input) ) THEN
	688
	689	DO ind_voc= 1, emt_att%nvoc
	690	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) ) THEN
	691	len_index = len_index + 1
	692	len_index_voc = len_index_voc + 1
	693
	694	match_spec_input(len_index) = ind_inp
	695	match_spec_model(len_index) = ind_mod
	696
	697	match_spec_voc_input(len_index_voc) = ind_voc
	698	match_spec_voc_model(len_index_voc) = ind_mod
	699	ENDIF
	700	END DO
	701	ENDIF
	702
	703	!
	704	!-- Other Species
	705
	706	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
	707	len_index = len_index + 1
	708	match_spec_input(len_index) = ind_inp
	709	match_spec_model(len_index) = ind_mod
	710	ENDIF
	711
	712	END DO ! ind_inp
	713	END DO ! ind_mod
	714
	715	ELSE ! if len_index_voc .le. 0
	716
	717	!
	718	!-- in case there are no species matching
	719
	720	message_string = 'Non of given emission species' // &
	721	' matches' // &
	722	' model chemical species:' // &
	723	' Emission routine is not called'
	724	CALL message( 'chem_emissions_matching', 'CM0438', 0, 0, 0, 6, 0 )
	725	ENDIF
	726
	727	!
	728	!-- Error check (no matching)
	729
	730	ELSE
	731
	732	!
	733	!-- either spc_names is zero or nspec_emis_inp is not allocated
	734	message_string = 'Array of Emission species not allocated:' // &
	735	' Either no emission species are provided as input or' // &
	736	' no chemical species are used by PALM:' // &
	737	' Emission routine is not called'
	738	CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 )
	739
	740	ENDIF
	741
	742	!
	743	!-- If emission module is switched on but mode_emis is not specified or it is given the wrong name
	744
	745	!
	746	!-- Error check (no species)
	747
	748	CASE DEFAULT
	749
	750	message_string = 'Emission Module switched ON, but' // &
	751	' either no emission mode specified or incorrectly given :' // &
	752	' please, pass the correct value to the namelist parameter "mode_emis"'
	753	CALL message( 'chem_emissions_matching', 'CM0445', 2, 2, 0, 6, 0 )
	754
	755	END SELECT
	756
	757	IF ( debug_output ) CALL debug_message( 'chem_emissions_match', 'end' )
	758
	759	END SUBROUTINE chem_emissions_match
	760
	761
	762	!------------------------------------------------------------------------------!
	763	! Description:
	764	! ------------
	765	!> Initialization:
	766	!> Netcdf reading, arrays allocation and first calculation of cssws
	767	!> fluxes at timestep 0
	768	!------------------------------------------------------------------------------!
	769	SUBROUTINE chem_emissions_init
	770
	771	USE netcdf_data_input_mod, &
	772	ONLY: chem_emis, chem_emis_att
	773
	774	IMPLICIT NONE
	775
	776	INTEGER(iwp) :: ispec !< running index
	777
	778	!
	779	!-- Actions for initial runs
	780	! IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	781	!-- ...
	782	!
	783	!
	784	!-- Actions for restart runs
	785	! ELSE
	786	!-- ...
	787	!
	788	! ENDIF
	789
	790
	791	IF ( debug_output ) CALL debug_message( 'chem_emissions_init', 'start' )
	792
	793	!
	794	!-- Matching
	795
	796	CALL chem_emissions_match( chem_emis_att, n_matched_vars )
	797
	798	IF ( n_matched_vars == 0 ) THEN
	799
	800	emission_output_required = .FALSE.
	801
	802	ELSE
	803
	804	emission_output_required = .TRUE.
	805
	806
	807	!
	808	!-- Set molecule masses (in kg/mol)
	809
	810	ALLOCATE( chem_emis_att%xm(n_matched_vars) )
	811
	812	DO ispec = 1, n_matched_vars
	813	SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
	814	CASE ( 'SO2' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 2
	815	CASE ( 'SO4' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 4
	816	CASE ( 'NO' ); chem_emis_att%xm(ispec) = xm_N + xm_O
	817	CASE ( 'NO2' ); chem_emis_att%xm(ispec) = xm_N + xm_O * 2
	818	CASE ( 'NH3' ); chem_emis_att%xm(ispec) = xm_N + xm_H * 3
	819	CASE ( 'CO' ); chem_emis_att%xm(ispec) = xm_C + xm_O
	820	CASE ( 'CO2' ); chem_emis_att%xm(ispec) = xm_C + xm_O * 2
	821	CASE ( 'CH4' ); chem_emis_att%xm(ispec) = xm_C + xm_H * 4
	822	CASE ( 'HNO3' ); chem_emis_att%xm(ispec) = xm_H + xm_N + xm_O*3
	823	CASE DEFAULT
	824	chem_emis_att%xm(ispec) = 1.0_wp
	825	END SELECT
	826	ENDDO
	827
	828
	829	!
	830	!-- Get emissions for the first time step base on LOD (if defined)
	831	!-- or emission mode (if no LOD defined)
	832
	833	!
	834	!-- NOTE - I could use a combined if ( lod = xxx .or. mode = 'XXX' )
	835	!-- type of decision structure but I think it is much better
	836	!-- to implement it this way (i.e., conditional on lod if it
	837	!-- is defined, and mode if not) as we can easily take out
	838	!-- the case structure for mode_emis later on.
	839	!-- (ecc 20140424)
	840
	841	IF ( emiss_lod < 0 ) THEN !-- no LOD defined (not likely)
	842
	843	SELECT CASE ( TRIM( mode_emis ) )
	844
	845	CASE ( 'PARAMETERIZED' ) ! LOD 0
	846
	847	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
	848	ALLOCATE( emis_distribution(1,0:ny,0:nx,n_matched_vars) )
	849	ENDIF
	850
	851	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars)
	852
	853	CASE ( 'DEFAULT' ) ! LOD 1
	854
	855	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
	856	ALLOCATE( emis_distribution(1,0:ny,0:nx,n_matched_vars) )
	857	ENDIF
	858
	859	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
	860
	861	CASE ( 'PRE-PROCESSED' ) ! LOD 2
	862
	863	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
	864	ALLOCATE( emis_distribution(nzb:nzt+1,0:ny,0:nx,n_matched_vars) )
	865	ENDIF
	866
	867	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
	868
	869	END SELECT
	870
	871	ELSE ! if LOD is defined
	872
	873	SELECT CASE ( emiss_lod )
	874
	875	CASE ( 0 ) ! parameterized mode
	876
	877	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
	878	ALLOCATE( emis_distribution(1,0:ny,0:nx,n_matched_vars) )
	879	ENDIF
	880
	881	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars)
	882
	883	CASE ( 1 ) ! default mode
	884
	885	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
	886	ALLOCATE( emis_distribution(1,0:ny,0:nx,n_matched_vars) )
	887	ENDIF
	888
	889	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
	890
	891	CASE ( 2 ) ! pre-processed mode
	892
	893	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
	894	ALLOCATE( emis_distribution(nzb:nzt+1,0:ny,0:nx,n_matched_vars) )
	895	ENDIF
	896
	897	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
	898
	899	END SELECT
	900
	901	ENDIF
	902
	903	ENDIF
	904
	905	IF ( debug_output ) CALL debug_message( 'chem_emissions_init', 'end' )
	906
	907	END SUBROUTINE chem_emissions_init
	908
	909
	910
	911	!------------------------------------------------------------------------------!
	912	! Description:
	913	! ------------
	914	!> Routine for Update of Emission values at each timestep
	915	!-------------------------------------------------------------------------------!
	916
	917	SUBROUTINE chem_emissions_setup( emt_att, emt, n_matched_vars )
	918
	919	USE surface_mod, &
	920	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
	921	USE netcdf_data_input_mod, &
	922	ONLY: street_type_f
	923	USE arrays_3d, &
	924	ONLY: hyp, pt
	925
	926
	927	IMPLICIT NONE
	928
	929
	930	TYPE(chem_emis_att_type), INTENT(INOUT) :: emt_att !< variable to store emission information
	931
	932	TYPE(chem_emis_val_type), INTENT(INOUT), ALLOCATABLE, DIMENSION(:) :: emt !< variable to store emission input values,
	933	!< depending on the considered species
	934
	935	INTEGER,INTENT(IN) :: n_matched_vars !< Output of matching routine with number
	936	!< of matched species
	937
	938	INTEGER(iwp) :: i !< running index for grid in x-direction
	939	INTEGER(iwp) :: j !< running index for grid in y-direction
	940	INTEGER(iwp) :: k !< running index for grid in z-direction
	941	INTEGER(iwp) :: m !< running index for horizontal surfaces
	942
	943	INTEGER(iwp) :: icat !< Index for number of categories
	944	INTEGER(iwp) :: ispec !< index for number of species
	945	INTEGER(iwp) :: i_pm_comp !< index for number of PM components
	946	INTEGER(iwp) :: ivoc !< Index for number of VOCs
	947
	948	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: delta_emis
	949	REAL(wp), ALLOCATABLE, DIMENSION(:) :: time_factor !< factor for time scaling of emissions
	950	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: emis
	951
	952	REAL(wp), DIMENSION(24) :: par_emis_time_factor !< time factors for the parameterized mode:
	953	!< fixed houlry profile for example day
	954	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: conv_to_ratio !< factor used for converting input
	955	!< to concentration ratio
	956	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: tmp_temp
	957
	958	!
	959	!-- CONVERSION FACTORS: TIME
	960	REAL(wp), PARAMETER :: s_per_hour = 3600.0 !< number of sec per hour (s)/(hour)
	961	REAL(wp), PARAMETER :: s_per_day = 86400.0 !< number of sec per day (s)/(day)
	962	REAL(wp), PARAMETER :: hour_per_year = 8760.0 !< number of hours in a year of 365 days
	963	REAL(wp), PARAMETER :: hour_per_day = 24.0 !< number of hours in a day
	964
	965	REAL(wp), PARAMETER :: hour_to_s = 1/s_per_hour !< conversion from hours to seconds (s/hour) ~ 0.2777778
	966	REAL(wp), PARAMETER :: day_to_s = 1/s_per_day !< conversion from day to seconds (s/day) ~ 1.157407e-05
	967	REAL(wp), PARAMETER :: year_to_s = 1/(s_per_hour*hour_per_year) !< conversion from year to sec (s/year) ~ 3.170979e-08
	968	!
	969	!-- CONVERSION FACTORS: WEIGHT
	970	REAL(wp), PARAMETER :: tons_to_kg = 100 !< Conversion from tons to kg (kg/tons)
	971	REAL(wp), PARAMETER :: g_to_kg = 0.001 !< Conversion from g to kg (kg/g)
	972	REAL(wp), PARAMETER :: miug_to_kg = 0.000000001 !< Conversion from g to kg (kg/g)
	973	!
	974	!-- CONVERSION FACTORS: fraction to ppm
	975	REAL(wp), PARAMETER :: ratio2ppm = 1.0e06
	976	!------------------------------------------------------
	977
	978	IF ( emission_output_required ) THEN
	979
	980	!
	981	!-- Set emis_dt
	982	IF ( call_chem_at_all_substeps ) THEN
	983
	984	dt_emis = dt_3d * weight_pres(intermediate_timestep_count)
	985
	986	ELSE
	987
	988	dt_emis = dt_3d
	989
	990	ENDIF
	991
	992	!
	993	!-- Conversion of units to the ones employed in PALM
	994	!-- In PARAMETERIZED mode no conversion is performed: in this case input units are fixed
	995
	996	IF ( TRIM( mode_emis ) == "DEFAULT" .OR. TRIM( mode_emis ) == "PRE-PROCESSED" ) THEN
	997
	998	SELECT CASE ( TRIM( emt_att%units ) )
	999
	1000	CASE ( 'kg/m2/s', 'KG/M2/S' ); con_factor = 1.0_wp ! kg
	1001	CASE ( 'kg/m2/hour', 'KG/M2/HOUR' ); con_factor = hour_to_s
	1002	CASE ( 'kg/m2/day', 'KG/M2/DAY' ); con_factor = day_to_s
	1003	CASE ( 'kg/m2/year', 'KG/M2/YEAR' ); con_factor = year_to_s
	1004
	1005	CASE ( 'ton/m2/s', 'TON/M2/S' ); con_factor = tons_to_kg ! tonnes
	1006	CASE ( 'ton/m2/hour', 'TON/M2/HOUR' ); con_factor = tons_to_kg*hour_to_s
	1007	CASE ( 'ton/m2/year', 'TON/M2/YEAR' ); con_factor = tons_to_kg*year_to_s
	1008
	1009	CASE ( 'g/m2/s', 'G/M2/S' ); con_factor = g_to_kg ! grams
	1010	CASE ( 'g/m2/hour', 'G/M2/HOUR' ); con_factor = g_to_kg*hour_to_s
	1011	CASE ( 'g/m2/year', 'G/M2/YEAR' ); con_factor = g_to_kg*year_to_s
	1012
	1013	CASE ( 'micrograms/m2/s', 'MICROGRAMS/M2/S' ); con_factor = miug_to_kg ! ug
	1014	CASE ( 'micrograms/m2/hour', 'MICROGRAMS/M2/HOUR' ); con_factor = miug_to_kg*hour_to_s
	1015	CASE ( 'micrograms/m2/year', 'MICROGRAMS/M2/YEAR' ); con_factor = miug_to_kg*year_to_s
	1016
	1017	!
	1018	!-- Error check (need units)
	1019
	1020	CASE DEFAULT
	1021	message_string = 'The Units of the provided emission input' // &
	1022	' are not the ones required by PALM-4U: please check ' // &
	1023	' emission module documentation.'
	1024	CALL message( 'chem_emissions_setup', 'CM0446', 2, 2, 0, 6, 0 )
	1025
	1026	END SELECT
	1027
	1028	ENDIF
	1029
	1030	!
	1031	!-- Conversion factor to convert kg/m**2/s to ppm/s
	1032
	1033	DO i = nxl, nxr
	1034	DO j = nys, nyn
	1035
	1036	!
	1037	!-- Derive Temperature from Potential Temperature
	1038
	1039	tmp_temp(nzb:nzt+1,j,i) = pt(nzb:nzt+1,j,i) * ( hyp(nzb:nzt+1) * pref_i )**r_cp
	1040
	1041	!
	1042	!-- We need to pass to cssws <- (ppm/s) * dz
	1043	!-- Input is Nmole/(m^2*s)
	1044	!-- To go to ppmdz multiply the input by (m2/N)dz
	1045	!-- (m*2/N)dz == V/N
	1046	!-- V/N = RT/P
	1047
	1048	! m*3/Nmole (J/mol)K^-1 K Pa
	1049	conv_to_ratio(nzb:nzt+1,j,i) = ( (Rgas * tmp_temp(nzb:nzt+1,j,i)) / ((hyp(nzb:nzt+1))) )
	1050
	1051	ENDDO
	1052	ENDDO
	1053
	1054
	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	con_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	con_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	con_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	con_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	con_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	con_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	con_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	con_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	con_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	con_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	! con_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	IF ( street_type_f%var(j,i) >= main_street_id .AND. &
	1467	street_type_f%var(j,i) < max_street_id ) THEN
	1468
	1469	!
	1470	!-- Cycle over matched species
	1471
	1472	DO ispec = 1, n_matched_vars
	1473
	1474	!
	1475	!-- PMs are already in kilograms
	1476
	1477	IF ( TRIM(spc_names(match_spec_model(ispec))) == "PM1" .OR. &
	1478	TRIM(spc_names(match_spec_model(ispec))) == "PM25" .OR. &
	1479	TRIM(spc_names(match_spec_model(ispec))) == "PM10" ) THEN
	1480
	1481	!
	1482	!-- kg/(m^2s) kg/m^3
	1483	surf_lsm_h%cssws(match_spec_model(ispec),m) = &
	1484	emiss_factor_main(match_spec_input(ispec)) * &
	1485	emis_distribution(1,j,i,ispec) * & ! kg/(m^2*s)
	1486	rho_air(k) ! kg/m^3
	1487
	1488	!
	1489	!-- Other Species
	1490	!-- Inputs are micromoles
	1491
	1492	ELSE
	1493
	1494	!
	1495	!-- ppm/s m kg/m^3
	1496	surf_lsm_h%cssws(match_spec_model(ispec),m) = &
	1497	emiss_factor_main(match_spec_input(ispec)) * &
	1498	emis_distribution(1,j,i,ispec) * & ! micromoles/(m^2*s)
	1499	conv_to_ratio(k,j,i) * & ! m^3/Nmole
	1500	rho_air(k) ! kg/m^3
	1501
	1502	ENDIF
	1503
	1504	ENDDO ! ispec
	1505
	1506
	1507	ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. &
	1508	street_type_f%var(j,i) < main_street_id ) THEN
	1509
	1510	!
	1511	!-- Cycle over matched species
	1512
	1513	DO ispec = 1, n_matched_vars
	1514
	1515	!
	1516	!-- PMs are already in kilograms
	1517
	1518	IF ( TRIM(spc_names(match_spec_model(ispec))) == "PM1" .OR. &
	1519	TRIM(spc_names(match_spec_model(ispec))) == "PM25" .OR. &
	1520	TRIM(spc_names(match_spec_model(ispec))) == "PM10" ) THEN
	1521
	1522	!
	1523	!-- kg/(m^2s) kg/m^3
	1524	surf_lsm_h%cssws(match_spec_model(ispec),m) = &
	1525	emiss_factor_side(match_spec_input(ispec)) * &
	1526	emis_distribution(1,j,i,ispec) * & ! kg/(m^2*s)
	1527	rho_air(k) ! kg/m^3
	1528	!
	1529	!-- Other species
	1530	!-- Inputs are micromoles
	1531
	1532	ELSE
	1533
	1534	!
	1535	!-- ppm/s m kg/m^3
	1536
	1537	surf_lsm_h%cssws(match_spec_model(ispec),m) = &
	1538	emiss_factor_side(match_spec_input(ispec)) * &
	1539	emis_distribution(1,j,i,ispec) * & ! micromoles/(m^2*s)
	1540	conv_to_ratio(k,j,i) * & ! m^3/Nmole
	1541	rho_air(k) ! kg/m^3
	1542
	1543	ENDIF
	1544
	1545	ENDDO ! ispec
	1546
	1547	!
	1548	!-- If no street type is defined, then assign zero emission to all the species
	1549
	1550	ELSE
	1551
	1552	surf_lsm_h%cssws(:,m) = 0.0_wp
	1553
	1554	ENDIF ! street type
	1555
	1556	ENDDO ! m
	1557
	1558	ENDIF ! street_type_f%from_file
	1559
	1560
	1561	!
	1562	!-- LOD 1 (DEFAULT) and LOD 2 (PRE-PROCESSED)
	1563
	1564
	1565	ELSE
	1566
	1567
	1568	DO ispec = 1, n_matched_vars
	1569
	1570	!
	1571	!-- Default surfaces
	1572
	1573	DO m = 1, surf_def_h(0)%ns
	1574
	1575	i = surf_def_h(0)%i(m)
	1576	j = surf_def_h(0)%j(m)
	1577
	1578	IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN
	1579
	1580	!
	1581	!-- PMs
	1582	IF ( TRIM(spc_names(match_spec_model(ispec))) == "PM1" .OR. &
	1583	TRIM(spc_names(match_spec_model(ispec))) == "PM25" .OR. &
	1584	TRIM(spc_names(match_spec_model(ispec))) == "PM10" ) THEN
	1585
	1586	surf_def_h(0)%cssws(match_spec_model(ispec),m) = & ! kg/m2/s * kg/m3
	1587	emis_distribution(1,j,i,ispec)* & ! kg/m2/s
	1588	rho_air(nzb) ! kg/m^3
	1589
	1590	ELSE
	1591
	1592	!
	1593	!-- VOCs
	1594	IF ( len_index_voc > 0 .AND. &
	1595	emt_att%species_name(match_spec_input(ispec)) == "VOC" ) THEN
	1596
	1597	surf_def_h(0)%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
	1598	emis_distribution(1,j,i,ispec) * & ! mole/m2/s
	1599	conv_to_ratio(nzb,j,i) * & ! m^3/mole
	1600	ratio2ppm * & ! ppm
	1601	rho_air(nzb) ! kg/m^3
	1602
	1603
	1604	!
	1605	!-- Other species
	1606
	1607	ELSE
	1608
	1609	surf_def_h(0)%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
	1610	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
	1611	( 1.0_wp / emt_att%xm(ispec) ) * & ! mole/kg
	1612	conv_to_ratio(nzb,j,i) * & ! m^3/mole
	1613	ratio2ppm * & ! ppm
	1614	rho_air(nzb) ! kg/m^3
	1615
	1616	ENDIF ! VOC
	1617
	1618	ENDIF ! PM
	1619
	1620	ENDIF ! emis_distribution > 0
	1621
	1622	ENDDO ! m
	1623
	1624	!
	1625	!-- LSM surfaces
	1626
	1627	DO m = 1, surf_lsm_h%ns
	1628
	1629	i = surf_lsm_h%i(m)
	1630	j = surf_lsm_h%j(m)
	1631	k = surf_lsm_h%k(m)
	1632
	1633	IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN
	1634
	1635	!
	1636	!-- PMs
	1637	IF ( TRIM(spc_names(match_spec_model(ispec))) == "PM1" .OR. &
	1638	TRIM(spc_names(match_spec_model(ispec))) == "PM25" .OR. &
	1639	TRIM(spc_names(match_spec_model(ispec))) == "PM10" ) THEN
	1640
	1641	surf_lsm_h%cssws(match_spec_model(ispec),m) = & ! kg/m2/s * kg/m3
	1642	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
	1643	rho_air(k) ! kg/m^3
	1644
	1645	ELSE
	1646
	1647	!
	1648	!-- VOCs
	1649
	1650	IF ( len_index_voc > 0 .AND. &
	1651	emt_att%species_name(match_spec_input(ispec)) == "VOC" ) THEN
	1652
	1653	surf_lsm_h%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
	1654	emis_distribution(1,j,i,ispec) * & ! mole/m2/s
	1655	conv_to_ratio(k,j,i) * & ! m^3/mole
	1656	ratio2ppm * & ! ppm
	1657	rho_air(k) ! kg/m^3
	1658
	1659	!
	1660	!-- Other species
	1661
	1662	ELSE
	1663
	1664	surf_lsm_h%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
	1665	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
	1666	( 1.0_wp / emt_att%xm(ispec) ) * & ! mole/kg
	1667	conv_to_ratio(k,j,i) * & ! m^3/mole
	1668	ratio2ppm * & ! ppm
	1669	rho_air(k) ! kg/m^3
	1670
	1671	ENDIF ! VOC
	1672
	1673	ENDIF ! PM
	1674
	1675	ENDIF ! emis_distribution
	1676
	1677	ENDDO ! m
	1678
	1679	!
	1680	!-- USM surfaces
	1681
	1682	DO m = 1, surf_usm_h%ns
	1683
	1684	i = surf_usm_h%i(m)
	1685	j = surf_usm_h%j(m)
	1686	k = surf_usm_h%k(m)
	1687
	1688	IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN
	1689
	1690	!
	1691	!-- PMs
	1692	IF ( TRIM(spc_names(match_spec_model(ispec))) == "PM1" .OR. &
	1693	TRIM(spc_names(match_spec_model(ispec))) == "PM25" .OR. &
	1694	TRIM(spc_names(match_spec_model(ispec))) == "PM10" ) THEN
	1695
	1696	surf_usm_h%cssws(match_spec_model(ispec),m) = & ! kg/m2/s * kg/m3
	1697	emis_distribution(1,j,i,ispec)* & ! kg/m2/s
	1698	rho_air(k) ! kg/m^3
	1699
	1700	ELSE
	1701
	1702	!
	1703	!-- VOCs
	1704	IF ( len_index_voc > 0 .AND. &
	1705	emt_att%species_name(match_spec_input(ispec)) == "VOC" ) THEN
	1706
	1707	surf_usm_h%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
	1708	emis_distribution(1,j,i,ispec) * & ! m2/s
	1709	conv_to_ratio(k,j,i) * & ! m^3/mole
	1710	ratio2ppm * & ! ppm
	1711	rho_air(k) ! kg/m^3
	1712
	1713	!
	1714	!-- Other species
	1715	ELSE
	1716
	1717	surf_usm_h%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
	1718	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
	1719	( 1.0_wp / emt_att%xm(ispec) ) * & ! mole/kg
	1720	conv_to_ratio(k,j,i) * & ! m^3/mole
	1721	ratio2ppm* & ! ppm
	1722	rho_air(k) ! kg/m^3
	1723
	1724
	1725	ENDIF ! VOC
	1726
	1727	ENDIF ! PM
	1728
	1729	ENDIF ! emis_distribution
	1730
	1731	ENDDO ! m
	1732
	1733	ENDDO
	1734
	1735	ENDIF
	1736
	1737	!
	1738	!-- Deallocate time_factor in case of DEFAULT mode)
	1739
	1740	IF ( ALLOCATED (time_factor) ) DEALLOCATE (time_factor)
	1741
	1742	ENDIF
	1743
	1744	END SUBROUTINE chem_emissions_setup
	1745
	1746	END MODULE chem_emissions_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |