Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |