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

Last change on this file since 4227 was 4227, checked in by gronemeier, 22 months ago

implement new palm_date_time_mod; replaced namelist parameters time_utc_init and day_of_year_init by origin_date_time

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