Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |