Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |