Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |