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

Last change on this file since 3483 was 3483, checked in by raasch, 3 years ago

bugfix: misplaced positions of cpp-directives for netCDF and MPI fixed; output format limited to a maximum line length of 80

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