Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |