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

Last change on this file since 3589 was 3589, checked in by suehring, 5 years ago

Remove erroneous UTF encoding; last commit documented

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