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

Last change on this file since 3780 was 3772, checked in by suehring, 6 years ago

Assure that emissions are only at natural (pavement) surfaces in case of parametrized emission mode, not on urban surfaces

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