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

Last change on this file since 3891 was 3885, checked in by kanani, 6 years ago

restructure/add location/debug messages

  • Property svn:keywords set to Id
File size: 64.2 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-2019 Leibniz Universitaet Hannover
18! Copyright 2018-2019 Freie Universitaet Berlin
19! Copyright 2018-2019 Karlsruhe Institute of Technology
20!--------------------------------------------------------------------------------!
21!
22! Current revisions:
23! ------------------
24!
25!
26! Former revisions:
27! -----------------
28! $Id: chem_emissions_mod.f90 3885 2019-04-11 11:29:34Z suehring $
29! Changes related to global restructuring of location messages and introduction
30! of additional debug messages
31!
32! 3831 2019-03-28 09:11:22Z forkel
33! added nvar to USE chem_gasphase_mod (chem_modules will not include nvar anymore)
34!
35! 3788 2019-03-07 11:40:09Z banzhafs
36! Removed unused variables from chem_emissions_mod
37!
38!3772 2019-02-28 15:51:57Z suehring
39! - In case of parametrized emissions, assure that emissions are only on natural
40!   surfaces (i.e. streets) and not on urban surfaces.
41! - some unnecessary if clauses removed
42!
43!3685 2019 -01-21 01:02:11Z knoop
44! Some interface calls moved to module_interface + cleanup
45!
46! 3611 2018-12-07 14:14:11Z banzhafs
47! Code update to comply PALM coding rules
48! removed unnecessary informative messsages/location
49! messages and corrected comments on PM units from to kg 
50! bug fix: spcs_name replaced by nvar in DO loops
51!
52! 3591 2018-11-30 17:37:32Z suehring
53! - Removed salsa dependency.
54! - Enabled PARAMETRIZED mode for default surfaces when LSM is not applied but
55!   salsa is (M. Kurppa)
56!
57! 3582 2018-11-29 19:16:36Z suehring
58! resler:
59! Break lines at 132 characters
60!
61! 3483 2018-11-02 14:19:26Z raasch
62! bugfix: wrong locations of netCDF directives fixed
63!
64! 3467 2018-10-30 19:05:21Z suehring
65! Enabled PARAMETRIZED mode for default surfaces when LSM is not applied but
66! salsa is used
67!
68! 3458 2018-10-30 14:51:23Z kanani
69! from chemistry branch r3443, banzhafs, Russo:
70! Additional correction for index of input file of pre-processed mode
71! Removed atomic and molecular weights as now available in chem_modules and
72! added par_emis_time_factor (formerly in netcdf_data_input_mod)
73! Added correction for index of input file of pre-processed mode
74! Added a condition for default mode necessary for information read from ncdf file
75! in pre-processed and default mode
76! Correction of multiplying factor necessary for scaling emission values in time
77! Introduced correction for scaling units in the case of DEFAULT emission mode
78!
79! 3373 2018-10-18 15:25:56Z kanani
80! Fix wrong location of __netcdf directive
81!
82! 3337 2018-10-12 15:17:09Z kanani
83! (from branch resler)
84! Formatting
85!
86! 3312 2018-10-06 14:15:46Z knoop
87! Initial revision
88!
89! 3286 2018-09-28 07:41:39Z forkel
90!
91! Authors:
92! --------
93! @author Emmanuele Russo (FU-Berlin)
94! @author Sabine Banzhaf  (FU-Berlin)
95! @author Martijn Schaap  (FU-Berlin, TNO Utrecht)
96!
97! Description:
98! ------------
99!>  MODULE for reading-in Chemistry Emissions data
100!>
101!> @todo Rename nspec to n_emis to avoid inteferece with nspec from chem_gasphase_mod
102!> @todo Check_parameters routine should be developed: for now it includes just one condition
103!> @todo Use of Restart files not contempled at the moment
104!> @todo revise indices of files read from the netcdf: preproc_emission_data and expert_emission_data
105!> @todo for now emission data may be passed on a singular vertical level: need to be more flexible
106!> @todo fill/activate restart option in chem_emissions_init
107!> @todo discuss dt_emis
108!> @note <Enter notes on the module>
109!> @bug  <Enter known bugs here>
110!------------------------------------------------------------------------------!
111
112 MODULE chem_emissions_mod
113
114    USE arrays_3d,                                                             &
115        ONLY:  rho_air
116
117    USE control_parameters,                                                    &
118        ONLY:  debug_output,                                                   &
119               end_time, message_string, initializing_actions,                 &
120               intermediate_timestep_count, dt_3d
121 
122    USE indices
123
124    USE kinds
125
126#if defined( __netcdf )
127    USE netcdf
128#endif
129
130    USE netcdf_data_input_mod,                                                  &
131        ONLY: chem_emis_att_type, chem_emis_val_type
132
133    USE date_and_time_mod,                                                      &
134        ONLY: day_of_month, hour_of_day,                                        &
135             index_mm, index_dd, index_hh,                                      &
136             month_of_year, hour_of_day,                                        &
137             time_default_indices, time_preprocessed_indices
138   
139    USE chem_gasphase_mod,                                                      &
140        ONLY: nvar, spc_names
141 
142    USE chem_modules
143
144    USE statistics,                                                             &   
145        ONLY:  weight_pres
146
147   
148    IMPLICIT NONE
149
150!
151!-- Declare all global variables within the module
152   
153    CHARACTER (LEN=80) ::  filename_emis             !< Variable for the name of the netcdf input file
154
155    INTEGER(iwp) ::  i                               !< index 1st selected dimension (some dims are not spatial)
156    INTEGER(iwp) ::  j                               !< index 2nd selected dimension
157    INTEGER(iwp) ::  i_start                         !< Index to start read variable from netcdf along one dims
158    INTEGER(iwp) ::  i_end                           !< Index to end read variable from netcdf in one dims
159    INTEGER(iwp) ::  j_start                         !< Index to start read variable from netcdf in additional dims
160    INTEGER(iwp) ::  j_end                           !< Index to end read variable from netcdf in additional dims
161    INTEGER(iwp) ::  z_start                         !< Index to start read variable from netcdf in additional dims
162    INTEGER(iwp) ::  z_end                           !< Index to end read variable from netcdf in additional dims
163    INTEGER(iwp) ::  dt_emis                         !< Time Step Emissions
164    INTEGER(iwp) ::  len_index                       !< length of index (used for several indices)
165    INTEGER(iwp) ::  len_index_voc                   !< length of voc index
166    INTEGER(iwp) ::  len_index_pm                    !< length of PMs index
167
168    REAL(wp) ::  con_factor                          !< Units Conversion Factor
169                           
170    REAL(wp), PARAMETER ::  Rgas = 8.3144                    !< gas constant in J/mol/K           
171    REAL(wp), PARAMETER ::  pref_i  = 1.0_wp / 100000.0_wp   !< Inverse Reference Pressure (1/Pa)
172    REAL(wp), PARAMETER ::  r_cp    = 0.286_wp               !< R / cp (exponent for potential temperature)
173
174
175    SAVE
176
177
178!   
179!-- Checks Input parameters
180    INTERFACE chem_emissions_check_parameters
181       MODULE PROCEDURE chem_emissions_check_parameters
182    END INTERFACE chem_emissions_check_parameters
183!
184!-- Matching Emissions actions 
185    INTERFACE chem_emissions_match
186       MODULE PROCEDURE chem_emissions_match
187    END INTERFACE chem_emissions_match
188!
189!-- Initialization actions 
190    INTERFACE chem_emissions_init
191       MODULE PROCEDURE chem_emissions_init
192    END INTERFACE chem_emissions_init
193!
194!-- Setup of Emissions
195    INTERFACE chem_emissions_setup
196       MODULE PROCEDURE chem_emissions_setup
197    END INTERFACE chem_emissions_setup
198
199
200   
201    PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup
202!
203!-- Public Variables
204    PUBLIC con_factor, len_index, len_index_pm, len_index_voc
205
206 CONTAINS
207
208!------------------------------------------------------------------------------!
209! Description:
210! ------------
211!> Routine for checking input parameters
212!------------------------------------------------------------------------------!
213 SUBROUTINE chem_emissions_check_parameters
214
215
216    IMPLICIT NONE
217
218    TYPE(chem_emis_att_type) ::  emt
219
220    !
221    !-- Check Emission Species Number equal to number of passed names for the chemistry species:
222    IF ( SIZE(emt%species_name) /= emt%nspec  )  THEN
223
224       message_string = 'Numbers of input emission species names and number of species'     //      &
225                         'for which emission values are given do not match'                 
226       CALL message( 'chem_emissions_check_parameters', 'CM0437', 2, 2, 0, 6, 0 ) 
227           
228    ENDIF
229   
230
231 END SUBROUTINE chem_emissions_check_parameters
232
233!------------------------------------------------------------------------------!
234! Description:
235! ------------
236!> Matching the chemical species indices. The routine checks what are the indices of the emission input species
237!> and the corresponding ones of the model species. The routine gives as output a vector containing the number
238!> of common species: it is important to note that while the model species are distinct, their values could be
239!> given to a single species in input: for example, in the case of NO2 and NO, values may be passed in input as
240!> NOx values.
241!------------------------------------------------------------------------------!
242SUBROUTINE chem_emissions_match( emt_att,len_index )   
243
244
245    INTEGER(iwp), INTENT(INOUT)             ::  len_index   !< Variable where to store the number of common species between the input dataset and the model species   
246
247    TYPE(chem_emis_att_type), INTENT(INOUT) ::  emt_att     !< Chemistry Emission Array containing information for all the input chemical emission species
248   
249    INTEGER(iwp) ::  ind_mod, ind_inp                       !< Parameters for cycling through chemical model and input species
250    INTEGER(iwp) ::  nspec_emis_inp                         !< Variable where to store the number of the emission species in input
251    INTEGER(iwp) ::  ind_voc                                !< Indices to check whether a split for voc should be done
252    INTEGER(iwp) ::  ispec                                  !< index for cycle over effective number of emission species
253
254
255    IF ( debug_output )  CALL debug_message( 'chem_emissions_match', 'start' )
256
257    !
258    !-- Number of input emission species.
259    nspec_emis_inp=emt_att%nspec 
260
261    !
262    !-- Check the emission mode: DEFAULT, PRE-PROCESSED or PARAMETERIZED
263    SELECT CASE( TRIM( mode_emis ) ) 
264
265       !
266       !-- PRE-PROCESSED mode
267       CASE ( "PRE-PROCESSED" )
268
269          len_index = 0
270          len_index_voc = 0
271
272          IF ( nvar > 0 .AND. (nspec_emis_inp > 0) )  THEN
273             !
274             !-- Cycle over model species
275             DO ind_mod = 1,  nvar 
276                !
277                !-- Cycle over input species 
278                DO ind_inp = 1, nspec_emis_inp
279
280                   !
281                   !-- Check for VOC Species 
282                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN       
283                      DO ind_voc = 1, emt_att%nvoc
284             
285                         IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
286                            len_index = len_index + 1
287                            len_index_voc = len_index_voc + 1
288                         ENDIF
289                      END DO
290                   ENDIF
291                   !
292                   !-- Other Species
293                   IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
294                      len_index = len_index + 1
295                   ENDIF
296                ENDDO
297             ENDDO
298
299             !
300             !-- Allocate array for storing the indices of the matched species
301             IF ( len_index > 0 )  THEN
302 
303                ALLOCATE( match_spec_input(len_index) ) 
304 
305                ALLOCATE( match_spec_model(len_index) )
306
307                IF ( len_index_voc > 0 )  THEN
308                   !
309                   !-- contains indices of the VOC model species
310                   ALLOCATE( match_spec_voc_model(len_index_voc) )
311                   !
312                   !-- contains the indices of different values of VOC composition of input variable VOC_composition
313                   ALLOCATE( match_spec_voc_input(len_index_voc) )
314
315                ENDIF
316
317                !
318                !-- pass the species indices to declared arrays
319                len_index = 0
320
321                !
322                !-- Cycle over model species
323                DO ind_mod = 1, nvar 
324                   !
325                   !-- Cycle over Input species 
326                   DO ind_inp = 1, nspec_emis_inp
327                      !
328                      !-- VOCs
329                      IF ( TRIM(emt_att%species_name(ind_inp) ) == "VOC" .AND.    &
330                           ALLOCATED(match_spec_voc_input) )  THEN
331                         
332                         DO ind_voc= 1, emt_att%nvoc
333                            IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
334                               len_index = len_index + 1
335                               len_index_voc = len_index_voc + 1
336                       
337                               match_spec_input(len_index) = ind_inp
338                               match_spec_model(len_index) = ind_mod
339
340                               match_spec_voc_input(len_index_voc) = ind_voc
341                               match_spec_voc_model(len_index_voc) = ind_mod                         
342                            ENDIF
343                         END DO
344                      ENDIF
345
346                      !
347                      !-- Other Species
348                      IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
349                         len_index = len_index + 1
350                         match_spec_input(len_index) = ind_inp
351                         match_spec_model(len_index) = ind_mod
352                      ENDIF
353                   END DO
354                END DO
355
356             ELSE
357                !
358                !-- in case there are no species matching
359                message_string = 'Non of given emission species'            //         &
360                                 ' matches'                                //          &
361                                 ' model chemical species:'                //          &
362                                 ' Emission routine is not called' 
363                CALL message( 'chem_emissions_matching', 'CM0438', 0, 0, 0, 6, 0 )
364             ENDIF
365 
366          ELSE
367
368             !
369             !-- either spc_names is zero or nspec_emis_inp is not allocated
370             message_string = 'Array of Emission species not allocated:'             //          &
371                              ' Either no emission species are provided as input or'  //         &
372                              ' no chemical species are used by PALM:'                //         &
373                              ' Emission routine is not called'                 
374             CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 ) 
375
376          ENDIF
377
378       !
379       !-- DEFAULT mode
380       CASE ("DEFAULT")
381
382          len_index = 0          !<  index for TOTAL number of species   
383          len_index_voc = 0      !<  index for TOTAL number of VOCs
384          len_index_pm = 3       !<  index for TOTAL number of PMs: PM1, PM2.5, PM10.
385 
386          IF ( nvar > 0 .AND. nspec_emis_inp > 0 )  THEN
387
388             !
389             !-- Cycle over model species
390             DO ind_mod = 1, nvar
391                !
392                !-- Cycle over input species
393                DO ind_inp = 1, nspec_emis_inp
394
395                   !
396                   !-- Check for VOC Species 
397                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN
398                      DO ind_voc= 1, emt_att%nvoc
399                       
400                         IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
401                            len_index = len_index + 1
402                            len_index_voc = len_index_voc + 1
403                         ENDIF
404                         
405                      END DO
406                   ENDIF
407
408                   !
409                   !-- PMs: There is one input species name for all PM
410                   !-- This variable has 3 dimensions, one for PM1, PM2.5 and PM10
411                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" )  THEN
412                      !
413                      !-- PM1
414                      IF ( TRIM( spc_names(ind_mod) ) == "PM1" )  THEN
415                         len_index = len_index + 1
416                      !
417                      !-- PM2.5
418                      ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" )  THEN
419                         len_index = len_index + 1
420                      !
421                      !-- PM10
422                      ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" )  THEN
423                         len_index = len_index + 1
424                      ENDIF
425                   ENDIF
426
427                   !
428                   !-- NOx: NO2 and NO   
429                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOx" )  THEN 
430                      !
431                      !-- NO
432                      IF ( TRIM( spc_names(ind_mod) ) == "NO" )  THEN
433                         len_index = len_index + 1
434                      !
435                      !-- NO2
436                      ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" )  THEN
437                         len_index = len_index + 1
438                      ENDIF
439                   ENDIF
440
441                   !
442                   !-- SOx: SO2 and SO4
443                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOx" )  THEN
444                      !
445                      !-- SO2
446                      IF ( TRIM( spc_names(ind_mod) ) == "SO2" )  THEN
447                         len_index = len_index + 1
448                      !
449                      !-- SO4
450                      ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" )  THEN
451                         len_index = len_index + 1
452                      ENDIF
453                   ENDIF
454
455                   !
456                   !-- Other Species
457                   IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
458                      len_index = len_index + 1
459                   ENDIF
460                END DO
461             END DO
462
463
464             !
465             !-- Allocate arrays
466             IF ( len_index > 0 )  THEN
467
468                ALLOCATE( match_spec_input(len_index) )
469                ALLOCATE( match_spec_model(len_index) )
470
471                IF ( len_index_voc > 0 )  THEN
472                   !
473                   !-- Contains indices of the VOC model species
474                   ALLOCATE( match_spec_voc_model(len_index_voc) )
475                   !
476                   !-- Contains the indices of different values of VOC composition
477                   !-- of input variable VOC_composition
478                   ALLOCATE( match_spec_voc_input(len_index_voc) )                                                 
479                ENDIF
480
481                !
482                !-- Pass the species indices to declared arrays
483                len_index = 0
484                len_index_voc = 0
485               
486                DO ind_mod = 1, nvar
487                   DO ind_inp = 1, nspec_emis_inp 
488
489                      !
490                      !-- VOCs
491                      IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" .AND.   &
492                         ALLOCATED(match_spec_voc_input) )  THEN     
493                         DO ind_voc= 1, emt_att%nvoc
494                            IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
495                               len_index = len_index + 1
496                               len_index_voc = len_index_voc + 1
497                       
498                               match_spec_input(len_index) = ind_inp
499                               match_spec_model(len_index) = ind_mod
500
501                               match_spec_voc_input(len_index_voc) = ind_voc
502                               match_spec_voc_model(len_index_voc) = ind_mod                         
503                            ENDIF
504                         END DO
505                      ENDIF
506
507                      !
508                      !-- PMs
509                      IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" )  THEN
510                         !
511                         !-- PM1
512                         IF ( TRIM( spc_names(ind_mod) ) == "PM1" )  THEN
513                            len_index = len_index + 1
514
515                            match_spec_input(len_index) = ind_inp
516                            match_spec_model(len_index) = ind_mod
517                         !
518                         !-- PM2.5
519                         ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" )  THEN
520                            len_index = len_index + 1
521
522                            match_spec_input(len_index) = ind_inp
523                            match_spec_model(len_index) = ind_mod
524                         !
525                         !-- PM10
526                         ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" )  THEN
527                            len_index = len_index + 1
528   
529                            match_spec_input(len_index) = ind_inp
530                            match_spec_model(len_index) = ind_mod
531 
532                         ENDIF
533                      ENDIF
534
535                      !
536                      !-- NOx
537                      IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOx" )  THEN
538                         !
539                         !-- NO
540                         IF ( TRIM( spc_names(ind_mod) ) == "NO" )  THEN
541                            len_index = len_index + 1
542
543                            match_spec_input(len_index) = ind_inp
544                            match_spec_model(len_index) = ind_mod
545                         !
546                         !-- NO2
547                         ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" )  THEN
548                            len_index = len_index + 1
549
550                            match_spec_input(len_index) = ind_inp
551                            match_spec_model(len_index) = ind_mod
552 
553                         ENDIF
554                      ENDIF
555
556                      !
557                      !-- SOx
558                      IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOx" ) THEN
559                         !
560                         !-- SO2
561                         IF ( TRIM( spc_names(ind_mod) ) == "SO2" )  THEN
562                            len_index = len_index + 1
563
564                            match_spec_input(len_index) = ind_inp
565                            match_spec_model(len_index) = ind_mod
566 
567                         !
568                         !-- SO4
569                         ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" )  THEN
570                            len_index = len_index + 1
571
572                            match_spec_input(len_index) = ind_inp
573                            match_spec_model(len_index) = ind_mod
574 
575                         ENDIF
576                      ENDIF
577
578                      !
579                      !-- Other Species
580                      IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
581                         len_index = len_index + 1
582                           
583                         match_spec_input(len_index) = ind_inp
584                         match_spec_model(len_index) = ind_mod
585                      ENDIF
586                   END DO
587                END DO
588
589             ELSE
590
591                message_string = 'Non of given Emission Species'            //         &
592                                 ' matches'                                //          &
593                                 ' model chemical species'                 //          &
594                                 ' Emission routine is not called'         
595                CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 ) 
596
597             ENDIF
598
599          ELSE
600
601             message_string = 'Array of Emission species not allocated: '            //          &
602                              ' Either no emission species are provided as input or'  //         &
603                              ' no chemical species are used by PALM:'                //         &
604                              ' Emission routine is not called'                                   
605             CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 ) 
606 
607          ENDIF
608 
609       !
610       !-- PARAMETERIZED mode
611       CASE ("PARAMETERIZED")
612
613          len_index = 0
614
615          IF ( nvar > 0 .AND. nspec_emis_inp > 0 )  THEN
616
617             !
618             !-- Cycle over model species
619             DO ind_mod = 1, nvar
620                ind_inp = 1
621                DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )    !< 'novalue' is the default 
622                   IF ( TRIM( surface_csflux_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
623                      len_index = len_index + 1
624                   ENDIF
625                   ind_inp = ind_inp + 1
626                ENDDO
627             ENDDO
628
629             IF ( len_index > 0 )  THEN
630
631                !
632                !-- Allocation of Arrays of the matched species
633                ALLOCATE( match_spec_input(len_index) ) 
634 
635                ALLOCATE( match_spec_model(len_index) )
636
637                !
638                !-- Pass the species indices to declared arrays   
639                len_index = 0
640
641                DO ind_mod = 1, nvar 
642                   ind_inp = 1
643                   DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )     
644                      IF ( TRIM( surface_csflux_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
645                         len_index = len_index + 1
646                         match_spec_input(len_index) = ind_inp
647                         match_spec_model(len_index) = ind_mod
648                      ENDIF
649                   ind_inp = ind_inp + 1
650                   END DO
651                END DO
652
653                !
654                !-- Check
655                DO ispec = 1, len_index
656
657                   IF ( emiss_factor_main(match_spec_input(ispec) ) < 0 .AND.    &
658                        emiss_factor_side(match_spec_input(ispec) ) < 0 )  THEN
659
660                      message_string = 'PARAMETERIZED emissions mode selected:'            //           &
661                                       ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES'        //          &
662                                       ' but values of scaling factors for street types'    //          &
663                                       ' emiss_factor_main AND emiss_factor_side'           //          &
664                                       ' not provided for each of the emissions species'    //          &
665                                       ' or not provided at all: PLEASE set a finite value' //          &
666                                       ' for these parameters in the chemistry namelist'         
667                      CALL message( 'chem_emissions_matching', 'CM0442', 2, 2, 0, 6, 0 ) 
668                   ENDIF
669                END DO
670
671
672             ELSE
673               
674                message_string = 'Non of given Emission Species'            //          &
675                                 ' matches'                                //           &
676                                 ' model chemical species'                 //           &
677                                 ' Emission routine is not called'         
678                CALL message( 'chem_emissions_matching', 'CM0443', 0, 0, 0, 6, 0 ) 
679             ENDIF
680
681          ELSE
682     
683             message_string = 'Array of Emission species not allocated: '            //           &
684                              ' Either no emission species are provided as input or'  //          &
685                              ' no chemical species are used by PALM.'                //          &
686                              ' Emission routine is not called'                                   
687             CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 ) 
688 
689          ENDIF             
690
691
692       !
693       !-- If emission module is switched on but mode_emis is not specified or it is given the wrong name
694       CASE DEFAULT
695
696          message_string = 'Emission Module switched ON, but'            //                         &
697                           ' either no emission mode specified or incorrectly given :'  //          &
698                           ' please, pass the correct value to the namelist parameter "mode_emis"'                 
699          CALL message( 'chem_emissions_matching', 'CM0445', 2, 2, 0, 6, 0 )             
700
701       END SELECT
702
703       IF ( debug_output )  CALL debug_message( 'chem_emissions_match', 'end' )
704
705 END SUBROUTINE chem_emissions_match
706
707 
708!------------------------------------------------------------------------------!
709! Description:
710! ------------
711!> Initialization:
712!> Netcdf reading, arrays allocation and first calculation of cssws
713!> fluxes at timestep 0
714!------------------------------------------------------------------------------!
715 SUBROUTINE chem_emissions_init
716
717    USE netcdf_data_input_mod,                                                 &
718        ONLY:  chem_emis, chem_emis_att
719   
720    IMPLICIT NONE
721 
722    INTEGER(iwp)                :: ispec                                            !< running index
723
724!   
725!-- Actions for initial runs
726!  IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
727!--    ...
728!   
729!
730!-- Actions for restart runs
731!  ELSE
732!--    ...
733!
734!  ENDIF
735
736
737  IF ( debug_output )  CALL debug_message( 'chem_emissions_init', 'start' )
738
739  !
740  !-- Matching
741  CALL  chem_emissions_match( chem_emis_att, nspec_out )
742 
743  IF ( nspec_out == 0 )  THEN
744 
745     emission_output_required = .FALSE.
746
747  ELSE
748
749     emission_output_required = .TRUE.
750
751
752     !
753     !-- Set molecule masses'
754     ALLOCATE( chem_emis_att%xm(nspec_out) )
755
756     DO ispec = 1, nspec_out
757        SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
758           CASE ( 'SO2' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 2        !< kg/mole
759           CASE ( 'SO4' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 4        !< kg/mole
760           CASE ( 'NO' ); chem_emis_att%xm(ispec) = xm_N + xm_O             !< kg/mole
761           CASE ( 'NO2' ); chem_emis_att%xm(ispec) = xm_N + xm_O * 2        !< kg/mole   
762           CASE ( 'NH3' ); chem_emis_att%xm(ispec) = xm_N + xm_H * 3        !< kg/mole
763           CASE ( 'CO'  ); chem_emis_att%xm(ispec) = xm_C + xm_O            !< kg/mole
764           CASE ( 'CO2' ); chem_emis_att%xm(ispec) = xm_C + xm_O * 2        !< kg/mole
765           CASE ( 'CH4' ); chem_emis_att%xm(ispec) = xm_C + xm_H * 4        !< kg/mole     
766           CASE ( 'HNO3' ); chem_emis_att%xm(ispec) = xm_H + xm_N + xm_O*3  !< kg/mole 
767           CASE DEFAULT
768              chem_emis_att%xm(ispec) = 1.0_wp
769        END SELECT
770     ENDDO
771
772   
773     !
774     !-- assign emission values
775     SELECT CASE ( TRIM( mode_emis ) )   
776
777
778        !
779        !-- PRE-PROCESSED case
780        CASE ( "PRE-PROCESSED" )
781
782           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution(nzb:nzt+1,0:ny,0:nx,nspec_out) ) 
783 
784           !
785           !-- Get emissions at the first time step
786           CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out )
787
788        !
789        !-- Default case
790        CASE ( "DEFAULT" )
791
792           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution(1,0:ny,0:nx,nspec_out) ) 
793
794           !
795           !-- Get emissions at the first time step
796           CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out )
797
798        !
799        !-- PARAMETERIZED case
800        CASE ( "PARAMETERIZED" )
801
802           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution(1,0:ny,0:nx,nspec_out) )
803
804           !
805           !-- Get emissions at the first time step
806           CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out)
807
808     END SELECT
809
810  ENDIF   
811
812  IF ( debug_output )  CALL debug_message( 'chem_emissions_init', 'end' )
813
814 END SUBROUTINE chem_emissions_init
815
816
817
818!------------------------------------------------------------------------------!
819! Description:
820! ------------
821!> Routine for Update of Emission values at each timestep
822!-------------------------------------------------------------------------------!
823
824 SUBROUTINE chem_emissions_setup( emt_att, emt, nspec_out )
825 
826   USE surface_mod,                                                  &
827       ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
828   USE netcdf_data_input_mod,                                        &
829       ONLY:  street_type_f
830   USE arrays_3d,                                                    &       
831       ONLY: hyp, pt 
832
833   
834 IMPLICIT NONE
835 
836
837    TYPE(chem_emis_att_type), INTENT(INOUT) ::  emt_att                         !< variable to store emission information                                                                           
838
839    TYPE(chem_emis_val_type), INTENT(INOUT), ALLOCATABLE, DIMENSION(:) ::  emt  !< variable to store emission input values,
840                                                                                !< depending on the considered species
841
842    INTEGER,INTENT(IN) ::  nspec_out                                            !< Output of matching routine with number
843                                                                                !< of matched species
844
845    INTEGER(iwp) ::  i                                                          !< running index for grid in x-direction
846    INTEGER(iwp) ::  j                                                          !< running index for grid in y-direction
847    INTEGER(iwp) ::  k                                                          !< running index for grid in z-direction
848    INTEGER(iwp) ::  m                                                          !< running index for horizontal surfaces
849   
850    INTEGER(iwp) ::  icat                                                       !< Index for number of categories
851    INTEGER(iwp) ::  ispec                                                      !< index for number of species
852    INTEGER(iwp) ::  i_pm_comp                                                  !< index for number of PM components
853    INTEGER(iwp) ::  ivoc                                                       !< Index for number of VOCs
854
855    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  delta_emis                       
856    REAL(wp), ALLOCATABLE, DIMENSION(:)   ::  time_factor                       !< factor for time scaling of emissions
857    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  emis
858
859    REAL(wp), DIMENSION(24) :: par_emis_time_factor                             !< time factors for the parameterized mode:
860                                                                               !< fixed houlry profile for example day
861    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  conv_to_ratio          !< factor used for converting input
862                                                                                !< to concentration ratio
863    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  tmp_temp
864
865    !
866    !-- CONVERSION FACTORS: TIME 
867    REAL(wp), PARAMETER ::  s_per_hour = 3600.0                       !< number of sec per hour (s)/(hour)   
868    REAL(wp), PARAMETER ::  s_per_day = 86400.0                       !< number of sec per day (s)/(day) 
869    REAL(wp), PARAMETER ::  hour_per_year = 8760.0                    !< number of hours in a year of 365 days 
870    REAL(wp), PARAMETER ::  hour_per_day = 24.0                       !< number of hours in a day
871   
872    REAL(wp), PARAMETER ::  hour_to_s = 1/s_per_hour                  !< conversion from hours to seconds (s/hour) ~ 0.2777778     
873    REAL(wp), PARAMETER ::  day_to_s = 1/s_per_day                    !< conversion from day to seconds (s/day) ~ 1.157407e-05
874    REAL(wp), PARAMETER ::  year_to_s = 1/(s_per_hour*hour_per_year)  !< conversion from year to sec (s/year) ~ 3.170979e-08
875    !
876    !-- CONVERSION FACTORS: WEIGHT 
877    REAL(wp), PARAMETER ::  tons_to_kg = 100                          !< Conversion from tons to kg (kg/tons)   
878    REAL(wp), PARAMETER ::  g_to_kg = 0.001                           !< Conversion from g to kg (kg/g)
879    REAL(wp), PARAMETER ::  miug_to_kg = 0.000000001                  !< Conversion from g to kg (kg/g)
880    !
881    !-- CONVERSION FACTORS: fraction to ppm
882    REAL(wp), PARAMETER ::  ratio2ppm  = 1.0e06 
883    !------------------------------------------------------   
884
885    IF ( emission_output_required )  THEN
886
887       !
888       !-- Set emis_dt 
889       IF ( call_chem_at_all_substeps )  THEN
890
891          dt_emis = dt_3d * weight_pres(intermediate_timestep_count)
892
893       ELSE
894
895          dt_emis = dt_3d
896
897       ENDIF
898
899
900       !
901       !-- Conversion of units to the ones employed in PALM
902       !-- In PARAMETERIZED mode no conversion is performed: in this case input units are fixed
903
904        IF ( TRIM( mode_emis ) == "DEFAULT" .OR. TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
905
906          SELECT CASE ( TRIM( emt_att%units ) )
907             !
908             !-- kilograms
909             CASE ( 'kg/m2/s', 'KG/M2/S' ) 
910               
911                con_factor=1
912
913             CASE ('kg/m2/hour', 'KG/M2/HOUR' )
914
915                con_factor=hour_to_s
916
917             CASE ( 'kg/m2/day', 'KG/M2/DAY' ) 
918               
919                con_factor=day_to_s
920
921             CASE ( 'kg/m2/year', 'KG/M2/YEAR' )
922             
923                con_factor=year_to_s
924
925             !
926             !-- Tons
927             CASE ( 'ton/m2/s', 'TON/M2/S' ) 
928               
929                con_factor=tons_to_kg
930
931             CASE ( 'ton/m2/hour', 'TON/M2/HOUR' ) 
932               
933                con_factor=tons_to_kg*hour_to_s
934
935             CASE ( 'ton/m2/year', 'TON/M2/YEAR' ) 
936               
937                con_factor=tons_to_kg*year_to_s
938
939             !
940             !-- Grams
941             CASE ( 'g/m2/s', 'G/M2/S' )
942         
943                con_factor=g_to_kg
944
945             CASE ( 'g/m2/hour', 'G/M2/HOUR' ) 
946
947                con_factor=g_to_kg*hour_to_s
948
949             CASE ( 'g/m2/year', 'G/M2/YEAR' )
950 
951                con_factor=g_to_kg*year_to_s
952
953             !
954             !-- Micrograms
955             CASE ( 'micrograms/m2/s', 'MICROGRAMS/M2/S' )
956 
957                con_factor=miug_to_kg
958
959             CASE ( 'micrograms/m2/hour', 'MICROGRAMS/M2/HOUR' ) 
960 
961                con_factor=miug_to_kg*hour_to_s
962
963             CASE ( 'micrograms/m2/year', 'MICROGRAMS/M2/YEAR' )
964 
965                con_factor=miug_to_kg*year_to_s
966
967             CASE DEFAULT 
968                message_string = 'The Units of the provided emission input' // &
969                                 ' are not the ones required by PALM-4U: please check '      // &
970                                 ' emission module documentation.'                                 
971                CALL message( 'chem_emissions_setup', 'CM0446', 2, 2, 0, 6, 0 )
972
973          END SELECT
974
975         
976       ENDIF
977
978       !
979       !-- Conversion factor to convert  kg/m**2/s to ppm/s
980       DO i = nxl, nxr
981          DO j = nys, nyn
982             !
983             !-- Derive Temperature from Potential Temperature
984             tmp_temp(nzb:nzt+1,j,i) = pt(nzb:nzt+1,j,i) * ( hyp(nzb:nzt+1) * pref_i )**r_cp
985             
986             
987             !>  We need to pass to cssws <- (ppm/s) * dz
988             !>  Input is Nmole/(m^2*s)
989             !>  To go to ppm*dz multiply the input by (m**2/N)*dz
990             !>  (m**2/N)*dz == V/N
991             !>  V/N = RT/P
992             !>  m**3/Nmole               (J/mol)*K^-1           K                      Pa         
993             conv_to_ratio(nzb:nzt+1,j,i) = ( (Rgas * tmp_temp(nzb:nzt+1,j,i)) / ((hyp(nzb:nzt+1))) ) 
994          ENDDO
995       ENDDO
996
997
998       !
999       !-- Initialize
1000       emis_distribution(:,nys:nyn,nxl:nxr,:) = 0.0_wp
1001
1002 
1003       !
1004       !-- PRE-PROCESSED MODE
1005       IF ( TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
1006
1007          !
1008          !-- Update time indices
1009          CALL time_preprocessed_indices( index_hh )
1010
1011       ELSEIF ( TRIM( mode_emis ) == "DEFAULT" )  THEN
1012
1013          !
1014          !-- Allocate array where to store temporary emission values     
1015          IF ( .NOT. ALLOCATED(emis) ) ALLOCATE( emis(nys:nyn,nxl:nxr) )
1016          !
1017          !-- Allocate time factor per category
1018          ALLOCATE( time_factor(emt_att%ncat) ) 
1019          !
1020          !-- Read-in hourly emission time factor
1021          IF ( TRIM( time_fac_type ) == "HOUR" )  THEN
1022
1023             !
1024             !-- Update time indices
1025             CALL time_default_indices( month_of_year, day_of_month, hour_of_day, index_hh )
1026             !
1027             !-- Check if the index is less or equal to the temporal dimension of HOURLY emission files             
1028             IF ( index_hh <= SIZE( emt_att%hourly_emis_time_factor(1,:) ) )  THEN
1029                !
1030                !-- Read-in the correspondant time factor             
1031                time_factor(:) = emt_att%hourly_emis_time_factor(:,index_hh)     
1032
1033             ELSE
1034
1035                message_string = 'The "HOUR" time-factors in the DEFAULT mode '           // &
1036                              ' are not provided for each hour of the total simulation time'     
1037                CALL message( 'chem_emissions_setup', 'CM0448', 2, 2, 0, 6, 0 ) 
1038
1039             ENDIF
1040          !
1041          !-- Read-in MDH emissions time factors
1042          ELSEIF ( TRIM( time_fac_type ) == "MDH" )  THEN
1043
1044             !
1045             !-- Update time indices     
1046             CALL time_default_indices( daytype_mdh, month_of_year, day_of_month,     &
1047                  hour_of_day, index_mm, index_dd,index_hh )
1048
1049             !
1050             !-- Check if the index is less or equal to the temporal dimension of MDH emission files             
1051             IF ( ( index_hh + index_dd + index_mm) <= SIZE( emt_att%mdh_emis_time_factor(1,:) ) )  THEN
1052                !
1053                !-- Read-in the correspondant time factor             
1054                time_factor(:) = emt_att%mdh_emis_time_factor(:,index_mm) * emt_att%mdh_emis_time_factor(:,index_dd) *   &
1055                                                                         emt_att%mdh_emis_time_factor(:,index_hh)
1056     
1057             ELSE
1058
1059                message_string = 'The "MDH" time-factors in the DEFAULT mode '           //  &
1060                              ' are not provided for each hour/day/month  of the total simulation time'     
1061                CALL message( 'chem_emissions_setup', 'CM0449', 2, 2, 0, 6, 0 )
1062
1063             ENDIF
1064
1065          ELSE
1066                 
1067             message_string = 'In the DEFAULT mode the time factor'           //  &
1068                              ' has to be defined in the NAMELIST'     
1069             CALL message( 'chem_emissions_setup', 'CM0450', 2, 2, 0, 6, 0 )
1070         
1071          ENDIF
1072
1073       !
1074       !-- PARAMETERIZED MODE
1075       ELSEIF ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
1076       
1077          !
1078          !-- assign constant values of time factors, diurnal time profile for traffic sector
1079          par_emis_time_factor( : ) =  &
1080            (/ 0.009, 0.004, 0.004, 0.009, 0.029, 0.039, 0.056, 0.053, 0.051, 0.051, 0.052, 0.055,  &
1081               0.059, 0.061, 0.064, 0.067, 0.069, 0.069, 0.049, 0.039, 0.039, 0.029, 0.024, 0.019 /)
1082         
1083          IF ( .NOT. ALLOCATED( time_factor ) ) ALLOCATE( time_factor(1) )
1084
1085          !
1086          !-- Get time-factor for specific hour
1087          index_hh = hour_of_day
1088
1089          time_factor(1) = par_emis_time_factor(index_hh)
1090
1091       ENDIF
1092
1093       
1094       !
1095       !--  Emission distribution calculation
1096
1097       !
1098       !-- PARAMETERIZED case
1099       IF ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
1100
1101          DO ispec = 1, nspec_out
1102
1103             !
1104             !-- Units are micromoles/m**2*day (or kilograms/m**2*day for PMs)
1105             emis_distribution(1,nys:nyn,nxl:nxr,ispec) = surface_csflux(match_spec_input(ispec)) *  &
1106                                                           time_factor(1) * hour_to_s
1107
1108          ENDDO 
1109
1110       !
1111       !-- PRE-PROCESSED case
1112       ELSEIF ( TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
1113
1114          !
1115          !-- Cycle over species:
1116          !-- nspec_out represents the number of species in common between the emission input data
1117          !-- and the chemistry mechanism used
1118          DO ispec=1,nspec_out 
1119   
1120             emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emt(match_spec_input(ispec))%                                 &
1121                                                            preproc_emission_data(index_hh,1,nys+1:nyn+1,nxl+1:nxr+1) * &
1122                                                                con_factor
1123         
1124          ENDDO
1125         
1126       !
1127       !-- DEFAULT case
1128       ELSEIF ( TRIM( mode_emis ) == "DEFAULT" )  THEN
1129
1130          !
1131          !-- Allocate array for the emission value corresponding to a specific category and time factor
1132          ALLOCATE( delta_emis(nys:nyn,nxl:nxr) ) 
1133
1134          !
1135          !-- Cycle over categories
1136          DO icat = 1, emt_att%ncat
1137 
1138             !
1139             !-- Cycle over Species
1140             !-- nspec_out represents the number of species in common between the emission input data
1141             !-- and the chemistry mechanism used
1142             DO ispec = 1, nspec_out
1143
1144                emis(nys:nyn,nxl:nxr) = emt(match_spec_input(ispec))%default_emission_data(icat,nys+1:nyn+1,nxl+1:nxr+1)
1145
1146
1147                !
1148                !-- NOx
1149                IF ( TRIM( spc_names(match_spec_model(ispec)) ) == "NO" )  THEN
1150               
1151                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &                 !<  kg/m2*s
1152                                                  emt_att%nox_comp(icat,1) * con_factor * hour_per_day
1153
1154                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
1155                                                                 delta_emis(nys:nyn,nxl:nxr)
1156
1157                ELSEIF ( TRIM( spc_names(match_spec_model(ispec)) ) == "NO2" )  THEN
1158
1159                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &                 !<  kg/m2*s
1160                                                  emt_att%nox_comp(icat,2) * con_factor * hour_per_day
1161
1162                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
1163                                                                 delta_emis(nys:nyn,nxl:nxr)
1164 
1165                !
1166                !-- SOx
1167                ELSEIF ( TRIM( spc_names(match_spec_model(ispec)) ) == "SO2" )  THEN
1168                 
1169                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &                 !<  kg/m2*s
1170                                                 emt_att%sox_comp(icat,1) * con_factor * hour_per_day
1171
1172                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
1173                                                                 delta_emis(nys:nyn,nxl:nxr)
1174
1175                ELSEIF ( TRIM( spc_names(match_spec_model(ispec)) ) == "SO4" )  THEN
1176                 
1177                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &                 !<  kg/m2*s
1178                                                 emt_att%sox_comp(icat,2) * con_factor * hour_per_day
1179
1180                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
1181                                                                 delta_emis(nys:nyn,nxl:nxr)
1182 
1183
1184                !
1185                !-- PMs
1186                !-- PM1
1187                ELSEIF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM1" )  THEN
1188
1189                   !
1190                   !-- Cycle over PM1 components
1191                   DO i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,1) )
1192
1193                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &               !<  kg/m2*s
1194                                                     emt_att%pm_comp(icat,i_pm_comp,1) * con_factor * hour_per_day 
1195
1196                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
1197                                                                    delta_emis(nys:nyn,nxl:nxr)
1198                   ENDDO
1199
1200                !
1201                !-- PM2.5
1202                ELSEIF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM25" )  THEN
1203
1204                   !
1205                   !-- Cycle over PM2.5 components
1206                   DO i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,2) )
1207
1208                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &               !<  kg/m2*s
1209                                                     emt_att%pm_comp(icat,i_pm_comp,2) * con_factor * hour_per_day 
1210
1211                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
1212                                                                    delta_emis(nys:nyn,nxl:nxr)
1213 
1214                   ENDDO
1215
1216                !
1217                !-- PM10
1218                ELSEIF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM10" )  THEN
1219
1220                   !
1221                   !-- Cycle over PM10 components
1222                   DO i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,3) ) 
1223
1224                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &               !<  kg/m2*s
1225                                                     emt_att%pm_comp(icat,i_pm_comp,3) * con_factor * hour_per_day 
1226
1227                      emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+ &
1228                                                                 delta_emis(nys:nyn,nxl:nxr) 
1229
1230                   ENDDO
1231
1232                !
1233                !-- VOCs
1234                ELSEIF ( SIZE( match_spec_voc_input ) > 0 )  THEN
1235
1236                   DO ivoc = 1, SIZE( match_spec_voc_input )
1237
1238                      IF ( TRIM( spc_names(match_spec_model(ispec)) ) == TRIM( emt_att%voc_name(ivoc) ) )  THEN   
1239
1240                         delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &
1241                                                        emt_att%voc_comp(icat,match_spec_voc_input(ivoc)) * &
1242                                                         con_factor * hour_per_day
1243
1244                         emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
1245                                                                       delta_emis(nys:nyn,nxl:nxr) 
1246
1247                      ENDIF                       
1248
1249                   ENDDO
1250               
1251                !
1252                !-- any other species
1253                ELSE
1254
1255                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &
1256                                                  con_factor * hour_per_day
1257 
1258                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
1259                                                                 delta_emis(nys:nyn,nxl:nxr) 
1260
1261                ENDIF
1262               
1263                emis(:,:)= 0 
1264               
1265             ENDDO
1266             
1267             delta_emis(:,:)=0 
1268         
1269          ENDDO
1270
1271       ENDIF
1272
1273       
1274!
1275!-- Cycle to transform x,y coordinates to the one of surface_mod and to assign emission values to cssws
1276!
1277!-- PARAMETERIZED mode
1278       !
1279       !-- Units of inputs are micromoles/(m**2*s)
1280       IF ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
1281
1282          IF ( street_type_f%from_file )  THEN
1283
1284!
1285!--          Streets are lsm surfaces, hence, no usm surface treatment required.
1286!--          However, urban surface may be initialized via default initialization
1287!--          in surface_mod, e.g. at horizontal urban walls that are at k == 0
1288!--          (building is lower than the first grid point). Hence, in order to
1289!--          have only emissions at streets, set the surfaces emissions to zero
1290!--          at urban walls.
1291             IF ( surf_usm_h%ns >=1 )  surf_usm_h%cssws = 0.0_wp
1292!
1293!--          Now, treat land-surfaces.
1294             DO  m = 1, surf_lsm_h%ns
1295                i = surf_lsm_h%i(m)
1296                j = surf_lsm_h%j(m)
1297                k = surf_lsm_h%k(m)
1298
1299                IF ( street_type_f%var(j,i) >= main_street_id  .AND. street_type_f%var(j,i) < max_street_id )  THEN
1300
1301                   !
1302                   !-- Cycle over matched species
1303                   DO  ispec = 1, nspec_out
1304
1305                      !
1306                      !-- PMs are already in kilograms
1307                      IF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM1 "  &
1308                          .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM25"  &
1309                          .OR. TRIM( spc_names(match_spec_model(ispec)) )=="PM10")  THEN
1310
1311                         !
1312                         !--           kg/(m^2*s) * kg/m^3
1313                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec)) * &
1314                                                                        emis_distribution(1,j,i,ispec) *            &  !< in kg/(m^2*s)
1315                                                                         rho_air(k)                                    !< in kg/m^3
1316                         
1317                      !
1318                      !-- Other Species
1319                      !-- Inputs are micromoles
1320                      ELSE
1321
1322                         !   
1323                         !--             ppm/s *m *kg/m^3               
1324                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec)) * &
1325                                                                        emis_distribution(1,j,i,ispec) *            &  !< in micromoles/(m^2*s)
1326                                                                         conv_to_ratio(k,j,i) *                     &  !< in m^3/Nmole     
1327                                                                          rho_air(k)                                   !< in kg/m^3
1328
1329                      ENDIF
1330
1331                   ENDDO
1332
1333
1334                ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND. street_type_f%var(j,i) < main_street_id )  THEN
1335
1336                   !
1337                   !-- Cycle over matched species
1338                   DO  ispec = 1, nspec_out
1339
1340                      !
1341                      !-- PMs are already in kilograms
1342                      IF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM1"  &
1343                          .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM25"  &
1344                          .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM10" )  THEN
1345
1346                         !
1347                         !--           kg/(m^2*s) * kg/m^3
1348                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * &
1349                                                                        emis_distribution(1,j,i,ispec) *            &  !< in kg/(m^2*s)
1350                                                                         rho_air(k)                                    !< in kg/m^3   
1351                      !
1352                      !-- Other species
1353                      !-- Inputs are micromoles
1354                      ELSE
1355                   
1356                         !   
1357                         !--             ppm/s *m *kg/m^3     
1358                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * &
1359                                                                        emis_distribution(1,j,i,ispec) *            &  !< in micromoles/(m^2*s)
1360                                                                         conv_to_ratio(k,j,i) *                     &  !< in m^3/Nmole       
1361                                                                          rho_air(k)                                   !< in kg/m^3   
1362                      ENDIF
1363
1364                   ENDDO
1365
1366                ELSE
1367
1368                   !
1369                   !-- If no street type is defined, then assign zero emission to all the species
1370                   surf_lsm_h%cssws(:,m) = 0.0_wp
1371
1372                ENDIF
1373
1374             ENDDO
1375
1376          ENDIF
1377
1378       !
1379       !-- For both DEFAULT and PRE-PROCESSED mode
1380       ELSE   
1381       
1382
1383          DO ispec = 1, nspec_out 
1384                   
1385!
1386!--          Default surfaces
1387             DO  m = 1, surf_def_h(0)%ns
1388
1389                i = surf_def_h(0)%i(m)
1390                j = surf_def_h(0)%j(m)
1391
1392                IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
1393
1394                   !
1395                   !-- PMs
1396                   IF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM1"  &
1397                        .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM25"  &
1398                          .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM10" )  THEN
1399               
1400                      !
1401                      !--            kg/(m^2*s) *kg/m^3                         kg/(m^2*s)                 
1402                      surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)*  &
1403                                                                        rho_air(nzb)                           !< in kg/m^3 
1404 
1405 
1406                   ELSE
1407
1408                      !
1409                      !-- VOCs
1410                      IF ( len_index_voc > 0 .AND. emt_att%species_name(match_spec_input(ispec)) == "VOC" )  THEN
1411                         !
1412                         !--           (ppm/s) * m * kg/m^3                        mole/(m^2/s)   
1413                         surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &
1414                                                                           conv_to_ratio(nzb,j,i) *         &  !< in m^3/mole
1415                                                                            ratio2ppm *                     &  !< in ppm
1416                                                                             rho_air(nzb)                      !< in kg/m^3 
1417
1418
1419                      !
1420                      !-- Other species
1421                      ELSE
1422
1423                         !
1424                         !--           (ppm/s) * m  * kg/m^3                    kg/(m^2/s)
1425                         surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &       
1426                                                                           ( 1.0_wp / emt_att%xm(ispec) ) * &  !< in mole/kg
1427                                                                             conv_to_ratio(nzb,j,i) *       &  !< in m^3/mole 
1428                                                                              ratio2ppm *                   &  !< in ppm
1429                                                                               rho_air(nzb)                    !< in  kg/m^3 
1430 
1431
1432                      ENDIF
1433
1434                   ENDIF
1435
1436                ENDIF
1437
1438             ENDDO
1439         
1440!
1441!--          LSM surfaces
1442             DO m = 1, surf_lsm_h%ns
1443
1444                i = surf_lsm_h%i(m)
1445                j = surf_lsm_h%j(m)
1446                k = surf_lsm_h%k(m)
1447
1448                IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
1449
1450                   !
1451                   !-- PMs
1452                   IF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM1"  &
1453                        .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM25"  &
1454                          .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM10" )  THEN
1455   
1456                      !
1457                      !--         kg/(m^2*s) * kg/m^3                           kg/(m^2*s)           
1458                      surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &
1459                                                                     rho_air(k)                                !< in kg/m^3
1460 
1461                   ELSE
1462
1463                      !
1464                      !-- VOCs
1465                      IF ( len_index_voc > 0 .AND. emt_att%species_name(match_spec_input(ispec)) == "VOC" )  THEN
1466                         !
1467                         !--          (ppm/s) * m * kg/m^3                        mole/(m^2/s)   
1468                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &     
1469                                                                        conv_to_ratio(k,j,i) *           &     !< in m^3/mole
1470                                                                         ratio2ppm *                     &     !< in ppm
1471                                                                          rho_air(k)                           !< in kg/m^3 
1472
1473                      !
1474                      !-- Other species
1475                      ELSE
1476                         !
1477                         !--           (ppm/s) * m  * kg/m^3                    kg/(m^2/s)
1478                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  & 
1479                                                                        ( 1.0_wp / emt_att%xm(ispec) ) * &     !< in mole/kg
1480                                                                         conv_to_ratio(k,j,i) *          &     !< in m^3/mole
1481                                                                          ratio2ppm *                    &     !< in ppm
1482                                                                           rho_air(k)                          !< in kg/m^3     
1483                                                   
1484                      ENDIF
1485
1486                   ENDIF
1487
1488                ENDIF
1489
1490             ENDDO
1491
1492!
1493!--          USM surfaces
1494             DO m = 1, surf_usm_h%ns
1495
1496                i = surf_usm_h%i(m)
1497                j = surf_usm_h%j(m)
1498                k = surf_usm_h%k(m)
1499
1500                IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
1501
1502                   !
1503                   !-- PMs
1504                   IF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM1" &
1505                        .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM25"  &
1506                          .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM10" )  THEN
1507                   
1508                      !
1509                      !--          kg/(m^2*s) *kg/m^3                             kg/(m^2*s)                     
1510                      surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)*  & 
1511                                                                     rho_air(k)                                !< in kg/m^3
1512
1513 
1514                   ELSE
1515                     
1516                      !
1517                      !-- VOCs
1518                      IF ( len_index_voc > 0 .AND. emt_att%species_name(match_spec_input(ispec)) == "VOC" ) THEN
1519                         !
1520                         !--          (ppm/s) * m * kg/m^3                        mole/(m^2/s)   
1521                         surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &
1522                                                                        conv_to_ratio(k,j,i) *           &     !< in m^3/mole
1523                                                                         ratio2ppm *                     &     !< in ppm
1524                                                                          rho_air(k)                           !< in kg/m^3   
1525
1526                      !
1527                      !-- Other species
1528                      ELSE
1529
1530                         !
1531                         !--          (ppm/s) * m * kg/m^3                        kg/(m^2/s)                     
1532                         surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &   
1533                                                                        ( 1.0_wp / emt_att%xm(ispec) ) * &     !< in mole/kg
1534                                                                         conv_to_ratio(k,j,i) *          &     !< in m^3/mole
1535                                                                          ratio2ppm*                     &     !< in ppm
1536                                                                           rho_air(k)                          !< in kg/m^3   
1537
1538
1539                      ENDIF
1540
1541                   ENDIF
1542
1543                ENDIF
1544 
1545             ENDDO
1546
1547          ENDDO
1548
1549       ENDIF
1550
1551       !
1552       !-- Deallocate time_factor in case of DEFAULT mode)
1553       IF ( ALLOCATED ( time_factor ) ) DEALLOCATE( time_factor )
1554
1555   ENDIF
1556
1557 END SUBROUTINE chem_emissions_setup
1558
1559 END MODULE chem_emissions_mod
Note: See TracBrowser for help on using the repository browser.