Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |