source: palm/trunk/SOURCE/chem_emissions_mod.f90

Last change on this file was 4887, checked in by banzhafs, 4 months ago

Unnecessary comments removed from chemistry_model_mod and chem_emissions_mod

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