Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |