Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |