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

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

Bugfixes in indoor model and in chemical emissions

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