Home

Context Navigation

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

Last change on this file since 4144 was 4144, checked in by raasch, 5 years ago
relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
Property svn:keywords set to `Id`
File size: 64.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |