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

Last change on this file since 4118 was 4055, checked in by suehring, 5 years ago

chem_emissions_mod: Formatting adjustments; initialization of arrays fixed; univsersal gas constant moved to basic_constants_and_equations_mod.f90

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