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

Last change on this file since 4218 was 4218, checked in by knoop, 5 years ago

Unused routine chem_emissions_check_parameters in chem_emissions_mod.f90 commented out due to uninitialized content

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