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

Last change on this file since 4154 was 4154, checked in by suehring, 22 months ago

Bugfix in chemistry decycling; Replace global arrays to setup chemical emissions by local ones

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