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

Last change on this file since 3831 was 3831, checked in by forkel, 2 years ago

added nvar to USE chem_gasphase_mod - nvar will not be included in chem_modules anymore soon

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