Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |