Home

Context Navigation

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

Last change on this file since 3924 was 3885, checked in by kanani, 6 years ago
restructure/add location/debug messages
Property svn:keywords set to `Id`
File size: 64.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |