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

Last change on this file since 4391 was 4356, checked in by suehring, 4 years ago

Bugfix in message calls for local checks; error messages in init_grid slightly revised; bugfix in time_integration (uninitialized emission time index); read_restart_data (change from J.Resler): use dynamic array allocation rather than automatic arrays to avoid problems with stack memory

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