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

Last change on this file since 4441 was 4403, checked in by banzhafs, 4 years ago

chemistry model: implemented on-demand emission read mode for LOD 2

  • Property svn:keywords set to Id
File size: 96.2 KB
Line 
1!> @file chem_emissions_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2018-2019 Leibniz Universitaet Hannover
18! Copyright 2018-2019 Freie Universitaet Berlin
19! Copyright 2018-2019 Karlsruhe Institute of Technology
20!--------------------------------------------------------------------------------!
21!
22! Current revisions:
23! ------------------
24!
25!
26! Former revisions:
27! -----------------
28! $Id: chem_emissions_mod.f90 4403 2020-02-12 13:08:46Z suehring $
29! Implemented on-demand read mode for LOD 2 (ECC)
30!  - added following module global variables
31!    - input_file_chem (namesake in netcdf_data_input_mod is local)
32!    - timestamps
33!  - added following public subroutines / interfaces
34!    - chem_emisisons_header_init
35!    - chem_emisisons_update_on_demand
36!  - added following local subroutines
37!    - chem_emisisons_header_init_lod2
38!    - chem_emisisons_update_on_demand_lod2
39!  - added following local auxiliary subroutines
40!    - chem_emissions_init_species ( )
41!    - chem_emissions_init_timestamps ( )
42!    - chem_emissions_assign_surface_flux ( )
43!  - added following local functions
44!    - chem_emisisons_convert_base_units ( )
45!    - chem_emissions_mass_2_molar_flux ( )
46!    - chem_emissions_locate_species ( )
47!    - chem_emissions_locate_timestep ( )
48!  - added following error messages
49!    - CM0468 - LOD mismatch (namelist / chemistry file)
50!    - CM0469 - Timestamps no in choronological order
51!  - depreciated unused module variable filename_emis
52!
53! 4356 2019-12-20 17:09:33Z suehring
54! Minor formatting adjustment
55!
56! 4242 2019-09-27 12:59:10Z suehring
57! Adjust index_hh access to new definition accompanied with new
58! palm_date_time_mod. Note, this is just a preleminary fix. (E Chan)
59!
60! 4230 2019-09-11 13:58:14Z suehring
61! Bugfix, consider that time_since_reference_point can be also negative when
62! time indices are determined.
63!
64! 4227 2019-09-10 18:04:34Z gronemeier
65! implement new palm_date_time_mod
66!
67! 4223 2019-09-10 09:20:47Z gronemeier
68! Unused routine chem_emissions_check_parameters commented out due to uninitialized content
69!
70! 4182 2019-08-22 15:20:23Z scharf
71! Corrected "Former revisions" section
72!
73! 4157 2019-08-14 09:19:12Z suehring
74! Replace global arrays also in mode_emis branch
75!
76! 4154 2019-08-13 13:35:59Z suehring
77! Replace global arrays for emissions by local ones.
78!
79! 4144 2019-08-06 09:11:47Z raasch
80! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
81!
82! 4055 2019-06-27 09:47:29Z suehring
83! - replaced local thermo. constants w/ module definitions in
84!   basic_constants_and_equations_mod (rgas_univ, p_0, r_d_cp)
85! - initialize array emis_distribution immediately following allocation
86! - lots of minor formatting changes based on review sesson in 20190325
87!   (E.C. Chan)
88!
89! 3968 2019-05-13 11:04:01Z suehring
90! - in subroutine chem_emissions_match replace all decision structures relating to
91!   mode_emis to emiss_lod
92! - in subroutine chem_check_parameters replace emt%nspec with emt%n_emiss_species
93! - spring cleaning (E.C. Chan)
94!
95! 3885 2019-04-11 11:29:34Z kanani
96! Changes related to global restructuring of location messages and introduction
97! of additional debug messages
98!
99! 3831 2019-03-28 09:11:22Z forkel
100! added nvar to USE chem_gasphase_mod (chem_modules will not include nvar anymore)
101!
102! 3788 2019-03-07 11:40:09Z banzhafs
103! Removed unused variables from chem_emissions_mod
104!
105! 3772 2019-02-28 15:51:57Z suehring
106! - In case of parametrized emissions, assure that emissions are only on natural
107!   surfaces (i.e. streets) and not on urban surfaces.
108! - some unnecessary if clauses removed
109!
110! 3685 2019 -01-21 01:02:11Z knoop
111! Some interface calls moved to module_interface + cleanup
112! 3286 2018-09-28 07:41:39Z forkel
113!
114! Authors:
115! --------
116! @author Emmanuele Russo (FU-Berlin)
117! @author Sabine Banzhaf  (FU-Berlin)
118! @author Martijn Schaap  (FU-Berlin, TNO Utrecht)
119!
120! Description:
121! ------------
122!>  MODULE for reading-in Chemistry Emissions data
123!>
124!> @todo Rename nspec to n_emis to avoid inteferece with nspec from chem_gasphase_mod
125!> @todo Check_parameters routine should be developed: for now it includes just one condition
126!> @todo Use of Restart files not contempled at the moment
127!> @todo revise indices of files read from the netcdf: preproc_emission_data and expert_emission_data
128!> @todo for now emission data may be passed on a singular vertical level: need to be more flexible
129!> @todo fill/activate restart option in chem_emissions_init
130!> @todo discuss dt_emis
131!> @note <Enter notes on the module>
132!> @bug  <Enter known bugs here>
133!------------------------------------------------------------------------------!
134
135 MODULE chem_emissions_mod
136
137    USE arrays_3d,                                                          &
138        ONLY:  rho_air
139
140    USE basic_constants_and_equations_mod,                                  &
141        ONLY:  rgas_univ, p_0, rd_d_cp
142
143    USE control_parameters,                                                 &
144        ONLY:  debug_output,                                                &
145               end_time, message_string, initializing_actions,              &
146               intermediate_timestep_count, dt_3d
147 
148    USE indices
149
150    USE kinds
151
152#if defined( __netcdf )
153    USE netcdf
154#endif
155
156    USE netcdf_data_input_mod,                                               &
157        ONLY: chem_emis_att_type, chem_emis_val_type
158
159    USE chem_gasphase_mod,                                                   &
160        ONLY: nvar, spc_names
161 
162    USE chem_modules
163
164    USE statistics,                                                          &   
165        ONLY:  weight_pres
166
167!
168!-- 20200203 (ECC)
169!-- Added new palm_date_time_mod for on-demand emission reading
170
171    USE palm_date_time_mod,                                                  &
172        ONLY:  get_date_time
173
174    IMPLICIT NONE
175
176!
177!-- Declare all global variables within the module
178
179!
180!-- 20200203 (ECC) - variable unused
181!    CHARACTER (LEN=80) ::  filename_emis             !< Variable for the name of the netcdf input file
182
183!
184!-- 20200203 (ECC) new variables for on-demand read mode
185
186    CHARACTER(LEN=512), ALLOCATABLE, DIMENSION(:) :: timestamps   !< timestamps in chemistry file
187   
188    CHARACTER(LEN=*),   PARAMETER   :: input_file_chem = 'PIDS_CHEM' !< chemistry file
189
190    INTEGER(iwp) ::  dt_emis         !< Time Step Emissions
191    INTEGER(iwp) ::  i               !< index 1st selected dimension (some dims are not spatial)
192    INTEGER(iwp) ::  j               !< index 2nd selected dimension
193    INTEGER(iwp) ::  i_start         !< Index to start read variable from netcdf along one dims
194    INTEGER(iwp) ::  i_end           !< Index to end read variable from netcdf in one dims
195    INTEGER(iwp) ::  j_start         !< Index to start read variable from netcdf in additional dims
196    INTEGER(iwp) ::  j_end           !< Index to end read variable from netcdf in additional dims
197    INTEGER(iwp) ::  len_index       !< length of index (used for several indices)
198    INTEGER(iwp) ::  len_index_pm    !< length of PMs index
199    INTEGER(iwp) ::  len_index_voc   !< length of voc index
200    INTEGER(iwp) ::  previous_timestamp_index !< index for current timestamp (20200203 ECC)
201    INTEGER(iwp) ::  z_start         !< Index to start read variable from netcdf in additional dims
202    INTEGER(iwp) ::  z_end           !< Index to end read variable from netcdf in additional dims
203
204    REAL(wp) ::  conversion_factor                   !< Units Conversion Factor
205
206    SAVE
207
208! !
209! !-- Checks Input parameters
210!     INTERFACE chem_emissions_check_parameters
211!        MODULE PROCEDURE chem_emissions_check_parameters
212!     END INTERFACE chem_emissions_check_parameters
213!
214!-- Matching Emissions actions 
215    INTERFACE chem_emissions_match
216       MODULE PROCEDURE chem_emissions_match
217    END INTERFACE chem_emissions_match
218!
219!-- Initialization actions 
220    INTERFACE chem_emissions_init
221       MODULE PROCEDURE chem_emissions_init
222    END INTERFACE chem_emissions_init
223!
224!-- Setup of Emissions
225    INTERFACE chem_emissions_setup
226       MODULE PROCEDURE chem_emissions_setup
227    END INTERFACE chem_emissions_setup
228
229!
230!-- 20200203 (ECC) new interfaces for on-demand mode
231
232!
233!-- initialization actions for on-demand mode
234    INTERFACE chem_emissions_header_init
235       MODULE PROCEDURE chem_emissions_header_init
236    END INTERFACE chem_emissions_header_init
237!
238!-- load emission data for on-demand mode
239    INTERFACE chem_emissions_update_on_demand
240       MODULE PROCEDURE chem_emissions_update_on_demand
241    END INTERFACE chem_emissions_update_on_demand
242
243!
244!-- 20200203 (ECC) update public routines
245
246!    PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup
247
248    PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup,    &
249           chem_emissions_header_init, chem_emissions_update_on_demand
250!
251!-- Public Variables
252    PUBLIC conversion_factor, len_index, len_index_pm, len_index_voc
253
254 CONTAINS
255
256! !------------------------------------------------------------------------------!
257! ! Description:
258! ! ------------
259! !> Routine for checking input parameters
260! !------------------------------------------------------------------------------!
261!  SUBROUTINE chem_emissions_check_parameters
262!
263!     IMPLICIT NONE
264!
265!     TYPE(chem_emis_att_type) ::  emt
266!
267! !
268! !-- Check if species count matches the number of names
269! !-- passed for the chemiscal species
270!
271!     IF  ( SIZE(emt%species_name) /= emt%n_emiss_species  )  THEN
272! !    IF  ( SIZE(emt%species_name) /= emt%nspec  )  THEN
273!
274!        message_string = 'Numbers of input emission species names and number of species'     //      &
275!                          'for which emission values are given do not match'
276!        CALL message( 'chem_emissions_check_parameters', 'CM0437', 2, 2, 0, 6, 0 )
277!
278!     ENDIF
279!
280!  END SUBROUTINE chem_emissions_check_parameters
281
282
283!------------------------------------------------------------------------------!
284! Description:
285! ------------
286!> Matching the chemical species indices. The routine checks what are the
287!> indices of the emission input species and the corresponding ones of the
288!> model species. The routine gives as output a vector containing the number
289!> of common species: it is important to note that while the model species
290!> are distinct, their values could be given to a single species in input.
291!> For example, in the case of NO2 and NO, values may be passed in input as
292!> NOX values.
293!------------------------------------------------------------------------------!
294
295SUBROUTINE chem_emissions_match( emt_att,len_index )   
296
297    INTEGER(iwp)  ::  ind_inp                    !< Parameters for cycling through chemical input species
298    INTEGER(iwp)  ::  ind_mod                    !< Parameters for cycling through chemical model species
299    INTEGER(iwp)  ::  ind_voc                    !< Indices to check whether a split for voc should be done
300    INTEGER(iwp)  ::  ispec                      !< index for cycle over effective number of emission species
301    INTEGER(iwp)  ::  nspec_emis_inp             !< Variable where to store # of emission species in input
302
303    INTEGER(iwp), INTENT(INOUT)  ::  len_index   !< number of common species between input dataset & model species
304
305    TYPE(chem_emis_att_type), INTENT(INOUT) ::  emt_att     !< Chemistry Emission Array (decl. netcdf_data_input.f90)
306
307
308    IF  ( debug_output )  CALL debug_message( 'chem_emissions_match', 'start' )
309
310!
311!-- Number of input emission species
312
313    nspec_emis_inp = emt_att%n_emiss_species
314!    nspec_emis_inp=emt_att%nspec
315
316!
317!-- Check the emission LOD: 0 (PARAMETERIZED), 1 (DEFAULT), 2 (PRE-PROCESSED)
318!
319    SELECT CASE (emiss_lod)
320
321!
322!-- LOD 0 (PARAMETERIZED mode)
323
324       CASE (0)
325
326          len_index = 0
327
328! number of species and number of matched species can be different
329! but call is only made if both are greater than zero
330
331          IF  ( nvar > 0  .AND.  nspec_emis_inp > 0 )  THEN
332
333!
334!-- Cycle over model species
335
336             DO  ind_mod = 1, nvar
337                ind_inp = 1
338                DO  WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )    !< 'novalue' is the default 
339                   IF  ( TRIM( surface_csflux_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
340                      len_index = len_index + 1
341                   ENDIF
342                   ind_inp = ind_inp + 1
343                ENDDO
344             ENDDO
345
346             IF  ( len_index > 0 )  THEN
347
348!
349!-- Allocation of Arrays of the matched species
350
351                ALLOCATE ( match_spec_input(len_index) ) 
352                ALLOCATE ( match_spec_model(len_index) )
353
354!
355!-- Pass species indices to declared arrays
356
357                len_index = 0
358
359                DO  ind_mod = 1, nvar 
360                   ind_inp = 1
361                   DO  WHILE  ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )
362                      IF  ( TRIM(surface_csflux_name(ind_inp)) ==          &
363                            TRIM(spc_names(ind_mod))            )  THEN
364                         len_index = len_index + 1
365                         match_spec_input(len_index) = ind_inp
366                         match_spec_model(len_index) = ind_mod
367                      ENDIF
368                   ind_inp = ind_inp + 1
369                   END DO
370                END DO
371
372!
373!-- Check
374
375                DO  ispec = 1, len_index
376
377                   IF  ( emiss_factor_main(match_spec_input(ispec) ) < 0    .AND.     &
378                         emiss_factor_side(match_spec_input(ispec) ) < 0 )  THEN
379
380                      message_string = 'PARAMETERIZED emissions mode selected:'             //          &
381                                       ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES'        //          &
382                                       ' but values of scaling factors for street types'    //          &
383                                       ' emiss_factor_main AND emiss_factor_side'           //          &
384                                       ' not provided for each of the emissions species'    //          &
385                                       ' or not provided at all: PLEASE set a finite value' //          &
386                                       ' for these parameters in the chemistry namelist'         
387                      CALL message( 'chem_emissions_matching', 'CM0442', 2, 2, 0, 6, 0 )
388 
389                   ENDIF
390
391                END DO
392
393
394             ELSE
395               
396                message_string = 'Non of given Emission Species'           //           &
397                                 ' matches'                                //           &
398                                 ' model chemical species'                 //           &
399                                 ' Emission routine is not called'         
400                CALL message( 'chem_emissions_matching', 'CM0443', 0, 0, 0, 6, 0 )
401 
402             ENDIF
403
404          ELSE
405     
406             message_string = 'Array of Emission species not allocated: '             //          &
407                              ' Either no emission species are provided as input or'  //          &
408                              ' no chemical species are used by PALM.'                //          &
409                              ' Emission routine is not called'                                   
410             CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 ) 
411 
412          ENDIF
413
414!
415!-- LOD 1 (DEFAULT mode)
416
417       CASE (1)
418
419          len_index = 0          ! total number of species (to be accumulated) 
420          len_index_voc = 0      ! total number of VOCs (to be accumulated)
421          len_index_pm = 3       ! total number of PMs: PM1, PM2.5, PM10.
422
423!
424!-- number of model species and input species could be different
425!-- but process this only when both are non-zero
426 
427          IF  ( nvar > 0  .AND.  nspec_emis_inp > 0 )  THEN
428
429!
430!-- Cycle over model species
431             DO  ind_mod = 1, nvar
432
433!
434!-- Cycle over input species
435
436                DO  ind_inp = 1, nspec_emis_inp
437
438!
439!-- Check for VOC Species
440
441                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN
442                      DO  ind_voc= 1, emt_att%nvoc
443                       
444                         IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
445                            len_index = len_index + 1
446                            len_index_voc = len_index_voc + 1
447                         ENDIF
448                         
449                      END DO
450                   ENDIF
451
452!
453!-- PMs: There is one input species name for all PM
454!-- This variable has 3 dimensions, one for PM1, PM2.5 and PM10
455
456                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "PM" )  THEN
457
458                      IF      ( TRIM( spc_names(ind_mod) ) == "PM1" )   THEN
459                         len_index = len_index + 1
460                      ELSEIF  ( TRIM( spc_names(ind_mod) ) == "PM25" )  THEN
461                         len_index = len_index + 1
462                      ELSEIF  ( TRIM( spc_names(ind_mod) ) == "PM10" )  THEN
463                         len_index = len_index + 1
464                      ENDIF
465
466                   ENDIF
467
468!
469!-- NOX: NO2 and NO
470
471                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "NOX" )  THEN
472
473                      IF     ( TRIM( spc_names(ind_mod) ) == "NO"  )  THEN
474                         len_index = len_index + 1
475                      ELSEIF  ( TRIM( spc_names(ind_mod) ) == "NO2" )  THEN
476                         len_index = len_index + 1
477                      ENDIF
478
479                   ENDIF
480
481!
482!-- SOX: SO2 and SO4
483
484                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "SOX" )  THEN
485
486                      IF      ( TRIM( spc_names(ind_mod) ) == "SO2" )  THEN
487                         len_index = len_index + 1
488                      ELSEIF  ( TRIM( spc_names(ind_mod) ) == "SO4" )  THEN
489                         len_index = len_index + 1
490                      ENDIF
491
492                   ENDIF
493
494!
495!-- Other Species
496
497                   IF  ( TRIM( emt_att%species_name(ind_inp) ) ==             &
498                        TRIM( spc_names(ind_mod) ) )  THEN
499                      len_index = len_index + 1
500                   ENDIF
501
502                END DO  ! ind_inp ...
503
504             END DO  ! ind_mod ...
505
506
507!
508!-- Allocate arrays
509
510             IF  ( len_index > 0 )  THEN
511
512                ALLOCATE ( match_spec_input(len_index) )
513                ALLOCATE ( match_spec_model(len_index) )
514
515                IF  ( len_index_voc > 0 )  THEN
516
517!
518!-- Contains indices of the VOC model species
519
520                   ALLOCATE( match_spec_voc_model(len_index_voc) )
521
522!
523!-- Contains the indices of different values of VOC composition
524!-- of input variable VOC_composition
525
526                   ALLOCATE( match_spec_voc_input(len_index_voc) )
527
528                ENDIF
529
530!
531!-- Pass the species indices to declared arrays
532
533                len_index = 0
534                len_index_voc = 0
535               
536                DO  ind_mod = 1, nvar
537                   DO  ind_inp = 1, nspec_emis_inp 
538
539!
540!-- VOCs
541
542                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC"  .AND.        &
543                            ALLOCATED (match_spec_voc_input) )  THEN
544
545                         DO  ind_voc = 1, emt_att%nvoc
546
547                            IF  ( TRIM( emt_att%voc_name(ind_voc) ) ==                  &
548                                 TRIM( spc_names(ind_mod) ) )  THEN
549
550                               len_index     = len_index + 1
551                               len_index_voc = len_index_voc + 1
552                       
553                               match_spec_input(len_index) = ind_inp
554                               match_spec_model(len_index) = ind_mod
555
556                               match_spec_voc_input(len_index_voc) = ind_voc
557                               match_spec_voc_model(len_index_voc) = ind_mod
558
559                            ENDIF
560
561                         END DO
562
563                      ENDIF
564
565!
566!-- PMs
567
568                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == "PM" )  THEN
569
570                         IF      ( TRIM( spc_names(ind_mod) ) == "PM1"  )  THEN
571                            len_index = len_index + 1
572                            match_spec_input(len_index) = ind_inp
573                            match_spec_model(len_index) = ind_mod
574                         ELSEIF  ( TRIM( spc_names(ind_mod) ) == "PM25" )  THEN
575                            len_index = len_index + 1
576                            match_spec_input(len_index) = ind_inp
577                            match_spec_model(len_index) = ind_mod
578                         ELSEIF  ( TRIM( spc_names(ind_mod) ) == "PM10" )  THEN
579                            len_index = len_index + 1
580                            match_spec_input(len_index) = ind_inp
581                            match_spec_model(len_index) = ind_mod
582                         ENDIF
583
584                      ENDIF
585
586!
587!-- NOX
588                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == "NOX" )  THEN
589
590                         IF      ( TRIM( spc_names(ind_mod) ) == "NO"  )  THEN
591                            len_index = len_index + 1
592
593                            match_spec_input(len_index) = ind_inp
594                            match_spec_model(len_index) = ind_mod
595
596                         ELSEIF  ( TRIM( spc_names(ind_mod) ) == "NO2" )  THEN
597                            len_index = len_index + 1
598
599                            match_spec_input(len_index) = ind_inp
600                            match_spec_model(len_index) = ind_mod
601 
602                         ENDIF
603
604                      ENDIF
605
606
607!
608!-- SOX
609
610                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == "SOX" ) THEN
611
612                         IF  ( TRIM( spc_names(ind_mod) ) == "SO2" )  THEN
613                            len_index = len_index + 1
614                            match_spec_input(len_index) = ind_inp
615                            match_spec_model(len_index) = ind_mod
616                         ELSEIF  ( TRIM( spc_names(ind_mod) ) == "SO4" )  THEN
617                            len_index = len_index + 1
618                            match_spec_input(len_index) = ind_inp
619                            match_spec_model(len_index) = ind_mod
620                         ENDIF
621
622                      ENDIF
623
624!
625!-- Other Species
626
627                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
628                         len_index = len_index + 1
629                         match_spec_input(len_index) = ind_inp
630                         match_spec_model(len_index) = ind_mod
631                      ENDIF
632
633                   END DO  ! inp_ind
634
635                END DO     ! inp_mod
636
637!
638!-- Error reporting (no matching)
639
640             ELSE
641
642                message_string = 'None of given Emission Species matches'           //   &
643                                 ' model chemical species'                          //   &
644                                 ' Emission routine is not called'         
645                CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 ) 
646
647             ENDIF
648
649!
650!-- Error reporting (no species)
651
652          ELSE
653
654             message_string = 'Array of Emission species not allocated: '             // &
655                              ' Either no emission species are provided as input or'  // &
656                              ' no chemical species are used by PALM:'                // &
657                              ' Emission routine is not called'                                   
658             CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 ) 
659 
660          ENDIF
661
662!
663!-- LOD 2 (PRE-PROCESSED mode)
664
665       CASE (2)
666
667          len_index = 0
668          len_index_voc = 0
669
670          IF  ( nvar > 0  .AND.  nspec_emis_inp > 0 )  THEN
671!
672!-- Cycle over model species
673             DO  ind_mod = 1, nvar 
674
675!
676!-- Cycle over input species 
677                DO  ind_inp = 1, nspec_emis_inp
678
679!
680!-- Check for VOC Species
681
682                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN       
683                      DO  ind_voc = 1, emt_att%nvoc
684                         IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
685                            len_index     = len_index + 1
686                            len_index_voc = len_index_voc + 1
687                         ENDIF
688                      END DO
689                   ENDIF
690
691!
692!-- Other Species
693
694                   IF  ( TRIM(emt_att%species_name(ind_inp)) == TRIM(spc_names(ind_mod)) )  THEN
695                      len_index = len_index + 1
696                   ENDIF
697                ENDDO
698             ENDDO
699
700!
701!-- Allocate array for storing the indices of the matched species
702
703             IF  ( len_index > 0 )  THEN
704 
705                ALLOCATE ( match_spec_input(len_index) ) 
706 
707                ALLOCATE ( match_spec_model(len_index) )
708
709                IF  ( len_index_voc > 0 )  THEN
710!
711!-- contains indices of the VOC model species
712                   ALLOCATE( match_spec_voc_model(len_index_voc) )
713!
714!-- contains the indices of different values of VOC composition of input variable VOC_composition
715                   ALLOCATE( match_spec_voc_input(len_index_voc) )
716
717                ENDIF
718
719!
720!-- pass the species indices to declared arrays
721
722                len_index = 0
723
724!
725!-- Cycle over model species
726
727                DO  ind_mod = 1, nvar
728 
729!
730!-- Cycle over Input species 
731
732                   DO  ind_inp = 1, nspec_emis_inp
733
734!
735!-- VOCs
736
737                      IF  ( TRIM(emt_att%species_name(ind_inp) ) == "VOC"  .AND.     &
738                           ALLOCATED(match_spec_voc_input) )  THEN
739                         
740                         DO  ind_voc= 1, emt_att%nvoc
741                            IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
742                               len_index = len_index + 1
743                               len_index_voc = len_index_voc + 1
744                       
745                               match_spec_input(len_index) = ind_inp
746                               match_spec_model(len_index) = ind_mod
747
748                               match_spec_voc_input(len_index_voc) = ind_voc
749                               match_spec_voc_model(len_index_voc) = ind_mod                         
750                            ENDIF
751                         END DO
752                      ENDIF
753
754!
755!-- Other Species
756
757                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
758                         len_index = len_index + 1
759                         match_spec_input(len_index) = ind_inp
760                         match_spec_model(len_index) = ind_mod
761                      ENDIF
762
763                   END DO  ! ind_inp
764                END DO     ! ind_mod
765
766             ELSE  ! if len_index_voc <= 0
767
768!
769!-- in case there are no species matching (just informational message)
770
771                message_string = 'Non of given emission species'            //         &
772                                 ' matches'                                //          &
773                                 ' model chemical species:'                //          &
774                                 ' Emission routine is not called' 
775                CALL message( 'chem_emissions_matching', 'CM0438', 0, 0, 0, 6, 0 )
776             ENDIF
777
778!
779!-- Error check (no matching)
780 
781          ELSE
782
783!
784!-- either spc_names is zero or nspec_emis_inp is not allocated
785             message_string = 'Array of Emission species not allocated:'              // &
786                              ' Either no emission species are provided as input or'  // &
787                              ' no chemical species are used by PALM:'                // &
788                              ' Emission routine is not called'                 
789             CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 ) 
790
791          ENDIF 
792
793!
794!-- If emission module is switched on but mode_emis is not specified or it is given the wrong name
795
796!
797!-- Error check (no species)
798
799       CASE DEFAULT
800
801          message_string = 'Emission Module switched ON, but'                           //         &
802                           ' either no emission mode specified or incorrectly given :'  //         &
803                           ' please, pass the correct value to the namelist parameter "mode_emis"'                 
804          CALL message( 'chem_emissions_matching', 'CM0445', 2, 2, 0, 6, 0 )             
805
806       END SELECT
807
808       IF  ( debug_output )  CALL debug_message( 'chem_emissions_match', 'end' )
809
810 END SUBROUTINE chem_emissions_match
811
812 
813!------------------------------------------------------------------------------!
814! Description:
815! ------------
816!> Initialization:
817!> Netcdf reading, arrays allocation and first calculation of cssws
818!> fluxes at timestep 0
819!------------------------------------------------------------------------------!
820
821 SUBROUTINE chem_emissions_init
822
823    USE netcdf_data_input_mod,                                              &
824        ONLY:  chem_emis, chem_emis_att
825   
826    IMPLICIT NONE
827 
828    INTEGER(iwp) :: ispec                        !< running index
829
830!   
831!-- Actions for initial runs
832!  IF  (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
833!--    ...
834!   
835!
836!-- Actions for restart runs
837!  ELSE
838!--    ...
839!
840!  ENDIF
841
842
843    IF  ( debug_output )  CALL debug_message( 'chem_emissions_init', 'start' )
844
845!
846!-- Matching
847
848    CALL  chem_emissions_match( chem_emis_att, n_matched_vars )
849 
850    IF  ( n_matched_vars == 0 )  THEN
851 
852       emission_output_required = .FALSE.
853
854    ELSE
855
856       emission_output_required = .TRUE.
857
858
859!
860!-- Set molecule masses  (in kg/mol)
861
862       ALLOCATE( chem_emis_att%xm(n_matched_vars) )
863
864       DO  ispec = 1, n_matched_vars
865          SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
866             CASE ( 'SO2'  );  chem_emis_att%xm(ispec) = xm_S + xm_O * 2
867             CASE ( 'SO4'  );  chem_emis_att%xm(ispec) = xm_S + xm_O * 4
868             CASE ( 'NO'   );  chem_emis_att%xm(ispec) = xm_N + xm_O
869             CASE ( 'NO2'  );  chem_emis_att%xm(ispec) = xm_N + xm_O * 2
870             CASE ( 'NH3'  );  chem_emis_att%xm(ispec) = xm_N + xm_H * 3
871             CASE ( 'CO'   );  chem_emis_att%xm(ispec) = xm_C + xm_O
872             CASE ( 'CO2'  );  chem_emis_att%xm(ispec) = xm_C + xm_O * 2
873             CASE ( 'CH4'  );  chem_emis_att%xm(ispec) = xm_C + xm_H * 4
874             CASE ( 'HNO3' );  chem_emis_att%xm(ispec) = xm_H + xm_N + xm_O*3
875             CASE DEFAULT
876                chem_emis_att%xm(ispec) = 1.0_wp
877          END SELECT
878       ENDDO
879
880   
881!
882!-- Get emissions for the first time step base on LOD (if defined)
883!-- or emission mode (if no LOD defined)
884
885!
886!-- NOTE - I could use a combined if ( lod = xxx .or. mode = 'XXX' )
887!--        type of decision structure but I think it is much better
888!--        to implement it this way (i.e., conditional on lod if it
889!--        is defined, and mode if not) as we can easily take out
890!--        the case structure for mode_emis later on.
891
892       IF   ( emiss_lod < 0 )  THEN  !-- no LOD defined (not likely)
893
894          SELECT CASE ( TRIM( mode_emis ) )   
895
896             CASE ( 'PARAMETERIZED' )     ! LOD 0
897
898                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
899                   ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
900                ENDIF
901
902                CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars)
903
904             CASE ( 'DEFAULT' )           ! LOD 1
905
906                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
907                   ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
908                ENDIF
909
910                CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
911
912             CASE ( 'PRE-PROCESSED' )     ! LOD 2
913
914                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
915!
916!--                Note, at the moment emissions are considered only by surface fluxes
917!--                rather than by volume sources. Therefore, no vertical dimension is
918!--                required and is thus allocated with 1. Later when volume sources
919!--                are considered, the vertical dimension will increase.                 
920                   !ALLOCATE( emis_distribution(nzb:nzt+1,nys:nyn,nxl:nxr,n_matched_vars) )
921                   ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
922                ENDIF
923 
924                CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
925
926          END SELECT
927
928       ELSE  ! if LOD is defined
929
930          SELECT CASE ( emiss_lod )
931
932             CASE ( 0 )     ! parameterized mode
933
934                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
935                   ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
936                ENDIF
937
938                CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars)
939
940             CASE ( 1 )     ! default mode
941
942                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
943                   ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
944                ENDIF
945
946                CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
947
948             CASE ( 2 )     ! pre-processed mode
949
950                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
951                   ALLOCATE( emis_distribution(nzb:nzt+1,nys:nyn,nxl:nxr,n_matched_vars) )
952                ENDIF
953 
954                CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
955
956          END SELECT
957
958       ENDIF
959
960!
961! -- initialize
962
963       emis_distribution = 0.0_wp
964
965    ENDIF
966
967    IF  ( debug_output )  CALL debug_message( 'chem_emissions_init', 'end' )
968
969 END SUBROUTINE chem_emissions_init
970
971
972
973!------------------------------------------------------------------------------!
974! Description:
975! ------------
976!> Routine for Update of Emission values at each timestep.
977!>
978!> @todo Clarify the correct usage of index_dd, index_hh and index_mm. Consider
979!>       renaming of these variables.
980!> @todo Clarify time used in emis_lod=2 mode. ATM, the used time seems strange.
981!-------------------------------------------------------------------------------!
982
983 SUBROUTINE chem_emissions_setup( emt_att, emt, n_matched_vars )
984 
985   USE surface_mod,                                               &
986       ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
987
988   USE netcdf_data_input_mod,                                     &
989       ONLY:  street_type_f
990
991   USE arrays_3d,                                                 &       
992       ONLY: hyp, pt 
993
994    USE control_parameters, &
995        ONLY: time_since_reference_point
996
997    USE palm_date_time_mod, &
998        ONLY: days_per_week, get_date_time, hours_per_day, months_per_year, seconds_per_day
999   
1000    IMPLICIT NONE
1001 
1002
1003    TYPE(chem_emis_att_type), INTENT(INOUT) ::  emt_att                         !< variable to store emission information
1004
1005    TYPE(chem_emis_val_type), INTENT(INOUT), ALLOCATABLE, DIMENSION(:) ::  emt  !< variable to store emission input values,
1006                                                                                !< depending on the considered species
1007
1008    INTEGER,INTENT(IN) ::  n_matched_vars                                       !< Output of matching routine with number
1009                                                                                !< of matched species
1010
1011    INTEGER(iwp) ::  i                                                          !< running index for grid in x-direction
1012    INTEGER(iwp) ::  i_pm_comp                                                  !< index for number of PM components
1013    INTEGER(iwp) ::  icat                                                       !< Index for number of categories
1014    INTEGER(iwp) ::  ispec                                                      !< index for number of species
1015    INTEGER(iwp) ::  ivoc                                                       !< Index for number of VOCs
1016    INTEGER(iwp) ::  j                                                          !< running index for grid in y-direction
1017    INTEGER(iwp) ::  k                                                          !< running index for grid in z-direction
1018    INTEGER(iwp) ::  m                                                          !< running index for horizontal surfaces
1019
1020    INTEGER(iwp) ::  day_of_month                                               !< day of the month
1021    INTEGER(iwp) ::  day_of_week                                                !< day of the week
1022    INTEGER(iwp) ::  day_of_year                                                !< day of the year
1023    INTEGER(iwp) ::  days_since_reference_point                                 !< days since reference point
1024    INTEGER(iwp) ::  hour_of_day                                                !< hour of the day
1025    INTEGER(iwp) ::  month_of_year                                              !< month of the year
1026    INTEGER(iwp) ::  index_dd                                                   !< index day
1027    INTEGER(iwp) ::  index_hh                                                   !< index hour
1028    INTEGER(iwp) ::  index_mm                                                   !< index month
1029
1030    REAL(wp) ::  time_utc_init                                                  !< second of day of initial date
1031
1032    !
1033    !-- CONVERSION FACTORS: TIME 
1034    REAL(wp), PARAMETER ::  hour_per_year =  8760.0_wp  !< number of hours in a year of 365 days 
1035    REAL(wp), PARAMETER ::  s_per_hour    =  3600.0_wp  !< number of sec per hour (s)/(hour)   
1036    REAL(wp), PARAMETER ::  s_per_day     = 86400.0_wp  !< number of sec per day (s)/(day) 
1037
1038    REAL(wp), PARAMETER ::  day_to_s      = 1.0_wp/s_per_day                   !< conversion day   -> sec
1039    REAL(wp), PARAMETER ::  hour_to_s     = 1.0_wp/s_per_hour                  !< conversion hours -> sec
1040    REAL(wp), PARAMETER ::  year_to_s     = 1.0_wp/(s_per_hour*hour_per_year)  !< conversion year  -> sec
1041    !
1042    !-- CONVERSION FACTORS: MASS
1043    REAL(wp), PARAMETER ::  g_to_kg       = 1.0E-03_wp     !< Conversion from g to kg (kg/g)
1044    REAL(wp), PARAMETER ::  miug_to_kg    = 1.0E-09_wp     !< Conversion from g to kg (kg/g)
1045    REAL(wp), PARAMETER ::  tons_to_kg    = 100.0_wp       !< Conversion from tons to kg (kg/tons)   
1046    !
1047    !-- CONVERSION FACTORS: PPM
1048    REAL(wp), PARAMETER ::  ratio2ppm     = 1.0E+06_wp 
1049 
1050    REAL(wp), DIMENSION(24)                        ::  par_emis_time_factor    !< time factors for the parameterized mode:
1051                                                                               !< fixed houlry profile for example day
1052    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  conv_to_ratio           !< factor used for converting input
1053                                                                               !< to concentration ratio
1054    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  tmp_temp                !< temporary variable for abs. temperature
1055
1056    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  delta_emis   !< incremental emission factor                     
1057    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  time_factor  !< factor for time scaling of emissions
1058    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  emis         !< emission factor
1059
1060    IF  ( emission_output_required )  THEN
1061
1062!
1063!-- Set emis_dt to be used - since chemistry ODEs can be stiff, the option
1064!-- to solve them at every RK substep is present to help improve stability
1065!-- should the need arises
1066 
1067       IF  ( call_chem_at_all_substeps )  THEN
1068
1069          dt_emis = dt_3d * weight_pres(intermediate_timestep_count)
1070
1071       ELSE
1072
1073          dt_emis = dt_3d
1074
1075       ENDIF
1076
1077 !
1078 !-- Conversion of units to the ones employed in PALM
1079 !-- In PARAMETERIZED mode no conversion is performed: in this case input units are fixed
1080
1081        IF  ( TRIM( mode_emis ) == "DEFAULT"  .OR.  TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
1082
1083          SELECT CASE ( TRIM( emt_att%units ) )
1084
1085             CASE ( 'kg/m2/s', 'KG/M2/S'    );     conversion_factor = 1.0_wp                   ! kg
1086             CASE ( 'kg/m2/hour', 'KG/M2/HOUR' );  conversion_factor = hour_to_s
1087             CASE ( 'kg/m2/day', 'KG/M2/DAY'  );   conversion_factor = day_to_s
1088             CASE ( 'kg/m2/year', 'KG/M2/YEAR' );  conversion_factor = year_to_s
1089
1090             CASE ( 'ton/m2/s', 'TON/M2/S'    );     conversion_factor = tons_to_kg             ! tonnes
1091             CASE ( 'ton/m2/hour', 'TON/M2/HOUR' );  conversion_factor = tons_to_kg*hour_to_s
1092             CASE ( 'ton/m2/year', 'TON/M2/YEAR' );  conversion_factor = tons_to_kg*year_to_s
1093
1094             CASE ( 'g/m2/s', 'G/M2/S'    );     conversion_factor = g_to_kg                    ! grams
1095             CASE ( 'g/m2/hour', 'G/M2/HOUR' );  conversion_factor = g_to_kg*hour_to_s
1096             CASE ( 'g/m2/year', 'G/M2/YEAR' );  conversion_factor = g_to_kg*year_to_s
1097
1098             CASE ( 'micrograms/m2/s', 'MICROGRAMS/M2/S'    );     conversion_factor = miug_to_kg            ! ug
1099             CASE ( 'micrograms/m2/hour', 'MICROGRAMS/M2/HOUR' );  conversion_factor = miug_to_kg*hour_to_s
1100             CASE ( 'micrograms/m2/year', 'MICROGRAMS/M2/YEAR' );  conversion_factor = miug_to_kg*year_to_s
1101
1102!
1103!-- Error check (need units)
1104
1105             CASE DEFAULT 
1106                message_string = 'The Units of the provided emission input'                 // &
1107                                 ' are not the ones required by PALM-4U: please check '     // &
1108                                 ' emission module documentation.'                                 
1109                CALL message( 'chem_emissions_setup', 'CM0446', 2, 2, 0, 6, 0 )
1110
1111          END SELECT
1112
1113       ENDIF
1114
1115!
1116!-- Conversion factor to convert  kg/m**2/s to ppm/s
1117
1118       DO  i = nxl, nxr
1119          DO  j = nys, nyn
1120
1121!
1122!-- Derive Temperature from Potential Temperature
1123
1124             tmp_temp(nzb:nzt+1,j,i) = pt(nzb:nzt+1,j,i) *                     &
1125                                       ( hyp(nzb:nzt+1) / p_0 )**rd_d_cp
1126 
1127!           
1128!-- We need to pass to cssws <- (ppm/s) * dz
1129!-- Input is Nmole/(m^2*s)
1130!-- To go to ppm*dz multiply the input by (m**2/N)*dz
1131!-- (m**2/N)*dz == V/N
1132!-- V/N = RT/P
1133
1134             conv_to_ratio(nzb:nzt+1,j,i) =  rgas_univ *                       &  ! J K-1 mol-1
1135                                             tmp_temp(nzb:nzt+1,j,i) /         &  ! K
1136                                             hyp(nzb:nzt+1)                       ! Pa
1137
1138! (ecc) for reference
1139!                   m**3/Nmole               (J/mol)*K^-1           K                      Pa         
1140!             conv_to_ratio(nzb:nzt+1,j,i) = ( (Rgas * tmp_temp(nzb:nzt+1,j,i)) / ((hyp(nzb:nzt+1))) ) 
1141
1142          ENDDO
1143       ENDDO
1144
1145
1146! (ecc) moved initialization immediately after allocation
1147!
1148!-- Initialize
1149
1150!       emis_distribution(:,nys:nyn,nxl:nxr,:) = 0.0_wp
1151
1152 
1153!
1154!-- LOD 2 (PRE-PROCESSED MODE)
1155
1156       IF  ( emiss_lod == 2 )  THEN
1157
1158! for reference (ecc)
1159!       IF  ( TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
1160
1161!
1162!-- Update time indices
1163
1164          CALL get_date_time( 0.0_wp, second_of_day=time_utc_init )
1165          CALL get_date_time( MAX( 0.0_wp, time_since_reference_point ),       &
1166                              hour=hour_of_day )
1167
1168          days_since_reference_point = INT( FLOOR( (                            &
1169                    time_utc_init + MAX( 0.0_wp, time_since_reference_point ) ) &
1170                                                   / seconds_per_day ) )
1171
1172          index_hh = days_since_reference_point * hours_per_day + hour_of_day
1173
1174!
1175!-- LOD 1 (DEFAULT MODE)
1176
1177        ELSEIF  ( emiss_lod == 1 )  THEN
1178
1179! for reference (ecc)
1180!       ELSEIF  ( TRIM( mode_emis ) == "DEFAULT" )  THEN
1181
1182!
1183!-- Allocate array where to store temporary emission values
1184
1185          IF  (  .NOT.  ALLOCATED(emis) ) ALLOCATE( emis(nys:nyn,nxl:nxr) )
1186
1187!
1188!-- Allocate time factor per category
1189
1190          ALLOCATE( time_factor(emt_att%ncat) )
1191
1192!
1193!-- Read-in hourly emission time factor
1194
1195          IF  ( TRIM(time_fac_type) == "HOUR" )  THEN
1196
1197!
1198!-- Update time indices
1199
1200             CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ),    &
1201                                 day_of_year=day_of_year, hour=hour_of_day )
1202             index_hh = ( day_of_year - 1_iwp ) * hour_of_day
1203
1204!
1205!-- Check if the index is less or equal to the temporal dimension of HOURLY emission files
1206
1207             IF  ( index_hh <= SIZE( emt_att%hourly_emis_time_factor(1,:) ) )  THEN
1208
1209!
1210!-- Read-in the correspondant time factor
1211
1212                time_factor(:) = emt_att%hourly_emis_time_factor(:,index_hh+1)     
1213
1214!
1215!-- Error check (time out of range)
1216
1217             ELSE
1218
1219                message_string = 'The "HOUR" time-factors in the DEFAULT mode '           // &
1220                              ' are not provided for each hour of the total simulation time'     
1221                CALL message( 'chem_emissions_setup', 'CM0448', 2, 2, 0, 6, 0 ) 
1222
1223             ENDIF
1224
1225!
1226!-- Read-in MDH emissions time factors
1227
1228          ELSEIF  ( TRIM( time_fac_type ) == "MDH" )  THEN
1229
1230!
1231!-- Update time indices
1232             CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ), &
1233                                 month=month_of_year,                       &
1234                                 day=day_of_month,                          &
1235                                 hour=hour_of_day,                          &
1236                                 day_of_week=day_of_week     )
1237             index_mm = month_of_year
1238             index_dd = months_per_year + day_of_week
1239             SELECT CASE(TRIM(daytype_mdh))
1240
1241                CASE ("workday")
1242                   index_hh = months_per_year + days_per_week + hour_of_day
1243
1244                CASE ("weekend")
1245                   index_hh = months_per_year + days_per_week + hours_per_day + hour_of_day
1246
1247                CASE ("holiday")
1248                   index_hh = months_per_year + days_per_week + 2*hours_per_day + hour_of_day
1249
1250             END SELECT
1251!
1252!-- Check if the index is less or equal to the temporal dimension of MDH emission files
1253
1254             IF  ( ( index_hh + index_dd + index_mm) <= SIZE( emt_att%mdh_emis_time_factor(1,:) ) )  THEN
1255
1256!
1257!-- Read in corresponding time factor
1258
1259                time_factor(:) = emt_att%mdh_emis_time_factor(:,index_mm) *    &
1260                                 emt_att%mdh_emis_time_factor(:,index_dd) *    &
1261                                 emt_att%mdh_emis_time_factor(:,index_hh+1)
1262
1263!
1264!-- Error check (MDH time factor not provided)
1265     
1266             ELSE
1267
1268                message_string = 'The "MDH" time-factors in the DEFAULT mode '           //  &
1269                              ' are not provided for each hour/day/month  of the total simulation time'     
1270                CALL message( 'chem_emissions_setup', 'CM0449', 2, 2, 0, 6, 0 )
1271
1272             ENDIF 
1273
1274!
1275!-- Error check (no time factor defined)
1276
1277          ELSE
1278                 
1279             message_string = 'In the DEFAULT mode the time factor'           //  &
1280                              ' has to be defined in the NAMELIST'     
1281             CALL message( 'chem_emissions_setup', 'CM0450', 2, 2, 0, 6, 0 )
1282         
1283          ENDIF
1284
1285!
1286!-- PARAMETERIZED MODE
1287
1288       ELSEIF ( emiss_lod == 0 )  THEN
1289
1290
1291! for reference (ecc)
1292!       ELSEIF  ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
1293       
1294!
1295!-- assign constant values of time factors, diurnal profile for traffic sector
1296
1297          par_emis_time_factor( : ) =  (/ 0.009, 0.004, 0.004, 0.009, 0.029, 0.039,   &
1298                                          0.056, 0.053, 0.051, 0.051, 0.052, 0.055,   &
1299                                          0.059, 0.061, 0.064, 0.067, 0.069, 0.069,   &
1300                                          0.049, 0.039, 0.039, 0.029, 0.024, 0.019 /)
1301         
1302          IF  (  .NOT.  ALLOCATED (time_factor) )  ALLOCATE (time_factor(1))
1303
1304!
1305!--       Get time-factor for specific hour
1306          CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ),              &
1307                              hour=hour_of_day )
1308
1309          index_hh = hour_of_day
1310          time_factor(1) = par_emis_time_factor(index_hh+1)
1311
1312       ENDIF  ! emiss_lod
1313
1314
1315!
1316!--  Emission distribution calculation
1317
1318!
1319!-- LOD 0 (PARAMETERIZED mode)
1320
1321       IF  ( emiss_lod == 0 )  THEN
1322
1323! for reference (ecc)
1324!       IF  ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
1325
1326          DO  ispec = 1, n_matched_vars
1327
1328!
1329!-- Units are micromoles/m**2*day (or kilograms/m**2*day for PMs)
1330
1331             emis_distribution(1,nys:nyn,nxl:nxr,ispec) =            &
1332                       surface_csflux(match_spec_input(ispec)) *     &
1333                       time_factor(1) * hour_to_s
1334
1335          ENDDO 
1336
1337         
1338!
1339!-- LOD 1 (DEFAULT mode)
1340
1341
1342       ELSEIF  ( emiss_lod == 1 )  THEN
1343
1344! for referene (ecc)
1345!       ELSEIF  ( TRIM( mode_emis ) == "DEFAULT" )  THEN
1346
1347!
1348!-- Allocate array for the emission value corresponding to a specific category and time factor
1349
1350          ALLOCATE (delta_emis(nys:nyn,nxl:nxr))
1351
1352!
1353!-- Cycle over categories
1354
1355          DO  icat = 1, emt_att%ncat
1356 
1357!
1358!-- Cycle over Species:  n_matched_vars represents the number of species
1359!-- in common between the emission input data and the chemistry mechanism used
1360
1361             DO  ispec = 1, n_matched_vars
1362
1363                emis(nys:nyn,nxl:nxr) =                    &
1364                       emt(match_spec_input(ispec))%       &
1365                           default_emission_data(icat,nys+1:nyn+1,nxl+1:nxr+1)
1366
1367!
1368!-- NO
1369
1370                IF  ( TRIM(spc_names(match_spec_model(ispec))) == "NO" )  THEN
1371               
1372                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *       &         ! kg/m2/s
1373                                                 time_factor(icat) *           &
1374                                                 emt_att%nox_comp(icat,1) *    &
1375                                                 conversion_factor * hours_per_day
1376
1377                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
1378                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +    &
1379                               delta_emis(nys:nyn,nxl:nxr)
1380!
1381!-- NO2
1382
1383                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "NO2" )  THEN
1384
1385                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *       &         ! kg/m2/s
1386                                                 time_factor(icat) *           &
1387                                                 emt_att%nox_comp(icat,2) *    &
1388                                                 conversion_factor * hours_per_day
1389
1390                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
1391                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +    &
1392                               delta_emis(nys:nyn,nxl:nxr)
1393 
1394!
1395!-- SO2
1396                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "SO2" )  THEN
1397                 
1398                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *       &         ! kg/m2/s
1399                                                 time_factor(icat) *           &
1400                                                 emt_att%sox_comp(icat,1) *    &
1401                                                 conversion_factor * hours_per_day
1402
1403                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
1404                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +    &
1405                               delta_emis(nys:nyn,nxl:nxr)
1406
1407!
1408!-- SO4
1409
1410                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "SO4" )  THEN
1411
1412                 
1413                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *       &         ! kg/m2/s
1414                                                 time_factor(icat) *           &
1415                                                 emt_att%sox_comp(icat,2) *    &
1416                                                 conversion_factor * hours_per_day
1417
1418                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
1419                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +    &
1420                               delta_emis(nys:nyn,nxl:nxr)
1421 
1422
1423!
1424!-- PM1
1425
1426                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "PM1" )  THEN
1427
1428                   DO  i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,1) )   ! cycle through components
1429
1430                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *              &     ! kg/m2/s
1431                                                    time_factor(icat) *                  &
1432                                                    emt_att%pm_comp(icat,i_pm_comp,1) *  &
1433                                                    conversion_factor * hours_per_day 
1434
1435                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                       &
1436                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +              &
1437                               delta_emis(nys:nyn,nxl:nxr)
1438
1439                   ENDDO
1440
1441!
1442!-- PM2.5
1443
1444                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "PM25" )  THEN
1445
1446                   DO  i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,2) )   ! cycle through components
1447
1448                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *              &     ! kg/m2/s
1449                                                    time_factor(icat) *                  &
1450                                                    emt_att%pm_comp(icat,i_pm_comp,2) *  &
1451                                                    conversion_factor * hours_per_day 
1452
1453                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                       &
1454                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +              &
1455                               delta_emis(nys:nyn,nxl:nxr)
1456 
1457                   ENDDO
1458
1459!
1460!-- PM10
1461
1462                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "PM10" )  THEN
1463
1464                   DO  i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,3) )   ! cycle through components
1465
1466                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *              &     ! kg/m2/s
1467                                                    time_factor(icat)     *              &
1468                                                    emt_att%pm_comp(icat,i_pm_comp,3) *  &
1469                                                    conversion_factor * hours_per_day 
1470
1471                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                       &
1472                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +              &
1473                               delta_emis(nys:nyn,nxl:nxr) 
1474
1475                   ENDDO
1476
1477!
1478!-- VOCs
1479
1480                ELSEIF  ( SIZE( match_spec_voc_input ) > 0 )  THEN
1481
1482                   DO  ivoc = 1, SIZE( match_spec_voc_input )          ! cycle through components
1483
1484                      IF  ( TRIM(spc_names(match_spec_model(ispec))) ==                &
1485                            TRIM(emt_att%voc_name(ivoc)) )  THEN   
1486
1487                         delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *           &
1488                                                       time_factor(icat) *               &
1489                                                       emt_att%voc_comp(icat,match_spec_voc_input(ivoc)) *   &
1490                                                       conversion_factor * hours_per_day
1491
1492                         emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                    &
1493                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +              &
1494                               delta_emis(nys:nyn,nxl:nxr) 
1495
1496                      ENDIF                       
1497
1498                   ENDDO
1499               
1500!
1501!-- any other species
1502
1503                ELSE
1504
1505                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *                 &
1506                                                 time_factor(icat) *                     &
1507                                                 conversion_factor * hours_per_day
1508 
1509                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                          &
1510                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +              &
1511                               delta_emis(nys:nyn,nxl:nxr) 
1512
1513                ENDIF  ! TRIM spc_names
1514               
1515                emis = 0 
1516               
1517             ENDDO
1518             
1519             delta_emis = 0 
1520         
1521          ENDDO
1522
1523!
1524!-- LOD 2 (PRE-PROCESSED mode)
1525
1526       ELSEIF  ( emiss_lod == 2 )  THEN
1527
1528! for reference (ecc)
1529!       ELSEIF  ( TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
1530
1531!
1532!-- Cycle over species: n_matched_vars represents the number of species
1533!-- in common between the emission input data and the chemistry mechanism used
1534
1535          DO  ispec = 1, n_matched_vars 
1536 
1537! (ecc)   
1538             emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                &
1539                       emt(match_spec_input(ispec))%                                     &
1540                           preproc_emission_data(index_hh+1,1,nys+1:nyn+1,nxl+1:nxr+1) * &
1541                       conversion_factor
1542
1543
1544!             emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                &
1545!                       emt(match_spec_input(ispec))%                                     &
1546!                           preproc_emission_data(index_hh,1,:,:) *   &
1547!                       conversion_factor
1548          ENDDO
1549
1550       ENDIF  ! emiss_lod
1551
1552       
1553!
1554!-- Cycle to transform x,y coordinates to the one of surface_mod and to assign emission values to cssws
1555
1556!
1557!-- LOD 0 (PARAMETERIZED mode)
1558!-- Units of inputs are micromoles/m2/s
1559
1560       IF  ( emiss_lod == 0 )  THEN
1561! for reference (ecc)
1562!       IF  ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
1563
1564          IF  (street_type_f%from_file)  THEN
1565
1566!
1567!-- Streets are lsm surfaces, hence, no usm surface treatment required.
1568!-- However, urban surface may be initialized via default initialization
1569!-- in surface_mod, e.g. at horizontal urban walls that are at k == 0
1570!-- (building is lower than the first grid point). Hence, in order to
1571!-- have only emissions at streets, set the surfaces emissions to zero
1572!-- at urban walls.
1573 
1574             IF  ( surf_usm_h%ns >=1 )  surf_usm_h%cssws = 0.0_wp
1575
1576!
1577!-- Treat land-surfaces.
1578
1579             DO  m = 1, surf_lsm_h%ns
1580
1581                i = surf_lsm_h%i(m)
1582                j = surf_lsm_h%j(m)
1583                k = surf_lsm_h%k(m)
1584
1585!
1586!-- set everything to zero then reassign according to street type
1587
1588                surf_lsm_h%cssws(:,m) = 0.0_wp
1589
1590                IF  ( street_type_f%var(j,i) >= main_street_id    .AND.                  &
1591                      street_type_f%var(j,i) < max_street_id   )  THEN
1592
1593!
1594!-- Cycle over matched species
1595
1596                   DO  ispec = 1, n_matched_vars
1597
1598!
1599!-- PMs are already in kilograms
1600
1601                      IF  ( TRIM(spc_names(match_spec_model(ispec))) == "PM1"     .OR.  &
1602                            TRIM(spc_names(match_spec_model(ispec))) == "PM25"    .OR.  &
1603                            TRIM(spc_names(match_spec_model(ispec))) == "PM10" )  THEN
1604
1605!
1606!-- kg/(m^2*s) * kg/m^3
1607                         surf_lsm_h%cssws(match_spec_model(ispec),m) =                   &
1608                                   emiss_factor_main(match_spec_input(ispec)) *          &
1609                                   emis_distribution(1,j,i,ispec) *                      &     ! kg/(m^2*s)
1610                                   rho_air(k)                                                  ! kg/m^3
1611                         
1612!
1613!-- Other Species
1614!-- Inputs are micromoles
1615
1616                      ELSE
1617
1618!   
1619!-- ppm/s *m *kg/m^3               
1620                         surf_lsm_h%cssws(match_spec_model(ispec),m) =                   &
1621                                   emiss_factor_main(match_spec_input(ispec)) *          &
1622                                   emis_distribution(1,j,i,ispec) *                      &     ! micromoles/(m^2*s)
1623                                   conv_to_ratio(k,j,i) *                                &     ! m^3/Nmole     
1624                                   rho_air(k)                                                  ! kg/m^3
1625
1626                      ENDIF
1627
1628                   ENDDO  ! ispec
1629
1630
1631                ELSEIF  ( street_type_f%var(j,i) >= side_street_id   .AND.               &
1632                          street_type_f%var(j,i) < main_street_id )  THEN
1633
1634!
1635!-- Cycle over matched species
1636
1637                   DO  ispec = 1, n_matched_vars
1638
1639!
1640!-- PMs are already in kilograms
1641
1642                      IF  ( TRIM(spc_names(match_spec_model(ispec))) == "PM1"      .OR.  &
1643                            TRIM(spc_names(match_spec_model(ispec))) == "PM25"     .OR.  &
1644                            TRIM(spc_names(match_spec_model(ispec))) == "PM10" )  THEN
1645
1646!
1647!-- kg/(m^2*s) * kg/m^3
1648                         surf_lsm_h%cssws(match_spec_model(ispec),m) =                   &
1649                                   emiss_factor_side(match_spec_input(ispec)) *          &
1650                                   emis_distribution(1,j,i,ispec) *                      &     ! kg/(m^2*s)
1651                                   rho_air(k)                                                  ! kg/m^3   
1652!
1653!-- Other species
1654!-- Inputs are micromoles
1655
1656                      ELSE
1657
1658!   
1659!-- ppm/s *m *kg/m^3
1660
1661                         surf_lsm_h%cssws(match_spec_model(ispec),m) =                   &
1662                                   emiss_factor_side(match_spec_input(ispec)) *          &
1663                                   emis_distribution(1,j,i,ispec) *                      &  ! micromoles/(m^2*s)
1664                                   conv_to_ratio(k,j,i) *                                &  ! m^3/Nmole       
1665                                   rho_air(k)                                               ! kg/m^3   
1666
1667                      ENDIF
1668
1669                   ENDDO  ! ispec
1670
1671!
1672!-- If no street type is defined, then assign zero emission to all the species
1673
1674! (ecc) moved to front (for reference)
1675!                ELSE
1676!
1677!                   surf_lsm_h%cssws(:,m) = 0.0_wp
1678
1679                ENDIF  ! street type
1680
1681             ENDDO  ! m
1682
1683          ENDIF  ! street_type_f%from_file
1684
1685
1686!
1687!-- LOD 1 (DEFAULT) and LOD 2 (PRE-PROCESSED)
1688
1689
1690       ELSE   
1691       
1692
1693          DO  ispec = 1, n_matched_vars 
1694                   
1695!
1696!--          Default surfaces
1697
1698             DO  m = 1, surf_def_h(0)%ns
1699
1700                i = surf_def_h(0)%i(m)
1701                j = surf_def_h(0)%j(m)
1702
1703                IF  ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
1704
1705!
1706!-- PMs
1707                   IF  ( TRIM(spc_names(match_spec_model(ispec))) == "PM1"     .OR.    &
1708                         TRIM(spc_names(match_spec_model(ispec))) == "PM25"    .OR.    &
1709                         TRIM(spc_names(match_spec_model(ispec))) == "PM10" )  THEN
1710               
1711                      surf_def_h(0)%cssws(match_spec_model(ispec),m) =         &     ! kg/m2/s * kg/m3
1712                               emis_distribution(1,j,i,ispec)*                 &     ! kg/m2/s
1713                               rho_air(nzb)                                          ! kg/m^3 
1714 
1715                   ELSE
1716
1717!
1718!-- VOCs
1719                      IF  ( len_index_voc > 0                                           .AND.      &
1720                            emt_att%species_name(match_spec_input(ispec)) == "VOC" )  THEN
1721
1722                         surf_def_h(0)%cssws(match_spec_model(ispec),m) =      &  ! ppm/s * m * kg/m3
1723                               emis_distribution(1,j,i,ispec) *                &  ! mole/m2/s
1724                               conv_to_ratio(nzb,j,i) *                        &  ! m^3/mole
1725                               ratio2ppm *                                     &  ! ppm
1726                               rho_air(nzb)                                       ! kg/m^3 
1727
1728
1729!
1730!-- Other species
1731
1732                      ELSE
1733
1734                         surf_def_h(0)%cssws(match_spec_model(ispec),m) =      &   ! ppm/s * m * kg/m3
1735                               emis_distribution(1,j,i,ispec) *                &   ! kg/m2/s
1736                               ( 1.0_wp / emt_att%xm(ispec) ) *                &   ! mole/kg
1737                               conv_to_ratio(nzb,j,i) *                        &   ! m^3/mole 
1738                               ratio2ppm *                                     &   ! ppm
1739                               rho_air(nzb)                                        !  kg/m^3 
1740 
1741                      ENDIF  ! VOC
1742
1743                   ENDIF  ! PM
1744
1745                ENDIF  ! emis_distribution > 0
1746
1747             ENDDO  ! m
1748         
1749!
1750!-- LSM surfaces
1751
1752
1753             DO  m = 1, surf_lsm_h%ns
1754
1755                i = surf_lsm_h%i(m)
1756                j = surf_lsm_h%j(m)
1757                k = surf_lsm_h%k(m)
1758
1759                IF  ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
1760
1761!
1762!-- PMs
1763                   IF  ( TRIM(spc_names(match_spec_model(ispec))) == "PM1"     .OR.    &
1764                         TRIM(spc_names(match_spec_model(ispec))) == "PM25"    .OR.    &
1765                         TRIM(spc_names(match_spec_model(ispec))) == "PM10" )  THEN
1766
1767                      surf_lsm_h%cssws(match_spec_model(ispec),m) =            &    ! kg/m2/s * kg/m3
1768                               emis_distribution(1,j,i,ispec) *                &    ! kg/m2/s
1769                               rho_air(k)                                           ! kg/m^3
1770 
1771                   ELSE
1772
1773!
1774!-- VOCs
1775
1776                      IF  ( len_index_voc > 0                                           .AND.      &
1777                            emt_att%species_name(match_spec_input(ispec)) == "VOC" )  THEN
1778
1779                         surf_lsm_h%cssws(match_spec_model(ispec),m) =         &   ! ppm/s * m * kg/m3
1780                               emis_distribution(1,j,i,ispec) *                &   ! mole/m2/s 
1781                               conv_to_ratio(k,j,i) *                          &   ! m^3/mole
1782                               ratio2ppm *                                     &   ! ppm
1783                               rho_air(k)                                          ! kg/m^3 
1784
1785!
1786!-- Other species
1787
1788                      ELSE
1789
1790                         surf_lsm_h%cssws(match_spec_model(ispec),m) =         &   ! ppm/s * m * kg/m3
1791                               emis_distribution(1,j,i,ispec) *                &   ! kg/m2/s
1792                               ( 1.0_wp / emt_att%xm(ispec) ) *                &   ! mole/kg
1793                               conv_to_ratio(k,j,i) *                          &   ! m^3/mole
1794                               ratio2ppm *                                     &   ! ppm
1795                               rho_air(k)                                          ! kg/m^3     
1796                                                   
1797                      ENDIF  ! VOC
1798
1799                   ENDIF  ! PM
1800
1801                ENDIF  ! emis_distribution
1802
1803             ENDDO  ! m
1804
1805
1806
1807!
1808!-- USM surfaces
1809
1810             DO  m = 1, surf_usm_h%ns
1811
1812                i = surf_usm_h%i(m)
1813                j = surf_usm_h%j(m)
1814                k = surf_usm_h%k(m)
1815
1816                IF  ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
1817
1818!
1819!-- PMs
1820                   IF  ( TRIM(spc_names(match_spec_model(ispec))) == "PM1"     .OR.    &
1821                         TRIM(spc_names(match_spec_model(ispec))) == "PM25"    .OR.    &
1822                         TRIM(spc_names(match_spec_model(ispec))) == "PM10" )  THEN
1823                   
1824                      surf_usm_h%cssws(match_spec_model(ispec),m) =            &    ! kg/m2/s * kg/m3
1825                               emis_distribution(1,j,i,ispec)*                 &    ! kg/m2/s
1826                               rho_air(k)                                           ! kg/m^3
1827
1828                   ELSE
1829                     
1830!
1831!-- VOCs
1832                      IF  ( len_index_voc > 0                                         .AND.      &
1833                            emt_att%species_name(match_spec_input(ispec)) == "VOC" )  THEN
1834
1835                         surf_usm_h%cssws(match_spec_model(ispec),m) =         &   ! ppm/s * m * kg/m3
1836                               emis_distribution(1,j,i,ispec) *                &   ! m2/s
1837                               conv_to_ratio(k,j,i) *                          &   ! m^3/mole
1838                               ratio2ppm *                                     &   ! ppm
1839                               rho_air(k)                                          ! kg/m^3   
1840
1841!
1842!-- Other species
1843                      ELSE
1844
1845                         surf_usm_h%cssws(match_spec_model(ispec),m) =         &   ! ppm/s * m * kg/m3
1846                               emis_distribution(1,j,i,ispec) *                &   ! kg/m2/s
1847                               ( 1.0_wp / emt_att%xm(ispec) ) *                &   ! mole/kg
1848                               conv_to_ratio(k,j,i) *                          &   ! m^3/mole
1849                               ratio2ppm*                                      &   ! ppm
1850                               rho_air(k)                                          ! kg/m^3   
1851
1852
1853                      ENDIF  ! VOC
1854
1855                   ENDIF  ! PM
1856
1857                ENDIF  ! emis_distribution
1858
1859             ENDDO  ! m
1860
1861          ENDDO
1862
1863       ENDIF
1864
1865!
1866!-- Deallocate time_factor in case of DEFAULT mode)
1867
1868       IF  ( ALLOCATED (time_factor) )  DEALLOCATE (time_factor)
1869
1870   ENDIF
1871
1872 END SUBROUTINE chem_emissions_setup
1873
1874
1875!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1876!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1877!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1878!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1879!!
1880!! 20200203 (ECC) - ON DEMAND EMISSION UPDATE MODE
1881!!
1882!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1883!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1884!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1885!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1886
1887
1888!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1889!!
1890!! WRAPPER / INTERFACE FUNCTIONS
1891!!
1892!! NOTE - I find using an explicity wrapper provides much better flow control
1893!!          over an interface block
1894!!
1895!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1896
1897!
1898!-- 20200203 (ECC)
1899!
1900!------------------------------------------------------------------------------!
1901! Description:
1902! ------------
1903!>  interface for initiation of emission arrays based on emission LOD
1904!
1905!------------------------------------------------------------------------------!
1906
1907 SUBROUTINE chem_emissions_header_init
1908
1909    IMPLICIT NONE
1910
1911    SELECT CASE ( emiss_lod )
1912       CASE ( 0 )
1913! do nothing at the moment
1914       CASE ( 1 )
1915! do nothing at the moment
1916       CASE ( 2 )
1917          CALL chem_emissions_header_init_lod2
1918    END SELECT
1919
1920 END SUBROUTINE chem_emissions_header_init
1921
1922
1923!
1924!-- 20200203 (ECC)
1925!
1926!------------------------------------------------------------------------------!
1927! Description:
1928! ------------
1929!>  interface for initiation of emission arrays based on emission LOD
1930!
1931!------------------------------------------------------------------------------!
1932
1933 SUBROUTINE chem_emissions_update_on_demand
1934
1935    IMPLICIT NONE
1936
1937    SELECT CASE ( emiss_lod )
1938       CASE ( 0 )
1939! do nothing at the moment
1940       CASE ( 1 )
1941! do nothing at the moment
1942       CASE ( 2 )
1943          CALL chem_emissions_update_on_demand_lod2
1944    END SELECT
1945
1946 END SUBROUTINE        ! chem_emisisons_update_on_demand
1947
1948 
1949!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1950!!
1951!! SUBROUTINES SPECIFIC FOR LOD 2
1952!!
1953!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1954
1955!
1956!-- 20200203 (ECC)
1957!
1958!------------------------------------------------------------------------------!
1959! Description:
1960! ------------
1961!>  Initiates header for emissions data attributes for LOD 2
1962!------------------------------------------------------------------------------!
1963
1964 SUBROUTINE chem_emissions_header_init_lod2
1965
1966    USE control_parameters,                                                 &
1967        ONLY: coupling_char, message_string
1968
1969    USE netcdf_data_input_mod,                                              &
1970        ONLY: chem_emis_att,                                                &
1971              open_read_file, close_input_file,                             &
1972              get_dimension_length, get_variable,                           &
1973              get_attribute
1974
1975    IMPLICIT NONE
1976
1977    INTEGER(iwp) ::  ncid      !< chemistry file netCDF handle
1978    INTEGER(iwp) ::  att_lod   !<  lod attribute in chemistry file
1979
1980    IF  ( debug_output )  &
1981       CALL debug_message( 'chem_emissions_header_init_lod2', 'start' )
1982
1983!
1984!-- opens _chemistry input file and obtain header information
1985
1986    CALL open_read_file ( TRIM(input_file_chem) // TRIM(coupling_char), ncid )
1987!
1988!-- check if LOD in chemistry file matches LOD in namelist
1989
1990    CALL get_attribute ( ncid, 'lod', att_lod, .TRUE. )
1991   
1992    IF  ( att_lod /= emiss_lod )  THEN
1993       message_string = ''   ! to get around unused variable warning / error
1994       WRITE ( message_string, * )                                           &
1995             'LOD mismatch between namelist (emiss_lod) and',                &
1996             CHAR(10), '                    ',                               &
1997             'chemistry input file (global attributes>lod)'
1998       CALL message( 'chem_emissions_header_init_lod2', 'CM0468', 1, 2, 0, 6, 0 )
1999    ENDIF
2000!
2001!-- obtain unit conversion factor
2002
2003    CALL get_attribute ( ncid, 'units', chem_emis_att%units, .FALSE., "emission_values" )
2004    conversion_factor = chem_emissions_convert_base_units ( chem_emis_att%units )
2005!
2006!-- obtain header attributes
2007
2008    CALL chem_emissions_init_species    ( ncid )
2009    CALL chem_emissions_init_timestamps ( ncid )
2010!
2011!-- done reading file
2012
2013    CALL close_input_file (ncid)
2014
2015!
2016!-- set previous timestamp index to something different
2017!-- to trigger a read event later on
2018
2019    previous_timestamp_index = -1
2020 
2021    IF  ( debug_output )  &
2022       CALL debug_message( 'chem_emissions_header_init_lod2', 'end' )
2023
2024 END SUBROUTINE chem_emissions_header_init_lod2
2025
2026!
2027!-- 20200203 (ECC)
2028!
2029!------------------------------------------------------------------------------!
2030! Description:
2031! ------------
2032!> Reads emission data on demand for LOD2
2033!------------------------------------------------------------------------------!
2034
2035 SUBROUTINE chem_emissions_update_on_demand_lod2
2036
2037    USE control_parameters,                                                  &
2038        ONLY: coupling_char,                                                 &
2039              time_since_reference_point
2040
2041    USE netcdf_data_input_mod,                                               &
2042        ONLY: chem_emis_att,                                                 &
2043              open_read_file, close_input_file, get_variable
2044
2045    USE arrays_3d,                                                           &
2046        ONLY: pt, hyp
2047
2048    USE surface_mod,                                                         &
2049        ONLY: surf_def_h, surf_lsm_h, surf_usm_h
2050
2051
2052    IMPLICIT NONE
2053
2054    CHARACTER(LEN=80) ::  this_timestamp   !< writes out timestamp
2055
2056    INTEGER(iwp) ::  i,j,k,m                !< generic counters
2057    INTEGER(iwp) ::  kmatch                 !< index of matched species
2058    INTEGER(iwp) ::  ncid                   !< netCDF file handle (chemistry file)
2059    INTEGER(iwp) ::  time_index_location    !< location of most recent timestamp
2060
2061    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:) ::  emissions_raw  !< raw emissions data
2062    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)     ::  emis_distrib   !< surface emissions
2063    REAL(wp), ALLOCATABLE, DIMENSION(:,:)       ::  mass2mole      !< conversion factor mass 2 molar (ppm) flux
2064    REAL(wp), ALLOCATABLE, DIMENSION(:,:)       ::  cssws_def_h    !< dummy default surface array
2065    REAL(wp), ALLOCATABLE, DIMENSION(:,:)       ::  cssws_lsm_h    !< dummy LSM surface array
2066    REAL(wp), ALLOCATABLE, DIMENSION(:,:)       ::  cssws_usm_h    !< dummy USM surface array
2067
2068    IF  ( debug_output )  &
2069       CALL debug_message ( 'chem_emissions_update_on_demand_lod2', 'start' )
2070!
2071!-- obtain current timestamp and locate index
2072!-- for most recent timestamp element
2073!-- end subroutine (RETURN) if it is still the same
2074!-- index as the existing time index
2075
2076    this_timestamp = ''   ! string must be initiated before using
2077    CALL get_date_time( time_since_reference_point, date_time_str=this_timestamp )
2078
2079    time_index_location = chem_emissions_locate_timestep                     & 
2080                             ( this_timestamp, timestamps,                   &
2081                               1, chem_emis_att%dt_emission )
2082
2083    IF  ( time_index_location == previous_timestamp_index )  RETURN
2084
2085!
2086!-- begin extract emission data for matched species from netCDF file
2087
2088    previous_timestamp_index = time_index_location
2089
2090    ALLOCATE ( emis_distrib(n_matched_vars,nys:nyn,nxl:nxr) )
2091    emis_distrib = 0.0_wp
2092
2093!
2094!-- open netCDF file and allocate temp memory
2095
2096    CALL open_read_file( TRIM(input_file_chem) // TRIM(coupling_char), ncid )
2097    ALLOCATE( emissions_raw(1,1,nys:nyn,nxl:nxr,1) )
2098
2099    DO k = 1, n_matched_vars
2100!
2101!-- get index for matching species
2102
2103       kmatch = chem_emissions_locate_species (                               &
2104                               spc_names(match_spec_model(k)),                &
2105                               chem_emis_att%species_name )
2106!
2107!-- extract variable as-is
2108!-- (note C index notations for nx and ny due to MPI and
2109!-- reversed index dimension order for netCDF Fortran API)
2110
2111       emissions_raw = 0.0_wp
2112
2113       CALL get_variable ( ncid, 'emission_values', emissions_raw,            &
2114                           kmatch, nxl+1, nys+1, 1, time_index_location,      &
2115                           1, nxr-nxl+1, nyn-nys+1, 1, 1, .FALSE. )
2116!
2117!-- transfer emission data
2118
2119       DO j = nys,nyn
2120         DO i = nxl,nxr
2121            emis_distrib(k,j,i) = emissions_raw(1,1,j,i,1) * conversion_factor
2122         ENDDO
2123       ENDDO
2124
2125    ENDDO  ! k = n_matched_vars
2126!
2127!-- netCDF handle and temp memory no longer needed
2128
2129    DEALLOCATE( emissions_raw )
2130    CALL close_input_file( ncid )
2131!
2132!-- Set emis_dt since chemistry ODEs can be stiff, the option
2133!-- to solve them at every RK substep is present to help improve
2134!-- stability should the need arises
2135 
2136    dt_emis = dt_3d
2137
2138    IF  ( call_chem_at_all_substeps )                                          &
2139       dt_emis = dt_emis * weight_pres(intermediate_timestep_count)
2140!
2141!-- calculate conversion factor from mass flux to molar flux (mixing ratio)
2142
2143    ALLOCATE ( mass2mole(nys:nyn,nxl:nxr) )
2144    mass2mole = 0.0_wp
2145
2146    DO i = nxl, nxr
2147       DO j = nys, nyn
2148          mass2mole(j,i) = mass_2_molar_flux ( hyp(nzb), pt(nzb,j,i) )
2149       ENDDO
2150    ENDDO
2151
2152!
2153!-- calculate surface fluxes
2154!-- NOTE - For some reason I cannot pass surf_xxx%cssws as output argument
2155!--        into subroutine assign_surface_flux ( ).  The contents got mixed up
2156!--        once the subroutine is finished.  I don't know why and I don't have
2157!--        time to investigate.  As workaround I declared dummy variables
2158!--        and reassign them one by one (i.e., in a loop)
2159!-- ECC 20200206
2160
2161!
2162!-- allocate and initialize dummy surface fluxes
2163
2164    ALLOCATE( cssws_def_h(n_matched_vars,surf_def_h(0)%ns) )
2165    cssws_def_h = 0.0_wp
2166
2167    ALLOCATE( cssws_lsm_h(n_matched_vars,surf_lsm_h%ns) )
2168    cssws_lsm_h = 0.0_wp
2169
2170    ALLOCATE( cssws_usm_h(n_matched_vars,surf_usm_h%ns) )
2171    cssws_usm_h = 0.0_wp
2172
2173!
2174!-- assign and transfer emissions as surface fluxes
2175
2176    CALL assign_surface_flux ( cssws_def_h, surf_def_h(0)%ns,                  &
2177                               surf_def_h(0)%j, surf_def_h(0)%i,               &
2178                               emis_distrib, mass2mole )
2179
2180
2181    CALL assign_surface_flux ( cssws_lsm_h, surf_lsm_h%ns,                     &
2182                               surf_lsm_h%j, surf_lsm_h%i,                     &
2183                               emis_distrib, mass2mole )
2184
2185
2186    CALL assign_surface_flux ( cssws_usm_h, surf_usm_h%ns,                     &
2187                               surf_usm_h%j, surf_usm_h%i,                     &
2188                               emis_distrib, mass2mole )
2189
2190    DO k = 1, n_matched_vars
2191
2192       DO m = 1, surf_def_h(0)%ns
2193          surf_def_h(0)%cssws(k,m) = cssws_def_h(k,m)
2194       ENDDO
2195
2196       DO m = 1, surf_lsm_h%ns
2197          surf_lsm_h%cssws(k,m)    = cssws_lsm_h(k,m)
2198       ENDDO
2199
2200       DO m = 1, surf_usm_h%ns
2201          surf_usm_h%cssws(k,m)    = cssws_usm_h(k,m)
2202       ENDDO
2203
2204    ENDDO
2205
2206!
2207!-- cleaning up
2208
2209    DEALLOCATE( cssws_def_h )
2210    DEALLOCATE( cssws_lsm_h )
2211    DEALLOCATE( cssws_usm_h )
2212
2213    DEALLOCATE ( emis_distrib )
2214    DEALLOCATE ( mass2mole )
2215
2216    IF  ( debug_output )  &
2217       CALL debug_message ( 'chem_emissions_update_on_demand_lod2', 'end' )
2218
2219 END SUBROUTINE        ! chem_emissions_update_on_demand_lod2
2220
2221 
2222!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2223!!
2224!! AUXILIARY SUBROUTINES
2225!!
2226!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2227
2228!
2229!-- 20200203 (ECC)
2230!
2231!------------------------------------------------------------------------------!
2232! Description:
2233! ------------
2234!> look for matched species between emissions attributes and selected
2235!> chemical mechanisms and determine corresponding molecular weights
2236!------------------------------------------------------------------------------!
2237
2238 SUBROUTINE chem_emissions_init_species ( ncid )
2239
2240    USE netcdf_data_input_mod,                                              &
2241        ONLY: chem_emis_att,                                                &
2242              open_read_file, close_input_file,                             &
2243              get_dimension_length, get_variable
2244
2245    IMPLICIT NONE
2246
2247    INTEGER(iwp) :: ispec  !< generic counter 4 species
2248
2249    INTEGER(iwp), INTENT(IN) :: ncid   !< netcdf file ID
2250
2251    IF  ( debug_output )  &
2252       CALL debug_message( 'chem_emissions_init_species', 'start' )
2253!
2254!- assign species
2255
2256    CALL get_dimension_length ( ncid, chem_emis_att%n_emiss_species, 'nspecies' )
2257    CALL get_variable ( ncid, 'emission_name', chem_emis_att%species_name,   &
2258                              chem_emis_att%n_emiss_species )
2259!
2260!- backward compatibility for salsa_mod (ECC)
2261    chem_emis_att%nspec = chem_emis_att%n_emiss_species
2262!
2263!-- get a list of matched species between emission_attributes and
2264!-- selected chemical mechanism
2265
2266    emission_output_required = .FALSE.
2267    CALL  chem_emissions_match( chem_emis_att, n_matched_vars )
2268
2269!
2270!-- if matched species found (at least 1)
2271!-- allocate memory for emission attributes
2272!-- assign molecular masses [kg/mol]
2273!-- see chemistry_model_mod.f90 for reference
2274
2275    IF  ( n_matched_vars > 0 )  THEN
2276
2277       emission_output_required = .TRUE.
2278
2279       ALLOCATE( chem_emis_att%xm(n_matched_vars) )
2280
2281       DO  ispec = 1, n_matched_vars
2282          chem_emis_att%xm(ispec) = 1.0_wp
2283          SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
2284             CASE ( 'SO2'  );  chem_emis_att%xm(ispec) = xm_S + xm_O * 2
2285             CASE ( 'SO4'  );  chem_emis_att%xm(ispec) = xm_S + xm_O * 4
2286             CASE ( 'NO'   );  chem_emis_att%xm(ispec) = xm_N + xm_O
2287             CASE ( 'NO2'  );  chem_emis_att%xm(ispec) = xm_N + xm_O * 2
2288             CASE ( 'NH3'  );  chem_emis_att%xm(ispec) = xm_N + xm_H * 3
2289             CASE ( 'CO'   );  chem_emis_att%xm(ispec) = xm_C + xm_O
2290             CASE ( 'CO2'  );  chem_emis_att%xm(ispec) = xm_C + xm_O * 2
2291             CASE ( 'CH4'  );  chem_emis_att%xm(ispec) = xm_C + xm_H * 4
2292             CASE ( 'HNO3' );  chem_emis_att%xm(ispec) = xm_H + xm_N + xm_O*3
2293          END SELECT
2294       ENDDO
2295   
2296    ENDIF  ! IF ( n_matched_vars > 0 )
2297
2298    IF  ( debug_output )  &
2299       CALL debug_message( 'chem_emissions_init_species', 'end' )
2300
2301 END SUBROUTINE chem_emissions_init_species
2302
2303 
2304!
2305!-- 20200203 (ECC)
2306!
2307!------------------------------------------------------------------------------!
2308! Description:
2309! ------------
2310!> extract timestamps from netCDF input
2311!------------------------------------------------------------------------------!
2312
2313 SUBROUTINE chem_emissions_init_timestamps ( ncid )
2314
2315    USE control_parameters,                                                 &
2316        ONLY: message_string
2317
2318    USE netcdf_data_input_mod,                                              &
2319        ONLY: chem_emis_att,                                                &
2320              open_read_file, close_input_file,                             &
2321              get_dimension_length, get_variable
2322
2323    IMPLICIT NONE
2324
2325    INTEGER(iwp) ::  fld_len   !< string field length
2326    INTEGER(iwp) ::  itime     !< generic counter (4 species)
2327
2328    INTEGER(iwp), INTENT(IN) :: ncid   !< netcdf file handle
2329
2330    IF  ( debug_output )  &
2331       CALL debug_message( 'chem_emissions_read_timestamps', 'start' )
2332!
2333!-- import timestamps from netCDF input
2334
2335    CALL get_dimension_length ( ncid, chem_emis_att%dt_emission, 'time' )
2336    CALL get_dimension_length ( ncid, fld_len, 'field_length' )
2337    CALL get_variable ( ncid, 'timestamp', timestamps, chem_emis_att%dt_emission, fld_len )
2338!
2339!-- throw error at first instance of timestamps
2340!-- not in listed in chronological order
2341
2342    DO itime = 2,chem_emis_att%dt_emission
2343
2344       IF ( timestamps(itime-1) > timestamps(itime) )  THEN
2345
2346           WRITE( message_string, * ) &
2347                'input timestamps not in chronological order for',             &
2348                CHAR(10), '                    ',                              &
2349                'index ', (itime-1), ' : ', TRIM(timestamps(itime-1)), ' and', &
2350                CHAR(10), '                    ',                              &
2351                'index ', (itime),   ' : ', TRIM(timestamps(itime)) 
2352
2353          CALL message( 'chem_emissions_read_timestamps', 'CM0469', 1, 2, 0, 6, 0 )
2354
2355       ENDIF
2356
2357    ENDDO
2358
2359    IF  ( debug_output )  &
2360       CALL debug_message( 'chem_emissions_read_timestamps', 'end' )
2361
2362 END SUBROUTINE chem_emissions_init_timestamps
2363
2364
2365!
2366!-- 20200203 (ECC)
2367!
2368!------------------------------------------------------------------------------!
2369! Description:
2370! ------------
2371!> assign emissions as surface fluxes
2372!
2373!> NOTE:  For arguments, I originally wanted to use unspecified dimensions,
2374!>        but I could not get this to work properly, hence the dimensioned
2375!>        array arguments
2376!------------------------------------------------------------------------------!
2377
2378SUBROUTINE assign_surface_flux ( surf_array, nsurfs, surf_j, surf_i,         &
2379                                 emis_dist, conv_mole )
2380
2381    USE arrays_3d,                                                           &
2382        ONLY: rho_air
2383
2384    USE netcdf_data_input_mod,                                               &
2385        ONLY: chem_emis_att
2386
2387    USE surface_mod   !< for surf_type
2388
2389    IMPLICIT NONE
2390!
2391!-- input arguments
2392
2393    INTEGER(iwp),                    INTENT(IN) ::  nsurfs   !< # surfaces in surf_array
2394    INTEGER(iwp), DIMENSION(nsurfs), INTENT(IN) ::  surf_i   !< i indices 4 surf. elements
2395    INTEGER(iwp), DIMENSION(nsurfs), INTENT(IN) ::  surf_j   !< j indices 4 surf. elements
2396
2397    REAL(wp), DIMENSION(nys:nyn,nxl:nxr),                INTENT(IN)    ::  conv_mole   !< conv. 2 molar flux
2398    REAL(wp), DIMENSION(n_matched_vars,nys:nyn,nxl:nxr), INTENT(IN)    ::  emis_dist   !< surf. emissions
2399
2400    REAL(wp), DIMENSION(n_matched_vars,nsurfs),          INTENT(INOUT) ::  surf_array  !< surface listing
2401
2402!
2403!-- parameters (magic numbers)
2404
2405    CHARACTER(LEN=2),  PARAMETER ::  sp_PM  = 'PM'   !< id string for all PMs
2406    CHARACTER(LEN=3),  PARAMETER ::  sp_VOC = 'VOC'  !< id string for VOC
2407
2408    REAL(wp), PARAMETER ::  mol2ppm = 1.0E+06_wp     !< conversion from mole 2 ppm
2409!
2410!-- local variables
2411
2412    CHARACTER(LEN=80) ::  this_species_name  !< matched species name
2413
2414    INTEGER(iwp) ::  i,j,k,m        !< generic counters
2415
2416    REAL(wp) ::  flux_conv_factor   !< conversion factor
2417
2418    IF  ( debug_output )  &
2419       CALL debug_message( 'chem_emissions_header_init_lod2', 'start' )
2420
2421    DO k = 1, n_matched_vars
2422
2423       this_species_name = spc_names(k)  !< species already matched
2424
2425       DO m = 1, nsurfs
2426
2427          j = surf_j(m)   ! get surface coordinates
2428          i = surf_i(m)
2429
2430!
2431!-- calculate conversion factor depending on emission species type
2432
2433          flux_conv_factor = rho_air(nzb)
2434!
2435!-- account for conversion to different types of emisison species
2436
2437          IF       ( TRIM(this_species_name(1:LEN(sp_PM)))  == sp_PM  )  THEN
2438
2439                ! do nothing (use mass flux directly)
2440
2441          ELSE IF  ( TRIM(this_species_name(1:LEN(sp_VOC))) == sp_VOC )  THEN
2442
2443             flux_conv_factor = flux_conv_factor *                           &
2444                                conv_mole(j,i) * mol2ppm
2445
2446          ELSE
2447
2448             flux_conv_factor = flux_conv_factor *                           &
2449                                conv_mole(j,i) * mol2ppm /                   &
2450                                chem_emis_att%xm(k)
2451
2452          ENDIF
2453!
2454!-- finally assign surface flux
2455
2456          surf_array(k,m) = emis_dist(k,j,i) * flux_conv_factor
2457
2458       ENDDO   ! m = 1, nsurfs
2459
2460    ENDDO   ! k = 1, n_matched_vars
2461
2462
2463    IF  ( debug_output )  &
2464       CALL debug_message( 'chem_emissions_header_init_lod2', 'end' )
2465   
2466END SUBROUTINE assign_surface_flux
2467
2468 
2469!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2470!!
2471!! AUXILIARY FUNCTIONS
2472!!
2473!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2474
2475!
2476!-- 20200203 (ECC)
2477!
2478!------------------------------------------------------------------------------!
2479! Description:
2480! ------------
2481!> given incoming flux units ( mass / area / time ) provide single-valued
2482!> conversion factor to ( kg / m2 / s )
2483!------------------------------------------------------------------------------!
2484
2485 FUNCTION chem_emissions_convert_base_units ( units_in ) RESULT ( conv_factor )
2486
2487    IMPLICIT NONE
2488!
2489!-- function arguments
2490
2491    REAL(wp) ::  conv_factor                       !< convertion factor
2492
2493    CHARACTER(LEN=*), INTENT(IN) ::  units_in      !< incoming units (ie emt_att%units)
2494!
2495!-- parameters (magic numbers)
2496
2497    INTEGER(iwp), PARAMETER ::  up2lo = 32         !< convert letter to lower case
2498!
2499!-- base unit conversion factors (should be self-explanatory)
2500
2501    REAL(wp), PARAMETER ::  g_to_kg       =     1.0E-03_wp
2502    REAL(wp), PARAMETER ::  miug_to_kg    =     1.0E-09_wp
2503    REAL(wp), PARAMETER ::  tons_to_kg    =   100.0_wp
2504
2505    REAL(wp), PARAMETER ::  hour_per_year =  8760.0_wp
2506    REAL(wp), PARAMETER ::  s_per_hour    =  3600.0_wp
2507    REAL(wp), PARAMETER ::  s_per_day     = 86400.0_wp
2508
2509    REAL(wp), PARAMETER ::  day_to_s      =     1.0_wp / s_per_day
2510    REAL(wp), PARAMETER ::  hour_to_s     =     1.0_wp / s_per_hour
2511    REAL(wp), PARAMETER ::  year_to_s     =     1.0_wp / s_per_hour / hour_per_year 
2512!
2513!-- local variables
2514
2515    CHARACTER(LEN=LEN(units_in)) ::  units_in_lo   !< units in lower case
2516
2517    INTEGER(iwp) ::  j,k                           !< generic counters
2518    INTEGER(iwp) ::  str_len                       !< length of unit string
2519!
2520!-- turn units string to lower case
2521 
2522    units_in_lo = ''
2523    str_len = LEN(TRIM(units_in))
2524
2525    DO k = 1,str_len
2526       j = IACHAR( units_in(k:k) )
2527       units_in_lo(k:k) = ACHAR(j)
2528       IF  ( (j>=IACHAR("A")) .AND. (j<=IACHAR("Z")) )   &
2529          units_in_lo(k:k) = ACHAR ( j + up2lo )
2530    ENDDO
2531
2532    conv_factor = 1.0_wp     !< default value (kg/m2/s)
2533
2534    SELECT CASE ( TRIM( units_in_lo ) )
2535       CASE ( 'kg/m2/s'            );   conv_factor = 1.0_wp
2536       CASE ( 'kg/m2/hour'         );   conv_factor = hour_to_s
2537       CASE ( 'kg/m2/day'          );   conv_factor = day_to_s
2538       CASE ( 'kg/m2/year'         );   conv_factor = year_to_s
2539       CASE ( 'ton/m2/s'           );   conv_factor = tons_to_kg
2540       CASE ( 'ton/m2/hour'        );   conv_factor = tons_to_kg * hour_to_s
2541       CASE ( 'ton/m2/day'         );   conv_factor = tons_to_kg * day_to_s
2542       CASE ( 'ton/m2/year'        );   conv_factor = tons_to_kg * year_to_s
2543       CASE ( 'g/m2/s'             );   conv_factor = g_to_kg
2544       CASE ( 'g/m2/hour'          );   conv_factor = g_to_kg * hour_to_s
2545       CASE ( 'g/m2/day'           );   conv_factor = g_to_kg * day_to_s
2546       CASE ( 'g/m2/year'          );   conv_factor = g_to_kg * year_to_s
2547       CASE ( 'micrograms/m2/s'    );   conv_factor = miug_to_kg
2548       CASE ( 'micrograms/m2/hour' );   conv_factor = miug_to_kg * hour_to_s
2549       CASE ( 'micrograms/m2/day'  );   conv_factor = miug_to_kg * day_to_s
2550       CASE ( 'micrograms/m2/year' );   conv_factor = miug_to_kg * year_to_s
2551       CASE DEFAULT
2552          message_string = ''   ! to get around unused variable warning / error
2553          WRITE ( message_string, * ) 'Specified emission units (',          &
2554                                      TRIM(units_in),                        &
2555                                      ') not recognized in PALM-4U'
2556          CALL message ( 'chem_emission_convert_units', 'CM0446', 2, 2, 0, 6, 0 )
2557    END SELECT
2558
2559 END FUNCTION chem_emissions_convert_base_units
2560
2561
2562!
2563!-- 20200203 (ECC)
2564!
2565!------------------------------------------------------------------------------!
2566! Description:
2567! ------------
2568!> calculates conversion factor from mass flux to ppm (molar flux)
2569!------------------------------------------------------------------------------!
2570
2571 FUNCTION mass_2_molar_flux ( rhogh, theta ) RESULT ( conv_factor )
2572
2573    USE basic_constants_and_equations_mod,                                   &
2574        ONLY: p_0, rgas_univ, rd_d_cp
2575
2576    IMPLICIT NONE
2577!
2578!-- function arguments
2579
2580    REAL(wp)             ::  conv_factor   !< conversion factor
2581    REAL(wp), INTENT(IN) ::  rhogh         !< hydrostatic pressure
2582    REAL(wp), INTENT(IN) ::  theta         !< potential temperature
2583
2584    conv_factor = ( rgas_univ / rhogh ) * theta * ( ( rhogh / p_0 ) ** rd_d_cp )
2585
2586 END FUNCTION mass_2_molar_flux
2587
2588
2589!
2590!-- 20200203 (ECC)
2591!
2592!------------------------------------------------------------------------------!
2593! Description:
2594! ------------
2595!> given target sepecies locate index in species array
2596!> returns 0 if none is found
2597!------------------------------------------------------------------------------!
2598
2599 FUNCTION chem_emissions_locate_species ( this_species, species_array )      &
2600          RESULT ( species_index )
2601
2602    IMPLICIT NONE
2603!
2604!-- function arguments
2605
2606    INTEGER(iwp) ::  species_index   !> index matching species
2607
2608    CHARACTER(LEN=*),  INTENT(IN) ::  this_species       !> target species
2609    CHARACTER(LEN=25), INTENT(IN) ::  species_array(:)   !> array of species
2610!
2611!-- local variables
2612
2613    INTEGER(iwp) :: k           !> generic counter
2614    INTEGER(iwp) :: n_species   !> number of species in species_array
2615
2616    n_species = SIZE( species_array, 1 )
2617
2618    DO k = 1, n_species
2619       IF ( TRIM(species_array(k)) == TRIM(this_species) )  EXIT
2620    ENDDO
2621
2622    species_index = 0   !> assume no matching index is found
2623
2624    IF ( TRIM(species_array(k)) == TRIM(this_species) )  specieS_index = k
2625
2626 END FUNCTION chem_emissions_locate_species
2627
2628
2629!
2630!-- 20200203 (ECC)
2631!
2632!------------------------------------------------------------------------------!
2633! Description:
2634! ------------
2635!> given target timestamp locate most recent timestep in timestamp array
2636!> using bisection search (since array is sorted)
2637!------------------------------------------------------------------------------!
2638
2639 RECURSIVE FUNCTION chem_emissions_locate_timestep                           &
2640                       ( this_timestamp, timestamp_array,                    &
2641                         lower_bound, upper_bound )                          &
2642                       RESULT ( timestamp_index )
2643
2644!
2645!-- function arguments
2646
2647    INTEGER(iwp) ::  timestamp_index   !> index for most recent timestamp in timestamp_array
2648
2649    CHARACTER(LEN=*),   INTENT(IN) ::  this_timestamp       !> target timestamp
2650    CHARACTER(LEN=512), INTENT(IN) ::  timestamp_array(:)   !> array of timestamps
2651
2652    INTEGER(iwp), INTENT(IN) ::  lower_bound, upper_bound   !> timestamp_array index bounds
2653
2654!
2655!-- local variables
2656
2657    INTEGER(iwp) :: k0,km,k1          !> lower, central, and upper index bounds
2658!
2659!-- assign bounds
2660
2661    k0 = lower_bound
2662    k1 = upper_bound
2663!
2664!-- make sure k1 is always not smaller than k0
2665
2666    IF  ( k0 > k1 )  THEN
2667       k0 = upper_bound
2668       k1 = lower_bound
2669    ENDIF
2670!
2671!-- make sure k0 and k1 stay within array bounds by timestamp_array
2672
2673    IF  ( k0 < 1                       )  k0 = 1
2674    IF  ( k1 > SIZE(timestamp_array,1) )  k1 = SIZE(timestamp_array,1)
2675!
2676!-- terminate if target is contained within 2 consecutive indices
2677!-- otherwise calculate central bound (km) and determine new
2678!-- index bounds for the next iteration
2679
2680    IF ( ( k1 - k0 ) > 1 )  THEN
2681       km = ( k0 + k1 ) / 2
2682       IF  ( TRIM(this_timestamp) > TRIM(timestamp_array(km)) )  THEN
2683          k0 = km
2684       ELSE
2685          k1 = km
2686       ENDIF
2687       timestamp_index = chem_emissions_locate_timestep &
2688                         ( this_timestamp, timestamp_array, k0, k1 )
2689    ELSE
2690       timestamp_index = k0
2691    ENDIF
2692
2693 END FUNCTION chem_emissions_locate_timestep
2694
2695
2696!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2697!!
2698!! END OF MODULE
2699!!
2700!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2701
2702 END MODULE chem_emissions_mod
Note: See TracBrowser for help on using the repository browser.