Home

Context Navigation

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

Last change on this file since 4154 was 4154, checked in by suehring, 5 years ago
Bugfix in chemistry decycling; Replace global arrays to setup chemical emissions by local ones
Property svn:keywords set to `Id`
File size: 65.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |