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

Last change on this file since 3351 was 3337, checked in by kanani, 6 years ago

reintegrate branch resler to trunk

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