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

Last change on this file since 4144 was 4144, checked in by raasch, 23 months ago

relational operators .EQ., .NE., etc. replaced by ==, /=, etc.

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