Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |