Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |