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

Last change on this file since 3685 was 3685, checked in by knoop, 3 years ago

Some interface calls moved to module_interface + cleanup

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