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

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

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

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