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

Last change on this file since 4598 was 4559, checked in by raasch, 11 months ago

files re-formatted to follow the PALM coding standard

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