Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |