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

Last change on this file since 3473 was 3467, checked in by suehring, 6 years ago

Branch salsa @3446 re-integrated into trunk

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