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

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

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

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