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

Last change on this file since 3312 was 3312, checked in by knoop, 5 years ago

Added proper exit code to PALM and fixed exit code handling by palmbuild and palmrun

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