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

Last change on this file since 4867 was 4828, checked in by Giersch, 4 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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