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

Last change on this file since 3984 was 3968, checked in by suehring, 5 years ago

Updates from chemistriy branched merged into trunk: code cleaning and formatting, code structure optimizations

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