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

Last change on this file since 3664 was 3611, checked in by banzhafs, 5 years ago

chem_emissions_mod and chem_modules update to comply PALM coding rules

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