Changeset 4559 for palm/trunk/SOURCE/chemistry_model_mod.f90

Timestamp:

Jun 11, 2020 8:51:48 AM (4 years ago)

Author:

raasch

Message:

files re-formatted to follow the PALM coding standard

File:

• palm/trunk/SOURCE/chemistry_model_mod.f90 (modified) (267 diffs)

palm/trunk/SOURCE/chemistry_model_mod.f90

@@ -1 +1 @@
 !> @file chemistry_model_mod.f90
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! This file is part of the PALM model system.
 !
-! PALM is free software: you can redistribute it and/or modify it under the
-! terms of the GNU General Public License as published by the Free Software
-! Foundation, either version 3 of the License, or (at your option) any later
-! version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
+! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
+! Public License as published by the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
+! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+! Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with PALM. If not, see
+! <http://www.gnu.org/licenses/>.
 !
 ! Copyright 2017-2020 Leibniz Universitaet Hannover
 ! Copyright 2017-2020 Karlsruhe Institute of Technology
 ! Copyright 2017-2020 Freie Universitaet Berlin
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 !
 ! Current revisions:
@@ -27 +26 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+! 4550 2020-05-29 15:22:13Z raasch
 ! bugfix for reading local restart data
-! 
+!
 ! 4544 2020-05-21 14:43:05Z raasch
 ! conc_av changed from pointer to allocatable array, array spec_conc_av removed
-! 
+!
 ! 4542 2020-05-19 15:45:12Z raasch
 ! redundant if statement removed
-! 
+!
 ! 4535 2020-05-15 12:07:23Z raasch
 ! bugfix for restart data format query
-! 
+!
 ! 4517 2020-05-03 14:29:30Z raasch
 ! added restart with MPI-IO
-! 
+!
 ! 4511 2020-04-30 12:20:40Z raasch
 ! decycling replaced by explicit setting of lateral boundary conditions
@@ -50 +52 @@
 ! use statement for exchange horiz added,
 ! bugfix for call of exchange horiz 2d
-! 
+!
 ! 4442 2020-03-04 19:21:13Z suehring
-! Change order of dimension in surface array %frac to allow for better 
+! Change order of dimension in surface array %frac to allow for better
 ! vectorization.
 !
 ! 4441 2020-03-04 19:20:35Z suehring
 ! in subroutine chem_init (ECC)
-! - allows different init paths emission data for legacy
-!   mode emission and on-demand mode
-! in subroutine chem_init_internal (ECC)
-! - reads netcdf file only when legacy mode is activated
-!   (i.e., emiss_read_legacy_mode = .TRUE.)
-!   otherwise file is read once at the beginning to obtain
-!   header information, and emission data are extracted on
-!   an on-demand basis
+! - allows different init paths emission data for legacy mode emission and on-demand mode in
+!   subroutine chem_init_internal (ECC)
+! - reads netcdf file only when legacy mode is activated (i.e., emiss_read_legacy_mode = .TRUE.)
+!   otherwise file is read once at the beginning to obtain header information, and emission data
+!   are extracted on an on-demand basis
 !
 ! 4372 2020-01-14 10:20:35Z banzhafs
-! chem_parin : added handler for new namelist item emiss_legacy_read_mode (ECC)
-! added messages
+! chem_parin : added handler for new namelist item emiss_legacy_read_mode (ECC) added messages
 ! CM0465 - legacy read mode selection
 ! CM0466 - legacy read mode force selection
 ! CM0467 - new read mode selection
 !
-! 4370 2020-01-10 14:00:44Z raasch 
+! 4370 2020-01-10 14:00:44Z raasch
 ! vector directives added to force vectorization on Intel19 compiler
 !
 ! 4346 2019-12-18 11:55:56Z motisi
-! Introduction of wall_flags_total_0, which currently sets bits based on static
-! topography information used in wall_flags_static_0
-! 
+! Introduction of wall_flags_total_0, which currently sets bits based on static topography
+! information used in wall_flags_static_0
+!
 ! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
-! 
+!
 ! 4306 2019-11-25 12:04:48Z banzhafs
 ! Corretion for r4304 commit
@@ -96 +94 @@
 !
 ! 4273 2019-10-24 13:40:54Z monakurppa
-! Add logical switches nesting_chem and nesting_offline_chem (both .TRUE.
-! by default)
-! 
+! Add logical switches nesting_chem and nesting_offline_chem (both .TRUE. by default)
+!
 ! 4272 2019-10-23 15:18:57Z schwenkel
-! Further modularization of boundary conditions: moved boundary conditions to
-! respective modules
+! Further modularization of boundary conditions: moved boundary conditions to respective modules
 !
 ! 4268 2019-10-17 11:29:38Z schwenkel
 ! Moving module specific boundary conditions from time_integration to module
-! 
+!
 ! 4230 2019-09-11 13:58:14Z suehring
-! Bugfix, initialize mean profiles also in restart runs. Also initialize 
-! array used for Runge-Kutta tendecies in restart runs.  
-! 
+! Bugfix, initialize mean profiles also in restart runs. Also initialize array used for Runge-Kutta
+! tendecies in restart runs.
+!
 ! 4227 2019-09-10 18:04:34Z gronemeier
 ! implement new palm_date_time_mod
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 4167 2019-08-16 11:01:48Z suehring
-! Changed behaviour of masked output over surface to follow terrain and ignore 
-! buildings (J.Resler, T.Gronemeier)
-! 
-! 
+! Changed behaviour of masked output over surface to follow terrain and ignore buildings
+! (J.Resler, T.Gronemeier)
+!
 ! 4166 2019-08-16 07:54:21Z resler
 ! Bugfix in decycling
-! 
+!
 ! 4115 2019-07-24 12:50:49Z suehring
 ! Fix faulty IF statement in decycling initialization
-! 
+!
 ! 4110 2019-07-22 17:05:21Z suehring
-! - Decycling boundary conditions are only set at the ghost points not on the
-!   prognostic grid points 
-! - Allocation and initialization of special advection flags cs_advc_flags_s
-!   used for chemistry. These are exclusively used for chemical species in 
-!   order to distinguish from the usually-used flags which might be different
-!   when decycling is applied in combination with cyclic boundary conditions. 
-!   Moreover, cs_advc_flags_s considers extended zones around buildings where 
-!   first-order upwind scheme is applied for the horizontal advection terms, 
-!   in order to overcome high concentration peaks due to stationary numerical 
-!   oscillations caused by horizontal advection discretization.
-! 
+! - Decycling boundary conditions are only set at the ghost points not on the prognostic grid points
+! - Allocation and initialization of special advection flags cs_advc_flags_s used for chemistry.
+!   These are exclusively used for chemical species in order to distinguish from the usually-used
+!   flags which might be different when decycling is applied in combination with cyclic boundary
+!   conditions.
+!   Moreover, cs_advc_flags_s considers extended zones around buildings where first-order upwind
+!   scheme is applied for the horizontal advection terms, in order to overcome high concentration
+!   peaks due to stationary numerical oscillations caused by horizontal advection discretization.
+!
 ! 4109 2019-07-22 17:00:34Z suehring
-! Slightly revise setting of boundary conditions at horizontal walls, use 
-! data-structure offset index instead of pre-calculate it for each facing
-! 
+! Slightly revise setting of boundary conditions at horizontal walls, use data-structure offset
+! index instead of pre-calculate it for each facing
+!
 ! 4080 2019-07-09 18:17:37Z suehring
 ! Restore accidantly removed limitation to positive values
-! 
+!
 ! 4079 2019-07-09 18:04:41Z suehring
-! Application of monotonic flux limiter for the vertical scalar advection 
-! up to the topography top (only for the cache-optimized version at the 
-! moment).
-! 
+! Application of monotonic flux limiter for the vertical scalar advection up to the topography top
+! (only for the cache-optimized version at the moment).
+!
 ! 4069 2019-07-01 14:05:51Z Giersch
-! Masked output running index mid has been introduced as a local variable to 
-! avoid runtime error (Loop variable has been modified) in time_integration
-! 
+! Masked output running index mid has been introduced as a local variable to avoid runtime error
+! (Loop variable has been modified) in time_integration
+!
 ! 4029 2019-06-14 14:04:35Z raasch
 ! nest_chemistry option removed
-! 
+!
 ! 4004 2019-05-24 11:32:38Z suehring
-! in subroutine chem_parin check emiss_lod / mod_emis only
-! when emissions_anthropogenic is activated in namelist (E.C. Chan)
-! 
+! in subroutine chem_parin check emiss_lod / mod_emis only when emissions_anthropogenic is activated
+! in namelist (E.C. Chan)
+!
 ! 3968 2019-05-13 11:04:01Z suehring
-! - added "emiss_lod" which serves the same function as "mode_emis"
-!   both will be synchronized with emiss_lod having pirority over
-!   mode_emis (see informational messages)
+! - added "emiss_lod" which serves the same function as "mode_emis" both will be synchronized with
+!   emiss_lod having pirority over mode_emis (see informational messages)
 ! - modified existing error message and introduced new informational messages
 !    - CM0436 - now also applies to invalid LOD settings
 !    - CM0463 - emiss_lod take precedence in case of conflict with mod_emis
 !    - CM0464 - emiss_lod valued assigned based on mode_emis if undefined
-! 
+!
 ! 3930 2019-04-24 14:57:18Z forkel
 ! Changed chem_non_transport_physics to chem_non_advective_processes
-! 
-! 
+!
 ! 3929 2019-04-24 12:52:08Z banzhafs
 ! Correct/complete module_interface introduction for chemistry model
@@ -182 +172 @@
 ! 3784 2019-03-05 14:16:20Z banzhafs
 ! 2D output of emission fluxes
-! 
+!
 ! 3784 2019-03-05 14:16:20Z banzhafs
-! Bugfix, uncomment erroneous commented variable used for dry deposition. 
-! Bugfix in 3D emission output. 
-! 
+! Bugfix, uncomment erroneous commented variable used for dry deposition.
+! Bugfix in 3D emission output.
+!
 ! 3784 2019-03-05 14:16:20Z banzhafs
-! Changes related to global restructuring of location messages and introduction 
+! Changes related to global restructuring of location messages and introduction
 ! of additional debug messages
-! 
+!
 ! 3784 2019-03-05 14:16:20Z banzhafs
 ! some formatting of the deposition code
@@ -196 +186 @@
 ! 3784 2019-03-05 14:16:20Z banzhafs
 ! some formatting
-! 
+!
 ! 3784 2019-03-05 14:16:20Z banzhafs
-! added cs_mech to USE chem_gasphase_mod 
-! 
+! added cs_mech to USE chem_gasphase_mod
+!
 ! 3784 2019-03-05 14:16:20Z banzhafs
 ! renamed get_mechanismname to get_mechanism_name
 ! renamed do_emiss to emissions_anthropogenic and do_depo to deposition_dry (ecc)
-! 
+!
 ! 3784 2019-03-05 14:16:20Z banzhafs
 ! Unused variables removed/taken care of
-! 
 !
 ! 3784 2019-03-05 14:16:20Z forkel
 ! Replaced READ from unit 10 by CALL get_mechanismname also in chem_header
 !
-!
 ! 3783 2019-03-05 13:23:50Z forkel
-! Removed forgotte write statements an some unused variables (did not touch the 
-! parts related to deposition)
-! 
-! 
+! Removed forgotte write statements an some unused variables (did not touch the parts related to
+! deposition)
+!
 ! 3780 2019-03-05 11:19:45Z forkel
 ! Removed READ from unit 10, added CALL get_mechanismname
-! 
-! 
+!
 ! 3767 2019-02-27 08:18:02Z raasch
 ! unused variable for file index removed from rrd-subroutines parameter list
-! 
+!
 ! 3738 2019-02-12 17:00:45Z suehring
 ! Clean-up debug prints
-! 
+!
 ! 3737 2019-02-12 16:57:06Z suehring
-! Enable mesoscale offline nesting for chemistry variables as well as 
-! initialization of chemistry via dynamic input file.
-! 
+! Enable mesoscale offline nesting for chemistry variables as well as initialization of chemistry
+! via dynamic input file.
+!
 ! 3719 2019-02-06 13:10:18Z kanani
-! Resolved cpu logpoint overlap with all progn.equations, moved cpu_log call
-! to prognostic_equations for better overview
-! 
+! Resolved cpu logpoint overlap with all progn.equations, moved cpu_log call to prognostic_equations
+! for better overview
+!
 ! 3700 2019-01-26 17:03:42Z knoop
 ! Some interface calls moved to module_interface + cleanup
-! 
+!
 ! 3664 2019-01-09 14:00:37Z forkel
 ! Replaced misplaced location message by @todo
-! 
-! 
+!
 ! 3654 2019-01-07 16:31:57Z suehring
 ! Disable misplaced location message
-! 
+!
 ! 3652 2019-01-07 15:29:59Z forkel
 ! Checks added for chemistry mechanism, parameter chem_mechanism added (basit)
-! 
+!
 ! 2718 2018-01-02 08:49:38Z maronga
 ! Initial revision
 !
-! 
+!
 !
 !
@@ -263 +248 @@
 !
 !
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
@@ -269 +254 @@
 !> @todo Extend chem_species type by nspec and nvar as addititional elements (RF)
 !> @todo Check possibility to reduce dimension of chem_species%conc from nspec to nvar (RF)
-!> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not 
-!>       allowed to use the chemistry model in a precursor run and additionally 
-!>       not using it in a main run
+!> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not allowed to use the
+!>       chemistry model in a precursor run and additionally not using it in a main run
 !> @todo Implement turbulent inflow of chem spcs in inflow_turbulence.
 !> @todo Separate boundary conditions for each chem spcs to be implemented
-!> @todo Currently only total concentration are calculated. Resolved, parameterized
-!>       and chemistry fluxes although partially and some completely coded but
-!>       are not operational/activated in this version. bK.
-!> @todo slight differences in passive scalar and chem spcs when chem reactions
-!>       turned off. Need to be fixed. bK
+!> @todo Currently only total concentration are calculated. Resolved, parameterized and chemistry
+!>       fluxes although partially and some completely coded but are not operational/activated in
+!>       this version. bK.
+!> @todo slight differences in passive scalar and chem spcs when chem reactions turned off. Need to
+!>       be fixed. bK
 !> @todo test nesting for chem spcs, was implemented by suehring (kanani)
 !> @todo chemistry error messages
-! 
-!------------------------------------------------------------------------------!
+!
+!--------------------------------------------------------------------------------------------------!
 
  MODULE chemistry_model_mod
@@ -313 +297 @@
 
     USE control_parameters,                                                                        &
-         ONLY:  bc_lr_cyc, bc_ns_cyc,                                                              &
+         ONLY:  air_chemistry,                                                                     &
                 bc_dirichlet_l,                                                                    &
                 bc_dirichlet_n,                                                                    &
@@ -319 +303 @@
                 bc_dirichlet_s,                                                                    &
                 bc_radiation_l,                                                                    &
+                bc_lr_cyc,                                                                         &
+                bc_ns_cyc,                                                                         &
                 bc_radiation_n,                                                                    &
                 bc_radiation_r,                                                                    &
                 bc_radiation_s,                                                                    &
                 debug_output,                                                                      &
-                dt_3d, humidity, initializing_actions, message_string,                             &
-                omega, tsc, intermediate_timestep_count, intermediate_timestep_count_max,          &
+                dt_3d,                                                                             &
+                humidity,                                                                          &
+                initializing_actions,                                                              &
+                intermediate_timestep_count,                                                       &
+                intermediate_timestep_count_max,                                                   &
                 max_pr_user,                                                                       &
+                message_string,                                                                    &
                 monotonic_limiter_z,                                                               &
+                omega,                                                                             &
                 restart_data_format_output,                                                        &
                 scalar_advec,                                                                      &
-                timestep_scheme, use_prescribed_profile_data, ws_scheme_sca, air_chemistry
+                timestep_scheme,                                                                   &
+                tsc,                                                                               &
+                use_prescribed_profile_data,                                                       &
+                ws_scheme_sca
 
     USE arrays_3d,                                                                                 &
@@ -359 +353 @@
     SAVE
 
+    INTEGER, DIMENSION(nkppctrl)                           ::  icntrl       !< 20 integer parameters for fine tuning KPP code
+
+    REAL(kind=wp), PUBLIC ::  cs_time_step = 0._wp
+
+    REAL(kind=wp), DIMENSION(nkppctrl)                     ::  rcntrl       !< 20 real parameters for fine tuning of KPP code
+                                                                            !< (e.g starting internal timestep of solver)
+
     REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_1  !< pointer for swapping of timelevels for conc
     REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_2  !< pointer for swapping of timelevels for conc
     REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_3  !< pointer for swapping of timelevels for conc
-    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  freq_1       !< pointer for phtolysis frequncies 
+    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  freq_1       !< pointer for phtolysis frequncies
                                                                             !< (only 1 timelevel required)
-    INTEGER, DIMENSION(nkppctrl)                           ::  icntrl       !< 20 integer parameters for fine tuning KPP code
                                                                             !< (e.g. solver type)
-    REAL(kind=wp), DIMENSION(nkppctrl)                     ::  rcntrl       !< 20 real parameters for fine tuning of KPP code
-                                                                            !< (e.g starting internal timestep of solver)
-
-    REAL(kind=wp), PUBLIC ::  cs_time_step = 0._wp
+
 
 !
 !-- Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010)
-    !
-    INTEGER(iwp), PARAMETER ::  nlu_dep = 15                   !< Number of DEPAC landuse classes (lu's)
-    INTEGER(iwp), PARAMETER ::  ncmp = 10                      !< Number of DEPAC gas components
-    INTEGER(iwp), PARAMETER ::  nposp = 69                     !< Number of possible species for deposition
-!
-!-- DEPAC landuse classes as defined in LOTOS-EUROS model v2.1                              
-    INTEGER(iwp) ::  ilu_grass              = 1        
-    INTEGER(iwp) ::  ilu_arable             = 2        
-    INTEGER(iwp) ::  ilu_permanent_crops    = 3         
-    INTEGER(iwp) ::  ilu_coniferous_forest  = 4         
-    INTEGER(iwp) ::  ilu_deciduous_forest   = 5         
-    INTEGER(iwp) ::  ilu_water_sea          = 6        
-    INTEGER(iwp) ::  ilu_urban              = 7        
-    INTEGER(iwp) ::  ilu_other              = 8  
-    INTEGER(iwp) ::  ilu_desert             = 9  
-    INTEGER(iwp) ::  ilu_ice                = 10 
-    INTEGER(iwp) ::  ilu_savanna            = 11 
-    INTEGER(iwp) ::  ilu_tropical_forest    = 12 
-    INTEGER(iwp) ::  ilu_water_inland       = 13 
-    INTEGER(iwp) ::  ilu_mediterrean_scrub  = 14 
-    INTEGER(iwp) ::  ilu_semi_natural_veg   = 15 
+    INTEGER(iwp), PARAMETER ::  nlu_dep = 15            !< Number of DEPAC landuse classes (lu's)
+    INTEGER(iwp), PARAMETER ::  ncmp = 10               !< Number of DEPAC gas components
+    INTEGER(iwp), PARAMETER ::  nposp = 69              !< Number of possible species for deposition
+!
+!-- DEPAC landuse classes as defined in LOTOS-EUROS model v2.1
+    INTEGER(iwp) ::  ilu_grass              = 1
+    INTEGER(iwp) ::  ilu_arable             = 2
+    INTEGER(iwp) ::  ilu_permanent_crops    = 3
+    INTEGER(iwp) ::  ilu_coniferous_forest  = 4
+    INTEGER(iwp) ::  ilu_deciduous_forest   = 5
+    INTEGER(iwp) ::  ilu_water_sea          = 6
+    INTEGER(iwp) ::  ilu_urban              = 7
+    INTEGER(iwp) ::  ilu_other              = 8
+    INTEGER(iwp) ::  ilu_desert             = 9
+    INTEGER(iwp) ::  ilu_ice                = 10
+    INTEGER(iwp) ::  ilu_savanna            = 11
+    INTEGER(iwp) ::  ilu_tropical_forest    = 12
+    INTEGER(iwp) ::  ilu_water_inland       = 13
+    INTEGER(iwp) ::  ilu_mediterrean_scrub  = 14
+    INTEGER(iwp) ::  ilu_semi_natural_veg   = 15
 
 !
@@ -405 +401 @@
 !
 !-- Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) to W m-2
-    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  alpha   =(/ 0.009, 0.009, 0.009, 0.006, 0.006, -999., -999., 0.009, -999.,      &
-         -999., 0.009, 0.006, -999., 0.009, 0.008/)*4.57
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  alpha   = (/ 0.009, 0.009, 0.009, 0.006, 0.006,    &
+                    -999.0, -999., 0.009, -999.0, -999.0, 0.009, 0.006, -999.0, 0.009, 0.008 /)*4.57
 !
 !-- Set temperatures per land use for f_temp
-    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmin = (/ 12.0, 12.0,  12.0,  0.0,  0.0, -999., -999., 12.0, -999., -999.,      &
-         12.0,  0.0, -999., 12.0,  8.0/)
-    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  topt = (/ 26.0, 26.0,  26.0, 18.0, 20.0, -999., -999., 26.0, -999., -999.,      &
-         26.0, 20.0, -999., 26.0, 24.0 /)
-    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmax = (/ 40.0, 40.0,  40.0, 36.0, 35.0, -999., -999., 40.0, -999., -999.,      &
-         40.0, 35.0, -999., 40.0, 39.0 /)
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmin = (/ 12.0,  12.0,  12.0,  0.0,  0.0, -999.0,  &
+                                      -999.0, 12.0, -999.0, -999.0, 12.0,  0.0, -999.0, 12.0,  8.0/)
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  topt = (/ 26.0, 26.0,  26.0, 18.0, 20.0, -999.0,   &
+                                     -999.0, 26.0, -999.0, -999.0, 26.0, 20.0, -999.0, 26.0, 24.0 /)
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmax = (/ 40.0, 40.0,  40.0, 36.0, 35.0, -999.0,   &
+                                    -999.0, 40.0, -999.0, -999.0,  40.0, 35.0, -999.0, 40.0, 39.0 /)
 !
 !-- Set f_min:
-    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,  &
-         0.1, -999., 0.01, 0.04/)
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  f_min = (/ 0.01, 0.01, 0.01, 0.1, 0.1, -999.0,     &
+                                       -999.0, 0.01, -999.0, -999.0, 0.01, 0.1, -999.0, 0.01, 0.04/)
 
 !
 !-- Set maximal conductance (m/s)
 !-- (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from  mmol O3/m2/s to m/s
-    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  g_max = (/ 270., 300., 300., 140., 150., -999., -999., 270., -999., -999.,      &
-         270., 150., -999., 300., 422./)/41000.
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  g_max = (/ 270.0, 300.0, 300.0, 140.0, 150.0,      &
+              -999.0, -999.0, 270.0, -999.0, -999.0, 270.0, 150.0, -999.0, 300.0, 422.0 /) / 41000.0
 !
 !-- Set max, min for vapour pressure deficit vpd
-    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,      &
-         1.0, -999., 0.9, 2.8/) 
-    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,     &
-         3.25, -999., 2.8, 4.5/) 
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_max = (/ 1.3, 0.9, 0.9, 0.5, 1.0, -999.0,      &
+                                           -999.0, 1.3, -999.0, -999.0, 1.3, 1.0, -999.0, 0.9, 2.8/)
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_min = (/ 3.0, 2.8, 2.8, 3.0, 3.25, -999.0,     &
+                                         -999.0, 3.0, -999.0, -999.0, 3.0,  3.25, -999.0, 2.8, 4.5/)
 
     PUBLIC nreact
@@ -436 +432 @@
     PUBLIC nvar                !< number of variable gas phase chemical species (nvar <= nspec)
     PUBLIC spc_names           !< names of gas phase chemical species (come from KPP) (come from KPP)
-    PUBLIC spec_conc_2  
-!    
+    PUBLIC spec_conc_2
+!
 !-- Interface section
     INTERFACE chem_3d_data_averaging
-       MODULE PROCEDURE chem_3d_data_averaging 
+       MODULE PROCEDURE chem_3d_data_averaging
     END INTERFACE chem_3d_data_averaging
 
@@ -508 +504 @@
        MODULE PROCEDURE chem_non_advective_processes_ij
     END INTERFACE chem_non_advective_processes
-    
+
     INTERFACE chem_exchange_horiz_bounds
        MODULE PROCEDURE chem_exchange_horiz_bounds
-    END INTERFACE chem_exchange_horiz_bounds    
+    END INTERFACE chem_exchange_horiz_bounds
 
     INTERFACE chem_prognostic_equations
@@ -532 +528 @@
 
     INTERFACE chem_wrd_local
-       MODULE PROCEDURE chem_wrd_local 
+       MODULE PROCEDURE chem_wrd_local
     END INTERFACE chem_wrd_local
 
     INTERFACE chem_depo
-       MODULE PROCEDURE chem_depo 
+       MODULE PROCEDURE chem_depo
     END INTERFACE chem_depo
 
     INTERFACE drydepos_gas_depac
-       MODULE PROCEDURE drydepos_gas_depac 
+       MODULE PROCEDURE drydepos_gas_depac
     END INTERFACE drydepos_gas_depac
 
     INTERFACE rc_special
-       MODULE PROCEDURE rc_special 
+       MODULE PROCEDURE rc_special
     END INTERFACE rc_special
 
     INTERFACE  rc_gw
-       MODULE PROCEDURE rc_gw 
+       MODULE PROCEDURE rc_gw
     END INTERFACE rc_gw
 
-    INTERFACE rw_so2 
-       MODULE PROCEDURE rw_so2  
+    INTERFACE rw_so2
+       MODULE PROCEDURE rw_so2
     END INTERFACE rw_so2
 
     INTERFACE rw_nh3_sutton
-       MODULE PROCEDURE rw_nh3_sutton 
+       MODULE PROCEDURE rw_nh3_sutton
     END INTERFACE rw_nh3_sutton
 
     INTERFACE rw_constant
-       MODULE PROCEDURE rw_constant 
+       MODULE PROCEDURE rw_constant
     END INTERFACE rw_constant
 
     INTERFACE rc_gstom
-       MODULE PROCEDURE rc_gstom 
+       MODULE PROCEDURE rc_gstom
     END INTERFACE rc_gstom
 
     INTERFACE rc_gstom_emb
-       MODULE PROCEDURE rc_gstom_emb 
+       MODULE PROCEDURE rc_gstom_emb
     END INTERFACE rc_gstom_emb
 
     INTERFACE par_dir_diff
-       MODULE PROCEDURE par_dir_diff 
+       MODULE PROCEDURE par_dir_diff
     END INTERFACE par_dir_diff
 
     INTERFACE rc_get_vpd
-       MODULE PROCEDURE rc_get_vpd 
+       MODULE PROCEDURE rc_get_vpd
     END INTERFACE rc_get_vpd
 
     INTERFACE rc_gsoil_eff
-       MODULE PROCEDURE rc_gsoil_eff 
+       MODULE PROCEDURE rc_gsoil_eff
     END INTERFACE rc_gsoil_eff
 
     INTERFACE rc_rinc
-       MODULE PROCEDURE rc_rinc 
+       MODULE PROCEDURE rc_rinc
     END INTERFACE rc_rinc
 
     INTERFACE rc_rctot
-       MODULE PROCEDURE rc_rctot  
+       MODULE PROCEDURE rc_rctot
     END INTERFACE rc_rctot
 
     INTERFACE drydepo_aero_zhang_vd
-       MODULE PROCEDURE drydepo_aero_zhang_vd 
+       MODULE PROCEDURE drydepo_aero_zhang_vd
     END INTERFACE drydepo_aero_zhang_vd
 
@@ -612 +608 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Subroutine for averaging 3D data of chemical species. Due to the fact that 
-!> the averaged chem arrays are allocated in chem_init, no if-query concerning
-!> the allocation is required (in any mode). Attention: If you just specify an 
-!> averaged output quantity in the _p3dr file during restarts the first output
-!> includes the time between the beginning of the restart run and the first
-!> output time (not necessarily the whole averaging_interval you have 
-!> specified in your _p3d/_p3dr file )
-!------------------------------------------------------------------------------!
+!> Subroutine for averaging 3D data of chemical species. Due to the fact that the averaged chem
+!> arrays are allocated in chem_init, no if-query concerning the allocation is required (in any
+!> mode). Attention: If you just specify an averaged output quantity in the _p3dr file during
+!> restarts the first output includes the time between the beginning of the restart run and the
+!> first output time (not necessarily the whole averaging_interval you have specified in your
+!> _p3d/_p3dr file ).
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_3d_data_averaging( mode, variable )
 
     USE control_parameters
 
-    USE exchange_horiz_mod,                                                    &
+    USE exchange_horiz_mod,                                                                        &
         ONLY:  exchange_horiz_2d
 
 
-    CHARACTER (LEN=*) ::  mode     !< 
-    CHARACTER (LEN=*) ::  variable !< 
+    CHARACTER (LEN=*) ::  mode     !<
+    CHARACTER (LEN=*) ::  variable !<
 
     LOGICAL ::  match_def  !< flag indicating default-type surface
@@ -641 +636 @@
     INTEGER(iwp) ::  j                  !< grid index y direction
     INTEGER(iwp) ::  k                  !< grid index z direction
+    INTEGER(iwp) ::  lsp                !< running index for chem spcs
     INTEGER(iwp) ::  m                  !< running index surface type
-    INTEGER(iwp) ::  lsp               !< running index for chem spcs
-
-    IF ( (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_')  )  THEN
+
+    IF ( ( variable(1:3) == 'kc_'  .OR.  variable(1:3) == 'em_' )  )  THEN
 
        IF ( mode == 'allocate' )  THEN
 
           DO  lsp = 1, nspec
-             IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
-                  TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
+             IF ( TRIM( variable(1:3) ) == 'kc_'  .AND.                                            &
+                  TRIM( variable(4:) )  == TRIM( chem_species(lsp)%name ) )  THEN
                 IF ( .NOT. ALLOCATED( chem_species(lsp)%conc_av ) )  THEN
                    ALLOCATE( chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
@@ -661 +656 @@
 
           DO  lsp = 1, nspec
-             IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
-                  TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
+             IF ( TRIM( variable(1:3) ) == 'kc_'  .AND.                                            &
+                  TRIM( variable(4:) )  == TRIM( chem_species(lsp)%name ) )  THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
-                         chem_species(lsp)%conc_av(k,j,i) =                              &
-                                           chem_species(lsp)%conc_av(k,j,i) +            &
-                                           chem_species(lsp)%conc(k,j,i)
+                         chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) +     &
+                                                            chem_species(lsp)%conc(k,j,i)
                       ENDDO
                    ENDDO
@@ -675 +669 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      match_def = surf_def_h(0)%start_index(j,i) <=                      &
-                           surf_def_h(0)%end_index(j,i)
-                      match_lsm = surf_lsm_h%start_index(j,i) <=                         &
-                           surf_lsm_h%end_index(j,i)
-                      match_usm = surf_usm_h%start_index(j,i) <=                         &
-                           surf_usm_h%end_index(j,i)
+                      match_def = surf_def_h(0)%start_index(j,i) <= surf_def_h(0)%end_index(j,i)
+                      match_lsm = surf_lsm_h%start_index(j,i)    <= surf_lsm_h%end_index(j,i)
+                      match_usm = surf_usm_h%start_index(j,i)    <= surf_usm_h%end_index(j,i)
 
                       IF ( match_def )  THEN
                          m = surf_def_h(0)%end_index(j,i)
-                         chem_species(lsp)%cssws_av(j,i) =                               &
-                              chem_species(lsp)%cssws_av(j,i) + &
-                              surf_def_h(0)%cssws(lsp,m)
+                         chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) +       &
+                                                           surf_def_h(0)%cssws(lsp,m)
                       ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
                          m = surf_lsm_h%end_index(j,i)
-                         chem_species(lsp)%cssws_av(j,i) =                               &
-                              chem_species(lsp)%cssws_av(j,i) + &
-                              surf_lsm_h%cssws(lsp,m)
+                         chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) +       &
+                                                           surf_lsm_h%cssws(lsp,m)
                       ELSEIF ( match_usm )  THEN
                          m = surf_usm_h%end_index(j,i)
-                         chem_species(lsp)%cssws_av(j,i) =                               &
-                              chem_species(lsp)%cssws_av(j,i) + &
-                              surf_usm_h%cssws(lsp,m)
+                         chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) +       &
+                                                           surf_usm_h%cssws(lsp,m)
                       ENDIF
                    ENDDO
@@ -706 +694 @@
 
           DO  lsp = 1, nspec
-             IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
-                  TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
+             IF ( TRIM( variable(1:3) ) == 'kc_'  .AND.                                            &
+                  TRIM( variable(4:) )  == TRIM( chem_species(lsp)%name ) )  THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
-                         chem_species(lsp)%conc_av(k,j,i) =                              &
-                             chem_species(lsp)%conc_av(k,j,i) /                          &
-                             REAL( average_count_3d, KIND=wp )
+                         chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) /     &
+                                                            REAL( average_count_3d, KIND = wp )
                       ENDDO
                    ENDDO
@@ -721 +708 @@
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
-                      chem_species(lsp)%cssws_av(j,i) =                                  &
-                      chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
+                      chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) /          &
+                                                        REAL( average_count_3d, KIND = wp )
                    ENDDO
                 ENDDO
@@ -734 +721 @@
  END SUBROUTINE chem_3d_data_averaging
 
-    
-!------------------------------------------------------------------------------!
+
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine to initialize and set all boundary conditions for chemical species
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_boundary_conditions( horizontal_conditions_only )
 
-    USE arrays_3d,                                                             &
+    USE arrays_3d,                                                                                 &
         ONLY:  dzu
 
-    USE surface_mod,                                                           &
-        ONLY:  bc_h                                                           
+    USE surface_mod,                                                                               &
+        ONLY:  bc_h
 
     INTEGER(iwp) ::  i                            !< grid index x direction.
@@ -752 +739 @@
     INTEGER(iwp) ::  k                            !< grid index z direction.
     INTEGER(iwp) ::  l                            !< running index boundary type, for up- and downward-facing walls.
+    INTEGER(iwp) ::  lsp                          !< running index for chem spcs.
     INTEGER(iwp) ::  m                            !< running index surface elements.
-    INTEGER(iwp) ::  lsp                          !< running index for chem spcs.
 
     LOGICAL, OPTIONAL ::  horizontal_conditions_only  !< switch to set horizontal bc only
@@ -770 +757 @@
                     j = bc_h(l)%j(m)
                     k = bc_h(l)%k(m)
-                   chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
-                                      chem_species(lsp)%conc(k+bc_h(l)%koff,j,i) 
+                   chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =                                  &
+                                                          chem_species(lsp)%conc(k+bc_h(l)%koff,j,i)
                 ENDDO
              ENDDO
@@ -785 +772 @@
                    j = bc_h(l)%j(m)
                    k = bc_h(l)%k(m)
-                   chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
-                                         chem_species(lsp)%conc_p(k,j,i)
+                   chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) = chem_species(lsp)%conc_p(k,j,i)
                 ENDDO
              ENDDO
           ENDIF
 
-       ENDDO    ! end lsp loop  
+       ENDDO    ! end lsp loop
 
 !
@@ -798 +784 @@
 !>           already runs over all species? (Siggi)
        DO  lsp = 1, nspec
-          IF ( ibc_cs_t == 0 )  THEN                    
+          IF ( ibc_cs_t == 0 )  THEN
              chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)
           ELSEIF ( ibc_cs_t == 1 )  THEN
@@ -904 +890 @@
 !
 !--    Lateral Neumann booundary conditions for timelevel t.
-!--    This branch is executed when routine is called after the non-advective processes /
-!--    before the prognostic equations.
+!--    This branch is executed when routine is called after the non-advective processes / before the
+!--    prognostic equations.
        IF ( bc_radiation_cs_s )  THEN
           DO  lsp = 1, nspec
@@ -998 +984 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine for checking data output for chemical species
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
 
@@ -1010 +996 @@
 
     INTEGER(iwp) ::  i
+    INTEGER(iwp) ::  ilen
     INTEGER(iwp) ::  lsp
-    INTEGER(iwp) ::  ilen
     INTEGER(iwp) ::  k
 
@@ -1026 +1012 @@
     IF ( TRIM( var(1:3) ) == 'em_' )  THEN
        DO  lsp=1,nspec
-          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
+          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
              unit = 'mol m-2 s-1'
           ENDIF
 !
-!--       It is possible to plant PM10 and PM25 into the gasphase chemistry code
-!--       as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
-!--       set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
+!--       It is possible to plant PM10 and PM25 into the gasphase chemistry code as passive species
+!--       (e.g. 'passive' in GASPHASE_PREPROC/mechanisms): set unit to micrograms per m**3 for PM10
+!--       and PM25 (PM2.5)
           IF (spec_name(1:2) == 'PM')  THEN
              unit = 'kg m-2 s-1'
@@ -1041 +1027 @@
 
        DO  lsp=1,nspec
-          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
+          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
              unit = 'ppm'
           ENDIF
 !
-!--            It is possible  to plant PM10 and PM25 into the gasphase chemistry code 
-!--            as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
-!--            set unit to kilograms per m**3 for PM10 and PM25 (PM2.5)
-          IF (spec_name(1:2) == 'PM')  THEN  
+!--       It is possible to plant PM10 and PM25 into the gasphase chemistry code as passive species
+!--       (e.g. 'passive' in GASPHASE_PREPROC/mechanisms): set unit to kilograms per m**3 for PM10
+!--       and PM25 (PM2.5)
+          IF (spec_name(1:2) == 'PM')  THEN
             unit = 'kg m-3'
           ENDIF
@@ -1054 +1040 @@
 
        DO  lsp=1,nphot
-          IF (TRIM( spec_name ) == TRIM( phot_frequen(lsp)%name ) )  THEN
+          IF ( TRIM( spec_name ) == TRIM( phot_frequen(lsp)%name ) )  THEN
              unit = 'sec-1'
           ENDIF
@@ -1065 +1051 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine for checking data output of profiles for chemistry model
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
 
     USE arrays_3d
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                        &
         ONLY:  data_output_pr, message_string
 
@@ -1082 +1068 @@
 
 
-    CHARACTER (LEN=*) ::  unit     !< 
-    CHARACTER (LEN=*) ::  variable !< 
-    CHARACTER (LEN=*) ::  dopr_unit
+    CHARACTER (LEN=*)  ::  unit     !<
+    CHARACTER (LEN=*)  ::  variable !<
+    CHARACTER (LEN=*)  ::  dopr_unit
     CHARACTER (LEN=16) ::  spec_name
 
@@ -1090 +1076 @@
 
 
-    spec_name = TRIM( variable(4:) )    
-
-       IF (  .NOT.  air_chemistry )  THEN
-          message_string = 'data_output_pr = ' //                        &
-          TRIM( data_output_pr(var_count) ) // ' is not imp' // &
-                      'lemented for air_chemistry = .FALSE.'
-          CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )              
+    spec_name = TRIM( variable(4:) )
+
+       IF ( .NOT. air_chemistry )  THEN
+          message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //             &
+                           ' is not imp' // 'lemented for air_chemistry = .FALSE.'
+          CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )
 
        ELSE
           DO  lsp = 1, nspec
-             IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN 
+             IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
                 cs_pr_count = cs_pr_count+1
                 cs_pr_index(cs_pr_count) = lsp
-                dopr_index(var_count) = pr_palm+cs_pr_count  
+                dopr_index(var_count) = pr_palm+cs_pr_count
                 dopr_unit  = 'ppm'
-                IF (spec_name(1:2) == 'PM')  THEN
+                IF ( spec_name(1:2) == 'PM')  THEN
                      dopr_unit = 'kg m-3'
                 ENDIF
                    hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
-                   unit = dopr_unit  
+                   unit = dopr_unit
              ENDIF
           ENDDO
@@ -1117 +1102 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Check parameters routine for chemistry_model_mod
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_check_parameters
 
@@ -1127 +1112 @@
         ONLY:  bc_lr, bc_ns
 
+    INTEGER (iwp) ::  lsp          !< running index for chem spcs.
+    INTEGER (iwp) ::  lsp_usr      !< running index for user defined chem spcs
+
     LOGICAL  ::  found
-    INTEGER (iwp) ::  lsp_usr      !< running index for user defined chem spcs 
-    INTEGER (iwp) ::  lsp          !< running index for chem spcs.
-!
-!-- check for chemical reactions status
+
+!
+!-- Check for chemical reactions status
     IF ( chem_gasphase_on )  THEN
        message_string = 'Chemical reactions: ON'
@@ -1140 +1127 @@
     ENDIF
 !
-!-- check for chemistry time-step
+!-- Check for chemistry time-step
     IF ( call_chem_at_all_substeps )  THEN
        message_string = 'Chemistry is calculated at all meteorology time-step'
@@ -1149 +1136 @@
     ENDIF
 !
-!-- check for photolysis scheme 
+!-- Check for photolysis scheme
     IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant')  )  THEN
        message_string = 'Incorrect photolysis scheme selected, please check spelling'
@@ -1155 +1142 @@
     ENDIF
 !
-!-- check for chemical mechanism used 
+!-- Check for chemical mechanism used
     CALL get_mechanism_name
     IF ( chem_mechanism /= TRIM( cs_mech ) )  THEN
-       message_string = 'Incorrect chemistry mechanism selected, check spelling in namelist and/or chem_gasphase_mod'
+       message_string =                                                                            &
+       'Incorrect chemistry mechanism selected, check spelling in namelist and/or chem_gasphase_mod'
        CALL message( 'chem_check_parameters', 'CM0462', 1, 2, 0, 6, 0 )
     ENDIF
@@ -1252 +1240 @@
 !
 !-- Cyclic conditions must be set identically at opposing boundaries
-    IF ( ( bc_cs_l == 'cyclic' .AND. bc_cs_r /= 'cyclic' )  .OR.                                 &
+    IF ( ( bc_cs_l == 'cyclic' .AND. bc_cs_r /= 'cyclic' )  .OR.                                   &
          ( bc_cs_r == 'cyclic' .AND. bc_cs_l /= 'cyclic' ) )  THEN
        message_string = 'boundary conditions bc_cs_l and bc_cs_r must both be cyclic or non-cyclic'
        CALL message( 'chem_check_parameters','PA0716', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( ( bc_cs_n == 'cyclic' .AND. bc_cs_s /= 'cyclic' )  .OR.                                 &
+    IF ( ( bc_cs_n == 'cyclic' .AND. bc_cs_s /= 'cyclic' )  .OR.                                   &
          ( bc_cs_s == 'cyclic' .AND. bc_cs_n /= 'cyclic' ) )  THEN
        message_string = 'boundary conditions bc_cs_n and bc_cs_s must both be cyclic or non-cyclic'
@@ -1298 +1286 @@
     CALL chem_init_internal
 !
-!-- check for initial chem species input
+!-- Check for initial chem species input
     lsp_usr = 1
     lsp     = 1
@@ -1307 +1295 @@
              found = .TRUE.
              EXIT
-          ENDIF 
+          ENDIF
        ENDDO
        IF ( .NOT.  found )  THEN
@@ -1317 +1305 @@
     ENDDO
 !
-!-- check for surface  emission flux chem species
+!-- Check for surface  emission flux chem species
     lsp_usr = 1
     lsp     = 1
@@ -1326 +1314 @@
              found = .TRUE.
              EXIT
-          ENDIF 
+          ENDIF
        ENDDO
        IF ( .NOT.  found )  THEN
@@ -1339 +1327 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine defining 2D output variables for chemical species
 !> @todo: Remove "mode" from argument list, not used.
-!------------------------------------------------------------------------------!
- SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf,   &
-                                  two_d, nzb_do, nzt_do, fill_value )
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do, &
+                                 fill_value )
 
 
     CHARACTER (LEN=*) ::  grid       !<
-    CHARACTER (LEN=*) ::  mode       !< 
-    CHARACTER (LEN=*) ::  variable   !< 
-    INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or time-averaged data
+    CHARACTER (LEN=*) ::  mode       !<
+    CHARACTER (LEN=*) ::  variable   !<
+
+    INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or
+                                     !< time-averaged data
     INTEGER(iwp) ::  nzb_do          !< lower limit of the domain (usually nzb)
     INTEGER(iwp) ::  nzt_do          !< upper limit of the domain (usually nzt+1)
-    LOGICAL      ::  found           !< 
-    LOGICAL      ::  two_d           !< flag parameter that indicates 2D variables (horizontal cross sections)
-    REAL(wp)     ::  fill_value 
-    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf  
+
+    LOGICAL      ::  found           !<
+    LOGICAL      ::  two_d           !< flag parameter that indicates 2D variables (horizontal cross
+                                     !< sections)
+
+    REAL(wp)     ::  fill_value
+
+    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf
 
 !
@@ -1378 +1372 @@
     spec_name = TRIM( variable(4:char_len-3) )
 !
-!-- Output of emission values, i.e. surface fluxes cssws. 
+!-- Output of emission values, i.e. surface fluxes cssws.
     IF ( variable(1:3) == 'em_' )  THEN
 
@@ -1386 +1380 @@
           IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
 !
-!--          No average output for now. 
+!--          No average output for now.
              DO  m = 1, surf_lsm_h%ns
-                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1) =              &
-                              local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1)  &
-                                            + surf_lsm_h%cssws(lsp,m)
+                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1) =                                  &
+                                                   local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1) &
+                                                   + surf_lsm_h%cssws(lsp,m)
              ENDDO
              DO  m = 1, surf_usm_h%ns
-                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1) =           &
-                              local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1)  &
-                                            + surf_usm_h%cssws(lsp,m)
+                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1) =                               &
+                                                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1) &
+                                                   + surf_usm_h%cssws(lsp,m)
              ENDDO
              grid = 'zu'
@@ -1405 +1399 @@
 
        DO  lsp=1,nspec
-          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name )  .AND.                           &
+          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name )  .AND.                          &
                 ( (variable(char_len-2:) == '_xy')  .OR.                                           &
                   (variable(char_len-2:) == '_xz')  .OR.                                           &
-                  (variable(char_len-2:) == '_yz') ) )  THEN             
+                  (variable(char_len-2:) == '_yz') ) )  THEN
 !
 !--   todo: remove or replace by "CALL message" mechanism (kanani)
 !                    IF(myid == 0)  WRITE(6,*) 'Output of species ' // TRIM( variable )  //       &
-!                                                             TRIM( chem_species(lsp)%name )       
+!                                                             TRIM( chem_species(lsp)%name )
              IF (av == 0)  THEN
                 DO  i = nxl, nxr
@@ -1437 +1431 @@
                 ENDDO
              ENDIF
-             grid = 'zu'           
+             grid = 'zu'
              found = .TRUE.
           ENDIF
@@ -1445 +1439 @@
     RETURN
 
- END SUBROUTINE chem_data_output_2d      
-
-
-!------------------------------------------------------------------------------!
+ END SUBROUTINE chem_data_output_2d
+
+
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine defining 3D output variables for chemical species
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
 
@@ -1459 +1453 @@
 
     CHARACTER (LEN=*)    ::  variable     !<
-    INTEGER(iwp)         ::  av           !<
-    INTEGER(iwp) ::  nzb_do               !< lower limit of the data output (usually 0)
-    INTEGER(iwp) ::  nzt_do               !< vertical upper limit of the data output (usually nz_do3d)
+
+    INTEGER(iwp)         ::  av         !<
+    INTEGER(iwp) ::  nzb_do             !< lower limit of the data output (usually 0)
+    INTEGER(iwp) ::  nzt_do             !< vertical upper limit of the data output (usually nz_do3d)
 
     LOGICAL      ::  found                !<
 
     REAL(wp)             ::  fill_value   !<
-    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf 
-!
-!-- local variables
+
+    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf
+!
+!-- Local variables
     CHARACTER(LEN=16)    ::  spec_name
+
     INTEGER(iwp)         ::  i
     INTEGER(iwp)         ::  j
@@ -1489 +1486 @@
        DO  lsp = 1, nvar   !!! cssws - nvar species, chem_species - nspec species !!!
           IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
-          
+
              local_pf = 0.0_wp
 !
 !--          no average for now
              DO  m = 1, surf_usm_h%ns
-                local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = &
-                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m)
+                local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) =                        &
+                                         local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) &
+                                         + surf_usm_h%cssws(lsp,m)
              ENDDO
              DO  m = 1, surf_lsm_h%ns
-                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = &
-                  local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) + surf_lsm_h%cssws(lsp,m)
+                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) =                        &
+                                         local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) &
+                                         + surf_lsm_h%cssws(lsp,m)
              ENDDO
              DO  l = 0, 3
                 DO  m = 1, surf_usm_v(l)%ns
-                   local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = &
-                     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)
+                   local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) =            &
+                                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)
                 ENDDO
                 DO  m = 1, surf_lsm_v(l)%ns
-                   local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = &
-                      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)
+                   local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) =            &
+                                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)
                 ENDDO
              ENDDO
@@ -1516 +1517 @@
     ELSE
       DO  lsp = 1, nspec
-         IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
+         IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
 !
 !--   todo: remove or replace by "CALL message" mechanism (kanani)
 !              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
-!                                                           TRIM( chem_species(lsp)%name )       
+!                                                           TRIM( chem_species(lsp)%name )
             IF (av == 0)  THEN
                DO  i = nxl, nxr
                   DO  j = nys, nyn
                      DO  k = nzb_do, nzt_do
-                         local_pf(i,j,k) = MERGE(                             &
-                                             chem_species(lsp)%conc(k,j,i),   &
-                                             REAL( fill_value, KIND = wp ),   &
+                         local_pf(i,j,k) = MERGE(                                                  &
+                                             chem_species(lsp)%conc(k,j,i),                        &
+                                             REAL( fill_value, KIND = wp ),                        &
                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
                      ENDDO
@@ -1538 +1539 @@
                   DO  j = nys, nyn
                      DO  k = nzb_do, nzt_do
-                         local_pf(i,j,k) = MERGE(                             &
-                                             chem_species(lsp)%conc_av(k,j,i),&
-                                             REAL( fill_value, KIND = wp ),   &
+                         local_pf(i,j,k) = MERGE(                                                  &
+                                             chem_species(lsp)%conc_av(k,j,i),                     &
+                                             REAL( fill_value, KIND = wp ),                        &
                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
                      ENDDO
@@ -1556 +1557 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine defining mask output variables for chemical species
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_data_output_mask( av, variable, found, local_pf, mid )
 
@@ -1569 +1570 @@
     CHARACTER(LEN=*)  ::  variable    !<
 
-    INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or time-averaged data
+    INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or
+                                     !< time-averaged data
+    INTEGER(iwp) ::  i               !< grid index along x-direction
+    INTEGER(iwp) ::  im              !< loop index for masked variables
+    INTEGER(iwp) ::  j               !< grid index along y-direction
+    INTEGER(iwp) ::  jm              !< loop index for masked variables
+    INTEGER(iwp) ::  k               !< grid index along z-direction
+    INTEGER(iwp) ::  kk              !< masked output index along z-direction
+    INTEGER(iwp) ::  ktt             !< k index of highest terrain surface
     INTEGER(iwp) ::  lsp
-    INTEGER(iwp) ::  i               !< grid index along x-direction
-    INTEGER(iwp) ::  j               !< grid index along y-direction
-    INTEGER(iwp) ::  k               !< grid index along z-direction
-    INTEGER(iwp) ::  im              !< loop index for masked variables
-    INTEGER(iwp) ::  jm              !< loop index for masked variables
-    INTEGER(iwp) ::  kk              !< masked output index along z-direction
     INTEGER(iwp) ::  mid             !< masked output running index
-    INTEGER(iwp) ::  ktt             !< k index of highest terrain surface
-    
+
     LOGICAL ::  found
-    
-    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
-              local_pf   !< 
+
+    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  local_pf   !<
 
     REAL(wp), PARAMETER ::  fill_value = -9999.0_wp    !< value for the _FillValue attribute
 
 !
-!-- local variables.
+!-- Local variables.
 
     spec_name = TRIM( variable(4:) )
@@ -1594 +1595 @@
 
     DO  lsp=1,nspec
-       IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN             
+       IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
 !
 !-- todo: remove or replace by "CALL message" mechanism (kanani)
 !              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
-!                                                        TRIM( chem_species(lsp)%name )       
+!                                                        TRIM( chem_species(lsp)%name )
           IF (av == 0)  THEN
              IF ( .NOT. mask_surface(mid) )  THEN
 
                 DO  i = 1, mask_size_l(mid,1)
-                   DO  j = 1, mask_size_l(mid,2) 
-                      DO  k = 1, mask_size(mid,3) 
-                          local_pf(i,j,k) = chem_species(lsp)%conc(  &
-                                               mask_k(mid,k),        &
-                                               mask_j(mid,j),        &
-                                               mask_i(mid,i)      )
+                   DO  j = 1, mask_size_l(mid,2)
+                      DO  k = 1, mask_size(mid,3)
+                          local_pf(i,j,k) = chem_species(lsp)%conc( mask_k(mid,k),                 &
+                                                                    mask_j(mid,j),                 &
+                                                                    mask_i(mid,i)      )
                       ENDDO
                    ENDDO
@@ -1614 +1614 @@
 
              ELSE
-!              
+!
 !--             Terrain-following masked output
                 DO  i = 1, mask_size_l(mid,1)
                    DO  j = 1, mask_size_l(mid,2)
+!
 !--                   Get k index of the highest terraing surface
                       im = mask_i(mid,i)
                       jm = mask_j(mid,j)
-                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), DIM = 1 ) - 1
+                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 ) ),        &
+                                    DIM = 1 ) - 1
                       DO  k = 1, mask_size_l(mid,3)
                          kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
 !--                      Set value if not in building
                          IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
@@ -1641 +1644 @@
                    DO  j = 1, mask_size_l(mid,2)
                       DO  k =  1, mask_size_l(mid,3)
-                          local_pf(i,j,k) = chem_species(lsp)%conc_av(  &
-                                               mask_k(mid,k),           &
-                                               mask_j(mid,j),           &
-                                               mask_i(mid,i)         )
+                          local_pf(i,j,k) = chem_species(lsp)%conc_av(  mask_k(mid,k),             &
+                                                                        mask_j(mid,j),             &
+                                                                        mask_i(mid,i)         )
                       ENDDO
                    ENDDO
@@ -1650 +1652 @@
 
              ELSE
-!              
+!
 !--             Terrain-following masked output
                 DO  i = 1, mask_size_l(mid,1)
                    DO  j = 1, mask_size_l(mid,2)
+!
 !--                   Get k index of the highest terraing surface
                       im = mask_i(mid,i)
                       jm = mask_j(mid,j)
-                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), DIM = 1 ) - 1
+                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )),         &
+                                    DIM = 1 ) - 1
                       DO  k = 1, mask_size_l(mid,3)
                          kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
 !--                      Set value if not in building
                          IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
@@ -1680 +1685 @@
     RETURN
 
- END SUBROUTINE chem_data_output_mask      
-
-
-!------------------------------------------------------------------------------!
+ END SUBROUTINE chem_data_output_mask
+
+
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine defining appropriate grid for netcdf variables.
 !> It is called out from subroutine netcdf.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
 
 
     CHARACTER (LEN=*), INTENT(IN)  ::  var          !<
-    LOGICAL, INTENT(OUT)           ::  found        !<
+
     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x       !<
     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y       !<
     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z       !<
 
+    LOGICAL, INTENT(OUT)           ::  found        !<
+
     found  = .TRUE.
 
-    IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' )  THEN                   !< always the same grid for chemistry variables
+    IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' )  THEN                !< always the same grid for
+                                                                         !< chemistry variables
        grid_x = 'x'
        grid_y = 'y'
-       grid_z = 'zu'                             
+       grid_z = 'zu'
     ELSE
        found  = .FALSE.
@@ -1715 +1723 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine defining header output for chemistry model
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_header( io )
 
-
-    INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
-    INTEGER(iwp)  :: lsp                       !< running index for chem spcs
-    INTEGER(iwp)  :: cs_fixed 
     CHARACTER (LEN=80)  :: docsflux_chr
     CHARACTER (LEN=80)  :: docsinit_chr
-!
-! Get name of chemical mechanism from chem_gasphase_mod
+
+    INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
+
+    INTEGER(iwp)  :: cs_fixed
+    INTEGER(iwp)  :: lsp                       !< running index for chem spcs
+
+!
+!-- Get name of chemical mechanism from chem_gasphase_mod
     CALL get_mechanism_name
 !
@@ -1736 +1746 @@
 !
 !-- Gasphase reaction status
-    IF ( chem_gasphase_on )  THEN 
+    IF ( chem_gasphase_on )  THEN
        WRITE( io, 2 )
     ELSE
@@ -1746 +1756 @@
 !
 !-- Emission mode info
-!-- At the moment the evaluation is done with both emiss_lod and mode_emis
-!-- but once salsa has been migrated to emiss_lod the .OR. mode_emis
-!-- conditions can be removed (ecc 20190513)
-    IF     ( (emiss_lod == 1) .OR. (mode_emis == 'DEFAULT') )        THEN
+!-- At the moment the evaluation is done with both emiss_lod and mode_emis but once salsa has been
+!-- migrated to emiss_lod the .OR. mode_emis conditions can be removed (ecc 20190513)
+    IF     ( ( emiss_lod == 1 )  .OR.  ( mode_emis == 'DEFAULT' ) )        THEN
         WRITE ( io, 5 )
-    ELSEIF ( (emiss_lod == 0) .OR. (mode_emis == 'PARAMETERIZED') )  THEN 
+    ELSEIF ( ( emiss_lod == 0 )  .OR.  ( mode_emis == 'PARAMETERIZED' ) )  THEN
         WRITE ( io, 6 )
-    ELSEIF ( (emiss_lod == 2) .OR. (mode_emis == 'PRE-PROCESSED') )  THEN 
+    ELSEIF ( ( emiss_lod == 2 )  .OR.  ( mode_emis == 'PRE-PROCESSED' ) )  THEN
         WRITE ( io, 7 )
     ENDIF
@@ -1759 +1768 @@
 !-- Photolysis scheme info
     IF ( photolysis_scheme == "simple" )  THEN
-       WRITE( io, 8 ) 
+       WRITE( io, 8 )
     ELSEIF (photolysis_scheme == "constant" )  THEN
        WRITE( io, 9 )
@@ -1766 +1775 @@
 !-- Emission flux info
     lsp = 1
-    docsflux_chr ='Chemical species for surface emission flux: ' 
+    docsflux_chr ='Chemical species for surface emission flux: '
     DO WHILE ( surface_csflux_name(lsp) /= 'novalue' )
-       docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ','  
+       docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ','
        IF ( LEN_TRIM( docsflux_chr ) >= 75 )  THEN
           WRITE ( io, 10 )  docsflux_chr
@@ -1780 +1789 @@
     ENDIF
 !
-!-- initializatoin of Surface and profile chemical species
+!-- Initialization of Surface and profile chemical species
     lsp = 1
-    docsinit_chr ='Chemical species for initial surface and profile emissions: ' 
+    docsinit_chr ='Chemical species for initial surface and profile emissions: '
     DO WHILE ( cs_name(lsp) /= 'novalue' )
-       docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ','  
+       docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ','
        IF ( LEN_TRIM( docsinit_chr ) >= 75 )  THEN
-        WRITE ( io, 11 )  docsinit_chr
-        docsinit_chr = '       '
+          WRITE ( io, 11 )  docsinit_chr
+          docsinit_chr = '       '
        ENDIF
        lsp = lsp + 1
@@ -1803 +1812 @@
 
 !
-!-- number of variable and fix chemical species and number of reactions
+!-- Number of variable and fix chemical species and number of reactions
     cs_fixed = nspec - nvar
-    WRITE ( io, * ) '   --> Chemical Mechanism        : ', cs_mech 
+    WRITE ( io, * ) '   --> Chemical Mechanism        : ', cs_mech
     WRITE ( io, * ) '   --> Chemical species, variable: ', nvar
     WRITE ( io, * ) '   --> Chemical species, fixed   : ', cs_fixed
@@ -1811 +1820 @@
 
 
-1   FORMAT (//' Chemistry model information:'/                                  &
-           ' ----------------------------'/)
+1   FORMAT (//' Chemistry model information:'/' ----------------------------'/)
 2   FORMAT ('    --> Chemical reactions are turned on')
 3   FORMAT ('    --> Chemical reactions are turned off')
@@ -1821 +1829 @@
 8   FORMAT ('    --> Photolysis scheme used =  simple ')
 9   FORMAT ('    --> Photolysis scheme used =  constant ')
-10  FORMAT (/'    ',A)  
-11  FORMAT (/'    ',A) 
+10  FORMAT (/'    ',A)
+11  FORMAT (/'    ',A)
 12  FORMAT (/'   Nesting for chemistry variables: ', L1 )
 13  FORMAT (/'   Offline nesting for chemistry variables: ', L1 )
@@ -1833 +1841 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine initializating chemistry_model_mod specific arrays
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_init_arrays
 !
@@ -1845 +1853 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine initializating chemistry_model_mod
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_init
 
@@ -1856 +1864 @@
 !-- introduced additional interfaces for on-demand emission update
 
-!    USE chem_emissions_mod,                                                    &
+!    USE chem_emissions_mod,                                                                        &
 !        ONLY:  chem_emissions_init
-        
-    USE chem_emissions_mod,                                                    &
-        ONLY:  chem_emissions_init, chem_emissions_header_init
-        
-    USE netcdf_data_input_mod,                                                 &
+
+    USE chem_emissions_mod,                                                                        &
+        ONLY:  chem_emissions_header_init, chem_emissions_init
+
+    USE netcdf_data_input_mod,                                                                     &
         ONLY:  init_3d
 
@@ -1874 +1882 @@
 !
 !-- Next statement is to avoid compiler warning about unused variables
-    IF ( ( ilu_arable + ilu_coniferous_forest + ilu_deciduous_forest + ilu_mediterrean_scrub + &
-           ilu_permanent_crops + ilu_savanna + ilu_semi_natural_veg + ilu_tropical_forest +    &
+    IF ( ( ilu_arable + ilu_coniferous_forest + ilu_deciduous_forest + ilu_mediterrean_scrub +     &
+           ilu_permanent_crops + ilu_savanna + ilu_semi_natural_veg + ilu_tropical_forest +        &
            ilu_urban ) == 0 )  CONTINUE
 
 !
 !-- 20200203 (ECC)
-!-- calls specific emisisons initialization subroutines
-!-- for legacy mode and on-demand mode         
+!-- Calls specific emisisons initialization subroutines for legacy mode and on-demand mode
 
 !    IF ( emissions_anthropogenic )  CALL chem_emissions_init
@@ -1897 +1904 @@
 
 !
-!-- Chemistry variables will be initialized if availabe from dynamic 
-!-- input file. Note, it is possible to initialize only part of the chemistry
-!-- variables from dynamic input. 
+!-- Chemistry variables will be initialized if availabe from dynamic input file. Note, it is
+!-- possible to initialize only part of the chemistry variables from dynamic input.
     IF ( INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
        DO  n = 1, nspec
@@ -1917 +1923 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine initializating chemistry_model_mod
 !> internal workaround for chem_species dependency in chem_check_parameters
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_init_internal
- 
+
     USE pegrid
 
-    USE netcdf_data_input_mod,                                                 &
-        ONLY:  chem_emis, chem_emis_att, input_pids_dynamic, init_3d,          &
+    USE netcdf_data_input_mod,                                                                     &
+        ONLY:  chem_emis, chem_emis_att, input_pids_dynamic, init_3d,                              &
                netcdf_data_input_chemistry_data
 
@@ -1935 +1941 @@
     INTEGER(iwp) ::  i                 !< running index for for horiz numerical grid points
     INTEGER(iwp) ::  j                 !< running index for for horiz numerical grid points
+    INTEGER(iwp) ::  lpr_lev           !< running index for chem spcs profile level
     INTEGER(iwp) ::  lsp               !< running index for chem spcs
-    INTEGER(iwp) ::  lpr_lev           !< running index for chem spcs profile level
 
     REAL(wp)     ::  flag              !< flag for masking topography/building grid points
@@ -1959 +1965 @@
     ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
     ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
-    ALLOCATE( phot_frequen(nphot) ) 
+    ALLOCATE( phot_frequen(nphot) )
     ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) )
     ALLOCATE( bc_cs_t_val(nspec) )
@@ -1976 +1982 @@
        chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_3 (:,:,:,lsp)
 
-       ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))                   
+       ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))
        chem_species(lsp)%cssws_av    = 0.0_wp
 !
 !--    The following block can be useful when emission module is not applied. &
-!--    if emission module is applied the following block will be overwritten.
-       ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))   
-       ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))   
-       ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1)) 
-       ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))   
-       chem_species(lsp)%flux_s_cs = 0.0_wp                                     
-       chem_species(lsp)%flux_l_cs = 0.0_wp                                     
-       chem_species(lsp)%diss_s_cs = 0.0_wp                                      
-       chem_species(lsp)%diss_l_cs = 0.0_wp                                     
-!
-!--   Allocate memory for initial concentration profiles
-!--   (concentration values come from namelist)
-!--   (@todo (FK): Because of this, chem_init is called in palm before
-!--               check_parameters, since conc_pr_init is used there.
-!--               We have to find another solution since chem_init should
-!--               eventually be called from init_3d_model!!)
+!--    If emission module is applied the following block will be overwritten.
+       ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))
+       ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))
+       ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
+       ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
+       chem_species(lsp)%flux_s_cs = 0.0_wp
+       chem_species(lsp)%flux_l_cs = 0.0_wp
+       chem_species(lsp)%diss_s_cs = 0.0_wp
+       chem_species(lsp)%diss_l_cs = 0.0_wp
+!
+!--   Allocate memory for initial concentration profiles (concentration values come from namelist)
+!--   (@todo (FK): Because of this, chem_init is called in palm before check_parameters, since
+!--                conc_pr_init is used there.
+!--                We have to find another solution since chem_init should eventually be called from
+!--                init_3d_model!!)
        ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) )
        chem_species(lsp)%conc_pr_init(:) = 0.0_wp
@@ -2002 +2007 @@
 
 !
-!-- For some passive scalars decycling may be enabled. This case, the lateral 
-!-- boundary conditions are non-cyclic for these scalars (chemical species
-!-- and aerosols), while the other scalars may have
-!-- cyclic boundary conditions. However, large gradients near the boundaries
-!-- may produce stationary numerical oscillations near the lateral boundaries
-!-- when a higher-order scheme is applied near these boundaries. 
-!-- To get rid-off this, set-up additional flags that control the order of the 
-!-- scalar advection scheme near the lateral boundaries for passive scalars 
-!-- with decycling.
+!-- For some passive scalars decycling may be enabled. This case, the lateral boundary conditions
+!-- are non-cyclic for these scalars (chemical species and aerosols), while the other scalars may
+!-- have cyclic boundary conditions. However, large gradients near the boundaries may produce
+!-- stationary numerical oscillations near the lateral boundaries when a higher-order scheme is
+!-- applied near these boundaries.
+!-- To get rid-off this, set-up additional flags that control the order of the scalar advection
+!-- scheme near the lateral boundaries for passive scalars with decycling.
     IF ( scalar_advec == 'ws-scheme' )  THEN
        ALLOCATE( cs_advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 !
-!--    In case of decyling, set Neumann boundary conditions for wall_flags_total_0 
-!--    bit 31 instead of cyclic boundary conditions.
-!--    Bit 31 is used to identify extended degradation zones (please see
-!--    following comment). 
-!--    Note, since several also other modules like Salsa or other future 
-!--    one may access this bit but may have other boundary conditions, the 
-!--    original value of wall_flags_total_0 bit 31 must not be modified. Hence, 
-!--    store the boundary conditions directly on cs_advc_flags_s. 
-!--    cs_advc_flags_s will be later overwritten in ws_init_flags_scalar and
-!--    bit 31 won't be used to control the numerical order. 
-!--    Initialize with flag 31 only. 
+!--    In case of decyling, set Neumann boundary conditions for wall_flags_total_0 bit 31 instead of
+!--    cyclic boundary conditions.
+!--    Bit 31 is used to identify extended degradation zones (please see following comment).
+!--    Note, since several also other modules like Salsa or other future one may access this bit but
+!--    may have other boundary conditions, the original value of wall_flags_total_0 bit 31 must not
+!--    be modified. Hence, store the boundary conditions directly on cs_advc_flags_s.
+!--    cs_advc_flags_s will be later overwritten in ws_init_flags_scalar and bit 31 won't be used to
+!--    control the numerical order.
+!--    Initialize with flag 31 only.
        cs_advc_flags_s = 0
        cs_advc_flags_s = MERGE( IBSET( cs_advc_flags_s, 31 ), 0, BTEST( wall_flags_total_0, 31 ) )
@@ -2031 +2032 @@
           IF ( nys == 0  )  THEN
              DO  i = 1, nbgp
-                cs_advc_flags_s(:,nys-i,:) = MERGE(                            &
-                                        IBSET( cs_advc_flags_s(:,nys,:), 31 ), &
-                                        IBCLR( cs_advc_flags_s(:,nys,:), 31 ), &
-                                        BTEST( cs_advc_flags_s(:,nys,:), 31 )  &
+                cs_advc_flags_s(:,nys-i,:) = MERGE(                                                &
+                                                    IBSET( cs_advc_flags_s(:,nys,:), 31 ),         &
+                                                    IBCLR( cs_advc_flags_s(:,nys,:), 31 ),         &
+                                                    BTEST( cs_advc_flags_s(:,nys,:), 31 )          &
                                                   )
              ENDDO
@@ -2040 +2041 @@
           IF ( nyn == ny )  THEN
              DO  i = 1, nbgp
-                cs_advc_flags_s(:,nyn+i,:) = MERGE(                            &
-                                        IBSET( cs_advc_flags_s(:,nyn,:), 31 ), &
-                                        IBCLR( cs_advc_flags_s(:,nyn,:), 31 ), &
-                                        BTEST( cs_advc_flags_s(:,nyn,:), 31 )  &
+                cs_advc_flags_s(:,nyn+i,:) = MERGE(                                                &
+                                                    IBSET( cs_advc_flags_s(:,nyn,:), 31 ),         &
+                                                    IBCLR( cs_advc_flags_s(:,nyn,:), 31 ),         &
+                                                    BTEST( cs_advc_flags_s(:,nyn,:), 31 )          &
                                                   )
              ENDDO
@@ -2051 +2052 @@
           IF ( nxl == 0  )  THEN
              DO  i = 1, nbgp
-                cs_advc_flags_s(:,:,nxl-i) = MERGE(                            &
-                                        IBSET( cs_advc_flags_s(:,:,nxl), 31 ), &
-                                        IBCLR( cs_advc_flags_s(:,:,nxl), 31 ), &
-                                        BTEST( cs_advc_flags_s(:,:,nxl), 31 )  &
+                cs_advc_flags_s(:,:,nxl-i) = MERGE(                                                &
+                                                    IBSET( cs_advc_flags_s(:,:,nxl), 31 ),         &
+                                                    IBCLR( cs_advc_flags_s(:,:,nxl), 31 ),         &
+                                                    BTEST( cs_advc_flags_s(:,:,nxl), 31 )          &
                                                   )
              ENDDO
@@ -2060 +2061 @@
           IF ( nxr == nx )  THEN
              DO  i = 1, nbgp
-                cs_advc_flags_s(:,:,nxr+i) = MERGE(                            &
-                                        IBSET( cs_advc_flags_s(:,:,nxr), 31 ), &
-                                        IBCLR( cs_advc_flags_s(:,:,nxr), 31 ), &
-                                        BTEST( cs_advc_flags_s(:,:,nxr), 31 )  &
+                cs_advc_flags_s(:,:,nxr+i) = MERGE(                                                &
+                                                    IBSET( cs_advc_flags_s(:,:,nxr), 31 ), &
+                                                    IBCLR( cs_advc_flags_s(:,:,nxr), 31 ), &
+                                                    BTEST( cs_advc_flags_s(:,:,nxr), 31 )  &
                                                   )
              ENDDO
@@ -2069 +2070 @@
        ENDIF
 !
-!--    To initialize advection flags appropriately, pass the boundary flags. 
-!--    The last argument indicates that a passive scalar is treated, where
-!--    the horizontal advection terms are degraded already 2 grid points before
-!--    the lateral boundary to avoid stationary oscillations at large-gradients. 
-!--    Also, extended degradation zones are applied, where horizontal advection of 
-!--    passive scalars is discretized by first-order scheme at all grid points 
-!--    that in the vicinity of buildings (<= 3 grid points). Even though no 
-!--    building is within the numerical stencil, first-order scheme is used. 
+!--    To initialize advection flags appropriately, pass the boundary flags.
+!--    The last argument indicates that a passive scalar is treated, where the horizontal advection
+!--    terms are degraded already 2 grid points before the lateral boundary to avoid stationary
+!--    oscillations at large-gradients.
+!--    Also, extended degradation zones are applied, where horizontal advection of passive scalars
+!--    is discretized by first-order scheme at all grid points that in the vicinity of buildings
+!--    (<= 3 grid points). Even though no building is within the numerical stencil, first-order
+!--    scheme is used.
 !--    At fourth and fifth grid point the order of the horizontal advection scheme
 !--    is successively upgraded.
-!--    These extended degradation zones are used to avoid stationary numerical 
-!--    oscillations, which are responsible for high concentration maxima that may
-!--    appear under shear-free stable conditions. 
+!--    These extended degradation zones are used to avoid stationary numerical oscillations, which
+!--    are responsible for high concentration maxima that may appear under shear-free stable
+!--    conditions.
        CALL ws_init_flags_scalar( bc_dirichlet_cs_l  .OR.  bc_radiation_cs_l,                      &
                                   bc_dirichlet_cs_n  .OR.  bc_radiation_cs_n,                      &
@@ -2089 +2090 @@
     ENDIF
 !
-!-- Initial concentration of profiles is prescribed by parameters cs_profile
-!-- and cs_heights in the namelist &chemistry_parameters
-
+!-- Initial concentration of profiles is prescribed by parameters cs_profile and cs_heights in the
+!-- namelist &chemistry_parameters.
     CALL chem_init_profiles
-!    
-!-- In case there is dynamic input file, create a list of names for chemistry
-!-- initial input files. Also, initialize array that indicates whether the
-!-- respective variable is on file or not. 
-    IF ( input_pids_dynamic )  THEN    
+!
+!-- In case there is dynamic input file, create a list of names for chemistry initial input files.
+!-- Also, initialize array that indicates whether the respective variable is on file or not.
+    IF ( input_pids_dynamic )  THEN
        ALLOCATE( init_3d%var_names_chem(1:nspec) )
        ALLOCATE( init_3d%from_file_chem(1:nspec) )
        init_3d%from_file_chem(:) = .FALSE.
-       
+
        DO  lsp = 1, nspec
           init_3d%var_names_chem(lsp) = init_3d%init_char // TRIM( chem_species(lsp)%name )
@@ -2108 +2107 @@
 !
 !-- Initialize model variables
-    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
+    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                                &
          TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
 !
@@ -2126 +2125 @@
                    ENDDO
                 ENDDO
-             ENDDO   
+             ENDDO
           ENDDO
 
-       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
-       THEN
-
-          DO  lsp = 1, nspec 
+       ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles') /= 0 )  THEN
+
+          DO  lsp = 1, nspec
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
@@ -2147 +2145 @@
 !
 !--    If required, change the surface chem spcs at the start of the 3D run
-       IF ( cs_surface_initial_change(1) /= 0.0_wp )  THEN            
-          DO  lsp = 1, nspec 
+       IF ( cs_surface_initial_change(1) /= 0.0_wp )  THEN
+          DO  lsp = 1, nspec
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
                    flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb,j,i), 0 ) )
-                   chem_species(lsp)%conc(nzb,j,i) = chem_species(lsp)%conc(nzb,j,i) +  &
+                   chem_species(lsp)%conc(nzb,j,i) = chem_species(lsp)%conc(nzb,j,i) +             &
                                                      cs_surface_initial_change(lsp) * flag
                 ENDDO
              ENDDO
           ENDDO
-       ENDIF 
+       ENDIF
 
     ENDIF
 !
-!-- Initiale old and new time levels. Note, this has to be done also in restart runs
+!-- Initial old and new time levels. Note, this has to be done also in restart runs
     DO  lsp = 1, nvar
-       chem_species(lsp)%tconc_m = 0.0_wp                      
-       chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
+       chem_species(lsp)%tconc_m = 0.0_wp
+       chem_species(lsp)%conc_p  = chem_species(lsp)%conc
     ENDDO
 
@@ -2173 +2171 @@
 !--       IF( myid == 0 )  THEN
 !--          WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,TRIM( phot_names(lsp) )
-!--       ENDIF 
+!--       ENDIF
        phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  =>  freq_1(:,:,:,lsp)
     ENDDO
@@ -2184 +2182 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Subroutine defining initial vertical profiles of chemical species (given by
-!> namelist parameters chem_profiles and chem_heights)  --> which should work
-!> analogue to parameters u_profile, v_profile and uv_heights)
-!------------------------------------------------------------------------------!
- SUBROUTINE chem_init_profiles              
+!> Subroutine defining initial vertical profiles of chemical species (given by namelist parameters
+!> chem_profiles and chem_heights)  --> which should work analogically to parameters u_profile,
+!> v_profile and uv_heights)
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE chem_init_profiles
 
     USE chem_modules
@@ -2197 +2195 @@
 !
 !-- Local variables
+    INTEGER ::  lpr_lev    !< running index for profile level for each chem spcs.
     INTEGER ::  lsp        !< running index for number of species in derived data type species_def
-    INTEGER ::  lsp_usr     !< running index for number of species (user defined)  in cs_names, cs_profiles etc
-    INTEGER ::  lpr_lev    !< running index for profile level for each chem spcs. 
+    INTEGER ::  lsp_usr    !< running index for number of species (user defined)  in cs_names,
+                           !< cs_profiles etc
     INTEGER ::  npr_lev    !< the next available profile lev
 !
@@ -2207 +2206 @@
 !-- "cs_heights" are prescribed, their values will!override the constant profile given by
 !-- "cs_surface".
-!     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN 
+!     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
     lsp_usr = 1
     DO  WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )   !'novalue' is the default
        DO  lsp = 1, nspec                                !
 !
-!--       create initial profile (conc_pr_init) for each chemical species
-          IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN   !
+!--       Create initial profile (conc_pr_init) for each chemical species
+          IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN
              IF ( cs_profile(lsp_usr,1) == 9999999.9_wp )  THEN
 !
-!--            set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
+!--            Set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of
+!--            each species
                 DO lpr_lev = 0, nzt+1
                    chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
@@ -2236 +2236 @@
                          IF ( npr_lev == 100 )  THEN
                             message_string = 'number of chem spcs exceeding the limit'
-                            CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )               
+                            CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )
                             EXIT
                          ENDIF
@@ -2242 +2242 @@
                    ENDIF
                    IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr,npr_lev+1) /= 9999999.9_wp )  THEN
-                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +       &
-                           ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                          &
-                           ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *  &
+                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +     &
+                           ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                        &
+                           ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *&
                            ( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
                    ELSE
@@ -2252 +2252 @@
              ENDIF
 !
-!--       If a profile is prescribed explicity using cs_profiles and cs_heights, then  
-!--       chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
-!--       on the cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:). 
+!--       If a profile is prescribed explicity using cs_profiles and cs_heights, then
+!--       chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based on the
+!--       cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:).
           ENDIF
 
@@ -2265 +2265 @@
  END SUBROUTINE chem_init_profiles
 
- 
-!------------------------------------------------------------------------------!
+
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine to integrate chemical species in the given chemical mechanism
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_integrate_ij( i, j )
 
-    USE statistics,                                                          &
+    USE statistics,                                                                                &
         ONLY:  weight_pres
 
-    USE control_parameters,                                                  &
+    USE control_parameters,                                                                        &
         ONLY:  dt_3d, intermediate_timestep_count, time_since_reference_point
 
@@ -2283 +2283 @@
     INTEGER,INTENT(IN)       :: j
 !
-!--   local variables
-    INTEGER(iwp) ::  lsp                                                     !< running index for chem spcs.
-    INTEGER(iwp) ::  lph                                                     !< running index for photolysis frequencies
-    INTEGER, DIMENSION(20)    :: istatus
+!-- Local variables
+    INTEGER(iwp) ::  lph                                 !< running index for photolysis frequencies
+    INTEGER(iwp) ::  lsp                                 !< running index for chem spcs.
+
+    INTEGER, DIMENSION(20)        :: istatus
+
+    INTEGER,DIMENSION(nzb+1:nzt)  :: nacc          !< Number of accepted steps
+    INTEGER,DIMENSION(nzb+1:nzt)  :: nrej          !< Number of rejected steps
+
+    REAL(wp)                         ::  conv                                !< conversion factor
+    REAL(kind=wp)                    ::  dt_chem
+
+    REAL(wp), PARAMETER              ::  fr2ppm  = 1.0e6_wp              !< Conversion factor
+                                                                         !< fraction to ppm
+!    REAL(wp), PARAMETER              ::  xm_air  = 28.96_wp              !< Mole mass of dry air
+!    REAL(wp), PARAMETER              ::  xm_h2o  = 18.01528_wp           !< Mole mass of water vapor
+    REAL(wp), PARAMETER              ::  p_std   = 101325.0_wp           !< standard pressure (Pa)
+    REAL(wp), PARAMETER              ::  ppm2fr  = 1.0e-6_wp             !< Conversion factor ppm to
+                                                                         !< fraction
+    REAL(wp), PARAMETER              ::  t_std   = 273.15_wp             !< standard pressure (Pa)
+    REAL(wp), PARAMETER              ::  vmolcm  = 22.414e3_wp           !< Mole volume (22.414 l)
+                                                                         !< in cm^3
+    REAL(wp), PARAMETER              ::  xna     = 6.022e23_wp           !< Avogadro number
+                                                                         !< (molecules/mol)
+
+    REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
+
+    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact
+    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact_i    !< conversion factor between
+                                                                              !< molecules cm^{-3} and ppm
+    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_qvap
+    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_temp
+
     REAL(kind=wp), DIMENSION(nzb+1:nzt,nspec)                :: tmp_conc
-    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_temp
-    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_qvap
     REAL(kind=wp), DIMENSION(nzb+1:nzt,nphot)                :: tmp_phot
-    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact
-    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact_i    !< conversion factor between 
-                                                                              !<    molecules cm^{-3} and ppm
-
-    INTEGER,DIMENSION(nzb+1:nzt)                            :: nacc          !< Number of accepted steps
-    INTEGER,DIMENSION(nzb+1:nzt)                            :: nrej          !< Number of rejected steps
-
-    REAL(wp)                         ::  conv                                !< conversion factor
-    REAL(wp), PARAMETER              ::  ppm2fr  = 1.0e-6_wp                 !< Conversion factor ppm to fraction
-    REAL(wp), PARAMETER              ::  fr2ppm  = 1.0e6_wp                  !< Conversion factor fraction to ppm
-!    REAL(wp), PARAMETER              ::  xm_air  = 28.96_wp                  !< Mole mass of dry air
-!    REAL(wp), PARAMETER              ::  xm_h2o  = 18.01528_wp               !< Mole mass of water vapor
-    REAL(wp), PARAMETER              ::  t_std   = 273.15_wp                 !< standard pressure (Pa)
-    REAL(wp), PARAMETER              ::  p_std   = 101325.0_wp               !< standard pressure (Pa)
-    REAL(wp), PARAMETER              ::  vmolcm  = 22.414e3_wp               !< Mole volume (22.414 l) in cm^3
-    REAL(wp), PARAMETER              ::  xna     = 6.022e23_wp               !< Avogadro number (molecules/mol)
-
-    REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
-
-    REAL(kind=wp)  :: dt_chem                                             
+
 !
 !-- Set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
@@ -2319 +2327 @@
        tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt)
 !
-!--    convert ppm to molecules/cm**3
-!--    tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp * 
-!--               hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp )  
+!--    Convert ppm to molecules/cm**3
+!--    tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp *
+!--               hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp )
        conv = ppm2fr * xna / vmolcm
        tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
@@ -2328 +2336 @@
        IF ( humidity )  THEN
           IF ( bulk_cloud_model )  THEN
-             tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) *  &
+             tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) *                              &
                              xm_air/xm_h2o * fr2ppm * tmp_fact(:)
           ELSE
@@ -2334 +2342 @@
           ENDIF
        ELSE
-          tmp_qvap(:) = 0.01 * xm_air/xm_h2o * fr2ppm * tmp_fact(:)          !< Constant value for q if water vapor is not computed
+          tmp_qvap(:) = 0.01 * xm_air/xm_h2o * fr2ppm * tmp_fact(:)   !< Constant value for q if
+                                                                      !< water vapor is not computed
        ENDIF
 
        DO  lsp = 1,nspec
-          tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:) 
+          tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:)
        ENDDO
 
        DO lph = 1,nphot
-          tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)               
+          tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)
        ENDDO
 !
@@ -2354 +2363 @@
        cs_time_step = dt_chem
 
-       IF(maxval(rcntrl) > 0.0)  THEN    ! Only if rcntrl is set
+       IF ( MAXVAL( rcntrl ) > 0.0 )  THEN    ! Only if rcntrl is set
           IF( time_since_reference_point <= 2*dt_3d)  THEN
              rcntrl_local = 0
@@ -2364 +2373 @@
        END IF
 
-       CALL chem_gasphase_integrate ( dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot, &
+       CALL chem_gasphase_integrate ( dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot,   &
             icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus )
 
@@ -2378 +2387 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine defining parin for &chemistry_parameters for chemistry model
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_parin
 
@@ -2392 +2401 @@
 
 
-    CHARACTER (LEN=80) ::  line                        !< dummy string that contains the current line of the parameter file
-
-    REAL(wp), DIMENSION(nmaxfixsteps) ::   my_steps    !< List of fixed timesteps   my_step(1) = 0.0 automatic stepping
-    INTEGER(iwp) ::  i                                 !< 
-    INTEGER(iwp) ::  max_pr_cs_tmp                     !< 
-
-
-    NAMELIST /chemistry_parameters/  bc_cs_b,                          &
-         bc_cs_l,                          &
-         bc_cs_n,                          &
-         bc_cs_r,                          &
-         bc_cs_s,                          &
-         bc_cs_t,                          &
-         call_chem_at_all_substeps,        &
-         chem_debug0,                      &
-         chem_debug1,                      &
-         chem_debug2,                      &
-         chem_gasphase_on,                 &
-         chem_mechanism,                   &          
-         cs_heights,                       &
-         cs_name,                          &
-         cs_profile,                       &
-         cs_surface,                       &
-         cs_surface_initial_change,        &
-         cs_vertical_gradient_level,       &
-         daytype_mdh,                      &
-         deposition_dry,                   &
-         emissions_anthropogenic,          &
-         emiss_lod,                        &
-         emiss_factor_main,                &
-         emiss_factor_side,                &
-         emiss_read_legacy_mode,           &
-         icntrl,                           &
-         main_street_id,                   &
-         max_street_id,                    &
-         mode_emis,                        &
-         my_steps,                         &
-         nesting_chem,                     &
-         nesting_offline_chem,             &
-         rcntrl,                           &
-         side_street_id,                   &
-         photolysis_scheme,                &
-         wall_csflux,                      &
-         cs_vertical_gradient,             &
-         top_csflux,                       & 
-         surface_csflux,                   &
-         surface_csflux_name,              &
+    CHARACTER (LEN=80) ::  line                        !< dummy string that contains the current
+                                                       !< line of the parameter file
+
+    INTEGER(iwp) ::  i                                 !<
+    INTEGER(iwp) ::  max_pr_cs_tmp                     !<
+
+    REAL(wp), DIMENSION(nmaxfixsteps) ::   my_steps    !< List of fixed timesteps   my_step(1) = 0.0
+                                                       !< automatic stepping
+
+
+    NAMELIST /chemistry_parameters/                                                                &
+         bc_cs_b,                                                                                  &
+         bc_cs_l,                                                                                  &
+         bc_cs_n,                                                                                  &
+         bc_cs_r,                                                                                  &
+         bc_cs_s,                                                                                  &
+         bc_cs_t,                                                                                  &
+         call_chem_at_all_substeps,                                                                &
+         chem_debug0,                                                                              &
+         chem_debug1,                                                                              &
+         chem_debug2,                                                                              &
+         chem_gasphase_on,                                                                         &
+         chem_mechanism,                                                                           &
+         cs_heights,                                                                               &
+         cs_name,                                                                                  &
+         cs_profile,                                                                               &
+         cs_surface,                                                                               &
+         cs_surface_initial_change,                                                                &
+         cs_vertical_gradient_level,                                                               &
+         daytype_mdh,                                                                              &
+         deposition_dry,                                                                           &
+         emissions_anthropogenic,                                                                  &
+         emiss_lod,                                                                                &
+         emiss_factor_main,                                                                        &
+         emiss_factor_side,                                                                        &
+         emiss_read_legacy_mode,                                                                   &
+         icntrl,                                                                                   &
+         main_street_id,                                                                           &
+         max_street_id,                                                                            &
+         mode_emis,                                                                                &
+         my_steps,                                                                                 &
+         nesting_chem,                                                                             &
+         nesting_offline_chem,                                                                     &
+         rcntrl,                                                                                   &
+         side_street_id,                                                                           &
+         photolysis_scheme,                                                                        &
+         wall_csflux,                                                                              &
+         cs_vertical_gradient,                                                                     &
+         top_csflux,                                                                               &
+         surface_csflux,                                                                           &
+         surface_csflux_name,                                                                      &
          time_fac_type
 !
-!-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
+!-- Analogically to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
 !-- so this way we could prescribe a specific flux value for each species
     !>  chemistry_parameters for initial profiles
     !>  cs_names = 'O3', 'NO2', 'NO', ...   to set initial profiles)
     !>  cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3)
-    !>  cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.) 
+    !>  cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.)
     !>  cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, .....  (chem spcs conc at height lvls chem_heights(1,:)) etc.
     !>  If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
     !>  then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
 !
-!--   Read chem namelist   
-
+!-- Read chem namelist
     CHARACTER(LEN=8)    :: solver_type
 
@@ -2462 +2474 @@
     rtol = 0.01_wp
 !
-!--   Try to find chemistry package
+!-- Try to find chemistry package
     REWIND ( 11 )
     line = ' '
@@ -2470 +2482 @@
     BACKSPACE ( 11 )
 !
-!--   Read chemistry namelist
-    READ ( 11, chemistry_parameters, ERR = 10, END = 20 )      
-!
-!--   Enable chemistry model
-    air_chemistry = .TRUE.                    
+!-- Read chemistry namelist
+    READ ( 11, chemistry_parameters, ERR = 10, END = 20 )
+!
+!-- Enable chemistry model
+    air_chemistry = .TRUE.
     GOTO 20
 
@@ -2486 +2498 @@
 
 !
-!-- synchronize emiss_lod and mod_emis only if emissions_anthropogenic
-!-- is activated in the namelist.  Otherwise their values are "don't care"
+!-- Synchronize emiss_lod and mod_emis only if emissions_anthropogenic is activated in the namelist.
+!-- Otherwise their values are "don't care"
     IF ( emissions_anthropogenic )  THEN
 
 !
-!--    check for emission mode for chem species
-
+!--    Check for emission mode for chem species
        IF ( emiss_lod < 0 )  THEN   !- if LOD not defined in namelist
-          IF ( ( mode_emis /= 'PARAMETERIZED'  )    .AND.      &
-               ( mode_emis /= 'DEFAULT'        )    .AND.      &
+          IF ( ( mode_emis /= 'PARAMETERIZED'  )    .AND.                                          &
+               ( mode_emis /= 'DEFAULT'        )    .AND.                                          &
                ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
              message_string = 'Incorrect mode_emiss  option select. Please check spelling'
@@ -2501 +2512 @@
           ENDIF
        ELSE
-          IF ( ( emiss_lod /= 0 )    .AND.         &
-               ( emiss_lod /= 1 )    .AND.         &
+          IF ( ( emiss_lod /= 0 )    .AND.                                                         &
+               ( emiss_lod /= 1 )    .AND.                                                         &
                ( emiss_lod /= 2 ) )  THEN
              message_string = 'Invalid value for emiss_lod (0, 1, or 2)'
@@ -2510 +2521 @@
 
 !
-! for reference (ecc)
+! For reference (ecc)
 !    IF ( (mode_emis /= 'PARAMETERIZED')  .AND. ( mode_emis /= 'DEFAULT' ) .AND. ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
 !       message_string = 'Incorrect mode_emiss  option select. Please check spelling'
@@ -2517 +2528 @@
 
 !
-!-- conflict resolution for emiss_lod and mode_emis
+!-- Conflict resolution for emiss_lod and mode_emis
 !-- 1) if emiss_lod is defined, have mode_emis assume same setting as emiss_lod
 !-- 2) if emiss_lod it not defined, have emiss_lod assuem same setting as mode_emis
-!-- this check is in place to retain backward compatibility with salsa until the
-!-- code is migrated completed to emiss_lod
-!-- note that 
-
+!-- this check is in place to retain backward compatibility with salsa until the code is migrated
+!-- completed to emiss_lod
+!-- note that
        IF  ( emiss_lod >= 0 ) THEN
 
@@ -2534 +2544 @@
                 mode_emis = 'PRE-PROCESSED'
           END SELECT
-       
-          message_string = 'Synchronizing mode_emis to defined emiss_lod'               //  &
-                           CHAR(10)  //  '                    '                         //  &
-                           'NOTE - mode_emis will be depreciated in future releases'    //  &
-                           CHAR(10)  //  '                    '                         //  &
+
+          message_string = 'Synchronizing mode_emis to defined emiss_lod'               //         &
+                           CHAR( 10 )  //  '                    '                       //         &
+                           'NOTE - mode_emis will be depreciated in future releases'    //         &
+                           CHAR( 10 )  //  '                    '                       //         &
                            'please use emiss_lod to define emission mode'
- 
+
           CALL message ( 'parin_chem', 'CM0463', 0, 0, 0, 6, 0 )
 
@@ -2554 +2564 @@
           END SELECT
 
-          message_string = 'emiss_lod undefined.  Using existing mod_emiss setting'     //  &
-                           CHAR(10)  //  '                    '                         //  &
-                           'NOTE - mode_emis will be depreciated in future releases'    //  &
-                           CHAR(10)  //  '                    '                         //  &
+          message_string = 'emiss_lod undefined.  Using existing mod_emiss setting'     //         &
+                           CHAR( 10 )  //  '                    '                       //         &
+                           'NOTE - mode_emis will be depreciated in future releases'    //         &
+                           CHAR( 10 )  //  '                    '                       //         &
                            'please use emiss_lod to define emission mode'
 
@@ -2576 +2586 @@
        IF ( emiss_read_legacy_mode )  THEN       !< notify legacy read mode
 
-          message_string = 'Legacy emission read mode activated'            // &
-                           CHAR(10)  //  '                    '             // &
-                           'All emissions data will be loaded '             // &
+          message_string = 'Legacy emission read mode activated'            //                     &
+                           CHAR( 10 )  //  '                    '           //                     &
+                           'All emissions data will be loaded '             //                     &
                            'prior to start of simulation'
           CALL message ( 'parin_chem', 'CM0465', 0, 0, 0, 6, 0 )
 
-       ELSE                                     !< if new read mode selected 
+       ELSE                                     !< if new read mode selected
 
           IF ( emiss_lod < 2 )  THEN            !< check LOD compatibility
 
-             message_string = 'New emission read mode '                    // &
-                              'currently unavailable for LODs 0 and 1.'    // &
-                              CHAR(10)  //  '                    '         // &
+             message_string = 'New emission read mode '                    //                      &
+                              'currently unavailable for LODs 0 and 1.'    //                      &
+                              CHAR( 10 )  //  '                    '       //                      &
                               'Reverting to legacy emission read mode'
              CALL message ( 'parin_chem', 'CM0466', 0, 0, 0, 6, 0 )
@@ -2596 +2606 @@
           ELSE                                  !< notify new read mode
 
-             message_string = 'New emission read mode activated'           // &
-                              CHAR(10)  //  '                    '         // &
-                              'LOD 2 emissions will be updated on-demand ' // &           
+             message_string = 'New emission read mode activated'           //                      &
+                              CHAR( 10 )  //  '                    '       //                      &
+                              'LOD 2 emissions will be updated on-demand ' //                      &
                               'according to indicated timestamps'
              CALL message ( 'parin_chem', 'CM0467', 0, 0, 0, 6, 0 )
@@ -2605 +2615 @@
 
        ENDIF ! if emiss_read_legacy_mode
- 
-                           
+
+
     ENDIF  ! if emissions_anthropengic
 
 
-    t_steps = my_steps          
-!
-!--    Determine the number of user-defined profiles and append them to the
-!--    standard data output (data_output_pr)
-    max_pr_cs_tmp = 0 
+    t_steps = my_steps
+!
+!-- Determine the number of user-defined profiles and append them to the standard data output
+!>> (data_output_pr)
+    max_pr_cs_tmp = 0
     i = 1
 
@@ -2628 +2638 @@
        max_pr_cs = max_pr_cs_tmp
     ENDIF
-
-    !      Set Solver Type
+!
+!-- Set Solver Type
     IF(icntrl(3) == 0)  THEN
        solver_type = 'rodas3'           !Default
@@ -2676 +2686 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for all grid points
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
     SUBROUTINE chem_actions( location )
 
@@ -2691 +2701 @@
 !
 !--       Chemical reactions and deposition
-          IF ( chem_gasphase_on ) THEN
+          IF ( chem_gasphase_on )  THEN
 !
 !--          If required, calculate photolysis frequencies -
@@ -2709 +2719 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for grid points i,j
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 
     SUBROUTINE chem_actions_ij( i, j, location )
 
+    CHARACTER (LEN=*), INTENT(IN) ::  location  !< call location string
+
+    INTEGER(iwp)  ::  dummy                     !< call location string
 
     INTEGER(iwp),      INTENT(IN) ::  i         !< grid index in x-direction
     INTEGER(iwp),      INTENT(IN) ::  j         !< grid index in y-direction
-    CHARACTER (LEN=*), INTENT(IN) ::  location  !< call location string
-    INTEGER(iwp)  ::  dummy  !< call location string
 
     IF ( air_chemistry    )   dummy = i + j
@@ -2736 +2747 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for all grid points
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
     SUBROUTINE chem_non_advective_processes()
 
@@ -2747 +2758 @@
       INTEGER(iwp) ::  j  !<
 
-      !
-      !-- Calculation of chemical reactions and deposition.
-
-
-      IF ( intermediate_timestep_count == 1 .OR. call_chem_at_all_substeps )  THEN
-
-         IF ( chem_gasphase_on ) THEN
+!
+!--   Calculation of chemical reactions and deposition.
+      IF ( intermediate_timestep_count == 1  .OR.  call_chem_at_all_substeps )  THEN
+
+         IF ( chem_gasphase_on )  THEN
             CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
             !$OMP PARALLEL PRIVATE (i,j)
@@ -2783 +2792 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for grid points i,j
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_non_advective_processes_ij( i, j )
 
@@ -2796 +2805 @@
 !
 !-- Calculation of chemical reactions and deposition.
-
-
-   IF ( intermediate_timestep_count == 1 .OR. call_chem_at_all_substeps )  THEN
-
-      IF ( chem_gasphase_on ) THEN
+   IF ( intermediate_timestep_count == 1  .OR.  call_chem_at_all_substeps )  THEN
+
+      IF ( chem_gasphase_on )  THEN
          CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
          CALL chem_integrate( i, j )
@@ -2817 +2824 @@
 
  END SUBROUTINE chem_non_advective_processes_ij
- 
-!------------------------------------------------------------------------------!
+
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> routine for exchange horiz of chemical quantities  
-!------------------------------------------------------------------------------! 
+!> routine for exchange horiz of chemical quantities
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_exchange_horiz_bounds
-  
-    USE exchange_horiz_mod,                                                    &
+
+    USE exchange_horiz_mod,                                                                        &
         ONLY:  exchange_horiz
 
    INTEGER(iwp) ::  lsp       !<
- 
-!
-!--    Loop over chemical species       
+
+!
+!--    Loop over chemical species
        CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' )
        DO  lsp = 1, nvar
@@ -2845 +2852 @@
  END SUBROUTINE chem_exchange_horiz_bounds
 
- 
-!------------------------------------------------------------------------------!
+
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Subroutine calculating prognostic equations for chemical species 
-!> (vector-optimized).
-!> Routine is called separately for each chemical species over a loop from
-!> prognostic_equations.
-!------------------------------------------------------------------------------!
+!> Subroutine calculating prognostic equations for chemical species (vector-optimized).
+!> Routine is called separately for each chemical species over a loop from prognostic_equations.
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_prognostic_equations()
 
@@ -2945 +2950 @@
                    DO  k = nzb+1, nzt
                       chem_species(ilsp)%tconc_m(k,j,i) = - 9.5625_wp * tend(k,j,i)                &
-                           + 5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
+                                                     + 5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
                    ENDDO
                 ENDDO
@@ -2959 +2964 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Subroutine calculating prognostic equations for chemical species 
-!> (cache-optimized).
-!> Routine is called separately for each chemical species over a loop from
-!> prognostic_equations.
-!------------------------------------------------------------------------------!
+!> Subroutine calculating prognostic equations for chemical species (cache-optimized).
+!> Routine is called separately for each chemical species over a loop from prognostic_equations.
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_prognostic_equations_ij( i, j, i_omp_start, tn )
 
 
     INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn
+
     INTEGER(iwp) :: ilsp
 !
@@ -3009 +3013 @@
 !
 !--    Diffusion terms (the last three arguments are zero)
-
        CALL diffusion_s( i, j, chem_species(ilsp)%conc,                                            &
             surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:),                              &
@@ -3033 +3036 @@
                )
 
-          IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp )  THEN                
-             chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)    !FKS6                
+          IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp )  THEN
+             chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)    !FKS6
           ENDIF
        ENDDO
@@ -3048 +3051 @@
              DO  k = nzb+1, nzt
                 chem_species(ilsp)%tconc_m(k,j,i) = -9.5625_wp * tend(k,j,i) +                     &
-                     5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
+                                                     5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
              ENDDO
           ENDIF
@@ -3058 +3061 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Read module-specific local restart data arrays (Fortran binary format).
-!------------------------------------------------------------------------------!
- SUBROUTINE chem_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,             &
-                                nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc,   &
-                                nys_on_file, tmp_3d, found )
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE chem_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync,   &
+                                nyn_on_file, nysf, nysc, nys_on_file, tmp_3d, found )
 
     USE control_parameters
 
-
+    INTEGER(iwp) ::  k               !<
     INTEGER(iwp) ::  lsp             !<
-    INTEGER(iwp) ::  k               !< 
-    INTEGER(iwp) ::  nxlc            !< 
-    INTEGER(iwp) ::  nxlf            !< 
-    INTEGER(iwp) ::  nxl_on_file     !<   
-    INTEGER(iwp) ::  nxrc            !< 
-    INTEGER(iwp) ::  nxrf            !< 
-    INTEGER(iwp) ::  nxr_on_file     !<   
-    INTEGER(iwp) ::  nync            !< 
-    INTEGER(iwp) ::  nynf            !< 
-    INTEGER(iwp) ::  nyn_on_file     !<   
-    INTEGER(iwp) ::  nysc            !< 
-    INTEGER(iwp) ::  nysf            !< 
-    INTEGER(iwp) ::  nys_on_file     !<   
-
-    LOGICAL, INTENT(OUT) :: found 
-
-    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
-
-
-    found = .FALSE.  
+    INTEGER(iwp) ::  nxlc            !<
+    INTEGER(iwp) ::  nxlf            !<
+    INTEGER(iwp) ::  nxl_on_file     !<
+    INTEGER(iwp) ::  nxrc            !<
+    INTEGER(iwp) ::  nxrf            !<
+    INTEGER(iwp) ::  nxr_on_file     !<
+    INTEGER(iwp) ::  nync            !<
+    INTEGER(iwp) ::  nynf            !<
+    INTEGER(iwp) ::  nyn_on_file     !<
+    INTEGER(iwp) ::  nysc            !<
+    INTEGER(iwp) ::  nysf            !<
+    INTEGER(iwp) ::  nys_on_file     !<
+
+    LOGICAL, INTENT(OUT) :: found
+
+    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
+
+
+    found = .FALSE.
 
 
@@ -3123 +3125 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Read module-specific local restart data arrays (Fortran binary format).
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_rrd_local_mpi
 
@@ -3154 +3156 @@
 
 
-!-------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 !> Description:
 !> Calculation of horizontally averaged profiles
@@ -3160 +3162 @@
 !> but at least for the region "total domain" (sr=0).
 !> quantities.
-!-------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_statistics( mode, sr, tn )
 
@@ -3169 +3171 @@
 
 
-    CHARACTER (LEN=*) ::  mode   !< 
+    CHARACTER (LEN=*) ::  mode   !<
 
     INTEGER(iwp) ::  i    !< running index on x-axis
     INTEGER(iwp) ::  j    !< running index on y-axis
     INTEGER(iwp) ::  k    !< vertical index counter
+    INTEGER(iwp) ::  lpr  !< running index chem spcs
     INTEGER(iwp) ::  sr   !< statistical region
-    INTEGER(iwp) ::  tn   !< thread number 
-    INTEGER(iwp) ::  lpr  !< running index chem spcs
+    INTEGER(iwp) ::  tn   !< thread number
 
     IF ( mode == 'profiles' )  THEN
        !
 !
-!--    Sample on how to calculate horizontally averaged profiles of user-
-!--    defined quantities. Each quantity is identified by the index
-!--    "pr_palm+#" where "#" is an integer starting from 1. These
-!--    user-profile-numbers must also be assigned to the respective strings
-!--    given by data_output_pr_cs in routine user_check_data_output_pr.
+!--    Sample on how to calculate horizontally averaged profiles of user-defined quantities. Each
+!--    quantity is identified by the index "pr_palm+#" where "#" is an integer starting from 1.
+!--    These user-profile-numbers must also be assigned to the respective strings given by
+!--    data_output_pr_cs in routine user_check_data_output_pr.
 !--    hom(:,:,:,:) =  dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
 !--                     w*pt*), dim-4 = statistical region.
@@ -3195 +3196 @@
                 DO lpr = 1, cs_pr_count
 
-                   sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) +    &
-                        chem_species(cs_pr_index(lpr))%conc(k,j,i) *       &
-                        rmask(j,i,sr)  *                                   &
-                        MERGE( 1.0_wp, 0.0_wp,                             &
+                   sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) +   &
+                        chem_species(cs_pr_index(lpr))%conc(k,j,i) *                               &
+                        rmask(j,i,sr)  *                                                           &
+                        MERGE( 1.0_wp, 0.0_wp,                                                     &
                         BTEST( wall_flags_total_0(k,j,i), 22 ) )
                 ENDDO
@@ -3211 +3212 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Subroutine for swapping of timelevels for chemical species
-!> called out from subroutine swap_timelevel
-!------------------------------------------------------------------------------!
+!> Subroutine for swapping of timelevels for chemical species called out from subroutine
+!> swap_timelevel
+!--------------------------------------------------------------------------------------------------!
 
 
@@ -3224 +3225 @@
     INTEGER(iwp), INTENT(IN) ::  level
 !
-!-- local variables
+!-- Local variables
     INTEGER(iwp)             ::  lsp
 
 
     IF ( level == 0 )  THEN
-       DO  lsp=1, nvar                                       
+       DO  lsp=1, nvar
           chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_1(:,:,:,lsp)
           chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_2(:,:,:,lsp)
        ENDDO
     ELSE
-       DO  lsp=1, nvar                                        
+       DO  lsp=1, nvar
           chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_2(:,:,:,lsp)
           chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_1(:,:,:,lsp)
@@ -3244 +3245 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine to write restart data for chemistry model
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_wrd_local
 
 
-    INTEGER(iwp) ::  lsp  !< running index for chem spcs. 
+    INTEGER(iwp) ::  lsp  !< running index for chem spcs.
 
     IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
@@ -3279 +3280 @@
 
 
-!-------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Subroutine to calculate the deposition of gases and PMs. For now deposition
-!> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT
-!> considered. The deposition of particles is derived following Zhang et al.,
-!> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010).
-!>     
-!> @TODO: Consider deposition on vertical surfaces    
+!> Subroutine to calculate the deposition of gases and PMs. For now deposition only takes place on
+!> lsm and usm horizontal surfaces. Default surfaces are NOT considered. The deposition of particles
+!> is derived following Zhang et al., 2001, gases are deposited using the DEPAC module
+!> (van Zanten et al., 2010).
+!>
+!> @TODO: Consider deposition on vertical surfaces
 !> @TODO: Consider overlaying horizontal surfaces
-!> @TODO: Consider resolved vegetation    
+!> @TODO: Consider resolved vegetation
 !> @TODO: Check error messages
-!-------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE chem_depo( i, j )
 
-    USE control_parameters,                                                 &   
-         ONLY:  dt_3d, intermediate_timestep_count, latitude,               &
-                time_since_reference_point
-
-    USE arrays_3d,                                                          &
+    USE control_parameters,                                                                        &
+         ONLY:  dt_3d, intermediate_timestep_count, latitude, time_since_reference_point
+
+    USE arrays_3d,                                                                                 &
          ONLY:  dzw, rho_air_zw
 
-    USE palm_date_time_mod,                                                 &
+    USE palm_date_time_mod,                                                                        &
          ONLY:  get_date_time
 
-    USE surface_mod,                                                        &
-         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_lsm_h,        &
-         surf_type, surf_usm_h
-
-    USE radiation_model_mod,                                                &
+    USE surface_mod,                                                                               &
+         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_lsm_h, surf_type, surf_usm_h
+
+    USE radiation_model_mod,                                                                       &
          ONLY:  cos_zenith
 
@@ -3314 +3313 @@
     INTEGER(iwp) ::  day_of_year                   !< current day of the year
     INTEGER(iwp), INTENT(IN) ::  i
+    INTEGER(iwp) ::  i_pspec                       !< index for matching depac gas component
     INTEGER(iwp), INTENT(IN) ::  j
     INTEGER(iwp) ::  k                             !< matching k to surface m at i,j
     INTEGER(iwp) ::  lsp                           !< running index for chem spcs.
-    INTEGER(iwp) ::  luv_palm                      !< index of PALM LSM vegetation_type at current surface element
-    INTEGER(iwp) ::  lup_palm                      !< index of PALM LSM pavement_type at current surface element
-    INTEGER(iwp) ::  luw_palm                      !< index of PALM LSM water_type at current surface element
-    INTEGER(iwp) ::  luu_palm                      !< index of PALM USM walls/roofs at current surface element
-    INTEGER(iwp) ::  lug_palm                      !< index of PALM USM green walls/roofs at current surface element
-    INTEGER(iwp) ::  lud_palm                      !< index of PALM USM windows at current surface element
+    INTEGER(iwp) ::  luv_palm                      !< index of PALM LSM vegetation_type at current
+                                                   !< surface element
+    INTEGER(iwp) ::  lup_palm                      !< index of PALM LSM pavement_type at current
+                                                   !< surface element
+    INTEGER(iwp) ::  luw_palm                      !< index of PALM LSM water_type at current
+                                                   !< surface element
+    INTEGER(iwp) ::  luu_palm                      !< index of PALM USM walls/roofs at current
+                                                   !< surface element
+    INTEGER(iwp) ::  lug_palm                      !< index of PALM USM green walls/roofs at current
+                                                   !< surface element
+    INTEGER(iwp) ::  lud_palm                      !< index of PALM USM windows at current surface
+                                                   !< element
     INTEGER(iwp) ::  luv_dep                       !< matching DEPAC LU to luv_palm
     INTEGER(iwp) ::  lup_dep                       !< matching DEPAC LU to lup_palm
@@ -3330 +3336 @@
     INTEGER(iwp) ::  lud_dep                       !< matching DEPAC LU to lud_palm
     INTEGER(iwp) ::  m                             !< index for horizontal surfaces
-
     INTEGER(iwp) ::  pspec                         !< running index
-    INTEGER(iwp) ::  i_pspec                       !< index for matching depac gas component
-!
-!-- Vegetation                                               !< Assign PALM classes to DEPAC land use classes
-    INTEGER(iwp) ::  ind_luv_user = 0                        !<  ERROR as no class given in PALM 
+!
+!-- Vegetation                                               !<  Assign PALM classes to DEPAC land
+                                                             !<  use classes
+    INTEGER(iwp) ::  ind_luv_user = 0                        !<  ERROR as no class given in PALM
     INTEGER(iwp) ::  ind_luv_b_soil = 1                      !<  assigned to ilu_desert
     INTEGER(iwp) ::  ind_luv_mixed_crops = 2                 !<  assigned to ilu_arable
@@ -3352 +3357 @@
     INTEGER(iwp) ::  ind_luv_ev_shrubs = 15                  !<  assigned to ilu_mediterrean_scrub
     INTEGER(iwp) ::  ind_luv_de_shrubs = 16                  !<  assigned to ilu_mediterrean_scrub
-    INTEGER(iwp) ::  ind_luv_mixed_forest = 17               !<  assigned to ilu_coniferous_forest (ave(decid+conif))
+    INTEGER(iwp) ::  ind_luv_mixed_forest = 17               !<  assigned to ilu_coniferous_forest(ave(decid+conif))
     INTEGER(iwp) ::  ind_luv_intrup_forest = 18              !<  assigned to ilu_other (ave(other+decid))
 !
 !-- Water
-    INTEGER(iwp) ::  ind_luw_user = 0                        !<  ERROR as no class given in PALM  
+    INTEGER(iwp) ::  ind_luw_user = 0                        !<  ERROR as no class given in PALM
     INTEGER(iwp) ::  ind_luw_lake = 1                        !<  assigned to ilu_water_inland
     INTEGER(iwp) ::  ind_luw_river = 2                       !<  assigned to ilu_water_inland
     INTEGER(iwp) ::  ind_luw_ocean = 3                       !<  assigned to ilu_water_sea
     INTEGER(iwp) ::  ind_luw_pond = 4                        !<  assigned to ilu_water_inland
-    INTEGER(iwp) ::  ind_luw_fountain = 5                    !<  assigned to ilu_water_inland 
+    INTEGER(iwp) ::  ind_luw_fountain = 5                    !<  assigned to ilu_water_inland
 !
 !-- Pavement
@@ -3369 +3374 @@
     INTEGER(iwp) ::  ind_lup_conc = 3                        !<  assigned to ilu_desert
     INTEGER(iwp) ::  ind_lup_sett = 4                        !<  assigned to ilu_desert
-    INTEGER(iwp) ::  ind_lup_pav_stones = 5                  !<  assigned to ilu_desert 
+    INTEGER(iwp) ::  ind_lup_pav_stones = 5                  !<  assigned to ilu_desert
     INTEGER(iwp) ::  ind_lup_cobblest = 6                    !<  assigned to ilu_desert
     INTEGER(iwp) ::  ind_lup_metal = 7                       !<  assigned to ilu_desert
     INTEGER(iwp) ::  ind_lup_wood = 8                        !<  assigned to ilu_desert
-    INTEGER(iwp) ::  ind_lup_gravel = 9                      !<  assigned to ilu_desert 
+    INTEGER(iwp) ::  ind_lup_gravel = 9                      !<  assigned to ilu_desert
     INTEGER(iwp) ::  ind_lup_f_gravel = 10                   !<  assigned to ilu_desert
     INTEGER(iwp) ::  ind_lup_pebblest = 11                   !<  assigned to ilu_desert
@@ -3385 +3390 @@
     INTEGER(iwp) ::  ind_p_dens = 2     !< index for rhopart in particle_pars
     INTEGER(iwp) ::  ind_p_slip = 3     !< index for slipcor in particle_pars
-
+    INTEGER(iwp) ::  nwet               !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     INTEGER(iwp) ::  part_type          !< index for particle type (PM10 or PM25) in particle_pars
 
-    INTEGER(iwp) ::  nwet               !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
-
-    REAL(wp) ::  dt_chem                !< length of chem time step
-    REAL(wp) ::  dh                     !< vertical grid size 
-    REAL(wp) ::  inv_dh                 !< inverse of vertical grid size 
-    REAL(wp) ::  dt_dh                  !< dt_chem/dh
-
-    REAL(wp) ::  dens              !< density at layer k at i,j  
+
+    REAL(wp) ::  conc_ijk_ugm3     !< concentration at i, j, k in ug/m3
+    REAL(wp) ::  dens              !< density at layer k at i,j
+    REAL(wp) ::  dh                !< vertical grid size
+    REAL(wp) ::  diffusivity       !< diffusivity
+    REAL(wp) ::  dt_chem           !< length of chem time step
+    REAL(wp) ::  dt_dh             !< dt_chem/dh
+    REAL(wp) ::  inv_dh            !< inverse of vertical grid size
+    REAL(wp) ::  lai               !< leaf area index at current surface element
+    REAL(wp) ::  ppm2ugm3          !< conversion factor from ppm to ug/m3
+    REAL(wp) ::  qv_tmp            !< surface mixing ratio at current surface element
     REAL(wp) ::  r_aero_surf       !< aerodynamic resistance (s/m) at current surface element
+    REAL(wp) ::  rb                !< quasi-laminar boundary layer resistance (s/m)
+    REAL(wp) ::  rc_tot            !< total canopy resistance (s/m)
+    REAL(wp) ::  rh_surf           !< relative humidity at current surface element
+    REAL(wp) ::  rs                !< Sedimentaion resistance (s/m)
+    REAL(wp) ::  sai               !< surface area index at current surface element assumed to be
+                                   !< lai + 1
+    REAL(wp) ::  slinnfac
+    REAL(wp) ::  solar_rad         !< solar radiation, direct and diffuse, at current surface
+                                   !< element
+    REAL(wp) ::  temp_tmp          !< temperatur at i,j,k
+    REAL(wp) ::  ts                !< surface temperatur in degrees celsius
     REAL(wp) ::  ustar_surf        !< ustar at current surface element
-    REAL(wp) ::  z0h_surf          !< roughness length for heat at current surface element 
-    REAL(wp) ::  solar_rad         !< solar radiation, direct and diffuse, at current surface element 
-    REAL(wp) ::  ppm2ugm3          !< conversion factor from ppm to ug/m3
-    REAL(wp) ::  rh_surf           !< relative humidity at current surface element
-    REAL(wp) ::  lai               !< leaf area index at current surface element
-    REAL(wp) ::  sai               !< surface area index at current surface element assumed to be lai + 1
-
-    REAL(wp) ::  slinnfac        
+    REAL(wp) ::  vd_lu             !< deposition velocity (m/s)
     REAL(wp) ::  visc              !< Viscosity
     REAL(wp) ::  vs                !< Sedimentation velocity
-    REAL(wp) ::  vd_lu             !< deposition velocity (m/s)
-    REAL(wp) ::  rs                !< Sedimentaion resistance (s/m)
-    REAL(wp) ::  rb                !< quasi-laminar boundary layer resistance (s/m)
-    REAL(wp) ::  rc_tot            !< total canopy resistance (s/m)
-
-    REAL(wp) ::  conc_ijk_ugm3     !< concentration at i, j, k in ug/m3
-    REAL(wp) ::  diffusivity       !< diffusivity
-
-
-    REAL(wp), DIMENSION(nspec) ::  bud_luv      !< budget for LSM vegetation type at current surface element
-    REAL(wp), DIMENSION(nspec) ::  bud_lup      !< budget for LSM pavement type at current surface element
-    REAL(wp), DIMENSION(nspec) ::  bud_luw      !< budget for LSM water type at current surface element
-    REAL(wp), DIMENSION(nspec) ::  bud_luu      !< budget for USM walls/roofs at current surface element
-    REAL(wp), DIMENSION(nspec) ::  bud_lug      !< budget for USM green surfaces at current surface element
+    REAL(wp) ::  z0h_surf          !< roughness length for heat at current surface element
+
+    REAL(wp), DIMENSION(nspec) ::  bud          !< overall budget at current surface element
     REAL(wp), DIMENSION(nspec) ::  bud_lud      !< budget for USM windows at current surface element
-    REAL(wp), DIMENSION(nspec) ::  bud          !< overall budget at current surface element
+    REAL(wp), DIMENSION(nspec) ::  bud_lug      !< budget for USM green surfaces at current surface
+                                                !< element
+    REAL(wp), DIMENSION(nspec) ::  bud_lup      !< budget for LSM pavement type at current surface
+                                                !< element
+    REAL(wp), DIMENSION(nspec) ::  bud_luu      !< budget for USM walls/roofs at current surface
+                                                !< element
+    REAL(wp), DIMENSION(nspec) ::  bud_luv      !< budget for LSM vegetation type at current surface
+                                                !< element
+    REAL(wp), DIMENSION(nspec) ::  bud_luw      !< budget for LSM water type at current surface
+                                                !< element
+    REAL(wp), DIMENSION(nspec) ::  ccomp_tot    !< total compensation point (ug/m3), for now kept to
+                                                !< zero for all species!
     REAL(wp), DIMENSION(nspec) ::  conc_ijk     !< concentration at i,j,k
-    REAL(wp), DIMENSION(nspec) ::  ccomp_tot    !< total compensation point (ug/m3), for now kept to zero for all species!
-                                                  
-
-    REAL(wp) ::  temp_tmp       !< temperatur at i,j,k
-    REAL(wp) ::  ts             !< surface temperatur in degrees celsius
-    REAL(wp) ::  qv_tmp         !< surface mixing ratio at current surface element
-!
-!-- Particle parameters (PM10 (1), PM25 (2))
-!-- partsize (diameter in m), rhopart (density in kg/m3), slipcor
-!-- (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3)
-    REAL(wp), DIMENSION(1:3,1:2), PARAMETER ::  particle_pars = RESHAPE( (/ &
-         8.0e-6_wp, 1.14e3_wp, 1.016_wp, &  !<  1
-         0.7e-6_wp, 1.14e3_wp, 1.082_wp &   !<  2
-         /), (/ 3, 2 /) )
-
+
+
+!
+!-- Particle parameters (PM10 (1), PM25 (2)) partsize (diameter in m), rhopart (density in kg/m3),
+!-- slipcor (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3)
     LOGICAL ::  match_lsm     !< flag indicating natural-type surface
     LOGICAL ::  match_usm     !< flag indicating urban-type surface
+
+    REAL(wp), DIMENSION(1:3,1:2), PARAMETER ::  particle_pars = RESHAPE( (/                        &
+                                                         8.0e-6_wp, 1.14e3_wp, 1.016_wp, &   !<  1
+                                                         0.7e-6_wp, 1.14e3_wp, 1.082_wp  &   !<  2
+                                                         /), (/ 3, 2 /) )
 !
 !-- List of names of possible tracers
-    CHARACTER(LEN=*), PARAMETER ::  pspecnames(nposp) = (/ &
+    CHARACTER(LEN=*), PARAMETER ::  pspecnames(nposp) = (/                                         &
          'NO2           ', &    !< NO2
          'NO            ', &    !< NO
@@ -3515 +3522 @@
          'O3_biascorr   ' /)    !< o3_biascorr
 !
-!-- tracer mole mass:
-    REAL(wp), PARAMETER ::  specmolm(nposp) = (/ &
+!-- Tracer mole mass:
+    REAL(wp), PARAMETER ::  specmolm(nposp) = (/                                                   &
          xm_O * 2 + xm_N, &                         !< NO2
          xm_O + xm_N, &                             !< NO
@@ -3600 +3607 @@
 !
 !--For LSM surfaces
-
     IF ( match_lsm )  THEN
 !
@@ -3612 +3618 @@
        r_aero_surf = surf_lsm_h%r_a(m)
        solar_rad   = surf_lsm_h%rad_sw_dir(m) + surf_lsm_h%rad_sw_dif(m)
-       
+
        lai = surf_lsm_h%lai(m)
        sai = lai + 1
@@ -3628 +3634 @@
        luw_palm = 0
        luw_dep = 0
-! 
+!
 !--    Initialize budgets
        bud_luv = 0.0_wp
@@ -3675 +3681 @@
              luv_dep = 4
           ELSEIF ( luv_palm == ind_luv_intrup_forest )  THEN
-             luv_dep = 8      
+             luv_dep = 8
           ENDIF
        ENDIF
@@ -3695 +3701 @@
              lup_dep = 9
           ELSEIF ( lup_palm == ind_lup_cobblest )  THEN
-             lup_dep = 9       
+             lup_dep = 9
           ELSEIF ( lup_palm == ind_lup_metal )  THEN
              lup_dep = 9
           ELSEIF ( lup_palm == ind_lup_wood )  THEN
-             lup_dep = 9    
+             lup_dep = 9
           ELSEIF ( lup_palm == ind_lup_gravel )  THEN
              lup_dep = 9
@@ -3718 +3724 @@
 
        IF ( surf_lsm_h%frac(m,ind_wat_win) > 0 )  THEN
-          luw_palm = surf_lsm_h%water_type(m)     
+          luw_palm = surf_lsm_h%water_type(m)
           IF ( luw_palm == ind_luw_user )  THEN
              message_string = 'No water type defined. Please define water type to enable deposition calculation'
@@ -3731 +3737 @@
              luw_dep = 13
           ELSEIF ( luw_palm == ind_luw_fountain )  THEN
-             luw_dep = 13 
+             luw_dep = 13
           ENDIF
        ENDIF
@@ -3742 +3748 @@
        ENDIF
 !
-!--    Compute length of time step 
+!--    Compute length of time step
        IF ( call_chem_at_all_substeps )  THEN
           dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
@@ -3770 +3776 @@
 !
 !--    Calculate relative humidity from specific humidity for DEPAC
-       qv_tmp = MAX(q(k,j,i),0.0_wp)
+       qv_tmp = MAX( q(k,j,i), 0.0_wp)
        rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
 !
@@ -3781 +3787 @@
 !
 !--       No vegetation on bare soil, desert or ice:
-          IF ( ( luv_palm == ind_luv_b_soil ) .OR. &
-                ( luv_palm == ind_luv_desert ) .OR. &
-                 ( luv_palm == ind_luv_ice ) ) THEN
+          IF ( ( luv_palm == ind_luv_b_soil )  .OR.  ( luv_palm == ind_luv_desert )  .OR.          &
+               ( luv_palm == ind_luv_ice ) ) THEN
 
              lai = 0.0_wp
@@ -3789 +3794 @@
 
           ENDIF
-          
+
           slinnfac = 1.0_wp
 !
@@ -3795 +3800 @@
           DO  lsp = 1, nvar
 !
-!--          Initialize 
+!--          Initialize
              vs = 0.0_wp
              vd_lu = 0.0_wp
@@ -3805 +3810 @@
 !
 !--             Sedimentation velocity
-                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     visc)
-
-                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
-                     vs, &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     nwet, temp_tmp, dens, visc, &
-                     luv_dep,  &
-                     r_aero_surf, ustar_surf )
-
-                bud_luv(lsp) = - conc_ijk(lsp) * &
-                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
+                                                        particle_pars(ind_p_size, part_type),      &
+                                                        particle_pars(ind_p_slip, part_type),      &
+                                                        visc)
+
+                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
+                                            vs,                                                    &
+                                            particle_pars(ind_p_size, part_type),                  &
+                                            particle_pars(ind_p_slip, part_type),                  &
+                                            nwet, temp_tmp, dens, visc,                            &
+                                            luv_dep,                                               &
+                                            r_aero_surf, ustar_surf )
+
+                bud_luv(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
 
 
@@ -3826 +3830 @@
 !
 !--             Sedimentation velocity
-                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     visc)
-
-                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
-                     vs, &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     nwet, temp_tmp, dens, visc, &
-                     luv_dep , &
-                     r_aero_surf, ustar_surf )
-
-                bud_luv(lsp) = - conc_ijk(lsp) * &
-                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
+                                                        particle_pars(ind_p_size, part_type),      &
+                                                        particle_pars(ind_p_slip, part_type),      &
+                                                        visc)
+
+                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
+                                            vs,                                                    &
+                                            particle_pars(ind_p_size, part_type),                  &
+                                            particle_pars(ind_p_slip, part_type),                  &
+                                            nwet, temp_tmp, dens, visc,                            &
+                                            luv_dep ,                                              &
+                                            r_aero_surf, ustar_surf )
+
+                bud_luv(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
 
              ELSE !< GASES
@@ -3866 +3869 @@
 !--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
 !--             Use density at k:
-
                 ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
 !
@@ -3874 +3876 @@
 !
 !--             Diffusivity for DEPAC relevant gases
-!--             Use default value 
+!--             Use default value
                 diffusivity            = 0.11e-4
 !
-!--             overwrite with known coefficients of diffusivity from Massman (1998)
-                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
+!--             Overwrite with known coefficients of diffusivity from Massman (1998)
+                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
                 IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
                 IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
@@ -3887 +3889 @@
 !
 !--             Get quasi-laminar boundary layer resistance rb:
-                CALL get_rb_cell( (luv_dep == ilu_water_sea) .OR. (luv_dep == ilu_water_inland), &
-                     z0h_surf, ustar_surf, diffusivity, &
-                     rb )
+                CALL get_rb_cell( ( luv_dep == ilu_water_sea )  .OR.                               &
+                                  ( luv_dep == ilu_water_inland ), z0h_surf, ustar_surf,           &
+                                  diffusivity, rb )
 !
 !--             Get rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf, solar_rad, cos_zenith, &
-                     rh_surf, lai, sai, nwet, luv_dep, 2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity, &
-                     r_aero_surf , rb )
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
+                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luv_dep,  &
+                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
+                                         diffusivity, r_aero_surf , rb )
 !
 !--             Calculate budget
@@ -3903 +3906 @@
                 ELSE
 
-                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
-
-                   bud_luv(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
-                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                   vd_lu = 1.0_wp / ( r_aero_surf + rb + rc_tot )
+
+                   bud_luv(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                           &
+                                    ( 1.0_wp - EXP( -vd_lu * dt_dh ) ) * dh
                 ENDIF
 
@@ -3935 +3938 @@
 !
 !--             Sedimentation velocity:
-                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     visc)
-
-                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
-                     vs, &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     nwet, temp_tmp, dens, visc, &
-                     lup_dep,  &
-                     r_aero_surf, ustar_surf )
-
-                bud_lup(lsp) = - conc_ijk(lsp) * &
-                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
+                                                        particle_pars(ind_p_size, part_type),      &
+                                                        particle_pars(ind_p_slip, part_type),      &
+                                                        visc)
+
+                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
+                                            vs,                                                    &
+                                            particle_pars(ind_p_size, part_type),                  &
+                                            particle_pars(ind_p_slip, part_type),                  &
+                                            nwet, temp_tmp, dens, visc,                            &
+                                            lup_dep,                                               &
+                                            r_aero_surf, ustar_surf )
+
+                bud_lup(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
 
 
@@ -3956 +3958 @@
 !
 !--             Sedimentation velocity:
-                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     visc)
-
-                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
-                     vs, &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     nwet, temp_tmp, dens, visc, &
-                     lup_dep, &
-                     r_aero_surf, ustar_surf )
-
-                bud_lup(lsp) = - conc_ijk(lsp) * &
-                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
+                                                        particle_pars(ind_p_size, part_type),      &
+                                                        particle_pars(ind_p_slip, part_type),      &
+                                                        visc)
+
+                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
+                                            vs,                                                    &
+                                            particle_pars(ind_p_size, part_type),                  &
+                                            particle_pars(ind_p_slip, part_type),                  &
+                                            nwet, temp_tmp, dens, visc,                            &
+                                            lup_dep,                                               &
+                                            r_aero_surf, ustar_surf )
+
+                bud_lup(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
 
              ELSE  !<GASES
 !
 !--             Read spc_name of current species for gas parameter
-
+                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
+                   i_pspec = 0
+                   DO  pspec = 1, nposp
+                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
+                         i_pspec = pspec
+                      END IF
+                   ENDDO
+
+                ELSE
+!
+!--                For now species not deposited
+                   CYCLE
+                ENDIF
+!
+!--             Factor used for conversion from ppb to ug/m3 :
+!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
+!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
+!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
+!--             thus:
+!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
+!--             Use density at lowest layer:
+                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
+!
+!--             Atmospheric concentration in DEPAC is requested in ug/m3:
+                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
+                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
+!
+!--             Diffusivity for DEPAC relevant gases
+!--             Use default value
+                diffusivity            = 0.11e-4
+!
+!--             Overwrite with known coefficients of diffusivity from Massman (1998)
+                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
+                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
+                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
+                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
+                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
+                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
+                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
+!
+!--             Get quasi-laminar boundary layer resistance rb:
+                CALL get_rb_cell( ( lup_dep == ilu_water_sea )  .OR.                               &
+                                  ( lup_dep == ilu_water_inland ), z0h_surf, ustar_surf,           &
+                                    diffusivity, rb )
+!
+!--             Get rc_tot
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
+                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lup_dep,  &
+                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
+                                         diffusivity, r_aero_surf , rb )
+!
+!--             Calculate budget
+                IF ( rc_tot <= 0.0 )  THEN
+                   bud_lup(lsp) = 0.0_wp
+                ELSE
+                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
+                   bud_lup(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                           &
+                                  ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
+                ENDIF
+
+             ENDIF
+          ENDDO
+       ENDIF
+!
+!--    Water
+       IF ( surf_lsm_h%frac(m,ind_wat_win) > 0 )  THEN
+!
+!--       No vegetation on water:
+          lai = 0.0_wp
+          sai = 0.0_wp
+          slinnfac = 1.0_wp
+!
+!--       Get vd
+          DO  lsp = 1, nvar
+!
+!--          Initialize
+             vs = 0.0_wp
+             vd_lu = 0.0_wp
+             rs = 0.0_wp
+             rb = 0.0_wp
+             rc_tot = 0.0_wp
+             IF ( spc_names(lsp) == 'PM10' )  THEN
+                part_type = 1
+!
+!--             Sedimentation velocity:
+                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
+                                                        particle_pars(ind_p_size, part_type),      &
+                                                        particle_pars(ind_p_slip, part_type),      &
+                                                        visc)
+
+                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
+                                            vs,                                                    &
+                                            particle_pars(ind_p_size, part_type),                  &
+                                            particle_pars(ind_p_slip, part_type),                  &
+                                            nwet, temp_tmp, dens, visc,                            &
+                                            luw_dep,                                               &
+                                            r_aero_surf, ustar_surf )
+
+                bud_luw(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
+
+             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
+                part_type = 2
+!
+!--             Sedimentation velocity:
+                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
+                                                        particle_pars(ind_p_size, part_type), &
+                                                        particle_pars(ind_p_slip, part_type), &
+                                                        visc)
+
+                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
+                                            vs,                                                    &
+                                            particle_pars(ind_p_size, part_type),                  &
+                                            particle_pars(ind_p_slip, part_type),                  &
+                                            nwet, temp_tmp, dens, visc,                            &
+                                            luw_dep,                                               &
+                                            r_aero_surf, ustar_surf )
+
+                bud_luw(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
+
+             ELSE  !<GASES
+!
+!--             Read spc_name of current species for gas PARAMETER
                 IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
                    i_pspec = 0
@@ -3997 +4119 @@
 !--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
 !--             Use density at lowest layer:
-
-                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
+                ppm2ugm3 = (dens/xm_air) * 0.001_wp  !< (mole air)/m3
 !
 !--             Atmospheric concentration in DEPAC is requested in ug/m3:
-                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
-                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
+!--                 ug/m3        ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
+                conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 *  specmolm(i_pspec)  ! in ug/m3
 !
 !--             Diffusivity for DEPAC relevant gases
-!--             Use default value 
+!--             Use default value
                 diffusivity            = 0.11e-4
 !
-!--             overwrite with known coefficients of diffusivity from Massman (1998)
-                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
+!--             Overwrite with known coefficients of diffusivity from Massman (1998)
+                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
                 IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
                 IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
@@ -4018 +4139 @@
 !
 !--             Get quasi-laminar boundary layer resistance rb:
-                CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland),   &
-                     z0h_surf, ustar_surf, diffusivity, rb )
-!
+                CALL get_rb_cell( ( luw_dep == ilu_water_sea )  .OR.                               &
+                                  ( luw_dep == ilu_water_inland ), z0h_surf, ustar_surf,           &
+                                    diffusivity, rb )
+
 !--             Get rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,      &
-                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lup_dep, 2,    &
-                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,              &
-                                         r_aero_surf , rb )
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
+                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luw_dep,  &
+                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
+                                         diffusivity, r_aero_surf , rb )
 !
 !--             Calculate budget
                 IF ( rc_tot <= 0.0 )  THEN
-                   bud_lup(lsp) = 0.0_wp
+
+                   bud_luw(lsp) = 0.0_wp
+
                 ELSE
+
                    vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
-                   bud_lup(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
-                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+
+                   bud_luw(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                           &
+                                  ( 1.0_wp - EXP(-vd_lu * dt_dh ) ) * dh
                 ENDIF
 
@@ -4039 +4165 @@
           ENDDO
        ENDIF
-!
-!--    Water
-       IF ( surf_lsm_h%frac(m,ind_wat_win) > 0 )  THEN
-!
-!--       No vegetation on water:
-          lai = 0.0_wp
-          sai = 0.0_wp
-          slinnfac = 1.0_wp
-!
-!--       Get vd
-          DO  lsp = 1, nvar
-!
-!--          Initialize
-             vs = 0.0_wp
-             vd_lu = 0.0_wp
-             rs = 0.0_wp
-             rb = 0.0_wp
-             rc_tot = 0.0_wp 
-             IF ( spc_names(lsp) == 'PM10' )  THEN
-                part_type = 1
-!
-!--             Sedimentation velocity:
-                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     visc)
-
-                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
-                     vs, &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     nwet, temp_tmp, dens, visc, &
-                     luw_dep, &
-                     r_aero_surf, ustar_surf )
-
-                bud_luw(lsp) = - conc_ijk(lsp) * &
-                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
-
-             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
-                part_type = 2
-!
-!--             Sedimentation velocity:
-                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     visc)
-
-                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
-                     vs, &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     nwet, temp_tmp, dens, visc, &
-                     luw_dep, &
-                     r_aero_surf, ustar_surf )
-
-                bud_luw(lsp) = - conc_ijk(lsp) * &
-                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
-
-             ELSE  !<GASES
-!
-!--             Read spc_name of current species for gas PARAMETER
-
-                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
-                   i_pspec = 0
-                   DO  pspec = 1, nposp
-                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
-                         i_pspec = pspec
-                      END IF
-                   ENDDO
-
-                ELSE
-!
-!--                For now species not deposited
-                   CYCLE
-                ENDIF
-!
-!--             Factor used for conversion from ppb to ug/m3 :
-!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
-!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
-!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
-!--             thus:
-!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
-!--             Use density at lowest layer:
-
-                ppm2ugm3 = (dens/xm_air) * 0.001_wp  !< (mole air)/m3
-!
-!--             Atmospheric concentration in DEPAC is requested in ug/m3:
-!--                 ug/m3        ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
-                conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 *  specmolm(i_pspec)  ! in ug/m3
-!
-!--             Diffusivity for DEPAC relevant gases
-!--             Use default value 
-                diffusivity            = 0.11e-4
-!
-!--             overwrite with known coefficients of diffusivity from Massman (1998)
-                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
-                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
-                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
-                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
-                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
-                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
-                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
-!
-!--             Get quasi-laminar boundary layer resistance rb:
-                CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland),  &
-                     z0h_surf, ustar_surf, diffusivity, rb )
-
-!--             Get rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         & 
-                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luw_dep, 2,    &
-                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
-                                         r_aero_surf , rb )
-!
-!--             Calculate budget 
-                IF ( rc_tot <= 0.0 )  THEN
-
-                   bud_luw(lsp) = 0.0_wp
-
-                ELSE
-
-                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
-
-                   bud_luw(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
-                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
-                ENDIF
-
-             ENDIF
-          ENDDO
-       ENDIF
 
 
@@ -4175 +4172 @@
        DO  lsp = 1, nspec
 
-          bud(lsp) = surf_lsm_h%frac(m,ind_veg_wall) * bud_luv(lsp) + &
-               surf_lsm_h%frac(m,ind_pav_green) * bud_lup(lsp) + &
-               surf_lsm_h%frac(m,ind_wat_win) * bud_luw(lsp)
+          bud(lsp) = surf_lsm_h%frac(m,ind_veg_wall)  * bud_luv(lsp) +                             &
+                     surf_lsm_h%frac(m,ind_pav_green) * bud_lup(lsp) +                             &
+                     surf_lsm_h%frac(m,ind_wat_win)   * bud_luw(lsp)
 !
 !--       Compute new concentration:
           conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
 
-          chem_species(lsp)%conc(k,j,i) = MAX(0.0_wp, conc_ijk(lsp))
+          chem_species(lsp)%conc(k,j,i) = MAX( 0.0_wp, conc_ijk(lsp) )
 
        ENDDO
@@ -4188 +4185 @@
     ENDIF
 !
-!-- For USM surfaces    
-
+!-- For USM surfaces
     IF ( match_usm )  THEN
 !
@@ -4225 +4221 @@
        IF ( surf_usm_h%frac(m,ind_pav_green) > 0 )  THEN
 !
-!--       For green urban surfaces (e.g. green roofs
-!--       assume LU short grass
+!--       For green urban surfaces (e.g. green roofs assume LU short grass
           lug_palm = ind_luv_s_grass
           IF ( lug_palm == ind_luv_user )  THEN
@@ -4266 +4261 @@
              lug_dep = 4
           ELSEIF ( lug_palm == ind_luv_intrup_forest )  THEN
-             lug_dep = 8      
+             lug_dep = 8
           ENDIF
        ENDIF
@@ -4273 +4268 @@
 !
 !--       For walls in USM assume concrete walls/roofs,
-!--       assumed LU class desert as also assumed for
-!--       pavements in LSM
+!--       assumed LU class desert as also assumed for pavements in LSM
           luu_palm = ind_lup_conc
           IF ( luu_palm == ind_lup_user )  THEN
@@ -4290 +4284 @@
              luu_dep = 9
           ELSEIF ( luu_palm == ind_lup_cobblest )  THEN
-             luu_dep = 9       
+             luu_dep = 9
           ELSEIF ( luu_palm == ind_lup_metal )  THEN
              luu_dep = 9
           ELSEIF ( luu_palm == ind_lup_wood )  THEN
-             luu_dep = 9    
+             luu_dep = 9
           ELSEIF ( luu_palm == ind_lup_gravel )  THEN
              luu_dep = 9
@@ -4314 +4308 @@
        IF ( surf_usm_h%frac(m,ind_wat_win) > 0 )  THEN
 !
-!--       For windows in USM assume metal as this is
-!--       as close as we get, assumed LU class desert
-!--       as also assumed for pavements in LSM
-          lud_palm = ind_lup_metal     
+!--       For windows in USM assume metal as this is as close as we get, assumed LU class desert as
+!--       also assumed for pavements in LSM.
+          lud_palm = ind_lup_metal
           IF ( lud_palm == ind_lup_user )  THEN
              message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
@@ -4332 +4325 @@
              lud_dep = 9
           ELSEIF ( lud_palm == ind_lup_cobblest )  THEN
-             lud_dep = 9       
+             lud_dep = 9
           ELSEIF ( lud_palm == ind_lup_metal )  THEN
              lud_dep = 9
           ELSEIF ( lud_palm == ind_lup_wood )  THEN
-             lud_dep = 9    
+             lud_dep = 9
           ELSEIF ( lud_palm == ind_lup_gravel )  THEN
              lud_dep = 9
@@ -4362 +4355 @@
        !ENDIF
 !
-!--    Compute length of time step 
+!--    Compute length of time step
        IF ( call_chem_at_all_substeps )  THEN
           dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
@@ -4371 +4364 @@
        dh = dzw(k)
        inv_dh = 1.0_wp / dh
-       dt_dh = dt_chem/dh
+       dt_dh = dt_chem / dh
 !
 !--    Concentration at i,j,k
@@ -4384 +4377 @@
 !
 !--    Viscosity of air
-       visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
+       visc = 1.496e-6 * temp_tmp**1.5 / ( temp_tmp + 120.0 )
 !
 !--    Air density at k
@@ -4390 +4383 @@
 !
 !--    Calculate relative humidity from specific humidity for DEPAC
-       qv_tmp = MAX(q(k,j,i),0.0_wp)
-       rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
-!
-!--    Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
-!--    for each surface fraction. Then derive overall budget taking into account the surface fractions.
-!
+       qv_tmp = MAX( q(k,j,i), 0.0_wp )
+       rh_surf = relativehumidity_from_specifichumidity( qv_tmp, temp_tmp, hyp(k) )
+!
+!--    Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget for
+!--    each surface fraction. Then derive overall budget taking into account the surface fractions.
+
 !--    Walls/roofs
        IF ( surf_usm_h%frac(m,ind_veg_wall) > 0 )  THEN
@@ -4408 +4401 @@
           DO  lsp = 1, nvar
 !
-!--          Initialize 
+!--          Initialize
              vs = 0.0_wp
              vd_lu = 0.0_wp
@@ -4418 +4411 @@
 !
 !--             Sedimentation velocity
-                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     visc)
-
-                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
-                     vs, &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     nwet, temp_tmp, dens, visc, &
-                     luu_dep,  &
-                     r_aero_surf, ustar_surf )
-
-                bud_luu(lsp) = - conc_ijk(lsp) * &
-                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
+                                                        particle_pars(ind_p_size, part_type),      &
+                                                        particle_pars(ind_p_slip, part_type),      &
+                                                        visc)
+
+                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
+                                            vs,                                                    &
+                                            particle_pars(ind_p_size, part_type),                  &
+                                            particle_pars(ind_p_slip, part_type),                  &
+                                            nwet, temp_tmp, dens, visc,                            &
+                                            luu_dep,                                               &
+                                            r_aero_surf, ustar_surf )
+
+                bud_luu(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
 
              ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
@@ -4438 +4430 @@
 !
 !--             Sedimentation velocity
-                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     visc)
-
-                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
-                     vs, &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     nwet, temp_tmp, dens, visc, &
-                     luu_dep , &
-                     r_aero_surf, ustar_surf )
-
-                bud_luu(lsp) = - conc_ijk(lsp) * &
-                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
+                                                        particle_pars(ind_p_size, part_type),      &
+                                                        particle_pars(ind_p_slip, part_type),      &
+                                                        visc)
+
+                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
+                                            vs,                                                    &
+                                            particle_pars(ind_p_size, part_type),                  &
+                                            particle_pars(ind_p_slip, part_type),                  &
+                                            nwet, temp_tmp, dens, visc,                            &
+                                            luu_dep ,                                              &
+                                            r_aero_surf, ustar_surf )
+
+                bud_luu(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
 
              ELSE  !< GASES
 !
 !--             Read spc_name of current species for gas parameter
-
                 IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
                    i_pspec = 0
@@ -4478 +4468 @@
 !--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
 !--             Use density at k:
-
                 ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
 
@@ -4487 +4476 @@
 !
 !--             Diffusivity for DEPAC relevant gases
-!--             Use default value 
+!--             Use default value
                 diffusivity            = 0.11e-4
 !
-!--             overwrite with known coefficients of diffusivity from Massman (1998)
-                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
+!--             Overwrite with known coefficients of diffusivity from Massman (1998)
+                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
                 IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
                 IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
@@ -4500 +4489 @@
 !
 !--             Get quasi-laminar boundary layer resistance rb:
-                CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland),   &
-                     z0h_surf, ustar_surf, diffusivity, &
-                     rb )
+                CALL get_rb_cell( ( luu_dep == ilu_water_sea )  .OR.                               &
+                                  ( luu_dep == ilu_water_inland ), z0h_surf, ustar_surf,           &
+                                    diffusivity, rb )
 !
 !--             Get rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,          &
-                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luu_dep, 2,     &
-                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,   &
-                                         r_aero_surf, rb )
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
+                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luu_dep,  &
+                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
+                                         diffusivity, r_aero_surf, rb )
 !
 !--             Calculate budget
@@ -4519 +4508 @@
                    vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
 
-                   bud_luu(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
-                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                   bud_luu(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                           &
+                                  ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
                 ENDIF
 
@@ -4532 +4521 @@
 !
 !--       No vegetation on bare soil, desert or ice:
-          IF ( ( lug_palm == ind_luv_b_soil ) .OR. &
-                ( lug_palm == ind_luv_desert ) .OR. &
-                 ( lug_palm == ind_luv_ice ) ) THEN
+          IF ( ( lug_palm == ind_luv_b_soil )  .OR.  ( lug_palm == ind_luv_desert ) .OR.           &
+               ( lug_palm == ind_luv_ice ) ) THEN
 
              lai = 0.0_wp
@@ -4541 +4529 @@
           ENDIF
 
-          
+
           slinnfac = 1.0_wp
 !
@@ -4557 +4545 @@
 !
 !--             Sedimentation velocity
-                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     visc)
-
-                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
-                     vs, &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     nwet, temp_tmp, dens, visc, &
-                     lug_dep,  &
-                     r_aero_surf, ustar_surf )
-
-                bud_lug(lsp) = - conc_ijk(lsp) * &
-                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
+                                                        particle_pars(ind_p_size, part_type),      &
+                                                        particle_pars(ind_p_slip, part_type),      &
+                                                        visc)
+
+                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
+                                            vs,                                                    &
+                                            particle_pars(ind_p_size, part_type),                  &
+                                            particle_pars(ind_p_slip, part_type),                  &
+                                            nwet, temp_tmp, dens, visc,                            &
+                                            lug_dep,                                               &
+                                            r_aero_surf, ustar_surf )
+
+                bud_lug(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
 
              ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
@@ -4577 +4564 @@
 !
 !--             Sedimentation velocity
-                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     visc)
-
-                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
-                     vs, &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     nwet, temp_tmp, dens, visc, &
-                     lug_dep, &
-                     r_aero_surf, ustar_surf )
-
-                bud_lug(lsp) = - conc_ijk(lsp) * &
-                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
+                                                        particle_pars(ind_p_size, part_type),      &
+                                                        particle_pars(ind_p_slip, part_type),      &
+                                                        visc)
+
+                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
+                                            vs,                                                    &
+                                            particle_pars(ind_p_size, part_type),                  &
+                                            particle_pars(ind_p_slip, part_type),                  &
+                                            nwet, temp_tmp, dens, visc,                            &
+                                            lug_dep,                                               &
+                                            r_aero_surf, ustar_surf )
+
+                bud_lug(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
 
              ELSE  !< GASES
 !
 !--             Read spc_name of current species for gas parameter
-
                 IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
                    i_pspec = 0
@@ -4617 +4602 @@
 !--                  c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
 !--             Use density at k:
-
                 ppm2ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
 !
@@ -4625 +4609 @@
 !
 !--             Diffusivity for DEPAC relevant gases
-!--             Use default value 
+!--             Use default value
                 diffusivity            = 0.11e-4
 !
-!--             overwrite with known coefficients of diffusivity from Massman (1998)
-                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
+!--             Overwrite with known coefficients of diffusivity from Massman (1998)
+                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
                 IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
                 IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
@@ -4638 +4622 @@
 !
 !--             Get quasi-laminar boundary layer resistance rb:
-                CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland),    &
-                     z0h_surf, ustar_surf, diffusivity, &
-                     rb )
+                CALL get_rb_cell( ( lug_dep == ilu_water_sea )  .OR.                               &
+                                  ( lug_dep == ilu_water_inland ), z0h_surf, ustar_surf,           &
+                                    diffusivity, rb )
 !
 !--             Get rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,           &
-                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lug_dep, 2,      &
-                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,    &
-                                         r_aero_surf , rb )
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
+                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lug_dep,  &
+                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
+                                         diffusivity, r_aero_surf , rb )
 !
 !--             Calculate budget
@@ -4655 +4639 @@
                 ELSE
 
-                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
-
-                   bud_lug(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
-                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                   vd_lu = 1.0_wp / ( r_aero_surf + rb + rc_tot )
+
+                   bud_lug(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                           &
+                                  ( 1.0_wp - EXP( - vd_lu * dt_dh ) )  * dh
                 ENDIF
 
@@ -4682 +4666 @@
              rs = 0.0_wp
              rb = 0.0_wp
-             rc_tot = 0.0_wp 
+             rc_tot = 0.0_wp
              IF ( spc_names(lsp) == 'PM10' )  THEN
                 part_type = 1
 !
 !--             Sedimentation velocity
-                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     visc)
-
-                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     nwet, temp_tmp, dens, visc,              &
-                     lud_dep, r_aero_surf, ustar_surf )
-
-                bud_lud(lsp) = - conc_ijk(lsp) * &
-                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
+                                                        particle_pars(ind_p_size, part_type),      &
+                                                        particle_pars(ind_p_slip, part_type),      &
+                                                        visc)
+
+                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs,                                         &
+                                            particle_pars(ind_p_size, part_type),                  &
+                                            particle_pars(ind_p_slip, part_type),                  &
+                                            nwet, temp_tmp, dens, visc,                            &
+                                            lud_dep, r_aero_surf, ustar_surf )
+
+                bud_lud(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
 
              ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
@@ -4705 +4688 @@
 !
 !--             Sedimentation velocity
-                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     visc)
-
-                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
-                     particle_pars(ind_p_size, part_type), &
-                     particle_pars(ind_p_slip, part_type), &
-                     nwet, temp_tmp, dens, visc, &
-                     lud_dep, &
-                     r_aero_surf, ustar_surf )
-
-                bud_lud(lsp) = - conc_ijk(lsp) * &
-                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
+                                                        particle_pars(ind_p_size, part_type),      &
+                                                        particle_pars(ind_p_slip, part_type),      &
+                                                        visc)
+
+                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs,                                         &
+                                            particle_pars(ind_p_size, part_type),                  &
+                                            particle_pars(ind_p_slip, part_type),                  &
+                                            nwet, temp_tmp, dens, visc,                            &
+                                            lud_dep,                                               &
+                                            r_aero_surf, ustar_surf )
+
+                bud_lud(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
 
              ELSE  !< GASES
 !
 !--             Read spc_name of current species for gas PARAMETER
-
                 IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
                    i_pspec = 0
@@ -4752 +4733 @@
 !
 !--             Diffusivity for DEPAC relevant gases
-!--             Use default value 
+!--             Use default value
                 diffusivity = 0.11e-4
 !
-!--             overwrite with known coefficients of diffusivity from Massman (1998)
-                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
+!--             Overwrite with known coefficients of diffusivity from Massman (1998)
+                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
                 IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
                 IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
@@ -4765 +4746 @@
 !
 !--             Get quasi-laminar boundary layer resistance rb:
-                CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland),   &
-                     z0h_surf, ustar_surf, diffusivity, rb )
+                CALL get_rb_cell( ( lud_dep == ilu_water_sea )  .OR.                               &
+                                  ( lud_dep == ilu_water_inland ), z0h_surf, ustar_surf,           &
+                                    diffusivity, rb )
 !
 !--             Get rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         &
-                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lud_dep, 2,    &
-                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
-                                         r_aero_surf , rb )
-!
-!--             Calculate budget 
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
+                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lud_dep,  &
+                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
+                                         diffusivity, r_aero_surf , rb )
+!
+!--             Calculate budget
                 IF ( rc_tot <= 0.0 )  THEN
 
@@ -4783 +4765 @@
                    vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
 
-                   bud_lud(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
-                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
+                   bud_lud(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                           &
+                                  ( 1.0_wp - EXP( - vd_lu * dt_dh ) )  * dh
                 ENDIF
 
@@ -4798 +4780 @@
 
 
-          bud(lsp) = surf_usm_h%frac(m,ind_veg_wall) * bud_luu(lsp) + &
-               surf_usm_h%frac(m,ind_pav_green) * bud_lug(lsp) + &
-               surf_usm_h%frac(m,ind_wat_win) * bud_lud(lsp)
+          bud(lsp) = surf_usm_h%frac(m,ind_veg_wall)  * bud_luu(lsp) +                             &
+                     surf_usm_h%frac(m,ind_pav_green) * bud_lug(lsp) +                             &
+                     surf_usm_h%frac(m,ind_wat_win)   * bud_lud(lsp)
 !
 !--       Compute new concentration
@@ -4815 +4797 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
@@ -4821 +4803 @@
 !>
 !> DEPAC:
-!> Code of the DEPAC routine and corresponding subroutines below from the DEPAC
-!> module of the LOTOS-EUROS model (Manders et al., 2017)
+!> Code of the DEPAC routine and corresponding subroutines below from the DEPAC module of the
+!> LOTOS-EUROS model (Manders et al., 2017)
 !>
-!> Original DEPAC routines by RIVM and TNO (2015), for Documentation see
-!> van Zanten et al., 2010.
-!------------------------------------------------------------------------------!
- SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, solar_rad, sinphi,    &
-      rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, conc_ijk_ugm3, diffusivity,  &
-      ra, rb )  
+!> Original DEPAC routines by RIVM and TNO (2015), for Documentation see van Zanten et al., 2010.
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, solar_rad, sinphi, rh, lai, sai,&
+                                nwet, lu, iratns, rc_tot, ccomp_tot, p, conc_ijk_ugm3, diffusivity,&
+                                ra, rb )
 !
 !--   Some of depac arguments are OPTIONAL:
 !--    A. compute Rc_tot without compensation points (ccomp_tot will be zero):
-!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi])
+!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot,&
+!--                    ccomp_tot, [smi])
 !--    B. compute Rc_tot with compensation points (used for LOTOS-EUROS):
-!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
-!--                c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water)
+!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot,&
+!--                    ccomp_tot, [smi], c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water)
 !--
 !--    C. compute effective Rc based on compensation points (used for OPS):
-!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
-!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
-!--               ra, rb, rc_eff)
-!--    X1. Extra (OPTIONAL) output variables: 
-!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
-!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
-!--               ra, rb, rc_eff, &
-!--               gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out)
-!--    X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves): 
-!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
-!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
-!--               ra, rb, rc_eff, &
-!--               gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, &
-!--               calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA)
+!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot,&
+!--                    ccomp_tot, [smi], c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, ra,  &
+!--                    rb, rc_eff)
+!--    X1. Extra (OPTIONAL) output variables:
+!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot,&
+!--                    ccomp_tot, [smi], c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, ra,  &
+!--                    rb, rc_eff, gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out,    &
+!--                    lai_out, sai_out)
+!--    X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves):
+!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot,&
+!--                    ccomp_tot, [smi], c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, ra,  &
+!--                    rb, rc_eff, gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out,    &
+!--                    lai_out, sai_out, calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA)
 
 
@@ -4858 +4839 @@
                                                      !< 'HNO3','NO','NO2','O3','SO2','NH3'
     INTEGER(iwp), INTENT(IN) ::  day_of_year         !< day of year, 1 ... 365 (366)
-    INTEGER(iwp), INTENT(IN) ::  nwet                !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
-    INTEGER(iwp), INTENT(IN) ::  lu                  !< land use type, lu = 1,...,nlu
     INTEGER(iwp), INTENT(IN) ::  iratns              !< index for NH3/SO2 ratio used for SO2:
                                                      !< iratns = 1: low NH3/SO2
                                                      !< iratns = 2: high NH3/SO2
                                                      !< iratns = 3: very low NH3/SO2
-    REAL(wp), INTENT(IN) ::  lat                     !< latitude Northern hemisphere (degrees) (S. hemisphere not possible)
+    INTEGER(iwp), INTENT(IN) ::  lu                  !< land use type, lu = 1,...,nlu
+    INTEGER(iwp), INTENT(IN) ::  nwet                !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
+
+    REAL(wp), INTENT(IN) ::  conc_ijk_ugm3           !< actual atmospheric concentration (ug/m3), in
+                                                     !< DEPAC=Catm
+    REAL(wp), INTENT(IN) ::  diffusivity             !< diffusivity
+    REAL(wp), INTENT(IN) ::  lai                     !< one-sidedleaf area index (-)
+    REAL(wp), INTENT(IN) ::  lat                     !< latitude Northern hemisphere (degrees) (S.
+                                                     !< hemisphere not possible)
+    REAL(wp), INTENT(IN) ::  p                       !< pressure (Pa)
+    REAL(wp), INTENT(IN) ::  ra                      !< aerodynamic resistance (s/m)
+    REAL(wp), INTENT(IN) ::  rb                      !< boundary layer resistance (s/m)
+    REAL(wp), INTENT(IN) ::  rh                      !< relative humidity (%)
+    REAL(wp), INTENT(IN) ::  sai                     !< surface area index (-) (lai + branches and
+                                                     !< stems)
+    REAL(wp), INTENT(IN) ::  sinphi                  !< sin of solar elevation angle
+    REAL(wp), INTENT(IN) ::  solar_rad               !< solar radiation, dirict+diffuse (W/m2)
     REAL(wp), INTENT(IN) ::  t                       !< temperature (C)
     REAL(wp), INTENT(IN) ::  ust                     !< friction velocity (m/s)
-    REAL(wp), INTENT(IN) ::  solar_rad               !< solar radiation, dirict+diffuse (W/m2)
-    REAL(wp), INTENT(IN) ::  sinphi                  !< sin of solar elevation angle
-    REAL(wp), INTENT(IN) ::  rh                      !< relative humidity (%)
-    REAL(wp), INTENT(IN) ::  lai                     !< one-sidedleaf area index (-)
-    REAL(wp), INTENT(IN) ::  sai                     !< surface area index (-) (lai + branches and stems)
-    REAL(wp), INTENT(IN) ::  diffusivity             !< diffusivity
-    REAL(wp), INTENT(IN) ::  p                       !< pressure (Pa)
-    REAL(wp), INTENT(IN) ::  conc_ijk_ugm3           !< actual atmospheric concentration (ug/m3), in DEPAC=Catm
-    REAL(wp), INTENT(IN) ::  ra                      !< aerodynamic resistance (s/m)
-    REAL(wp), INTENT(IN) ::  rb                      !< boundary layer resistance (s/m)
-
+
+    REAL(wp), INTENT(OUT) ::  ccomp_tot              !< total compensation point (ug/m3)
+!                                                    !< [= 0 for species that don't have a compensation
     REAL(wp), INTENT(OUT) ::  rc_tot                 !< total canopy resistance Rc (s/m)
-    REAL(wp), INTENT(OUT) ::  ccomp_tot              !< total compensation point (ug/m3)
-!                                                     !< [= 0 for species that don't have a compensation
+
 !-- Local variables:
 !
@@ -4899 +4885 @@
 
     LOGICAL ::  ready                                !< Rc has been set:
-    !< = 1 -> constant Rc
-    !< = 2 -> temperature dependent Rc
+                                                     !< = 1 -> constant Rc
+                                                     !< = 2 -> temperature dependent Rc
 !
 !-- Vegetation indicators:
     LOGICAL ::  lai_present                          !< leaves are present for current land use type
-    LOGICAL ::  sai_present                          !< vegetation is present for current land use type
-
+    LOGICAL ::  sai_present                          !< vegetation is present for current land use
+                                                     !< type
+
+    REAL(wp) ::  csoil                               !< soil compensation point (ug/m3)
+    REAL(wp) ::  cstom                               !< stomatal compensation point (ug/m3)
+    REAL(wp) ::  cw                                  !< external leaf surface compensation point
+                                                     !< (ug/m3)
+    REAL(wp) ::  gc_tot                              !< total canopy conductance (m/s)
+    REAL(wp) ::  gsoil_eff                           !< effective soil conductance (m/s)
+    REAL(wp) ::  gstom                               !< stomatal conductance (m/s)
+    REAL(wp) ::  gw                                  !< external leaf conductance (m/s)
 !    REAL(wp) ::  laimax                              !< maximum leaf area index (-)
-    REAL(wp) ::  gw                                  !< external leaf conductance (m/s)
-    REAL(wp) ::  gstom                               !< stomatal conductance (m/s)
-    REAL(wp) ::  gsoil_eff                           !< effective soil conductance (m/s)
-    REAL(wp) ::  gc_tot                              !< total canopy conductance (m/s)
-    REAL(wp) ::  cw                                  !< external leaf surface compensation point (ug/m3)
-    REAL(wp) ::  cstom                               !< stomatal compensation point (ug/m3)
-    REAL(wp) ::  csoil                               !< soil compensation point (ug/m3)
 !
 !-- Next statement is just to avoid compiler warning about unused variable
@@ -4931 +4919 @@
 
     CASE ( 'NO', 'no' )
-       icmp = icmp_no 
+       icmp = icmp_no
 
     CASE ( 'NH3', 'nh3' )
@@ -4979 +4967 @@
 !
 !--    External conductance:
-       CALL rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai,gw )         
+       CALL rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai,gw )
 !
 !--    Stomatal conductance:
-       CALL rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity, gstom, p )
+       CALL rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity,  &
+                      gstom, p )
 !
 !--    Effective soil conductance:
@@ -5000 +4989 @@
 !
 !--    Effective Rc based on compensation points:
-!--        IF ( present(rc_eff) ) then
+!--        IF ( PRESENT( rc_eff ) )  THEN
 !--          check on required arguments:
-!--           IF ( (.not. present(conc_ijk_ugm3)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then
-!--              stop 'output argument rc_eff requires input arguments conc_ijk_ugm3, ra and rb'
-!--           END IF
-!
-!--       Compute rc_eff :
-       !      CALL rc_comp_point_rc_eff(ccomp_tot,conc_ijk_ugm3,ra,rb,rc_tot,rc_eff)
-       !   ENDIF
+!--           IF ( ( .NOT. PRESENT( conc_ijk_ugm3 ) )  .OR.  ( .NOT. PRESENT( ra ) )  .OR.         &
+!--                ( .NOT. PRESENT( rb ) ) )  THEN
+!--              STOP 'output argument rc_eff requires input arguments conc_ijk_ugm3, ra and rb'
+!--           ENDIF
+!
+!--           Compute rc_eff :
+!--           CALL rc_comp_point_rc_eff(ccomp_tot,conc_ijk_ugm3,ra,rb,rc_tot,rc_eff)
+!--        ENDIF
     ENDIF
 
@@ -5014 +5004 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine to compute total canopy resistance in special cases
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE rc_special( icmp, compnam, lu, t, nwet, rc_tot, ready, ccomp_tot )
 
-   
+
     CHARACTER(LEN=*), INTENT(IN)  ::  compnam     !< component name
 
-    INTEGER(iwp), INTENT(IN)  ::  icmp            !< component index      
+    INTEGER(iwp), INTENT(IN)  ::  icmp            !< component index
     INTEGER(iwp), INTENT(IN)  ::  lu              !< land use type, lu = 1,...,nlu
     INTEGER(iwp), INTENT(IN)  ::  nwet            !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
 
-    REAL(wp), INTENT(IN)  ::  t                   !< temperature (C)
-
-    REAL(wp), INTENT(OUT) ::  rc_tot             !< total canopy resistance Rc (s/m)
-    REAL(wp), INTENT(OUT) ::  ccomp_tot          !< total compensation point (ug/m3)
-
     LOGICAL, INTENT(OUT) ::  ready               !< Rc has been set
                                                  !< = 1 -> constant Rc
+
+    REAL(wp), INTENT(IN)  ::  t                   !< temperature (C)
+
+    REAL(wp), INTENT(OUT) ::  ccomp_tot          !< total compensation point (ug/m3)
+    REAL(wp), INTENT(OUT) ::  rc_tot             !< total canopy resistance Rc (s/m)
+
 !
 !-- Next line is to avoid compiler warning about unused variable
@@ -5055 +5046 @@
        ELSE
 !
-!--       all other circumstances:
+!--       All other circumstances:
           rc_tot = 10.0_wp
        ENDIF
@@ -5085 +5076 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine to compute external conductance
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai, gw )
 
@@ -5095 +5086 @@
 !-- Input/output variables:
     CHARACTER(LEN=*), INTENT(IN) ::  compnam      !< component name
-
-    INTEGER(iwp), INTENT(IN) ::  nwet             !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     INTEGER(iwp), INTENT(IN) ::  iratns           !< index for NH3/SO2 ratio;
                                                   !< iratns = 1: low NH3/SO2
                                                   !< iratns = 2: high NH3/SO2
                                                   !< iratns = 3: very low NH3/SO2
+
+    INTEGER(iwp), INTENT(IN) ::  nwet             !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
+
     LOGICAL, INTENT(IN) ::  sai_present
 
-    REAL(wp), INTENT(IN) ::  t                    !< temperature (C)
     REAL(wp), INTENT(IN) ::  rh                   !< relative humidity (%)
     REAL(wp), INTENT(IN) ::  sai                  !< one-sided leaf area index (-)
+    REAL(wp), INTENT(IN) ::  t                    !< temperature (C)
 
     REAL(wp), INTENT(OUT) ::  gw                  !< external leaf conductance (m/s)
@@ -5126 +5118 @@
        CALL rw_nh3_sutton( t, rh, sai_present, gw )
 !
-!--    conversion from leaf resistance to canopy resistance by multiplying with sai:
+!--    Conversion from leaf resistance to canopy resistance by multiplying with sai:
        gw = sai * gw
 
@@ -5137 +5129 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine to compute external leaf conductance for SO2
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE rw_so2( t, nwet, rh, iratns, sai_present, gw )
 
 !
 !-- Input/output variables:
-    INTEGER(iwp), INTENT(IN) ::  nwet        !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     INTEGER(iwp), INTENT(IN) ::  iratns      !< index for NH3/SO2 ratio:
                                              !< iratns = 1: low NH3/SO2
                                              !< iratns = 2: high NH3/SO2
                                              !< iratns = 3: very low NH3/SO2
+    INTEGER(iwp), INTENT(IN) ::  nwet        !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
+
     LOGICAL, INTENT(IN) ::  sai_present
 
+    REAL(wp), INTENT(IN) ::  rh              !< relative humidity (%)
     REAL(wp), INTENT(IN) ::  t               !< temperature (C)
-    REAL(wp), INTENT(IN) ::  rh              !< relative humidity (%)   
 
     REAL(wp), INTENT(OUT) ::  gw             !< external leaf conductance (m/s)
@@ -5172 +5165 @@
           IF ( t > -1.0_wp )  THEN
              IF ( rh < 81.3_wp )  THEN
-                rw = 25000.0_wp * exp( -0.0693_wp * rh )
+                rw = 25000.0_wp * EXP( -0.0693_wp * rh )
              ELSE
-                rw = 0.58e12 * exp( -0.278_wp * rh ) + 10.0_wp
+                rw = 0.58e12 * EXP( -0.278_wp * rh ) + 10.0_wp
              ENDIF
           ELSE
@@ -5193 +5186 @@
        ENDIF
 !
-!--    very low NH3/SO2 ratio:
+!--    Very low NH3/SO2 ratio:
        IF ( iratns == iratns_very_low ) rw = rw + 50.0_wp
 !
@@ -5200 +5193 @@
     ELSE
 !
-!--      no vegetation:
+!--      No vegetation:
        gw = 0.0_wp
     ENDIF
@@ -5207 +5200 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Subroutine to compute external leaf conductance for NH3,
-!>                  following Sutton & Fowler, 1993
-!------------------------------------------------------------------------------!
+!> Subroutine to compute external leaf conductance for NH3, following Sutton & Fowler, 1993
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE rw_nh3_sutton( tsurf, rh,sai_present, gw )
 
@@ -5219 +5211 @@
     LOGICAL, INTENT(IN) ::  sai_present
 
+    REAL(wp), INTENT(IN) ::  rh             !< relative humidity (%)
     REAL(wp), INTENT(IN) ::  tsurf          !< surface temperature (C)
-    REAL(wp), INTENT(IN) ::  rh             !< relative humidity (%)
 
     REAL(wp), INTENT(OUT) ::  gw            !< external leaf conductance (m/s)
@@ -5239 +5231 @@
 !
 !--    External resistance according to Sutton & Fowler, 1993
-       rw = 2.0_wp * exp( ( 100.0_wp - rh ) / 12.0_wp )
+       rw = 2.0_wp * EXP( ( 100.0_wp - rh ) / 12.0_wp )
        rw = sai_grass_haarweg * rw
 !
@@ -5255 +5247 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Subroutine to compute external leaf conductance 
-!------------------------------------------------------------------------------!
+!> Subroutine to compute external leaf conductance
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE rw_constant( rw_val, sai_present, gw )
 
@@ -5266 +5258 @@
     LOGICAL, INTENT(IN) ::  sai_present
 
-    REAL(wp), INTENT(IN) ::  rw_val       !< constant value of Rw    
+    REAL(wp), INTENT(IN) ::  rw_val       !< constant value of Rw
 
     REAL(wp), INTENT(OUT) ::  gw          !< wernal leaf conductance (m/s)
 !
 !-- Compute conductance:
-    IF ( sai_present .AND. .NOT.missing(rw_val) )  THEN
+    IF ( sai_present  .AND.  .NOT. missing(rw_val) )  THEN
        gw = 1.0_wp / rw_val
     ELSE
@@ -5280 +5272 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Subroutine to compute stomatal conductance 
-!------------------------------------------------------------------------------!
- SUBROUTINE rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity, gstom, p ) 
+!> Subroutine to compute stomatal conductance
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity,  &
+                      gstom, p )
 
 !
@@ -5296 +5289 @@
     LOGICAL, INTENT(IN) ::  lai_present
 
+    REAL(wp), INTENT(IN) ::  diffusivity           !< diffusion coefficient of the gas involved
     REAL(wp), INTENT(IN) ::  lai                   !< one-sided leaf area index
+    REAL(wp), INTENT(IN) ::  rh                    !< relative humidity (%)
+    REAL(wp), INTENT(IN) ::  sinphi                !< sin of solar elevation angle
     REAL(wp), INTENT(IN) ::  solar_rad             !< solar radiation, dirict+diffuse (W/m2)
-    REAL(wp), INTENT(IN) ::  sinphi                !< sin of solar elevation angle
     REAL(wp), INTENT(IN) ::  t                     !< temperature (C)
-    REAL(wp), INTENT(IN) ::  rh                    !< relative humidity (%)
-    REAL(wp), INTENT(IN) ::  diffusivity           !< diffusion coefficient of the gas involved
 
     REAL(wp), OPTIONAL,INTENT(IN) :: p             !< pressure (Pa)
@@ -5324 +5317 @@
     CASE( 'NO2', 'O3', 'SO2', 'NH3' )
 !
-!--    if vegetation present:
+!--    If vegetation present:
        IF ( lai_present )  THEN
 
@@ -5330 +5323 @@
              CALL rc_get_vpd( t, rh, vpd )
              CALL rc_gstom_emb( lu, solar_rad, t, vpd, lai_present, lai, sinphi, gstom, p )
-             gstom = gstom * diffusivity / dO3       !< Gstom of Emberson is derived for ozone 
+             gstom = gstom * diffusivity / dO3       !< Gstom of Emberson is derived for ozone
           ELSE
              gstom = 0.0_wp
@@ -5336 +5329 @@
        ELSE
 !
-!--       no vegetation; zero conductance (infinite resistance):
+!--       No vegetation; zero conductance (infinite resistance):
           gstom = 0.0_wp
        ENDIF
@@ -5348 +5341 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine to compute stomatal conductance according to Emberson
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE rc_gstom_emb( lu, solar_rad, T, vpd, lai_present, lai, sinp, Gsto, p )
 !
@@ -5360 +5353 @@
 !>     Updated and extended.
 !>   2009-09, Arjo Segers, TNO
-!>     Limitted temperature influence to range to avoid
-!>     floating point exceptions.
+!>     Limitted temperature influence to range to avoid floating point exceptions.
 
 !> Method
@@ -5391 +5383 @@
     LOGICAL, INTENT(IN) ::  lai_present
 
+    REAL(wp), INTENT(IN) ::  lai                !< one-sided leaf area index
+    REAL(wp), INTENT(IN) ::  sinp               !< sin of solar elevation angle
     REAL(wp), INTENT(IN) ::  solar_rad          !< solar radiation, dirict+diffuse (W/m2)
     REAL(wp), INTENT(IN) ::  t                  !< temperature (C)
     REAL(wp), INTENT(IN) ::  vpd                !< vapour pressure deficit (kPa)
 
-    REAL(wp), INTENT(IN) ::  lai                !< one-sided leaf area index
-    REAL(wp), INTENT(IN) ::  sinp               !< sin of solar elevation angle
 
     REAL(wp), OPTIONAL, INTENT(IN) ::  p        !< pressure (Pa)
@@ -5403 +5395 @@
 !
 !-- Local variables:
+    REAL(wp) ::  bt
+    REAL(wp) ::  f_env
     REAL(wp) ::  f_light
     REAL(wp) ::  f_phen
     REAL(wp) ::  f_temp
+    REAL(wp) ::  f_swp
     REAL(wp) ::  f_vpd
-    REAL(wp) ::  f_swp
-    REAL(wp) ::  bt
-    REAL(wp) ::  f_env
+    REAL(wp) ::  laishade
+    REAL(wp) ::  laisun
+    REAL(wp) ::  pardiff
     REAL(wp) ::  pardir
-    REAL(wp) ::  pardiff
     REAL(wp) ::  parshade
     REAL(wp) ::  parsun
-    REAL(wp) ::  laisun
-    REAL(wp) ::  laishade
+    REAL(wp) ::  pres
     REAL(wp) ::  sinphi
-    REAL(wp) ::  pres
+
     REAL(wp), PARAMETER ::  p_sealevel = 1.01325e05    !< Pa
 !
@@ -5423 +5416 @@
     IF ( lai_present )  THEN
 
-       ! calculation of correction factors for stomatal conductance
-       IF ( sinp <= 0.0_wp )  THEN  
+       ! Calculation of correction factors for stomatal conductance
+       IF ( sinp <= 0.0_wp )  THEN
           sinphi = 0.0001_wp
        ELSE
@@ -5430 +5423 @@
        END IF
 !
-!--    ratio between actual and sea-level pressure is used
-!--    to correct for height in the computation of par;
-!--    should not exceed sea-level pressure therefore ...
-       IF (  present(p) )  THEN
-          pres = min( p, p_sealevel )
-       ELSE 
+!--    Ratio between actual and sea-level pressure is used to correct for height in the computation
+!--    of par; should not exceed sea-level pressure therefore ...
+       IF (  PRESENT( p ) )  THEN
+          pres = MIN( p, p_sealevel )
+       ELSE
           pres = p_sealevel
        ENDIF
 !
-!--    direct and diffuse par, Photoactive (=visible) radiation:
+!--    Direct and diffuse par, Photoactive (=visible) radiation:
        CALL par_dir_diff( solar_rad, sinphi, pres, p_sealevel, pardir, pardiff )
 !
-!--    par for shaded leaves (canopy averaged):
-       parshade = pardiff * exp( -0.5 * lai**0.7 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) * exp( -sinphi )     !< Norman,1982
-       IF ( solar_rad > 200.0_wp .AND. lai > 2.5_wp )  THEN
-          parshade = pardiff * exp( -0.5 * lai**0.8 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) * exp( -sinphi )  !< Zhang et al., 2001
+!--    Par for shaded leaves (canopy averaged):
+       parshade = pardiff * EXP( -0.5 * lai**0.7 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) *         &
+                  EXP( -sinphi )     !< Norman,1982
+       IF ( solar_rad > 200.0_wp  .AND.  lai > 2.5_wp )  THEN
+          parshade = pardiff * EXP( -0.5 * lai**0.8 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) *      &
+                     EXP( -sinphi )  !< Zhang et al., 2001
        END IF
 !
-!--    par for sunlit leaves (canopy averaged):
+!--    Par for sunlit leaves (canopy averaged):
 !--    alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5
        parsun = pardir * 0.5/sinphi + parshade             !< Norman, 1982
-       IF ( solar_rad > 200.0_wp .AND. lai > 2.5_wp )  THEN
+       IF ( solar_rad > 200.0_wp  .AND.  lai > 2.5_wp )  THEN
           parsun = pardir**0.8 * 0.5 / sinphi + parshade   !< Zhang et al., 2001
        END IF
 !
-!--    leaf area index for sunlit and shaded leaves:
+!--    Leaf area index for sunlit and shaded leaves:
        IF ( sinphi > 0 )  THEN
-          laisun = 2 * sinphi * ( 1 - exp( -0.5 * lai / sinphi ) )
+          laisun = 2 * sinphi * ( 1 - EXP( -0.5 * lai / sinphi ) )
           laishade = lai - laisun
        ELSE
@@ -5464 +5458 @@
        END IF
 
-       f_light = ( laisun * ( 1 - exp( -1.0_wp * alpha(lu) * parsun ) ) + &
-            laishade * ( 1 - exp( -1.0_wp * alpha(lu) * parshade ) ) ) / lai 
+       f_light = ( laisun * ( 1 - EXP( -1.0_wp * alpha(lu) * parsun ) ) +                          &
+                   laishade * ( 1 - EXP( -1.0_wp * alpha(lu) * parshade ) ) ) / lai
 
        f_light = MAX(f_light,f_min(lu))
 !
-!--    temperature influence; only non-zero within range [tmin,tmax]:
-       IF ( ( tmin(lu) < t ) .AND. ( t < tmax(lu) ) )  THEN
+!--    Temperature influence; only non-zero within range [tmin,tmax]:
+       IF ( ( tmin(lu) < t )  .AND.  ( t < tmax(lu) ) )  THEN
           bt = ( tmax(lu) - topt(lu) ) / ( topt(lu) - tmin(lu) )
-          f_temp = ( ( t - tmin(lu) ) / ( topt(lu) - tmin(lu) ) ) * ( ( tmax(lu) - t ) / ( tmax(lu) - topt(lu) ) )**bt
+          f_temp = ( ( t - tmin(lu) ) / ( topt(lu) - tmin(lu) ) ) *                                &
+                   ( ( tmax(lu) - t ) / ( tmax(lu) - topt(lu) ) )**bt
        ELSE
           f_temp = 0.0_wp
@@ -5478 +5473 @@
        f_temp = MAX( f_temp, f_min(lu) )
 !
-!--      vapour pressure deficit influence
-       f_vpd = MIN( 1.0_wp, ( ( 1.0_wp - f_min(lu) ) * ( vpd_min(lu) - vpd ) / ( vpd_min(lu) - vpd_max(lu) ) + f_min(lu) ) )
+!--    Vapour pressure deficit influence
+       f_vpd = MIN( 1.0_wp, ( ( 1.0_wp - f_min(lu) ) * ( vpd_min(lu) - vpd ) /                     &
+                              ( vpd_min(lu) - vpd_max(lu) ) + f_min(lu) ) )
        f_vpd = MAX( f_vpd, f_min(lu) )
 
        f_swp = 1.0_wp
 !
-!--      influence of phenology on stom. conductance
-!--      ignored for now in DEPAC since influence of f_phen on lu classes in use is negligible.
-!--      When other EMEP classes (e.g. med. broadleaf) are used f_phen might be too important to ignore
+!--    Influence of phenology on stom. conductance
+!--    Ignored for now in DEPAC since influence of f_phen on lu classes in use is negligible.
+!--    When other EMEP classes (e.g. med. broadleaf) are used f_phen might be too important to
+!--    ignore.
        f_phen = 1.0_wp
 !
-!--      evaluate total stomatal conductance
+!--    Evaluate total stomatal conductance
        f_env = f_temp * f_vpd * f_swp
        f_env = MAX( f_env,f_min(lu) )
        gsto = g_max(lu) * f_light * f_phen * f_env
 !
-!--      gstom expressed per m2 leafarea;
-!--      this is converted with lai to m2 surface.
+!--    gstom expressed per m2 leafarea;
+!--    This is converted with lai to m2 surface.
        gsto = lai * gsto    ! in m/s
 
@@ -5505 +5502 @@
 
 
- !-------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------!
  !> par_dir_diff
- !>     Weiss, A., Norman, J.M. (1985) Partitioning solar radiation into direct and
- !>     diffuse, visible and near-infrared components. Agric. Forest Meteorol.
- !>     34, 205-213.
+ !>     Weiss, A., Norman, J.M. (1985) Partitioning solar radiation into direct and diffuse, visible
+ !>     and near-infrared components. Agric. Forest Meteorol. 34, 205-213.
  !>     From a SUBROUTINE obtained from Leiming Zhang,
  !>     Meteorological Service of Canada
@@ -5515 +5511 @@
  !>     Willem Asman set it to global radiation (here defined as solar radiation, dirict+diffuse)
  !>
- !>     @todo Check/connect/replace with radiation_model_mod variables    
- !-------------------------------------------------------------------
+ !>     @todo Check/connect/replace with radiation_model_mod variables
+ !-------------------------------------------------------------------------------------------------!
  SUBROUTINE par_dir_diff( solar_rad, sinphi, pres, pres_0, par_dir, par_diff )
 
 
-    REAL(wp), INTENT(IN) ::  solar_rad       !< solar radiation, dirict+diffuse (W m-2)
-    REAL(wp), INTENT(IN) ::  sinphi          !< sine of the solar elevation
     REAL(wp), INTENT(IN) ::  pres            !< actual pressure (to correct for height) (Pa)
     REAL(wp), INTENT(IN) ::  pres_0          !< pressure at sea level (Pa)
-
-    REAL(wp), INTENT(OUT) ::  par_dir        !< par direct : visible (photoactive) direct beam radiation (W m-2)
-    REAL(wp), INTENT(OUT) ::  par_diff       !< par diffuse: visible (photoactive) diffuse radiation (W m-2)
-
-
-    REAL(wp) ::  sv                          !< total visible radiation
+    REAL(wp), INTENT(IN) ::  sinphi          !< sine of the solar elevation
+    REAL(wp), INTENT(IN) ::  solar_rad       !< solar radiation, dirict+diffuse (W m-2)
+
+    REAL(wp), INTENT(OUT) ::  par_diff       !< par diffuse: visible (photoactive) diffuse radiation
+                                             !< (W m-2)
+    REAL(wp), INTENT(OUT) ::  par_dir        !< par direct : visible (photoactive) direct beam
+                                             !< radiation (W m-2)
+
     REAL(wp) ::  fv                          !< par direct beam fraction (dimensionless)
-    REAL(wp) ::  ratio                       !< ratio measured to potential solar radiation (dimensionless)
-    REAL(wp) ::  rdm                         !< potential direct beam near-infrared radiation (W m-2); "potential" means clear-sky
+    REAL(wp) ::  ratio                       !< ratio measured to potential solar radiation
+                                             !< (dimensionless)
+    REAL(wp) ::  rdm                         !< potential direct beam near-infrared radiation
+                                             !< (W m-2); "potential" means clear-sky
     REAL(wp) ::  rdn                         !< potential diffuse near-infrared radiation (W m-2)
     REAL(wp) ::  rdu                         !< visible (par) direct beam radiation (W m-2)
@@ -5538 +5536 @@
     REAL(wp) ::  rn                          !< near-infrared radiation (W m-2)
     REAL(wp) ::  rv                          !< visible radiation (W m-2)
-    REAL(wp) ::  ww                          !< water absorption in the near infrared for 10 mm of precipitable water
+    REAL(wp) ::  sv                          !< total visible radiation
+    REAL(wp) ::  ww                          !< water absorption in the near infrared for 10 mm of
+                                             !< precipitable water
 
 !
 !-- Calculate visible (PAR) direct beam radiation
-!-- 600 W m-2 represents average amount of par (400-700 nm wavelength)
-!-- at the top of the atmosphere; this is roughly 0.45*solar constant (solar constant=1320 Wm-2)
-    rdu = 600.0_wp* exp( -0.185_wp * ( pres / pres_0 ) / sinphi ) * sinphi
+!-- 600 W m-2 represents average amount of par (400-700 nm wavelength) at the top of the atmosphere;
+!-- this is roughly 0.45*solar constant (solar constant=1320 Wm-2)
+    rdu = 600.0_wp* EXP( -0.185_wp * ( pres / pres_0 ) / sinphi ) * sinphi
 !
 !-- Calculate potential visible diffuse radiation
@@ -5550 +5550 @@
 !
 !-- Calculate the water absorption in the-near infrared
-    ww = 1320 * 10**( -1.195_wp + 0.4459_wp * log10( 1.0_wp / sinphi ) - 0.0345_wp * ( log10( 1.0_wp / sinphi ) )**2 )
+    ww = 1320 * 10**( -1.195_wp + 0.4459_wp * LOG10( 1.0_wp / sinphi ) - 0.0345_wp *               &
+                      ( LOG10( 1.0_wp / sinphi ) )**2 )
 !
 !-- Calculate potential direct beam near-infrared radiation
-    rdm = (720.0_wp * exp(-0.06_wp * (pres / pres_0) / sinphi ) - ww ) * sinphi     !< 720 = solar constant - 600
+    rdm = ( 720.0_wp * EXP( -0.06_wp * ( pres / pres_0) / sinphi ) - ww ) * sinphi  !< 720 = solar
+                                                                                    !< constant - 600
 !
 !-- Calculate potential diffuse near-infrared radiation
@@ -5563 +5565 @@
 !
 !-- Compute ratio between input global radiation (here defined as solar radiation, dirict+diffuse)
-!-- and total radiation computed here
+!-- and total radiation computed here.
     ratio = MIN( 0.89_wp, solar_rad / ( rv + rn ) )
 !
@@ -5571 +5573 @@
 !-- Calculate fraction of par in the direct beam
     fv = MIN( 0.99_wp, ( 0.9_wp - ratio ) / 0.7_wp )              !< help variable
-    fv = MAX( 0.01_wp, rdu / rv * ( 1.0_wp - fv**0.6667_wp ) )    !< fraction of par in the direct beam
+    fv = MAX( 0.01_wp, rdu / rv * ( 1.0_wp - fv**0.6667_wp ) )    !< fraction of par in the direct
+                                                                  !< beam
 !
 !-- Compute direct and diffuse parts of par
@@ -5579 +5582 @@
  END SUBROUTINE par_dir_diff
 
- 
- !-------------------------------------------------------------------
- !> rc_get_vpd: get vapour pressure deficit (kPa)
- !-------------------------------------------------------------------
+
+!--------------------------------------------------------------------------------------------------!
+!> rc_get_vpd: get vapour pressure deficit (kPa)
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE rc_get_vpd( temp, rh, vpd )
 
 !
 !-- Input/output variables:
+    REAL(wp), INTENT(IN) ::  rh    !< relative humidity (%)
     REAL(wp), INTENT(IN) ::  temp    !< temperature (C)
-    REAL(wp), INTENT(IN) ::  rh    !< relative humidity (%)
 
     REAL(wp), INTENT(OUT) ::  vpd    !< vapour pressure deficit (kPa)
@@ -5595 +5598 @@
     REAL(wp) ::  esat
 !
-!-- fit parameters:
+!-- Fit parameters:
     REAL(wp), PARAMETER ::  a1 = 6.113718e-01
     REAL(wp), PARAMETER ::  a2 = 4.43839e-02
@@ -5603 +5606 @@
     REAL(wp), PARAMETER ::  a6 = 3.0e-09
 !
-!-- esat is saturation vapour pressure (kPa) at temp(C) following Monteith(1973)
+!-- esat is saturation vapour pressure (kPa) at temp(C) following Monteith (1973)
     esat = a1 + a2 * temp + a3 * temp**2 + a4 * temp**3 + a5 * temp**4 + a6 * temp**5
     vpd  = esat * ( 1 - rh / 100 )
@@ -5610 +5613 @@
 
 
- !-------------------------------------------------------------------
- !> rc_gsoil_eff: compute effective soil conductance
- !-------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------!
+!> rc_gsoil_eff: compute effective soil conductance
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE rc_gsoil_eff( icmp, lu, sai, ust, nwet, t, gsoil_eff )
 
@@ -5618 +5621 @@
 !-- Input/output variables:
     INTEGER(iwp), INTENT(IN) ::  icmp          !< component index
-    INTEGER(iwp), INTENT(IN) ::  lu            !< land use type, lu = 1,..., nlu
     INTEGER(iwp), INTENT(IN) ::  nwet          !< index for wetness
                                                !< nwet = 0 -> dry; nwet = 1 -> wet; nwet = 9 -> snow
                                                !< N.B. this routine cannot be called with nwet = 9,
                                                !< nwet = 9 should be handled outside this routine.
+    INTEGER(iwp), INTENT(IN) ::  lu            !< land use type, lu = 1,..., nlu
+
     REAL(wp), INTENT(IN) ::  sai               !< surface area index
+    REAL(wp), INTENT(IN) ::  t                 !< temperature (C)
     REAL(wp), INTENT(IN) ::  ust               !< friction velocity (m/s)
-    REAL(wp), INTENT(IN) ::  t                 !< temperature (C)
+
     REAL(wp), INTENT(OUT) ::  gsoil_eff        !< effective soil conductance (m/s)
 !
@@ -5634 +5639 @@
 !-- Soil resistance (numbers matched with lu_classes and component numbers)
     !     grs    ara    crp    cnf    dec    wat    urb   oth    des    ice    sav    trf    wai    med    sem
-    REAL(wp), PARAMETER ::  rsoil(nlu_dep,ncmp) = reshape( (/ &
+    REAL(wp), PARAMETER ::  rsoil(nlu_dep,ncmp) = RESHAPE( (/                                      &
          1000.,  200.,  200.,  200.,  200., 2000.,  400., 1000., 2000., 2000., 1000.,  200., 2000.,  200.,  400., &    !< O3
          1000., 1000., 1000., 1000., 1000.,   10., 1000., 1000., 1000.,  500., 1000., 1000.,   10., 1000., 1000., &    !< SO2
@@ -5640 +5645 @@
          -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., &    !< NO
          100.,  100.,  100.,  100.,  100.,   10.,  100.,  100.,  100., 1000.,  100.,  100.,   10.,  100.,  100.,  &    !< NH3
-         -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., &    !< CO 
-         -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., &    !< NO3 
+         -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., &    !< CO
+         -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., &    !< NO3
          -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., &    !< HNO3
          -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., &    !< N2O5
-         -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999. /),&  !< H2O2   
+         -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999. /),&  !< H2O2
          (/nlu_dep,ncmp/) )
 !
 !-- For                                          o3    so2   no2     no    nh3     co     no3    hno3   n2o5   h2o2
-    REAL(wp), PARAMETER ::  rsoil_wet(ncmp)    = (/2000., 10. , 2000., -999., 10.  , -999., -999., -999., -999., -999./)
-    REAL(wp), PARAMETER ::  rsoil_frozen(ncmp) = (/2000., 500., 2000., -999., 1000., -999., -999., -999., -999., -999./)
+    REAL(wp), PARAMETER ::  rsoil_frozen(ncmp) = (/ 2000., 500., 2000., -999., 1000., -999., -999., -999., -999., -999. /)
+    REAL(wp), PARAMETER ::  rsoil_wet(ncmp)    = (/ 2000., 10. , 2000., -999., 10.  , -999., -999., -999., -999., -999. /)
 !
 !-- Compute in canopy (in crop) resistance:
@@ -5700 +5705 @@
 
 
- !-------------------------------------------------------------------
- !> rc_rinc: compute in canopy (or in crop) resistance
- !> van Pul and Jacobs, 1993, BLM
- !-------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------!
+!> rc_rinc: compute in canopy (or in crop) resistance van Pul and Jacobs, 1993, BLM
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE rc_rinc( lu, sai, ust, rinc )
 
@@ -5715 +5719 @@
     REAL(wp), INTENT(OUT) ::  rinc           !< in canopy resistance (s/m)
 !
-!-- b = empirical constant for computation of rinc (in canopy resistance) (= 14 m-1 or -999 if not applicable)
+!-- b = empirical constant for computation of rinc (in canopy resistance) (= 14 m-1 or -999 if not
+!--     applicable)
 !-- h = vegetation height (m)                     gra  ara crop con dec wat   urb   oth   des   ice   sav   trf  wai  med semi
     REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  b = (/ -999, 14, 14, 14, 14, -999, -999, -999, -999, -999, -999, 14, -999,  &
@@ -5724 +5729 @@
 !-- Compute Rinc only for arable land, perm. crops, forest; otherwise Rinc = 0:
     IF ( b(lu) > 0.0_wp )  THEN
-!       !
+!
 !--    Check for u* > 0 (otherwise denominator = 0):
        IF ( ust > 0.0_wp )  THEN
@@ -5732 +5737 @@
        ENDIF
     ELSE
-       IF ( lu == ilu_grass .OR. lu == ilu_other  )  THEN
+       IF ( lu == ilu_grass  .OR.  lu == ilu_other )  THEN
           rinc = -999.0_wp     !< no deposition path for grass, other, and semi-natural
        ELSE
@@ -5742 +5747 @@
 
 
- !-------------------------------------------------------------------
- !> rc_rctot: compute total canopy (or surface) resistance Rc
- !-------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------!
+!> rc_rctot: compute total canopy (or surface) resistance Rc
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE rc_rctot( gstom, gsoil_eff, gw, gc_tot, rc_tot )
 
 !
 !-- Input/output variables:
+    REAL(wp), INTENT(IN) ::  gsoil_eff     !< effective soil conductance (s/m)
     REAL(wp), INTENT(IN) ::  gstom         !< stomatal conductance (s/m)
-    REAL(wp), INTENT(IN) ::  gsoil_eff     !< effective soil conductance (s/m)
     REAL(wp), INTENT(IN) ::  gw            !< external leaf conductance (s/m)
 
@@ -5769 +5774 @@
 
 
- !-------------------------------------------------------------------
- !> rc_comp_point_rc_eff: calculate the effective resistance Rc
- !> based on one or more compensation points
- !-------------------------------------------------------------------
- !> NH3rc (see depac v3.6 is based on Avero workshop Marc Sutton. p. 173.
- !> Sutton 1998 AE 473-480)
- !>
- !> Documentation by Ferd Sauter, 2008; see also documentation block in header of depac subroutine.
- !> FS 2009-01-29: variable names made consistent with DEPAC 
- !> FS 2009-03-04: use total compensation point
- !>
- !> C: with total compensation point   ! D: approximation of C
- !>                                    !    with classical approach
- !>  zr --------- Catm                 !  zr --------- Catm    
- !>         |                          !         |       
- !>         Ra                         !         Ra      
- !>         |                          !         |       
- !>         Rb                         !         Rb      
- !>         |                          !         |       
- !>  z0 --------- Cc                   !  z0 --------- Cc
- !>         |                          !         |              
- !>        Rc                          !        Rc_eff          
- !>         |                          !         |              
- !>     --------- Ccomp_tot            !     --------- C=0
- !>
- !>
- !> The effective Rc is defined such that instead of using
- !>
- !>   F = -vd*[Catm - Ccomp_tot]                                    (1)
- !>
- !> we can use the 'normal' flux formula
- !>
- !>   F = -vd'*Catm,                                                (2)
- !>
- !> with vd' = 1/(Ra + Rb + Rc')                                    (3)
- !>
- !> and Rc' the effective Rc (rc_eff).
- !>                                                (Catm - Ccomp_tot)
- !> vd'*Catm = vd*(Catm - Ccomp_tot) <=> vd' = vd* ------------------
- !>                                                      Catm
- !>
- !>                                        (Catm - Ccomp_tot)
- !> 1/(Ra + Rb + Rc') = (1/Ra + Rb + Rc) * ------------------
- !>                                              Catm
- !>
- !>                                          Catm
- !> (Ra + Rb + Rc') = (Ra + Rb + Rc) * ------------------
- !>                                     (Catm - Ccomp_tot)
- !>
- !>                              Catm
- !> Rc' = (Ra + Rb + Rc) * ------------------ - Ra - Rb
- !>                        (Catm - Ccomp_tot)
- !>
- !>                        Catm                           Catm
- !> Rc' = (Ra + Rb) [------------------ - 1 ] + Rc * ------------------
- !>                  (Catm - Ccomp_tot)              (Catm - Ccomp_tot)
- !>
- !> Rc' = [(Ra + Rb)*Ccomp_tot + Rc*Catm ] / (Catm - Ccomp_tot)
- !>
- ! -------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------!
+!> rc_comp_point_rc_eff: calculate the effective resistance Rc
+!> based on one or more compensation points
+!--------------------------------------------------------------------------------------------------!
+!> NH3rc (see depac v3.6 is based on Avero workshop Marc Sutton. p. 173. Sutton 1998 AE 473-480)
+!>
+!> Documentation by Ferd Sauter, 2008; see also documentation block in header of depac subroutine.
+!> FS 2009-01-29: variable names made consistent with DEPAC
+!> FS 2009-03-04: use total compensation point
+!>
+!> C: with total compensation point   ! D: approximation of C
+!>                                    !    with classical approach
+!>  zr --------- Catm                 !  zr --------- Catm
+!>         |                          !         |
+!>         Ra                         !         Ra
+!>         |                          !         |
+!>         Rb                         !         Rb
+!>         |                          !         |
+!>  z0 --------- Cc                   !  z0 --------- Cc
+!>         |                          !         |
+!>        Rc                          !        Rc_eff
+!>         |                          !         |
+!>     --------- Ccomp_tot            !     --------- C=0
+!>
+!>
+!> The effective Rc is defined such that instead of using
+!>
+!>   F = -vd*[Catm - Ccomp_tot]                                    (1)
+!>
+!> we can use the 'normal' flux formula
+!>
+!>   F = -vd'*Catm,                                                (2)
+!>
+!> with vd' = 1/(Ra + Rb + Rc')                                    (3)
+!>
+!> and Rc' the effective Rc (rc_eff).
+!>                                                (Catm - Ccomp_tot)
+!> vd'*Catm = vd*(Catm - Ccomp_tot) <=> vd' = vd* ------------------
+!>                                                      Catm
+!>
+!>                                        (Catm - Ccomp_tot)
+!> 1/(Ra + Rb + Rc') = (1/Ra + Rb + Rc) * ------------------
+!>                                              Catm
+!>
+!>                                          Catm
+!> (Ra + Rb + Rc') = (Ra + Rb + Rc) * ------------------
+!>                                     (Catm - Ccomp_tot)
+!>
+!>                              Catm
+!> Rc' = (Ra + Rb + Rc) * ------------------ - Ra - Rb
+!>                        (Catm - Ccomp_tot)
+!>
+!>                        Catm                           Catm
+!> Rc' = (Ra + Rb) [------------------ - 1 ] + Rc * ------------------
+!>                  (Catm - Ccomp_tot)              (Catm - Ccomp_tot)
+!>
+!> Rc' = [(Ra + Rb)*Ccomp_tot + Rc*Catm ] / (Catm - Ccomp_tot)
+!>
+! -------------------------------------------------------------------------------------------
 ! SUBROUTINE rc_comp_point_rc_eff( ccomp_tot, conc_ijk_ugm3, ra, rb, rc_tot, rc_eff )
 !
@@ -5843 +5847 @@
 !    !
 !!-- Compute effective resistance:
-!    IF (  ccomp_tot == 0.0_wp )  THEN
+!    IF ( ccomp_tot == 0.0_wp )  THEN
 !       !
-!!--    trace with no compensiation point ( or compensation point equal to zero)
+!!--    Trace with no compensiation point ( or compensation point equal to zero)
 !       rc_eff = rc_tot
 !
-!    ELSE IF ( ccomp_tot > 0.0_wp .AND. ( abs( conc_ijk_ugm3 - ccomp_tot ) < 1.e-8 ) )  THEN
+!    ELSE IF ( ccomp_tot > 0.0_wp  .AND.  ( abs( conc_ijk_ugm3 - ccomp_tot ) < 1.e-8 ) )  THEN
 !       !
-!!--   surface concentration (almost) equal to atmospheric concentration
+!!--   Surface concentration (almost) equal to atmospheric concentration
 !!--    no exchange between surface and atmosphere, infinite RC --> vd=0
 !       rc_eff = 9999999999.0_wp
@@ -5856 +5860 @@
 !    ELSE IF ( ccomp_tot > 0.0_wp )  THEN
 !       !
-!!--    compensation point available, calculate effective resistance
+!!--    Compensation point available, calculate effective resistance
 !       rc_eff = ( ( ra + rb ) * ccomp_tot + rc_tot * conc_ijk_ugm3 ) / ( conc_ijk_ugm3 - ccomp_tot )
 !
@@ -5866 +5870 @@
 !
 !    RETURN
-!    
+!
 ! END SUBROUTINE rc_comp_point_rc_eff
 
 
- !-------------------------------------------------------------------
- !> missing: check for data that correspond with a missing deposition path
- !>          this data is represented by -999
- !-------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------!
+!> missing: check for data that correspond with a missing deposition path this data is represented
+!>          by -999
+!--------------------------------------------------------------------------------------------------!
  LOGICAL function missing( x )
 
@@ -5879 +5883 @@
 
 !
-!-- bandwidth for checking (in)equalities of floats
+!-- Bandwidth for checking (in)equalities of floats
     REAL(wp), PARAMETER :: eps = 1.0e-5
 
-    missing = (abs(x + 999.0_wp) <= eps)
+    missing = ( ABS( x + 999.0_wp ) <= eps )
 
  END function missing
@@ -5891 +5895 @@
 !
 !-- in/out
-
     REAL(wp), INTENT(IN) ::  rhopart                 !< particle density (kg/m3)
     REAL(wp), INTENT(IN) ::  partsize                !< particle size (m)
@@ -5907 +5910 @@
  END FUNCTION sedimentation_velocity
 
- 
- !------------------------------------------------------------------------
+
+ !-------------------------------------------------------------------------------------------------!
  !> Boundary-layer deposition resistance following Zhang (2001)
- !------------------------------------------------------------------------
- SUBROUTINE drydepo_aero_zhang_vd( vd, rs, vs1, partsize, slipcor, nwet, tsurf, dens1, viscos1, &
-      luc, ftop_lu, ustar )
-
-!
-!-- in/out 
-
+ !-------------------------------------------------------------------------------------------------!
+ SUBROUTINE drydepo_aero_zhang_vd( vd, rs, vs1, partsize, slipcor, nwet, tsurf, dens1, viscos1,    &
+                                   luc, ftop_lu, ustar )
+
+!
+!-- in/out
+    INTEGER(iwp), INTENT(IN) ::  luc         !< DEPAC LU
     INTEGER(iwp), INTENT(IN) ::  nwet        !< 1=rain, 9=snowcover
-    INTEGER(iwp), INTENT(IN) ::  luc         !< DEPAC LU
-
-    REAL(wp), INTENT(IN) ::  vs1             !< sedimentation velocity in lowest layer
+
+    REAL(wp), INTENT(IN) ::  dens1           !< air density (kg/m3) in lowest layer
+    REAL(wp), INTENT(IN) ::  ftop_lu         !< atmospheric resistnace Ra
     REAL(wp), INTENT(IN) ::  partsize        !< particle diameter (m)
     REAL(wp), INTENT(IN) ::  slipcor         !< slip correction factor
     REAL(wp), INTENT(IN) ::  tsurf           !< surface temperature (K)
-    REAL(wp), INTENT(IN) ::  dens1           !< air density (kg/m3) in lowest layer
+    REAL(wp), INTENT(IN) ::  ustar           !< friction velocity u*
     REAL(wp), INTENT(IN) ::  viscos1         !< air viscosity in lowest layer
-    REAL(wp), INTENT(IN) ::  ftop_lu         !< atmospheric resistnace Ra
-    REAL(wp), INTENT(IN) ::  ustar           !< friction velocity u*    
-
+    REAL(wp), INTENT(IN) ::  vs1             !< sedimentation velocity in lowest layer
+
+    REAL(wp), INTENT(OUT) ::  rs             !< sedimentaion resistance (s/m)
     REAL(wp), INTENT(OUT) ::  vd             !< deposition velocity (m/s)
-    REAL(wp), INTENT(OUT) ::  rs             !< sedimentaion resistance (s/m)
-!
-!-- constants 
-
+!
+!-- Constants
     REAL(wp), PARAMETER ::  grav     = 9.80665_wp             !< acceleration of gravity (m/s2)
-
     REAL(wp), PARAMETER ::  epsilon0 = 3.0_wp
     REAL(wp), PARAMETER ::  kb       = 1.38066E-23_wp
     REAL(wp), PARAMETER ::  pi       = 3.141592654_wp      !< pi
 
-    REAL(wp), PARAMETER :: alfa_lu(nlu_dep) = & 
-         (/1.2_wp,  1.2_wp,   1.2_wp,  1.0_wp,  1.0_wp,   100.0_wp, 1.5_wp,  1.2_wp, 50.0_wp, 100.0_wp, &
-               1.2_wp, 1.0_wp, 100.0_wp, 1.2_wp, 50.0_wp/)   
+    REAL(wp), PARAMETER :: alfa_lu(nlu_dep) = &
+         (/ 1.2_wp,  1.2_wp,   1.2_wp,  1.0_wp,  1.0_wp,   100.0_wp, 1.5_wp,  1.2_wp, 50.0_wp, 100.0_wp, &
+              1.2_wp, 1.0_wp, 100.0_wp, 1.2_wp, 50.0_wp /)
     REAL(wp), PARAMETER :: gamma_lu(nlu_dep) = &
-         (/0.54_wp, 0.54_wp,  0.54_wp, 0.56_wp, 0.56_wp,  0.50_wp,  0.56_wp, 0.54_wp, 0.58_wp, 0.50_wp, &
-              0.54_wp, 0.56_wp, 0.50_wp, 0.54_wp, 0.54_wp/)   
-    REAL(wp), PARAMETER ::A_lu(nlu_dep) = &   
-         (/3.0_wp,  3.0_wp,   2.0_wp,  2.0_wp,  7.0_wp, -99.0_wp, 10.0_wp, 3.0_wp, -99.0_wp, -99.0_wp,  &
-              3.0_wp, 7.0_wp, -99.0_wp, 2.0_wp, -99.0_wp/)
-!
-!--   grass  arabl crops conif decid  water  urba  othr  desr  ice   sav  trf   wai  med   sem     
+         (/ 0.54_wp, 0.54_wp,  0.54_wp, 0.56_wp, 0.56_wp,  0.50_wp,  0.56_wp, 0.54_wp, 0.58_wp, 0.50_wp, &
+              0.54_wp, 0.56_wp, 0.50_wp, 0.54_wp, 0.54_wp /)
+    REAL(wp), PARAMETER ::A_lu(nlu_dep) = &
+         (/ 3.0_wp,  3.0_wp,   2.0_wp,  2.0_wp,  7.0_wp, -99.0_wp, 10.0_wp, 3.0_wp, -99.0_wp, -99.0_wp,  &
+              3.0_wp, 7.0_wp, -99.0_wp, 2.0_wp, -99.0_wp /)
+!
+!--   grass  arabl crops conif decid  water  urba  othr  desr  ice   sav  trf   wai  med   sem
 !
 !-- local
+    REAL(wp) ::  diff_part
+    REAL(wp) ::  ebrown
+    REAL(wp) ::  eimpac
+    REAL(wp) ::  einterc
     REAL(wp) ::  kinvisc
-    REAL(wp) ::  diff_part
+    REAL(wp) ::  reffic
     REAL(wp) ::  schmidt
     REAL(wp) ::  stokes
-    REAL(wp) ::  Ebrown
-    REAL(wp) ::  Eimpac
-    REAL(wp) ::  Einterc
-    REAL(wp) ::  Reffic
-!
-!-- kinetic viscosity & diffusivity
+!
+!-- Kinetic viscosity & diffusivity
     kinvisc = viscos1 / dens1    !< only needed at surface
 
@@ -5970 +5970 @@
     schmidt = kinvisc / diff_part
 !
-!-- calculate collection efficiencie E
+!-- Calculate collection efficiencie E
     Ebrown = Schmidt**( -gamma_lu(luc) )    !< Brownian diffusion
 !
-!-- determine Stokes number, interception efficiency 
-!-- and sticking efficiency R (1 = no rebound)
-    IF ( luc == ilu_ice .OR. nwet==9 .OR. luc == ilu_water_sea .OR. luc == ilu_water_inland )  THEN
+!-- Determine Stokes number, interception efficiency and sticking efficiency R (1 = no rebound)
+    IF ( luc == ilu_ice  .OR.  nwet == 9  .OR.  luc == ilu_water_sea  .OR.                         &
+         luc == ilu_water_inland )  THEN
        stokes = vs1 * ustar**2 / ( grav * kinvisc )
-       Einterc = 0.0_wp
-       Reffic = 1.0_wp
-    ELSE IF ( luc == ilu_other .OR. luc == ilu_desert )  THEN     !<tundra of desert
+       einterc = 0.0_wp
+       reffic = 1.0_wp
+    ELSE IF ( luc == ilu_other  .OR.  luc == ilu_desert )  THEN     !<tundra of desert
        stokes = vs1 * ustar**2 / ( grav * kinvisc )
-       Einterc = 0.0_wp
-       Reffic = exp( -Stokes**0.5_wp )
+       einterc = 0.0_wp
+       reffic = EXP( - stokes**0.5_wp )
     ELSE
-       stokes = vs1 * ustar / ( grav * A_lu(luc) * 1.0E-3_wp )
-       Einterc = 0.5_wp * ( partsize / (A_lu(luc) * 1.0E-3_wp ) )**2
-       Reffic = exp( -Stokes**0.5_wp )
+       stokes = vs1 * ustar / ( grav * a_lu(luc) * 1.0E-3_wp )
+       einterc = 0.5_wp * ( partsize / ( a_lu(luc) * 1.0E-3_wp ) )**2
+       reffic = EXP( - stokes**0.5_wp )
     END IF
 !
-!-- when surface is wet all particles do not rebound:
-    IF ( nwet==1 )  Reffic = 1.0_wp
-!
-!-- determine impaction efficiency:
-    Eimpac = ( stokes / ( alfa_lu(luc) + stokes ) )**2
-!
-!-- sedimentation resistance:
-    rs = 1.0_wp / ( epsilon0 * MAX( 1.0E-5_wp, ustar ) * ( Ebrown + Eimpac + Einterc ) * Reffic )
-
-!-- deposition velocity according to Seinfeld and Pandis (2006; eq 19.7):
-!--  
+!-- When surface is wet all particles do not rebound:
+    IF ( nwet == 1 )  reffic = 1.0_wp
+!
+!-- Determine impaction efficiency:
+    eimpac = ( stokes / ( alfa_lu(luc) + stokes ) )**2
+!
+!-- Sedimentation resistance:
+    rs = 1.0_wp / ( epsilon0 * MAX( 1.0E-5_wp, ustar ) * ( ebrown + eimpac + einterc ) * reffic )
+
+!-- Deposition velocity according to Seinfeld and Pandis (2006; eq 19.7):
+!--
 !--              1
 !--      vd = ------------------ + vs
 !--           Ra + Rs + Ra*Rs*vs
-!--  
+!--
 !-- where: Rs = Rb (in Seinfeld and Pandis, 2006)
 
@@ -6012 +6012 @@
 
 
- !-------------------------------------------------------------------------------------
+ !-------------------------------------------------------------------------------------------------!
  !> Compute quasi-laminar boundary layer resistance as a function of landuse and tracer
- !> Original EMEP formulation by (Simpson et al, 2003) is used 
- !-------------------------------------------------------------------------------------
- SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffusivity, rb )    
-
-!
-!-- in/out 
-
+ !> Original EMEP formulation by (Simpson et al, 2003) is used.
+ !-------------------------------------------------------------------------------------------------!
+ SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffusivity, rb )
+
+!
+!-- in/out
     LOGICAL , INTENT(IN) ::  is_water
 
+    REAL(wp), INTENT(IN) ::  diffusivity          !< coefficient of diffusivity
+    REAL(wp), INTENT(IN) ::  ustar                !< friction velocity
     REAL(wp), INTENT(IN) ::  z0h                  !< roughness length for heat
-    REAL(wp), INTENT(IN) ::  ustar                !< friction velocity
-    REAL(wp), INTENT(IN) ::  diffusivity          !< coefficient of diffusivity
 
     REAL(wp), INTENT(OUT) ::  rb                  !< boundary layer resistance
@@ -6031 +6030 @@
 !-- const
 
+    REAL(wp), PARAMETER ::  kappa_stab = 0.35     !< von Karman constant
     REAL(wp), PARAMETER ::  thk = 0.19e-4         !< thermal diffusivity of dry air 20 C
-    REAL(wp), PARAMETER ::  kappa_stab = 0.35     !< von Karman constant
 !
 !-- Next line is to avoid compiler warning about unused variable
     IF ( is_water  .OR.  ( z0h + kappa_stab ) > 0.0_wp )  CONTINUE
 !
-!-- Use Simpson et al. (2003) 
+!-- Use Simpson et al. (2003)
 !-- @TODO: Check rb over water calculation, until then leave commented lines
 !--  IF ( is_water )  THEN
-!--   org: rb = 1.0_wp / (kappa_stab*MAX(0.01_wp,ustar)) * log(z0h/diffusivity*kappa_stab*MAX(0.01_wp,ustar))
-!--        rb = 1.0_wp / (kappa_stab*MAX(0.1_wp,ustar)) * log(z0h/diffusivity*kappa_stab*MAX(0.1_wp,ustar))
+!--   org: rb = 1.0_wp / ( kappa_stab * MAX( 0.01_wp, ustar ) ) * LOG( z0h / diffusivity * kappa_stab * MAX( 0.01_wp, ustar ) )
+!--        rb = 1.0_wp / ( kappa_stab * MAX( 0.1_wp, ustar ) ) *  LOG( z0h / diffusivity * kappa_stab * MAX( 0.1_wp, ustar ) )
 !--  ELSE
     rb = 5.0_wp / MAX( 0.01_wp, ustar ) * ( thk / diffusivity )**0.67_wp
@@ -6049 +6048 @@
 
 
- !-----------------------------------------------------------------
- !>  Compute water vapor partial pressure (e_w)
- !>  given specific humidity Q [(kg water)/(kg air)].
- !>
- !>  Use that gas law for volume V with temperature T
- !>  holds for the total mixture as well as the water part:
- !>
- !>    R T / V = p_air / n_air = p_water / n_water
- !>
- !>  thus:
- !>
- !>    p_water = p_air n_water / n_air
- !>
- !>  Use:
- !>    n_air =   m_air   /        xm_air
- !>            [kg air]  /  [(kg air)/(mole air)]
- !>  and:
- !>    n_water =  m_air * Q  /     xm_water
- !>              [kg water]  /  [(kg water)/(mole water)]
- !>  thus:
- !>    p_water = p_air Q / (xm_water/xm_air)
- !------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------!
+!>  Compute water vapor partial pressure (e_w) given specific humidity Q [(kg water)/(kg air)].
+!>
+!>  Use that gas law for volume V with temperature T holds for the total mixture as well as the
+!>  water part:
+!>
+!>    R T / V = p_air / n_air = p_water / n_water
+!>
+!>  thus:
+!>
+!>    p_water = p_air n_water / n_air
+!>
+!>  Use:
+!>    n_air =   m_air   /        xm_air
+!>            [kg air]  /  [(kg air)/(mole air)]
+!>  and:
+!>    n_water =  m_air * Q  /     xm_water
+!>              [kg water]  /  [(kg water)/(mole water)]
+!>  thus:
+!>    p_water = p_air Q / (xm_water/xm_air)
+!--------------------------------------------------------------------------------------------------!
 
  ELEMENTAL FUNCTION watervaporpartialpressure( q, p ) RESULT( p_w )
 
 !
-!-- in/out 
-
+!-- in/out
+    REAL(wp), INTENT(IN) ::  p                      !< air pressure [Pa]
     REAL(wp), INTENT(IN) ::  q                      !< specific humidity [(kg water)/(kg air)]
-    REAL(wp), INTENT(IN) ::  p                      !< air pressure [Pa]
 
     REAL(wp) ::  p_w                                !< water vapor partial pressure [Pa]
 !
-!-- const 
-
+!-- Const
     REAL(wp), PARAMETER  ::  eps = xm_h2o / xm_air  !< mole mass ratio ~ 0.622
 !
-!-- partial pressure of water vapor:
+!-- Partial pressure of water vapor:
     p_w = p * q / eps
 
- END function watervaporpartialpressure
-
-
- !------------------------------------------------------------------    
- !>  Saturation vapor pressure.
- !>  From (Stull 1988, eq. 7.5.2d):
- !>
- !>      e_sat = p0 exp( 17.67 * (T-273.16) / (T-29.66) )     [Pa]
- !>
- !>  where:
- !>      p0 = 611.2 [Pa]   : reference pressure
- !>
- !>  Arguments:
- !>      T  [K]  : air temperature
- !>  Result:
- !>      e_sat_w  [Pa]  : saturation vapor pressure
- !>
- !>  References:
- !>      Roland B. Stull, 1988
- !>      An introduction to boundary layer meteorology.
- !-----------------------------------------------------------------
-
+ END FUNCTION watervaporpartialpressure
+
+
+!--------------------------------------------------------------------------------------------------!
+!>  Saturation vapor pressure.
+!>  From (Stull 1988, eq. 7.5.2d):
+!>
+!>      e_sat = p0 exp( 17.67 * (T-273.16) / (T-29.66) )     [Pa]
+!>
+!>  where:
+!>      p0 = 611.2 [Pa]   : reference pressure
+!>
+!>  Arguments:
+!>      T  [K]  : air temperature
+!>  Result:
+!>      e_sat_w  [Pa]  : saturation vapor pressure
+!>
+!>  References:
+!>      Roland B. Stull, 1988
+!>      An introduction to boundary layer meteorology.
+!--------------------------------------------------------------------------------------------------!
  ELEMENTAL FUNCTION saturationvaporpressure( t ) RESULT( e_sat_w )
 
 !
 !-- in/out
-
     REAL(wp), INTENT(IN) ::  t            !< temperature [K]
 
     REAL(wp) ::  e_sat_w                  !< saturation vapor pressure  [Pa]
 !
-!-- const 
+!-- Const
     REAL(wp), PARAMETER ::  p0 = 611.2   !< base pressure [Pa]
 !
-!-- saturation vapor pressure:
-    e_sat_w = p0 * exp( 17.67_wp * ( t - 273.16_wp ) / ( t - 29.66_wp ) )    !< [Pa]
+!-- Saturation vapor pressure:
+    e_sat_w = p0 * EXP( 17.67_wp * ( t - 273.16_wp ) / ( t - 29.66_wp ) )    !< [Pa]
 
  END FUNCTION saturationvaporpressure
 
 
- !------------------------------------------------------------------------
- !>  Relative humidity RH [%] is by definition:
- !>
- !>           e_w             water vapor partial pressure
- !>    Rh = -------- * 100
- !>         e_sat_w           saturation vapor pressure
- !------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------!
+!>  Relative humidity RH [%] is by definition:
+!>
+!>           e_w             water vapor partial pressure
+!>    Rh = -------- * 100
+!>         e_sat_w           saturation vapor pressure
+!--------------------------------------------------------------------------------------------------!
 
  ELEMENTAL FUNCTION relativehumidity_from_specifichumidity( q, t, p ) RESULT( rh )
 
 !
-!-- in/out 
-
+!-- in/out
+    REAL(wp), INTENT(IN) ::  p    !< air pressure [Pa]
     REAL(wp), INTENT(IN) ::  q    !< specific humidity [(kg water)/(kg air)]
     REAL(wp), INTENT(IN) ::  t    !< temperature [K]
-    REAL(wp), INTENT(IN) ::  p    !< air pressure [Pa]
 
     REAL(wp) ::  rh               !< relative humidity [%]
 !
-!-- relative humidity:
+!-- Relative humidity:
     rh = watervaporpartialpressure( q, p ) / saturationvaporpressure( t ) * 100.0_wp
 
  END FUNCTION relativehumidity_from_specifichumidity
 
-     
+
  END MODULE chemistry_model_mod