Home

Context Navigation

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

Last change on this file since 3311 was 3298, checked in by kanani, 6 years ago
Merge chemistry branch at r3297 to trunk
Property svn:keywords set to `Id`
File size: 76.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |