Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |