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

Last change on this file since 4230 was 4230, checked in by suehring, 5 years ago

Several bugfixes: profile initialization of chemical species in restart runs; Runge-Kutta tendency array not initialized in chemistry model in restart runs; fixed determination of time indices for chemical emissions (introduced with commit -4227); Update chemistry profiles in offline nesting; initialize canopy resistances for greened building walls (even if green fraction is zero)

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