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

Last change on this file since 3298 was 3298, checked in by kanani, 3 years ago

Merge chemistry branch at r3297 to trunk

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