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

Last change on this file since 3819 was 3788, checked in by banzhafs, 6 years ago

Removed unused variables from chem_emissions_mod

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