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 3483 2018-11-02 14:19:26Z eckhard $ |
---|
29 | ! bugfix: wrong locations of netCDF directives fixed |
---|
30 | ! |
---|
31 | ! 3467 2018-10-30 19:05:21Z suehring |
---|
32 | ! Enabled PARAMETRIZED mode for default surfaces when LSM is not applied but |
---|
33 | ! salsa is used |
---|
34 | ! |
---|
35 | ! 3458 2018-10-30 14:51:23Z kanani |
---|
36 | ! from chemistry branch r3443, banzhafs, Russo: |
---|
37 | ! Additional correction for index of input file of pre-processed mode |
---|
38 | ! Removed atomic and molecular weights as now available in chem_modules and |
---|
39 | ! added par_emis_time_factor (formerly in netcdf_data_input_mod) |
---|
40 | ! Added correction for index of input file of pre-processed mode |
---|
41 | ! Added a condition for default mode necessary for information read from ncdf file |
---|
42 | ! in pre-processed and default mode |
---|
43 | ! Correction of multiplying factor necessary for scaling emission values in time |
---|
44 | ! Introduced correction for scaling units in the case of DEFAULT emission mode |
---|
45 | ! |
---|
46 | ! 3373 2018-10-18 15:25:56Z kanani |
---|
47 | ! Fix wrong location of __netcdf directive |
---|
48 | ! |
---|
49 | ! 3337 2018-10-12 15:17:09Z kanani |
---|
50 | ! (from branch resler) |
---|
51 | ! Formatting |
---|
52 | ! |
---|
53 | ! 3312 2018-10-06 14:15:46Z knoop |
---|
54 | ! Initial revision |
---|
55 | ! |
---|
56 | ! 3286 2018-09-28 07:41:39Z forkel |
---|
57 | ! |
---|
58 | ! Authors: |
---|
59 | ! -------- |
---|
60 | ! @author Emmanuele Russo (Fu-Berlin) |
---|
61 | ! @author Sabine Banzhaf (FU-Berlin) |
---|
62 | ! @author Martijn Schaap (FU-Berlin, TNO Utrecht) |
---|
63 | ! |
---|
64 | ! Description: |
---|
65 | ! ------------ |
---|
66 | !> MODULE for reading-in Chemistry Emissions data |
---|
67 | !> |
---|
68 | !> @todo Check_parameters routine should be developed: for now it includes just one condition |
---|
69 | !> @todo Use of Restart files not contempled at the moment |
---|
70 | !> @todo revise indices of files read from the netcdf: preproc_emission_data and expert_emission_data |
---|
71 | !> @todo for now emission data may be passed on a singular vertical level: need to be more flexible |
---|
72 | !> @note <Enter notes on the module> |
---|
73 | !> @bug <Enter known bugs here> |
---|
74 | !------------------------------------------------------------------------------! |
---|
75 | |
---|
76 | MODULE chem_emissions_mod |
---|
77 | |
---|
78 | USE arrays_3d, & |
---|
79 | ONLY: rho_air |
---|
80 | |
---|
81 | USE control_parameters, & |
---|
82 | ONLY: initializing_actions, end_time, message_string, & |
---|
83 | intermediate_timestep_count, dt_3d |
---|
84 | |
---|
85 | USE indices |
---|
86 | |
---|
87 | USE kinds |
---|
88 | |
---|
89 | #if defined ( __netcdf ) |
---|
90 | USE NETCDF |
---|
91 | #endif |
---|
92 | |
---|
93 | USE netcdf_data_input_mod, & |
---|
94 | ONLY: chem_emis_att_type, chem_emis_val_type |
---|
95 | |
---|
96 | USE date_and_time_mod, & |
---|
97 | ONLY: time_default_indices, time_preprocessed_indices, & |
---|
98 | month_of_year, day_of_month, hour_of_day, & |
---|
99 | index_mm, index_dd, index_hh |
---|
100 | |
---|
101 | USE chem_gasphase_mod, & |
---|
102 | ONLY: spc_names |
---|
103 | |
---|
104 | USE chem_modules |
---|
105 | |
---|
106 | USE statistics, & |
---|
107 | ONLY: weight_pres |
---|
108 | |
---|
109 | IMPLICIT NONE |
---|
110 | |
---|
111 | !-- Declare all global variables within the module |
---|
112 | |
---|
113 | CHARACTER (LEN=80) :: filename_emis !> Variable for the name of the netcdf input file |
---|
114 | |
---|
115 | INTEGER(iwp) :: i !> index 1st selected dimension (some dims are not spatial) |
---|
116 | INTEGER(iwp) :: j !> index 2nd selected dimension |
---|
117 | INTEGER(iwp) :: i_start !> Index to start read variable from netcdf along one dims |
---|
118 | INTEGER(iwp) :: i_end !> Index to end read variable from netcdf in one dims |
---|
119 | INTEGER(iwp) :: j_start !> Index to start read variable from netcdf in additional dims |
---|
120 | INTEGER(iwp) :: j_end !> Index to end read variable from netcdf in additional dims |
---|
121 | INTEGER(iwp) :: z_start !> Index to start read variable from netcdf in additional dims |
---|
122 | INTEGER(iwp) :: z_end !> Index to end read variable from netcdf in additional dims |
---|
123 | INTEGER(iwp) :: dt_emis !> Time Step Emissions |
---|
124 | INTEGER(iwp) :: len_index !> length of index (used for several indices) |
---|
125 | INTEGER(iwp) :: len_index_voc !> length of voc index |
---|
126 | INTEGER(iwp) :: len_index_pm !> length of PMs index |
---|
127 | REAL(wp) :: con_factor !> Units Conversion Factor |
---|
128 | |
---|
129 | |
---|
130 | ! --------------------------------------------------------------- |
---|
131 | ! Set Values of molecules, mols, etc |
---|
132 | ! --------------------------------------------------------------- |
---|
133 | |
---|
134 | !> Avogadro number |
---|
135 | REAL, PARAMETER :: Avog = 6.02205e23 ! mlc/mol |
---|
136 | |
---|
137 | !> Dobson units: |
---|
138 | REAL, PARAMETER :: Dobs = 2.68668e16 ! (mlc/cm2) / DU |
---|
139 | |
---|
140 | !> sesalt composition: |
---|
141 | ! (Seinfeld and Pandis, "Atmospheric Chemistry and Physics", |
---|
142 | ! table 7.8 "Composition of Sea-Salt", p. 444) |
---|
143 | REAL, PARAMETER :: massfrac_Cl_in_seasalt = 0.5504 ! (kg Cl )/(kg seasalt) |
---|
144 | REAL, PARAMETER :: massfrac_Na_in_seasalt = 0.3061 ! (kg Na )/(kg seasalt) |
---|
145 | REAL, PARAMETER :: massfrac_SO4_in_seasalt = 0.0768 ! (kg SO4)/(kg seasalt) |
---|
146 | |
---|
147 | !> other numbers |
---|
148 | REAL, PARAMETER :: xm_seasalt = 74.947e-3 ! kg/mol : NaCl, SO4, .. |
---|
149 | REAL, PARAMETER :: rho_seasalt = 2.2e3 ! kg/m3 |
---|
150 | |
---|
151 | !> * amonium sulphate |
---|
152 | |
---|
153 | REAL, PARAMETER :: xm_NH4HSO4 = xm_NH4 + xm_H + xm_SO4 ! kg/mol |
---|
154 | REAL, PARAMETER :: rho_NH4HSO4a = 1.8e3 ! kg/m3 |
---|
155 | |
---|
156 | ! --------------------------------------------------------------- |
---|
157 | ! gas |
---|
158 | ! --------------------------------------------------------------- |
---|
159 | |
---|
160 | !> gas constant |
---|
161 | REAL, PARAMETER :: Rgas = 8.3144 ! J/mol/K |
---|
162 | |
---|
163 | !> gas constant for dry air |
---|
164 | REAL, PARAMETER :: Rgas_air = Rgas / xm_air ! J/kg/K |
---|
165 | |
---|
166 | !> water vapour |
---|
167 | REAL, PARAMETER :: Rgas_h2o = Rgas / xm_h2o ! J/kg/K |
---|
168 | |
---|
169 | !> standard pressure |
---|
170 | REAL, PARAMETER :: p0 = 1.0e5 ! Pa |
---|
171 | |
---|
172 | !> sea-level pressure |
---|
173 | REAL, PARAMETER :: p_sealevel = 1.01325e5 ! Pa <-- suggestion Bram Bregman |
---|
174 | |
---|
175 | !> global mean pressure |
---|
176 | REAL, PARAMETER :: p_global = 98500.0 ! Pa |
---|
177 | |
---|
178 | !> specific heat of dry air at constant pressure |
---|
179 | REAL, PARAMETER :: cp_air = 1004.0 ! J/kg/K |
---|
180 | |
---|
181 | !> Latent heat of evaporation |
---|
182 | REAL, PARAMETER :: lvap = 2.5e6 ! [J kg-1] |
---|
183 | |
---|
184 | !> Latent heat of condensation at 0 deg Celcius |
---|
185 | ! (heat (J) necesarry to evaporate 1 kg of water) |
---|
186 | REAL, PARAMETER :: Lc = 22.6e5 ! J/kg |
---|
187 | |
---|
188 | !> kappa = R/cp = 2/7 |
---|
189 | REAL, PARAMETER :: kappa = 2.0/7.0 |
---|
190 | |
---|
191 | !> Von Karman constant (dry_dep) |
---|
192 | REAL, PARAMETER :: vkarman = 0.4 |
---|
193 | |
---|
194 | !> Boltzmann constant: |
---|
195 | REAL(wp), PARAMETER :: kbolz = 1.38066e-23_wp ! J/K |
---|
196 | |
---|
197 | !> Inverse Reference Pressure (1/Pa) |
---|
198 | REAL(wp), PARAMETER :: pref_i = 1.0_wp / 100000.0_wp |
---|
199 | |
---|
200 | !> |
---|
201 | REAL(wp), PARAMETER :: r_cp = 0.286_wp !< R / cp (exponent for potential temperature) |
---|
202 | |
---|
203 | |
---|
204 | SAVE |
---|
205 | |
---|
206 | !-- Interfaces Part |
---|
207 | |
---|
208 | !-- Checks Input parameters |
---|
209 | INTERFACE chem_emissions_check_parameters |
---|
210 | MODULE PROCEDURE chem_emissions_check_parameters |
---|
211 | END INTERFACE chem_emissions_check_parameters |
---|
212 | |
---|
213 | !-- Reading of NAMELIST parameters |
---|
214 | ! INTERFACE chem_emissions_parin |
---|
215 | ! MODULE PROCEDURE chem_emissions_parin |
---|
216 | ! END INTERFACE chem_emissions_parin |
---|
217 | |
---|
218 | !-- Output of information to the header file |
---|
219 | ! INTERFACE chem_emissions_header |
---|
220 | ! MODULE PROCEDURE chem_emissions_header |
---|
221 | ! END INTERFACE chem_emissions_header |
---|
222 | |
---|
223 | !-- Matching Emissions actions |
---|
224 | INTERFACE chem_emissions_match |
---|
225 | MODULE PROCEDURE chem_emissions_match |
---|
226 | END INTERFACE chem_emissions_match |
---|
227 | |
---|
228 | !-- Initialization actions |
---|
229 | INTERFACE chem_emissions_init |
---|
230 | MODULE PROCEDURE chem_emissions_init |
---|
231 | END INTERFACE chem_emissions_init |
---|
232 | |
---|
233 | !-- Setup of Emissions |
---|
234 | INTERFACE chem_emissions_setup |
---|
235 | MODULE PROCEDURE chem_emissions_setup |
---|
236 | END INTERFACE chem_emissions_setup |
---|
237 | |
---|
238 | |
---|
239 | !-- Public Interfaces |
---|
240 | PUBLIC chem_emissions_init,chem_emissions_match,chem_emissions_setup |
---|
241 | |
---|
242 | !-- Public Variables |
---|
243 | |
---|
244 | PUBLIC con_factor, len_index,len_index_voc,len_index_pm |
---|
245 | |
---|
246 | CONTAINS |
---|
247 | |
---|
248 | !------------------------------------------------------------------------------! |
---|
249 | ! Description: |
---|
250 | ! ------------ |
---|
251 | !> Routine for checking input parameters |
---|
252 | !------------------------------------------------------------------------------! |
---|
253 | SUBROUTINE chem_emissions_check_parameters |
---|
254 | |
---|
255 | |
---|
256 | !TBD: Where Should we put the call to chem_emissions_check_parameters? In chem_init or in check_parameters? |
---|
257 | |
---|
258 | IMPLICIT NONE |
---|
259 | |
---|
260 | INTEGER(iwp) :: tmp |
---|
261 | |
---|
262 | TYPE(chem_emis_att_type) :: emt |
---|
263 | ! |
---|
264 | !-- Call routine for controlling chem_emissions namelist |
---|
265 | ! CALL chem_emissions_parin |
---|
266 | |
---|
267 | !TBD: In case the given emission values do not coincide with the passed names we should think of a solution: |
---|
268 | ! I would personally do that if the name of the species differ from the number of emission values, I would |
---|
269 | ! print a warning that says that in correspondance of that particular species the value is zero. |
---|
270 | ! An alternative would be to initialize the cs_surface_flux values to a negative number. |
---|
271 | |
---|
272 | !-- Check Emission Species Number equal to number of passed names for the chemistry species: |
---|
273 | IF ( SIZE(emt%species_name) /= emt%nspec ) THEN |
---|
274 | |
---|
275 | message_string = 'Numbers of input emission species names and number of species' // & |
---|
276 | 'for which emission values are given do not match' |
---|
277 | CALL message( 'chem_emissions_check_parameters', 'CM0437', 2, 2, 0, 6, 0 ) |
---|
278 | |
---|
279 | |
---|
280 | ENDIF |
---|
281 | |
---|
282 | !-- Check Emission Species Number equals to number of passed names for the species |
---|
283 | !IF ( SIZE(emt%species_name) /= SIZE(emt%species_index) ) THEN |
---|
284 | ! message_string = 'Number of input emission species names and ' // & |
---|
285 | ! ' number of passed species indices do not match' |
---|
286 | ! CALL message( 'chem_emissions_check_parameters', 'CM0101', 2, 2, 0, 6, 0 ) |
---|
287 | |
---|
288 | !ENDIF |
---|
289 | |
---|
290 | !-- Check Emission Categories |
---|
291 | ! IF ( SIZE(chem_emis%cat_name) /= SIZE(chem_emis%cat_index) ) THEN |
---|
292 | ! WRITE( message_string, * ) & |
---|
293 | ! 'Number of input emission categories name and categories index does not match\\' |
---|
294 | ! CALL message( 'chem_emissions_check_parameters', 'CM0101', 1, 2, 0, 6, 0 ) |
---|
295 | ! ENDIF |
---|
296 | |
---|
297 | |
---|
298 | |
---|
299 | !TBD: Check which other consistency controls should be included |
---|
300 | |
---|
301 | !TBD: Include check for spatial dimension of netcdf file. If they do not match with the ones |
---|
302 | ! of the simulation, the model should print an error. |
---|
303 | |
---|
304 | END SUBROUTINE chem_emissions_check_parameters |
---|
305 | |
---|
306 | !------------------------------------------------------------------------------! |
---|
307 | ! Description: |
---|
308 | ! ------------ |
---|
309 | !> Matching the chemical species indices. The routine checks what are the indices of the emission input species |
---|
310 | ! and the corresponding ones of the model species. The routine gives as output a vector containing the number |
---|
311 | ! of common species: it is important to note that while the model species are distinct, their values could be |
---|
312 | ! 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. |
---|
313 | !------------------------------------------------------------------------------! |
---|
314 | SUBROUTINE chem_emissions_match(emt_att,len_index) |
---|
315 | |
---|
316 | |
---|
317 | INTEGER(iwp), INTENT(INOUT) :: len_index !< Variable where to store the number of common species between the input dataset and the model species |
---|
318 | |
---|
319 | TYPE(chem_emis_att_type), INTENT(INOUT) :: emt_att !< Chemistry Emission Array containing information for all the input chemical emission species |
---|
320 | |
---|
321 | INTEGER(iwp) :: ind_mod, ind_inp !< Parameters for cycling through chemical model and input species |
---|
322 | INTEGER(iwp) :: nspec_emis_inp !< Variable where to store the number of the emission species in input |
---|
323 | |
---|
324 | INTEGER(iwp) :: ind_voc !< Indices to check whether a split for voc should be done |
---|
325 | |
---|
326 | INTEGER(iwp) :: ispec !> index for cycle over effective number of emission species |
---|
327 | |
---|
328 | |
---|
329 | !> Tell the user where we are!! |
---|
330 | CALL location_message( 'Matching input emissions and model chemistry species', .FALSE. ) |
---|
331 | |
---|
332 | !> Number of input emission species. |
---|
333 | |
---|
334 | nspec_emis_inp=emt_att%nspec |
---|
335 | |
---|
336 | !> Check the emission mode: DEFAULT, PRE-PROCESSED or PARAMETERIZED !TBD: Add option for capital or small letters |
---|
337 | SELECT CASE(TRIM(mode_emis)) |
---|
338 | |
---|
339 | !> 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. |
---|
340 | CASE ("PRE-PROCESSED") |
---|
341 | |
---|
342 | CALL location_message( 'chem_emissions PRE-PROCESSED_mode_matching', .FALSE. ) |
---|
343 | |
---|
344 | len_index=0 |
---|
345 | len_index_voc=0 |
---|
346 | |
---|
347 | !> The first condition is that both the number of model and input emissions species are not null |
---|
348 | IF ( (SIZE(spc_names) .GT. 0) .AND. (nspec_emis_inp .GT. 0)) THEN |
---|
349 | |
---|
350 | !> Cycle over model species |
---|
351 | DO ind_mod = 1, SIZE(spc_names) |
---|
352 | !> Cycle over Input species |
---|
353 | DO ind_inp = 1, nspec_emis_inp |
---|
354 | |
---|
355 | !> 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 |
---|
356 | ! > Check for VOC Species: In input in this case we only have one species (VOC) |
---|
357 | IF (TRIM(emt_att%species_name(ind_inp))=="VOC") THEN |
---|
358 | !> 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 |
---|
359 | DO ind_voc= 1,emt_att%nvoc |
---|
360 | !> Determine the indices of the coinciding species in the two cases and assign them to matching arrays |
---|
361 | IF (TRIM(emt_att%voc_name(ind_voc))==TRIM(spc_names(ind_mod))) THEN |
---|
362 | len_index=len_index+1 |
---|
363 | len_index_voc=len_index_voc+1 |
---|
364 | ENDIF |
---|
365 | END DO |
---|
366 | ENDIF |
---|
367 | !> Other Species |
---|
368 | IF (TRIM(emt_att%species_name(ind_inp))==TRIM(spc_names(ind_mod))) THEN |
---|
369 | len_index=len_index+1 |
---|
370 | ENDIF |
---|
371 | ENDDO |
---|
372 | ENDDO |
---|
373 | |
---|
374 | !> Allocate array for storing the indices of the matched species |
---|
375 | IF (len_index>0) THEN |
---|
376 | |
---|
377 | ALLOCATE(match_spec_input(len_index)) |
---|
378 | |
---|
379 | ALLOCATE(match_spec_model(len_index)) |
---|
380 | |
---|
381 | IF (len_index_voc>0) THEN |
---|
382 | |
---|
383 | ALLOCATE(match_spec_voc_model(len_index_voc)) !> Contains indices of the VOC model species |
---|
384 | |
---|
385 | 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 |
---|
386 | |
---|
387 | ENDIF |
---|
388 | |
---|
389 | !> Repeat the same cycle of above but passing the species indices to the newly declared arrays |
---|
390 | len_index=0 |
---|
391 | |
---|
392 | !> Cycle over model species |
---|
393 | DO ind_mod = 1, SIZE(spc_names) |
---|
394 | !> Cycle over Input species |
---|
395 | DO ind_inp = 1, nspec_emis_inp |
---|
396 | |
---|
397 | !> Determine the indices of the coinciding species in the two cases and assign them to matching arrays |
---|
398 | |
---|
399 | ! > VOCs |
---|
400 | IF ( TRIM(emt_att%species_name(ind_inp))=="VOC" .AND. ALLOCATED(match_spec_voc_input) ) THEN |
---|
401 | DO ind_voc= 1,emt_att%nvoc |
---|
402 | IF (TRIM(emt_att%voc_name(ind_voc))==TRIM(spc_names(ind_mod))) THEN |
---|
403 | len_index=len_index+1 |
---|
404 | len_index_voc=len_index_voc+1 |
---|
405 | |
---|
406 | match_spec_input(len_index)=ind_inp |
---|
407 | match_spec_model(len_index)=ind_mod |
---|
408 | |
---|
409 | match_spec_voc_input(len_index_voc)=ind_voc |
---|
410 | match_spec_voc_model(len_index_voc)=ind_mod |
---|
411 | ENDIF |
---|
412 | END DO |
---|
413 | ENDIF |
---|
414 | |
---|
415 | !>Other Species |
---|
416 | IF (TRIM(emt_att%species_name(ind_inp))==TRIM(spc_names(ind_mod))) THEN |
---|
417 | len_index=len_index+1 |
---|
418 | match_spec_input(len_index)=ind_inp |
---|
419 | match_spec_model(len_index)=ind_mod |
---|
420 | ENDIF |
---|
421 | END DO |
---|
422 | END DO |
---|
423 | |
---|
424 | !> In the case there are no species matching, the emission module should not be called |
---|
425 | ELSE |
---|
426 | |
---|
427 | message_string = 'Given Chemistry Emission Species' // & |
---|
428 | ' DO NOT MATCH' // & |
---|
429 | ' model chemical species:' // & |
---|
430 | ' Chemistry Emissions routine is not called' |
---|
431 | CALL message( 'chem_emissions_matching', 'CM0438', 0, 0, 0, 6, 0 ) |
---|
432 | ENDIF !> IF (len_index>0) |
---|
433 | |
---|
434 | ELSE |
---|
435 | |
---|
436 | !In this case either spc_names is null or nspec_emis_inp is not allocated |
---|
437 | |
---|
438 | message_string = 'Array of Emission species not allocated:' // & |
---|
439 | ' Either no emission species are provided as input or' // & |
---|
440 | ' no chemical species are used by PALM:' // & |
---|
441 | ' Chemistry Emissions routine is not called' |
---|
442 | CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 ) |
---|
443 | |
---|
444 | ENDIF !> IF ( (SIZE(spc_names) .GT. 0) .AND. ALLOCATED(nspec_emis_inp)) |
---|
445 | |
---|
446 | !> ------------------------------------------------------------------ |
---|
447 | !> DEFAULT case |
---|
448 | |
---|
449 | CASE ("DEFAULT") |
---|
450 | |
---|
451 | CALL location_message( 'chem_emissions DEFAULT_mode_matching', .FALSE. ) |
---|
452 | |
---|
453 | len_index=0 !> index for TOTAL number of species |
---|
454 | len_index_voc=0 !> index for TOTAL number of VOCs |
---|
455 | 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. |
---|
456 | |
---|
457 | !> The first condition is that both the number of model and input emissions species are not null |
---|
458 | IF ( (SIZE(spc_names) .GT. 0) .AND. (nspec_emis_inp .GT. 0) ) THEN |
---|
459 | |
---|
460 | !> Cycle over model species |
---|
461 | DO ind_mod = 1, SIZE(spc_names) |
---|
462 | !> Cycle over input species |
---|
463 | DO ind_inp = 1, nspec_emis_inp |
---|
464 | |
---|
465 | ! > Check for VOC Species: In input in this case we only have one species (VOC) |
---|
466 | IF (TRIM(emt_att%species_name(ind_inp))=="VOC") THEN |
---|
467 | !> 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 |
---|
468 | DO ind_voc= 1,emt_att%nvoc |
---|
469 | !> Determine the indices of the coinciding species in the two cases and assign them to matching arrays |
---|
470 | IF (TRIM(emt_att%voc_name(ind_voc))==TRIM(spc_names(ind_mod))) THEN |
---|
471 | len_index=len_index+1 |
---|
472 | len_index_voc=len_index_voc+1 |
---|
473 | ENDIF |
---|
474 | END DO |
---|
475 | ENDIF |
---|
476 | |
---|
477 | !> 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. |
---|
478 | IF (TRIM(emt_att%species_name(ind_inp))=="PM") THEN |
---|
479 | !>pm1 |
---|
480 | IF (TRIM(spc_names(ind_mod))=="PM1") THEN |
---|
481 | len_index=len_index+1 |
---|
482 | !>pm2.5 |
---|
483 | ELSE IF (TRIM(spc_names(ind_mod))=="PM25") THEN |
---|
484 | len_index=len_index+1 |
---|
485 | !>pm10 |
---|
486 | ELSE IF (TRIM(spc_names(ind_mod))=="PM10") THEN |
---|
487 | len_index=len_index+1 |
---|
488 | ENDIF |
---|
489 | ENDIF |
---|
490 | |
---|
491 | !> 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 |
---|
492 | IF (TRIM(emt_att%species_name(ind_inp))=="NOx") THEN |
---|
493 | !>no |
---|
494 | IF (TRIM(spc_names(ind_mod))=="NO") THEN |
---|
495 | len_index=len_index+1 |
---|
496 | !>no2 |
---|
497 | ELSE IF (TRIM(spc_names(ind_mod))=="NO2") THEN |
---|
498 | len_index=len_index+1 |
---|
499 | ENDIF |
---|
500 | ENDIF |
---|
501 | |
---|
502 | !>SOX: same as for NOx, but with SO2 and SO4 |
---|
503 | IF (TRIM(emt_att%species_name(ind_inp))=="SOx") THEN |
---|
504 | !>no |
---|
505 | IF (TRIM(spc_names(ind_mod))=="SO2") THEN |
---|
506 | len_index=len_index+1 |
---|
507 | !>no2 |
---|
508 | ELSE IF (TRIM(spc_names(ind_mod))=="SO4") THEN |
---|
509 | len_index=len_index+1 |
---|
510 | ENDIF |
---|
511 | ENDIF |
---|
512 | |
---|
513 | !> Other Species |
---|
514 | IF (TRIM(emt_att%species_name(ind_inp))==TRIM(spc_names(ind_mod))) THEN |
---|
515 | len_index=len_index+1 |
---|
516 | ENDIF |
---|
517 | END DO !> number of emission input species |
---|
518 | END DO !> number of model species |
---|
519 | |
---|
520 | |
---|
521 | !> Allocate Arrays |
---|
522 | IF (len_index>0) THEN |
---|
523 | |
---|
524 | ALLOCATE(match_spec_input(len_index)) |
---|
525 | ALLOCATE(match_spec_model(len_index)) |
---|
526 | |
---|
527 | IF (len_index_voc>0) THEN |
---|
528 | ALLOCATE(match_spec_voc_model(len_index_voc)) !> contains indices of the VOC model species |
---|
529 | ALLOCATE(match_spec_voc_input(len_index_voc)) !> In input there is only one value for VOCs in the DEFAULT mode. |
---|
530 | ! This array contains the indices of the different values of VOC compositions of the input variable VOC_composition |
---|
531 | ENDIF |
---|
532 | |
---|
533 | !> ASSIGN INDICES |
---|
534 | !> Repeat the same cycles of above, but passing the species indices to the newly declared arrays |
---|
535 | len_index=0 |
---|
536 | len_index_voc=0 |
---|
537 | |
---|
538 | DO ind_mod = 1, SIZE(spc_names) |
---|
539 | DO ind_inp = 1, nspec_emis_inp |
---|
540 | |
---|
541 | ! > VOCs |
---|
542 | IF ( TRIM(emt_att%species_name(ind_inp))=="VOC" .AND. ALLOCATED(match_spec_voc_input) ) THEN |
---|
543 | DO ind_voc= 1,emt_att%nvoc |
---|
544 | IF (TRIM(emt_att%voc_name(ind_voc))==TRIM(spc_names(ind_mod))) THEN |
---|
545 | len_index=len_index+1 |
---|
546 | len_index_voc=len_index_voc+1 |
---|
547 | |
---|
548 | match_spec_input(len_index)=ind_inp |
---|
549 | match_spec_model(len_index)=ind_mod |
---|
550 | |
---|
551 | match_spec_voc_input(len_index_voc)=ind_voc |
---|
552 | match_spec_voc_model(len_index_voc)=ind_mod |
---|
553 | ENDIF |
---|
554 | END DO |
---|
555 | ENDIF |
---|
556 | |
---|
557 | !> 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. |
---|
558 | IF (TRIM(emt_att%species_name(ind_inp))=="PM") THEN |
---|
559 | !>pm1 |
---|
560 | IF (TRIM(spc_names(ind_mod))=="PM1") THEN |
---|
561 | len_index=len_index+1 |
---|
562 | |
---|
563 | match_spec_input(len_index)=ind_inp |
---|
564 | match_spec_model(len_index)=ind_mod |
---|
565 | |
---|
566 | !match_spec_pm(1)=ind_mod |
---|
567 | !>pm2.5 |
---|
568 | ELSE IF (TRIM(spc_names(ind_mod))=="PM25") THEN |
---|
569 | len_index=len_index+1 |
---|
570 | |
---|
571 | match_spec_input(len_index)=ind_inp |
---|
572 | match_spec_model(len_index)=ind_mod |
---|
573 | |
---|
574 | !match_spec_pm(2)=ind_mod |
---|
575 | !>pm10 |
---|
576 | ELSE IF (TRIM(spc_names(ind_mod))=="PM10") THEN |
---|
577 | len_index=len_index+1 |
---|
578 | |
---|
579 | match_spec_input(len_index)=ind_inp |
---|
580 | match_spec_model(len_index)=ind_mod |
---|
581 | |
---|
582 | !match_spec_pm(3)=ind_mod |
---|
583 | ENDIF |
---|
584 | ENDIF |
---|
585 | |
---|
586 | !> NOx - The same as for PMs but here the model species are only 2: NO and NO2 |
---|
587 | IF (TRIM(emt_att%species_name(ind_inp))=="NOx") THEN |
---|
588 | !>no |
---|
589 | IF (TRIM(spc_names(ind_mod))=="NO") THEN |
---|
590 | len_index=len_index+1 |
---|
591 | |
---|
592 | match_spec_input(len_index)=ind_inp |
---|
593 | match_spec_model(len_index)=ind_mod |
---|
594 | |
---|
595 | !match_spec_nox(1)=ind_mod |
---|
596 | !>no2 |
---|
597 | ELSE IF (TRIM(spc_names(ind_mod))=="NO2") THEN |
---|
598 | len_index=len_index+1 |
---|
599 | |
---|
600 | match_spec_input(len_index)=ind_inp |
---|
601 | match_spec_model(len_index)=ind_mod |
---|
602 | |
---|
603 | ! match_spec_nox(2)=ind_mod |
---|
604 | ENDIF |
---|
605 | ENDIF |
---|
606 | |
---|
607 | !> SOx |
---|
608 | IF (TRIM(emt_att%species_name(ind_inp))=="SOx") THEN |
---|
609 | !>so2 |
---|
610 | IF (TRIM(spc_names(ind_mod))=="SO2") THEN |
---|
611 | len_index=len_index+1 |
---|
612 | |
---|
613 | match_spec_input(len_index)=ind_inp |
---|
614 | match_spec_model(len_index)=ind_mod |
---|
615 | |
---|
616 | ! match_spec_sox(1)=ind_mod |
---|
617 | !>so4 |
---|
618 | ELSE IF (TRIM(spc_names(ind_mod))=="SO4") THEN |
---|
619 | len_index=len_index+1 |
---|
620 | |
---|
621 | match_spec_input(len_index)=ind_inp |
---|
622 | match_spec_model(len_index)=ind_mod |
---|
623 | |
---|
624 | ! match_spec_sox(2)=ind_mod |
---|
625 | ENDIF |
---|
626 | ENDIF |
---|
627 | |
---|
628 | !> Other Species |
---|
629 | IF (TRIM(emt_att%species_name(ind_inp))==TRIM(spc_names(ind_mod))) THEN |
---|
630 | len_index=len_index+1 |
---|
631 | |
---|
632 | match_spec_input(len_index)=ind_inp |
---|
633 | match_spec_model(len_index)=ind_mod |
---|
634 | ENDIF |
---|
635 | END DO |
---|
636 | END DO |
---|
637 | |
---|
638 | ELSE |
---|
639 | |
---|
640 | message_string = 'Given Chemistry Emission Species' // & |
---|
641 | ' DO NOT MATCH' // & |
---|
642 | ' model chemical species' // & |
---|
643 | ' Chemistry Emissions routine is not called' |
---|
644 | CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 ) |
---|
645 | |
---|
646 | ENDIF |
---|
647 | |
---|
648 | ELSE |
---|
649 | |
---|
650 | message_string = 'Array of Emission species not allocated: ' // & |
---|
651 | ' Either no emission species are provided as input or' // & |
---|
652 | ' no chemical species are used by PALM:' // & |
---|
653 | ' Chemistry Emissions routine is not called' |
---|
654 | CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 ) |
---|
655 | |
---|
656 | ENDIF |
---|
657 | |
---|
658 | !-- 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. |
---|
659 | |
---|
660 | CASE ("PARAMETERIZED") |
---|
661 | |
---|
662 | len_index=0 |
---|
663 | |
---|
664 | !spc_names have to be greater than zero for proceeding |
---|
665 | IF ( (SIZE(spc_names) .GT. 0) .AND. (nspec_emis_inp .GT. 0) ) THEN |
---|
666 | |
---|
667 | |
---|
668 | !cycle over model species |
---|
669 | DO ind_mod = 1, SIZE(spc_names) |
---|
670 | ind_inp=1 |
---|
671 | DO WHILE ( TRIM( surface_csflux_name( ind_inp ) ) /= 'novalue' ) !<'novalue' is the default |
---|
672 | IF (TRIM(surface_csflux_name( ind_inp ))==TRIM(spc_names(ind_mod))) THEN |
---|
673 | len_index=len_index+1 |
---|
674 | ENDIF |
---|
675 | ind_inp=ind_inp+1 |
---|
676 | |
---|
677 | ENDDO |
---|
678 | ENDDO |
---|
679 | |
---|
680 | !Proceed only if there are matched species |
---|
681 | IF ( len_index .GT. 0 ) THEN |
---|
682 | |
---|
683 | !Allocation of Arrays of the matched species |
---|
684 | ALLOCATE(match_spec_input(len_index)) |
---|
685 | |
---|
686 | ALLOCATE(match_spec_model(len_index)) |
---|
687 | |
---|
688 | !Assign corresponding indices of input and model matched species |
---|
689 | len_index=0 |
---|
690 | |
---|
691 | DO ind_mod = 1, SIZE(spc_names) |
---|
692 | ind_inp = 1 |
---|
693 | DO WHILE ( TRIM( surface_csflux_name( ind_inp ) ) /= 'novalue' ) !<'novalue' is the default |
---|
694 | IF (TRIM( surface_csflux_name( ind_inp ) ) == TRIM(spc_names(ind_mod))) THEN |
---|
695 | len_index=len_index+1 |
---|
696 | match_spec_input(len_index)=ind_inp |
---|
697 | match_spec_model(len_index)=ind_mod |
---|
698 | ENDIF |
---|
699 | ind_inp=ind_inp+1 |
---|
700 | END DO |
---|
701 | END DO |
---|
702 | |
---|
703 | ! also, Add AN if to check that when the surface_csflux are passed to the namelist, also the street factor values are provided |
---|
704 | DO ispec = 1 , len_index |
---|
705 | |
---|
706 | IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND. emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN |
---|
707 | |
---|
708 | message_string = 'PARAMETERIZED emissions mode selected:' // & |
---|
709 | ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES' // & |
---|
710 | ' but values of scaling factors for street types' // & |
---|
711 | ' emiss_factor_main AND emiss_factor_side' // & |
---|
712 | ' not provided for each of the emissions species' // & |
---|
713 | ' or not provided at all: PLEASE set a finite value' // & |
---|
714 | ' for these parameters in the chemistry namelist' |
---|
715 | CALL message( 'chem_emissions_matching', 'CM0442', 2, 2, 0, 6, 0 ) |
---|
716 | ENDIF |
---|
717 | END DO |
---|
718 | |
---|
719 | |
---|
720 | ELSE |
---|
721 | |
---|
722 | message_string = 'Given Chemistry Emission Species' // & |
---|
723 | ' DO NOT MATCH' // & |
---|
724 | ' model chemical species' // & |
---|
725 | ' Chemistry Emissions routine is not called' |
---|
726 | CALL message( 'chem_emissions_matching', 'CM0443', 0, 0, 0, 6, 0 ) |
---|
727 | ENDIF |
---|
728 | |
---|
729 | ELSE |
---|
730 | |
---|
731 | message_string = 'Array of Emission species not allocated: ' // & |
---|
732 | ' Either no emission species are provided as input or' // & |
---|
733 | ' no chemical species are used by PALM.' // & |
---|
734 | ' Chemistry Emissions routine is not called' |
---|
735 | CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 ) |
---|
736 | |
---|
737 | ENDIF |
---|
738 | |
---|
739 | |
---|
740 | !-- CASE when emission module is switched on but mode_emis is not specified or it is given the wrong name |
---|
741 | CASE DEFAULT |
---|
742 | |
---|
743 | message_string = 'Chemistry Emissions Module switched ON, but' // & |
---|
744 | ' either no emission mode specified or incorrectly given :' // & |
---|
745 | ' please, pass the correct value to the namelist parameter "mode_emis"' |
---|
746 | CALL message( 'chem_emissions_matching', 'CM0445', 2, 2, 0, 6, 0 ) |
---|
747 | |
---|
748 | END SELECT |
---|
749 | |
---|
750 | END SUBROUTINE chem_emissions_match |
---|
751 | |
---|
752 | !------------------------------------------------------------------------------! |
---|
753 | ! Description: |
---|
754 | ! ------------ |
---|
755 | !> Initialization: |
---|
756 | !> Netcdf reading, arrays allocation and first calculation of cssws fluxes at timestep 0 |
---|
757 | !> |
---|
758 | !------------------------------------------------------------------------------! |
---|
759 | SUBROUTINE chem_emissions_init(emt_att,emt,nspec_out) |
---|
760 | |
---|
761 | USE surface_mod, & |
---|
762 | ONLY: surf_lsm_h,surf_def_h,surf_usm_h |
---|
763 | |
---|
764 | IMPLICIT NONE |
---|
765 | |
---|
766 | TYPE(chem_emis_att_type),INTENT(INOUT) :: emt_att !> variable where to store all the information of |
---|
767 | ! emission inputs whose values do not depend |
---|
768 | ! on the considered species |
---|
769 | |
---|
770 | TYPE(chem_emis_val_type),INTENT(INOUT), ALLOCATABLE, DIMENSION(:) :: emt !> variable where to store emission inputs values, |
---|
771 | ! depending on the considered species |
---|
772 | |
---|
773 | INTEGER(iwp),INTENT(INOUT) :: nspec_out !> number of outputs of chemistry emission routines |
---|
774 | |
---|
775 | CHARACTER (LEN=80) :: units !> units of chemistry inputs |
---|
776 | |
---|
777 | INTEGER(iwp) :: ispec !> Index to go through the emission chemical species |
---|
778 | |
---|
779 | |
---|
780 | !-- Actions for initial runs : TBD: needs to be updated |
---|
781 | IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN |
---|
782 | !-- ... |
---|
783 | |
---|
784 | ! |
---|
785 | !-- Actions for restart runs |
---|
786 | ELSE |
---|
787 | !-- ... |
---|
788 | |
---|
789 | ENDIF |
---|
790 | |
---|
791 | |
---|
792 | CALL location_message( 'Starting initialization of chemistry emissions arrays', .FALSE. ) |
---|
793 | |
---|
794 | |
---|
795 | !-- Match Input and Model emissions |
---|
796 | CALL chem_emissions_match(emt_att,nspec_out) |
---|
797 | |
---|
798 | !> If nspec_out==0, then do not use emission module: The corresponding message is already printed in the matching routine |
---|
799 | IF ( nspec_out == 0 ) THEN |
---|
800 | |
---|
801 | emission_output_required=.FALSE. |
---|
802 | |
---|
803 | ELSE |
---|
804 | |
---|
805 | emission_output_required=.TRUE. |
---|
806 | |
---|
807 | |
---|
808 | !----------------------------------------------------------------- |
---|
809 | !> Set molecule masses' |
---|
810 | ALLOCATE(emt_att%xm(nspec_out)) |
---|
811 | ! loop over emisisons: |
---|
812 | |
---|
813 | DO ispec = 1, nspec_out |
---|
814 | ! switch: |
---|
815 | SELECT CASE ( TRIM(spc_names(match_spec_model(ispec))) ) |
---|
816 | CASE ( 'SO2','SO4' ) ; emt_att%xm(ispec) = xm_S + xm_O * 2 ! kg/mole |
---|
817 | CASE ( 'NO','NO2' ) ; emt_att%xm(ispec) = xm_N + xm_O * 2 ! kg/mole NO2 equivalent |
---|
818 | CASE ( 'NH3' ) ; emt_att%xm(ispec) = xm_N + xm_H * 3 ! kg/mole |
---|
819 | CASE ( 'CO' ) ; emt_att%xm(ispec) = xm_C + xm_O ! kg/mole |
---|
820 | CASE ( 'CO2' ) ; emt_att%xm(ispec) = xm_C + xm_O * 2 ! kg/mole |
---|
821 | CASE ( 'CH4' ) ; emt_att%xm(ispec) = xm_C + xm_H * 4 ! kg/mole |
---|
822 | CASE ( 'HNO3' ); emt_att%xm(ispec) = xm_H + xm_N + xm_O*3 !kg/mole |
---|
823 | CASE DEFAULT |
---|
824 | !-- TBD: check if this hase to be removed |
---|
825 | !emt_att%xm(ispec) = 1.0_wp |
---|
826 | END SELECT |
---|
827 | ENDDO |
---|
828 | |
---|
829 | |
---|
830 | !> ASSIGN EMISSION VALUES for the different emission modes |
---|
831 | SELECT CASE(TRIM(mode_emis)) !TBD: Add the option for CApital or small letters |
---|
832 | |
---|
833 | |
---|
834 | !> PRE-PROCESSED case |
---|
835 | CASE ("PRE-PROCESSED") |
---|
836 | |
---|
837 | IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(nzb:nzt+1,0:ny,0:nx,nspec_out)) |
---|
838 | |
---|
839 | CALL location_message( 'emis_distribution array allocated in PRE-PROCESSED mode', .FALSE. ) |
---|
840 | |
---|
841 | !> Calculate the values of the emissions at the first time step |
---|
842 | CALL chem_emissions_setup(emt_att,emt,nspec_out) |
---|
843 | |
---|
844 | !> Default case |
---|
845 | CASE ("DEFAULT") |
---|
846 | |
---|
847 | IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution( 1, 0:ny, 0:nx, nspec_out ) ) |
---|
848 | |
---|
849 | CALL location_message( 'emis_distribution array allocated in DEFAULT 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 | !> PARAMETERIZED case |
---|
855 | CASE ("PARAMETERIZED") |
---|
856 | |
---|
857 | CALL location_message( 'emis_distribution array allocated in PARAMETERIZED mode', .FALSE. ) |
---|
858 | |
---|
859 | ! For now for PAR and DEF values only, first vertical level of emis_distribution is allocated, while for PRE-PROCESSED all. |
---|
860 | IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(1,0:ny,0:nx,nspec_out)) |
---|
861 | |
---|
862 | !> Calculate the values of the emissions at the first time step |
---|
863 | CALL chem_emissions_setup(emt_att,emt,nspec_out) |
---|
864 | |
---|
865 | END SELECT |
---|
866 | |
---|
867 | ENDIF |
---|
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 salsa_mod, & |
---|
888 | ONLY: salsa |
---|
889 | USE surface_mod, & |
---|
890 | ONLY: surf_lsm_h,surf_def_h,surf_usm_h |
---|
891 | USE netcdf_data_input_mod, & |
---|
892 | ONLY: street_type_f |
---|
893 | USE arrays_3d, & |
---|
894 | ONLY: hyp, pt |
---|
895 | |
---|
896 | IMPLICIT NONE |
---|
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 ppm*dz basically we need to multiply the input by (m**2/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**2*s). Units in this case are always micromoles/m**2*day (or kilograms/m**2*day 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)*(24*3600)=value/day ---- (value/day)*time_factor=value/hour ---(value/hour)/(3600)=value/sec ) |
---|
1247 | ! ((value/sec)*(24*3600)*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/m2*s kg/m2*s |
---|
1252 | delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,1)*con_factor*hour_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_factor*hour_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/m2*s kg/m2*s |
---|
1266 | delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,1)*con_factor*hour_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/m2*s kg/m2*s |
---|
1272 | delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,2)*con_factor*hour_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_factor*hour_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_factor*hour_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_factor*hour_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**2*s but the composition is provided in mole/kg: a conversion factor for the input that could be eventually provided in, for example, tons/(m**2*s) 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_factor*hour_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**2*s). 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^2*s) *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^2*s) *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 | ELSEIF ( salsa ) THEN |
---|
1449 | DO m = 1, surf_def_h(0)%ns |
---|
1450 | i = surf_def_h(0)%i(m) |
---|
1451 | j = surf_def_h(0)%j(m) |
---|
1452 | k = surf_def_h(0)%k(m) |
---|
1453 | |
---|
1454 | |
---|
1455 | IF ( street_type_f%var(j,i) >= main_street_id .AND. & |
---|
1456 | street_type_f%var(j,i) < max_street_id ) & |
---|
1457 | THEN |
---|
1458 | |
---|
1459 | !> Cycle over already matched species |
---|
1460 | DO ispec=1,nspec_out |
---|
1461 | |
---|
1462 | !> PMs are already in mass units:micrograms: have to be converted to kilograms |
---|
1463 | IF ( TRIM(spc_names(match_spec_model(ispec)))=="PM1" & |
---|
1464 | .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & |
---|
1465 | .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10")& |
---|
1466 | THEN |
---|
1467 | |
---|
1468 | surf_def_h(0)%cssws(match_spec_model(ispec),m) = & |
---|
1469 | emiss_factor_main(match_spec_input(ispec)) * & |
---|
1470 | emis_distribution(1,j,i,ispec) * rho_air(k) /& |
---|
1471 | time_factor(1) |
---|
1472 | ELSE |
---|
1473 | |
---|
1474 | !> Other Species: inputs are micromoles: have to be converted |
---|
1475 | surf_def_h(0)%cssws(match_spec_model(ispec),m) = & |
---|
1476 | emiss_factor_main(match_spec_input(ispec)) * & |
---|
1477 | emis_distribution(1,j,i,ispec) * & |
---|
1478 | conv_to_ratio(k,j,i) * rho_air(k) / time_factor(1) |
---|
1479 | ENDIF |
---|
1480 | ENDDO |
---|
1481 | |
---|
1482 | ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. & |
---|
1483 | street_type_f%var(j,i) < main_street_id ) & |
---|
1484 | THEN |
---|
1485 | |
---|
1486 | !> Cycle over already matched species |
---|
1487 | DO ispec=1,nspec_out |
---|
1488 | |
---|
1489 | !> PMs are already in mass units: micrograms |
---|
1490 | IF ( TRIM(spc_names(match_spec_model(ispec)))=="PM1" & |
---|
1491 | .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & |
---|
1492 | .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10")& |
---|
1493 | THEN |
---|
1494 | |
---|
1495 | surf_def_h(0)%cssws(match_spec_model(ispec),m) = & |
---|
1496 | emiss_factor_side(match_spec_input(ispec)) * & |
---|
1497 | emis_distribution(1,j,i,ispec) * rho_air(k) / & |
---|
1498 | time_factor(1) |
---|
1499 | ELSE |
---|
1500 | |
---|
1501 | surf_def_h(0)%cssws(match_spec_model(ispec),m) = & |
---|
1502 | emiss_factor_side(match_spec_input(ispec)) * & |
---|
1503 | emis_distribution(1,j,i,ispec) * & |
---|
1504 | conv_to_ratio(k,j,i) * rho_air(k) / time_factor(1) |
---|
1505 | ENDIF |
---|
1506 | |
---|
1507 | ENDDO |
---|
1508 | |
---|
1509 | ELSE |
---|
1510 | |
---|
1511 | !> If no street type is defined, then assign null emissions to all the species |
---|
1512 | surf_def_h(0)%cssws(:,m) = 0.0_wp |
---|
1513 | |
---|
1514 | ENDIF |
---|
1515 | |
---|
1516 | ENDDO |
---|
1517 | |
---|
1518 | ENDIF |
---|
1519 | |
---|
1520 | ENDIF |
---|
1521 | |
---|
1522 | !> For both DEFAULT and PRE-PROCESSED |
---|
1523 | ELSE |
---|
1524 | |
---|
1525 | |
---|
1526 | DO ispec=1,nspec_out |
---|
1527 | !TBD: for the PRE-PROCESSED mode in the future, values at different heights should be included! |
---|
1528 | !> Default surface type |
---|
1529 | IF (surf_def_h(0)%ns .GT. 0) THEN |
---|
1530 | |
---|
1531 | DO m = 1, surf_def_h(0)%ns |
---|
1532 | |
---|
1533 | i = surf_def_h(0)%i(m) |
---|
1534 | j = surf_def_h(0)%j(m) |
---|
1535 | |
---|
1536 | IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN |
---|
1537 | |
---|
1538 | |
---|
1539 | !> Distinguish between PMs (no needing conversion in ppms), |
---|
1540 | ! VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and |
---|
1541 | ! other Species (still expressed in Kg/(m**2*s) at this point) |
---|
1542 | |
---|
1543 | !> PMs |
---|
1544 | IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & |
---|
1545 | .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN |
---|
1546 | |
---|
1547 | ! kg/(m^2*s) *kg/m^3 kg/(m^2*s) |
---|
1548 | surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)* & |
---|
1549 | ! kg/m^3 |
---|
1550 | rho_air(nzb) |
---|
1551 | |
---|
1552 | |
---|
1553 | ELSE |
---|
1554 | |
---|
1555 | !> VOCs |
---|
1556 | IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN |
---|
1557 | ! ( ppm/s) * m * kg/m^3 mole/(m^2/s) |
---|
1558 | surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & |
---|
1559 | ! m^3/mole ppm |
---|
1560 | conv_to_ratio(nzb,j,i)*ratio2ppm * & |
---|
1561 | ! kg/m^3 |
---|
1562 | rho_air(nzb) |
---|
1563 | |
---|
1564 | |
---|
1565 | !> OTHER SPECIES |
---|
1566 | ELSE |
---|
1567 | |
---|
1568 | ! ( ppm/s )*m * kg/m^3 kg/(m^2/s) |
---|
1569 | surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & |
---|
1570 | ! mole/Kg |
---|
1571 | (1/emt_att%xm(ispec))* & |
---|
1572 | ! m^3/mole ppm |
---|
1573 | conv_to_ratio(nzb,j,i)*ratio2ppm* & |
---|
1574 | ! kg/m^3 |
---|
1575 | rho_air(nzb) |
---|
1576 | |
---|
1577 | |
---|
1578 | ENDIF |
---|
1579 | |
---|
1580 | ENDIF |
---|
1581 | |
---|
1582 | ENDIF |
---|
1583 | |
---|
1584 | ENDDO |
---|
1585 | |
---|
1586 | END IF |
---|
1587 | |
---|
1588 | !-- Land Surface Mode surface type |
---|
1589 | IF (surf_lsm_h%ns .GT. 0) THEN |
---|
1590 | |
---|
1591 | DO m = 1, surf_lsm_h%ns |
---|
1592 | |
---|
1593 | i = surf_lsm_h%i(m) |
---|
1594 | j = surf_lsm_h%j(m) |
---|
1595 | k = surf_lsm_h%k(m) |
---|
1596 | |
---|
1597 | IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN |
---|
1598 | |
---|
1599 | !> Distinguish between PMs (no needing conversion in ppms), |
---|
1600 | ! VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and |
---|
1601 | ! other Species (still expressed in Kg/(m**2*s) at this point) |
---|
1602 | |
---|
1603 | !> PMs |
---|
1604 | IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & |
---|
1605 | .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN |
---|
1606 | |
---|
1607 | ! kg/(m^2*s) * kg/m^3 kg/(m^2*s) |
---|
1608 | surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & |
---|
1609 | ! kg/m^3 |
---|
1610 | rho_air(k) |
---|
1611 | |
---|
1612 | ELSE |
---|
1613 | |
---|
1614 | !> VOCs |
---|
1615 | IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN |
---|
1616 | ! ( ppm/s) * m * kg/m^3 mole/(m^2/s) |
---|
1617 | surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & |
---|
1618 | ! m^3/mole ppm |
---|
1619 | conv_to_ratio(k,j,i)*ratio2ppm* & |
---|
1620 | ! kg/m^3 |
---|
1621 | rho_air(k) |
---|
1622 | |
---|
1623 | !> OTHER SPECIES |
---|
1624 | ELSE |
---|
1625 | |
---|
1626 | ! ( ppm/s) * m * kg/m^3 kg/(m^2/s) |
---|
1627 | surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & |
---|
1628 | ! mole/Kg |
---|
1629 | (1/emt_att%xm(ispec))* & |
---|
1630 | ! m^3/mole ppm |
---|
1631 | conv_to_ratio(k,j,i)*ratio2ppm* & |
---|
1632 | ! kg/m^3 |
---|
1633 | rho_air(k) |
---|
1634 | |
---|
1635 | ENDIF |
---|
1636 | |
---|
1637 | ENDIF |
---|
1638 | |
---|
1639 | ENDIF |
---|
1640 | |
---|
1641 | ENDDO |
---|
1642 | |
---|
1643 | END IF |
---|
1644 | |
---|
1645 | !-- Urban Surface Mode surface type |
---|
1646 | IF (surf_usm_h%ns .GT. 0) THEN |
---|
1647 | |
---|
1648 | |
---|
1649 | DO m = 1, surf_usm_h%ns |
---|
1650 | |
---|
1651 | i = surf_usm_h%i(m) |
---|
1652 | j = surf_usm_h%j(m) |
---|
1653 | k = surf_usm_h%k(m) |
---|
1654 | |
---|
1655 | IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN |
---|
1656 | |
---|
1657 | !> PMs |
---|
1658 | IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & |
---|
1659 | .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN |
---|
1660 | |
---|
1661 | ! kg/(m^2*s) *kg/m^3 kg/(m^2*s) |
---|
1662 | surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)* & |
---|
1663 | ! kg/m^3 |
---|
1664 | rho_air(k) |
---|
1665 | |
---|
1666 | |
---|
1667 | ELSE |
---|
1668 | |
---|
1669 | !> VOCs |
---|
1670 | IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN |
---|
1671 | ! ( ppm/s) * m * kg/m^3 mole/(m^2/s) |
---|
1672 | surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & |
---|
1673 | ! m^3/mole ppm |
---|
1674 | conv_to_ratio(k,j,i)*ratio2ppm* & |
---|
1675 | ! kg/m^3 |
---|
1676 | rho_air(k) |
---|
1677 | |
---|
1678 | !> OTHER SPECIES |
---|
1679 | ELSE |
---|
1680 | |
---|
1681 | |
---|
1682 | ! ( ppm/s) * m * kg/m^3 kg/(m^2/s) |
---|
1683 | surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & |
---|
1684 | ! mole/Kg |
---|
1685 | (1/emt_att%xm(ispec))* & |
---|
1686 | ! m^3/mole ppm |
---|
1687 | conv_to_ratio(k,j,i)*ratio2ppm* & |
---|
1688 | ! kg/m^3 |
---|
1689 | rho_air(k) |
---|
1690 | |
---|
1691 | |
---|
1692 | ENDIF |
---|
1693 | |
---|
1694 | ENDIF |
---|
1695 | |
---|
1696 | ENDIF |
---|
1697 | |
---|
1698 | ENDDO |
---|
1699 | END IF |
---|
1700 | ENDDO |
---|
1701 | |
---|
1702 | ENDIF |
---|
1703 | |
---|
1704 | |
---|
1705 | !> At the end of each call to chem_emissions setup, the time_factor is deallocated (ALLOCATED ONLY in the DEFAULT mode) |
---|
1706 | |
---|
1707 | IF ( ALLOCATED ( time_factor ) ) DEALLOCATE( time_factor ) |
---|
1708 | |
---|
1709 | ENDIF !> emis_output_required |
---|
1710 | |
---|
1711 | |
---|
1712 | END SUBROUTINE chem_emissions_setup |
---|
1713 | |
---|
1714 | END MODULE chem_emissions_mod |
---|