Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |