Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |