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

Last change on this file since 3572 was 3570, checked in by kanani, 5 years ago

Fix too long lines (chemistry_model_mod, chem_emissions_mod), correct terminal message (palmrun)

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