Home

Context Navigation

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

Last change on this file since 3775 was 3772, checked in by suehring, 6 years ago
Assure that emissions are only at natural (pavement) surfaces in case of parametrized emission mode, not on urban surfaces
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-2018 Leibniz Universitaet Hannover
18	! Copyright 2018-2018 Freie Universitaet Berlin
19	! Copyright 2018-2018 Karlsruhe Institute of Technology
20	!--------------------------------------------------------------------------------!
21	!
22	! Current revisions:
23	! ------------------
24	!
25	!
26	! Former revisions:
27	! -----------------
28	! $Id: chem_emissions_mod.f90 3772 2019-02-28 15:51:57Z gronemeier $
29	! - In case of parametrized emissions, assure that emissions are only on natural
30	! surfaces (i.e. streets) and not on urban surfaces.
31	! - some unnecessary if clauses removed
32	!
33	! 3685 2019-01-21 01:02:11Z knoop
34	! Some interface calls moved to module_interface + cleanup
35	!
36	! 3611 2018-12-07 14:14:11Z banzhafs
37	! Code update to comply PALM coding rules
38	! removed unnecessary informative messsages/location
39	! messages and corrected comments on PM units from to kg
40	! bug fix: spcs_name replaced by nvar in DO loops
41	!
42	! 3591 2018-11-30 17:37:32Z suehring
43	! - Removed salsa dependency.
44	! - Enabled PARAMETRIZED mode for default surfaces when LSM is not applied but
45	! salsa is (M. Kurppa)
46	!
47	! 3582 2018-11-29 19:16:36Z suehring
48	! resler:
49	! Break lines at 132 characters
50	!
51	! 3483 2018-11-02 14:19:26Z raasch
52	! bugfix: wrong locations of netCDF directives fixed
53	!
54	! 3467 2018-10-30 19:05:21Z suehring
55	! Enabled PARAMETRIZED mode for default surfaces when LSM is not applied but
56	! salsa is used
57	!
58	! 3458 2018-10-30 14:51:23Z kanani
59	! from chemistry branch r3443, banzhafs, Russo:
60	! Additional correction for index of input file of pre-processed mode
61	! Removed atomic and molecular weights as now available in chem_modules and
62	! added par_emis_time_factor (formerly in netcdf_data_input_mod)
63	! Added correction for index of input file of pre-processed mode
64	! Added a condition for default mode necessary for information read from ncdf file
65	! in pre-processed and default mode
66	! Correction of multiplying factor necessary for scaling emission values in time
67	! Introduced correction for scaling units in the case of DEFAULT emission mode
68	!
69	! 3373 2018-10-18 15:25:56Z kanani
70	! Fix wrong location of __netcdf directive
71	!
72	! 3337 2018-10-12 15:17:09Z kanani
73	! (from branch resler)
74	! Formatting
75	!
76	! 3312 2018-10-06 14:15:46Z knoop
77	! Initial revision
78	!
79	! 3286 2018-09-28 07:41:39Z forkel
80	!
81	! Authors:
82	! --------
83	! @author Emmanuele Russo (FU-Berlin)
84	! @author Sabine Banzhaf (FU-Berlin)
85	! @author Martijn Schaap (FU-Berlin, TNO Utrecht)
86	!
87	! Description:
88	! ------------
89	!> MODULE for reading-in Chemistry Emissions data
90	!>
91	!> @todo Check_parameters routine should be developed: for now it includes just one condition
92	!> @todo Use of Restart files not contempled at the moment
93	!> @todo revise indices of files read from the netcdf: preproc_emission_data and expert_emission_data
94	!> @todo for now emission data may be passed on a singular vertical level: need to be more flexible
95	!> @todo fill/activate restart option in chem_emissions_init
96	!> @todo discuss dt_emis
97	!> @note <Enter notes on the module>
98	!> @bug <Enter known bugs here>
99	!------------------------------------------------------------------------------!
100
101	MODULE chem_emissions_mod
102
103	USE arrays_3d, &
104	ONLY: rho_air
105
106	USE control_parameters, &
107	ONLY: end_time, message_string, initializing_actions, &
108	intermediate_timestep_count, dt_3d
109
110	USE indices
111
112	USE kinds
113
114	#if defined( __netcdf )
115	USE netcdf
116	#endif
117
118	USE netcdf_data_input_mod, &
119	ONLY: chem_emis_att_type, chem_emis_val_type
120
121	USE date_and_time_mod, &
122	ONLY: day_of_month, hour_of_day, &
123	index_mm, index_dd, index_hh, &
124	month_of_year, hour_of_day, &
125	time_default_indices, time_preprocessed_indices
126
127	USE chem_gasphase_mod, &
128	ONLY: spc_names
129
130	USE chem_modules
131
132	USE statistics, &
133	ONLY: weight_pres
134
135
136	IMPLICIT NONE
137
138	!
139	!-- Declare all global variables within the module
140
141	CHARACTER (LEN=80) :: filename_emis !< Variable for the name of the netcdf input file
142
143	INTEGER(iwp) :: i !< index 1st selected dimension (some dims are not spatial)
144	INTEGER(iwp) :: j !< index 2nd selected dimension
145	INTEGER(iwp) :: i_start !< Index to start read variable from netcdf along one dims
146	INTEGER(iwp) :: i_end !< Index to end read variable from netcdf in one dims
147	INTEGER(iwp) :: j_start !< Index to start read variable from netcdf in additional dims
148	INTEGER(iwp) :: j_end !< Index to end read variable from netcdf in additional dims
149	INTEGER(iwp) :: 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	INTEGER(iwp) :: dt_emis !< Time Step Emissions
152	INTEGER(iwp) :: len_index !< length of index (used for several indices)
153	INTEGER(iwp) :: len_index_voc !< length of voc index
154	INTEGER(iwp) :: len_index_pm !< length of PMs index
155
156	REAL(wp) :: con_factor !< Units Conversion Factor
157
158	REAL(wp), PARAMETER :: Rgas = 8.3144 !< gas constant in J/mol/K
159	REAL(wp), PARAMETER :: pref_i = 1.0_wp / 100000.0_wp !< Inverse Reference Pressure (1/Pa)
160	REAL(wp), PARAMETER :: r_cp = 0.286_wp !< R / cp (exponent for potential temperature)
161
162
163	SAVE
164
165
166	!
167	!-- Checks Input parameters
168	INTERFACE chem_emissions_check_parameters
169	MODULE PROCEDURE chem_emissions_check_parameters
170	END INTERFACE chem_emissions_check_parameters
171	!
172	!-- Matching Emissions actions
173	INTERFACE chem_emissions_match
174	MODULE PROCEDURE chem_emissions_match
175	END INTERFACE chem_emissions_match
176	!
177	!-- Initialization actions
178	INTERFACE chem_emissions_init
179	MODULE PROCEDURE chem_emissions_init
180	END INTERFACE chem_emissions_init
181	!
182	!-- Setup of Emissions
183	INTERFACE chem_emissions_setup
184	MODULE PROCEDURE chem_emissions_setup
185	END INTERFACE chem_emissions_setup
186
187
188
189	PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup
190	!
191	!-- Public Variables
192	PUBLIC con_factor, len_index, len_index_pm, len_index_voc
193
194	CONTAINS
195
196	!------------------------------------------------------------------------------!
197	! Description:
198	! ------------
199	!> Routine for checking input parameters
200	!------------------------------------------------------------------------------!
201	SUBROUTINE chem_emissions_check_parameters
202
203
204	IMPLICIT NONE
205
206	INTEGER(iwp) :: tmp
207
208	TYPE(chem_emis_att_type) :: emt
209
210	!
211	!-- Check Emission Species Number equal to number of passed names for the chemistry species:
212	IF ( SIZE(emt%species_name) /= emt%nspec ) THEN
213
214	message_string = 'Numbers of input emission species names and number of species' // &
215	'for which emission values are given do not match'
216	CALL message( 'chem_emissions_check_parameters', 'CM0437', 2, 2, 0, 6, 0 )
217
218	ENDIF
219
220
221	END SUBROUTINE chem_emissions_check_parameters
222
223	!------------------------------------------------------------------------------!
224	! Description:
225	! ------------
226	!> Matching the chemical species indices. The routine checks what are the indices of the emission input species
227	!> and the corresponding ones of the model species. The routine gives as output a vector containing the number
228	!> of common species: it is important to note that while the model species are distinct, their values could be
229	!> given to a single species in input: for example, in the case of NO2 and NO, values may be passed in input as
230	!> NOx values.
231	!------------------------------------------------------------------------------!
232	SUBROUTINE chem_emissions_match( emt_att,len_index )
233
234
235	INTEGER(iwp), INTENT(INOUT) :: len_index !< Variable where to store the number of common species between the input dataset and the model species
236
237	TYPE(chem_emis_att_type), INTENT(INOUT) :: emt_att !< Chemistry Emission Array containing information for all the input chemical emission species
238
239	INTEGER(iwp) :: ind_mod, ind_inp !< Parameters for cycling through chemical model and input species
240	INTEGER(iwp) :: nspec_emis_inp !< Variable where to store the number of the emission species in input
241	INTEGER(iwp) :: ind_voc !< Indices to check whether a split for voc should be done
242	INTEGER(iwp) :: ispec !< index for cycle over effective number of emission species
243
244
245	CALL location_message( 'Matching input emissions and model chemistry species', .FALSE. )
246
247	!
248	!-- Number of input emission species.
249	nspec_emis_inp=emt_att%nspec
250
251	!
252	!-- Check the emission mode: DEFAULT, PRE-PROCESSED or PARAMETERIZED
253	SELECT CASE( TRIM( mode_emis ) )
254
255	!
256	!-- PRE-PROCESSED mode
257	CASE ( "PRE-PROCESSED" )
258
259	len_index = 0
260	len_index_voc = 0
261
262	IF ( nvar > 0 .AND. (nspec_emis_inp > 0) ) THEN
263	!
264	!-- Cycle over model species
265	DO ind_mod = 1, nvar
266	!
267	!-- Cycle over input species
268	DO ind_inp = 1, nspec_emis_inp
269
270	!
271	!-- Check for VOC Species
272	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" ) THEN
273	DO ind_voc = 1, emt_att%nvoc
274
275	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) ) THEN
276	len_index = len_index + 1
277	len_index_voc = len_index_voc + 1
278	ENDIF
279	END DO
280	ENDIF
281	!
282	!-- Other Species
283	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
284	len_index = len_index + 1
285	ENDIF
286	ENDDO
287	ENDDO
288
289	!
290	!-- Allocate array for storing the indices of the matched species
291	IF ( len_index > 0 ) THEN
292
293	ALLOCATE( match_spec_input(len_index) )
294
295	ALLOCATE( match_spec_model(len_index) )
296
297	IF ( len_index_voc > 0 ) THEN
298	!
299	!-- contains indices of the VOC model species
300	ALLOCATE( match_spec_voc_model(len_index_voc) )
301	!
302	!-- contains the indices of different values of VOC composition of input variable VOC_composition
303	ALLOCATE( match_spec_voc_input(len_index_voc) )
304
305	ENDIF
306
307	!
308	!-- pass the species indices to declared arrays
309	len_index = 0
310
311	!
312	!-- Cycle over model species
313	DO ind_mod = 1, nvar
314	!
315	!-- Cycle over Input species
316	DO ind_inp = 1, nspec_emis_inp
317	!
318	!-- VOCs
319	IF ( TRIM(emt_att%species_name(ind_inp) ) == "VOC" .AND. &
320	ALLOCATED(match_spec_voc_input) ) THEN
321
322	DO ind_voc= 1, emt_att%nvoc
323	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) ) THEN
324	len_index = len_index + 1
325	len_index_voc = len_index_voc + 1
326
327	match_spec_input(len_index) = ind_inp
328	match_spec_model(len_index) = ind_mod
329
330	match_spec_voc_input(len_index_voc) = ind_voc
331	match_spec_voc_model(len_index_voc) = ind_mod
332	ENDIF
333	END DO
334	ENDIF
335
336	!
337	!-- Other Species
338	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
339	len_index = len_index + 1
340	match_spec_input(len_index) = ind_inp
341	match_spec_model(len_index) = ind_mod
342	ENDIF
343	END DO
344	END DO
345
346	ELSE
347	!
348	!-- in case there are no species matching
349	message_string = 'Non of given emission species' // &
350	' matches' // &
351	' model chemical species:' // &
352	' Emission routine is not called'
353	CALL message( 'chem_emissions_matching', 'CM0438', 0, 0, 0, 6, 0 )
354	ENDIF
355
356	ELSE
357
358	!
359	!-- either spc_names is zero or nspec_emis_inp is not allocated
360	message_string = 'Array of Emission species not allocated:' // &
361	' Either no emission species are provided as input or' // &
362	' no chemical species are used by PALM:' // &
363	' Emission routine is not called'
364	CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 )
365
366	ENDIF
367
368	!
369	!-- DEFAULT mode
370	CASE ("DEFAULT")
371
372	len_index = 0 !< index for TOTAL number of species
373	len_index_voc = 0 !< index for TOTAL number of VOCs
374	len_index_pm = 3 !< index for TOTAL number of PMs: PM1, PM2.5, PM10.
375
376	IF ( nvar > 0 .AND. nspec_emis_inp > 0 ) THEN
377
378	!
379	!-- Cycle over model species
380	DO ind_mod = 1, nvar
381	!
382	!-- Cycle over input species
383	DO ind_inp = 1, nspec_emis_inp
384
385	!
386	!-- Check for VOC Species
387	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" ) THEN
388	DO ind_voc= 1, emt_att%nvoc
389
390	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) ) THEN
391	len_index = len_index + 1
392	len_index_voc = len_index_voc + 1
393	ENDIF
394
395	END DO
396	ENDIF
397
398	!
399	!-- PMs: There is one input species name for all PM
400	!-- This variable has 3 dimensions, one for PM1, PM2.5 and PM10
401	IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" ) THEN
402	!
403	!-- PM1
404	IF ( TRIM( spc_names(ind_mod) ) == "PM1" ) THEN
405	len_index = len_index + 1
406	!
407	!-- PM2.5
408	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" ) THEN
409	len_index = len_index + 1
410	!
411	!-- PM10
412	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" ) THEN
413	len_index = len_index + 1
414	ENDIF
415	ENDIF
416
417	!
418	!-- NOx: NO2 and NO
419	IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOx" ) THEN
420	!
421	!-- NO
422	IF ( TRIM( spc_names(ind_mod) ) == "NO" ) THEN
423	len_index = len_index + 1
424	!
425	!-- NO2
426	ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" ) THEN
427	len_index = len_index + 1
428	ENDIF
429	ENDIF
430
431	!
432	!-- SOx: SO2 and SO4
433	IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOx" ) THEN
434	!
435	!-- SO2
436	IF ( TRIM( spc_names(ind_mod) ) == "SO2" ) THEN
437	len_index = len_index + 1
438	!
439	!-- SO4
440	ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" ) THEN
441	len_index = len_index + 1
442	ENDIF
443	ENDIF
444
445	!
446	!-- Other Species
447	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
448	len_index = len_index + 1
449	ENDIF
450	END DO
451	END DO
452
453
454	!
455	!-- Allocate arrays
456	IF ( len_index > 0 ) THEN
457
458	ALLOCATE( match_spec_input(len_index) )
459	ALLOCATE( match_spec_model(len_index) )
460
461	IF ( len_index_voc > 0 ) THEN
462	!
463	!-- Contains indices of the VOC model species
464	ALLOCATE( match_spec_voc_model(len_index_voc) )
465	!
466	!-- Contains the indices of different values of VOC composition
467	!-- of input variable VOC_composition
468	ALLOCATE( match_spec_voc_input(len_index_voc) )
469	ENDIF
470
471	!
472	!-- Pass the species indices to declared arrays
473	len_index = 0
474	len_index_voc = 0
475
476	DO ind_mod = 1, nvar
477	DO ind_inp = 1, nspec_emis_inp
478
479	!
480	!-- VOCs
481	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" .AND. &
482	ALLOCATED(match_spec_voc_input) ) THEN
483	DO ind_voc= 1, emt_att%nvoc
484	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) ) THEN
485	len_index = len_index + 1
486	len_index_voc = len_index_voc + 1
487
488	match_spec_input(len_index) = ind_inp
489	match_spec_model(len_index) = ind_mod
490
491	match_spec_voc_input(len_index_voc) = ind_voc
492	match_spec_voc_model(len_index_voc) = ind_mod
493	ENDIF
494	END DO
495	ENDIF
496
497	!
498	!-- PMs
499	IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" ) THEN
500	!
501	!-- PM1
502	IF ( TRIM( spc_names(ind_mod) ) == "PM1" ) THEN
503	len_index = len_index + 1
504
505	match_spec_input(len_index) = ind_inp
506	match_spec_model(len_index) = ind_mod
507	!
508	!-- PM2.5
509	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" ) THEN
510	len_index = len_index + 1
511
512	match_spec_input(len_index) = ind_inp
513	match_spec_model(len_index) = ind_mod
514	!
515	!-- PM10
516	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" ) THEN
517	len_index = len_index + 1
518
519	match_spec_input(len_index) = ind_inp
520	match_spec_model(len_index) = ind_mod
521
522	ENDIF
523	ENDIF
524
525	!
526	!-- NOx
527	IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOx" ) THEN
528	!
529	!-- NO
530	IF ( TRIM( spc_names(ind_mod) ) == "NO" ) THEN
531	len_index = len_index + 1
532
533	match_spec_input(len_index) = ind_inp
534	match_spec_model(len_index) = ind_mod
535	!
536	!-- NO2
537	ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" ) THEN
538	len_index = len_index + 1
539
540	match_spec_input(len_index) = ind_inp
541	match_spec_model(len_index) = ind_mod
542
543	ENDIF
544	ENDIF
545
546	!
547	!-- SOx
548	IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOx" ) THEN
549	!
550	!-- SO2
551	IF ( TRIM( spc_names(ind_mod) ) == "SO2" ) THEN
552	len_index = len_index + 1
553
554	match_spec_input(len_index) = ind_inp
555	match_spec_model(len_index) = ind_mod
556
557	!
558	!-- SO4
559	ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" ) THEN
560	len_index = len_index + 1
561
562	match_spec_input(len_index) = ind_inp
563	match_spec_model(len_index) = ind_mod
564
565	ENDIF
566	ENDIF
567
568	!
569	!-- Other Species
570	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
571	len_index = len_index + 1
572
573	match_spec_input(len_index) = ind_inp
574	match_spec_model(len_index) = ind_mod
575	ENDIF
576	END DO
577	END DO
578
579	ELSE
580
581	message_string = 'Non of given Emission Species' // &
582	' matches' // &
583	' model chemical species' // &
584	' Emission routine is not called'
585	CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 )
586
587	ENDIF
588
589	ELSE
590
591	message_string = 'Array of Emission species not allocated: ' // &
592	' Either no emission species are provided as input or' // &
593	' no chemical species are used by PALM:' // &
594	' Emission routine is not called'
595	CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 )
596
597	ENDIF
598
599	!
600	!-- PARAMETERIZED mode
601	CASE ("PARAMETERIZED")
602
603	len_index = 0
604
605	IF ( nvar > 0 .AND. nspec_emis_inp > 0 ) THEN
606
607	!
608	!-- Cycle over model species
609	DO ind_mod = 1, nvar
610	ind_inp = 1
611	DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' ) !< 'novalue' is the default
612	IF ( TRIM( surface_csflux_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
613	len_index = len_index + 1
614	ENDIF
615	ind_inp = ind_inp + 1
616	ENDDO
617	ENDDO
618
619	IF ( len_index > 0 ) THEN
620
621	!
622	!-- Allocation of Arrays of the matched species
623	ALLOCATE( match_spec_input(len_index) )
624
625	ALLOCATE( match_spec_model(len_index) )
626
627	!
628	!-- Pass the species indices to declared arrays
629	len_index = 0
630
631	DO ind_mod = 1, nvar
632	ind_inp = 1
633	DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )
634	IF ( TRIM( surface_csflux_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
635	len_index = len_index + 1
636	match_spec_input(len_index) = ind_inp
637	match_spec_model(len_index) = ind_mod
638	ENDIF
639	ind_inp = ind_inp + 1
640	END DO
641	END DO
642
643	!
644	!-- Check
645	DO ispec = 1, len_index
646
647	IF ( emiss_factor_main(match_spec_input(ispec) ) < 0 .AND. &
648	emiss_factor_side(match_spec_input(ispec) ) < 0 ) THEN
649
650	message_string = 'PARAMETERIZED emissions mode selected:' // &
651	' EMISSIONS POSSIBLE ONLY ON STREET SURFACES' // &
652	' but values of scaling factors for street types' // &
653	' emiss_factor_main AND emiss_factor_side' // &
654	' not provided for each of the emissions species' // &
655	' or not provided at all: PLEASE set a finite value' // &
656	' for these parameters in the chemistry namelist'
657	CALL message( 'chem_emissions_matching', 'CM0442', 2, 2, 0, 6, 0 )
658	ENDIF
659	END DO
660
661
662	ELSE
663
664	message_string = 'Non of given Emission Species' // &
665	' matches' // &
666	' model chemical species' // &
667	' Emission routine is not called'
668	CALL message( 'chem_emissions_matching', 'CM0443', 0, 0, 0, 6, 0 )
669	ENDIF
670
671	ELSE
672
673	message_string = 'Array of Emission species not allocated: ' // &
674	' Either no emission species are provided as input or' // &
675	' no chemical species are used by PALM.' // &
676	' Emission routine is not called'
677	CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 )
678
679	ENDIF
680
681
682	!
683	!-- If emission module is switched on but mode_emis is not specified or it is given the wrong name
684	CASE DEFAULT
685
686	message_string = 'Emission Module switched ON, but' // &
687	' either no emission mode specified or incorrectly given :' // &
688	' please, pass the correct value to the namelist parameter "mode_emis"'
689	CALL message( 'chem_emissions_matching', 'CM0445', 2, 2, 0, 6, 0 )
690
691	END SELECT
692
693	END SUBROUTINE chem_emissions_match
694
695
696	!------------------------------------------------------------------------------!
697	! Description:
698	! ------------
699	!> Initialization:
700	!> Netcdf reading, arrays allocation and first calculation of cssws
701	!> fluxes at timestep 0
702	!------------------------------------------------------------------------------!
703	SUBROUTINE chem_emissions_init
704
705	USE surface_mod, &
706	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
707
708	USE netcdf_data_input_mod, &
709	ONLY: chem_emis, chem_emis_att
710
711	IMPLICIT NONE
712
713	CHARACTER (LEN=80) :: units !< units of inputs
714
715	INTEGER(iwp) :: ispec !< running index
716
717	!
718	!-- Actions for initial runs
719	! IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
720	!-- ...
721	!
722	!
723	!-- Actions for restart runs
724	! ELSE
725	!-- ...
726	!
727	! ENDIF
728
729
730	CALL location_message( 'Starting initialization of emission arrays', .FALSE. )
731
732	!
733	!-- Matching
734	CALL chem_emissions_match( chem_emis_att, nspec_out )
735
736	IF ( nspec_out == 0 ) THEN
737
738	emission_output_required = .FALSE.
739
740	ELSE
741
742	emission_output_required = .TRUE.
743
744
745	!
746	!-- Set molecule masses'
747	ALLOCATE( chem_emis_att%xm(nspec_out) )
748
749	DO ispec = 1, nspec_out
750	SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
751	CASE ( 'SO2' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 2 !< kg/mole
752	CASE ( 'SO4' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 4 !< kg/mole
753	CASE ( 'NO' ); chem_emis_att%xm(ispec) = xm_N + xm_O !< kg/mole
754	CASE ( 'NO2' ); chem_emis_att%xm(ispec) = xm_N + xm_O * 2 !< kg/mole
755	CASE ( 'NH3' ); chem_emis_att%xm(ispec) = xm_N + xm_H * 3 !< kg/mole
756	CASE ( 'CO' ); chem_emis_att%xm(ispec) = xm_C + xm_O !< kg/mole
757	CASE ( 'CO2' ); chem_emis_att%xm(ispec) = xm_C + xm_O * 2 !< kg/mole
758	CASE ( 'CH4' ); chem_emis_att%xm(ispec) = xm_C + xm_H * 4 !< kg/mole
759	CASE ( 'HNO3' ); chem_emis_att%xm(ispec) = xm_H + xm_N + xm_O*3 !< kg/mole
760	CASE DEFAULT
761	chem_emis_att%xm(ispec) = 1.0_wp
762	END SELECT
763	ENDDO
764
765
766	!
767	!-- assign emission values
768	SELECT CASE ( TRIM( mode_emis ) )
769
770
771	!
772	!-- PRE-PROCESSED case
773	CASE ( "PRE-PROCESSED" )
774
775	IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution(nzb:nzt+1,0:ny,0:nx,nspec_out) )
776
777	!
778	!-- Get emissions at the first time step
779	CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out )
780
781	!
782	!-- Default case
783	CASE ( "DEFAULT" )
784
785	IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution(1,0:ny,0:nx,nspec_out) )
786
787	!
788	!-- Get emissions at the first time step
789	CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out )
790
791	!
792	!-- PARAMETERIZED case
793	CASE ( "PARAMETERIZED" )
794
795	IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution(1,0:ny,0:nx,nspec_out) )
796
797	!
798	!-- Get emissions at the first time step
799	CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out)
800
801	END SELECT
802
803	ENDIF
804
805	END SUBROUTINE chem_emissions_init
806
807
808
809	!------------------------------------------------------------------------------!
810	! Description:
811	! ------------
812	!> Routine for Update of Emission values at each timestep
813	!-------------------------------------------------------------------------------!
814
815	SUBROUTINE chem_emissions_setup( emt_att, emt, nspec_out )
816
817	USE arrays_3d, &
818	ONLY: dzw
819	USE grid_variables, &
820	ONLY: dx, dy
821	USE indices, &
822	ONLY: nnx, nny, nnz
823	USE surface_mod, &
824	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
825	USE netcdf_data_input_mod, &
826	ONLY: street_type_f
827	USE arrays_3d, &
828	ONLY: hyp, pt
829
830
831	IMPLICIT NONE
832
833
834	TYPE(chem_emis_att_type), INTENT(INOUT) :: emt_att !< variable to store emission information
835
836	TYPE(chem_emis_val_type), INTENT(INOUT), ALLOCATABLE, DIMENSION(:) :: emt !< variable to store emission input values,
837	!< depending on the considered species
838
839	INTEGER,INTENT(IN) :: nspec_out !< Output of matching routine with number
840	!< of matched species
841
842	CHARACTER(LEN=80) :: units !< Units of the emission data
843
844	INTEGER(iwp) :: i !< running index for grid in x-direction
845	INTEGER(iwp) :: j !< running index for grid in y-direction
846	INTEGER(iwp) :: k !< running index for grid in z-direction
847	INTEGER(iwp) :: m !< running index for horizontal surfaces
848
849	INTEGER(iwp) :: icat !< Index for number of categories
850	INTEGER(iwp) :: ispec !< index for number of species
851	INTEGER(iwp) :: i_pm_comp !< index for number of PM components
852	INTEGER(iwp) :: ivoc !< Index for number of VOCs
853	INTEGER(iwp) :: time_index !< Index for time
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 |