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

Last change on this file since 4173 was 4157, checked in by suehring, 5 years ago

chem_emissions: Replace global arrays also in mode_emis branch; diagnostic output: restructure initialization in order to work also when data output during spin-up is enabled; radiation: give informative message on raytracing distance only by core zero not by all cores

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