Home

Context Navigation

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

Last change on this file since 3593 was 3589, checked in by suehring, 6 years ago
Remove erroneous UTF encoding; last commit documented
Property svn:keywords set to `Id`
File size: 77.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |