Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |