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

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

Some interface calls moved to module_interface + cleanup

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