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

Last change on this file since 4180 was 4180, checked in by scharf, 5 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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