Ignore:
Timestamp:
Jun 11, 2020 8:51:48 AM (12 months ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/chemistry_model_mod.f90

    r4550 r4559  
    11!> @file chemistry_model_mod.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    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/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 2017-2020 Leibniz Universitaet Hannover
    1817! Copyright 2017-2020 Karlsruhe Institute of Technology
    1918! Copyright 2017-2020 Freie Universitaet Berlin
    20 !------------------------------------------------------------------------------!
     19!--------------------------------------------------------------------------------------------------!
    2120!
    2221! Current revisions:
     
    2726! -----------------
    2827! $Id$
     28! file re-formatted to follow the PALM coding standard
     29!
     30! 4550 2020-05-29 15:22:13Z raasch
    2931! bugfix for reading local restart data
    30 ! 
     32!
    3133! 4544 2020-05-21 14:43:05Z raasch
    3234! conc_av changed from pointer to allocatable array, array spec_conc_av removed
    33 ! 
     35!
    3436! 4542 2020-05-19 15:45:12Z raasch
    3537! redundant if statement removed
    36 ! 
     38!
    3739! 4535 2020-05-15 12:07:23Z raasch
    3840! bugfix for restart data format query
    39 ! 
     41!
    4042! 4517 2020-05-03 14:29:30Z raasch
    4143! added restart with MPI-IO
    42 ! 
     44!
    4345! 4511 2020-04-30 12:20:40Z raasch
    4446! decycling replaced by explicit setting of lateral boundary conditions
     
    5052! use statement for exchange horiz added,
    5153! bugfix for call of exchange horiz 2d
    52 ! 
     54!
    5355! 4442 2020-03-04 19:21:13Z suehring
    54 ! Change order of dimension in surface array %frac to allow for better 
     56! Change order of dimension in surface array %frac to allow for better
    5557! vectorization.
    5658!
    5759! 4441 2020-03-04 19:20:35Z suehring
    5860! in subroutine chem_init (ECC)
    59 ! - allows different init paths emission data for legacy
    60 !   mode emission and on-demand mode
    61 ! in subroutine chem_init_internal (ECC)
    62 ! - reads netcdf file only when legacy mode is activated
    63 !   (i.e., emiss_read_legacy_mode = .TRUE.)
    64 !   otherwise file is read once at the beginning to obtain
    65 !   header information, and emission data are extracted on
    66 !   an on-demand basis
     61! - allows different init paths emission data for legacy mode emission and on-demand mode in
     62!   subroutine chem_init_internal (ECC)
     63! - reads netcdf file only when legacy mode is activated (i.e., emiss_read_legacy_mode = .TRUE.)
     64!   otherwise file is read once at the beginning to obtain header information, and emission data
     65!   are extracted on an on-demand basis
    6766!
    6867! 4372 2020-01-14 10:20:35Z banzhafs
    69 ! chem_parin : added handler for new namelist item emiss_legacy_read_mode (ECC)
    70 ! added messages
     68! chem_parin : added handler for new namelist item emiss_legacy_read_mode (ECC) added messages
    7169! CM0465 - legacy read mode selection
    7270! CM0466 - legacy read mode force selection
    7371! CM0467 - new read mode selection
    7472!
    75 ! 4370 2020-01-10 14:00:44Z raasch 
     73! 4370 2020-01-10 14:00:44Z raasch
    7674! vector directives added to force vectorization on Intel19 compiler
    7775!
    7876! 4346 2019-12-18 11:55:56Z motisi
    79 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    80 ! topography information used in wall_flags_static_0
    81 ! 
     77! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     78! information used in wall_flags_static_0
     79!
    8280! 4329 2019-12-10 15:46:36Z motisi
    8381! Renamed wall_flags_0 to wall_flags_static_0
    84 ! 
     82!
    8583! 4306 2019-11-25 12:04:48Z banzhafs
    8684! Corretion for r4304 commit
     
    9694!
    9795! 4273 2019-10-24 13:40:54Z monakurppa
    98 ! Add logical switches nesting_chem and nesting_offline_chem (both .TRUE.
    99 ! by default)
    100 !
     96! Add logical switches nesting_chem and nesting_offline_chem (both .TRUE. by default)
     97!
    10198! 4272 2019-10-23 15:18:57Z schwenkel
    102 ! Further modularization of boundary conditions: moved boundary conditions to
    103 ! respective modules
     99! Further modularization of boundary conditions: moved boundary conditions to respective modules
    104100!
    105101! 4268 2019-10-17 11:29:38Z schwenkel
    106102! Moving module specific boundary conditions from time_integration to module
    107 ! 
     103!
    108104! 4230 2019-09-11 13:58:14Z suehring
    109 ! Bugfix, initialize mean profiles also in restart runs. Also initialize
    110 ! array used for Runge-Kutta tendecies in restart runs. 
    111 ! 
     105! Bugfix, initialize mean profiles also in restart runs. Also initialize array used for Runge-Kutta
     106! tendecies in restart runs.
     107!
    112108! 4227 2019-09-10 18:04:34Z gronemeier
    113109! implement new palm_date_time_mod
    114 ! 
     110!
    115111! 4182 2019-08-22 15:20:23Z scharf
    116112! Corrected "Former revisions" section
    117 ! 
     113!
    118114! 4167 2019-08-16 11:01:48Z suehring
    119 ! Changed behaviour of masked output over surface to follow terrain and ignore
    120 ! buildings (J.Resler, T.Gronemeier)
    121 !
    122 !
     115! Changed behaviour of masked output over surface to follow terrain and ignore buildings
     116! (J.Resler, T.Gronemeier)
     117!
    123118! 4166 2019-08-16 07:54:21Z resler
    124119! Bugfix in decycling
    125 ! 
     120!
    126121! 4115 2019-07-24 12:50:49Z suehring
    127122! Fix faulty IF statement in decycling initialization
    128 ! 
     123!
    129124! 4110 2019-07-22 17:05:21Z suehring
    130 ! - Decycling boundary conditions are only set at the ghost points not on the
    131 !   prognostic grid points
    132 ! - Allocation and initialization of special advection flags cs_advc_flags_s
    133 !   used for chemistry. These are exclusively used for chemical species in
    134 !   order to distinguish from the usually-used flags which might be different
    135 !   when decycling is applied in combination with cyclic boundary conditions.
    136 !   Moreover, cs_advc_flags_s considers extended zones around buildings where
    137 !   first-order upwind scheme is applied for the horizontal advection terms,
    138 !   in order to overcome high concentration peaks due to stationary numerical
    139 !   oscillations caused by horizontal advection discretization.
    140 !
     125! - Decycling boundary conditions are only set at the ghost points not on the prognostic grid points
     126! - Allocation and initialization of special advection flags cs_advc_flags_s used for chemistry.
     127!   These are exclusively used for chemical species in order to distinguish from the usually-used
     128!   flags which might be different when decycling is applied in combination with cyclic boundary
     129!   conditions.
     130!   Moreover, cs_advc_flags_s considers extended zones around buildings where first-order upwind
     131!   scheme is applied for the horizontal advection terms, in order to overcome high concentration
     132!   peaks due to stationary numerical oscillations caused by horizontal advection discretization.
     133!
    141134! 4109 2019-07-22 17:00:34Z suehring
    142 ! Slightly revise setting of boundary conditions at horizontal walls, use
    143 ! data-structure offset index instead of pre-calculate it for each facing
    144 ! 
     135! Slightly revise setting of boundary conditions at horizontal walls, use data-structure offset
     136! index instead of pre-calculate it for each facing
     137!
    145138! 4080 2019-07-09 18:17:37Z suehring
    146139! Restore accidantly removed limitation to positive values
    147 ! 
     140!
    148141! 4079 2019-07-09 18:04:41Z suehring
    149 ! Application of monotonic flux limiter for the vertical scalar advection
    150 ! up to the topography top (only for the cache-optimized version at the
    151 ! moment).
    152 !
     142! Application of monotonic flux limiter for the vertical scalar advection up to the topography top
     143! (only for the cache-optimized version at the moment).
     144!
    153145! 4069 2019-07-01 14:05:51Z Giersch
    154 ! Masked output running index mid has been introduced as a local variable to
    155 ! avoid runtime error (Loop variable has been modified) in time_integration
    156 ! 
     146! Masked output running index mid has been introduced as a local variable to avoid runtime error
     147! (Loop variable has been modified) in time_integration
     148!
    157149! 4029 2019-06-14 14:04:35Z raasch
    158150! nest_chemistry option removed
    159 ! 
     151!
    160152! 4004 2019-05-24 11:32:38Z suehring
    161 ! in subroutine chem_parin check emiss_lod / mod_emis only
    162 ! when emissions_anthropogenic is activated in namelist (E.C. Chan)
    163 ! 
     153! in subroutine chem_parin check emiss_lod / mod_emis only when emissions_anthropogenic is activated
     154! in namelist (E.C. Chan)
     155!
    164156! 3968 2019-05-13 11:04:01Z suehring
    165 ! - added "emiss_lod" which serves the same function as "mode_emis"
    166 !   both will be synchronized with emiss_lod having pirority over
    167 !   mode_emis (see informational messages)
     157! - added "emiss_lod" which serves the same function as "mode_emis" both will be synchronized with
     158!   emiss_lod having pirority over mode_emis (see informational messages)
    168159! - modified existing error message and introduced new informational messages
    169160!    - CM0436 - now also applies to invalid LOD settings
    170161!    - CM0463 - emiss_lod take precedence in case of conflict with mod_emis
    171162!    - CM0464 - emiss_lod valued assigned based on mode_emis if undefined
    172 ! 
     163!
    173164! 3930 2019-04-24 14:57:18Z forkel
    174165! Changed chem_non_transport_physics to chem_non_advective_processes
    175 !
    176 !
     166!
    177167! 3929 2019-04-24 12:52:08Z banzhafs
    178168! Correct/complete module_interface introduction for chemistry model
     
    182172! 3784 2019-03-05 14:16:20Z banzhafs
    183173! 2D output of emission fluxes
    184 ! 
     174!
    185175! 3784 2019-03-05 14:16:20Z banzhafs
    186 ! Bugfix, uncomment erroneous commented variable used for dry deposition. 
    187 ! Bugfix in 3D emission output. 
    188 ! 
     176! Bugfix, uncomment erroneous commented variable used for dry deposition.
     177! Bugfix in 3D emission output.
     178!
    189179! 3784 2019-03-05 14:16:20Z banzhafs
    190 ! Changes related to global restructuring of location messages and introduction 
     180! Changes related to global restructuring of location messages and introduction
    191181! of additional debug messages
    192 ! 
     182!
    193183! 3784 2019-03-05 14:16:20Z banzhafs
    194184! some formatting of the deposition code
     
    196186! 3784 2019-03-05 14:16:20Z banzhafs
    197187! some formatting
    198 ! 
     188!
    199189! 3784 2019-03-05 14:16:20Z banzhafs
    200 ! added cs_mech to USE chem_gasphase_mod 
    201 ! 
     190! added cs_mech to USE chem_gasphase_mod
     191!
    202192! 3784 2019-03-05 14:16:20Z banzhafs
    203193! renamed get_mechanismname to get_mechanism_name
    204194! renamed do_emiss to emissions_anthropogenic and do_depo to deposition_dry (ecc)
    205 ! 
     195!
    206196! 3784 2019-03-05 14:16:20Z banzhafs
    207197! Unused variables removed/taken care of
    208 !
    209198!
    210199! 3784 2019-03-05 14:16:20Z forkel
    211200! Replaced READ from unit 10 by CALL get_mechanismname also in chem_header
    212201!
    213 !
    214202! 3783 2019-03-05 13:23:50Z forkel
    215 ! Removed forgotte write statements an some unused variables (did not touch the
    216 ! parts related to deposition)
    217 !
    218 !
     203! Removed forgotte write statements an some unused variables (did not touch the parts related to
     204! deposition)
     205!
    219206! 3780 2019-03-05 11:19:45Z forkel
    220207! Removed READ from unit 10, added CALL get_mechanismname
    221 !
    222 !
     208!
    223209! 3767 2019-02-27 08:18:02Z raasch
    224210! unused variable for file index removed from rrd-subroutines parameter list
    225 ! 
     211!
    226212! 3738 2019-02-12 17:00:45Z suehring
    227213! Clean-up debug prints
    228 ! 
     214!
    229215! 3737 2019-02-12 16:57:06Z suehring
    230 ! Enable mesoscale offline nesting for chemistry variables as well as
    231 ! initialization of chemistry via dynamic input file.
    232 ! 
     216! Enable mesoscale offline nesting for chemistry variables as well as initialization of chemistry
     217! via dynamic input file.
     218!
    233219! 3719 2019-02-06 13:10:18Z kanani
    234 ! Resolved cpu logpoint overlap with all progn.equations, moved cpu_log call
    235 ! to prognostic_equations for better overview
    236 ! 
     220! Resolved cpu logpoint overlap with all progn.equations, moved cpu_log call to prognostic_equations
     221! for better overview
     222!
    237223! 3700 2019-01-26 17:03:42Z knoop
    238224! Some interface calls moved to module_interface + cleanup
    239 ! 
     225!
    240226! 3664 2019-01-09 14:00:37Z forkel
    241227! Replaced misplaced location message by @todo
    242 !
    243 !
     228!
    244229! 3654 2019-01-07 16:31:57Z suehring
    245230! Disable misplaced location message
    246 ! 
     231!
    247232! 3652 2019-01-07 15:29:59Z forkel
    248233! Checks added for chemistry mechanism, parameter chem_mechanism added (basit)
    249 ! 
     234!
    250235! 2718 2018-01-02 08:49:38Z maronga
    251236! Initial revision
    252237!
    253 ! 
     238!
    254239!
    255240!
     
    263248!
    264249!
    265 !------------------------------------------------------------------------------!
     250!--------------------------------------------------------------------------------------------------!
    266251! Description:
    267252! ------------
     
    269254!> @todo Extend chem_species type by nspec and nvar as addititional elements (RF)
    270255!> @todo Check possibility to reduce dimension of chem_species%conc from nspec to nvar (RF)
    271 !> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not
    272 !>       allowed to use the chemistry model in a precursor run and additionally
    273 !>       not using it in a main run
     256!> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not allowed to use the
     257!>       chemistry model in a precursor run and additionally not using it in a main run
    274258!> @todo Implement turbulent inflow of chem spcs in inflow_turbulence.
    275259!> @todo Separate boundary conditions for each chem spcs to be implemented
    276 !> @todo Currently only total concentration are calculated. Resolved, parameterized
    277 !>       and chemistry fluxes although partially and some completely coded but
    278 !>       are not operational/activated in this version. bK.
    279 !> @todo slight differences in passive scalar and chem spcs when chem reactions
    280 !>       turned off. Need to be fixed. bK
     260!> @todo Currently only total concentration are calculated. Resolved, parameterized and chemistry
     261!>       fluxes although partially and some completely coded but are not operational/activated in
     262!>       this version. bK.
     263!> @todo slight differences in passive scalar and chem spcs when chem reactions turned off. Need to
     264!>       be fixed. bK
    281265!> @todo test nesting for chem spcs, was implemented by suehring (kanani)
    282266!> @todo chemistry error messages
    283 ! 
    284 !------------------------------------------------------------------------------!
     267!
     268!--------------------------------------------------------------------------------------------------!
    285269
    286270 MODULE chemistry_model_mod
     
    313297
    314298    USE control_parameters,                                                                        &
    315          ONLY:  bc_lr_cyc, bc_ns_cyc,                                                              &
     299         ONLY:  air_chemistry,                                                                     &
    316300                bc_dirichlet_l,                                                                    &
    317301                bc_dirichlet_n,                                                                    &
     
    319303                bc_dirichlet_s,                                                                    &
    320304                bc_radiation_l,                                                                    &
     305                bc_lr_cyc,                                                                         &
     306                bc_ns_cyc,                                                                         &
    321307                bc_radiation_n,                                                                    &
    322308                bc_radiation_r,                                                                    &
    323309                bc_radiation_s,                                                                    &
    324310                debug_output,                                                                      &
    325                 dt_3d, humidity, initializing_actions, message_string,                             &
    326                 omega, tsc, intermediate_timestep_count, intermediate_timestep_count_max,          &
     311                dt_3d,                                                                             &
     312                humidity,                                                                          &
     313                initializing_actions,                                                              &
     314                intermediate_timestep_count,                                                       &
     315                intermediate_timestep_count_max,                                                   &
    327316                max_pr_user,                                                                       &
     317                message_string,                                                                    &
    328318                monotonic_limiter_z,                                                               &
     319                omega,                                                                             &
    329320                restart_data_format_output,                                                        &
    330321                scalar_advec,                                                                      &
    331                 timestep_scheme, use_prescribed_profile_data, ws_scheme_sca, air_chemistry
     322                timestep_scheme,                                                                   &
     323                tsc,                                                                               &
     324                use_prescribed_profile_data,                                                       &
     325                ws_scheme_sca
    332326
    333327    USE arrays_3d,                                                                                 &
     
    359353    SAVE
    360354
     355    INTEGER, DIMENSION(nkppctrl)                           ::  icntrl       !< 20 integer parameters for fine tuning KPP code
     356
     357    REAL(kind=wp), PUBLIC ::  cs_time_step = 0._wp
     358
     359    REAL(kind=wp), DIMENSION(nkppctrl)                     ::  rcntrl       !< 20 real parameters for fine tuning of KPP code
     360                                                                            !< (e.g starting internal timestep of solver)
     361
    361362    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_1  !< pointer for swapping of timelevels for conc
    362363    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_2  !< pointer for swapping of timelevels for conc
    363364    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_3  !< pointer for swapping of timelevels for conc
    364     REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  freq_1       !< pointer for phtolysis frequncies 
     365    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  freq_1       !< pointer for phtolysis frequncies
    365366                                                                            !< (only 1 timelevel required)
    366     INTEGER, DIMENSION(nkppctrl)                           ::  icntrl       !< 20 integer parameters for fine tuning KPP code
    367367                                                                            !< (e.g. solver type)
    368     REAL(kind=wp), DIMENSION(nkppctrl)                     ::  rcntrl       !< 20 real parameters for fine tuning of KPP code
    369                                                                             !< (e.g starting internal timestep of solver)
    370 
    371     REAL(kind=wp), PUBLIC ::  cs_time_step = 0._wp
     368
    372369
    373370!
    374371!-- Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010)
    375     !
    376     INTEGER(iwp), PARAMETER ::  nlu_dep = 15                   !< Number of DEPAC landuse classes (lu's)
    377     INTEGER(iwp), PARAMETER ::  ncmp = 10                      !< Number of DEPAC gas components
    378     INTEGER(iwp), PARAMETER ::  nposp = 69                     !< Number of possible species for deposition
    379 !
    380 !-- DEPAC landuse classes as defined in LOTOS-EUROS model v2.1                             
    381     INTEGER(iwp) ::  ilu_grass              = 1       
    382     INTEGER(iwp) ::  ilu_arable             = 2       
    383     INTEGER(iwp) ::  ilu_permanent_crops    = 3         
    384     INTEGER(iwp) ::  ilu_coniferous_forest  = 4         
    385     INTEGER(iwp) ::  ilu_deciduous_forest   = 5         
    386     INTEGER(iwp) ::  ilu_water_sea          = 6       
    387     INTEGER(iwp) ::  ilu_urban              = 7       
    388     INTEGER(iwp) ::  ilu_other              = 8 
    389     INTEGER(iwp) ::  ilu_desert             = 9 
    390     INTEGER(iwp) ::  ilu_ice                = 10
    391     INTEGER(iwp) ::  ilu_savanna            = 11
    392     INTEGER(iwp) ::  ilu_tropical_forest    = 12
    393     INTEGER(iwp) ::  ilu_water_inland       = 13
    394     INTEGER(iwp) ::  ilu_mediterrean_scrub  = 14
    395     INTEGER(iwp) ::  ilu_semi_natural_veg   = 15
     372    INTEGER(iwp), PARAMETER ::  nlu_dep = 15            !< Number of DEPAC landuse classes (lu's)
     373    INTEGER(iwp), PARAMETER ::  ncmp = 10               !< Number of DEPAC gas components
     374    INTEGER(iwp), PARAMETER ::  nposp = 69              !< Number of possible species for deposition
     375!
     376!-- DEPAC landuse classes as defined in LOTOS-EUROS model v2.1
     377    INTEGER(iwp) ::  ilu_grass              = 1
     378    INTEGER(iwp) ::  ilu_arable             = 2
     379    INTEGER(iwp) ::  ilu_permanent_crops    = 3
     380    INTEGER(iwp) ::  ilu_coniferous_forest  = 4
     381    INTEGER(iwp) ::  ilu_deciduous_forest   = 5
     382    INTEGER(iwp) ::  ilu_water_sea          = 6
     383    INTEGER(iwp) ::  ilu_urban              = 7
     384    INTEGER(iwp) ::  ilu_other              = 8
     385    INTEGER(iwp) ::  ilu_desert             = 9
     386    INTEGER(iwp) ::  ilu_ice                = 10
     387    INTEGER(iwp) ::  ilu_savanna            = 11
     388    INTEGER(iwp) ::  ilu_tropical_forest    = 12
     389    INTEGER(iwp) ::  ilu_water_inland       = 13
     390    INTEGER(iwp) ::  ilu_mediterrean_scrub  = 14
     391    INTEGER(iwp) ::  ilu_semi_natural_veg   = 15
    396392
    397393!
     
    405401!
    406402!-- Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) to W m-2
    407     REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  alpha   =(/ 0.009, 0.009, 0.009, 0.006, 0.006, -999., -999., 0.009, -999.,      &
    408          -999., 0.009, 0.006, -999., 0.009, 0.008/)*4.57
     403    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  alpha   = (/ 0.009, 0.009, 0.009, 0.006, 0.006,    &
     404                    -999.0, -999., 0.009, -999.0, -999.0, 0.009, 0.006, -999.0, 0.009, 0.008 /)*4.57
    409405!
    410406!-- Set temperatures per land use for f_temp
    411     REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmin = (/ 12.0, 12.0,  12.0,  0.0,  0.0, -999., -999., 12.0, -999., -999.,      &
    412          12.0,  0.0, -999., 12.0,  8.0/)
    413     REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  topt = (/ 26.0, 26.0,  26.0, 18.0, 20.0, -999., -999., 26.0, -999., -999.,      &
    414          26.0, 20.0, -999., 26.0, 24.0 /)
    415     REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmax = (/ 40.0, 40.0,  40.0, 36.0, 35.0, -999., -999., 40.0, -999., -999.,      &
    416          40.0, 35.0, -999., 40.0, 39.0 /)
     407    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmin = (/ 12.0,  12.0,  12.0,  0.0,  0.0, -999.0,  &
     408                                      -999.0, 12.0, -999.0, -999.0, 12.0,  0.0, -999.0, 12.0,  8.0/)
     409    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  topt = (/ 26.0, 26.0,  26.0, 18.0, 20.0, -999.0,   &
     410                                     -999.0, 26.0, -999.0, -999.0, 26.0, 20.0, -999.0, 26.0, 24.0 /)
     411    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmax = (/ 40.0, 40.0,  40.0, 36.0, 35.0, -999.0,   &
     412                                    -999.0, 40.0, -999.0, -999.0,  40.0, 35.0, -999.0, 40.0, 39.0 /)
    417413!
    418414!-- Set f_min:
    419     REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  f_min = (/ 0.01, 0.01, 0.01, 0.1, 0.1, -999., -999., 0.01, -999., -999., 0.01,  &
    420          0.1, -999., 0.01, 0.04/)
     415    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  f_min = (/ 0.01, 0.01, 0.01, 0.1, 0.1, -999.0,     &
     416                                       -999.0, 0.01, -999.0, -999.0, 0.01, 0.1, -999.0, 0.01, 0.04/)
    421417
    422418!
    423419!-- Set maximal conductance (m/s)
    424420!-- (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from  mmol O3/m2/s to m/s
    425     REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  g_max = (/ 270., 300., 300., 140., 150., -999., -999., 270., -999., -999.,      &
    426          270., 150., -999., 300., 422./)/41000.
     421    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  g_max = (/ 270.0, 300.0, 300.0, 140.0, 150.0,      &
     422              -999.0, -999.0, 270.0, -999.0, -999.0, 270.0, 150.0, -999.0, 300.0, 422.0 /) / 41000.0
    427423!
    428424!-- Set max, min for vapour pressure deficit vpd
    429     REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_max = (/1.3, 0.9, 0.9, 0.5, 1.0, -999., -999., 1.3, -999., -999., 1.3,      &
    430          1.0, -999., 0.9, 2.8/)
    431     REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_min = (/3.0, 2.8, 2.8, 3.0, 3.25, -999., -999., 3.0, -999., -999., 3.0,     &
    432          3.25, -999., 2.8, 4.5/)
     425    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_max = (/ 1.3, 0.9, 0.9, 0.5, 1.0, -999.0,      &
     426                                           -999.0, 1.3, -999.0, -999.0, 1.3, 1.0, -999.0, 0.9, 2.8/)
     427    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_min = (/ 3.0, 2.8, 2.8, 3.0, 3.25, -999.0,     &
     428                                         -999.0, 3.0, -999.0, -999.0, 3.0,  3.25, -999.0, 2.8, 4.5/)
    433429
    434430    PUBLIC nreact
     
    436432    PUBLIC nvar                !< number of variable gas phase chemical species (nvar <= nspec)
    437433    PUBLIC spc_names           !< names of gas phase chemical species (come from KPP) (come from KPP)
    438     PUBLIC spec_conc_2 
    439 !   
     434    PUBLIC spec_conc_2
     435!
    440436!-- Interface section
    441437    INTERFACE chem_3d_data_averaging
    442        MODULE PROCEDURE chem_3d_data_averaging 
     438       MODULE PROCEDURE chem_3d_data_averaging
    443439    END INTERFACE chem_3d_data_averaging
    444440
     
    508504       MODULE PROCEDURE chem_non_advective_processes_ij
    509505    END INTERFACE chem_non_advective_processes
    510    
     506
    511507    INTERFACE chem_exchange_horiz_bounds
    512508       MODULE PROCEDURE chem_exchange_horiz_bounds
    513     END INTERFACE chem_exchange_horiz_bounds   
     509    END INTERFACE chem_exchange_horiz_bounds
    514510
    515511    INTERFACE chem_prognostic_equations
     
    532528
    533529    INTERFACE chem_wrd_local
    534        MODULE PROCEDURE chem_wrd_local 
     530       MODULE PROCEDURE chem_wrd_local
    535531    END INTERFACE chem_wrd_local
    536532
    537533    INTERFACE chem_depo
    538        MODULE PROCEDURE chem_depo 
     534       MODULE PROCEDURE chem_depo
    539535    END INTERFACE chem_depo
    540536
    541537    INTERFACE drydepos_gas_depac
    542        MODULE PROCEDURE drydepos_gas_depac 
     538       MODULE PROCEDURE drydepos_gas_depac
    543539    END INTERFACE drydepos_gas_depac
    544540
    545541    INTERFACE rc_special
    546        MODULE PROCEDURE rc_special 
     542       MODULE PROCEDURE rc_special
    547543    END INTERFACE rc_special
    548544
    549545    INTERFACE  rc_gw
    550        MODULE PROCEDURE rc_gw 
     546       MODULE PROCEDURE rc_gw
    551547    END INTERFACE rc_gw
    552548
    553     INTERFACE rw_so2 
    554        MODULE PROCEDURE rw_so2 
     549    INTERFACE rw_so2
     550       MODULE PROCEDURE rw_so2
    555551    END INTERFACE rw_so2
    556552
    557553    INTERFACE rw_nh3_sutton
    558        MODULE PROCEDURE rw_nh3_sutton 
     554       MODULE PROCEDURE rw_nh3_sutton
    559555    END INTERFACE rw_nh3_sutton
    560556
    561557    INTERFACE rw_constant
    562        MODULE PROCEDURE rw_constant 
     558       MODULE PROCEDURE rw_constant
    563559    END INTERFACE rw_constant
    564560
    565561    INTERFACE rc_gstom
    566        MODULE PROCEDURE rc_gstom 
     562       MODULE PROCEDURE rc_gstom
    567563    END INTERFACE rc_gstom
    568564
    569565    INTERFACE rc_gstom_emb
    570        MODULE PROCEDURE rc_gstom_emb 
     566       MODULE PROCEDURE rc_gstom_emb
    571567    END INTERFACE rc_gstom_emb
    572568
    573569    INTERFACE par_dir_diff
    574        MODULE PROCEDURE par_dir_diff 
     570       MODULE PROCEDURE par_dir_diff
    575571    END INTERFACE par_dir_diff
    576572
    577573    INTERFACE rc_get_vpd
    578        MODULE PROCEDURE rc_get_vpd 
     574       MODULE PROCEDURE rc_get_vpd
    579575    END INTERFACE rc_get_vpd
    580576
    581577    INTERFACE rc_gsoil_eff
    582        MODULE PROCEDURE rc_gsoil_eff 
     578       MODULE PROCEDURE rc_gsoil_eff
    583579    END INTERFACE rc_gsoil_eff
    584580
    585581    INTERFACE rc_rinc
    586        MODULE PROCEDURE rc_rinc 
     582       MODULE PROCEDURE rc_rinc
    587583    END INTERFACE rc_rinc
    588584
    589585    INTERFACE rc_rctot
    590        MODULE PROCEDURE rc_rctot 
     586       MODULE PROCEDURE rc_rctot
    591587    END INTERFACE rc_rctot
    592588
    593589    INTERFACE drydepo_aero_zhang_vd
    594        MODULE PROCEDURE drydepo_aero_zhang_vd 
     590       MODULE PROCEDURE drydepo_aero_zhang_vd
    595591    END INTERFACE drydepo_aero_zhang_vd
    596592
     
    612608
    613609
    614 !------------------------------------------------------------------------------!
     610!--------------------------------------------------------------------------------------------------!
    615611! Description:
    616612! ------------
    617 !> Subroutine for averaging 3D data of chemical species. Due to the fact that
    618 !> the averaged chem arrays are allocated in chem_init, no if-query concerning
    619 !> the allocation is required (in any mode). Attention: If you just specify an
    620 !> averaged output quantity in the _p3dr file during restarts the first output
    621 !> includes the time between the beginning of the restart run and the first
    622 !> output time (not necessarily the whole averaging_interval you have
    623 !> specified in your _p3d/_p3dr file )
    624 !------------------------------------------------------------------------------!
     613!> Subroutine for averaging 3D data of chemical species. Due to the fact that the averaged chem
     614!> arrays are allocated in chem_init, no if-query concerning the allocation is required (in any
     615!> mode). Attention: If you just specify an averaged output quantity in the _p3dr file during
     616!> restarts the first output includes the time between the beginning of the restart run and the
     617!> first output time (not necessarily the whole averaging_interval you have specified in your
     618!> _p3d/_p3dr file ).
     619!--------------------------------------------------------------------------------------------------!
    625620 SUBROUTINE chem_3d_data_averaging( mode, variable )
    626621
    627622    USE control_parameters
    628623
    629     USE exchange_horiz_mod,                                                    &
     624    USE exchange_horiz_mod,                                                                        &
    630625        ONLY:  exchange_horiz_2d
    631626
    632627
    633     CHARACTER (LEN=*) ::  mode     !< 
    634     CHARACTER (LEN=*) ::  variable !< 
     628    CHARACTER (LEN=*) ::  mode     !<
     629    CHARACTER (LEN=*) ::  variable !<
    635630
    636631    LOGICAL ::  match_def  !< flag indicating default-type surface
     
    641636    INTEGER(iwp) ::  j                  !< grid index y direction
    642637    INTEGER(iwp) ::  k                  !< grid index z direction
     638    INTEGER(iwp) ::  lsp                !< running index for chem spcs
    643639    INTEGER(iwp) ::  m                  !< running index surface type
    644     INTEGER(iwp) ::  lsp               !< running index for chem spcs
    645 
    646     IF ( (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_')  )  THEN
     640
     641    IF ( ( variable(1:3) == 'kc_'  .OR.  variable(1:3) == 'em_' )  )  THEN
    647642
    648643       IF ( mode == 'allocate' )  THEN
    649644
    650645          DO  lsp = 1, nspec
    651              IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
    652                   TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
     646             IF ( TRIM( variable(1:3) ) == 'kc_'  .AND.                                            &
     647                  TRIM( variable(4:) )  == TRIM( chem_species(lsp)%name ) )  THEN
    653648                IF ( .NOT. ALLOCATED( chem_species(lsp)%conc_av ) )  THEN
    654649                   ALLOCATE( chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     
    661656
    662657          DO  lsp = 1, nspec
    663              IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
    664                   TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
     658             IF ( TRIM( variable(1:3) ) == 'kc_'  .AND.                                            &
     659                  TRIM( variable(4:) )  == TRIM( chem_species(lsp)%name ) )  THEN
    665660                DO  i = nxlg, nxrg
    666661                   DO  j = nysg, nyng
    667662                      DO  k = nzb, nzt+1
    668                          chem_species(lsp)%conc_av(k,j,i) =                              &
    669                                            chem_species(lsp)%conc_av(k,j,i) +            &
    670                                            chem_species(lsp)%conc(k,j,i)
     663                         chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) +     &
     664                                                            chem_species(lsp)%conc(k,j,i)
    671665                      ENDDO
    672666                   ENDDO
     
    675669                DO  i = nxl, nxr
    676670                   DO  j = nys, nyn
    677                       match_def = surf_def_h(0)%start_index(j,i) <=                      &
    678                            surf_def_h(0)%end_index(j,i)
    679                       match_lsm = surf_lsm_h%start_index(j,i) <=                         &
    680                            surf_lsm_h%end_index(j,i)
    681                       match_usm = surf_usm_h%start_index(j,i) <=                         &
    682                            surf_usm_h%end_index(j,i)
     671                      match_def = surf_def_h(0)%start_index(j,i) <= surf_def_h(0)%end_index(j,i)
     672                      match_lsm = surf_lsm_h%start_index(j,i)    <= surf_lsm_h%end_index(j,i)
     673                      match_usm = surf_usm_h%start_index(j,i)    <= surf_usm_h%end_index(j,i)
    683674
    684675                      IF ( match_def )  THEN
    685676                         m = surf_def_h(0)%end_index(j,i)
    686                          chem_species(lsp)%cssws_av(j,i) =                               &
    687                               chem_species(lsp)%cssws_av(j,i) + &
    688                               surf_def_h(0)%cssws(lsp,m)
     677                         chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) +       &
     678                                                           surf_def_h(0)%cssws(lsp,m)
    689679                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
    690680                         m = surf_lsm_h%end_index(j,i)
    691                          chem_species(lsp)%cssws_av(j,i) =                               &
    692                               chem_species(lsp)%cssws_av(j,i) + &
    693                               surf_lsm_h%cssws(lsp,m)
     681                         chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) +       &
     682                                                           surf_lsm_h%cssws(lsp,m)
    694683                      ELSEIF ( match_usm )  THEN
    695684                         m = surf_usm_h%end_index(j,i)
    696                          chem_species(lsp)%cssws_av(j,i) =                               &
    697                               chem_species(lsp)%cssws_av(j,i) + &
    698                               surf_usm_h%cssws(lsp,m)
     685                         chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) +       &
     686                                                           surf_usm_h%cssws(lsp,m)
    699687                      ENDIF
    700688                   ENDDO
     
    706694
    707695          DO  lsp = 1, nspec
    708              IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
    709                   TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
     696             IF ( TRIM( variable(1:3) ) == 'kc_'  .AND.                                            &
     697                  TRIM( variable(4:) )  == TRIM( chem_species(lsp)%name ) )  THEN
    710698                DO  i = nxlg, nxrg
    711699                   DO  j = nysg, nyng
    712700                      DO  k = nzb, nzt+1
    713                          chem_species(lsp)%conc_av(k,j,i) =                              &
    714                              chem_species(lsp)%conc_av(k,j,i) /                          &
    715                              REAL( average_count_3d, KIND=wp )
     701                         chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) /     &
     702                                                            REAL( average_count_3d, KIND = wp )
    716703                      ENDDO
    717704                   ENDDO
     
    721708                DO  i = nxlg, nxrg
    722709                   DO  j = nysg, nyng
    723                       chem_species(lsp)%cssws_av(j,i) =                                  &
    724                       chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
     710                      chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) /          &
     711                                                        REAL( average_count_3d, KIND = wp )
    725712                   ENDDO
    726713                ENDDO
     
    734721 END SUBROUTINE chem_3d_data_averaging
    735722
    736    
    737 !------------------------------------------------------------------------------!
     723
     724!--------------------------------------------------------------------------------------------------!
    738725! Description:
    739726! ------------
    740727!> Subroutine to initialize and set all boundary conditions for chemical species
    741 !------------------------------------------------------------------------------!
     728!--------------------------------------------------------------------------------------------------!
    742729 SUBROUTINE chem_boundary_conditions( horizontal_conditions_only )
    743730
    744     USE arrays_3d,                                                             &
     731    USE arrays_3d,                                                                                 &
    745732        ONLY:  dzu
    746733
    747     USE surface_mod,                                                           &
    748         ONLY:  bc_h                                                           
     734    USE surface_mod,                                                                               &
     735        ONLY:  bc_h
    749736
    750737    INTEGER(iwp) ::  i                            !< grid index x direction.
     
    752739    INTEGER(iwp) ::  k                            !< grid index z direction.
    753740    INTEGER(iwp) ::  l                            !< running index boundary type, for up- and downward-facing walls.
     741    INTEGER(iwp) ::  lsp                          !< running index for chem spcs.
    754742    INTEGER(iwp) ::  m                            !< running index surface elements.
    755     INTEGER(iwp) ::  lsp                          !< running index for chem spcs.
    756743
    757744    LOGICAL, OPTIONAL ::  horizontal_conditions_only  !< switch to set horizontal bc only
     
    770757                    j = bc_h(l)%j(m)
    771758                    k = bc_h(l)%k(m)
    772                    chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
    773                                       chem_species(lsp)%conc(k+bc_h(l)%koff,j,i)
     759                   chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =                                  &
     760                                                          chem_species(lsp)%conc(k+bc_h(l)%koff,j,i)
    774761                ENDDO
    775762             ENDDO
     
    785772                   j = bc_h(l)%j(m)
    786773                   k = bc_h(l)%k(m)
    787                    chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
    788                                          chem_species(lsp)%conc_p(k,j,i)
     774                   chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) = chem_species(lsp)%conc_p(k,j,i)
    789775                ENDDO
    790776             ENDDO
    791777          ENDIF
    792778
    793        ENDDO    ! end lsp loop 
     779       ENDDO    ! end lsp loop
    794780
    795781!
     
    798784!>           already runs over all species? (Siggi)
    799785       DO  lsp = 1, nspec
    800           IF ( ibc_cs_t == 0 )  THEN                   
     786          IF ( ibc_cs_t == 0 )  THEN
    801787             chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)
    802788          ELSEIF ( ibc_cs_t == 1 )  THEN
     
    904890!
    905891!--    Lateral Neumann booundary conditions for timelevel t.
    906 !--    This branch is executed when routine is called after the non-advective processes /
    907 !--    before the prognostic equations.
     892!--    This branch is executed when routine is called after the non-advective processes / before the
     893!--    prognostic equations.
    908894       IF ( bc_radiation_cs_s )  THEN
    909895          DO  lsp = 1, nspec
     
    998984
    999985
    1000 !------------------------------------------------------------------------------!
     986!--------------------------------------------------------------------------------------------------!
    1001987! Description:
    1002988! ------------
    1003989!> Subroutine for checking data output for chemical species
    1004 !------------------------------------------------------------------------------!
     990!--------------------------------------------------------------------------------------------------!
    1005991 SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
    1006992
     
    1010996
    1011997    INTEGER(iwp) ::  i
     998    INTEGER(iwp) ::  ilen
    1012999    INTEGER(iwp) ::  lsp
    1013     INTEGER(iwp) ::  ilen
    10141000    INTEGER(iwp) ::  k
    10151001
     
    10261012    IF ( TRIM( var(1:3) ) == 'em_' )  THEN
    10271013       DO  lsp=1,nspec
    1028           IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
     1014          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
    10291015             unit = 'mol m-2 s-1'
    10301016          ENDIF
    10311017!
    1032 !--       It is possible to plant PM10 and PM25 into the gasphase chemistry code
    1033 !--       as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
    1034 !--       set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
     1018!--       It is possible to plant PM10 and PM25 into the gasphase chemistry code as passive species
     1019!--       (e.g. 'passive' in GASPHASE_PREPROC/mechanisms): set unit to micrograms per m**3 for PM10
     1020!--       and PM25 (PM2.5)
    10351021          IF (spec_name(1:2) == 'PM')  THEN
    10361022             unit = 'kg m-2 s-1'
     
    10411027
    10421028       DO  lsp=1,nspec
    1043           IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
     1029          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
    10441030             unit = 'ppm'
    10451031          ENDIF
    10461032!
    1047 !--            It is possible  to plant PM10 and PM25 into the gasphase chemistry code
    1048 !--            as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
    1049 !--            set unit to kilograms per m**3 for PM10 and PM25 (PM2.5)
    1050           IF (spec_name(1:2) == 'PM')  THEN 
     1033!--       It is possible to plant PM10 and PM25 into the gasphase chemistry code as passive species
     1034!--       (e.g. 'passive' in GASPHASE_PREPROC/mechanisms): set unit to kilograms per m**3 for PM10
     1035!--       and PM25 (PM2.5)
     1036          IF (spec_name(1:2) == 'PM')  THEN
    10511037            unit = 'kg m-3'
    10521038          ENDIF
     
    10541040
    10551041       DO  lsp=1,nphot
    1056           IF (TRIM( spec_name ) == TRIM( phot_frequen(lsp)%name ) )  THEN
     1042          IF ( TRIM( spec_name ) == TRIM( phot_frequen(lsp)%name ) )  THEN
    10571043             unit = 'sec-1'
    10581044          ENDIF
     
    10651051
    10661052
    1067 !------------------------------------------------------------------------------!
     1053!--------------------------------------------------------------------------------------------------!
    10681054! Description:
    10691055! ------------
    10701056!> Subroutine for checking data output of profiles for chemistry model
    1071 !------------------------------------------------------------------------------!
     1057!--------------------------------------------------------------------------------------------------!
    10721058 SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
    10731059
    10741060    USE arrays_3d
    10751061
    1076     USE control_parameters,                                                    &
     1062    USE control_parameters,                                                                        &
    10771063        ONLY:  data_output_pr, message_string
    10781064
     
    10821068
    10831069
    1084     CHARACTER (LEN=*) ::  unit     !<
    1085     CHARACTER (LEN=*) ::  variable !<
    1086     CHARACTER (LEN=*) ::  dopr_unit
     1070    CHARACTER (LEN=*)  ::  unit     !<
     1071    CHARACTER (LEN=*)  ::  variable !<
     1072    CHARACTER (LEN=*)  ::  dopr_unit
    10871073    CHARACTER (LEN=16) ::  spec_name
    10881074
     
    10901076
    10911077
    1092     spec_name = TRIM( variable(4:) )   
    1093 
    1094        IF (  .NOT.  air_chemistry )  THEN
    1095           message_string = 'data_output_pr = ' //                        &
    1096           TRIM( data_output_pr(var_count) ) // ' is not imp' // &
    1097                       'lemented for air_chemistry = .FALSE.'
    1098           CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )             
     1078    spec_name = TRIM( variable(4:) )
     1079
     1080       IF ( .NOT. air_chemistry )  THEN
     1081          message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //             &
     1082                           ' is not imp' // 'lemented for air_chemistry = .FALSE.'
     1083          CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )
    10991084
    11001085       ELSE
    11011086          DO  lsp = 1, nspec
    1102              IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
     1087             IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
    11031088                cs_pr_count = cs_pr_count+1
    11041089                cs_pr_index(cs_pr_count) = lsp
    1105                 dopr_index(var_count) = pr_palm+cs_pr_count 
     1090                dopr_index(var_count) = pr_palm+cs_pr_count
    11061091                dopr_unit  = 'ppm'
    1107                 IF (spec_name(1:2) == 'PM')  THEN
     1092                IF ( spec_name(1:2) == 'PM')  THEN
    11081093                     dopr_unit = 'kg m-3'
    11091094                ENDIF
    11101095                   hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
    1111                    unit = dopr_unit 
     1096                   unit = dopr_unit
    11121097             ENDIF
    11131098          ENDDO
     
    11171102
    11181103
    1119 !------------------------------------------------------------------------------!
     1104!--------------------------------------------------------------------------------------------------!
    11201105! Description:
    11211106! ------------
    11221107!> Check parameters routine for chemistry_model_mod
    1123 !------------------------------------------------------------------------------!
     1108!--------------------------------------------------------------------------------------------------!
    11241109 SUBROUTINE chem_check_parameters
    11251110
     
    11271112        ONLY:  bc_lr, bc_ns
    11281113
     1114    INTEGER (iwp) ::  lsp          !< running index for chem spcs.
     1115    INTEGER (iwp) ::  lsp_usr      !< running index for user defined chem spcs
     1116
    11291117    LOGICAL  ::  found
    1130     INTEGER (iwp) ::  lsp_usr      !< running index for user defined chem spcs
    1131     INTEGER (iwp) ::  lsp          !< running index for chem spcs.
    1132 !
    1133 !-- check for chemical reactions status
     1118
     1119!
     1120!-- Check for chemical reactions status
    11341121    IF ( chem_gasphase_on )  THEN
    11351122       message_string = 'Chemical reactions: ON'
     
    11401127    ENDIF
    11411128!
    1142 !-- check for chemistry time-step
     1129!-- Check for chemistry time-step
    11431130    IF ( call_chem_at_all_substeps )  THEN
    11441131       message_string = 'Chemistry is calculated at all meteorology time-step'
     
    11491136    ENDIF
    11501137!
    1151 !-- check for photolysis scheme
     1138!-- Check for photolysis scheme
    11521139    IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant')  )  THEN
    11531140       message_string = 'Incorrect photolysis scheme selected, please check spelling'
     
    11551142    ENDIF
    11561143!
    1157 !-- check for chemical mechanism used
     1144!-- Check for chemical mechanism used
    11581145    CALL get_mechanism_name
    11591146    IF ( chem_mechanism /= TRIM( cs_mech ) )  THEN
    1160        message_string = 'Incorrect chemistry mechanism selected, check spelling in namelist and/or chem_gasphase_mod'
     1147       message_string =                                                                            &
     1148       'Incorrect chemistry mechanism selected, check spelling in namelist and/or chem_gasphase_mod'
    11611149       CALL message( 'chem_check_parameters', 'CM0462', 1, 2, 0, 6, 0 )
    11621150    ENDIF
     
    12521240!
    12531241!-- Cyclic conditions must be set identically at opposing boundaries
    1254     IF ( ( bc_cs_l == 'cyclic' .AND. bc_cs_r /= 'cyclic' )  .OR.                                 &
     1242    IF ( ( bc_cs_l == 'cyclic' .AND. bc_cs_r /= 'cyclic' )  .OR.                                   &
    12551243         ( bc_cs_r == 'cyclic' .AND. bc_cs_l /= 'cyclic' ) )  THEN
    12561244       message_string = 'boundary conditions bc_cs_l and bc_cs_r must both be cyclic or non-cyclic'
    12571245       CALL message( 'chem_check_parameters','PA0716', 1, 2, 0, 6, 0 )
    12581246    ENDIF
    1259     IF ( ( bc_cs_n == 'cyclic' .AND. bc_cs_s /= 'cyclic' )  .OR.                                 &
     1247    IF ( ( bc_cs_n == 'cyclic' .AND. bc_cs_s /= 'cyclic' )  .OR.                                   &
    12601248         ( bc_cs_s == 'cyclic' .AND. bc_cs_n /= 'cyclic' ) )  THEN
    12611249       message_string = 'boundary conditions bc_cs_n and bc_cs_s must both be cyclic or non-cyclic'
     
    12981286    CALL chem_init_internal
    12991287!
    1300 !-- check for initial chem species input
     1288!-- Check for initial chem species input
    13011289    lsp_usr = 1
    13021290    lsp     = 1
     
    13071295             found = .TRUE.
    13081296             EXIT
    1309           ENDIF 
     1297          ENDIF
    13101298       ENDDO
    13111299       IF ( .NOT.  found )  THEN
     
    13171305    ENDDO
    13181306!
    1319 !-- check for surface  emission flux chem species
     1307!-- Check for surface  emission flux chem species
    13201308    lsp_usr = 1
    13211309    lsp     = 1
     
    13261314             found = .TRUE.
    13271315             EXIT
    1328           ENDIF 
     1316          ENDIF
    13291317       ENDDO
    13301318       IF ( .NOT.  found )  THEN
     
    13391327
    13401328
    1341 !------------------------------------------------------------------------------!
     1329!--------------------------------------------------------------------------------------------------!
    13421330! Description:
    13431331! ------------
    13441332!> Subroutine defining 2D output variables for chemical species
    13451333!> @todo: Remove "mode" from argument list, not used.
    1346 !------------------------------------------------------------------------------!
    1347  SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf,   &
    1348                                   two_d, nzb_do, nzt_do, fill_value )
     1334!--------------------------------------------------------------------------------------------------!
     1335 SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do, &
     1336                                 fill_value )
    13491337
    13501338
    13511339    CHARACTER (LEN=*) ::  grid       !<
    1352     CHARACTER (LEN=*) ::  mode       !<
    1353     CHARACTER (LEN=*) ::  variable   !<
    1354     INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or time-averaged data
     1340    CHARACTER (LEN=*) ::  mode       !<
     1341    CHARACTER (LEN=*) ::  variable   !<
     1342
     1343    INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or
     1344                                     !< time-averaged data
    13551345    INTEGER(iwp) ::  nzb_do          !< lower limit of the domain (usually nzb)
    13561346    INTEGER(iwp) ::  nzt_do          !< upper limit of the domain (usually nzt+1)
    1357     LOGICAL      ::  found           !<
    1358     LOGICAL      ::  two_d           !< flag parameter that indicates 2D variables (horizontal cross sections)
    1359     REAL(wp)     ::  fill_value
    1360     REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf 
     1347
     1348    LOGICAL      ::  found           !<
     1349    LOGICAL      ::  two_d           !< flag parameter that indicates 2D variables (horizontal cross
     1350                                     !< sections)
     1351
     1352    REAL(wp)     ::  fill_value
     1353
     1354    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf
    13611355
    13621356!
     
    13781372    spec_name = TRIM( variable(4:char_len-3) )
    13791373!
    1380 !-- Output of emission values, i.e. surface fluxes cssws. 
     1374!-- Output of emission values, i.e. surface fluxes cssws.
    13811375    IF ( variable(1:3) == 'em_' )  THEN
    13821376
     
    13861380          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
    13871381!
    1388 !--          No average output for now. 
     1382!--          No average output for now.
    13891383             DO  m = 1, surf_lsm_h%ns
    1390                 local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1) =              &
    1391                               local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1) &
    1392                                             + surf_lsm_h%cssws(lsp,m)
     1384                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1) =                                  &
     1385                                                   local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1) &
     1386                                                   + surf_lsm_h%cssws(lsp,m)
    13931387             ENDDO
    13941388             DO  m = 1, surf_usm_h%ns
    1395                    local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1) =           &
    1396                               local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1) &
    1397                                             + surf_usm_h%cssws(lsp,m)
     1389                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1) =                               &
     1390                                                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1) &
     1391                                                   + surf_usm_h%cssws(lsp,m)
    13981392             ENDDO
    13991393             grid = 'zu'
     
    14051399
    14061400       DO  lsp=1,nspec
    1407           IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name )  .AND.                           &
     1401          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name )  .AND.                          &
    14081402                ( (variable(char_len-2:) == '_xy')  .OR.                                           &
    14091403                  (variable(char_len-2:) == '_xz')  .OR.                                           &
    1410                   (variable(char_len-2:) == '_yz') ) )  THEN             
     1404                  (variable(char_len-2:) == '_yz') ) )  THEN
    14111405!
    14121406!--   todo: remove or replace by "CALL message" mechanism (kanani)
    14131407!                    IF(myid == 0)  WRITE(6,*) 'Output of species ' // TRIM( variable )  //       &
    1414 !                                                             TRIM( chem_species(lsp)%name )       
     1408!                                                             TRIM( chem_species(lsp)%name )
    14151409             IF (av == 0)  THEN
    14161410                DO  i = nxl, nxr
     
    14371431                ENDDO
    14381432             ENDIF
    1439              grid = 'zu'           
     1433             grid = 'zu'
    14401434             found = .TRUE.
    14411435          ENDIF
     
    14451439    RETURN
    14461440
    1447  END SUBROUTINE chem_data_output_2d     
    1448 
    1449 
    1450 !------------------------------------------------------------------------------!
     1441 END SUBROUTINE chem_data_output_2d
     1442
     1443
     1444!--------------------------------------------------------------------------------------------------!
    14511445! Description:
    14521446! ------------
    14531447!> Subroutine defining 3D output variables for chemical species
    1454 !------------------------------------------------------------------------------!
     1448!--------------------------------------------------------------------------------------------------!
    14551449 SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
    14561450
     
    14591453
    14601454    CHARACTER (LEN=*)    ::  variable     !<
    1461     INTEGER(iwp)         ::  av           !<
    1462     INTEGER(iwp) ::  nzb_do               !< lower limit of the data output (usually 0)
    1463     INTEGER(iwp) ::  nzt_do               !< vertical upper limit of the data output (usually nz_do3d)
     1455
     1456    INTEGER(iwp)         ::  av         !<
     1457    INTEGER(iwp) ::  nzb_do             !< lower limit of the data output (usually 0)
     1458    INTEGER(iwp) ::  nzt_do             !< vertical upper limit of the data output (usually nz_do3d)
    14641459
    14651460    LOGICAL      ::  found                !<
    14661461
    14671462    REAL(wp)             ::  fill_value   !<
    1468     REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf
    1469 !
    1470 !-- local variables
     1463
     1464    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf
     1465!
     1466!-- Local variables
    14711467    CHARACTER(LEN=16)    ::  spec_name
     1468
    14721469    INTEGER(iwp)         ::  i
    14731470    INTEGER(iwp)         ::  j
     
    14891486       DO  lsp = 1, nvar   !!! cssws - nvar species, chem_species - nspec species !!!
    14901487          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
    1491          
     1488
    14921489             local_pf = 0.0_wp
    14931490!
    14941491!--          no average for now
    14951492             DO  m = 1, surf_usm_h%ns
    1496                 local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = &
    1497                    local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m)
     1493                local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) =                        &
     1494                                         local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) &
     1495                                         + surf_usm_h%cssws(lsp,m)
    14981496             ENDDO
    14991497             DO  m = 1, surf_lsm_h%ns
    1500                 local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = &
    1501                   local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) + surf_lsm_h%cssws(lsp,m)
     1498                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) =                        &
     1499                                         local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) &
     1500                                         + surf_lsm_h%cssws(lsp,m)
    15021501             ENDDO
    15031502             DO  l = 0, 3
    15041503                DO  m = 1, surf_usm_v(l)%ns
    1505                    local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = &
    1506                      local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) + surf_usm_v(l)%cssws(lsp,m)
     1504                   local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) =            &
     1505                                local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) &
     1506                                + surf_usm_v(l)%cssws(lsp,m)
    15071507                ENDDO
    15081508                DO  m = 1, surf_lsm_v(l)%ns
    1509                    local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = &
    1510                       local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) + surf_lsm_v(l)%cssws(lsp,m)
     1509                   local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) =            &
     1510                                local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) &
     1511                                + surf_lsm_v(l)%cssws(lsp,m)
    15111512                ENDDO
    15121513             ENDDO
     
    15161517    ELSE
    15171518      DO  lsp = 1, nspec
    1518          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
     1519         IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
    15191520!
    15201521!--   todo: remove or replace by "CALL message" mechanism (kanani)
    15211522!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
    1522 !                                                           TRIM( chem_species(lsp)%name )       
     1523!                                                           TRIM( chem_species(lsp)%name )
    15231524            IF (av == 0)  THEN
    15241525               DO  i = nxl, nxr
    15251526                  DO  j = nys, nyn
    15261527                     DO  k = nzb_do, nzt_do
    1527                          local_pf(i,j,k) = MERGE(                             &
    1528                                              chem_species(lsp)%conc(k,j,i),   &
    1529                                              REAL( fill_value, KIND = wp ),   &
     1528                         local_pf(i,j,k) = MERGE(                                                  &
     1529                                             chem_species(lsp)%conc(k,j,i),                        &
     1530                                             REAL( fill_value, KIND = wp ),                        &
    15301531                                             BTEST( wall_flags_total_0(k,j,i), 0 ) )
    15311532                     ENDDO
     
    15381539                  DO  j = nys, nyn
    15391540                     DO  k = nzb_do, nzt_do
    1540                          local_pf(i,j,k) = MERGE(                             &
    1541                                              chem_species(lsp)%conc_av(k,j,i),&
    1542                                              REAL( fill_value, KIND = wp ),   &
     1541                         local_pf(i,j,k) = MERGE(                                                  &
     1542                                             chem_species(lsp)%conc_av(k,j,i),                     &
     1543                                             REAL( fill_value, KIND = wp ),                        &
    15431544                                             BTEST( wall_flags_total_0(k,j,i), 0 ) )
    15441545                     ENDDO
     
    15561557
    15571558
    1558 !------------------------------------------------------------------------------!
     1559!--------------------------------------------------------------------------------------------------!
    15591560! Description:
    15601561! ------------
    15611562!> Subroutine defining mask output variables for chemical species
    1562 !------------------------------------------------------------------------------!
     1563!--------------------------------------------------------------------------------------------------!
    15631564 SUBROUTINE chem_data_output_mask( av, variable, found, local_pf, mid )
    15641565
     
    15691570    CHARACTER(LEN=*)  ::  variable    !<
    15701571
    1571     INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or time-averaged data
     1572    INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or
     1573                                     !< time-averaged data
     1574    INTEGER(iwp) ::  i               !< grid index along x-direction
     1575    INTEGER(iwp) ::  im              !< loop index for masked variables
     1576    INTEGER(iwp) ::  j               !< grid index along y-direction
     1577    INTEGER(iwp) ::  jm              !< loop index for masked variables
     1578    INTEGER(iwp) ::  k               !< grid index along z-direction
     1579    INTEGER(iwp) ::  kk              !< masked output index along z-direction
     1580    INTEGER(iwp) ::  ktt             !< k index of highest terrain surface
    15721581    INTEGER(iwp) ::  lsp
    1573     INTEGER(iwp) ::  i               !< grid index along x-direction
    1574     INTEGER(iwp) ::  j               !< grid index along y-direction
    1575     INTEGER(iwp) ::  k               !< grid index along z-direction
    1576     INTEGER(iwp) ::  im              !< loop index for masked variables
    1577     INTEGER(iwp) ::  jm              !< loop index for masked variables
    1578     INTEGER(iwp) ::  kk              !< masked output index along z-direction
    15791582    INTEGER(iwp) ::  mid             !< masked output running index
    1580     INTEGER(iwp) ::  ktt             !< k index of highest terrain surface
    1581    
     1583
    15821584    LOGICAL ::  found
    1583    
    1584     REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
    1585               local_pf   !<
     1585
     1586    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  local_pf   !<
    15861587
    15871588    REAL(wp), PARAMETER ::  fill_value = -9999.0_wp    !< value for the _FillValue attribute
    15881589
    15891590!
    1590 !-- local variables.
     1591!-- Local variables.
    15911592
    15921593    spec_name = TRIM( variable(4:) )
     
    15941595
    15951596    DO  lsp=1,nspec
    1596        IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN             
     1597       IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
    15971598!
    15981599!-- todo: remove or replace by "CALL message" mechanism (kanani)
    15991600!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
    1600 !                                                        TRIM( chem_species(lsp)%name )       
     1601!                                                        TRIM( chem_species(lsp)%name )
    16011602          IF (av == 0)  THEN
    16021603             IF ( .NOT. mask_surface(mid) )  THEN
    16031604
    16041605                DO  i = 1, mask_size_l(mid,1)
    1605                    DO  j = 1, mask_size_l(mid,2)
    1606                       DO  k = 1, mask_size(mid,3)
    1607                           local_pf(i,j,k) = chem_species(lsp)%conc(  &
    1608                                                mask_k(mid,k),        &
    1609                                                mask_j(mid,j),        &
    1610                                                mask_i(mid,i)      )
     1606                   DO  j = 1, mask_size_l(mid,2)
     1607                      DO  k = 1, mask_size(mid,3)
     1608                          local_pf(i,j,k) = chem_species(lsp)%conc( mask_k(mid,k),                 &
     1609                                                                    mask_j(mid,j),                 &
     1610                                                                    mask_i(mid,i)      )
    16111611                      ENDDO
    16121612                   ENDDO
     
    16141614
    16151615             ELSE
    1616 !             
     1616!
    16171617!--             Terrain-following masked output
    16181618                DO  i = 1, mask_size_l(mid,1)
    16191619                   DO  j = 1, mask_size_l(mid,2)
     1620!
    16201621!--                   Get k index of the highest terraing surface
    16211622                      im = mask_i(mid,i)
    16221623                      jm = mask_j(mid,j)
    1623                       ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), DIM = 1 ) - 1
     1624                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 ) ),        &
     1625                                    DIM = 1 ) - 1
    16241626                      DO  k = 1, mask_size_l(mid,3)
    16251627                         kk = MIN( ktt+mask_k(mid,k), nzt+1 )
     1628!
    16261629!--                      Set value if not in building
    16271630                         IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
     
    16411644                   DO  j = 1, mask_size_l(mid,2)
    16421645                      DO  k =  1, mask_size_l(mid,3)
    1643                           local_pf(i,j,k) = chem_species(lsp)%conc_av(  &
    1644                                                mask_k(mid,k),           &
    1645                                                mask_j(mid,j),           &
    1646                                                mask_i(mid,i)         )
     1646                          local_pf(i,j,k) = chem_species(lsp)%conc_av(  mask_k(mid,k),             &
     1647                                                                        mask_j(mid,j),             &
     1648                                                                        mask_i(mid,i)         )
    16471649                      ENDDO
    16481650                   ENDDO
     
    16501652
    16511653             ELSE
    1652 !             
     1654!
    16531655!--             Terrain-following masked output
    16541656                DO  i = 1, mask_size_l(mid,1)
    16551657                   DO  j = 1, mask_size_l(mid,2)
     1658!
    16561659!--                   Get k index of the highest terraing surface
    16571660                      im = mask_i(mid,i)
    16581661                      jm = mask_j(mid,j)
    1659                       ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), DIM = 1 ) - 1
     1662                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )),         &
     1663                                    DIM = 1 ) - 1
    16601664                      DO  k = 1, mask_size_l(mid,3)
    16611665                         kk = MIN( ktt+mask_k(mid,k), nzt+1 )
     1666!
    16621667!--                      Set value if not in building
    16631668                         IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
     
    16801685    RETURN
    16811686
    1682  END SUBROUTINE chem_data_output_mask     
    1683 
    1684 
    1685 !------------------------------------------------------------------------------!
     1687 END SUBROUTINE chem_data_output_mask
     1688
     1689
     1690!--------------------------------------------------------------------------------------------------!
    16861691! Description:
    16871692! ------------
    16881693!> Subroutine defining appropriate grid for netcdf variables.
    16891694!> It is called out from subroutine netcdf.
    1690 !------------------------------------------------------------------------------!
     1695!--------------------------------------------------------------------------------------------------!
    16911696 SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
    16921697
    16931698
    16941699    CHARACTER (LEN=*), INTENT(IN)  ::  var          !<
    1695     LOGICAL, INTENT(OUT)           ::  found        !<
     1700
    16961701    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x       !<
    16971702    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y       !<
    16981703    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z       !<
    16991704
     1705    LOGICAL, INTENT(OUT)           ::  found        !<
     1706
    17001707    found  = .TRUE.
    17011708
    1702     IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' )  THEN                   !< always the same grid for chemistry variables
     1709    IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' )  THEN                !< always the same grid for
     1710                                                                         !< chemistry variables
    17031711       grid_x = 'x'
    17041712       grid_y = 'y'
    1705        grid_z = 'zu'                             
     1713       grid_z = 'zu'
    17061714    ELSE
    17071715       found  = .FALSE.
     
    17151723
    17161724
    1717 !------------------------------------------------------------------------------!
     1725!--------------------------------------------------------------------------------------------------!
    17181726! Description:
    17191727! ------------
    17201728!> Subroutine defining header output for chemistry model
    1721 !------------------------------------------------------------------------------!
     1729!--------------------------------------------------------------------------------------------------!
    17221730 SUBROUTINE chem_header( io )
    17231731
    1724 
    1725     INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
    1726     INTEGER(iwp)  :: lsp                       !< running index for chem spcs
    1727     INTEGER(iwp)  :: cs_fixed
    17281732    CHARACTER (LEN=80)  :: docsflux_chr
    17291733    CHARACTER (LEN=80)  :: docsinit_chr
    1730 !
    1731 ! Get name of chemical mechanism from chem_gasphase_mod
     1734
     1735    INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
     1736
     1737    INTEGER(iwp)  :: cs_fixed
     1738    INTEGER(iwp)  :: lsp                       !< running index for chem spcs
     1739
     1740!
     1741!-- Get name of chemical mechanism from chem_gasphase_mod
    17321742    CALL get_mechanism_name
    17331743!
     
    17361746!
    17371747!-- Gasphase reaction status
    1738     IF ( chem_gasphase_on )  THEN 
     1748    IF ( chem_gasphase_on )  THEN
    17391749       WRITE( io, 2 )
    17401750    ELSE
     
    17461756!
    17471757!-- Emission mode info
    1748 !-- At the moment the evaluation is done with both emiss_lod and mode_emis
    1749 !-- but once salsa has been migrated to emiss_lod the .OR. mode_emis
    1750 !-- conditions can be removed (ecc 20190513)
    1751     IF     ( (emiss_lod == 1) .OR. (mode_emis == 'DEFAULT') )        THEN
     1758!-- At the moment the evaluation is done with both emiss_lod and mode_emis but once salsa has been
     1759!-- migrated to emiss_lod the .OR. mode_emis conditions can be removed (ecc 20190513)
     1760    IF     ( ( emiss_lod == 1 )  .OR.  ( mode_emis == 'DEFAULT' ) )        THEN
    17521761        WRITE ( io, 5 )
    1753     ELSEIF ( (emiss_lod == 0) .OR. (mode_emis == 'PARAMETERIZED') )  THEN
     1762    ELSEIF ( ( emiss_lod == 0 )  .OR.  ( mode_emis == 'PARAMETERIZED' ) )  THEN
    17541763        WRITE ( io, 6 )
    1755     ELSEIF ( (emiss_lod == 2) .OR. (mode_emis == 'PRE-PROCESSED') )  THEN
     1764    ELSEIF ( ( emiss_lod == 2 )  .OR.  ( mode_emis == 'PRE-PROCESSED' ) )  THEN
    17561765        WRITE ( io, 7 )
    17571766    ENDIF
     
    17591768!-- Photolysis scheme info
    17601769    IF ( photolysis_scheme == "simple" )  THEN
    1761        WRITE( io, 8 ) 
     1770       WRITE( io, 8 )
    17621771    ELSEIF (photolysis_scheme == "constant" )  THEN
    17631772       WRITE( io, 9 )
     
    17661775!-- Emission flux info
    17671776    lsp = 1
    1768     docsflux_chr ='Chemical species for surface emission flux: ' 
     1777    docsflux_chr ='Chemical species for surface emission flux: '
    17691778    DO WHILE ( surface_csflux_name(lsp) /= 'novalue' )
    1770        docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ',' 
     1779       docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ','
    17711780       IF ( LEN_TRIM( docsflux_chr ) >= 75 )  THEN
    17721781          WRITE ( io, 10 )  docsflux_chr
     
    17801789    ENDIF
    17811790!
    1782 !-- initializatoin of Surface and profile chemical species
     1791!-- Initialization of Surface and profile chemical species
    17831792    lsp = 1
    1784     docsinit_chr ='Chemical species for initial surface and profile emissions: ' 
     1793    docsinit_chr ='Chemical species for initial surface and profile emissions: '
    17851794    DO WHILE ( cs_name(lsp) /= 'novalue' )
    1786        docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ',' 
     1795       docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ','
    17871796       IF ( LEN_TRIM( docsinit_chr ) >= 75 )  THEN
    1788         WRITE ( io, 11 )  docsinit_chr
    1789         docsinit_chr = '       '
     1797          WRITE ( io, 11 )  docsinit_chr
     1798          docsinit_chr = '       '
    17901799       ENDIF
    17911800       lsp = lsp + 1
     
    18031812
    18041813!
    1805 !-- number of variable and fix chemical species and number of reactions
     1814!-- Number of variable and fix chemical species and number of reactions
    18061815    cs_fixed = nspec - nvar
    1807     WRITE ( io, * ) '   --> Chemical Mechanism        : ', cs_mech 
     1816    WRITE ( io, * ) '   --> Chemical Mechanism        : ', cs_mech
    18081817    WRITE ( io, * ) '   --> Chemical species, variable: ', nvar
    18091818    WRITE ( io, * ) '   --> Chemical species, fixed   : ', cs_fixed
     
    18111820
    18121821
    1813 1   FORMAT (//' Chemistry model information:'/                                  &
    1814            ' ----------------------------'/)
     18221   FORMAT (//' Chemistry model information:'/' ----------------------------'/)
    181518232   FORMAT ('    --> Chemical reactions are turned on')
    181618243   FORMAT ('    --> Chemical reactions are turned off')
     
    182118298   FORMAT ('    --> Photolysis scheme used =  simple ')
    182218309   FORMAT ('    --> Photolysis scheme used =  constant ')
    1823 10  FORMAT (/'    ',A) 
    1824 11  FORMAT (/'    ',A) 
     183110  FORMAT (/'    ',A)
     183211  FORMAT (/'    ',A)
    1825183312  FORMAT (/'   Nesting for chemistry variables: ', L1 )
    1826183413  FORMAT (/'   Offline nesting for chemistry variables: ', L1 )
     
    18331841
    18341842
    1835 !------------------------------------------------------------------------------!
     1843!--------------------------------------------------------------------------------------------------!
    18361844! Description:
    18371845! ------------
    18381846!> Subroutine initializating chemistry_model_mod specific arrays
    1839 !------------------------------------------------------------------------------!
     1847!--------------------------------------------------------------------------------------------------!
    18401848 SUBROUTINE chem_init_arrays
    18411849!
     
    18451853
    18461854
    1847 !------------------------------------------------------------------------------!
     1855!--------------------------------------------------------------------------------------------------!
    18481856! Description:
    18491857! ------------
    18501858!> Subroutine initializating chemistry_model_mod
    1851 !------------------------------------------------------------------------------!
     1859!--------------------------------------------------------------------------------------------------!
    18521860 SUBROUTINE chem_init
    18531861
     
    18561864!-- introduced additional interfaces for on-demand emission update
    18571865
    1858 !    USE chem_emissions_mod,                                                    &
     1866!    USE chem_emissions_mod,                                                                        &
    18591867!        ONLY:  chem_emissions_init
    1860        
    1861     USE chem_emissions_mod,                                                    &
    1862         ONLY:  chem_emissions_init, chem_emissions_header_init
    1863        
    1864     USE netcdf_data_input_mod,                                                 &
     1868
     1869    USE chem_emissions_mod,                                                                        &
     1870        ONLY:  chem_emissions_header_init, chem_emissions_init
     1871
     1872    USE netcdf_data_input_mod,                                                                     &
    18651873        ONLY:  init_3d
    18661874
     
    18741882!
    18751883!-- Next statement is to avoid compiler warning about unused variables
    1876     IF ( ( ilu_arable + ilu_coniferous_forest + ilu_deciduous_forest + ilu_mediterrean_scrub + &
    1877            ilu_permanent_crops + ilu_savanna + ilu_semi_natural_veg + ilu_tropical_forest +    &
     1884    IF ( ( ilu_arable + ilu_coniferous_forest + ilu_deciduous_forest + ilu_mediterrean_scrub +     &
     1885           ilu_permanent_crops + ilu_savanna + ilu_semi_natural_veg + ilu_tropical_forest +        &
    18781886           ilu_urban ) == 0 )  CONTINUE
    18791887
    18801888!
    18811889!-- 20200203 (ECC)
    1882 !-- calls specific emisisons initialization subroutines
    1883 !-- for legacy mode and on-demand mode         
     1890!-- Calls specific emisisons initialization subroutines for legacy mode and on-demand mode
    18841891
    18851892!    IF ( emissions_anthropogenic )  CALL chem_emissions_init
     
    18971904
    18981905!
    1899 !-- Chemistry variables will be initialized if availabe from dynamic
    1900 !-- input file. Note, it is possible to initialize only part of the chemistry
    1901 !-- variables from dynamic input.
     1906!-- Chemistry variables will be initialized if availabe from dynamic input file. Note, it is
     1907!-- possible to initialize only part of the chemistry variables from dynamic input.
    19021908    IF ( INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
    19031909       DO  n = 1, nspec
     
    19171923
    19181924
    1919 !------------------------------------------------------------------------------!
     1925!--------------------------------------------------------------------------------------------------!
    19201926! Description:
    19211927! ------------
    19221928!> Subroutine initializating chemistry_model_mod
    19231929!> internal workaround for chem_species dependency in chem_check_parameters
    1924 !------------------------------------------------------------------------------!
     1930!--------------------------------------------------------------------------------------------------!
    19251931 SUBROUTINE chem_init_internal
    1926  
     1932
    19271933    USE pegrid
    19281934
    1929     USE netcdf_data_input_mod,                                                 &
    1930         ONLY:  chem_emis, chem_emis_att, input_pids_dynamic, init_3d,          &
     1935    USE netcdf_data_input_mod,                                                                     &
     1936        ONLY:  chem_emis, chem_emis_att, input_pids_dynamic, init_3d,                              &
    19311937               netcdf_data_input_chemistry_data
    19321938
     
    19351941    INTEGER(iwp) ::  i                 !< running index for for horiz numerical grid points
    19361942    INTEGER(iwp) ::  j                 !< running index for for horiz numerical grid points
     1943    INTEGER(iwp) ::  lpr_lev           !< running index for chem spcs profile level
    19371944    INTEGER(iwp) ::  lsp               !< running index for chem spcs
    1938     INTEGER(iwp) ::  lpr_lev           !< running index for chem spcs profile level
    19391945
    19401946    REAL(wp)     ::  flag              !< flag for masking topography/building grid points
     
    19591965    ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
    19601966    ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
    1961     ALLOCATE( phot_frequen(nphot) ) 
     1967    ALLOCATE( phot_frequen(nphot) )
    19621968    ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) )
    19631969    ALLOCATE( bc_cs_t_val(nspec) )
     
    19761982       chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_3 (:,:,:,lsp)
    19771983
    1978        ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))                   
     1984       ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))
    19791985       chem_species(lsp)%cssws_av    = 0.0_wp
    19801986!
    19811987!--    The following block can be useful when emission module is not applied. &
    1982 !--    if emission module is applied the following block will be overwritten.
    1983        ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))   
    1984        ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))   
    1985        ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
    1986        ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))   
    1987        chem_species(lsp)%flux_s_cs = 0.0_wp                                     
    1988        chem_species(lsp)%flux_l_cs = 0.0_wp                                     
    1989        chem_species(lsp)%diss_s_cs = 0.0_wp                                     
    1990        chem_species(lsp)%diss_l_cs = 0.0_wp                                     
    1991 !
    1992 !--   Allocate memory for initial concentration profiles
    1993 !--   (concentration values come from namelist)
    1994 !--   (@todo (FK): Because of this, chem_init is called in palm before
    1995 !--               check_parameters, since conc_pr_init is used there.
    1996 !--               We have to find another solution since chem_init should
    1997 !--               eventually be called from init_3d_model!!)
     1988!--    If emission module is applied the following block will be overwritten.
     1989       ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))
     1990       ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))
     1991       ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
     1992       ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
     1993       chem_species(lsp)%flux_s_cs = 0.0_wp
     1994       chem_species(lsp)%flux_l_cs = 0.0_wp
     1995       chem_species(lsp)%diss_s_cs = 0.0_wp
     1996       chem_species(lsp)%diss_l_cs = 0.0_wp
     1997!
     1998!--   Allocate memory for initial concentration profiles (concentration values come from namelist)
     1999!--   (@todo (FK): Because of this, chem_init is called in palm before check_parameters, since
     2000!--                conc_pr_init is used there.
     2001!--                We have to find another solution since chem_init should eventually be called from
     2002!--                init_3d_model!!)
    19982003       ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) )
    19992004       chem_species(lsp)%conc_pr_init(:) = 0.0_wp
     
    20022007
    20032008!
    2004 !-- For some passive scalars decycling may be enabled. This case, the lateral
    2005 !-- boundary conditions are non-cyclic for these scalars (chemical species
    2006 !-- and aerosols), while the other scalars may have
    2007 !-- cyclic boundary conditions. However, large gradients near the boundaries
    2008 !-- may produce stationary numerical oscillations near the lateral boundaries
    2009 !-- when a higher-order scheme is applied near these boundaries.
    2010 !-- To get rid-off this, set-up additional flags that control the order of the
    2011 !-- scalar advection scheme near the lateral boundaries for passive scalars
    2012 !-- with decycling.
     2009!-- For some passive scalars decycling may be enabled. This case, the lateral boundary conditions
     2010!-- are non-cyclic for these scalars (chemical species and aerosols), while the other scalars may
     2011!-- have cyclic boundary conditions. However, large gradients near the boundaries may produce
     2012!-- stationary numerical oscillations near the lateral boundaries when a higher-order scheme is
     2013!-- applied near these boundaries.
     2014!-- To get rid-off this, set-up additional flags that control the order of the scalar advection
     2015!-- scheme near the lateral boundaries for passive scalars with decycling.
    20132016    IF ( scalar_advec == 'ws-scheme' )  THEN
    20142017       ALLOCATE( cs_advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    20152018!
    2016 !--    In case of decyling, set Neumann boundary conditions for wall_flags_total_0
    2017 !--    bit 31 instead of cyclic boundary conditions.
    2018 !--    Bit 31 is used to identify extended degradation zones (please see
    2019 !--    following comment).
    2020 !--    Note, since several also other modules like Salsa or other future
    2021 !--    one may access this bit but may have other boundary conditions, the
    2022 !--    original value of wall_flags_total_0 bit 31 must not be modified. Hence,
    2023 !--    store the boundary conditions directly on cs_advc_flags_s.
    2024 !--    cs_advc_flags_s will be later overwritten in ws_init_flags_scalar and
    2025 !--    bit 31 won't be used to control the numerical order.
    2026 !--    Initialize with flag 31 only.
     2019!--    In case of decyling, set Neumann boundary conditions for wall_flags_total_0 bit 31 instead of
     2020!--    cyclic boundary conditions.
     2021!--    Bit 31 is used to identify extended degradation zones (please see following comment).
     2022!--    Note, since several also other modules like Salsa or other future one may access this bit but
     2023!--    may have other boundary conditions, the original value of wall_flags_total_0 bit 31 must not
     2024!--    be modified. Hence, store the boundary conditions directly on cs_advc_flags_s.
     2025!--    cs_advc_flags_s will be later overwritten in ws_init_flags_scalar and bit 31 won't be used to
     2026!--    control the numerical order.
     2027!--    Initialize with flag 31 only.
    20272028       cs_advc_flags_s = 0
    20282029       cs_advc_flags_s = MERGE( IBSET( cs_advc_flags_s, 31 ), 0, BTEST( wall_flags_total_0, 31 ) )
     
    20312032          IF ( nys == 0  )  THEN
    20322033             DO  i = 1, nbgp
    2033                 cs_advc_flags_s(:,nys-i,:) = MERGE(                            &
    2034                                         IBSET( cs_advc_flags_s(:,nys,:), 31 ), &
    2035                                         IBCLR( cs_advc_flags_s(:,nys,:), 31 ), &
    2036                                         BTEST( cs_advc_flags_s(:,nys,:), 31 )  &
     2034                cs_advc_flags_s(:,nys-i,:) = MERGE(                                                &
     2035                                                    IBSET( cs_advc_flags_s(:,nys,:), 31 ),        &
     2036                                                    IBCLR( cs_advc_flags_s(:,nys,:), 31 ),        &
     2037                                                    BTEST( cs_advc_flags_s(:,nys,:), 31 )          &
    20372038                                                  )
    20382039             ENDDO
     
    20402041          IF ( nyn == ny )  THEN
    20412042             DO  i = 1, nbgp
    2042                 cs_advc_flags_s(:,nyn+i,:) = MERGE(                            &
    2043                                         IBSET( cs_advc_flags_s(:,nyn,:), 31 ), &
    2044                                         IBCLR( cs_advc_flags_s(:,nyn,:), 31 ), &
    2045                                         BTEST( cs_advc_flags_s(:,nyn,:), 31 )  &
     2043                cs_advc_flags_s(:,nyn+i,:) = MERGE(                                                &
     2044                                                    IBSET( cs_advc_flags_s(:,nyn,:), 31 ),        &
     2045                                                    IBCLR( cs_advc_flags_s(:,nyn,:), 31 ),        &
     2046                                                    BTEST( cs_advc_flags_s(:,nyn,:), 31 )          &
    20462047                                                  )
    20472048             ENDDO
     
    20512052          IF ( nxl == 0  )  THEN
    20522053             DO  i = 1, nbgp
    2053                 cs_advc_flags_s(:,:,nxl-i) = MERGE(                            &
    2054                                         IBSET( cs_advc_flags_s(:,:,nxl), 31 ), &
    2055                                         IBCLR( cs_advc_flags_s(:,:,nxl), 31 ), &
    2056                                         BTEST( cs_advc_flags_s(:,:,nxl), 31 )  &
     2054                cs_advc_flags_s(:,:,nxl-i) = MERGE(                                                &
     2055                                                    IBSET( cs_advc_flags_s(:,:,nxl), 31 ),        &
     2056                                                    IBCLR( cs_advc_flags_s(:,:,nxl), 31 ),        &
     2057                                                    BTEST( cs_advc_flags_s(:,:,nxl), 31 )          &
    20572058                                                  )
    20582059             ENDDO
     
    20602061          IF ( nxr == nx )  THEN
    20612062             DO  i = 1, nbgp
    2062                 cs_advc_flags_s(:,:,nxr+i) = MERGE(                            &
    2063                                         IBSET( cs_advc_flags_s(:,:,nxr), 31 ), &
    2064                                         IBCLR( cs_advc_flags_s(:,:,nxr), 31 ), &
    2065                                         BTEST( cs_advc_flags_s(:,:,nxr), 31 )  &
     2063                cs_advc_flags_s(:,:,nxr+i) = MERGE(                                                &
     2064                                                    IBSET( cs_advc_flags_s(:,:,nxr), 31 ), &
     2065                                                    IBCLR( cs_advc_flags_s(:,:,nxr), 31 ), &
     2066                                                    BTEST( cs_advc_flags_s(:,:,nxr), 31 )  &
    20662067                                                  )
    20672068             ENDDO
     
    20692070       ENDIF
    20702071!
    2071 !--    To initialize advection flags appropriately, pass the boundary flags. 
    2072 !--    The last argument indicates that a passive scalar is treated, where
    2073 !--    the horizontal advection terms are degraded already 2 grid points before
    2074 !--    the lateral boundary to avoid stationary oscillations at large-gradients.
    2075 !--    Also, extended degradation zones are applied, where horizontal advection of
    2076 !--    passive scalars is discretized by first-order scheme at all grid points
    2077 !--    that in the vicinity of buildings (<= 3 grid points). Even though no
    2078 !--    building is within the numerical stencil, first-order scheme is used.
     2072!--    To initialize advection flags appropriately, pass the boundary flags.
     2073!--    The last argument indicates that a passive scalar is treated, where the horizontal advection
     2074!--    terms are degraded already 2 grid points before the lateral boundary to avoid stationary
     2075!--    oscillations at large-gradients.
     2076!--    Also, extended degradation zones are applied, where horizontal advection of passive scalars
     2077!--    is discretized by first-order scheme at all grid points that in the vicinity of buildings
     2078!--    (<= 3 grid points). Even though no building is within the numerical stencil, first-order
     2079!--    scheme is used.
    20792080!--    At fourth and fifth grid point the order of the horizontal advection scheme
    20802081!--    is successively upgraded.
    2081 !--    These extended degradation zones are used to avoid stationary numerical
    2082 !--    oscillations, which are responsible for high concentration maxima that may
    2083 !--    appear under shear-free stable conditions.
     2082!--    These extended degradation zones are used to avoid stationary numerical oscillations, which
     2083!--    are responsible for high concentration maxima that may appear under shear-free stable
     2084!--    conditions.
    20842085       CALL ws_init_flags_scalar( bc_dirichlet_cs_l  .OR.  bc_radiation_cs_l,                      &
    20852086                                  bc_dirichlet_cs_n  .OR.  bc_radiation_cs_n,                      &
     
    20892090    ENDIF
    20902091!
    2091 !-- Initial concentration of profiles is prescribed by parameters cs_profile
    2092 !-- and cs_heights in the namelist &chemistry_parameters
    2093 
     2092!-- Initial concentration of profiles is prescribed by parameters cs_profile and cs_heights in the
     2093!-- namelist &chemistry_parameters.
    20942094    CALL chem_init_profiles
    2095 !   
    2096 !-- In case there is dynamic input file, create a list of names for chemistry
    2097 !-- initial input files. Also, initialize array that indicates whether the
    2098 !-- respective variable is on file or not.
    2099     IF ( input_pids_dynamic )  THEN   
     2095!
     2096!-- In case there is dynamic input file, create a list of names for chemistry initial input files.
     2097!-- Also, initialize array that indicates whether the respective variable is on file or not.
     2098    IF ( input_pids_dynamic )  THEN
    21002099       ALLOCATE( init_3d%var_names_chem(1:nspec) )
    21012100       ALLOCATE( init_3d%from_file_chem(1:nspec) )
    21022101       init_3d%from_file_chem(:) = .FALSE.
    2103        
     2102
    21042103       DO  lsp = 1, nspec
    21052104          init_3d%var_names_chem(lsp) = init_3d%init_char // TRIM( chem_species(lsp)%name )
     
    21082107!
    21092108!-- Initialize model variables
    2110     IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
     2109    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                                &
    21112110         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
    21122111!
     
    21262125                   ENDDO
    21272126                ENDDO
    2128              ENDDO   
     2127             ENDDO
    21292128          ENDDO
    21302129
    2131        ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
    2132        THEN
    2133 
    2134           DO  lsp = 1, nspec
     2130       ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles') /= 0 )  THEN
     2131
     2132          DO  lsp = 1, nspec
    21352133             DO  i = nxlg, nxrg
    21362134                DO  j = nysg, nyng
     
    21472145!
    21482146!--    If required, change the surface chem spcs at the start of the 3D run
    2149        IF ( cs_surface_initial_change(1) /= 0.0_wp )  THEN           
    2150           DO  lsp = 1, nspec 
     2147       IF ( cs_surface_initial_change(1) /= 0.0_wp )  THEN
     2148          DO  lsp = 1, nspec
    21512149             DO  i = nxlg, nxrg
    21522150                DO  j = nysg, nyng
    21532151                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb,j,i), 0 ) )
    2154                    chem_species(lsp)%conc(nzb,j,i) = chem_species(lsp)%conc(nzb,j,i) +  &
     2152                   chem_species(lsp)%conc(nzb,j,i) = chem_species(lsp)%conc(nzb,j,i) +             &
    21552153                                                     cs_surface_initial_change(lsp) * flag
    21562154                ENDDO
    21572155             ENDDO
    21582156          ENDDO
    2159        ENDIF 
     2157       ENDIF
    21602158
    21612159    ENDIF
    21622160!
    2163 !-- Initiale old and new time levels. Note, this has to be done also in restart runs
     2161!-- Initial old and new time levels. Note, this has to be done also in restart runs
    21642162    DO  lsp = 1, nvar
    2165        chem_species(lsp)%tconc_m = 0.0_wp                     
    2166        chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
     2163       chem_species(lsp)%tconc_m = 0.0_wp
     2164       chem_species(lsp)%conc_p  = chem_species(lsp)%conc
    21672165    ENDDO
    21682166
     
    21732171!--       IF( myid == 0 )  THEN
    21742172!--          WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,TRIM( phot_names(lsp) )
    2175 !--       ENDIF 
     2173!--       ENDIF
    21762174       phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  =>  freq_1(:,:,:,lsp)
    21772175    ENDDO
     
    21842182
    21852183
    2186 !------------------------------------------------------------------------------!
     2184!--------------------------------------------------------------------------------------------------!
    21872185! Description:
    21882186! ------------
    2189 !> Subroutine defining initial vertical profiles of chemical species (given by
    2190 !> namelist parameters chem_profiles and chem_heights)  --> which should work
    2191 !> analogue to parameters u_profile, v_profile and uv_heights)
    2192 !------------------------------------------------------------------------------!
    2193  SUBROUTINE chem_init_profiles             
     2187!> Subroutine defining initial vertical profiles of chemical species (given by namelist parameters
     2188!> chem_profiles and chem_heights)  --> which should work analogically to parameters u_profile,
     2189!> v_profile and uv_heights)
     2190!--------------------------------------------------------------------------------------------------!
     2191 SUBROUTINE chem_init_profiles
    21942192
    21952193    USE chem_modules
     
    21972195!
    21982196!-- Local variables
     2197    INTEGER ::  lpr_lev    !< running index for profile level for each chem spcs.
    21992198    INTEGER ::  lsp        !< running index for number of species in derived data type species_def
    2200     INTEGER ::  lsp_usr     !< running index for number of species (user defined)  in cs_names, cs_profiles etc
    2201     INTEGER ::  lpr_lev    !< running index for profile level for each chem spcs.
     2199    INTEGER ::  lsp_usr    !< running index for number of species (user defined)  in cs_names,
     2200                           !< cs_profiles etc
    22022201    INTEGER ::  npr_lev    !< the next available profile lev
    22032202!
     
    22072206!-- "cs_heights" are prescribed, their values will!override the constant profile given by
    22082207!-- "cs_surface".
    2209 !     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN 
     2208!     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
    22102209    lsp_usr = 1
    22112210    DO  WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )   !'novalue' is the default
    22122211       DO  lsp = 1, nspec                                !
    22132212!
    2214 !--       create initial profile (conc_pr_init) for each chemical species
    2215           IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN   !
     2213!--       Create initial profile (conc_pr_init) for each chemical species
     2214          IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN
    22162215             IF ( cs_profile(lsp_usr,1) == 9999999.9_wp )  THEN
    22172216!
    2218 !--            set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
     2217!--            Set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of
     2218!--            each species
    22192219                DO lpr_lev = 0, nzt+1
    22202220                   chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
     
    22362236                         IF ( npr_lev == 100 )  THEN
    22372237                            message_string = 'number of chem spcs exceeding the limit'
    2238                             CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )               
     2238                            CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )
    22392239                            EXIT
    22402240                         ENDIF
     
    22422242                   ENDIF
    22432243                   IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr,npr_lev+1) /= 9999999.9_wp )  THEN
    2244                       chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +       &
    2245                            ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                          &
    2246                            ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *  &
     2244                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +     &
     2245                           ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                        &
     2246                           ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *&
    22472247                           ( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
    22482248                   ELSE
     
    22522252             ENDIF
    22532253!
    2254 !--       If a profile is prescribed explicity using cs_profiles and cs_heights, then 
    2255 !--       chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
    2256 !--       on the cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:).
     2254!--       If a profile is prescribed explicity using cs_profiles and cs_heights, then
     2255!--       chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based on the
     2256!--       cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:).
    22572257          ENDIF
    22582258
     
    22652265 END SUBROUTINE chem_init_profiles
    22662266
    2267  
    2268 !------------------------------------------------------------------------------!
     2267
     2268!--------------------------------------------------------------------------------------------------!
    22692269! Description:
    22702270! ------------
    22712271!> Subroutine to integrate chemical species in the given chemical mechanism
    2272 !------------------------------------------------------------------------------!
     2272!--------------------------------------------------------------------------------------------------!
    22732273 SUBROUTINE chem_integrate_ij( i, j )
    22742274
    2275     USE statistics,                                                          &
     2275    USE statistics,                                                                                &
    22762276        ONLY:  weight_pres
    22772277
    2278     USE control_parameters,                                                  &
     2278    USE control_parameters,                                                                        &
    22792279        ONLY:  dt_3d, intermediate_timestep_count, time_since_reference_point
    22802280
     
    22832283    INTEGER,INTENT(IN)       :: j
    22842284!
    2285 !--   local variables
    2286     INTEGER(iwp) ::  lsp                                                     !< running index for chem spcs.
    2287     INTEGER(iwp) ::  lph                                                     !< running index for photolysis frequencies
    2288     INTEGER, DIMENSION(20)    :: istatus
     2285!-- Local variables
     2286    INTEGER(iwp) ::  lph                                 !< running index for photolysis frequencies
     2287    INTEGER(iwp) ::  lsp                                 !< running index for chem spcs.
     2288
     2289    INTEGER, DIMENSION(20)        :: istatus
     2290
     2291    INTEGER,DIMENSION(nzb+1:nzt)  :: nacc          !< Number of accepted steps
     2292    INTEGER,DIMENSION(nzb+1:nzt)  :: nrej          !< Number of rejected steps
     2293
     2294    REAL(wp)                         ::  conv                                !< conversion factor
     2295    REAL(kind=wp)                    ::  dt_chem
     2296
     2297    REAL(wp), PARAMETER              ::  fr2ppm  = 1.0e6_wp              !< Conversion factor
     2298                                                                         !< fraction to ppm
     2299!    REAL(wp), PARAMETER              ::  xm_air  = 28.96_wp              !< Mole mass of dry air
     2300!    REAL(wp), PARAMETER              ::  xm_h2o  = 18.01528_wp           !< Mole mass of water vapor
     2301    REAL(wp), PARAMETER              ::  p_std   = 101325.0_wp           !< standard pressure (Pa)
     2302    REAL(wp), PARAMETER              ::  ppm2fr  = 1.0e-6_wp             !< Conversion factor ppm to
     2303                                                                         !< fraction
     2304    REAL(wp), PARAMETER              ::  t_std   = 273.15_wp             !< standard pressure (Pa)
     2305    REAL(wp), PARAMETER              ::  vmolcm  = 22.414e3_wp           !< Mole volume (22.414 l)
     2306                                                                         !< in cm^3
     2307    REAL(wp), PARAMETER              ::  xna     = 6.022e23_wp           !< Avogadro number
     2308                                                                         !< (molecules/mol)
     2309
     2310    REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
     2311
     2312    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact
     2313    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact_i    !< conversion factor between
     2314                                                                              !< molecules cm^{-3} and ppm
     2315    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_qvap
     2316    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_temp
     2317
    22892318    REAL(kind=wp), DIMENSION(nzb+1:nzt,nspec)                :: tmp_conc
    2290     REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_temp
    2291     REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_qvap
    22922319    REAL(kind=wp), DIMENSION(nzb+1:nzt,nphot)                :: tmp_phot
    2293     REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact
    2294     REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact_i    !< conversion factor between
    2295                                                                               !<    molecules cm^{-3} and ppm
    2296 
    2297     INTEGER,DIMENSION(nzb+1:nzt)                            :: nacc          !< Number of accepted steps
    2298     INTEGER,DIMENSION(nzb+1:nzt)                            :: nrej          !< Number of rejected steps
    2299 
    2300     REAL(wp)                         ::  conv                                !< conversion factor
    2301     REAL(wp), PARAMETER              ::  ppm2fr  = 1.0e-6_wp                 !< Conversion factor ppm to fraction
    2302     REAL(wp), PARAMETER              ::  fr2ppm  = 1.0e6_wp                  !< Conversion factor fraction to ppm
    2303 !    REAL(wp), PARAMETER              ::  xm_air  = 28.96_wp                  !< Mole mass of dry air
    2304 !    REAL(wp), PARAMETER              ::  xm_h2o  = 18.01528_wp               !< Mole mass of water vapor
    2305     REAL(wp), PARAMETER              ::  t_std   = 273.15_wp                 !< standard pressure (Pa)
    2306     REAL(wp), PARAMETER              ::  p_std   = 101325.0_wp               !< standard pressure (Pa)
    2307     REAL(wp), PARAMETER              ::  vmolcm  = 22.414e3_wp               !< Mole volume (22.414 l) in cm^3
    2308     REAL(wp), PARAMETER              ::  xna     = 6.022e23_wp               !< Avogadro number (molecules/mol)
    2309 
    2310     REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
    2311 
    2312     REAL(kind=wp)  :: dt_chem                                             
     2320
    23132321!
    23142322!-- Set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
     
    23192327       tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt)
    23202328!
    2321 !--    convert ppm to molecules/cm**3
    2322 !--    tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp * 
    2323 !--               hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp ) 
     2329!--    Convert ppm to molecules/cm**3
     2330!--    tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp *
     2331!--               hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp )
    23242332       conv = ppm2fr * xna / vmolcm
    23252333       tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
     
    23282336       IF ( humidity )  THEN
    23292337          IF ( bulk_cloud_model )  THEN
    2330              tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) *  &
     2338             tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) *                              &
    23312339                             xm_air/xm_h2o * fr2ppm * tmp_fact(:)
    23322340          ELSE
     
    23342342          ENDIF
    23352343       ELSE
    2336           tmp_qvap(:) = 0.01 * xm_air/xm_h2o * fr2ppm * tmp_fact(:)          !< Constant value for q if water vapor is not computed
     2344          tmp_qvap(:) = 0.01 * xm_air/xm_h2o * fr2ppm * tmp_fact(:)   !< Constant value for q if
     2345                                                                      !< water vapor is not computed
    23372346       ENDIF
    23382347
    23392348       DO  lsp = 1,nspec
    2340           tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:) 
     2349          tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:)
    23412350       ENDDO
    23422351
    23432352       DO lph = 1,nphot
    2344           tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)               
     2353          tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)
    23452354       ENDDO
    23462355!
     
    23542363       cs_time_step = dt_chem
    23552364
    2356        IF(maxval(rcntrl) > 0.0)  THEN    ! Only if rcntrl is set
     2365       IF ( MAXVAL( rcntrl ) > 0.0 )  THEN    ! Only if rcntrl is set
    23572366          IF( time_since_reference_point <= 2*dt_3d)  THEN
    23582367             rcntrl_local = 0
     
    23642373       END IF
    23652374
    2366        CALL chem_gasphase_integrate ( dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot, &
     2375       CALL chem_gasphase_integrate ( dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot,   &
    23672376            icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus )
    23682377
     
    23782387
    23792388
    2380 !------------------------------------------------------------------------------!
     2389!--------------------------------------------------------------------------------------------------!
    23812390! Description:
    23822391! ------------
    23832392!> Subroutine defining parin for &chemistry_parameters for chemistry model
    2384 !------------------------------------------------------------------------------!
     2393!--------------------------------------------------------------------------------------------------!
    23852394 SUBROUTINE chem_parin
    23862395
     
    23922401
    23932402
    2394     CHARACTER (LEN=80) ::  line                        !< dummy string that contains the current line of the parameter file
    2395 
    2396     REAL(wp), DIMENSION(nmaxfixsteps) ::   my_steps    !< List of fixed timesteps   my_step(1) = 0.0 automatic stepping
    2397     INTEGER(iwp) ::  i                                 !<
    2398     INTEGER(iwp) ::  max_pr_cs_tmp                     !<
    2399 
    2400 
    2401     NAMELIST /chemistry_parameters/  bc_cs_b,                          &
    2402          bc_cs_l,                          &
    2403          bc_cs_n,                          &
    2404          bc_cs_r,                          &
    2405          bc_cs_s,                          &
    2406          bc_cs_t,                          &
    2407          call_chem_at_all_substeps,        &
    2408          chem_debug0,                      &
    2409          chem_debug1,                      &
    2410          chem_debug2,                      &
    2411          chem_gasphase_on,                 &
    2412          chem_mechanism,                   &         
    2413          cs_heights,                       &
    2414          cs_name,                          &
    2415          cs_profile,                       &
    2416          cs_surface,                       &
    2417          cs_surface_initial_change,        &
    2418          cs_vertical_gradient_level,       &
    2419          daytype_mdh,                      &
    2420          deposition_dry,                   &
    2421          emissions_anthropogenic,          &
    2422          emiss_lod,                        &
    2423          emiss_factor_main,                &
    2424          emiss_factor_side,                &
    2425          emiss_read_legacy_mode,           &
    2426          icntrl,                           &
    2427          main_street_id,                   &
    2428          max_street_id,                    &
    2429          mode_emis,                        &
    2430          my_steps,                         &
    2431          nesting_chem,                     &
    2432          nesting_offline_chem,             &
    2433          rcntrl,                           &
    2434          side_street_id,                   &
    2435          photolysis_scheme,                &
    2436          wall_csflux,                      &
    2437          cs_vertical_gradient,             &
    2438          top_csflux,                       &
    2439          surface_csflux,                   &
    2440          surface_csflux_name,              &
     2403    CHARACTER (LEN=80) ::  line                        !< dummy string that contains the current
     2404                                                       !< line of the parameter file
     2405
     2406    INTEGER(iwp) ::  i                                 !<
     2407    INTEGER(iwp) ::  max_pr_cs_tmp                     !<
     2408
     2409    REAL(wp), DIMENSION(nmaxfixsteps) ::   my_steps    !< List of fixed timesteps   my_step(1) = 0.0
     2410                                                       !< automatic stepping
     2411
     2412
     2413    NAMELIST /chemistry_parameters/                                                                &
     2414         bc_cs_b,                                                                                  &
     2415         bc_cs_l,                                                                                  &
     2416         bc_cs_n,                                                                                  &
     2417         bc_cs_r,                                                                                  &
     2418         bc_cs_s,                                                                                  &
     2419         bc_cs_t,                                                                                  &
     2420         call_chem_at_all_substeps,                                                                &
     2421         chem_debug0,                                                                              &
     2422         chem_debug1,                                                                              &
     2423         chem_debug2,                                                                              &
     2424         chem_gasphase_on,                                                                         &
     2425         chem_mechanism,                                                                           &
     2426         cs_heights,                                                                               &
     2427         cs_name,                                                                                  &
     2428         cs_profile,                                                                               &
     2429         cs_surface,                                                                               &
     2430         cs_surface_initial_change,                                                                &
     2431         cs_vertical_gradient_level,                                                               &
     2432         daytype_mdh,                                                                              &
     2433         deposition_dry,                                                                           &
     2434         emissions_anthropogenic,                                                                  &
     2435         emiss_lod,                                                                                &
     2436         emiss_factor_main,                                                                        &
     2437         emiss_factor_side,                                                                        &
     2438         emiss_read_legacy_mode,                                                                   &
     2439         icntrl,                                                                                   &
     2440         main_street_id,                                                                           &
     2441         max_street_id,                                                                            &
     2442         mode_emis,                                                                                &
     2443         my_steps,                                                                                 &
     2444         nesting_chem,                                                                             &
     2445         nesting_offline_chem,                                                                     &
     2446         rcntrl,                                                                                   &
     2447         side_street_id,                                                                           &
     2448         photolysis_scheme,                                                                        &
     2449         wall_csflux,                                                                              &
     2450         cs_vertical_gradient,                                                                     &
     2451         top_csflux,                                                                               &
     2452         surface_csflux,                                                                           &
     2453         surface_csflux_name,                                                                      &
    24412454         time_fac_type
    24422455!
    2443 !-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
     2456!-- Analogically to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
    24442457!-- so this way we could prescribe a specific flux value for each species
    24452458    !>  chemistry_parameters for initial profiles
    24462459    !>  cs_names = 'O3', 'NO2', 'NO', ...   to set initial profiles)
    24472460    !>  cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3)
    2448     !>  cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.) 
     2461    !>  cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.)
    24492462    !>  cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, .....  (chem spcs conc at height lvls chem_heights(1,:)) etc.
    24502463    !>  If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
    24512464    !>  then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
    24522465!
    2453 !--   Read chem namelist   
    2454 
     2466!-- Read chem namelist
    24552467    CHARACTER(LEN=8)    :: solver_type
    24562468
     
    24622474    rtol = 0.01_wp
    24632475!
    2464 !--   Try to find chemistry package
     2476!-- Try to find chemistry package
    24652477    REWIND ( 11 )
    24662478    line = ' '
     
    24702482    BACKSPACE ( 11 )
    24712483!
    2472 !--   Read chemistry namelist
    2473     READ ( 11, chemistry_parameters, ERR = 10, END = 20 )     
    2474 !
    2475 !--   Enable chemistry model
    2476     air_chemistry = .TRUE.                   
     2484!-- Read chemistry namelist
     2485    READ ( 11, chemistry_parameters, ERR = 10, END = 20 )
     2486!
     2487!-- Enable chemistry model
     2488    air_chemistry = .TRUE.
    24772489    GOTO 20
    24782490
     
    24862498
    24872499!
    2488 !-- synchronize emiss_lod and mod_emis only if emissions_anthropogenic
    2489 !-- is activated in the namelist.  Otherwise their values are "don't care"
     2500!-- Synchronize emiss_lod and mod_emis only if emissions_anthropogenic is activated in the namelist.
     2501!-- Otherwise their values are "don't care"
    24902502    IF ( emissions_anthropogenic )  THEN
    24912503
    24922504!
    2493 !--    check for emission mode for chem species
    2494 
     2505!--    Check for emission mode for chem species
    24952506       IF ( emiss_lod < 0 )  THEN   !- if LOD not defined in namelist
    2496           IF ( ( mode_emis /= 'PARAMETERIZED'  )    .AND.      &
    2497                ( mode_emis /= 'DEFAULT'        )    .AND.      &
     2507          IF ( ( mode_emis /= 'PARAMETERIZED'  )    .AND.                                          &
     2508               ( mode_emis /= 'DEFAULT'        )    .AND.                                          &
    24982509               ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
    24992510             message_string = 'Incorrect mode_emiss  option select. Please check spelling'
     
    25012512          ENDIF
    25022513       ELSE
    2503           IF ( ( emiss_lod /= 0 )    .AND.         &
    2504                ( emiss_lod /= 1 )    .AND.         &
     2514          IF ( ( emiss_lod /= 0 )    .AND.                                                         &
     2515               ( emiss_lod /= 1 )    .AND.                                                         &
    25052516               ( emiss_lod /= 2 ) )  THEN
    25062517             message_string = 'Invalid value for emiss_lod (0, 1, or 2)'
     
    25102521
    25112522!
    2512 ! for reference (ecc)
     2523! For reference (ecc)
    25132524!    IF ( (mode_emis /= 'PARAMETERIZED')  .AND. ( mode_emis /= 'DEFAULT' ) .AND. ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
    25142525!       message_string = 'Incorrect mode_emiss  option select. Please check spelling'
     
    25172528
    25182529!
    2519 !-- conflict resolution for emiss_lod and mode_emis
     2530!-- Conflict resolution for emiss_lod and mode_emis
    25202531!-- 1) if emiss_lod is defined, have mode_emis assume same setting as emiss_lod
    25212532!-- 2) if emiss_lod it not defined, have emiss_lod assuem same setting as mode_emis
    2522 !-- this check is in place to retain backward compatibility with salsa until the
    2523 !-- code is migrated completed to emiss_lod
    2524 !-- note that
    2525 
     2533!-- this check is in place to retain backward compatibility with salsa until the code is migrated
     2534!-- completed to emiss_lod
     2535!-- note that
    25262536       IF  ( emiss_lod >= 0 ) THEN
    25272537
     
    25342544                mode_emis = 'PRE-PROCESSED'
    25352545          END SELECT
    2536        
    2537           message_string = 'Synchronizing mode_emis to defined emiss_lod'               //  &
    2538                            CHAR(10)  //  '                    '                         //  &
    2539                            'NOTE - mode_emis will be depreciated in future releases'    //  &
    2540                            CHAR(10)  //  '                    '                         //  &
     2546
     2547          message_string = 'Synchronizing mode_emis to defined emiss_lod'               //         &
     2548                           CHAR( 10 )  //  '                    '                       //         &
     2549                           'NOTE - mode_emis will be depreciated in future releases'    //         &
     2550                           CHAR( 10 )  //  '                    '                       //         &
    25412551                           'please use emiss_lod to define emission mode'
    2542  
     2552
    25432553          CALL message ( 'parin_chem', 'CM0463', 0, 0, 0, 6, 0 )
    25442554
     
    25542564          END SELECT
    25552565
    2556           message_string = 'emiss_lod undefined.  Using existing mod_emiss setting'     //  &
    2557                            CHAR(10)  //  '                    '                         //  &
    2558                            'NOTE - mode_emis will be depreciated in future releases'    //  &
    2559                            CHAR(10)  //  '                    '                         //  &
     2566          message_string = 'emiss_lod undefined.  Using existing mod_emiss setting'     //         &
     2567                           CHAR( 10 )  //  '                    '                       //         &
     2568                           'NOTE - mode_emis will be depreciated in future releases'    //         &
     2569                           CHAR( 10 )  //  '                    '                       //         &
    25602570                           'please use emiss_lod to define emission mode'
    25612571
     
    25762586       IF ( emiss_read_legacy_mode )  THEN       !< notify legacy read mode
    25772587
    2578           message_string = 'Legacy emission read mode activated'            // &
    2579                            CHAR(10)  //  '                    '             // &
    2580                            'All emissions data will be loaded '             // &
     2588          message_string = 'Legacy emission read mode activated'            //                     &
     2589                           CHAR( 10 )  //  '                    '           //                    &
     2590                           'All emissions data will be loaded '             //                     &
    25812591                           'prior to start of simulation'
    25822592          CALL message ( 'parin_chem', 'CM0465', 0, 0, 0, 6, 0 )
    25832593
    2584        ELSE                                     !< if new read mode selected 
     2594       ELSE                                     !< if new read mode selected
    25852595
    25862596          IF ( emiss_lod < 2 )  THEN            !< check LOD compatibility
    25872597
    2588              message_string = 'New emission read mode '                    // &
    2589                               'currently unavailable for LODs 0 and 1.'    // &
    2590                               CHAR(10)  //  '                    '         // &
     2598             message_string = 'New emission read mode '                    //                      &
     2599                              'currently unavailable for LODs 0 and 1.'    //                      &
     2600                              CHAR( 10 )  //  '                    '       //                      &
    25912601                              'Reverting to legacy emission read mode'
    25922602             CALL message ( 'parin_chem', 'CM0466', 0, 0, 0, 6, 0 )
     
    25962606          ELSE                                  !< notify new read mode
    25972607
    2598              message_string = 'New emission read mode activated'           // &
    2599                               CHAR(10)  //  '                    '         // &
    2600                               'LOD 2 emissions will be updated on-demand ' // &           
     2608             message_string = 'New emission read mode activated'           //                      &
     2609                              CHAR( 10 )  //  '                    '       //                      &
     2610                              'LOD 2 emissions will be updated on-demand ' //                      &
    26012611                              'according to indicated timestamps'
    26022612             CALL message ( 'parin_chem', 'CM0467', 0, 0, 0, 6, 0 )
     
    26052615
    26062616       ENDIF ! if emiss_read_legacy_mode
    2607  
    2608                            
     2617
     2618
    26092619    ENDIF  ! if emissions_anthropengic
    26102620
    26112621
    2612     t_steps = my_steps         
    2613 !
    2614 !--    Determine the number of user-defined profiles and append them to the
    2615 !--    standard data output (data_output_pr)
    2616     max_pr_cs_tmp = 0 
     2622    t_steps = my_steps
     2623!
     2624!-- Determine the number of user-defined profiles and append them to the standard data output
     2625!>> (data_output_pr)
     2626    max_pr_cs_tmp = 0
    26172627    i = 1
    26182628
     
    26282638       max_pr_cs = max_pr_cs_tmp
    26292639    ENDIF
    2630 
    2631     !      Set Solver Type
     2640!
     2641!-- Set Solver Type
    26322642    IF(icntrl(3) == 0)  THEN
    26332643       solver_type = 'rodas3'           !Default
     
    26762686
    26772687
    2678 !------------------------------------------------------------------------------!
     2688!--------------------------------------------------------------------------------------------------!
    26792689! Description:
    26802690! ------------
    26812691!> Call for all grid points
    2682 !------------------------------------------------------------------------------!
     2692!--------------------------------------------------------------------------------------------------!
    26832693    SUBROUTINE chem_actions( location )
    26842694
     
    26912701!
    26922702!--       Chemical reactions and deposition
    2693           IF ( chem_gasphase_on ) THEN
     2703          IF ( chem_gasphase_on )  THEN
    26942704!
    26952705!--          If required, calculate photolysis frequencies -
     
    27092719
    27102720
    2711 !------------------------------------------------------------------------------!
     2721!--------------------------------------------------------------------------------------------------!
    27122722! Description:
    27132723! ------------
    27142724!> Call for grid points i,j
    2715 !------------------------------------------------------------------------------!
     2725!--------------------------------------------------------------------------------------------------!
    27162726
    27172727    SUBROUTINE chem_actions_ij( i, j, location )
    27182728
     2729    CHARACTER (LEN=*), INTENT(IN) ::  location  !< call location string
     2730
     2731    INTEGER(iwp)  ::  dummy                     !< call location string
    27192732
    27202733    INTEGER(iwp),      INTENT(IN) ::  i         !< grid index in x-direction
    27212734    INTEGER(iwp),      INTENT(IN) ::  j         !< grid index in y-direction
    2722     CHARACTER (LEN=*), INTENT(IN) ::  location  !< call location string
    2723     INTEGER(iwp)  ::  dummy  !< call location string
    27242735
    27252736    IF ( air_chemistry    )   dummy = i + j
     
    27362747
    27372748
    2738 !------------------------------------------------------------------------------!
     2749!--------------------------------------------------------------------------------------------------!
    27392750! Description:
    27402751! ------------
    27412752!> Call for all grid points
    2742 !------------------------------------------------------------------------------!
     2753!--------------------------------------------------------------------------------------------------!
    27432754    SUBROUTINE chem_non_advective_processes()
    27442755
     
    27472758      INTEGER(iwp) ::  j  !<
    27482759
    2749       !
    2750       !-- Calculation of chemical reactions and deposition.
    2751 
    2752 
    2753       IF ( intermediate_timestep_count == 1 .OR. call_chem_at_all_substeps )  THEN
    2754 
    2755          IF ( chem_gasphase_on ) THEN
     2760!
     2761!--   Calculation of chemical reactions and deposition.
     2762      IF ( intermediate_timestep_count == 1  .OR.  call_chem_at_all_substeps )  THEN
     2763
     2764         IF ( chem_gasphase_on )  THEN
    27562765            CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
    27572766            !$OMP PARALLEL PRIVATE (i,j)
     
    27832792
    27842793
    2785 !------------------------------------------------------------------------------!
     2794!--------------------------------------------------------------------------------------------------!
    27862795! Description:
    27872796! ------------
    27882797!> Call for grid points i,j
    2789 !------------------------------------------------------------------------------!
     2798!--------------------------------------------------------------------------------------------------!
    27902799 SUBROUTINE chem_non_advective_processes_ij( i, j )
    27912800
     
    27962805!
    27972806!-- Calculation of chemical reactions and deposition.
    2798 
    2799 
    2800    IF ( intermediate_timestep_count == 1 .OR. call_chem_at_all_substeps )  THEN
    2801 
    2802       IF ( chem_gasphase_on ) THEN
     2807   IF ( intermediate_timestep_count == 1  .OR.  call_chem_at_all_substeps )  THEN
     2808
     2809      IF ( chem_gasphase_on )  THEN
    28032810         CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
    28042811         CALL chem_integrate( i, j )
     
    28172824
    28182825 END SUBROUTINE chem_non_advective_processes_ij
    2819  
    2820 !------------------------------------------------------------------------------!
     2826
     2827!--------------------------------------------------------------------------------------------------!
    28212828! Description:
    28222829! ------------
    2823 !> routine for exchange horiz of chemical quantities 
    2824 !------------------------------------------------------------------------------!
     2830!> routine for exchange horiz of chemical quantities
     2831!--------------------------------------------------------------------------------------------------!
    28252832 SUBROUTINE chem_exchange_horiz_bounds
    2826  
    2827     USE exchange_horiz_mod,                                                    &
     2833
     2834    USE exchange_horiz_mod,                                                                        &
    28282835        ONLY:  exchange_horiz
    28292836
    28302837   INTEGER(iwp) ::  lsp       !<
    2831  
    2832 !
    2833 !--    Loop over chemical species       
     2838
     2839!
     2840!--    Loop over chemical species
    28342841       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' )
    28352842       DO  lsp = 1, nvar
     
    28452852 END SUBROUTINE chem_exchange_horiz_bounds
    28462853
    2847  
    2848 !------------------------------------------------------------------------------!
     2854
     2855!--------------------------------------------------------------------------------------------------!
    28492856! Description:
    28502857! ------------
    2851 !> Subroutine calculating prognostic equations for chemical species
    2852 !> (vector-optimized).
    2853 !> Routine is called separately for each chemical species over a loop from
    2854 !> prognostic_equations.
    2855 !------------------------------------------------------------------------------!
     2858!> Subroutine calculating prognostic equations for chemical species (vector-optimized).
     2859!> Routine is called separately for each chemical species over a loop from prognostic_equations.
     2860!--------------------------------------------------------------------------------------------------!
    28562861 SUBROUTINE chem_prognostic_equations()
    28572862
     
    29452950                   DO  k = nzb+1, nzt
    29462951                      chem_species(ilsp)%tconc_m(k,j,i) = - 9.5625_wp * tend(k,j,i)                &
    2947                            + 5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
     2952                                                     + 5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
    29482953                   ENDDO
    29492954                ENDDO
     
    29592964
    29602965
    2961 !------------------------------------------------------------------------------!
     2966!--------------------------------------------------------------------------------------------------!
    29622967! Description:
    29632968! ------------
    2964 !> Subroutine calculating prognostic equations for chemical species
    2965 !> (cache-optimized).
    2966 !> Routine is called separately for each chemical species over a loop from
    2967 !> prognostic_equations.
    2968 !------------------------------------------------------------------------------!
     2969!> Subroutine calculating prognostic equations for chemical species (cache-optimized).
     2970!> Routine is called separately for each chemical species over a loop from prognostic_equations.
     2971!--------------------------------------------------------------------------------------------------!
    29692972 SUBROUTINE chem_prognostic_equations_ij( i, j, i_omp_start, tn )
    29702973
    29712974
    29722975    INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn
     2976
    29732977    INTEGER(iwp) :: ilsp
    29742978!
     
    30093013!
    30103014!--    Diffusion terms (the last three arguments are zero)
    3011 
    30123015       CALL diffusion_s( i, j, chem_species(ilsp)%conc,                                            &
    30133016            surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:),                              &
     
    30333036               )
    30343037
    3035           IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp )  THEN               
    3036              chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)    !FKS6               
     3038          IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp )  THEN
     3039             chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)    !FKS6
    30373040          ENDIF
    30383041       ENDDO
     
    30483051             DO  k = nzb+1, nzt
    30493052                chem_species(ilsp)%tconc_m(k,j,i) = -9.5625_wp * tend(k,j,i) +                     &
    3050                      5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
     3053                                                     5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
    30513054             ENDDO
    30523055          ENDIF
     
    30583061
    30593062
    3060 !------------------------------------------------------------------------------!
     3063!--------------------------------------------------------------------------------------------------!
    30613064! Description:
    30623065! ------------
    30633066!> Read module-specific local restart data arrays (Fortran binary format).
    3064 !------------------------------------------------------------------------------!
    3065  SUBROUTINE chem_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,             &
    3066                                 nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc,   &
    3067                                 nys_on_file, tmp_3d, found )
     3067!--------------------------------------------------------------------------------------------------!
     3068 SUBROUTINE chem_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync,   &
     3069                                nyn_on_file, nysf, nysc, nys_on_file, tmp_3d, found )
    30683070
    30693071    USE control_parameters
    30703072
    3071 
     3073    INTEGER(iwp) ::  k               !<
    30723074    INTEGER(iwp) ::  lsp             !<
    3073     INTEGER(iwp) ::  k               !<
    3074     INTEGER(iwp) ::  nxlc            !<
    3075     INTEGER(iwp) ::  nxlf            !<
    3076     INTEGER(iwp) ::  nxl_on_file     !<   
    3077     INTEGER(iwp) ::  nxrc            !<
    3078     INTEGER(iwp) ::  nxrf            !<
    3079     INTEGER(iwp) ::  nxr_on_file     !<   
    3080     INTEGER(iwp) ::  nync            !<
    3081     INTEGER(iwp) ::  nynf            !<
    3082     INTEGER(iwp) ::  nyn_on_file     !<   
    3083     INTEGER(iwp) ::  nysc            !<
    3084     INTEGER(iwp) ::  nysf            !<
    3085     INTEGER(iwp) ::  nys_on_file     !<   
    3086 
    3087     LOGICAL, INTENT(OUT) :: found
    3088 
    3089     REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !< 3D array to temp store data
    3090 
    3091 
    3092     found = .FALSE. 
     3075    INTEGER(iwp) ::  nxlc            !<
     3076    INTEGER(iwp) ::  nxlf            !<
     3077    INTEGER(iwp) ::  nxl_on_file     !<
     3078    INTEGER(iwp) ::  nxrc            !<
     3079    INTEGER(iwp) ::  nxrf            !<
     3080    INTEGER(iwp) ::  nxr_on_file     !<
     3081    INTEGER(iwp) ::  nync            !<
     3082    INTEGER(iwp) ::  nynf            !<
     3083    INTEGER(iwp) ::  nyn_on_file     !<
     3084    INTEGER(iwp) ::  nysc            !<
     3085    INTEGER(iwp) ::  nysf            !<
     3086    INTEGER(iwp) ::  nys_on_file     !<
     3087
     3088    LOGICAL, INTENT(OUT) :: found
     3089
     3090    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) &
     3091                :: tmp_3d   !< 3D array to temp store data
     3092
     3093
     3094    found = .FALSE.
    30933095
    30943096
     
    31233125
    31243126
    3125 !------------------------------------------------------------------------------!
     3127!--------------------------------------------------------------------------------------------------!
    31263128! Description:
    31273129! ------------
    31283130!> Read module-specific local restart data arrays (Fortran binary format).
    3129 !------------------------------------------------------------------------------!
     3131!--------------------------------------------------------------------------------------------------!
    31303132 SUBROUTINE chem_rrd_local_mpi
    31313133
     
    31543156
    31553157
    3156 !-------------------------------------------------------------------------------!
     3158!--------------------------------------------------------------------------------------------------!
    31573159!> Description:
    31583160!> Calculation of horizontally averaged profiles
     
    31603162!> but at least for the region "total domain" (sr=0).
    31613163!> quantities.
    3162 !-------------------------------------------------------------------------------!
     3164!--------------------------------------------------------------------------------------------------!
    31633165 SUBROUTINE chem_statistics( mode, sr, tn )
    31643166
     
    31693171
    31703172
    3171     CHARACTER (LEN=*) ::  mode   !< 
     3173    CHARACTER (LEN=*) ::  mode   !<
    31723174
    31733175    INTEGER(iwp) ::  i    !< running index on x-axis
    31743176    INTEGER(iwp) ::  j    !< running index on y-axis
    31753177    INTEGER(iwp) ::  k    !< vertical index counter
     3178    INTEGER(iwp) ::  lpr  !< running index chem spcs
    31763179    INTEGER(iwp) ::  sr   !< statistical region
    3177     INTEGER(iwp) ::  tn   !< thread number
    3178     INTEGER(iwp) ::  lpr  !< running index chem spcs
     3180    INTEGER(iwp) ::  tn   !< thread number
    31793181
    31803182    IF ( mode == 'profiles' )  THEN
    31813183       !
    31823184!
    3183 !--    Sample on how to calculate horizontally averaged profiles of user-
    3184 !--    defined quantities. Each quantity is identified by the index
    3185 !--    "pr_palm+#" where "#" is an integer starting from 1. These
    3186 !--    user-profile-numbers must also be assigned to the respective strings
    3187 !--    given by data_output_pr_cs in routine user_check_data_output_pr.
     3185!--    Sample on how to calculate horizontally averaged profiles of user-defined quantities. Each
     3186!--    quantity is identified by the index "pr_palm+#" where "#" is an integer starting from 1.
     3187!--    These user-profile-numbers must also be assigned to the respective strings given by
     3188!--    data_output_pr_cs in routine user_check_data_output_pr.
    31883189!--    hom(:,:,:,:) =  dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
    31893190!--                     w*pt*), dim-4 = statistical region.
     
    31953196                DO lpr = 1, cs_pr_count
    31963197
    3197                    sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) +    &
    3198                         chem_species(cs_pr_index(lpr))%conc(k,j,i) *       &
    3199                         rmask(j,i,sr)  *                                   &
    3200                         MERGE( 1.0_wp, 0.0_wp,                             &
     3198                   sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) +   &
     3199                        chem_species(cs_pr_index(lpr))%conc(k,j,i) *                               &
     3200                        rmask(j,i,sr)  *                                                           &
     3201                        MERGE( 1.0_wp, 0.0_wp,                                                     &
    32013202                        BTEST( wall_flags_total_0(k,j,i), 22 ) )
    32023203                ENDDO
     
    32113212
    32123213
    3213 !------------------------------------------------------------------------------!
     3214!--------------------------------------------------------------------------------------------------!
    32143215! Description:
    32153216! ------------
    3216 !> Subroutine for swapping of timelevels for chemical species
    3217 !> called out from subroutine swap_timelevel
    3218 !------------------------------------------------------------------------------!
     3217!> Subroutine for swapping of timelevels for chemical species called out from subroutine
     3218!> swap_timelevel
     3219!--------------------------------------------------------------------------------------------------!
    32193220
    32203221
     
    32243225    INTEGER(iwp), INTENT(IN) ::  level
    32253226!
    3226 !-- local variables
     3227!-- Local variables
    32273228    INTEGER(iwp)             ::  lsp
    32283229
    32293230
    32303231    IF ( level == 0 )  THEN
    3231        DO  lsp=1, nvar                                       
     3232       DO  lsp=1, nvar
    32323233          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_1(:,:,:,lsp)
    32333234          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_2(:,:,:,lsp)
    32343235       ENDDO
    32353236    ELSE
    3236        DO  lsp=1, nvar                                       
     3237       DO  lsp=1, nvar
    32373238          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_2(:,:,:,lsp)
    32383239          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_1(:,:,:,lsp)
     
    32443245
    32453246
    3246 !------------------------------------------------------------------------------!
     3247!--------------------------------------------------------------------------------------------------!
    32473248! Description:
    32483249! ------------
    32493250!> Subroutine to write restart data for chemistry model
    3250 !------------------------------------------------------------------------------!
     3251!--------------------------------------------------------------------------------------------------!
    32513252 SUBROUTINE chem_wrd_local
    32523253
    32533254
    3254     INTEGER(iwp) ::  lsp  !< running index for chem spcs. 
     3255    INTEGER(iwp) ::  lsp  !< running index for chem spcs.
    32553256
    32563257    IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
     
    32793280
    32803281
    3281 !-------------------------------------------------------------------------------!
     3282!--------------------------------------------------------------------------------------------------!
    32823283! Description:
    32833284! ------------
    3284 !> Subroutine to calculate the deposition of gases and PMs. For now deposition
    3285 !> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT
    3286 !> considered. The deposition of particles is derived following Zhang et al.,
    3287 !> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010).
    3288 !>     
    3289 !> @TODO: Consider deposition on vertical surfaces   
     3285!> Subroutine to calculate the deposition of gases and PMs. For now deposition only takes place on
     3286!> lsm and usm horizontal surfaces. Default surfaces are NOT considered. The deposition of particles
     3287!> is derived following Zhang et al., 2001, gases are deposited using the DEPAC module
     3288!> (van Zanten et al., 2010).
     3289!>
     3290!> @TODO: Consider deposition on vertical surfaces
    32903291!> @TODO: Consider overlaying horizontal surfaces
    3291 !> @TODO: Consider resolved vegetation   
     3292!> @TODO: Consider resolved vegetation
    32923293!> @TODO: Check error messages
    3293 !-------------------------------------------------------------------------------!
     3294!--------------------------------------------------------------------------------------------------!
    32943295 SUBROUTINE chem_depo( i, j )
    32953296
    3296     USE control_parameters,                                                 &   
    3297          ONLY:  dt_3d, intermediate_timestep_count, latitude,               &
    3298                 time_since_reference_point
    3299 
    3300     USE arrays_3d,                                                          &
     3297    USE control_parameters,                                                                        &
     3298         ONLY:  dt_3d, intermediate_timestep_count, latitude, time_since_reference_point
     3299
     3300    USE arrays_3d,                                                                                 &
    33013301         ONLY:  dzw, rho_air_zw
    33023302
    3303     USE palm_date_time_mod,                                                 &
     3303    USE palm_date_time_mod,                                                                        &
    33043304         ONLY:  get_date_time
    33053305
    3306     USE surface_mod,                                                        &
    3307          ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_lsm_h,        &
    3308          surf_type, surf_usm_h
    3309 
    3310     USE radiation_model_mod,                                                &
     3306    USE surface_mod,                                                                               &
     3307         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_lsm_h, surf_type, surf_usm_h
     3308
     3309    USE radiation_model_mod,                                                                       &
    33113310         ONLY:  cos_zenith
    33123311
     
    33143313    INTEGER(iwp) ::  day_of_year                   !< current day of the year
    33153314    INTEGER(iwp), INTENT(IN) ::  i
     3315    INTEGER(iwp) ::  i_pspec                       !< index for matching depac gas component
    33163316    INTEGER(iwp), INTENT(IN) ::  j
    33173317    INTEGER(iwp) ::  k                             !< matching k to surface m at i,j
    33183318    INTEGER(iwp) ::  lsp                           !< running index for chem spcs.
    3319     INTEGER(iwp) ::  luv_palm                      !< index of PALM LSM vegetation_type at current surface element
    3320     INTEGER(iwp) ::  lup_palm                      !< index of PALM LSM pavement_type at current surface element
    3321     INTEGER(iwp) ::  luw_palm                      !< index of PALM LSM water_type at current surface element
    3322     INTEGER(iwp) ::  luu_palm                      !< index of PALM USM walls/roofs at current surface element
    3323     INTEGER(iwp) ::  lug_palm                      !< index of PALM USM green walls/roofs at current surface element
    3324     INTEGER(iwp) ::  lud_palm                      !< index of PALM USM windows at current surface element
     3319    INTEGER(iwp) ::  luv_palm                      !< index of PALM LSM vegetation_type at current
     3320                                                   !< surface element
     3321    INTEGER(iwp) ::  lup_palm                      !< index of PALM LSM pavement_type at current
     3322                                                   !< surface element
     3323    INTEGER(iwp) ::  luw_palm                      !< index of PALM LSM water_type at current
     3324                                                   !< surface element
     3325    INTEGER(iwp) ::  luu_palm                      !< index of PALM USM walls/roofs at current
     3326                                                   !< surface element
     3327    INTEGER(iwp) ::  lug_palm                      !< index of PALM USM green walls/roofs at current
     3328                                                   !< surface element
     3329    INTEGER(iwp) ::  lud_palm                      !< index of PALM USM windows at current surface
     3330                                                   !< element
    33253331    INTEGER(iwp) ::  luv_dep                       !< matching DEPAC LU to luv_palm
    33263332    INTEGER(iwp) ::  lup_dep                       !< matching DEPAC LU to lup_palm
     
    33303336    INTEGER(iwp) ::  lud_dep                       !< matching DEPAC LU to lud_palm
    33313337    INTEGER(iwp) ::  m                             !< index for horizontal surfaces
    3332 
    33333338    INTEGER(iwp) ::  pspec                         !< running index
    3334     INTEGER(iwp) ::  i_pspec                       !< index for matching depac gas component
    3335 !
    3336 !-- Vegetation                                               !< Assign PALM classes to DEPAC land use classes
    3337     INTEGER(iwp) ::  ind_luv_user = 0                        !<  ERROR as no class given in PALM 
     3339!
     3340!-- Vegetation                                               !<  Assign PALM classes to DEPAC land
     3341                                                             !< use classes
     3342    INTEGER(iwp) ::  ind_luv_user = 0                        !<  ERROR as no class given in PALM
    33383343    INTEGER(iwp) ::  ind_luv_b_soil = 1                      !<  assigned to ilu_desert
    33393344    INTEGER(iwp) ::  ind_luv_mixed_crops = 2                 !<  assigned to ilu_arable
     
    33523357    INTEGER(iwp) ::  ind_luv_ev_shrubs = 15                  !<  assigned to ilu_mediterrean_scrub
    33533358    INTEGER(iwp) ::  ind_luv_de_shrubs = 16                  !<  assigned to ilu_mediterrean_scrub
    3354     INTEGER(iwp) ::  ind_luv_mixed_forest = 17               !<  assigned to ilu_coniferous_forest (ave(decid+conif))
     3359    INTEGER(iwp) ::  ind_luv_mixed_forest = 17               !<  assigned to ilu_coniferous_forest(ave(decid+conif))
    33553360    INTEGER(iwp) ::  ind_luv_intrup_forest = 18              !<  assigned to ilu_other (ave(other+decid))
    33563361!
    33573362!-- Water
    3358     INTEGER(iwp) ::  ind_luw_user = 0                        !<  ERROR as no class given in PALM 
     3363    INTEGER(iwp) ::  ind_luw_user = 0                        !<  ERROR as no class given in PALM
    33593364    INTEGER(iwp) ::  ind_luw_lake = 1                        !<  assigned to ilu_water_inland
    33603365    INTEGER(iwp) ::  ind_luw_river = 2                       !<  assigned to ilu_water_inland
    33613366    INTEGER(iwp) ::  ind_luw_ocean = 3                       !<  assigned to ilu_water_sea
    33623367    INTEGER(iwp) ::  ind_luw_pond = 4                        !<  assigned to ilu_water_inland
    3363     INTEGER(iwp) ::  ind_luw_fountain = 5                    !<  assigned to ilu_water_inland 
     3368    INTEGER(iwp) ::  ind_luw_fountain = 5                    !<  assigned to ilu_water_inland
    33643369!
    33653370!-- Pavement
     
    33693374    INTEGER(iwp) ::  ind_lup_conc = 3                        !<  assigned to ilu_desert
    33703375    INTEGER(iwp) ::  ind_lup_sett = 4                        !<  assigned to ilu_desert
    3371     INTEGER(iwp) ::  ind_lup_pav_stones = 5                  !<  assigned to ilu_desert 
     3376    INTEGER(iwp) ::  ind_lup_pav_stones = 5                  !<  assigned to ilu_desert
    33723377    INTEGER(iwp) ::  ind_lup_cobblest = 6                    !<  assigned to ilu_desert
    33733378    INTEGER(iwp) ::  ind_lup_metal = 7                       !<  assigned to ilu_desert
    33743379    INTEGER(iwp) ::  ind_lup_wood = 8                        !<  assigned to ilu_desert
    3375     INTEGER(iwp) ::  ind_lup_gravel = 9                      !<  assigned to ilu_desert 
     3380    INTEGER(iwp) ::  ind_lup_gravel = 9                      !<  assigned to ilu_desert
    33763381    INTEGER(iwp) ::  ind_lup_f_gravel = 10                   !<  assigned to ilu_desert
    33773382    INTEGER(iwp) ::  ind_lup_pebblest = 11                   !<  assigned to ilu_desert
     
    33853390    INTEGER(iwp) ::  ind_p_dens = 2     !< index for rhopart in particle_pars
    33863391    INTEGER(iwp) ::  ind_p_slip = 3     !< index for slipcor in particle_pars
    3387 
     3392    INTEGER(iwp) ::  nwet               !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
    33883393    INTEGER(iwp) ::  part_type          !< index for particle type (PM10 or PM25) in particle_pars
    33893394
    3390     INTEGER(iwp) ::  nwet               !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
    3391 
    3392     REAL(wp) ::  dt_chem                !< length of chem time step
    3393     REAL(wp) ::  dh                     !< vertical grid size
    3394     REAL(wp) ::  inv_dh                 !< inverse of vertical grid size
    3395     REAL(wp) ::  dt_dh                  !< dt_chem/dh
    3396 
    3397     REAL(wp) ::  dens              !< density at layer k at i,j 
     3395
     3396    REAL(wp) ::  conc_ijk_ugm3     !< concentration at i, j, k in ug/m3
     3397    REAL(wp) ::  dens              !< density at layer k at i,j
     3398    REAL(wp) ::  dh                !< vertical grid size
     3399    REAL(wp) ::  diffusivity       !< diffusivity
     3400    REAL(wp) ::  dt_chem           !< length of chem time step
     3401    REAL(wp) ::  dt_dh             !< dt_chem/dh
     3402    REAL(wp) ::  inv_dh            !< inverse of vertical grid size
     3403    REAL(wp) ::  lai               !< leaf area index at current surface element
     3404    REAL(wp) ::  ppm2ugm3          !< conversion factor from ppm to ug/m3
     3405    REAL(wp) ::  qv_tmp            !< surface mixing ratio at current surface element
    33983406    REAL(wp) ::  r_aero_surf       !< aerodynamic resistance (s/m) at current surface element
     3407    REAL(wp) ::  rb                !< quasi-laminar boundary layer resistance (s/m)
     3408    REAL(wp) ::  rc_tot            !< total canopy resistance (s/m)
     3409    REAL(wp) ::  rh_surf           !< relative humidity at current surface element
     3410    REAL(wp) ::  rs                !< Sedimentaion resistance (s/m)
     3411    REAL(wp) ::  sai               !< surface area index at current surface element assumed to be
     3412                                   !< lai + 1
     3413    REAL(wp) ::  slinnfac
     3414    REAL(wp) ::  solar_rad         !< solar radiation, direct and diffuse, at current surface
     3415                                   !< element
     3416    REAL(wp) ::  temp_tmp          !< temperatur at i,j,k
     3417    REAL(wp) ::  ts                !< surface temperatur in degrees celsius
    33993418    REAL(wp) ::  ustar_surf        !< ustar at current surface element
    3400     REAL(wp) ::  z0h_surf          !< roughness length for heat at current surface element
    3401     REAL(wp) ::  solar_rad         !< solar radiation, direct and diffuse, at current surface element
    3402     REAL(wp) ::  ppm2ugm3          !< conversion factor from ppm to ug/m3
    3403     REAL(wp) ::  rh_surf           !< relative humidity at current surface element
    3404     REAL(wp) ::  lai               !< leaf area index at current surface element
    3405     REAL(wp) ::  sai               !< surface area index at current surface element assumed to be lai + 1
    3406 
    3407     REAL(wp) ::  slinnfac       
     3419    REAL(wp) ::  vd_lu             !< deposition velocity (m/s)
    34083420    REAL(wp) ::  visc              !< Viscosity
    34093421    REAL(wp) ::  vs                !< Sedimentation velocity
    3410     REAL(wp) ::  vd_lu             !< deposition velocity (m/s)
    3411     REAL(wp) ::  rs                !< Sedimentaion resistance (s/m)
    3412     REAL(wp) ::  rb                !< quasi-laminar boundary layer resistance (s/m)
    3413     REAL(wp) ::  rc_tot            !< total canopy resistance (s/m)
    3414 
    3415     REAL(wp) ::  conc_ijk_ugm3     !< concentration at i, j, k in ug/m3
    3416     REAL(wp) ::  diffusivity       !< diffusivity
    3417 
    3418 
    3419     REAL(wp), DIMENSION(nspec) ::  bud_luv      !< budget for LSM vegetation type at current surface element
    3420     REAL(wp), DIMENSION(nspec) ::  bud_lup      !< budget for LSM pavement type at current surface element
    3421     REAL(wp), DIMENSION(nspec) ::  bud_luw      !< budget for LSM water type at current surface element
    3422     REAL(wp), DIMENSION(nspec) ::  bud_luu      !< budget for USM walls/roofs at current surface element
    3423     REAL(wp), DIMENSION(nspec) ::  bud_lug      !< budget for USM green surfaces at current surface element
     3422    REAL(wp) ::  z0h_surf          !< roughness length for heat at current surface element
     3423
     3424    REAL(wp), DIMENSION(nspec) ::  bud          !< overall budget at current surface element
    34243425    REAL(wp), DIMENSION(nspec) ::  bud_lud      !< budget for USM windows at current surface element
    3425     REAL(wp), DIMENSION(nspec) ::  bud          !< overall budget at current surface element
     3426    REAL(wp), DIMENSION(nspec) ::  bud_lug      !< budget for USM green surfaces at current surface
     3427                                                !< element
     3428    REAL(wp), DIMENSION(nspec) ::  bud_lup      !< budget for LSM pavement type at current surface
     3429                                                !< element
     3430    REAL(wp), DIMENSION(nspec) ::  bud_luu      !< budget for USM walls/roofs at current surface
     3431                                                !< element
     3432    REAL(wp), DIMENSION(nspec) ::  bud_luv      !< budget for LSM vegetation type at current surface
     3433                                                !< element
     3434    REAL(wp), DIMENSION(nspec) ::  bud_luw      !< budget for LSM water type at current surface
     3435                                                !< element
     3436    REAL(wp), DIMENSION(nspec) ::  ccomp_tot    !< total compensation point (ug/m3), for now kept to
     3437                                                !< zero for all species!
    34263438    REAL(wp), DIMENSION(nspec) ::  conc_ijk     !< concentration at i,j,k
    3427     REAL(wp), DIMENSION(nspec) ::  ccomp_tot    !< total compensation point (ug/m3), for now kept to zero for all species!
    3428                                                  
    3429 
    3430     REAL(wp) ::  temp_tmp       !< temperatur at i,j,k
    3431     REAL(wp) ::  ts             !< surface temperatur in degrees celsius
    3432     REAL(wp) ::  qv_tmp         !< surface mixing ratio at current surface element
    3433 !
    3434 !-- Particle parameters (PM10 (1), PM25 (2))
    3435 !-- partsize (diameter in m), rhopart (density in kg/m3), slipcor
    3436 !-- (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3)
    3437     REAL(wp), DIMENSION(1:3,1:2), PARAMETER ::  particle_pars = RESHAPE( (/ &
    3438          8.0e-6_wp, 1.14e3_wp, 1.016_wp, &  !<  1
    3439          0.7e-6_wp, 1.14e3_wp, 1.082_wp &   !<  2
    3440          /), (/ 3, 2 /) )
    3441 
     3439
     3440
     3441!
     3442!-- Particle parameters (PM10 (1), PM25 (2)) partsize (diameter in m), rhopart (density in kg/m3),
     3443!-- slipcor (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3)
    34423444    LOGICAL ::  match_lsm     !< flag indicating natural-type surface
    34433445    LOGICAL ::  match_usm     !< flag indicating urban-type surface
     3446
     3447    REAL(wp), DIMENSION(1:3,1:2), PARAMETER ::  particle_pars = RESHAPE( (/                        &
     3448                                                         8.0e-6_wp, 1.14e3_wp, 1.016_wp, &   !<  1
     3449                                                         0.7e-6_wp, 1.14e3_wp, 1.082_wp  &   !<  2
     3450                                                         /), (/ 3, 2 /) )
    34443451!
    34453452!-- List of names of possible tracers
    3446     CHARACTER(LEN=*), PARAMETER ::  pspecnames(nposp) = (/ &
     3453    CHARACTER(LEN=*), PARAMETER ::  pspecnames(nposp) = (/                                         &
    34473454         'NO2           ', &    !< NO2
    34483455         'NO            ', &    !< NO
     
    35153522         'O3_biascorr   ' /)    !< o3_biascorr
    35163523!
    3517 !-- tracer mole mass:
    3518     REAL(wp), PARAMETER ::  specmolm(nposp) = (/ &
     3524!-- Tracer mole mass:
     3525    REAL(wp), PARAMETER ::  specmolm(nposp) = (/                                                   &
    35193526         xm_O * 2 + xm_N, &                         !< NO2
    35203527         xm_O + xm_N, &                             !< NO
     
    36003607!
    36013608!--For LSM surfaces
    3602 
    36033609    IF ( match_lsm )  THEN
    36043610!
     
    36123618       r_aero_surf = surf_lsm_h%r_a(m)
    36133619       solar_rad   = surf_lsm_h%rad_sw_dir(m) + surf_lsm_h%rad_sw_dif(m)
    3614        
     3620
    36153621       lai = surf_lsm_h%lai(m)
    36163622       sai = lai + 1
     
    36283634       luw_palm = 0
    36293635       luw_dep = 0
    3630 ! 
     3636!
    36313637!--    Initialize budgets
    36323638       bud_luv = 0.0_wp
     
    36753681             luv_dep = 4
    36763682          ELSEIF ( luv_palm == ind_luv_intrup_forest )  THEN
    3677              luv_dep = 8     
     3683             luv_dep = 8
    36783684          ENDIF
    36793685       ENDIF
     
    36953701             lup_dep = 9
    36963702          ELSEIF ( lup_palm == ind_lup_cobblest )  THEN
    3697              lup_dep = 9       
     3703             lup_dep = 9
    36983704          ELSEIF ( lup_palm == ind_lup_metal )  THEN
    36993705             lup_dep = 9
    37003706          ELSEIF ( lup_palm == ind_lup_wood )  THEN
    3701              lup_dep = 9   
     3707             lup_dep = 9
    37023708          ELSEIF ( lup_palm == ind_lup_gravel )  THEN
    37033709             lup_dep = 9
     
    37183724
    37193725       IF ( surf_lsm_h%frac(m,ind_wat_win) > 0 )  THEN
    3720           luw_palm = surf_lsm_h%water_type(m)     
     3726          luw_palm = surf_lsm_h%water_type(m)
    37213727          IF ( luw_palm == ind_luw_user )  THEN
    37223728             message_string = 'No water type defined. Please define water type to enable deposition calculation'
     
    37313737             luw_dep = 13
    37323738          ELSEIF ( luw_palm == ind_luw_fountain )  THEN
    3733              luw_dep = 13 
     3739             luw_dep = 13
    37343740          ENDIF
    37353741       ENDIF
     
    37423748       ENDIF
    37433749!
    3744 !--    Compute length of time step 
     3750!--    Compute length of time step
    37453751       IF ( call_chem_at_all_substeps )  THEN
    37463752          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
     
    37703776!
    37713777!--    Calculate relative humidity from specific humidity for DEPAC
    3772        qv_tmp = MAX(q(k,j,i),0.0_wp)
     3778       qv_tmp = MAX( q(k,j,i), 0.0_wp)
    37733779       rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
    37743780!
     
    37813787!
    37823788!--       No vegetation on bare soil, desert or ice:
    3783           IF ( ( luv_palm == ind_luv_b_soil ) .OR. &
    3784                 ( luv_palm == ind_luv_desert ) .OR. &
    3785                  ( luv_palm == ind_luv_ice ) ) THEN
     3789          IF ( ( luv_palm == ind_luv_b_soil )  .OR.  ( luv_palm == ind_luv_desert )  .OR.          &
     3790               ( luv_palm == ind_luv_ice ) ) THEN
    37863791
    37873792             lai = 0.0_wp
     
    37893794
    37903795          ENDIF
    3791          
     3796
    37923797          slinnfac = 1.0_wp
    37933798!
     
    37953800          DO  lsp = 1, nvar
    37963801!
    3797 !--          Initialize 
     3802!--          Initialize
    37983803             vs = 0.0_wp
    37993804             vd_lu = 0.0_wp
     
    38053810!
    38063811!--             Sedimentation velocity
    3807                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    3808                      particle_pars(ind_p_size, part_type), &
    3809                      particle_pars(ind_p_slip, part_type), &
    3810                      visc)
    3811 
    3812                 CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    3813                      vs, &
    3814                      particle_pars(ind_p_size, part_type), &
    3815                      particle_pars(ind_p_slip, part_type), &
    3816                      nwet, temp_tmp, dens, visc, &
    3817                      luv_dep,  &
    3818                      r_aero_surf, ustar_surf )
    3819 
    3820                 bud_luv(lsp) = - conc_ijk(lsp) * &
    3821                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3812                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
     3813                                                        particle_pars(ind_p_size, part_type),      &
     3814                                                        particle_pars(ind_p_slip, part_type),      &
     3815                                                        visc)
     3816
     3817                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
     3818                                            vs,                                                    &
     3819                                            particle_pars(ind_p_size, part_type),                  &
     3820                                            particle_pars(ind_p_slip, part_type),                  &
     3821                                            nwet, temp_tmp, dens, visc,                            &
     3822                                            luv_dep,                                               &
     3823                                            r_aero_surf, ustar_surf )
     3824
     3825                bud_luv(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
    38223826
    38233827
     
    38263830!
    38273831!--             Sedimentation velocity
    3828                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    3829                      particle_pars(ind_p_size, part_type), &
    3830                      particle_pars(ind_p_slip, part_type), &
    3831                      visc)
    3832 
    3833                 CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    3834                      vs, &
    3835                      particle_pars(ind_p_size, part_type), &
    3836                      particle_pars(ind_p_slip, part_type), &
    3837                      nwet, temp_tmp, dens, visc, &
    3838                      luv_dep , &
    3839                      r_aero_surf, ustar_surf )
    3840 
    3841                 bud_luv(lsp) = - conc_ijk(lsp) * &
    3842                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3832                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
     3833                                                        particle_pars(ind_p_size, part_type),      &
     3834                                                        particle_pars(ind_p_slip, part_type),      &
     3835                                                        visc)
     3836
     3837                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
     3838                                            vs,                                                    &
     3839                                            particle_pars(ind_p_size, part_type),                  &
     3840                                            particle_pars(ind_p_slip, part_type),                  &
     3841                                            nwet, temp_tmp, dens, visc,                            &
     3842                                            luv_dep ,                                              &
     3843                                            r_aero_surf, ustar_surf )
     3844
     3845                bud_luv(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
    38433846
    38443847             ELSE !< GASES
     
    38663869!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
    38673870!--             Use density at k:
    3868 
    38693871                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
    38703872!
     
    38743876!
    38753877!--             Diffusivity for DEPAC relevant gases
    3876 !--             Use default value 
     3878!--             Use default value
    38773879                diffusivity            = 0.11e-4
    38783880!
    3879 !--             overwrite with known coefficients of diffusivity from Massman (1998)
    3880                 IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
     3881!--             Overwrite with known coefficients of diffusivity from Massman (1998)
     3882                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
    38813883                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
    38823884                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
     
    38873889!
    38883890!--             Get quasi-laminar boundary layer resistance rb:
    3889                 CALL get_rb_cell( (luv_dep == ilu_water_sea) .OR. (luv_dep == ilu_water_inland), &
    3890                      z0h_surf, ustar_surf, diffusivity, &
    3891                      rb )
     3891                CALL get_rb_cell( ( luv_dep == ilu_water_sea )  .OR.                              &
     3892                                  ( luv_dep == ilu_water_inland ), z0h_surf, ustar_surf,          &
     3893                                  diffusivity, rb )
    38923894!
    38933895!--             Get rc_tot
    3894                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf, solar_rad, cos_zenith, &
    3895                      rh_surf, lai, sai, nwet, luv_dep, 2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity, &
    3896                      r_aero_surf , rb )
     3896                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
     3897                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luv_dep,  &
     3898                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
     3899                                         diffusivity, r_aero_surf , rb )
    38973900!
    38983901!--             Calculate budget
     
    39033906                ELSE
    39043907
    3905                    vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
    3906 
    3907                    bud_luv(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
    3908                         (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3908                   vd_lu = 1.0_wp / ( r_aero_surf + rb + rc_tot )
     3909
     3910                   bud_luv(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                          &
     3911                                    ( 1.0_wp - EXP( -vd_lu * dt_dh ) ) * dh
    39093912                ENDIF
    39103913
     
    39353938!
    39363939!--             Sedimentation velocity:
    3937                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    3938                      particle_pars(ind_p_size, part_type), &
    3939                      particle_pars(ind_p_slip, part_type), &
    3940                      visc)
    3941 
    3942                 CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    3943                      vs, &
    3944                      particle_pars(ind_p_size, part_type), &
    3945                      particle_pars(ind_p_slip, part_type), &
    3946                      nwet, temp_tmp, dens, visc, &
    3947                      lup_dep,  &
    3948                      r_aero_surf, ustar_surf )
    3949 
    3950                 bud_lup(lsp) = - conc_ijk(lsp) * &
    3951                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3940                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
     3941                                                        particle_pars(ind_p_size, part_type),      &
     3942                                                        particle_pars(ind_p_slip, part_type),      &
     3943                                                        visc)
     3944
     3945                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
     3946                                            vs,                                                    &
     3947                                            particle_pars(ind_p_size, part_type),                  &
     3948                                            particle_pars(ind_p_slip, part_type),                  &
     3949                                            nwet, temp_tmp, dens, visc,                            &
     3950                                            lup_dep,                                               &
     3951                                            r_aero_surf, ustar_surf )
     3952
     3953                bud_lup(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
    39523954
    39533955
     
    39563958!
    39573959!--             Sedimentation velocity:
    3958                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    3959                      particle_pars(ind_p_size, part_type), &
    3960                      particle_pars(ind_p_slip, part_type), &
    3961                      visc)
    3962 
    3963                 CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    3964                      vs, &
    3965                      particle_pars(ind_p_size, part_type), &
    3966                      particle_pars(ind_p_slip, part_type), &
    3967                      nwet, temp_tmp, dens, visc, &
    3968                      lup_dep, &
    3969                      r_aero_surf, ustar_surf )
    3970 
    3971                 bud_lup(lsp) = - conc_ijk(lsp) * &
    3972                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3960                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
     3961                                                        particle_pars(ind_p_size, part_type),      &
     3962                                                        particle_pars(ind_p_slip, part_type),      &
     3963                                                        visc)
     3964
     3965                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
     3966                                            vs,                                                    &
     3967                                            particle_pars(ind_p_size, part_type),                  &
     3968                                            particle_pars(ind_p_slip, part_type),                  &
     3969                                            nwet, temp_tmp, dens, visc,                            &
     3970                                            lup_dep,                                               &
     3971                                            r_aero_surf, ustar_surf )
     3972
     3973                bud_lup(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
    39733974
    39743975             ELSE  !<GASES
    39753976!
    39763977!--             Read spc_name of current species for gas parameter
    3977 
     3978                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
     3979                   i_pspec = 0
     3980                   DO  pspec = 1, nposp
     3981                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
     3982                         i_pspec = pspec
     3983                      END IF
     3984                   ENDDO
     3985
     3986                ELSE
     3987!
     3988!--                For now species not deposited
     3989                   CYCLE
     3990                ENDIF
     3991!
     3992!--             Factor used for conversion from ppb to ug/m3 :
     3993!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     3994!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     3995!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
     3996!--             thus:
     3997!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     3998!--             Use density at lowest layer:
     3999                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
     4000!
     4001!--             Atmospheric concentration in DEPAC is requested in ug/m3:
     4002                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     4003                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
     4004!
     4005!--             Diffusivity for DEPAC relevant gases
     4006!--             Use default value
     4007                diffusivity            = 0.11e-4
     4008!
     4009!--             Overwrite with known coefficients of diffusivity from Massman (1998)
     4010                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
     4011                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
     4012                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
     4013                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
     4014                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
     4015                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
     4016                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
     4017!
     4018!--             Get quasi-laminar boundary layer resistance rb:
     4019                CALL get_rb_cell( ( lup_dep == ilu_water_sea )  .OR.                               &
     4020                                  ( lup_dep == ilu_water_inland ), z0h_surf, ustar_surf,           &
     4021                                    diffusivity, rb )
     4022!
     4023!--             Get rc_tot
     4024                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
     4025                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lup_dep,  &
     4026                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
     4027                                         diffusivity, r_aero_surf , rb )
     4028!
     4029!--             Calculate budget
     4030                IF ( rc_tot <= 0.0 )  THEN
     4031                   bud_lup(lsp) = 0.0_wp
     4032                ELSE
     4033                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
     4034                   bud_lup(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                           &
     4035                                  ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
     4036                ENDIF
     4037
     4038             ENDIF
     4039          ENDDO
     4040       ENDIF
     4041!
     4042!--    Water
     4043       IF ( surf_lsm_h%frac(m,ind_wat_win) > 0 )  THEN
     4044!
     4045!--       No vegetation on water:
     4046          lai = 0.0_wp
     4047          sai = 0.0_wp
     4048          slinnfac = 1.0_wp
     4049!
     4050!--       Get vd
     4051          DO  lsp = 1, nvar
     4052!
     4053!--          Initialize
     4054             vs = 0.0_wp
     4055             vd_lu = 0.0_wp
     4056             rs = 0.0_wp
     4057             rb = 0.0_wp
     4058             rc_tot = 0.0_wp
     4059             IF ( spc_names(lsp) == 'PM10' )  THEN
     4060                part_type = 1
     4061!
     4062!--             Sedimentation velocity:
     4063                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
     4064                                                        particle_pars(ind_p_size, part_type),      &
     4065                                                        particle_pars(ind_p_slip, part_type),      &
     4066                                                        visc)
     4067
     4068                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
     4069                                            vs,                                                    &
     4070                                            particle_pars(ind_p_size, part_type),                  &
     4071                                            particle_pars(ind_p_slip, part_type),                  &
     4072                                            nwet, temp_tmp, dens, visc,                            &
     4073                                            luw_dep,                                               &
     4074                                            r_aero_surf, ustar_surf )
     4075
     4076                bud_luw(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
     4077
     4078             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
     4079                part_type = 2
     4080!
     4081!--             Sedimentation velocity:
     4082                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
     4083                                                        particle_pars(ind_p_size, part_type), &
     4084                                                        particle_pars(ind_p_slip, part_type), &
     4085                                                        visc)
     4086
     4087                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
     4088                                            vs,                                                    &
     4089                                            particle_pars(ind_p_size, part_type),                  &
     4090                                            particle_pars(ind_p_slip, part_type),                  &
     4091                                            nwet, temp_tmp, dens, visc,                            &
     4092                                            luw_dep,                                               &
     4093                                            r_aero_surf, ustar_surf )
     4094
     4095                bud_luw(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
     4096
     4097             ELSE  !<GASES
     4098!
     4099!--             Read spc_name of current species for gas PARAMETER
    39784100                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
    39794101                   i_pspec = 0
     
    39974119!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
    39984120!--             Use density at lowest layer:
    3999 
    4000                 ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
     4121                ppm2ugm3 = (dens/xm_air) * 0.001_wp  !< (mole air)/m3
    40014122!
    40024123!--             Atmospheric concentration in DEPAC is requested in ug/m3:
    4003                 !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    4004                 conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
     4124!--                 ug/m3        ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     4125                conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 *  specmolm(i_pspec)  ! in ug/m3
    40054126!
    40064127!--             Diffusivity for DEPAC relevant gases
    4007 !--             Use default value 
     4128!--             Use default value
    40084129                diffusivity            = 0.11e-4
    40094130!
    4010 !--             overwrite with known coefficients of diffusivity from Massman (1998)
    4011                 IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
     4131!--             Overwrite with known coefficients of diffusivity from Massman (1998)
     4132                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
    40124133                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
    40134134                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
     
    40184139!
    40194140!--             Get quasi-laminar boundary layer resistance rb:
    4020                 CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland),   &
    4021                      z0h_surf, ustar_surf, diffusivity, rb )
    4022 !
     4141                CALL get_rb_cell( ( luw_dep == ilu_water_sea )  .OR.                               &
     4142                                  ( luw_dep == ilu_water_inland ), z0h_surf, ustar_surf,           &
     4143                                    diffusivity, rb )
     4144
    40234145!--             Get rc_tot
    4024                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,      &
    4025                                          solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lup_dep, 2,    &
    4026                                          rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,              &
    4027                                          r_aero_surf , rb )
     4146                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
     4147                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luw_dep,  &
     4148                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
     4149                                         diffusivity, r_aero_surf , rb )
    40284150!
    40294151!--             Calculate budget
    40304152                IF ( rc_tot <= 0.0 )  THEN
    4031                    bud_lup(lsp) = 0.0_wp
     4153
     4154                   bud_luw(lsp) = 0.0_wp
     4155
    40324156                ELSE
     4157
    40334158                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
    4034                    bud_lup(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
    4035                         (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     4159
     4160                   bud_luw(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                           &
     4161                                  ( 1.0_wp - EXP(-vd_lu * dt_dh ) ) * dh
    40364162                ENDIF
    40374163
     
    40394165          ENDDO
    40404166       ENDIF
    4041 !
    4042 !--    Water
    4043        IF ( surf_lsm_h%frac(m,ind_wat_win) > 0 )  THEN
    4044 !
    4045 !--       No vegetation on water:
    4046           lai = 0.0_wp
    4047           sai = 0.0_wp
    4048           slinnfac = 1.0_wp
    4049 !
    4050 !--       Get vd
    4051           DO  lsp = 1, nvar
    4052 !
    4053 !--          Initialize
    4054              vs = 0.0_wp
    4055              vd_lu = 0.0_wp
    4056              rs = 0.0_wp
    4057              rb = 0.0_wp
    4058              rc_tot = 0.0_wp
    4059              IF ( spc_names(lsp) == 'PM10' )  THEN
    4060                 part_type = 1
    4061 !
    4062 !--             Sedimentation velocity:
    4063                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    4064                      particle_pars(ind_p_size, part_type), &
    4065                      particle_pars(ind_p_slip, part_type), &
    4066                      visc)
    4067 
    4068                 CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    4069                      vs, &
    4070                      particle_pars(ind_p_size, part_type), &
    4071                      particle_pars(ind_p_slip, part_type), &
    4072                      nwet, temp_tmp, dens, visc, &
    4073                      luw_dep, &
    4074                      r_aero_surf, ustar_surf )
    4075 
    4076                 bud_luw(lsp) = - conc_ijk(lsp) * &
    4077                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    4078 
    4079              ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
    4080                 part_type = 2
    4081 !
    4082 !--             Sedimentation velocity:
    4083                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    4084                      particle_pars(ind_p_size, part_type), &
    4085                      particle_pars(ind_p_slip, part_type), &
    4086                      visc)
    4087 
    4088                 CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    4089                      vs, &
    4090                      particle_pars(ind_p_size, part_type), &
    4091                      particle_pars(ind_p_slip, part_type), &
    4092                      nwet, temp_tmp, dens, visc, &
    4093                      luw_dep, &
    4094                      r_aero_surf, ustar_surf )
    4095 
    4096                 bud_luw(lsp) = - conc_ijk(lsp) * &
    4097                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    4098 
    4099              ELSE  !<GASES
    4100 !
    4101 !--             Read spc_name of current species for gas PARAMETER
    4102 
    4103                 IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
    4104                    i_pspec = 0
    4105                    DO  pspec = 1, nposp
    4106                       IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
    4107                          i_pspec = pspec
    4108                       END IF
    4109                    ENDDO
    4110 
    4111                 ELSE
    4112 !
    4113 !--                For now species not deposited
    4114                    CYCLE
    4115                 ENDIF
    4116 !
    4117 !--             Factor used for conversion from ppb to ug/m3 :
    4118 !--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
    4119 !--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
    4120 !--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
    4121 !--             thus:
    4122 !--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
    4123 !--             Use density at lowest layer:
    4124 
    4125                 ppm2ugm3 = (dens/xm_air) * 0.001_wp  !< (mole air)/m3
    4126 !
    4127 !--             Atmospheric concentration in DEPAC is requested in ug/m3:
    4128 !--                 ug/m3        ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    4129                 conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 *  specmolm(i_pspec)  ! in ug/m3
    4130 !
    4131 !--             Diffusivity for DEPAC relevant gases
    4132 !--             Use default value
    4133                 diffusivity            = 0.11e-4
    4134 !
    4135 !--             overwrite with known coefficients of diffusivity from Massman (1998)
    4136                 IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
    4137                 IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
    4138                 IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
    4139                 IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
    4140                 IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
    4141                 IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
    4142                 IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
    4143 !
    4144 !--             Get quasi-laminar boundary layer resistance rb:
    4145                 CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland),  &
    4146                      z0h_surf, ustar_surf, diffusivity, rb )
    4147 
    4148 !--             Get rc_tot
    4149                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         &
    4150                                          solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luw_dep, 2,    &
    4151                                          rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
    4152                                          r_aero_surf , rb )
    4153 !
    4154 !--             Calculate budget
    4155                 IF ( rc_tot <= 0.0 )  THEN
    4156 
    4157                    bud_luw(lsp) = 0.0_wp
    4158 
    4159                 ELSE
    4160 
    4161                    vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
    4162 
    4163                    bud_luw(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
    4164                         (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    4165                 ENDIF
    4166 
    4167              ENDIF
    4168           ENDDO
    4169        ENDIF
    41704167
    41714168
     
    41754172       DO  lsp = 1, nspec
    41764173
    4177           bud(lsp) = surf_lsm_h%frac(m,ind_veg_wall) * bud_luv(lsp) + &
    4178                surf_lsm_h%frac(m,ind_pav_green) * bud_lup(lsp) + &
    4179                surf_lsm_h%frac(m,ind_wat_win) * bud_luw(lsp)
     4174          bud(lsp) = surf_lsm_h%frac(m,ind_veg_wall)  * bud_luv(lsp) +                            &
     4175                     surf_lsm_h%frac(m,ind_pav_green) * bud_lup(lsp) +                            &
     4176                     surf_lsm_h%frac(m,ind_wat_win)  * bud_luw(lsp)
    41804177!
    41814178!--       Compute new concentration:
    41824179          conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
    41834180
    4184           chem_species(lsp)%conc(k,j,i) = MAX(0.0_wp, conc_ijk(lsp))
     4181          chem_species(lsp)%conc(k,j,i) = MAX( 0.0_wp, conc_ijk(lsp) )
    41854182
    41864183       ENDDO
     
    41884185    ENDIF
    41894186!
    4190 !-- For USM surfaces   
    4191 
     4187!-- For USM surfaces
    41924188    IF ( match_usm )  THEN
    41934189!
     
    42254221       IF ( surf_usm_h%frac(m,ind_pav_green) > 0 )  THEN
    42264222!
    4227 !--       For green urban surfaces (e.g. green roofs
    4228 !--       assume LU short grass
     4223!--       For green urban surfaces (e.g. green roofs assume LU short grass
    42294224          lug_palm = ind_luv_s_grass
    42304225          IF ( lug_palm == ind_luv_user )  THEN
     
    42664261             lug_dep = 4
    42674262          ELSEIF ( lug_palm == ind_luv_intrup_forest )  THEN
    4268              lug_dep = 8     
     4263             lug_dep = 8
    42694264          ENDIF
    42704265       ENDIF
     
    42734268!
    42744269!--       For walls in USM assume concrete walls/roofs,
    4275 !--       assumed LU class desert as also assumed for
    4276 !--       pavements in LSM
     4270!--       assumed LU class desert as also assumed for pavements in LSM
    42774271          luu_palm = ind_lup_conc
    42784272          IF ( luu_palm == ind_lup_user )  THEN
     
    42904284             luu_dep = 9
    42914285          ELSEIF ( luu_palm == ind_lup_cobblest )  THEN
    4292              luu_dep = 9       
     4286             luu_dep = 9
    42934287          ELSEIF ( luu_palm == ind_lup_metal )  THEN
    42944288             luu_dep = 9
    42954289          ELSEIF ( luu_palm == ind_lup_wood )  THEN
    4296              luu_dep = 9   
     4290             luu_dep = 9
    42974291          ELSEIF ( luu_palm == ind_lup_gravel )  THEN
    42984292             luu_dep = 9
     
    43144308       IF ( surf_usm_h%frac(m,ind_wat_win) > 0 )  THEN
    43154309!
    4316 !--       For windows in USM assume metal as this is
    4317 !--       as close as we get, assumed LU class desert
    4318 !--       as also assumed for pavements in LSM
    4319           lud_palm = ind_lup_metal     
     4310!--       For windows in USM assume metal as this is as close as we get, assumed LU class desert as
     4311!--       also assumed for pavements in LSM.
     4312          lud_palm = ind_lup_metal
    43204313          IF ( lud_palm == ind_lup_user )  THEN
    43214314             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
     
    43324325             lud_dep = 9
    43334326          ELSEIF ( lud_palm == ind_lup_cobblest )  THEN
    4334              lud_dep = 9       
     4327             lud_dep = 9
    43354328          ELSEIF ( lud_palm == ind_lup_metal )  THEN
    43364329             lud_dep = 9
    43374330          ELSEIF ( lud_palm == ind_lup_wood )  THEN
    4338              lud_dep = 9   
     4331             lud_dep = 9
    43394332          ELSEIF ( lud_palm == ind_lup_gravel )  THEN
    43404333             lud_dep = 9
     
    43624355       !ENDIF
    43634356!
    4364 !--    Compute length of time step 
     4357!--    Compute length of time step
    43654358       IF ( call_chem_at_all_substeps )  THEN
    43664359          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
     
    43714364       dh = dzw(k)
    43724365       inv_dh = 1.0_wp / dh
    4373        dt_dh = dt_chem/dh
     4366       dt_dh = dt_chem / dh
    43744367!
    43754368!--    Concentration at i,j,k
     
    43844377!
    43854378!--    Viscosity of air
    4386        visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
     4379       visc = 1.496e-6 * temp_tmp**1.5 / ( temp_tmp + 120.0 )
    43874380!
    43884381!--    Air density at k
     
    43904383!
    43914384!--    Calculate relative humidity from specific humidity for DEPAC
    4392        qv_tmp = MAX(q(k,j,i),0.0_wp)
    4393        rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
    4394 !
    4395 !--    Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
    4396 !--    for each surface fraction. Then derive overall budget taking into account the surface fractions.
    4397 !
     4385       qv_tmp = MAX( q(k,j,i), 0.0_wp )
     4386       rh_surf = relativehumidity_from_specifichumidity( qv_tmp, temp_tmp, hyp(k) )
     4387!
     4388!--    Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget for
     4389!--    each surface fraction. Then derive overall budget taking into account the surface fractions.
     4390
    43984391!--    Walls/roofs
    43994392       IF ( surf_usm_h%frac(m,ind_veg_wall) > 0 )  THEN
     
    44084401          DO  lsp = 1, nvar
    44094402!
    4410 !--          Initialize 
     4403!--          Initialize
    44114404             vs = 0.0_wp
    44124405             vd_lu = 0.0_wp
     
    44184411!
    44194412!--             Sedimentation velocity
    4420                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    4421                      particle_pars(ind_p_size, part_type), &
    4422                      particle_pars(ind_p_slip, part_type), &
    4423                      visc)
    4424 
    4425                 CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    4426                      vs, &
    4427                      particle_pars(ind_p_size, part_type), &
    4428                      particle_pars(ind_p_slip, part_type), &
    4429                      nwet, temp_tmp, dens, visc, &
    4430                      luu_dep,  &
    4431                      r_aero_surf, ustar_surf )
    4432 
    4433                 bud_luu(lsp) = - conc_ijk(lsp) * &
    4434                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     4413                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
     4414                                                        particle_pars(ind_p_size, part_type),      &
     4415                                                        particle_pars(ind_p_slip, part_type),      &
     4416                                                        visc)
     4417
     4418                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
     4419                                            vs,                                                    &
     4420                                            particle_pars(ind_p_size, part_type),                  &
     4421                                            particle_pars(ind_p_slip, part_type),                  &
     4422                                            nwet, temp_tmp, dens, visc,                            &
     4423                                            luu_dep,                                               &
     4424                                            r_aero_surf, ustar_surf )
     4425
     4426                bud_luu(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
    44354427
    44364428             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
     
    44384430!
    44394431!--             Sedimentation velocity
    4440                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    4441                      particle_pars(ind_p_size, part_type), &
    4442                      particle_pars(ind_p_slip, part_type), &
    4443                      visc)
    4444 
    4445                 CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    4446                      vs, &
    4447                      particle_pars(ind_p_size, part_type), &
    4448                      particle_pars(ind_p_slip, part_type), &
    4449                      nwet, temp_tmp, dens, visc, &
    4450                      luu_dep , &
    4451                      r_aero_surf, ustar_surf )
    4452 
    4453                 bud_luu(lsp) = - conc_ijk(lsp) * &
    4454                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     4432                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
     4433                                                        particle_pars(ind_p_size, part_type),      &
     4434                                                        particle_pars(ind_p_slip, part_type),      &
     4435                                                        visc)
     4436
     4437                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
     4438                                            vs,                                                    &
     4439                                            particle_pars(ind_p_size, part_type),                  &
     4440                                            particle_pars(ind_p_slip, part_type),                  &
     4441                                            nwet, temp_tmp, dens, visc,                            &
     4442                                            luu_dep ,                                              &
     4443                                            r_aero_surf, ustar_surf )
     4444
     4445                bud_luu(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
    44554446
    44564447             ELSE  !< GASES
    44574448!
    44584449!--             Read spc_name of current species for gas parameter
    4459 
    44604450                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
    44614451                   i_pspec = 0
     
    44784468!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
    44794469!--             Use density at k:
    4480 
    44814470                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
    44824471
     
    44874476!
    44884477!--             Diffusivity for DEPAC relevant gases
    4489 !--             Use default value 
     4478!--             Use default value
    44904479                diffusivity            = 0.11e-4
    44914480!
    4492 !--             overwrite with known coefficients of diffusivity from Massman (1998)
    4493                 IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
     4481!--             Overwrite with known coefficients of diffusivity from Massman (1998)
     4482                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
    44944483                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
    44954484                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
     
    45004489!
    45014490!--             Get quasi-laminar boundary layer resistance rb:
    4502                 CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland),   &
    4503                      z0h_surf, ustar_surf, diffusivity, &
    4504                      rb )
     4491                CALL get_rb_cell( ( luu_dep == ilu_water_sea )  .OR.                               &
     4492                                  ( luu_dep == ilu_water_inland ), z0h_surf, ustar_surf,          &
     4493                                    diffusivity, rb )
    45054494!
    45064495!--             Get rc_tot
    4507                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,          &
    4508                                          solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luu_dep, 2,    &
    4509                                          rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,   &
    4510                                          r_aero_surf, rb )
     4496                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
     4497                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luu_dep, &
     4498                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
     4499                                         diffusivity, r_aero_surf, rb )
    45114500!
    45124501!--             Calculate budget
     
    45194508                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
    45204509
    4521                    bud_luu(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
    4522                         (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     4510                   bud_luu(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                          &
     4511                                  ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
    45234512                ENDIF
    45244513
     
    45324521!
    45334522!--       No vegetation on bare soil, desert or ice:
    4534           IF ( ( lug_palm == ind_luv_b_soil ) .OR. &
    4535                 ( lug_palm == ind_luv_desert ) .OR. &
    4536                  ( lug_palm == ind_luv_ice ) ) THEN
     4523          IF ( ( lug_palm == ind_luv_b_soil )  .OR.  ( lug_palm == ind_luv_desert ) .OR.           &
     4524               ( lug_palm == ind_luv_ice ) ) THEN
    45374525
    45384526             lai = 0.0_wp
     
    45414529          ENDIF
    45424530
    4543          
     4531
    45444532          slinnfac = 1.0_wp
    45454533!
     
    45574545!
    45584546!--             Sedimentation velocity
    4559                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    4560                      particle_pars(ind_p_size, part_type), &
    4561                      particle_pars(ind_p_slip, part_type), &
    4562                      visc)
    4563 
    4564                 CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    4565                      vs, &
    4566                      particle_pars(ind_p_size, part_type), &
    4567                      particle_pars(ind_p_slip, part_type), &
    4568                      nwet, temp_tmp, dens, visc, &
    4569                      lug_dep,  &
    4570                      r_aero_surf, ustar_surf )
    4571 
    4572                 bud_lug(lsp) = - conc_ijk(lsp) * &
    4573                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     4547                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
     4548                                                        particle_pars(ind_p_size, part_type),      &
     4549                                                        particle_pars(ind_p_slip, part_type),      &
     4550                                                        visc)
     4551
     4552                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
     4553                                            vs,                                                    &
     4554                                            particle_pars(ind_p_size, part_type),                  &
     4555                                            particle_pars(ind_p_slip, part_type),                  &
     4556                                            nwet, temp_tmp, dens, visc,                            &
     4557                                            lug_dep,                                               &
     4558                                            r_aero_surf, ustar_surf )
     4559
     4560                bud_lug(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
    45744561
    45754562             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
     
    45774564!
    45784565!--             Sedimentation velocity
    4579                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    4580                      particle_pars(ind_p_size, part_type), &
    4581                      particle_pars(ind_p_slip, part_type), &
    4582                      visc)
    4583 
    4584                 CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    4585                      vs, &
    4586                      particle_pars(ind_p_size, part_type), &
    4587                      particle_pars(ind_p_slip, part_type), &
    4588                      nwet, temp_tmp, dens, visc, &
    4589                      lug_dep, &
    4590                      r_aero_surf, ustar_surf )
    4591 
    4592                 bud_lug(lsp) = - conc_ijk(lsp) * &
    4593                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     4566                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
     4567                                                        particle_pars(ind_p_size, part_type),      &
     4568                                                        particle_pars(ind_p_slip, part_type),      &
     4569                                                        visc)
     4570
     4571                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
     4572                                            vs,                                                    &
     4573                                            particle_pars(ind_p_size, part_type),                  &
     4574                                            particle_pars(ind_p_slip, part_type),                  &
     4575                                            nwet, temp_tmp, dens, visc,                            &
     4576                                            lug_dep,                                               &
     4577                                            r_aero_surf, ustar_surf )
     4578
     4579                bud_lug(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
    45944580
    45954581             ELSE  !< GASES
    45964582!
    45974583!--             Read spc_name of current species for gas parameter
    4598 
    45994584                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
    46004585                   i_pspec = 0
     
    46174602!--                  c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
    46184603!--             Use density at k:
    4619 
    46204604                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
    46214605!
     
    46254609!
    46264610!--             Diffusivity for DEPAC relevant gases
    4627 !--             Use default value 
     4611!--             Use default value
    46284612                diffusivity            = 0.11e-4
    46294613!
    4630 !--             overwrite with known coefficients of diffusivity from Massman (1998)
    4631                 IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
     4614!--             Overwrite with known coefficients of diffusivity from Massman (1998)
     4615                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
    46324616                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
    46334617                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
     
    46384622!
    46394623!--             Get quasi-laminar boundary layer resistance rb:
    4640                 CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland),    &
    4641                      z0h_surf, ustar_surf, diffusivity, &
    4642                      rb )
     4624                CALL get_rb_cell( ( lug_dep == ilu_water_sea )  .OR.                               &
     4625                                  ( lug_dep == ilu_water_inland ), z0h_surf, ustar_surf,          &
     4626                                    diffusivity, rb )
    46434627!
    46444628!--             Get rc_tot
    4645                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,           &
    4646                                          solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lug_dep, 2,      &
    4647                                          rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,    &
    4648                                          r_aero_surf , rb )
     4629                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
     4630                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lug_dep, &
     4631                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
     4632                                         diffusivity, r_aero_surf , rb )
    46494633!
    46504634!--             Calculate budget
     
    46554639                ELSE
    46564640
    4657                    vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
    4658 
    4659                    bud_lug(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
    4660                         (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     4641                   vd_lu = 1.0_wp / ( r_aero_surf + rb + rc_tot )
     4642
     4643                   bud_lug(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                          &
     4644                                  ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
    46614645                ENDIF
    46624646
     
    46824666             rs = 0.0_wp
    46834667             rb = 0.0_wp
    4684              rc_tot = 0.0_wp 
     4668             rc_tot = 0.0_wp
    46854669             IF ( spc_names(lsp) == 'PM10' )  THEN
    46864670                part_type = 1
    46874671!
    46884672!--             Sedimentation velocity
    4689                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    4690                      particle_pars(ind_p_size, part_type), &
    4691                      particle_pars(ind_p_slip, part_type), &
    4692                      visc)
    4693 
    4694                 CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
    4695                      particle_pars(ind_p_size, part_type), &
    4696                      particle_pars(ind_p_slip, part_type), &
    4697                      nwet, temp_tmp, dens, visc,              &
    4698                      lud_dep, r_aero_surf, ustar_surf )
    4699 
    4700                 bud_lud(lsp) = - conc_ijk(lsp) * &
    4701                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     4673                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
     4674                                                        particle_pars(ind_p_size, part_type),      &
     4675                                                        particle_pars(ind_p_slip, part_type),      &
     4676                                                        visc)
     4677
     4678                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs,                                         &
     4679                                            particle_pars(ind_p_size, part_type),                  &
     4680                                            particle_pars(ind_p_slip, part_type),                  &
     4681                                            nwet, temp_tmp, dens, visc,                            &
     4682                                            lud_dep, r_aero_surf, ustar_surf )
     4683
     4684                bud_lud(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
    47024685
    47034686             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
     
    47054688!
    47064689!--             Sedimentation velocity
    4707                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    4708                      particle_pars(ind_p_size, part_type), &
    4709                      particle_pars(ind_p_slip, part_type), &
    4710                      visc)
    4711 
    4712                 CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
    4713                      particle_pars(ind_p_size, part_type), &
    4714                      particle_pars(ind_p_slip, part_type), &
    4715                      nwet, temp_tmp, dens, visc, &
    4716                      lud_dep, &
    4717                      r_aero_surf, ustar_surf )
    4718 
    4719                 bud_lud(lsp) = - conc_ijk(lsp) * &
    4720                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     4690                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
     4691                                                        particle_pars(ind_p_size, part_type),      &
     4692                                                        particle_pars(ind_p_slip, part_type),      &
     4693                                                        visc)
     4694
     4695                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs,                                         &
     4696                                            particle_pars(ind_p_size, part_type),                  &
     4697                                            particle_pars(ind_p_slip, part_type),                  &
     4698                                            nwet, temp_tmp, dens, visc,                            &
     4699                                            lud_dep,                                               &
     4700                                            r_aero_surf, ustar_surf )
     4701
     4702                bud_lud(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
    47214703
    47224704             ELSE  !< GASES
    47234705!
    47244706!--             Read spc_name of current species for gas PARAMETER
    4725 
    47264707                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
    47274708                   i_pspec = 0
     
    47524733!
    47534734!--             Diffusivity for DEPAC relevant gases
    4754 !--             Use default value 
     4735!--             Use default value
    47554736                diffusivity = 0.11e-4
    47564737!
    4757 !--             overwrite with known coefficients of diffusivity from Massman (1998)
    4758                 IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
     4738!--             Overwrite with known coefficients of diffusivity from Massman (1998)
     4739                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
    47594740                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
    47604741                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
     
    47654746!
    47664747!--             Get quasi-laminar boundary layer resistance rb:
    4767                 CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland),   &
    4768                      z0h_surf, ustar_surf, diffusivity, rb )
     4748                CALL get_rb_cell( ( lud_dep == ilu_water_sea )  .OR.                               &
     4749                                  ( lud_dep == ilu_water_inland ), z0h_surf, ustar_surf,           &
     4750                                    diffusivity, rb )
    47694751!
    47704752!--             Get rc_tot
    4771                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         &
    4772                                          solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lud_dep, 2,    &
    4773                                          rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
    4774                                          r_aero_surf , rb )
    4775 !
    4776 !--             Calculate budget 
     4753                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
     4754                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lud_dep, &
     4755                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
     4756                                         diffusivity, r_aero_surf , rb )
     4757!
     4758!--             Calculate budget
    47774759                IF ( rc_tot <= 0.0 )  THEN
    47784760
     
    47834765                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
    47844766
    4785                    bud_lud(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
    4786                         (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     4767                   bud_lud(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                          &
     4768                                  ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
    47874769                ENDIF
    47884770
     
    47984780
    47994781
    4800           bud(lsp) = surf_usm_h%frac(m,ind_veg_wall) * bud_luu(lsp) + &
    4801                surf_usm_h%frac(m,ind_pav_green) * bud_lug(lsp) + &
    4802                surf_usm_h%frac(m,ind_wat_win) * bud_lud(lsp)
     4782          bud(lsp) = surf_usm_h%frac(m,ind_veg_wall)  * bud_luu(lsp) +                            &
     4783                     surf_usm_h%frac(m,ind_pav_green) * bud_lug(lsp) +                            &
     4784                     surf_usm_h%frac(m,ind_wat_win)  * bud_lud(lsp)
    48034785!
    48044786!--       Compute new concentration
     
    48154797
    48164798
    4817 !------------------------------------------------------------------------------!
     4799!--------------------------------------------------------------------------------------------------!
    48184800! Description:
    48194801! ------------
     
    48214803!>
    48224804!> DEPAC:
    4823 !> Code of the DEPAC routine and corresponding subroutines below from the DEPAC
    4824 !> module of the LOTOS-EUROS model (Manders et al., 2017)
     4805!> Code of the DEPAC routine and corresponding subroutines below from the DEPAC module of the
     4806!> LOTOS-EUROS model (Manders et al., 2017)
    48254807!>
    4826 !> Original DEPAC routines by RIVM and TNO (2015), for Documentation see
    4827 !> van Zanten et al., 2010.
    4828 !------------------------------------------------------------------------------!
    4829  SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, solar_rad, sinphi,    &
    4830       rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, conc_ijk_ugm3, diffusivity,  &
    4831       ra, rb ) 
     4808!> Original DEPAC routines by RIVM and TNO (2015), for Documentation see van Zanten et al., 2010.
     4809!--------------------------------------------------------------------------------------------------!
     4810 SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, solar_rad, sinphi, rh, lai, sai,&
     4811                                nwet, lu, iratns, rc_tot, ccomp_tot, p, conc_ijk_ugm3, diffusivity,&
     4812                                ra, rb )
    48324813!
    48334814!--   Some of depac arguments are OPTIONAL:
    48344815!--    A. compute Rc_tot without compensation points (ccomp_tot will be zero):
    4835 !--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi])
     4816!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot,&
     4817!--                    ccomp_tot, [smi])
    48364818!--    B. compute Rc_tot with compensation points (used for LOTOS-EUROS):
    4837 !--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
    4838 !--                c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water)
     4819!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot,&
     4820!--                    ccomp_tot, [smi], c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water)
    48394821!--
    48404822!--    C. compute effective Rc based on compensation points (used for OPS):
    4841 !--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
    4842 !--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
    4843 !--               ra, rb, rc_eff)
    4844 !--    X1. Extra (OPTIONAL) output variables:
    4845 !--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
    4846 !--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
    4847 !--               ra, rb, rc_eff, &
    4848 !--               gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out)
    4849 !--    X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves):
    4850 !--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
    4851 !--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
    4852 !--               ra, rb, rc_eff, &
    4853 !--               gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, &
    4854 !--               calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA)
     4823!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot,&
     4824!--                    ccomp_tot, [smi], c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, ra,  &
     4825!--                    rb, rc_eff)
     4826!--    X1. Extra (OPTIONAL) output variables:
     4827!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot,&
     4828!--                    ccomp_tot, [smi], c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, ra,  &
     4829!--                    rb, rc_eff, gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out,    &
     4830!--                    lai_out, sai_out)
     4831!--    X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves):
     4832!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot,&
     4833!--                    ccomp_tot, [smi], c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, ra,  &
     4834!--                    rb, rc_eff, gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out,    &
     4835!--                    lai_out, sai_out, calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA)
    48554836
    48564837
     
    48584839                                                     !< 'HNO3','NO','NO2','O3','SO2','NH3'
    48594840    INTEGER(iwp), INTENT(IN) ::  day_of_year         !< day of year, 1 ... 365 (366)
    4860     INTEGER(iwp), INTENT(IN) ::  nwet                !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
    4861     INTEGER(iwp), INTENT(IN) ::  lu                  !< land use type, lu = 1,...,nlu
    48624841    INTEGER(iwp), INTENT(IN) ::  iratns              !< index for NH3/SO2 ratio used for SO2:
    48634842                                                     !< iratns = 1: low NH3/SO2
    48644843                                                     !< iratns = 2: high NH3/SO2
    48654844                                                     !< iratns = 3: very low NH3/SO2
    4866     REAL(wp), INTENT(IN) ::  lat                     !< latitude Northern hemisphere (degrees) (S. hemisphere not possible)
     4845    INTEGER(iwp), INTENT(IN) ::  lu                  !< land use type, lu = 1,...,nlu
     4846    INTEGER(iwp), INTENT(IN) ::  nwet                !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     4847
     4848    REAL(wp), INTENT(IN) ::  conc_ijk_ugm3           !< actual atmospheric concentration (ug/m3), in
     4849                                                     !< DEPAC=Catm
     4850    REAL(wp), INTENT(IN) ::  diffusivity             !< diffusivity
     4851    REAL(wp), INTENT(IN) ::  lai                     !< one-sidedleaf area index (-)
     4852    REAL(wp), INTENT(IN) ::  lat                     !< latitude Northern hemisphere (degrees) (S.
     4853                                                     !< hemisphere not possible)
     4854    REAL(wp), INTENT(IN) ::  p                       !< pressure (Pa)
     4855    REAL(wp), INTENT(IN) ::  ra                      !< aerodynamic resistance (s/m)
     4856    REAL(wp), INTENT(IN) ::  rb                      !< boundary layer resistance (s/m)
     4857    REAL(wp), INTENT(IN) ::  rh                      !< relative humidity (%)
     4858    REAL(wp), INTENT(IN) ::  sai                     !< surface area index (-) (lai + branches and
     4859                                                     !< stems)
     4860    REAL(wp), INTENT(IN) ::  sinphi                  !< sin of solar elevation angle
     4861    REAL(wp), INTENT(IN) ::  solar_rad               !< solar radiation, dirict+diffuse (W/m2)
    48674862    REAL(wp), INTENT(IN) ::  t                       !< temperature (C)
    48684863    REAL(wp), INTENT(IN) ::  ust                     !< friction velocity (m/s)
    4869     REAL(wp), INTENT(IN) ::  solar_rad               !< solar radiation, dirict+diffuse (W/m2)
    4870     REAL(wp), INTENT(IN) ::  sinphi                  !< sin of solar elevation angle
    4871     REAL(wp), INTENT(IN) ::  rh                      !< relative humidity (%)
    4872     REAL(wp), INTENT(IN) ::  lai                     !< one-sidedleaf area index (-)
    4873     REAL(wp), INTENT(IN) ::  sai                     !< surface area index (-) (lai + branches and stems)
    4874     REAL(wp), INTENT(IN) ::  diffusivity             !< diffusivity
    4875     REAL(wp), INTENT(IN) ::  p                       !< pressure (Pa)
    4876     REAL(wp), INTENT(IN) ::  conc_ijk_ugm3           !< actual atmospheric concentration (ug/m3), in DEPAC=Catm
    4877     REAL(wp), INTENT(IN) ::  ra                      !< aerodynamic resistance (s/m)
    4878     REAL(wp), INTENT(IN) ::  rb                      !< boundary layer resistance (s/m)
    4879 
     4864
     4865    REAL(wp), INTENT(OUT) ::  ccomp_tot              !< total compensation point (ug/m3)
     4866!                                                    !< [= 0 for species that don't have a compensation
    48804867    REAL(wp), INTENT(OUT) ::  rc_tot                 !< total canopy resistance Rc (s/m)
    4881     REAL(wp), INTENT(OUT) ::  ccomp_tot              !< total compensation point (ug/m3)
    4882 !                                                     !< [= 0 for species that don't have a compensation
     4868
    48834869!-- Local variables:
    48844870!
     
    48994885
    49004886    LOGICAL ::  ready                                !< Rc has been set:
    4901     !< = 1 -> constant Rc
    4902     !< = 2 -> temperature dependent Rc
     4887                                                     !< = 1 -> constant Rc
     4888                                                     !< = 2 -> temperature dependent Rc
    49034889!
    49044890!-- Vegetation indicators:
    49054891    LOGICAL ::  lai_present                          !< leaves are present for current land use type
    4906     LOGICAL ::  sai_present                          !< vegetation is present for current land use type
    4907 
     4892    LOGICAL ::  sai_present                          !< vegetation is present for current land use
     4893                                                     !< type
     4894
     4895    REAL(wp) ::  csoil                               !< soil compensation point (ug/m3)
     4896    REAL(wp) ::  cstom                               !< stomatal compensation point (ug/m3)
     4897    REAL(wp) ::  cw                                  !< external leaf surface compensation point
     4898                                                     !< (ug/m3)
     4899    REAL(wp) ::  gc_tot                              !< total canopy conductance (m/s)
     4900    REAL(wp) ::  gsoil_eff                           !< effective soil conductance (m/s)
     4901    REAL(wp) ::  gstom                               !< stomatal conductance (m/s)
     4902    REAL(wp) ::  gw                                  !< external leaf conductance (m/s)
    49084903!    REAL(wp) ::  laimax                              !< maximum leaf area index (-)
    4909     REAL(wp) ::  gw                                  !< external leaf conductance (m/s)
    4910     REAL(wp) ::  gstom                               !< stomatal conductance (m/s)
    4911     REAL(wp) ::  gsoil_eff                           !< effective soil conductance (m/s)
    4912     REAL(wp) ::  gc_tot                              !< total canopy conductance (m/s)
    4913     REAL(wp) ::  cw                                  !< external leaf surface compensation point (ug/m3)
    4914     REAL(wp) ::  cstom                               !< stomatal compensation point (ug/m3)
    4915     REAL(wp) ::  csoil                               !< soil compensation point (ug/m3)
    49164904!
    49174905!-- Next statement is just to avoid compiler warning about unused variable
     
    49314919
    49324920    CASE ( 'NO', 'no' )
    4933        icmp = icmp_no 
     4921       icmp = icmp_no
    49344922
    49354923    CASE ( 'NH3', 'nh3' )
     
    49794967!
    49804968!--    External conductance:
    4981        CALL rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai,gw )         
     4969       CALL rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai,gw )
    49824970!
    49834971!--    Stomatal conductance:
    4984        CALL rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity, gstom, p )
     4972       CALL rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity,  &
     4973                      gstom, p )
    49854974!
    49864975!--    Effective soil conductance:
     
    50004989!
    50014990!--    Effective Rc based on compensation points:
    5002 !--        IF ( present(rc_eff) ) then
     4991!--        IF ( PRESENT( rc_eff ) )  THEN
    50034992!--          check on required arguments:
    5004 !--           IF ( (.not. present(conc_ijk_ugm3)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then
    5005 !--              stop 'output argument rc_eff requires input arguments conc_ijk_ugm3, ra and rb'
    5006 !--           END IF
    5007 !
    5008 !--       Compute rc_eff :
    5009        !&nb