Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |