Index: palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3847)
+++ palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3848)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id: chemistry_model_mod.f90 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 
 ! 
@@ -260,12 +263,9 @@
 ! ------------
 !> Chemistry model for PALM-4U
-!> @todo Extend chem_species type by nspec and nvar as addititional elements
+!> @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 Update/clean-up todo list! (FK)
-!> @todo Add routine chem_check_parameters, add checks for inconsistent or
-!>       unallowed parameter settings. 
-!>       CALL of chem_check_parameters from check_parameters. (FK)
 !> @todo Implement turbulent inflow of chem spcs in inflow_turbulence.
 !> @todo Separate boundary conditions for each chem spcs to be implemented
@@ -277,7 +277,4 @@
 !> @todo test nesting for chem spcs, was implemented by suehring (kanani)
 !> @todo chemistry error messages
-!> @todo Format this module according to PALM coding standard (see e.g. module
-!>       template under http://palm.muk.uni-hannover.de/mosaik/downloads/8 or
-!>       D3_coding_standard.pdf under https://palm.muk.uni-hannover.de/trac/downloads/16)
 ! 
 !------------------------------------------------------------------------------!
@@ -318,31 +315,27 @@
 
     LOGICAL ::  nest_chemistry = .TRUE.  !< flag for nesting mode of chemical species, independent on parent or not
-    !
-    !- Define chemical variables
-    TYPE   species_def
-       CHARACTER(LEN=15)                            ::  name
-       CHARACTER(LEN=15)                            ::  unit
-       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  conc
-       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  conc_av
-       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  conc_p
-       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  tconc_m
-       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   ::  cssws_av           
-       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   ::  flux_s_cs
-       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   ::  diss_s_cs
-       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) ::  flux_l_cs
-       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) ::  diss_l_cs
-       REAL(kind=wp), ALLOCATABLE, DIMENSION(:)     ::  conc_pr_init
+!
+!-  Define chemical variables within chem_species
+    TYPE species_def
+       CHARACTER(LEN=15)                            ::  name         !< name of chemical species
+       CHARACTER(LEN=15)                            ::  unit         !< unit (ppm for gases, kg m^-3 for aerosol tracers)
+       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  conc         !< concentrations of trace gases
+       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  conc_av      !< averaged concentrations
+       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  conc_p       !< conc at prognostic time level
+       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  tconc_m      !< weighted tendency of conc for previous sub-timestep (Runge-Kutta)
+       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   ::  cssws_av     !< averaged fluxes of trace gases at surface
+       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   ::  flux_s_cs    !< 6th-order advective flux at south face of grid box of chemical species (='cs')
+       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   ::  diss_s_cs    !< artificial numerical dissipation flux at south face of grid box of chemical species
+       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) ::  flux_l_cs    !< 6th-order advective flux at left face of grid box of chemical species (='cs')
+       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) ::  diss_l_cs    !< artificial numerical dissipation flux at left face of grid box of chemical species
+       REAL(kind=wp), ALLOCATABLE, DIMENSION(:)     ::  conc_pr_init !< initial profile of chemical species
     END TYPE species_def
-
-
-    TYPE   photols_def                                                            
-       CHARACTER(LEN=15)                            :: name
-       CHARACTER(LEN=15)                            :: unit
-       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     :: freq
+!
+!-- Define photolysis frequencies in phot_frequen
+    TYPE photols_def                                                            
+       CHARACTER(LEN=15)                            :: name          !< name of pgotolysis frequency
+       CHARACTER(LEN=15)                            :: unit          !< unit (1/s)
+       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     :: freq          !< photolysis frequency
     END TYPE photols_def
-
-
-    PUBLIC  species_def
-    PUBLIC  photols_def
 
 
@@ -350,14 +343,22 @@
     TYPE(photols_def), ALLOCATABLE, DIMENSION(:), TARGET, PUBLIC ::  phot_frequen  
 
-    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_1
-    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_2
-    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_3
-    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_av       
-    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  freq_1
-
-    INTEGER, DIMENSION(nkppctrl)                           ::  icntrl                            !< Fine tuning kpp
-    REAL(kind=wp), DIMENSION(nkppctrl)                     ::  rcntrl                            !< Fine tuning kpp
-    LOGICAL ::  decycle_chem_lr    = .FALSE.
-    LOGICAL ::  decycle_chem_ns    = .FALSE.
+    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 ::  spec_conc_av !< averaged concentrations of chemical species       
+    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)
+!
+!-- Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not
+!-- approproate for chemicals compounds since they may accumulate too much. 
+!-- If no proper boundary conditions from a DYNAMIC input file are available,
+!-- de-cycling applies the initial profiles at the inflow boundaries for 
+!-- Dirichlet boundary conditions
+    LOGICAL ::  decycle_chem_lr    = .FALSE.    !< switch for setting decycling in left-right direction
+    LOGICAL ::  decycle_chem_ns    = .FALSE.    !< switch for setting decycling in south-norht direction
     CHARACTER (LEN=20), DIMENSION(4) ::  decycle_method = &
          (/'dirichlet','dirichlet','dirichlet','dirichlet'/)
@@ -376,12 +377,12 @@
 
 
-    !
-    !-- Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010)
+!
+!-- 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                              
+!
+!-- DEPAC landuse classes as defined in LOTOS-EUROS model v2.1                              
     INTEGER(iwp) ::  ilu_grass              = 1        
     INTEGER(iwp) ::  ilu_arable             = 2        
@@ -400,18 +401,18 @@
     INTEGER(iwp) ::  ilu_semi_natural_veg   = 15 
 
-    !
-    !-- NH3/SO2 ratio regimes:
+!
+!-- NH3/SO2 ratio regimes:
     INTEGER(iwp), PARAMETER ::  iratns_low      = 1       !< low ratio NH3/SO2
     INTEGER(iwp), PARAMETER ::  iratns_high     = 2       !< high ratio NH3/SO2
     INTEGER(iwp), PARAMETER ::  iratns_very_low = 3       !< very low ratio NH3/SO2
-    !
-    !-- Default:
+!
+!-- Default:
     INTEGER, PARAMETER ::  iratns_default = iratns_low
-    !
-    !-- Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) to W m-2
+!
+!-- 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
-    !
-    !-- Set temperatures per land use for f_temp
+!
+!-- 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/)
@@ -420,16 +421,16 @@
     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 /)
-    !
-    !-- Set f_min:
+!
+!-- 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/)
 
-    !
-    !-- 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
+!
+!-- 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
-    !
-    !-- Set max, min for vapour pressure deficit vpd
+         270., 150., -999., 300., 422./)/41000.
+!
+!-- 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/) 
@@ -437,14 +438,14 @@
          3.25, -999., 2.8, 4.5/) 
 
-
     PUBLIC nest_chemistry
-    PUBLIC nspec
     PUBLIC nreact
-    PUBLIC nvar      
-    PUBLIC spc_names
+    PUBLIC nspec               !< number of gas phase chemical species including constant compound (e.g. N2)
+    PUBLIC nvar                !< number of variable gas phase chemical species (nvar <= nspec)
+    PUBLIC photols_def         !< type defining phot_frequen
+    PUBLIC species_def         !< type defining chem_species
+    PUBLIC spc_names           !< names of gas phase chemical species (come from KPP) (come from KPP)
     PUBLIC spec_conc_2  
-
-    !    
-    !-- Interface section
+!    
+!-- Interface section
     INTERFACE chem_3d_data_averaging
        MODULE PROCEDURE chem_3d_data_averaging 
@@ -607,22 +608,18 @@
          chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo 
 
-
-
  CONTAINS
 
- !
- !------------------------------------------------------------------------------!
- !
- ! 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 )
- !------------------------------------------------------------------------------!
-
+
+!------------------------------------------------------------------------------!
+! 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 chem_3d_data_averaging( mode, variable )
 
@@ -638,10 +635,10 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  mode    !< 
-    CHARACTER (LEN=*) :: variable !< 
-
-    LOGICAL      ::  match_def !< flag indicating natural-type surface
-    LOGICAL      ::  match_lsm !< flag indicating natural-type surface
-    LOGICAL      ::  match_usm !< flag indicating urban-type surface
+    CHARACTER (LEN=*) ::  mode     !< 
+    CHARACTER (LEN=*) ::  variable !< 
+
+    LOGICAL ::  match_def  !< flag indicating default-type surface
+    LOGICAL ::  match_lsm  !< flag indicating natural-type surface
+    LOGICAL ::  match_usm  !< flag indicating urban-type surface
 
     INTEGER(iwp) ::  i                  !< grid index x direction
@@ -649,96 +646,99 @@
     INTEGER(iwp) ::  k                  !< grid index z direction
     INTEGER(iwp) ::  m                  !< running index surface type
-    INTEGER(iwp) :: lsp                 !< running index for chem spcs
+    INTEGER(iwp) ::  lsp               !< running index for chem spcs
 
     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
-             chem_species(lsp)%conc_av = 0.0_wp
-          ENDIF
-       ENDDO
-
-    ELSEIF ( mode == 'sum' )  THEN
-
-       DO lsp = 1, nspec
-          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)
+       IF ( mode == 'allocate' )  THEN
+
+          DO  lsp = 1, nspec
+             IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
+                  TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
+                chem_species(lsp)%conc_av = 0.0_wp
+             ENDIF
+          ENDDO
+
+       ELSEIF ( mode == 'sum' )  THEN
+
+          DO  lsp = 1, nspec
+             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)
+                      ENDDO
                    ENDDO
                 ENDDO
-             ENDDO
-          ELSEIF ( TRIM(variable(4:)) == TRIM('cssws*') )        THEN
-             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)
-
-                   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)
-                   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)
-                   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)
-                   ENDIF
-                ENDDO
-             ENDDO
-          ENDIF
-       ENDDO
-
-    ELSEIF ( mode == 'average' )  THEN
-       DO lsp = 1, nspec
-          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 )
+             ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) )  THEN
+                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)
+
+                      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)
+                      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)
+                      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)
+                      ENDIF
                    ENDDO
                 ENDDO
-             ENDDO
-
-          ELSEIF (TRIM(variable(4:)) == TRIM('cssws*') )            THEN
-             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 )
+             ENDIF
+          ENDDO
+
+       ELSEIF ( mode == 'average' )  THEN
+
+          DO  lsp = 1, nspec
+             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 )
+                      ENDDO
+                   ENDDO
                 ENDDO
-             ENDDO
-             CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )                 
-          ENDIF
-       ENDDO
+
+             ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) )  THEN
+                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 )
+                   ENDDO
+                ENDDO
+                CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )                 
+             ENDIF
+          ENDDO
+       ENDIF
 
     ENDIF
 
-    ENDIF
-
  END SUBROUTINE chem_3d_data_averaging
 
-!    
-!------------------------------------------------------------------------------!
-!
+    
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine to initialize and set all boundary conditions for chemical species
 !------------------------------------------------------------------------------!
-
  SUBROUTINE chem_boundary_conds( mode )                                           
 
@@ -755,5 +755,5 @@
         ONLY:  bc_h                                                           
 
-    CHARACTER (len=*), INTENT(IN) ::  mode
+    CHARACTER (LEN=*), INTENT(IN) ::  mode
     INTEGER(iwp) ::  i                            !< grid index x direction.
     INTEGER(iwp) ::  j                            !< grid index y direction.
@@ -768,5 +768,5 @@
        CASE ( 'init' )
 
-          IF ( bc_cs_b == 'dirichlet' )    THEN
+          IF ( bc_cs_b == 'dirichlet' )  THEN
              ibc_cs_b = 0 
           ELSEIF ( bc_cs_b == 'neumann' )  THEN
@@ -774,16 +774,16 @@
           ELSE
              message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"'  
-             CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 )     !<
+             CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 )
           ENDIF                                                                  
 !
-!--          Set Integer flags and check for possible erroneous settings for top
-!--          boundary condition. 
-          IF ( bc_cs_t == 'dirichlet' )             THEN
+!--       Set Integer flags and check for possible erroneous settings for top
+!--       boundary condition. 
+          IF ( bc_cs_t == 'dirichlet' )  THEN
              ibc_cs_t = 0 
-          ELSEIF ( bc_cs_t == 'neumann' )           THEN
+          ELSEIF ( bc_cs_t == 'neumann' )  THEN
              ibc_cs_t = 1
           ELSEIF ( bc_cs_t == 'initial_gradient' )  THEN
              ibc_cs_t = 2
-          ELSEIF ( bc_cs_t == 'nested' )            THEN          
+          ELSEIF ( bc_cs_t == 'nested' )  THEN          
              ibc_cs_t = 3
           ELSE
@@ -792,10 +792,11 @@
           ENDIF
 
-
        CASE ( 'set_bc_bottomtop' )                   
+!
 !--          Bottom boundary condtions for chemical species      
-          DO lsp = 1, nspec                                                      
-             IF ( ibc_cs_b == 0 ) THEN                   
-                DO l = 0, 1 
+          DO  lsp = 1, nspec                                                      
+             IF ( ibc_cs_b == 0 )  THEN                   
+                DO  l = 0, 1 
+!
 !--                Set index kb: For upward-facing surfaces (l=0), kb=-1, i.e.
 !--                the chem_species(nspec)%conc_p value at the topography top (k-1)
@@ -805,15 +806,16 @@
                    kb = MERGE( -1, 1, l == 0 )
                    !$OMP PARALLEL DO PRIVATE( i, j, k )
-                   DO m = 1, bc_h(l)%ns
-                      i = bc_h(l)%i(m)            
-                      j = bc_h(l)%j(m)
-                      k = bc_h(l)%k(m)
+                   DO  m = 1, bc_h(l)%ns
+                       i = bc_h(l)%i(m)            
+                       j = bc_h(l)%j(m)
+                       k = bc_h(l)%k(m)
                       chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i) 
                    ENDDO                                        
                 ENDDO                                       
 
-             ELSEIF ( ibc_cs_b == 1 ) THEN
+             ELSEIF ( ibc_cs_b == 1 )  THEN
+!
 !--             in boundary_conds there is som extra loop over m here for passive tracer
-                DO l = 0, 1
+                DO  l = 0, 1
                    kb = MERGE( -1, 1, l == 0 )
                    !$OMP PARALLEL DO PRIVATE( i, j, k )                                           
@@ -828,21 +830,22 @@
              ENDIF
        ENDDO    ! end lsp loop  
-
+!
 !--    Top boundary conditions for chemical species - Should this not be done for all species?
           IF ( ibc_cs_t == 0 )  THEN                    
-             DO lsp = 1, nspec
+             DO  lsp = 1, nspec
                 chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)        
              ENDDO
           ELSEIF ( ibc_cs_t == 1 )  THEN
-             DO lsp = 1, nspec
+             DO  lsp = 1, nspec
                 chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
              ENDDO
           ELSEIF ( ibc_cs_t == 2 )  THEN
-             DO lsp = 1, nspec
+             DO  lsp = 1, nspec
                 chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)
              ENDDO
           ENDIF
-!
+
        CASE ( 'set_bc_lateral' )                       
+!
 !--             Lateral boundary conditions for chem species at inflow boundary
 !--             are automatically set when chem_species concentration is
@@ -854,17 +857,17 @@
 
           IF ( bc_radiation_s )  THEN
-             DO lsp = 1, nspec
+             DO  lsp = 1, nspec
                 chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
              ENDDO
           ELSEIF ( bc_radiation_n )  THEN
-             DO lsp = 1, nspec
+             DO  lsp = 1, nspec
                 chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
              ENDDO
           ELSEIF ( bc_radiation_l )  THEN
-             DO lsp = 1, nspec
+             DO  lsp = 1, nspec
                 chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
              ENDDO
           ELSEIF ( bc_radiation_r )  THEN
-             DO lsp = 1, nspec
+             DO  lsp = 1, nspec
                 chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
              ENDDO
@@ -875,5 +878,5 @@
  END SUBROUTINE chem_boundary_conds
 
-!
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -881,6 +884,12 @@
 !> Boundary conditions for prognostic variables in chemistry: decycling in the 
 !> x-direction
-!-----------------------------------------------------------------------------
+!------------------------------------------------------------------------------!
  SUBROUTINE chem_boundary_conds_decycle( cs_3d, cs_pr_init )
+!
+!-- Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not
+!-- approproate for chemicals compounds since they may accumulate too much.
+!-- If no proper boundary conditions from a DYNAMIC input file are available,
+!-- de-cycling applies the initial profiles at the inflow boundaries for
+!-- Dirichlet boundary conditions
 
     IMPLICIT NONE 
@@ -897,6 +906,5 @@
 
     flag = 0.0_wp
-
-
+!
 !-- Left and right boundaries
     IF ( decycle_chem_lr  .AND.  bc_lr_cyc )  THEN
@@ -962,6 +970,5 @@
        ENDDO
     ENDIF
-
-
+!
 !-- South and north boundaries
     IF ( decycle_chem_ns  .AND.  bc_ns_cyc )  THEN
@@ -1029,14 +1036,12 @@
     ENDIF
  END SUBROUTINE chem_boundary_conds_decycle
-   !
-!------------------------------------------------------------------------------!
-!
+
+
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine for checking data output for chemical species
 !------------------------------------------------------------------------------!
-
  SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
-
 
     IMPLICIT NONE
@@ -1050,5 +1055,5 @@
     INTEGER(iwp) ::  k
 
-    CHARACTER(len=16)    ::  spec_name
+    CHARACTER(LEN=16)    ::  spec_name
 
 !
@@ -1058,16 +1063,16 @@
     unit = 'illegal'
 
-    spec_name = TRIM(var(4:))             !< var 1:3 is 'kc_' or 'em_'.
-
-    IF ( TRIM(var(1:3)) == 'em_' )  THEN
-       DO lsp=1,nspec
-          IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
+    spec_name = TRIM( var(4:) )             !< var 1:3 is 'kc_' or 'em_'.
+
+    IF ( TRIM( var(1:3) ) == 'em_' )  THEN
+       DO  lsp=1,nspec
+          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)
-          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 micrograms per m**3 for PM10 and PM25 (PM2.5)
+          IF (spec_name(1:2) == 'PM')  THEN
              unit = 'kg m-2 s-1'
           ENDIF
@@ -1076,6 +1081,6 @@
     ELSE
 
-       DO lsp=1,nspec
-          IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
+       DO  lsp=1,nspec
+          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
              unit = 'ppm'
           ENDIF
@@ -1084,11 +1089,11 @@
 !--            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  
+          IF (spec_name(1:2) == 'PM')  THEN  
             unit = 'kg m-3'
           ENDIF
        ENDDO
 
-       DO lsp=1,nphot
-          IF (TRIM(spec_name) == TRIM(phot_frequen(lsp)%name))   THEN
+       DO  lsp=1,nphot
+          IF (TRIM( spec_name ) == TRIM( phot_frequen(lsp)%name ) )  THEN
              unit = 'sec-1'
           ENDIF
@@ -1099,12 +1104,11 @@
     RETURN
  END SUBROUTINE chem_check_data_output
-!
-!------------------------------------------------------------------------------!
-!
+
+
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine for checking data output of profiles for chemistry model
 !------------------------------------------------------------------------------!
-
  SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
 
@@ -1122,10 +1126,10 @@
     CHARACTER (LEN=*) ::  variable !< 
     CHARACTER (LEN=*) ::  dopr_unit
-    CHARACTER(len=16) ::  spec_name
+    CHARACTER (LEN=16) ::  spec_name
 
     INTEGER(iwp) ::  var_count, lsp  !<
 
 
-    spec_name = TRIM(variable(4:))    
+    spec_name = TRIM( variable(4:) )    
 
        IF (  .NOT.  air_chemistry )  THEN
@@ -1136,11 +1140,11 @@
 
        ELSE
-          DO lsp = 1, nspec
-             IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN 
+          DO  lsp = 1, nspec
+             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_unit  = 'ppm'
-                IF (spec_name(1:2) == 'PM')   THEN
+                IF (spec_name(1:2) == 'PM')  THEN
                      dopr_unit = 'kg m-3'
                 ENDIF
@@ -1153,5 +1157,5 @@
  END SUBROUTINE chem_check_data_output_pr
 
-!
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -1163,5 +1167,5 @@
     IMPLICIT NONE
 
-    LOGICAL       ::  found
+    LOGICAL  ::  found
     INTEGER (iwp) ::  lsp_usr      !< running index for user defined chem spcs 
     INTEGER (iwp) ::  lsp          !< running index for chem spcs.
@@ -1171,18 +1175,18 @@
        message_string = 'Chemical reactions: ON'
        CALL message( 'chem_check_parameters', 'CM0421', 0, 0, 0, 6, 0 )
-    ELSEIF ( .not. (chem_gasphase_on) ) THEN
+    ELSEIF ( .NOT. (chem_gasphase_on) )  THEN
        message_string = 'Chemical reactions: OFF'
        CALL message( 'chem_check_parameters', 'CM0422', 0, 0, 0, 6, 0 )
     ENDIF
-
+!
 !-- check for chemistry time-step
     IF ( call_chem_at_all_substeps )  THEN
        message_string = 'Chemistry is calculated at all meteorology time-step'
        CALL message( 'chem_check_parameters', 'CM0423', 0, 0, 0, 6, 0 )
-    ELSEIF ( .not. (call_chem_at_all_substeps) ) THEN
+    ELSEIF ( .not. (call_chem_at_all_substeps) )  THEN
        message_string = 'Sub-time-steps are skipped for chemistry time-steps'
        CALL message( 'chem_check_parameters', 'CM0424', 0, 0, 0, 6, 0 )
     ENDIF
-
+!
 !-- check for photolysis scheme 
     IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant')  )  THEN
@@ -1190,26 +1194,25 @@
        CALL message( 'chem_check_parameters', 'CM0425', 1, 2, 0, 6, 0 )
     ENDIF
-
+!
 !-- check for  decycling of chem species
-    IF ( (.not. any(decycle_method == 'neumann') ) .AND. (.not. any(decycle_method == 'dirichlet') ) )   THEN
+    IF ( (.NOT. any(decycle_method == 'neumann') ) .AND. (.NOT. any(decycle_method == 'dirichlet') ) )  THEN
        message_string = 'Incorrect boundary conditions. Only neumann or '   &
                 // 'dirichlet &available for decycling chemical species '
        CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 )
     ENDIF
+!
 !-- check for chemical mechanism used 
     CALL get_mechanism_name
-    IF (chem_mechanism /= trim(cs_mech) )  THEN
+    IF ( chem_mechanism /= TRIM( cs_mech ) )  THEN
        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
-
-
-!---------------------
+!
 !-- chem_check_parameters is called before the array chem_species is allocated!
 !-- temporary switch of this part of the check
 !    RETURN                !bK commented
+
     CALL chem_init_internal
-!---------------------
-
+!
 !-- check for initial chem species input
     lsp_usr = 1
@@ -1217,32 +1220,32 @@
     DO WHILE ( cs_name (lsp_usr) /= 'novalue')
        found = .FALSE.
-       DO lsp = 1, nvar
-          IF ( TRIM(cs_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
+       DO  lsp = 1, nvar
+          IF ( TRIM( cs_name (lsp_usr) ) == TRIM( chem_species(lsp)%name) )  THEN
              found = .TRUE.
              EXIT
           ENDIF 
        ENDDO
-       IF ( .not. found ) THEN
-          message_string = 'Unused/incorrect input for initial surface value: ' // TRIM(cs_name(lsp_usr))
+       IF ( .NOT.  found )  THEN
+          message_string = 'Unused/incorrect input for initial surface value: ' //     &
+                            TRIM( cs_name(lsp_usr) )
           CALL message( 'chem_check_parameters', 'CM0427', 1, 2, 0, 6, 0 )
        ENDIF
        lsp_usr = lsp_usr + 1
     ENDDO
-
+!
 !-- check for surface  emission flux chem species
-
     lsp_usr = 1
     lsp     = 1
     DO WHILE ( surface_csflux_name (lsp_usr) /= 'novalue')
        found = .FALSE.
-       DO lsp = 1, nvar
-          IF ( TRIM(surface_csflux_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
+       DO  lsp = 1, nvar
+          IF ( TRIM( surface_csflux_name (lsp_usr) ) == TRIM( chem_species(lsp)%name ) )  THEN
              found = .TRUE.
              EXIT
           ENDIF 
        ENDDO
-       IF ( .not. found ) THEN
+       IF ( .NOT.  found )  THEN
           message_string = 'Unused/incorrect input of chemical species for surface emission fluxes: '  &
-                            // TRIM(surface_csflux_name(lsp_usr))
+                            // TRIM( surface_csflux_name(lsp_usr) )
           CALL message( 'chem_check_parameters', 'CM0428', 1, 2, 0, 6, 0 )
        ENDIF
@@ -1252,7 +1255,6 @@
  END SUBROUTINE chem_check_parameters
 
-!
-!------------------------------------------------------------------------------!
-!
+
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
@@ -1260,5 +1262,4 @@
 !> @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 )
@@ -1281,7 +1282,7 @@
     REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf  
 
-    !
-    !-- local variables.
-    CHARACTER(len=16)    ::  spec_name
+!
+!-- local variables.
+    CHARACTER(LEN=16)    ::  spec_name
     INTEGER(iwp) ::  lsp
     INTEGER(iwp) ::  i               !< grid index along x-direction
@@ -1289,5 +1290,4 @@
     INTEGER(iwp) ::  k               !< grid index along z-direction
     INTEGER(iwp) ::  char_len        !< length of a character string
-
 !
 !-- Next statement is to avoid compiler warnings about unused variables
@@ -1295,18 +1295,18 @@
 
     found = .FALSE.
-    char_len  = LEN_TRIM(variable)
+    char_len  = LEN_TRIM( variable )
 
     spec_name = TRIM( variable(4:char_len-3) )
 
-       DO lsp=1,nspec
-          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             
+       DO  lsp=1,nspec
+          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             
 !
 !--   todo: remove or replace by "CALL message" mechanism (kanani)
-!                    IF(myid == 0)  WRITE(6,*) 'Output of species ' // TRIM(variable)  //               &
-!                                                             TRIM(chem_species(lsp)%name)       
-             IF (av == 0) THEN
+!                    IF(myid == 0)  WRITE(6,*) 'Output of species ' // TRIM( variable )  //       &
+!                                                             TRIM( chem_species(lsp)%name )       
+             IF (av == 0)  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -1316,6 +1316,4 @@
                                               REAL( fill_value, KIND = wp ),                       &
                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-
-
                       ENDDO
                    ENDDO
@@ -1334,6 +1332,6 @@
                 ENDDO
              ENDIF
-              grid = 'zu'           
-              found = .TRUE.
+             grid = 'zu'           
+             found = .TRUE.
           ENDIF
        ENDDO
@@ -1343,14 +1341,11 @@
  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 )
-
 
     USE indices
@@ -1372,8 +1367,7 @@
     REAL(wp)             ::  fill_value   !<
     REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf 
-
-
-    !-- local variables
-    CHARACTER(len=16)    ::  spec_name
+!
+!-- local variables
+    CHARACTER(LEN=16)    ::  spec_name
     INTEGER(iwp)         ::  i
     INTEGER(iwp)         ::  j
@@ -1385,47 +1379,47 @@
 
     found = .FALSE.
-    IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) ) THEN
+    IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) )  THEN
        RETURN
     ENDIF
 
-    spec_name = TRIM(variable(4:))
-
-    IF ( variable(1:3) == 'em_' ) THEN
-
-      local_pf = 0.0_wp
-
-      DO lsp = 1, nvar   !!! cssws - nvar species, chem_species - nspec species !!!
-         IF ( TRIM(spec_name) == TRIM(chem_species(lsp)%name) )   THEN
-            !
-            !-- 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)
-            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)) = &
+    spec_name = TRIM( variable(4:) )
+
+    IF ( variable(1:3) == 'em_' )  THEN
+
+       local_pf = 0.0_wp
+
+       DO  lsp = 1, nvar   !!! cssws - nvar species, chem_species - nspec species !!!
+          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
+!
+!--          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)
+             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)
-            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)) = &
+             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)
-               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)
-               ENDDO
-            ENDDO
-            found = .TRUE.
-         ENDIF
-      ENDDO
+                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)
+                ENDDO
+             ENDDO
+             found = .TRUE.
+          ENDIF
+       ENDDO
     ELSE
-      DO lsp=1,nspec
-         IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
+      DO  lsp = 1, nspec
+         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)       
-            IF (av == 0) THEN
+!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
+!                                                           TRIM( chem_species(lsp)%name )       
+            IF (av == 0)  THEN
                DO  i = nxl, nxr
                   DO  j = nys, nyn
@@ -1440,4 +1434,5 @@
 
             ELSE
+
                DO  i = nxl, nxr
                   DO  j = nys, nyn
@@ -1459,12 +1454,11 @@
 
  END SUBROUTINE chem_data_output_3d
-!
-!------------------------------------------------------------------------------!
-!
+
+
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine defining mask output variables for chemical species
 !------------------------------------------------------------------------------!
-    
  SUBROUTINE chem_data_output_mask( av, variable, found, local_pf )
 
@@ -1477,15 +1471,13 @@
     IMPLICIT NONE
 
-    CHARACTER(LEN=5) ::  grid        !< flag to distinquish between staggered grids
-
-    CHARACTER (LEN=*)::  variable    !< 
-    INTEGER(iwp)     ::  av          !< flag to control data output of instantaneous or time-averaged data
-    LOGICAL          ::  found       !< 
+    CHARACTER(LEN=5)  ::  grid        !< flag to distinquish between staggered grids
+    CHARACTER(LEN=*)  ::  variable    !< 
+    INTEGER(iwp)  ::  av              !< flag to control data output of instantaneous or time-averaged data
+    LOGICAL ::  found
     REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
               local_pf   !< 
-
-
-    !-- local variables.
-    CHARACTER(len=16)    ::  spec_name
+!
+!-- local variables.
+    CHARACTER(LEN=16)  ::  spec_name
     INTEGER(iwp) ::  lsp
     INTEGER(iwp) ::  i               !< grid index along x-direction
@@ -1495,15 +1487,15 @@
 
     found = .TRUE.
-    grid  = 's'
+    grid = 's'
 
     spec_name = TRIM( variable(4:) )
 
-    DO lsp=1,nspec
-       IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name) )               THEN             
+    DO  lsp=1,nspec
+       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)       
-          IF (av == 0) THEN
+!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
+!                                                        TRIM( chem_species(lsp)%name )       
+          IF (av == 0)  THEN
              IF ( .NOT. mask_surface(mid) )  THEN
 
@@ -1582,5 +1574,4 @@
              ENDIF
 
-
           ENDIF
           found = .FALSE.
@@ -1592,7 +1583,6 @@
  END SUBROUTINE chem_data_output_mask      
 
-!
-!------------------------------------------------------------------------------!
-!
+
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
@@ -1612,5 +1602,5 @@
     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'
@@ -1625,7 +1615,7 @@
 
  END SUBROUTINE chem_define_netcdf_grid
-!
-!------------------------------------------------------------------------------!
-!
+
+
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
@@ -1637,15 +1627,15 @@
 
     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
+    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
     CALL get_mechanism_name
-
+!
 !-- Write chemistry model  header
     WRITE( io, 1 )
-
+!
 !-- Gasphase reaction status
     IF ( chem_gasphase_on )  THEN 
@@ -1654,9 +1644,8 @@
        WRITE( io, 3 )
     ENDIF
-
 !
 !-- Chemistry time-step
     WRITE ( io, 4 ) cs_time_step
-
+!
 !-- Emission mode info
     IF ( mode_emis == "DEFAULT" )  THEN 
@@ -1667,12 +1656,12 @@
        WRITE( io, 7 )
     ENDIF 
-
+!
 !-- Photolysis scheme info
-    IF ( photolysis_scheme == "simple" )      THEN
+    IF ( photolysis_scheme == "simple" )  THEN
        WRITE( io, 8 ) 
-    ELSEIF (photolysis_scheme == "constant" ) THEN
+    ELSEIF (photolysis_scheme == "constant" )  THEN
        WRITE( io, 9 )
     ENDIF
-
+!
 !-- Emission flux info
     lsp = 1
@@ -1681,6 +1670,6 @@
        docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ','  
        IF ( LEN_TRIM( docsflux_chr ) >= 75 )  THEN
-        WRITE ( io, 10 )  docsflux_chr
-        docsflux_chr = '       '
+          WRITE ( io, 10 )  docsflux_chr
+          docsflux_chr = '       '
        ENDIF
        lsp = lsp + 1
@@ -1690,6 +1679,5 @@
        WRITE ( io, 10 )  docsflux_chr
     ENDIF
-
-
+!
 !-- initializatoin of Surface and profile chemical species
 
@@ -1708,5 +1696,5 @@
        WRITE ( io, 11 )  docsinit_chr
     ENDIF
-
+!
 !-- number of variable and fix chemical species and number of reactions
     cs_fixed = nspec - nvar
@@ -1734,7 +1722,6 @@
  END SUBROUTINE chem_header
 
-!
-!------------------------------------------------------------------------------!
-!
+
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
@@ -1742,12 +1729,11 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE chem_init_arrays
-
+!
 !-- Please use this place to allocate required arrays
 
  END SUBROUTINE chem_init_arrays
 
-!
-!------------------------------------------------------------------------------!
-!
+
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
@@ -1767,5 +1753,4 @@
     INTEGER(iwp) ::  j !< running index y dimension
     INTEGER(iwp) ::  n !< running index for chemical species
-
 !
 !-- Next statement is to avoid compiler warning about unused variables
@@ -1774,6 +1759,4 @@
            ilu_urban ) == 0 )  CONTINUE
          
-
-
     IF ( emissions_anthropogenic )  CALL chem_emissions_init
 !
@@ -1796,7 +1779,6 @@
  END SUBROUTINE chem_init
 
-!
-!------------------------------------------------------------------------------!
-!
+
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
@@ -1813,5 +1795,4 @@
 
     IMPLICIT NONE
-
 !
 !-- Local variables
@@ -1821,5 +1802,5 @@
     INTEGER(iwp) ::  lpr_lev           !< running index for chem spcs profile level
 
-    IF ( emissions_anthropogenic ) THEN
+    IF ( emissions_anthropogenic )  THEN
        CALL netcdf_data_input_chemistry_data( chem_emis_att, chem_emis )
     ENDIF
@@ -1874,8 +1855,8 @@
 
     ENDDO
-
 !
 !-- Initial concentration of profiles is prescribed by parameters cs_profile
 !-- and cs_heights in the namelist &chemistry_parameters
+
     CALL chem_init_profiles
 !    
@@ -1892,13 +1873,9 @@
        ENDDO
     ENDIF
-
-
-
 !
 !-- Initialize model variables
     IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
          TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
-
-
+!
 !--    First model run of a possible job queue.
 !--    Initial profiles of the variables must be computed.
@@ -1907,5 +1884,5 @@
 !
 !--       Transfer initial profiles to the arrays of the 3D model
-          DO lsp = 1, nspec
+          DO  lsp = 1, nspec
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
@@ -1921,5 +1898,5 @@
           CALL location_message( 'initializing with constant chemistry profiles', .FALSE. )
 
-          DO lsp = 1, nspec 
+          DO  lsp = 1, nspec 
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
@@ -1930,9 +1907,8 @@
 
        ENDIF
-
 !
 !--    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 
              chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) +  &
                                                cs_surface_initial_change(lsp)
@@ -1941,5 +1917,5 @@
 !
 !--    Initiale old and new time levels. 
-       DO lsp = 1, nvar
+       DO  lsp = 1, nvar
           chem_species(lsp)%tconc_m = 0.0_wp                      
           chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
@@ -1948,15 +1924,12 @@
     ENDIF
 
-
-
-!--- new code add above this line
-    DO lsp = 1, nphot
+    DO  lsp = 1, nphot
        phot_frequen(lsp)%name = phot_names(lsp)
 !
 !-- todo: remove or replace by "CALL message" mechanism (kanani)
-!         IF( myid == 0 )  THEN
-!            WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,TRIM(phot_names(lsp))
-!         ENDIF 
-          phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => freq_1(:,:,:,lsp)
+!--       IF( myid == 0 )  THEN
+!--          WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,TRIM( phot_names(lsp) )
+!--       ENDIF 
+          phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  =>  freq_1(:,:,:,lsp)
     ENDDO
 
@@ -1967,7 +1940,6 @@
  END SUBROUTINE chem_init_internal
 
-!
-!------------------------------------------------------------------------------!
-!
+
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
@@ -1976,32 +1948,34 @@
 !> analogue to parameters u_profile, v_profile and uv_heights)
 !------------------------------------------------------------------------------!
- SUBROUTINE chem_init_profiles              !< SUBROUTINE is called from chem_init in case of 
-   !< TRIM( initializing_actions ) /= 'read_restart_data'
-   !< We still need to see what has to be done in case of restart run
+ SUBROUTINE chem_init_profiles              
+!
+!-- SUBROUTINE is called from chem_init in case of TRIM( initializing_actions ) /= 'read_restart_data'
+!< We still need to see what has to be done in case of restart run
+
     USE chem_modules
 
     IMPLICIT NONE
-
-    !-- Local variables
+!
+!-- Local variables
     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 ::  npr_lev    !< the next available profile lev
-
-    !-----------------
-    !-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles
-    !-- and heights. If parameter "cs_profile" is not prescribed then initial surface values
-    !-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and
-    !-- "cs_heights" are prescribed, their values will!override the constant profile given by
-    !-- "cs_surface".
-    ! 
+!
+!-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles
+!-- and heights. If parameter "cs_profile" is not prescribed then initial surface values
+!-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and
+!-- "cs_heights" are prescribed, their values will!override the constant profile given by
+!-- "cs_surface".
     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
+!
+!--          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
+                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
                    DO lpr_lev = 0, nzt+1
                       chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
@@ -2016,5 +1990,5 @@
 
                    npr_lev = 1
-                   !                     chem_species(lsp)%conc_pr_init(0) = 0.0_wp
+!                   chem_species(lsp)%conc_pr_init(0) = 0.0_wp
                    DO  lpr_lev = 1, nz+1
                       IF ( npr_lev < 100 )  THEN
@@ -2028,5 +2002,5 @@
                          ENDDO
                       ENDIF
-                      IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr, npr_lev + 1) /= 9999999.9_wp )  THEN
+                      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) ) /                          &
@@ -2038,7 +2012,8 @@
                    ENDDO
                 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
           ENDDO
@@ -2048,12 +2023,11 @@
 
  END SUBROUTINE chem_init_profiles
- !
- !------------------------------------------------------------------------------!
- !
- ! Description:
- ! ------------
- !> Subroutine to integrate chemical species in the given chemical mechanism
- !------------------------------------------------------------------------------!
-
+
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine to integrate chemical species in the given chemical mechanism
+!------------------------------------------------------------------------------!
  SUBROUTINE chem_integrate_ij( i, j )
 
@@ -2068,6 +2042,6 @@
     INTEGER,INTENT(IN)       :: i
     INTEGER,INTENT(IN)       :: j
-
-    !--   local variables
+!
+!--   local variables
     INTEGER(iwp) ::  lsp                                                     !< running index for chem spcs.
     INTEGER(iwp) ::  lph                                                     !< running index for photolysis frequencies
@@ -2096,22 +2070,21 @@
     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
-    IF (chem_gasphase_on) THEN
+    IF (chem_gasphase_on)  THEN
        nacc = 0
        nrej = 0
 
        tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt)
-       !
-       !--       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)
        tmp_fact_i = 1.0_wp/tmp_fact
 
-       IF ( humidity ) THEN
+       IF ( humidity )  THEN
           IF ( bulk_cloud_model )  THEN
              tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) *  &
@@ -2124,5 +2097,5 @@
        ENDIF
 
-       DO lsp = 1,nspec
+       DO  lsp = 1,nspec
           tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:) 
        ENDDO
@@ -2131,11 +2104,6 @@
           tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)               
        ENDDO
-       !
-       !--   todo: remove (kanani)
-       !           IF(myid == 0 .AND. chem_debug0 ) THEN
-       !              IF (i == 10 .AND. j == 10) WRITE(0,*) 'begin chemics step ',dt_3d
-       !           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)
@@ -2146,5 +2114,5 @@
        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
@@ -2159,5 +2127,5 @@
             icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus )
 
-       DO lsp = 1,nspec
+       DO  lsp = 1,nspec
           chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:)
        ENDDO
@@ -2168,11 +2136,11 @@
     RETURN
  END SUBROUTINE chem_integrate_ij
- !
- !------------------------------------------------------------------------------!
- !
- ! Description:
- ! ------------
- !> Subroutine defining parin for &chemistry_parameters for chemistry model
- !------------------------------------------------------------------------------!
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine defining parin for &chemistry_parameters for chemistry model
+!------------------------------------------------------------------------------!
  SUBROUTINE chem_parin
 
@@ -2230,7 +2198,7 @@
          surface_csflux_name,              &
          time_fac_type
-
-    !-- analogue 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
+!
+!-- analogue 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)
@@ -2240,7 +2208,6 @@
     !>  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
@@ -2252,6 +2219,6 @@
     atol = 1.0_wp
     rtol = 0.01_wp
-    !
-    !--   Try to find chemistry package
+!
+!--   Try to find chemistry package
     REWIND ( 11 )
     line = ' '
@@ -2260,9 +2227,9 @@
     ENDDO
     BACKSPACE ( 11 )
-    !
-    !--   Read chemistry namelist
+!
+!--   Read chemistry namelist
     READ ( 11, chemistry_parameters, ERR = 10, END = 20 )      
-    !
-    !--   Enable chemistry model
+!
+!--   Enable chemistry model
     air_chemistry = .TRUE.                    
     GOTO 20
@@ -2273,7 +2240,6 @@
 
  20 CONTINUE
-
-    !
-    !--    check for  emission mode for chem species
+!
+!--    check for emission mode for chem species
     IF ( (mode_emis /= 'PARAMETERIZED')  .AND. ( mode_emis /= 'DEFAULT' ) .AND. ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
        message_string = 'Incorrect mode_emiss  option select. Please check spelling'
@@ -2282,12 +2248,12 @@
 
     t_steps = my_steps          
-
-    !--    Determine the number of user-defined profiles and append them to the
-    !--    standard data output (data_output_pr)
+!
+!--    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
 
     DO  WHILE ( data_output_pr(i)  /= ' '  .AND.  i <= 100 )
-       IF ( TRIM( data_output_pr(i)(1:3)) == 'kc_' ) THEN
+       IF ( TRIM( data_output_pr(i)(1:3) ) == 'kc_' )  THEN
           max_pr_cs_tmp = max_pr_cs_tmp+1
        ENDIF
@@ -2295,5 +2261,5 @@
     ENDDO
 
-    IF ( max_pr_cs_tmp > 0 ) THEN
+    IF ( max_pr_cs_tmp > 0 )  THEN
        cs_pr_namelist_found = .TRUE.
        max_pr_cs = max_pr_cs_tmp
@@ -2301,17 +2267,17 @@
 
     !      Set Solver Type
-    IF(icntrl(3) == 0)   THEN
+    IF(icntrl(3) == 0)  THEN
        solver_type = 'rodas3'           !Default
-    ELSE IF(icntrl(3) == 1)   THEN
+    ELSE IF(icntrl(3) == 1)  THEN
        solver_type = 'ros2'
-    ELSE IF(icntrl(3) == 2)   THEN
+    ELSE IF(icntrl(3) == 2)  THEN
        solver_type = 'ros3'
-    ELSE IF(icntrl(3) == 3)   THEN
+    ELSE IF(icntrl(3) == 3)  THEN
        solver_type = 'ro4'
-    ELSE IF(icntrl(3) == 4)   THEN
+    ELSE IF(icntrl(3) == 4)  THEN
        solver_type = 'rodas3'
-    ELSE IF(icntrl(3) == 5)   THEN
+    ELSE IF(icntrl(3) == 5)  THEN
        solver_type = 'rodas4'
-    ELSE IF(icntrl(3) == 6)   THEN
+    ELSE IF(icntrl(3) == 6)  THEN
        solver_type = 'Rang3'
     ELSE
@@ -2320,26 +2286,26 @@
     END IF
 
-    !
-    !--   todo: remove or replace by "CALL message" mechanism (kanani)
-    !       write(text,*) 'gas_phase chemistry: solver_type = ',TRIM(solver_type)
-    !kk    Has to be changed to right calling sequence
-    !kk       CALL location_message( TRIM(text), .FALSE. )
-    !        IF(myid == 0)   THEN
-    !           write(9,*) ' '
-    !           write(9,*) 'kpp setup '
-    !           write(9,*) ' '
-    !           write(9,*) '    gas_phase chemistry: solver_type = ',TRIM(solver_type)
-    !           write(9,*) ' '
-    !           write(9,*) '    Hstart  = ',rcntrl(3)
-    !           write(9,*) '    FacMin  = ',rcntrl(4)
-    !           write(9,*) '    FacMax  = ',rcntrl(5)
-    !           write(9,*) ' '
-    !           IF(vl_dim > 1)    THEN
-    !              write(9,*) '    Vector mode                   vektor length = ',vl_dim
-    !           ELSE
-    !              write(9,*) '    Scalar mode'
-    !           ENDIF
-    !           write(9,*) ' '
-    !        END IF
+!
+!--   todo: remove or replace by "CALL message" mechanism (kanani)
+!       write(text,*) 'gas_phase chemistry: solver_type = ',TRIM( solver_type )
+!kk    Has to be changed to right calling sequence
+!kk       CALL location_message( TRIM( text ), .FALSE. )
+!        IF(myid == 0)  THEN
+!           write(9,*) ' '
+!           write(9,*) 'kpp setup '
+!           write(9,*) ' '
+!           write(9,*) '    gas_phase chemistry: solver_type = ',TRIM( solver_type )
+!           write(9,*) ' '
+!           write(9,*) '    Hstart  = ',rcntrl(3)
+!           write(9,*) '    FacMin  = ',rcntrl(4)
+!           write(9,*) '    FacMax  = ',rcntrl(5)
+!           write(9,*) ' '
+!           IF(vl_dim > 1)  THEN
+!              write(9,*) '    Vector mode                   vektor length = ',vl_dim
+!           ELSE
+!              write(9,*) '    Scalar mode'
+!           ENDIF
+!           write(9,*) ' '
+!        END IF
 
     RETURN
@@ -2347,14 +2313,13 @@
  END SUBROUTINE chem_parin
 
- !
- !------------------------------------------------------------------------------!
- !
- ! Description:
- ! ------------
- !> Subroutine calculating prognostic equations for chemical species 
- !> (vector-optimized).
- !> Routine is called separately for each chemical species over a loop from
- !> prognostic_equations.
- !------------------------------------------------------------------------------!
+ 
+!------------------------------------------------------------------------------!
+! 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 chem_prognostic_equations( cs_scalar_p, cs_scalar, tcs_scalar_m,   &
       pr_init_cs, ilsp )
@@ -2390,9 +2355,9 @@
 
 
-    !
-    !-- Tendency terms for chemical species
+!
+!-- Tendency terms for chemical species
     tend = 0.0_wp
-    !    
-    !-- Advection terms
+!    
+!-- Advection terms
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( ws_scheme_sca )  THEN
@@ -2404,6 +2369,6 @@
        CALL advec_s_up( cs_scalar )
     ENDIF
-    !
-    !-- Diffusion terms  (the last three arguments are zero)
+!
+!-- Diffusion terms  (the last three arguments are zero)
     CALL diffusion_s( cs_scalar,                                               &
          surf_def_h(0)%cssws(ilsp,:),                             &
@@ -2424,6 +2389,6 @@
          surf_usm_v(2)%cssws(ilsp,:),                             &
          surf_usm_v(3)%cssws(ilsp,:) )
-    !    
-    !-- Prognostic equation for chemical species
+!    
+!-- Prognostic equation for chemical species
     DO  i = nxl, nxr
        DO  j = nys, nyn    
@@ -2443,6 +2408,6 @@
        ENDDO
     ENDDO
-    !
-    !-- Calculate tendencies for the next Runge-Kutta step
+!
+!-- Calculate tendencies for the next Runge-Kutta step
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( intermediate_timestep_count == 1 )  THEN
@@ -2469,15 +2434,15 @@
  END SUBROUTINE chem_prognostic_equations
 
- !------------------------------------------------------------------------------!
- !
- ! 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 chem_prognostic_equations_ij( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs,     &
-      i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs,                  &
+
+!------------------------------------------------------------------------------!
+! 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 chem_prognostic_equations_ij( cs_scalar_p, cs_scalar, tcs_scalar_m,                 &
+      pr_init_cs, i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs,                            &
       flux_l_cs, diss_l_cs )
     USE pegrid         
@@ -2507,13 +2472,13 @@
     REAL(wp), DIMENSION(0:nz+1)                                  :: pr_init_cs  !<
 
-    !
-    !-- local variables
+!
+!-- local variables
 
     INTEGER :: k
-    !
-    !--    Tendency-terms for chem spcs.
+!
+!--    Tendency-terms for chem spcs.
     tend(:,j,i) = 0.0_wp
-    !    
-    !-- Advection terms
+!    
+!-- Advection terms
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( ws_scheme_sca )  THEN
@@ -2526,8 +2491,6 @@
        CALL advec_s_up( i, j, cs_scalar )
     ENDIF
-
-    !
-
-    !-- Diffusion terms (the last three arguments are zero)
+!
+!-- Diffusion terms (the last three arguments are zero)
 
     CALL diffusion_s( i, j, cs_scalar,                                                 &
@@ -2541,8 +2504,7 @@
          surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:),        &
          surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
-
-    !    
-    !-- Prognostic equation for chem spcs
-    DO k = nzb+1, nzt
+!    
+!-- Prognostic equation for chem spcs
+    DO  k = nzb+1, nzt
        cs_scalar_p(k,j,i) = cs_scalar(k,j,i) + ( dt_3d  *                       &
             ( tsc(2) * tend(k,j,i) +         &
@@ -2557,7 +2519,6 @@
        IF ( cs_scalar_p(k,j,i) < 0.0_wp )  cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i)    !FKS6
     ENDDO
-
-    !
-    !-- Calculate tendencies for the next Runge-Kutta step
+!
+!-- Calculate tendencies for the next Runge-Kutta step
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( intermediate_timestep_count == 1 )  THEN
@@ -2576,12 +2537,10 @@
  END SUBROUTINE chem_prognostic_equations_ij
 
- !
- !------------------------------------------------------------------------------!
- !
- ! Description:
- ! ------------
- !> Subroutine to read restart data of chemical species
- !------------------------------------------------------------------------------!
-
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine to read restart data of chemical species
+!------------------------------------------------------------------------------!
  SUBROUTINE chem_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,             &
                             nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc,   &
@@ -2623,15 +2582,15 @@
     IF ( ALLOCATED(chem_species) )  THEN
 
-       DO lsp = 1, nspec
+       DO  lsp = 1, nspec
 
           !< for time-averaged chemical conc.
-          spc_name_av  =  TRIM(chem_species(lsp)%name)//'_av'
-
-          IF ( restart_string(1:length) == TRIM(chem_species(lsp)%name) )    &
-               THEN
+          spc_name_av  =  TRIM( chem_species(lsp)%name )//'_av'
+
+          IF ( restart_string(1:length) == TRIM( chem_species(lsp)%name) )    &
+             THEN
              !< read data into tmp_3d
              IF ( k == 1 )  READ ( 13 )  tmp_3d  
              !< fill ..%conc in the restart run    
-             chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,                  &
+             chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,                    &
                   nxlc-nbgp:nxrc+nbgp) =                  & 
                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
@@ -2639,5 +2598,5 @@
           ELSEIF (restart_string(1:length) == spc_name_av )  THEN
              IF ( k == 1 )  READ ( 13 )  tmp_3d
-             chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,               &
+             chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,                 &
                   nxlc-nbgp:nxrc+nbgp) =               &
                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
@@ -2651,12 +2610,13 @@
 
  END SUBROUTINE chem_rrd_local
- !
- !-------------------------------------------------------------------------------!
- !> Description:
- !> Calculation of horizontally averaged profiles
- !> This routine is called for every statistic region (sr) defined by the user,
- !> but at least for the region "total domain" (sr=0).
- !> quantities.
- !------------------------------------------------------------------------------!
+
+
+!-------------------------------------------------------------------------------!
+!> Description:
+!> Calculation of horizontally averaged profiles
+!> This routine is called for every statistic region (sr) defined by the user,
+!> but at least for the region "total domain" (sr=0).
+!> quantities.
+!-------------------------------------------------------------------------------!
  SUBROUTINE chem_statistics( mode, sr, tn )
 
@@ -2681,13 +2641,14 @@
     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.
-       !--    hom(:,:,:,:) =  dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
-       !                       w*pt*), dim-4 = statistical region.
-
-       !$OMP DO
+!
+!--    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.
+
+!$OMP DO
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -2710,12 +2671,12 @@
  END SUBROUTINE chem_statistics
 
- !
- !------------------------------------------------------------------------------!
- !
- ! Description:
- ! ------------
- !> Subroutine for swapping of timelevels for chemical species
- !> called out from subroutine swap_timelevel
- !------------------------------------------------------------------------------!
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine for swapping of timelevels for chemical species
+!> called out from subroutine swap_timelevel
+!------------------------------------------------------------------------------!
+
 
  SUBROUTINE chem_swap_timelevel( level )
@@ -2724,16 +2685,16 @@
 
     INTEGER(iwp), INTENT(IN) ::  level
-    !
-    !-- local variables
+!
+!-- local variables
     INTEGER(iwp)             ::  lsp
 
 
-    IF ( level == 0 ) THEN
-       DO lsp=1, nvar                                       
+    IF ( level == 0 )  THEN
+       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)
@@ -2743,12 +2704,13 @@
     RETURN
  END SUBROUTINE chem_swap_timelevel
- !
- !------------------------------------------------------------------------------!
- !
- ! Description:
- ! ------------
- !> Subroutine to write restart data for chemistry model
- !------------------------------------------------------------------------------!
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine to write restart data for chemistry model
+!------------------------------------------------------------------------------!
  SUBROUTINE chem_wrd_local
+
 
     IMPLICIT NONE 
@@ -2766,20 +2728,17 @@
 
 
-
- !-------------------------------------------------------------------------------!
- !
- ! 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    
- !> @TODO: Consider overlaying horizontal surfaces
- !> @TODO: Consider resolved vegetation    
- !> @TODO: Check error messages
- !-------------------------------------------------------------------------------!
-
+!-------------------------------------------------------------------------------!
+! 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    
+!> @TODO: Consider overlaying horizontal surfaces
+!> @TODO: Consider resolved vegetation    
+!> @TODO: Check error messages
+!-------------------------------------------------------------------------------!
  SUBROUTINE chem_depo( i, j )
 
@@ -2801,13 +2760,11 @@
 
 
-
     IMPLICIT NONE
 
     INTEGER(iwp), INTENT(IN) ::  i
     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) ::  lu                            !< running index for landuse classes
+!    INTEGER(iwp) ::  lu                           !< running index for landuse classes
     INTEGER(iwp) ::  lu_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
@@ -2826,7 +2783,6 @@
     INTEGER(iwp) ::  pspec                         !< running index
     INTEGER(iwp) ::  i_pspec                       !< index for matching depac gas component
-
-    !
-    !-- Vegetation                                              !<Match to DEPAC
+!
+!-- Vegetation                                              !<Match to DEPAC
     INTEGER(iwp) ::  ind_lu_user = 0                         !< --> ERROR  
     INTEGER(iwp) ::  ind_lu_b_soil = 1                       !< --> ilu_desert
@@ -2848,7 +2804,6 @@
     INTEGER(iwp) ::  ind_lu_mixed_forest = 17                !< --> ilu_coniferous_forest (ave(decid+conif))
     INTEGER(iwp) ::  ind_lu_intrup_forest = 18               !< --> ilu_other (ave(other+decid))
-
-    !
-    !-- Water
+!
+!-- Water
     INTEGER(iwp) ::  ind_luw_user = 0                        !< --> ERROR
     INTEGER(iwp) ::  ind_luw_lake = 1                        !< --> ilu_water_inland
@@ -2857,7 +2812,6 @@
     INTEGER(iwp) ::  ind_luw_pond = 4                        !< --> ilu_water_inland
     INTEGER(iwp) ::  ind_luw_fountain = 5                    !< --> ilu_water_inland 
-
-    !
-    !-- Pavement
+!
+!-- Pavement
     INTEGER(iwp) ::  ind_lup_user = 0                        !< --> ERROR
     INTEGER(iwp) ::  ind_lup_asph_conc = 1                   !< --> ilu_desert
@@ -2876,8 +2830,6 @@
     INTEGER(iwp) ::  ind_lup_art_turf = 14                   !< --> ilu_desert
     INTEGER(iwp) ::  ind_lup_clay = 15                       !< --> ilu_desert
-
-
-    !
-    !-- Particle parameters according to the respective aerosol classes (PM25, PM10)
+!
+!-- Particle parameters according to the respective aerosol classes (PM25, PM10)
     INTEGER(iwp) ::  ind_p_size = 1     !< index for partsize in particle_pars
     INTEGER(iwp) ::  ind_p_dens = 2     !< index for rhopart in particle_pars
@@ -2929,9 +2881,8 @@
     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)
+!
+!-- 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
@@ -2939,12 +2890,9 @@
          /), (/ 3, 2 /) )
 
-
     LOGICAL ::  match_lsm     !< flag indicating natural-type surface
     LOGICAL ::  match_usm     !< flag indicating urban-type surface
-
-
-    !
-    !-- List of names of possible tracers
-    CHARACTER(len=*), PARAMETER ::  pspecnames(nposp) = (/ &
+!
+!-- List of names of possible tracers
+    CHARACTER(LEN=*), PARAMETER ::  pspecnames(nposp) = (/ &
          'NO2           ', &    !< NO2
          'NO            ', &    !< NO
@@ -3016,8 +2964,6 @@
          'tpm25_biascorr', &    !< tpm25_biascorr
          'O3_biascorr   ' /)    !< o3_biascorr
-
-
-    !
-    !-- tracer mole mass:
+!
+!-- tracer mole mass:
     REAL(wp), PARAMETER ::  specmolm(nposp) = (/ &
          xm_O * 2 + xm_N, &                         !< NO2
@@ -3090,29 +3036,23 @@
          xm_dummy, &                                !< tpm25_biascorr
          xm_O * 3 /)                                !< o3_biascorr
-
-
-    !
-    !-- Initialize surface element m
+!
+!-- Initialize surface element m
     m = 0
     k = 0
-    !
-    !-- LSM or USM surface present at i,j:
-    !-- Default surfaces are NOT considered for deposition
+!
+!-- LSM or USM surface present at i,j:
+!-- Default surfaces are NOT considered for deposition
     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)
-
-
-    !
-    !--For LSM surfaces
+!
+!--For LSM surfaces
 
     IF ( match_lsm )  THEN
-
-       !
-       !-- Get surface element information at i,j:
+!
+!--    Get surface element information at i,j:
        m = surf_lsm_h%start_index(j,i)
        k = surf_lsm_h%k(m)
-
-       !
-       !-- Get needed variables for surface element m
+!
+!--    Get needed variables for surface element m
        ustar_lu  = surf_lsm_h%us(m)
        z0h_lu    = surf_lsm_h%z0h(m)
@@ -3121,13 +3061,11 @@
        lai = surf_lsm_h%lai(m)
        sai = lai + 1
-
-       !
-       !-- For small grid spacing neglect R_a
+!
+!--    For small grid spacing neglect R_a
        IF ( dzw(k) <= 1.0 )  THEN
           r_aero_lu = 0.0_wp
        ENDIF
-
-       !
-       !--Initialize lu's
+!
+!--    Initialize lu's
        lu_palm = 0
        lu_dep = 0
@@ -3136,14 +3074,11 @@
        luw_palm = 0
        luw_dep = 0
-
-       ! 
-       !--Initialize budgets
+! 
+!--    Initialize budgets
        bud_now_lu  = 0.0_wp
        bud_now_lup = 0.0_wp
        bud_now_luw = 0.0_wp
-
-
-       !
-       !-- Get land use for i,j and assign to DEPAC lu
+!
+!--    Get land use for i,j and assign to DEPAC lu
        IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
           lu_palm = surf_lsm_h%vegetation_type(m)
@@ -3245,8 +3180,6 @@
           ENDIF
        ENDIF
-
-
-       !
-       !-- Set wetness indicator to dry or wet for lsm vegetation or pavement
+!
+!--    Set wetness indicator to dry or wet for lsm vegetation or pavement
        IF ( surf_lsm_h%c_liq(m) > 0 )  THEN
           nwet = 1
@@ -3254,7 +3187,6 @@
           nwet = 0
        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)
@@ -3263,53 +3195,40 @@
        ENDIF
 
-
        dh = dzw(k)
        inv_dh = 1.0_wp / dh
        dt_dh = dt_chem/dh
-
-       !
-       !-- Concentration at i,j,k
+!
+!--    Concentration at i,j,k
        DO  lsp = 1, nspec
           cc(lsp) = chem_species(lsp)%conc(k,j,i)
        ENDDO
 
-
-       !
-       !-- Temperature at i,j,k
+!--    Temperature at i,j,k
        ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
 
        ts       = ttemp - 273.15  !< in degrees celcius
-
-       !
-       !-- Viscosity of air
+!
+!--    Viscosity of air
        visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
-
-       !
-       !-- Air density at k
+!
+!--    Air density at k
        dens = rho_air_zw(k)
-
-       !
-       !-- Calculate relative humidity from specific humidity for DEPAC
+!
+!--    Calculate relative humidity from specific humidity for DEPAC
        qv_tmp = MAX(q(k,j,i),0.0_wp)
        relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, 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.
-
-       !
-       !-- Vegetation
+!
+!-- 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.
+!
+!--    Vegetation
        IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
 
-
           slinnfac = 1.0_wp
-
-          !
-          !-- Get deposition velocity vd
+!
+!--       Get deposition velocity vd
           DO  lsp = 1, nvar
-             !
-             !-- Initialize 
+!
+!--          Initialize 
              vs = 0.0_wp
              vd_lu = 0.0_wp
@@ -3319,6 +3238,6 @@
              IF ( spc_names(lsp) == 'PM10' )  THEN
                 part_type = 1
-                !
-                !-- Sedimentation velocity
+!
+!--             Sedimentation velocity
                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
                      particle_pars(ind_p_size, part_type), &
@@ -3340,11 +3259,10 @@
              ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
                 part_type = 2
-                !
-                !-- Sedimentation velocity
+!
+!--             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, &
@@ -3359,10 +3277,7 @@
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
-
              ELSE !< GASES
-
-                !
-                !-- Read spc_name of current species for gas parameter
-
+!
+!--             Read spc_name of current species for gas parameter
                 IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
                    i_pspec = 0
@@ -3374,31 +3289,28 @@
 
                 ELSE
-                   !
-                   !-- For now species not deposited
+!
+!--             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 k:
+!
+!--             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 k:
 
                 ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
-
-                !
-                !-- Atmospheric concentration in DEPAC is requested in ug/m3:
+!
+!--             Atmospheric concentration in DEPAC is requested in ug/m3:
                 !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
                 catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
-
-                !
-                !-- Diffusivity for DEPAC relevant gases
-                !-- Use default value 
+!
+!--             Diffusivity for DEPAC relevant gases
+!--             Use default value 
                 diffc            = 0.11e-4
-                !
-                !-- overwrite with known coefficients of diffusivity from Massman (1998)
+!
+!--             overwrite with known coefficients of diffusivity from Massman (1998)
                 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
                 IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
@@ -3408,21 +3320,16 @@
                 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
                 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
-
-
-                !
-                !-- Get quasi-laminar boundary layer resistance Rb:
+!
+!--             Get quasi-laminar boundary layer resistance Rb:
                 CALL get_rb_cell( (lu_dep == ilu_water_sea) .OR. (lu_dep == ilu_water_inland), &
                      z0h_lu, ustar_lu, diffc, &
                      Rb )
-
-                !
-                !-- Get Rc_tot
+!
+!--             Get Rc_tot
                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
                      relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
                      r_aero_lu , Rb )
-
-
-                !
-                !-- Calculate budget
+!
+!--             Calculate budget
                 IF ( Rc_tot <= 0.0 )  THEN
 
@@ -3440,23 +3347,18 @@
           ENDDO
        ENDIF
-
-
-       !
-       !-- Pavement
+!
+!--    Pavement
        IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 )  THEN
-
-
-          !
-          !-- No vegetation on pavements:
+!
+!--       No vegetation on pavements:
           lai = 0.0_wp
           sai = 0.0_wp
 
           slinnfac = 1.0_wp
-
-          !
-          !-- Get vd
+!
+!--       Get vd
           DO  lsp = 1, nvar
-             !
-             !-- Initialize
+!
+!--       Initialize
              vs = 0.0_wp
              vd_lu = 0.0_wp
@@ -3466,6 +3368,6 @@
              IF ( spc_names(lsp) == 'PM10' )  THEN
                 part_type = 1
-                !
-                !-- Sedimentation velocity:
+!
+!--             Sedimentation velocity:
                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
                      particle_pars(ind_p_size, part_type), &
@@ -3487,11 +3389,10 @@
              ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
                 part_type = 2
-                !
-                !-- Sedimentation velocity:
+!
+!--             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, &
@@ -3506,9 +3407,7 @@
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
-
              ELSE  !<GASES
-
-                !
-                !-- Read spc_name of current species for gas parameter
+!
+!--             Read spc_name of current species for gas parameter
 
                 IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
@@ -3521,30 +3420,28 @@
 
                 ELSE
-                   !
-                   !-- For now species not deposited
+!
+!--                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:
+!
+!--             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:
 
                 ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
-
-                !
-                !-- Atmospheric concentration in DEPAC is requested in ug/m3:
+!
+!--             Atmospheric concentration in DEPAC is requested in ug/m3:
                 !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
                 catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
-
-                !
-                !-- Diffusivity for DEPAC relevant gases
-                !-- Use default value 
+!
+!--             Diffusivity for DEPAC relevant gases
+!--             Use default value 
                 diffc            = 0.11e-4
-                !
-                !-- overwrite with known coefficients of diffusivity from Massman (1998)
+!
+!--             overwrite with known coefficients of diffusivity from Massman (1998)
                 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
                 IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
@@ -3554,57 +3451,40 @@
                 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
                 IF ( spc_names(lsp) == 'NH3' ) diffc = 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_lu, ustar_lu, diffc, &
-                     Rb )
-
-
-                !
-                !-- Get Rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
-                     relh, lai, sai, nwet, lup_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
-                     r_aero_lu , Rb )
-
-
-                !
-                !-- Calculate budget
+!
+!--             Get quasi-laminar boundary layer resistance Rb:
+                CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland),   &
+                     z0h_lu, ustar_lu, diffc, Rb )
+!
+!--             Get Rc_tot
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,      &
+                                         glrad, cos_zenith, relh, lai, sai, nwet, lup_dep, 2,      &
+                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,             &
+                                         r_aero_lu , Rb )
+!
+!--             Calculate budget
                 IF ( Rc_tot <= 0.0 )  THEN
-
                    bud_now_lup(lsp) = 0.0_wp
-
                 ELSE
-
                    vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
-
                    bud_now_lup(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
                         (1.0_wp - exp(-vd_lu * dt_dh )) * dh
                 ENDIF
 
-
              ENDIF
           ENDDO
        ENDIF
-
-
-       !
-       !-- Water
+!
+!--    Water
        IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 )  THEN
-
-
-          !
-          !-- No vegetation on water:
+!
+!--       No vegetation on water:
           lai = 0.0_wp
           sai = 0.0_wp
-
           slinnfac = 1.0_wp
-
-          !
-          !-- Get vd
+!
+!--       Get vd
           DO  lsp = 1, nvar
-             !
-             !-- Initialize
+!
+!--          Initialize
              vs = 0.0_wp
              vd_lu = 0.0_wp
@@ -3614,6 +3494,6 @@
              IF ( spc_names(lsp) == 'PM10' )  THEN
                 part_type = 1
-                !
-                !-- Sedimentation velocity:
+!
+!--             Sedimentation velocity:
                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
                      particle_pars(ind_p_size, part_type), &
@@ -3632,14 +3512,12 @@
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
-
              ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
                 part_type = 2
-                !
-                !-- Sedimentation velocity:
+!
+!--             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, &
@@ -3654,9 +3532,7 @@
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
-
              ELSE  !<GASES
-
-                !
-                !-- Read spc_name of current species for gas PARAMETER
+!
+!--             Read spc_name of current species for gas PARAMETER
 
                 IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
@@ -3669,30 +3545,28 @@
 
                 ELSE
-                   !
-                   !-- For now species not deposited
+!
+!--                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:
-
-                ppm_to_ugm3 =  (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
-                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
-
-                !
-                !-- Diffusivity for DEPAC relevant gases
-                !-- Use default value 
+!
+!--             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:
+
+                ppm_to_ugm3 = (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
+                catm = cc(lsp) * ppm_to_ugm3 *  specmolm(i_pspec)  ! in ug/m3
+!
+!--             Diffusivity for DEPAC relevant gases
+!--             Use default value 
                 diffc            = 0.11e-4
-                !
-                !-- overwrite with known coefficients of diffusivity from Massman (1998)
+!
+!--             overwrite with known coefficients of diffusivity from Massman (1998)
                 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
                 IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
@@ -3702,20 +3576,16 @@
                 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
                 IF ( spc_names(lsp) == 'NH3' ) diffc = 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_lu, ustar_lu, diffc, &
-                     Rb )
-
-                !
-                !-- Get Rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
-                     relh, lai, sai, nwet, luw_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
-                     r_aero_lu , Rb )
-
-
-                ! Calculate budget 
+!
+!--             Get quasi-laminar boundary layer resistance Rb:
+                CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland),  &
+                     z0h_lu, ustar_lu, diffc, rb )
+
+!--             Get Rc_tot
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,      &
+                                         glrad, cos_zenith, relh, lai, sai, nwet, luw_dep, 2,      &
+                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,            &
+                                          r_aero_lu , Rb )
+!
+!--             Calculate budget 
                 IF ( Rc_tot <= 0.0 )  THEN
 
@@ -3736,15 +3606,13 @@
 
        bud_now = 0.0_wp
-       !
-       !-- Calculate overall budget for surface m and adapt concentration
+!
+!--    Calculate overall budget for surface m and adapt concentration
        DO  lsp = 1, nspec
-
 
           bud_now(lsp) = surf_lsm_h%frac(ind_veg_wall,m) * bud_now_lu(lsp) + &
                surf_lsm_h%frac(ind_pav_green,m) * bud_now_lup(lsp) + &
                surf_lsm_h%frac(ind_wat_win,m) * bud_now_luw(lsp)
-
-          !
-          !-- Compute new concentration:
+!
+!--       Compute new concentration:
           cc(lsp) = cc(lsp) + bud_now(lsp) * inv_dh
 
@@ -3754,20 +3622,14 @@
 
     ENDIF
-
-
-
-
-    !
-    !-- For USM surfaces    
+!
+!-- For USM surfaces    
 
     IF ( match_usm )  THEN
-
-       !
-       !-- Get surface element information at i,j:
+!
+!--    Get surface element information at i,j:
        m = surf_usm_h%start_index(j,i)
        k = surf_usm_h%k(m)
-
-       !
-       !-- Get needed variables for surface element m
+!
+!--    Get needed variables for surface element m
        ustar_lu  = surf_usm_h%us(m)
        z0h_lu    = surf_usm_h%z0h(m)
@@ -3776,13 +3638,11 @@
        lai = surf_usm_h%lai(m)
        sai = lai + 1
-
-       !
-       !-- For small grid spacing neglect R_a
+!
+!--    For small grid spacing neglect R_a
        IF ( dzw(k) <= 1.0 )  THEN
           r_aero_lu = 0.0_wp
        ENDIF
-
-       !
-       !-- Initialize lu's
+!
+!--    Initialize lu's
        luu_palm = 0
        luu_dep = 0
@@ -3791,18 +3651,15 @@
        lud_palm = 0
        lud_dep = 0
-
-       !
-       !-- Initialize budgets
+!
+!--    Initialize budgets
        bud_now_luu  = 0.0_wp
        bud_now_lug = 0.0_wp
        bud_now_lud = 0.0_wp
-
-
-       !
-       !-- Get land use for i,j and assign to DEPAC lu
+!
+!--    Get land use for i,j and assign to DEPAC lu
        IF ( surf_usm_h%frac(ind_pav_green,m) > 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_lu_s_grass
           IF ( lug_palm == ind_lu_user )  THEN
@@ -3849,8 +3706,8 @@
 
        IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 )  THEN
-          !
-          !-- For walls in USM assume concrete walls/roofs,
-          !-- assumed LU class desert as also assumed for
-          !-- pavements in LSM
+!
+!--       For walls in USM assume concrete walls/roofs,
+!--       assumed LU class desert as also assumed for
+!--       pavements in LSM
           luu_palm = ind_lup_conc
           IF ( luu_palm == ind_lup_user )  THEN
@@ -3891,8 +3748,8 @@
 
        IF ( surf_usm_h%frac(ind_wat_win,m) > 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
+!
+!--       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
@@ -3931,9 +3788,7 @@
           ENDIF
        ENDIF
-
-
-       !
-       !-- @TODO: Activate these lines as soon as new ebsolver branch is merged:
-       !-- Set wetness indicator to dry or wet for usm vegetation or pavement
+!
+!--    @TODO: Activate these lines as soon as new ebsolver branch is merged:
+!--    Set wetness indicator to dry or wet for usm vegetation or pavement
        !IF ( surf_usm_h%c_liq(m) > 0 )  THEN
        !   nwet = 1
@@ -3941,7 +3796,6 @@
        nwet = 0
        !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)
@@ -3950,57 +3804,44 @@
        ENDIF
 
-
        dh = dzw(k)
        inv_dh = 1.0_wp / dh
        dt_dh = dt_chem/dh
-
-       !
-       !-- Concentration at i,j,k
+!
+!--    Concentration at i,j,k
        DO  lsp = 1, nspec
           cc(lsp) = chem_species(lsp)%conc(k,j,i)
        ENDDO
-
-       !
-       !-- Temperature at i,j,k
+!
+!--    Temperature at i,j,k
        ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
 
        ts       = ttemp - 273.15  !< in degrees celcius
-
-       !
-       !-- Viscosity of air
+!
+!--    Viscosity of air
        visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
-
-       !
-       !-- Air density at k
+!
+!--    Air density at k
        dens = rho_air_zw(k)
-
-       !
-       !-- Calculate relative humidity from specific humidity for DEPAC
+!
+!--    Calculate relative humidity from specific humidity for DEPAC
        qv_tmp = MAX(q(k,j,i),0.0_wp)
        relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(nzb) )
-
-
-
-       !
-       !-- 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
+!
+!--    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(ind_veg_wall,m) > 0 )  THEN
-
-
-          !
-          !-- No vegetation on non-green walls:
+!
+!--       No vegetation on non-green walls:
           lai = 0.0_wp
           sai = 0.0_wp
 
           slinnfac = 1.0_wp
-
-          !
-          !-- Get vd
+!
+!--       Get vd
           DO  lsp = 1, nvar
-             !
-             !-- Initialize 
+!
+!--          Initialize 
              vs = 0.0_wp
              vd_lu = 0.0_wp
@@ -4010,6 +3851,6 @@
              IF (spc_names(lsp) == 'PM10' )  THEN
                 part_type = 1
-                !
-                !-- Sedimentation velocity
+!
+!--             Sedimentation velocity
                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
                      particle_pars(ind_p_size, part_type), &
@@ -4028,14 +3869,12 @@
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
-
              ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
                 part_type = 2
-                !
-                !-- Sedimentation velocity
+!
+!--             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, &
@@ -4050,9 +3889,7 @@
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
-
              ELSE  !< GASES
-
-                !
-                !-- Read spc_name of current species for gas parameter
+!
+!--             Read spc_name of current species for gas parameter
 
                 IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
@@ -4064,31 +3901,29 @@
                    ENDDO
                 ELSE
-                   !
-                   !-- For now species not deposited
+!
+!--                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 k:
+
+                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
 
                 !
-                !-- 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 k:
-
-                ppm_to_ugm3 =  (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
-                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
-
-                !
-                !-- Diffusivity for DEPAC relevant gases
-                !-- Use default value 
+!--             Atmospheric concentration in DEPAC is requested in ug/m3:
+!--                 ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
+!--             catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
+!
+!--             Diffusivity for DEPAC relevant gases
+!--             Use default value 
                 diffc            = 0.11e-4
-                !
-                !-- overwrite with known coefficients of diffusivity from Massman (1998)
+!
+!--             overwrite with known coefficients of diffusivity from Massman (1998)
                 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
                 IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
@@ -4098,21 +3933,17 @@
                 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
                 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
-
-
-                !
-                !-- Get quasi-laminar boundary layer resistance Rb:
-                CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland), &
+!
+!--             Get quasi-laminar boundary layer resistance Rb:
+                CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland),   &
                      z0h_lu, ustar_lu, diffc, &
                      Rb )
-
-                !
-                !-- Get Rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
-                     relh, lai, sai, nwet, luu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
-                     r_aero_lu , Rb )
-
-
-                !
-                !-- Calculate budget
+!
+!--             Get Rc_tot
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,       &
+                                         glrad, cos_zenith, relh, lai, sai, nwet, luu_dep, 2,       &
+                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,             &
+                                         r_aero_lu, rb )
+!
+!--             Calculate budget
                 IF ( Rc_tot <= 0.0 )  THEN
 
@@ -4130,18 +3961,14 @@
           ENDDO
        ENDIF
-
-
-       !
-       !-- Green usm surfaces
+!
+!--    Green usm surfaces
        IF ( surf_usm_h%frac(ind_pav_green,m) > 0 )  THEN
 
-
           slinnfac = 1.0_wp
-
-          !
-          !-- Get vd
+!
+!--       Get vd
           DO  lsp = 1, nvar
-             !
-             !-- Initialize
+!
+!--          Initialize
              vs = 0.0_wp
              vd_lu = 0.0_wp
@@ -4151,5 +3978,6 @@
              IF ( spc_names(lsp) == 'PM10' )  THEN
                 part_type = 1
-                ! Sedimentation velocity
+!
+!--             Sedimentation velocity
                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
                      particle_pars(ind_p_size, part_type), &
@@ -4168,13 +3996,12 @@
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
-
              ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
                 part_type = 2
-                ! Sedimentation velocity
+!
+!--             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, &
@@ -4189,9 +4016,7 @@
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
-
              ELSE  !< GASES
-
-                !
-                !-- Read spc_name of current species for gas parameter
+!
+!--             Read spc_name of current species for gas parameter
 
                 IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
@@ -4203,31 +4028,28 @@
                    ENDDO
                 ELSE
-                   !
-                   !-- For now species not deposited
+!
+!--                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 k:
+!
+!--             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 k:
 
                 ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
-
-                !
-                !-- Atmospheric concentration in DEPAC is requested in ug/m3:
+!
+!--             Atmospheric concentration in DEPAC is requested in ug/m3:
                 !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
                 catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
-
-                !
-                !-- Diffusivity for DEPAC relevant gases
-                !-- Use default value 
+!
+!--             Diffusivity for DEPAC relevant gases
+!--             Use default value 
                 diffc            = 0.11e-4
-                !
-                !-- overwrite with known coefficients of diffusivity from Massman (1998)
+!
+!--             overwrite with known coefficients of diffusivity from Massman (1998)
                 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
                 IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
@@ -4237,22 +4059,17 @@
                 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
                 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
-
-
-                !
-                !-- Get quasi-laminar boundary layer resistance Rb:
-                CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland), &
+!
+!--             Get quasi-laminar boundary layer resistance Rb:
+                CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland),    &
                      z0h_lu, ustar_lu, diffc, &
                      Rb )
-
-
-                !
-                !-- Get Rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
-                     relh, lai, sai, nwet, lug_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
-                     r_aero_lu , Rb )
-
-
-                !
-                !-- Calculate budget
+!
+!--             Get Rc_tot
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,        &
+                                         glrad, cos_zenith, relh, lai, sai, nwet, lug_dep, 2,        &
+                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,              &
+                                         r_aero_lu , Rb )
+!
+!--             Calculate budget
                 IF ( Rc_tot <= 0.0 )  THEN
 
@@ -4267,27 +4084,21 @@
                 ENDIF
 
-
              ENDIF
           ENDDO
        ENDIF
-
-
-       !
-       !-- Windows
+!
+!--    Windows
        IF ( surf_usm_h%frac(ind_wat_win,m) > 0 )  THEN
-
-
-          !
-          !-- No vegetation on windows:
+!
+!--       No vegetation on windows:
           lai = 0.0_wp
           sai = 0.0_wp
 
           slinnfac = 1.0_wp
-
-          !
-          !-- Get vd
+!
+!--       Get vd
           DO  lsp = 1, nvar
-             !
-             !-- Initialize
+!
+!--          Initialize
              vs = 0.0_wp
              vd_lu = 0.0_wp
@@ -4297,6 +4108,6 @@
              IF ( spc_names(lsp) == 'PM10' )  THEN
                 part_type = 1
-                !
-                !-- Sedimentation velocity
+!
+!--             Sedimentation velocity
                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
                      particle_pars(ind_p_size, part_type), &
@@ -4304,6 +4115,23 @@
                      visc)
 
-                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
-                     vs, &
+                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
+                     particle_pars(ind_p_size, part_type), &
+                     particle_pars(ind_p_slip, part_type), &
+                     nwet, ttemp, dens, visc,              &
+                     lud_dep, r_aero_lu, ustar_lu )
+
+                bud_now_lud(lsp) = - cc(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), &
@@ -4315,31 +4143,7 @@
                      (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, ttemp, dens, visc, &
-                     lud_dep, &
-                     r_aero_lu, ustar_lu )
-
-                bud_now_lud(lsp) = - cc(lsp) * &
-                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
-
-
              ELSE  !< GASES
-
-                !
-                !-- Read spc_name of current species for gas PARAMETER
+!
+!--             Read spc_name of current species for gas PARAMETER
 
                 IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
@@ -4351,30 +4155,28 @@
                    ENDDO
                 ELSE
-                   ! For now species not deposited
+!
+!--                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 k:
+!
+!--             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 k:
 
                 ppm_to_ugm3 =  (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
+!
+!--             Atmospheric concentration in DEPAC is requested in ug/m3:
+!--                 ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
                 catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
-
-                !
-                !-- Diffusivity for DEPAC relevant gases
-                !-- Use default value 
-                diffc            = 0.11e-4
-                !
-                !-- overwrite with known coefficients of diffusivity from Massman (1998)
+!
+!--             Diffusivity for DEPAC relevant gases
+!--             Use default value 
+                diffc = 0.11e-4
+!
+!--             overwrite with known coefficients of diffusivity from Massman (1998)
                 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
                 IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
@@ -4384,21 +4186,16 @@
                 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
                 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
-
-
-                !
-                !-- Get quasi-laminar boundary layer resistance Rb:
-                CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland), &
-                     z0h_lu, ustar_lu, diffc, &
-                     Rb )
-
-                !
-                !-- Get Rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
-                     relh, lai, sai, nwet, lud_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
-                     r_aero_lu , Rb )
-
-
-                !
-                !-- Calculate budget 
+!
+!--             Get quasi-laminar boundary layer resistance Rb:
+                CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland),   &
+                     z0h_lu, ustar_lu, diffc, rb )
+!
+!--             Get Rc_tot
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,      &
+                                         glrad, cos_zenith, relh, lai, sai, nwet, lud_dep, 2,      &
+                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,            &
+                                         r_aero_lu , rb )
+!
+!--             Calculate budget 
                 IF ( Rc_tot <= 0.0 )  THEN
 
@@ -4419,6 +4216,6 @@
 
        bud_now = 0.0_wp
-       !
-       !-- Calculate overall budget for surface m and adapt concentration
+!
+!--    Calculate overall budget for surface m and adapt concentration
        DO  lsp = 1, nspec
 
@@ -4427,7 +4224,6 @@
                surf_usm_h%frac(ind_pav_green,m) * bud_now_lug(lsp) + &
                surf_usm_h%frac(ind_wat_win,m) * bud_now_lud(lsp)
-
-          !
-          !-- Compute new concentration
+!
+!--       Compute new concentration
           cc(lsp) = cc(lsp) + bud_now(lsp) * inv_dh
 
@@ -4439,59 +4235,49 @@
 
 
-
  END SUBROUTINE chem_depo
 
 
-
- !----------------------------------------------------------------------------------
- !
- !> DEPAC:
- !> 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.
- !---------------------------------------------------------------------------------
- !    
- !----------------------------------------------------------------------------------    
- !
- !> depac: compute total canopy (or surface) resistance Rc for gases
- !----------------------------------------------------------------------------------
-
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine to compute total canopy (or surface) resistance Rc for gases
+!>
+!> DEPAC:
+!> 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, glrad, sinphi,  &
       rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, catm, diffc,           &
       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])
-   !
-   !-- 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)
-   !
-   !-- 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)
-
+!
+!--   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])
+!--    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)
+!--
+!--    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)
 
     IMPLICIT NONE
 
-    CHARACTER(len=*), INTENT(IN) ::  compnam         !< component name
+    CHARACTER(LEN=*), INTENT(IN) ::  compnam         !< component name
                                                      !< 'HNO3','NO','NO2','O3','SO2','NH3'
     INTEGER(iwp), INTENT(IN) ::  day_of_year         !< day of year, 1 ... 365 (366)
@@ -4518,13 +4304,11 @@
     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:
-    ! 
-    !-- Component number taken from component name, paramteres matched with include files
+!                                                     !< [= 0 for species that don't have a compensation
+!-- Local variables:
+!
+!-- Component number taken from component name, paramteres matched with include files
     INTEGER(iwp) ::  icmp
-
-    !
-    !-- Component numbers:
+!
+!-- Component numbers:
     INTEGER(iwp), PARAMETER ::  icmp_o3   = 1
     INTEGER(iwp), PARAMETER ::  icmp_so2  = 2
@@ -4541,6 +4325,6 @@
     !< = 1 -> constant Rc
     !< = 2 -> temperature dependent Rc
-    !
-    !-- Vegetation indicators:
+!
+!-- Vegetation indicators:
     LOGICAL ::  lai_present                          !< leaves are present for current land use type
     LOGICAL ::  sai_present                          !< vegetation is present for current land use type
@@ -4554,13 +4338,9 @@
     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
     IF ( day_of_year == 0  .OR.  ( catm + lat + ra + rb ) > 0.0_wp )  CONTINUE
-
-
-    !
-    !-- Define component number
+!
+!-- Define component number
     SELECT CASE ( TRIM( compnam ) )
 
@@ -4596,12 +4376,12 @@
 
     CASE default
-       !
-       !-- Component not part of DEPAC --> not deposited
+!
+!--    Component not part of DEPAC --> not deposited
        RETURN
 
     END SELECT
 
-    !
-    !-- Inititalize
+!
+!-- Inititalize
     gw        = 0.0_wp
     gstom     = 0.0_wp
@@ -4611,51 +4391,44 @@
     cstom     = 0.0_wp
     csoil     = 0.0_wp
-
-
-    !
-    !-- Check whether vegetation is present:
+!
+!-- Check whether vegetation is present:
     lai_present = ( lai > 0.0 )
     sai_present = ( sai > 0.0 )
-
-    !
-    !-- Set Rc (i.e. rc_tot) in special cases:
+!
+!-- Set Rc (i.e. rc_tot) in special cases:
     CALL rc_special( icmp, compnam, lu, t, nwet, rc_tot, ready, ccomp_tot )
-
-    !
-    !-- If Rc is not set:
+!
+!-- If Rc is not set:
     IF ( .NOT. ready ) then
-
-       !
-       !-- External conductance:
+!
+!--    External conductance:
        CALL rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai,gw )         
-
-       !
-       !-- Stomatal conductance:
+!
+!--    Stomatal conductance:
        CALL rc_gstom( icmp, compnam, lu, lai_present, lai, glrad, sinphi, t, rh, diffc, gstom, p )
-       !
-       !-- Effective soil conductance:
+!
+!--    Effective soil conductance:
        CALL rc_gsoil_eff( icmp, lu, sai, ust, nwet, t, gsoil_eff )
-
-       !
-       !-- Total canopy conductance (gc_tot) and resistance Rc (rc_tot):
+!
+!--    Total canopy conductance (gc_tot) and resistance Rc (rc_tot):
        CALL rc_rctot( gstom, gsoil_eff, gw, gc_tot, rc_tot )
-
-       !
-       !-- Needed to include compensation point for NH3
-       !-- Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3):
-       !-- CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,&
-       !     lai_present, sai_present, &
-       !     ccomp_tot, &
-       !     catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, &
-       !     c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, &
-       !     tsea=tsea,cw=cw,cstom=cstom,csoil=csoil )
-       !
-       !-- Effective Rc based on compensation points:
-       !   IF ( present(rc_eff) ) then
-       !     check on required arguments:
-       !      IF ( (.not. present(catm)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then
-       !         stop 'output argument rc_eff requires input arguments catm, ra and rb'
-       !      END IF
-       !--    Compute rc_eff :
+!
+!--    Needed to include compensation point for NH3
+!--    Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3):
+!--    CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,&
+!--          lai_present, sai_present, &
+!--          ccomp_tot, &
+!--          catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, &
+!--          c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, &
+!--          tsea=tsea,cw=cw,cstom=cstom,csoil=csoil )
+!
+!--    Effective Rc based on compensation points:
+!--        IF ( present(rc_eff) ) then
+!--          check on required arguments:
+!--           IF ( (.not. present(catm)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then
+!--              stop 'output argument rc_eff requires input arguments catm, ra and rb'
+!--           END IF
+!
+!--       Compute rc_eff :
        !      CALL rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff)
        !   ENDIF
@@ -4665,19 +4438,20 @@
 
 
-
- !-------------------------------------------------------------------
- !> rc_special: compute total canopy resistance in special CASEs
- !-------------------------------------------------------------------
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine to compute total canopy resistance in special cases
+!------------------------------------------------------------------------------!
  SUBROUTINE rc_special( icmp, compnam, lu, t, nwet, rc_tot, ready, ccomp_tot )
 
     IMPLICIT NONE
    
-    CHARACTER(len=*), INTENT(IN) ::  compnam     !< component name
-
-    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)
+    CHARACTER(LEN=*), INTENT(IN)  ::  compnam     !< component name
+
+    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)
@@ -4686,29 +4460,25 @@
     LOGICAL, INTENT(OUT) ::  ready               !< Rc has been set
                                                  !< = 1 -> constant Rc
-                                                 !< = 2 -> temperature dependent Rc
-
 !
 !-- Next line is to avoid compiler warning about unused variable
     IF ( icmp == 0 )  CONTINUE
-
-    !
-    !-- rc_tot is not yet set:
+!
+!-- rc_tot is not yet set:
     ready = .FALSE.
-
-    !
-    !-- Default compensation point in special CASEs = 0:
+!
+!-- Default compensation point in special CASEs = 0:
     ccomp_tot = 0.0_wp
 
     SELECT CASE( TRIM( compnam ) )
     CASE( 'HNO3', 'N2O5', 'NO3', 'H2O2' )
-       !
-       !-- No separate resistances for HNO3; just one total canopy resistance:
+!
+!--    No separate resistances for HNO3; just one total canopy resistance:
        IF ( t < -5.0_wp .AND. nwet == 9 )  THEN
-          !
-          !-- T < 5 C and snow:
+!
+!--       T < 5 C and snow:
           rc_tot = 50.0_wp
        ELSE
-          !
-          !-- all other circumstances:
+!
+!--       all other circumstances:
           rc_tot = 10.0_wp
        ENDIF
@@ -4724,9 +4494,9 @@
        ENDIF
     CASE( 'NO2', 'O3', 'SO2', 'NH3' )
-       !
-       !-- snow surface:
+!
+!--    snow surface:
        IF ( nwet == 9 )  THEN
-          !
-          !-- To be activated when snow is implemented
+!
+!--       To be activated when snow is implemented
           !CALL rc_snow(ipar_snow(icmp),t,rc_tot)
           ready = .TRUE.
@@ -4740,15 +4510,15 @@
 
 
-
- !-------------------------------------------------------------------
- !> rc_gw: compute external conductance
- !-------------------------------------------------------------------
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine to compute external conductance
+!------------------------------------------------------------------------------!
  SUBROUTINE rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai, gw )
 
     IMPLICIT NONE
-
-    !
-    !-- Input/output variables:
-    CHARACTER(len=*), INTENT(IN) ::  compnam      !< component name
+!
+!-- 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
@@ -4757,5 +4527,4 @@
                                                   !< iratns = 2: high NH3/SO2
                                                   !< iratns = 3: very low NH3/SO2
-
     LOGICAL, INTENT(IN) ::  sai_present
 
@@ -4766,5 +4535,4 @@
     REAL(wp), INTENT(OUT) ::  gw                  !< external leaf conductance (m/s)
 
-
     SELECT CASE( TRIM( compnam ) )
 
@@ -4783,11 +4551,10 @@
     CASE( 'NH3' )
        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
 
     CASE default
-       message_string = 'Component '// TRIM(compnam) // ' not supported'
+       message_string = 'Component '// TRIM( compnam ) // ' not supported'
        CALL message( 'rc_gw', 'CM0458', 1, 2, 0, 6, 0 )
     END SELECT
@@ -4796,14 +4563,14 @@
 
 
-
- !-------------------------------------------------------------------
- !> rw_so2: compute external leaf conductance for SO2
- !-------------------------------------------------------------------
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine to compute external leaf conductance for SO2
+!------------------------------------------------------------------------------!
  SUBROUTINE rw_so2( t, nwet, rh, iratns, sai_present, gw )
 
     IMPLICIT NONE
-
-    !
-    !-- Input/output variables:
+!
+!-- 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:
@@ -4817,19 +4584,17 @@
 
     REAL(wp), INTENT(OUT) ::  gw             !< external leaf conductance (m/s)
-
-    !
-    !-- Local variables:
+!
+!-- Local variables:
     REAL(wp) ::  rw                          !< external leaf resistance (s/m)
-
-
-    !
-    !-- Check if vegetation present:
+!
+!-- Check if vegetation present:
     IF ( sai_present )  THEN
 
        IF ( nwet == 0 )  THEN
-          !--------------------------
-          ! dry surface
-          !--------------------------
-          ! T > -1 C
+!
+!--   ------------------------
+!--         dry surface
+!--   ------------------------
+!--         T > -1 C
           IF ( t > -1.0_wp )  THEN
              IF ( rh < 81.3_wp )  THEN
@@ -4848,14 +4613,65 @@
           ENDIF
        ELSE
-          !--------------------------
-          ! wet surface
-          !--------------------------
+!
+!--   ------------------------
+!--         wet surface
+!--   ------------------------
           rw = 10.0_wp !see Table 5, Erisman et al, 1994 Atm. Environment, 0 is impl. as 10
        ENDIF
-
-       ! very low NH3/SO2 ratio:
+!
+!--    very low NH3/SO2 ratio:
        IF ( iratns == iratns_very_low ) rw = rw + 50.0_wp
-
-       ! Conductance:
+!
+!--      Conductance:
+       gw = 1.0_wp / rw
+    ELSE
+!
+!--      no vegetation:
+       gw = 0.0_wp
+    ENDIF
+
+ END SUBROUTINE rw_so2
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine to compute external leaf conductance for NH3,
+!>                  following Sutton & Fowler, 1993
+!------------------------------------------------------------------------------!
+ SUBROUTINE rw_nh3_sutton( tsurf, rh,sai_present, gw )
+
+    IMPLICIT NONE
+!
+!-- Input/output variables:
+    LOGICAL, INTENT(IN) ::  sai_present
+
+    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)
+!
+!-- Local variables:
+    REAL(wp) ::  rw                         !< external leaf resistance (s/m)
+    REAL(wp) ::  sai_grass_haarweg          !< surface area index at experimental site Haarweg
+!
+!-- Fix sai_grass at value valid for Haarweg data for which gamma_w parametrization is derived
+    sai_grass_haarweg = 3.5_wp
+!
+!-- Calculation rw:
+!--                    100 - rh
+!--    rw = 2.0 * exp(----------)
+!--                      12
+
+    IF ( sai_present )  THEN
+!
+!--    External resistance according to Sutton & Fowler, 1993
+       rw = 2.0_wp * exp( ( 100.0_wp - rh ) / 12.0_wp )
+       rw = sai_grass_haarweg * rw
+!
+!--    Frozen soil (from Depac v1):
+       IF ( tsurf < 0.0_wp ) rw = 200.0_wp
+!
+!--    Conductance:
        gw = 1.0_wp / rw
     ELSE
@@ -4864,77 +4680,24 @@
     ENDIF
 
- END SUBROUTINE rw_so2
-
-
-
- !-------------------------------------------------------------------
- !> rw_nh3_sutton: compute external leaf conductance for NH3,
- !>                  following Sutton & Fowler, 1993
- !-------------------------------------------------------------------
- SUBROUTINE rw_nh3_sutton( tsurf, rh,sai_present, gw )
+ END SUBROUTINE rw_nh3_sutton
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine to compute external leaf conductance 
+!------------------------------------------------------------------------------!
+ SUBROUTINE rw_constant( rw_val, sai_present, gw )
 
     IMPLICIT NONE
-
-    !
-    !-- Input/output variables:
+!
+!-- Input/output variables:
     LOGICAL, INTENT(IN) ::  sai_present
 
-    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)
-
-    !
-    !-- Local variables:
-    REAL(wp) ::  rw                         !< external leaf resistance (s/m)
-    REAL(wp) ::  sai_grass_haarweg          !< surface area index at experimental site Haarweg
-
-    !
-    !-- Fix sai_grass at value valid for Haarweg data for which gamma_w parametrization is derived
-    sai_grass_haarweg = 3.5_wp
-
-    !
-    !-- Calculation rw:
-    !                  100 - rh
-    !  rw = 2.0 * exp(----------)
-    !                    12
-
-    IF ( sai_present )  THEN
-
-       ! External resistance according to Sutton & Fowler, 1993
-       rw = 2.0_wp * exp( ( 100.0_wp - rh ) / 12.0_wp )
-       rw = sai_grass_haarweg * rw
-
-       ! Frozen soil (from Depac v1):
-       IF ( tsurf < 0.0_wp ) rw = 200.0_wp
-
-       ! Conductance:
-       gw = 1.0_wp / rw
-    ELSE
-       ! no vegetation:
-       gw = 0.0_wp
-    ENDIF
-
- END SUBROUTINE rw_nh3_sutton
-
-
-
- !-------------------------------------------------------------------
- !> rw_constant: compute constant external leaf conductance
- !-------------------------------------------------------------------
- SUBROUTINE rw_constant( rw_val, sai_present, gw )
-
-    IMPLICIT NONE
-
-    !
-    !-- Input/output variables:
-    LOGICAL, INTENT(IN) ::  sai_present
-
     REAL(wp), INTENT(IN) ::  rw_val       !< constant value of Rw    
 
     REAL(wp), INTENT(OUT) ::  gw          !< wernal leaf conductance (m/s)
-
-    !
-    !-- Compute conductance:
+!
+!-- Compute conductance:
     IF ( sai_present .AND. .NOT.missing(rw_val) )  THEN
        gw = 1.0_wp / rw_val
@@ -4946,15 +4709,15 @@
 
 
-
- !-------------------------------------------------------------------
- !> rc_gstom: compute stomatal conductance
- !-------------------------------------------------------------------
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine to compute stomatal conductance 
+!------------------------------------------------------------------------------!
  SUBROUTINE rc_gstom( icmp, compnam, lu, lai_present, lai, glrad, sinphi, t, rh, diffc, gstom, p ) 
 
     IMPLICIT NONE
-
-    !
-    !-- input/output variables:
-    CHARACTER(len=*), INTENT(IN) ::  compnam       !< component name
+!
+!-- input/output variables:
+    CHARACTER(LEN=*), INTENT(IN) ::  compnam       !< component name
 
     INTEGER(iwp), INTENT(IN) ::  icmp              !< component index
@@ -4973,26 +4736,23 @@
 
     REAL(wp), INTENT(OUT) ::  gstom                !< stomatal conductance (m/s)
-
-    !
-    !-- Local variables
+!
+!-- Local variables
     REAL(wp) ::  vpd                               !< vapour pressure deficit (kPa)
 
     REAL(wp), PARAMETER ::  dO3 = 0.13e-4          !< diffusion coefficient of ozon (m2/s)
-
 !
 !-- Next line is to avoid compiler warning about unused variables
     IF ( icmp == 0 )  CONTINUE
 
-    SELECT CASE( TRIM(compnam) )
+    SELECT CASE( TRIM( compnam ) )
 
     CASE( 'NO', 'CO' )
-       !
-       !-- For no stomatal uptake is neglected:
+!
+!--    For no stomatal uptake is neglected:
        gstom = 0.0_wp
 
     CASE( 'NO2', 'O3', 'SO2', 'NH3' )
-
-       !
-       !-- if vegetation present:
+!
+!--    if vegetation present:
        IF ( lai_present )  THEN
 
@@ -5005,11 +4765,11 @@
           ENDIF
        ELSE
-          !
-          !--no vegetation; zero conductance (infinite resistance):
+!
+!--       no vegetation; zero conductance (infinite resistance):
           gstom = 0.0_wp
        ENDIF
 
     CASE default
-       message_string = 'Component '// TRIM(compnam) // ' not supported'
+       message_string = 'Component '// TRIM( compnam ) // ' not supported'
        CALL message( 'rc_gstom', 'CM0459', 1, 2, 0, 6, 0 )
     END SELECT
@@ -5018,47 +4778,45 @@
 
 
-
- !-------------------------------------------------------------------
- !> rc_gstom_emb: stomatal conductance according to Emberson
- !-------------------------------------------------------------------
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine to compute stomatal conductance according to Emberson
+!------------------------------------------------------------------------------!
  SUBROUTINE rc_gstom_emb( lu, glrad, T, vpd, lai_present, lai, sinp, Gsto, p )
-    !
-    !-- History
-    !>   Original code from Lotos-Euros, TNO, M. Schaap
-    !>   2009-08, M.C. van Zanten, Rivm
-    !>     Updated and extended.
-    !>   2009-09, Arjo Segers, TNO
-    !>     Limitted temperature influence to range to avoid
-    !>     floating point exceptions.
-    !
-    !> Method
-    !
-    !>   Code based on Emberson et al, 2000, Env. Poll., 403-413
-    !>   Notation conform Unified EMEP Model Description Part 1, ch 8
-    !
-    !>   In the calculation of f_light the modification of L. Zhang 2001, AE to the PARshade and PARsun
-    !>   parametrizations of Norman 1982 are applied
-    !
-    !>   f_phen and f_SWP are set to 1
-    !
-    !>   Land use types DEPAC versus Emberson (Table 5.1, EMEP model description)
-    !>   DEPAC                     Emberson
-    !>     1 = grass                 GR = grassland
-    !>     2 = arable land           TC = temperate crops ( lai according to RC = rootcrops)
-    !<     3 = permanent crops       TC = temperate crops ( lai according to RC = rootcrops)
-    !<     4 = coniferous forest     CF = tempareate/boREAL(wp) coniferous forest
-    !>     5 = deciduous forest      DF = temperate/boREAL(wp) deciduous forest
-    !>     6 = water                 W  = water
-    !>     7 = urban                 U  = urban
-    !>     8 = other                 GR = grassland
-    !>     9 = desert                DE = desert
+!
+!>  History
+!>   Original code from Lotos-Euros, TNO, M. Schaap
+!>   2009-08, M.C. van Zanten, Rivm
+!>     Updated and extended.
+!>   2009-09, Arjo Segers, TNO
+!>     Limitted temperature influence to range to avoid
+!>     floating point exceptions.
+
+!> Method
+
+!>   Code based on Emberson et al, 2000, Env. Poll., 403-413
+!>   Notation conform Unified EMEP Model Description Part 1, ch 8
+!
+!>   In the calculation of f_light the modification of L. Zhang 2001, AE to the PARshade and PARsun
+!>   parametrizations of Norman 1982 are applied
+!>   f_phen and f_SWP are set to 1
+!
+!>   Land use types DEPAC versus Emberson (Table 5.1, EMEP model description)
+!>   DEPAC                     Emberson
+!>     1 = grass                 GR = grassland
+!>     2 = arable land           TC = temperate crops ( lai according to RC = rootcrops)
+!>     3 = permanent crops       TC = temperate crops ( lai according to RC = rootcrops)
+!>     4 = coniferous forest     CF = tempareate/boREAL(wp) coniferous forest
+!>     5 = deciduous forest      DF = temperate/boREAL(wp) deciduous forest
+!>     6 = water                 W  = water
+!>     7 = urban                 U  = urban
+!>     8 = other                 GR = grassland
+!>     9 = desert                DE = desert
 
     IMPLICIT NONE
-
-    !
-    !-- Emberson specific declarations
-
-    !
-    !-- Input/output variables:
+!
+!-- Emberson specific declarations
+!
+!-- Input/output variables:
     INTEGER(iwp), INTENT(IN) ::  lu             !< land use type, lu = 1,...,nlu
 
@@ -5075,7 +4833,6 @@
 
     REAL(wp), INTENT(OUT) ::  gsto              !< stomatal conductance (m/s)
-
-    !
-    !-- Local variables:
+!
+!-- Local variables:
     REAL(wp) ::  f_light
     REAL(wp) ::  f_phen
@@ -5094,8 +4851,6 @@
     REAL(wp) ::  pres
     REAL(wp), PARAMETER ::  p_sealevel = 1.01325e05    !< Pa
-
-
-    !
-    !-- Check whether vegetation is present:
+!
+!-- Check whether vegetation is present:
     IF ( lai_present )  THEN
 
@@ -5106,9 +4861,8 @@
           sinphi = sinp
        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 ...
+!
+!--    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 )
@@ -5116,26 +4870,22 @@
           pres = p_sealevel
        ENDIF
-
-       !
-       !-- direct and diffuse par, Photoactive (=visible) radiation:
+!
+!--    direct and diffuse par, Photoactive (=visible) radiation:
        CALL par_dir_diff( glrad, sinphi, pres, p_sealevel, pardir, pardiff )
-
-       !
-       !-- par for shaded leaves (canopy averaged):
+!
+!--    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 ( glrad > 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):
-       !-- alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5
+!
+!--    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 ( glrad > 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 ) )
@@ -5150,7 +4900,6 @@
 
        f_light = MAX(f_light,f_min(lu))
-
-       !
-       !-- temperature influence; only non-zero within range [tmin,tmax]:
+!
+!--    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) )
@@ -5160,23 +4909,23 @@
        END IF
        f_temp = MAX( f_temp, f_min(lu) )
-
-       ! vapour pressure deficit influence
+!
+!--      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
 
@@ -5188,15 +4937,11 @@
 
 
-
  !-------------------------------------------------------------------
  !> 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.
- !
  !>     From a SUBROUTINE obtained from Leiming Zhang,
  !>     Meteorological Service of Canada
- !
  !>     Leiming uses solar irradiance. This should be equal to global radiation and
  !>     Willem Asman set it to global radiation
@@ -5228,48 +4973,37 @@
     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)
+!
+!-- 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
-
-    !
-    !-- Calculate potential visible diffuse radiation
+!
+!-- Calculate potential visible diffuse radiation
     rdv = 0.4_wp * ( 600.0_wp - rdu ) * sinphi
-
-    !
-    !-- Calculate the water absorption in the-near infrared
+!
+!-- 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 )
-
-    !
-    !-- Calculate potential direct beam near-infrared radiation
+!
+!-- 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
-
-    !
-    !-- Calculate potential diffuse near-infrared radiation
+!
+!-- Calculate potential diffuse near-infrared radiation
     rdn = 0.6_wp * ( 720 - rdm - ww ) * sinphi
-
-    !
-    !-- Compute visible and near-infrared radiation
+!
+!-- Compute visible and near-infrared radiation
     rv = MAX( 0.1_wp, rdu + rdv )
     rn = MAX( 0.01_wp, rdm + rdn )
-
-    !
-    !-- Compute ratio between input global radiation and total radiation computed here
+!
+!-- Compute ratio between input global radiation and total radiation computed here
     ratio = MIN( 0.89_wp, glrad / ( rv + rn ) )
-
-    !
-    !-- Calculate total visible radiation
+!
+!-- Calculate total visible radiation
     sv = ratio * rv
-
-    !
-    !-- Calculate fraction of par in the direct beam
+!
+!-- 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
-
-    !
-    !-- Compute direct and diffuse parts of par
+!
+!-- Compute direct and diffuse parts of par
     par_dir = fv * sv
     par_diff = sv - par_dir
@@ -5284,18 +5018,15 @@
 
     IMPLICIT NONE
-
-    !
-    !-- Input/output variables:
+!
+!-- Input/output variables:
     REAL(wp), INTENT(IN) ::  temp    !< temperature (C)
     REAL(wp), INTENT(IN) ::  relh    !< relative humidity (%)
 
     REAL(wp), INTENT(OUT) ::  vpd    !< vapour pressure deficit (kPa)
-
-    !
-    !-- Local variables:
+!
+!-- Local variables:
     REAL(wp) ::  esat
-
-    !
-    !-- fit parameters:
+!
+!-- fit parameters:
     REAL(wp), PARAMETER ::  a1 = 6.113718e-01
     REAL(wp), PARAMETER ::  a2 = 4.43839e-02
@@ -5304,12 +5035,10 @@
     REAL(wp), PARAMETER ::  a5 = 2.16e-07
     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 - relh / 100 )
 
  END SUBROUTINE rc_get_vpd
-
 
 
@@ -5320,7 +5049,6 @@
 
     IMPLICIT NONE
-
-    !
-    !-- Input/output variables:
+!
+!-- Input/output variables:
     INTEGER(iwp), INTENT(IN) ::  icmp          !< component index
     INTEGER(iwp), INTENT(IN) ::  lu            !< land use type, lu = 1,..., nlu
@@ -5329,19 +5057,14 @@
                                                !< N.B. this routine cannot be called with nwet = 9,
                                                !< nwet = 9 should be handled outside this routine.
-
     REAL(wp), INTENT(IN) ::  sai               !< surface area index
     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)
-
-    !
-    !-- local variables:
+!
+!-- local variables:
     REAL(wp) ::  rinc                          !< in canopy resistance  (s/m)
     REAL(wp) ::  rsoil_eff                     !< effective soil resistance (s/m)
-
-
-    !
-    !-- Soil resistance (numbers matched with lu_classes and component numbers)
+!
+!-- 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( (/ &
@@ -5357,23 +5080,18 @@
          -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
+!
+!-- 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./)
-
-
-
-    !
-    !-- Compute in canopy (in crop) resistance:
+!
+!-- Compute in canopy (in crop) resistance:
     CALL rc_rinc( lu, sai, ust, rinc )
-
-    !
-    !-- Check for missing deposition path:
+!
+!-- Check for missing deposition path:
     IF ( missing(rinc) )  THEN
        rsoil_eff = -9999.0_wp
     ELSE
-       !
-       !-- Frozen soil (temperature below 0):
+!
+!--    Frozen soil (temperature below 0):
        IF ( t < 0.0_wp )  THEN
           IF ( missing( rsoil_frozen( icmp ) ) )  THEN
@@ -5383,6 +5101,6 @@
           ENDIF
        ELSE
-          !
-          !-- Non-frozen soil; dry:
+!
+!--       Non-frozen soil; dry:
           IF ( nwet == 0 )  THEN
              IF ( missing( rsoil( lu, icmp ) ) )  THEN
@@ -5391,7 +5109,6 @@
                 rsoil_eff = rsoil( lu, icmp ) + rinc
              ENDIF
-
-             !
-             !-- Non-frozen soil; wet:
+!
+!--       Non-frozen soil; wet:
           ELSEIF ( nwet == 1 )  THEN
              IF ( missing( rsoil_wet( icmp ) ) )  THEN
@@ -5406,7 +5123,6 @@
        ENDIF
     ENDIF
-
-    !
-    !-- Compute conductance:
+!
+!-- Compute conductance:
     IF ( rsoil_eff > 0.0_wp )  THEN
        gsoil_eff = 1.0_wp / rsoil_eff
@@ -5426,6 +5142,6 @@
     IMPLICIT NONE
 
-    !
-    !-- Input/output variables:
+!
+!-- Input/output variables:
     INTEGER(iwp), INTENT(IN) ::  lu          !< land use class, lu = 1, ..., nlu
 
@@ -5434,18 +5150,16 @@
 
     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)
-    !-- h = vegetation height (m)                     gra  ara crop con dec wat   urb   oth   des   ice   sav   trf  wai  med semi
+!
+!-- 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,  &
          14, 14 /)
     REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  h = (/ -999, 1,  1,  20, 20, -999, -999, -999, -999, -999, -999, 20, -999,  &
          1 ,  1 /)
-
-    !
-    !-- Compute Rinc only for arable land, perm. crops, forest; otherwise Rinc = 0:
+!
+!-- 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):
+!       !
+!--    Check for u* > 0 (otherwise denominator = 0):
        IF ( ust > 0.0_wp )  THEN
           rinc = b(lu) * h(lu) * sai/ust
@@ -5464,5 +5178,4 @@
 
 
-
  !-------------------------------------------------------------------
  !> rc_rctot: compute total canopy (or surface) resistance Rc
@@ -5472,6 +5185,6 @@
     IMPLICIT NONE
 
-    !
-    !-- Input/output variables:
+!
+!-- Input/output variables:
     REAL(wp), INTENT(IN) ::  gstom         !< stomatal conductance (s/m)
     REAL(wp), INTENT(IN) ::  gsoil_eff     !< effective soil conductance (s/m)
@@ -5480,11 +5193,9 @@
     REAL(wp), INTENT(OUT) ::  gc_tot       !< total canopy conductance (m/s)
     REAL(wp), INTENT(OUT) ::  rc_tot       !< total canopy resistance Rc (s/m)
-
-    !
-    !-- Total conductance:
+!
+!-- Total conductance:
     gc_tot = gstom + gsoil_eff + gw
-
-    !
-    !-- Total resistance (note: gw can be negative, but no total emission allowed here):
+!
+!-- Total resistance (note: gw can be negative, but no total emission allowed here):
     IF ( gc_tot <= 0.0_wp .OR. gw < 0.0_wp )  THEN
        rc_tot = -9999.0_wp
@@ -5496,10 +5207,8 @@
 
 
-
  !-------------------------------------------------------------------
  !> 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)
@@ -5558,11 +5267,9 @@
  !>
  ! -------------------------------------------------------------------------------------------
-
 ! SUBROUTINE rc_comp_point_rc_eff( ccomp_tot, catm, ra, rb, rc_tot, rc_eff )
 !
 !    IMPLICIT NONE
 !
-!    !
-!    !-- Input/output variables:
+!!-- Input/output variables:
 !    REAL(wp), INTENT(IN) ::  ccomp_tot    !< total compensation point (weighed average of separate compensation points) (ug/m3)
 !    REAL(wp), INTENT(IN) ::  catm         !< atmospheric concentration (ug/m3)
@@ -5574,19 +5281,19 @@
 !
 !    !
-!    !-- Compute effective resistance:
+!!-- Compute effective resistance:
 !    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( catm - ccomp_tot ) < 1.e-8 ) )  THEN
 !       !
-!       !--surface concentration (almost) equal to atmospheric concentration
-!       !-- no exchange between surface and atmosphere, infinite RC --> vd=0
+!!--   surface concentration (almost) equal to atmospheric concentration
+!!--    no exchange between surface and atmosphere, infinite RC --> vd=0
 !       rc_eff = 9999999999.0_wp
 !
 !    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 * catm ) / ( catm - ccomp_tot )
 !
@@ -5602,16 +5309,14 @@
 
 
-
  !-------------------------------------------------------------------
  !> missing: check for data that correspond with a missing deposition path
  !>          this data is represented by -999
  !-------------------------------------------------------------------
-
  LOGICAL function missing( x )
 
     REAL(wp), INTENT(IN) ::  x
 
-    !
-    !-- bandwidth for checking (in)equalities of floats
+!
+!-- bandwidth for checking (in)equalities of floats
     REAL(wp), PARAMETER :: eps = 1.0e-5
 
@@ -5621,12 +5326,10 @@
 
 
-
  ELEMENTAL FUNCTION sedimentation_velocity( rhopart, partsize, slipcor, visc ) RESULT( vs )
 
-
     IMPLICIT NONE
 
-    !
-    !-- in/out
+!
+!-- in/out
 
     REAL(wp), INTENT(IN) ::  rhopart                 !< particle density (kg/m3)
@@ -5636,13 +5339,9 @@
 
     REAL(wp) ::  vs
-
-    !
-    !-- acceleration of gravity:
+!
+!-- acceleration of gravity:
     REAL(wp), PARAMETER         ::  grav = 9.80665   !< m/s2
 
-
-
-    !
-    !-- sedimentation velocity
+!-- sedimentation velocity
     vs = rhopart * ( partsize**2.0_wp ) * grav * slipcor / ( 18.0_wp * visc )
 
@@ -5656,8 +5355,8 @@
       luc, ftop_lu, ustar_lu )
 
-
     IMPLICIT NONE  
 
-    ! -- in/out 
+!
+!-- in/out 
 
     INTEGER(iwp), INTENT(IN) ::  nwet        !< 1=rain, 9=snowcover
@@ -5675,8 +5374,6 @@
     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             !< acceleration of gravity (m/s2)
@@ -5693,10 +5390,8 @@
     REAL(wp), PARAMETER ::A_lu(nlu_dep) = &   
          (/3.0,  3.0,   2.0,  2.0,  7.0,  -99.,   10.0, 3.0, -99., -99.,  3.0, 7.0, -99., 2.0, -99./)
-    !
-    !--   grass  arabl crops conif decid  water  urba  othr  desr  ice   sav  trf   wai  med   sem     
-
-    !
-    !-- local
-
+!
+!--   grass  arabl crops conif decid  water  urba  othr  desr  ice   sav  trf   wai  med   sem     
+!
+!-- local
     REAL(wp) ::  kinvisc
     REAL(wp) ::  diff_part
@@ -5707,22 +5402,18 @@
     REAL(wp) ::  Einterc
     REAL(wp) ::  Reffic
-
-
-    !
-    !-- kinetic viscosity & diffusivity
+!
+!-- kinetic viscosity & diffusivity
     kinvisc = viscos1 / dens1    !< only needed at surface
 
     diff_part = kb * tsurf * slipcor / ( 3 * pi * viscos1 * partsize )
-
-    !
-    !-- Schmidt number
+!
+!-- Schmidt number
     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)
+!
+!-- 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_lu**2 / ( grav * kinvisc )
@@ -5738,25 +5429,21 @@
        Reffic = exp( -Stokes**0.5_wp )
     END IF
-    !
-    !-- when surface is wet all particles do not rebound:
+!
+!-- when surface is wet all particles do not rebound:
     IF ( nwet==1 )  Reffic = 1.0_wp
-    !
-    !-- determine impaction efficiency:
+!
+!-- determine impaction efficiency:
     Eimpac = ( stokes / ( alfa_lu(luc) + stokes ) )**beta
-
-    !
-    !-- sedimentation resistance:
+!
+!-- sedimentation resistance:
     Rs = 1.0_wp / ( epsilon0 * ustar_lu * ( 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)
-    !
-    !
+!-- 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)
 
     vd = 1.0_wp / ( ftop_lu + Rs + ftop_lu * Rs * vs1) + vs1
@@ -5772,8 +5459,9 @@
  SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffc, Rb )    
 
-   
+
     IMPLICIT NONE
 
-    !-- in/out 
+!
+!-- in/out 
 
     LOGICAL , INTENT(IN) ::  is_water
@@ -5784,24 +5472,21 @@
 
     REAL(wp), INTENT(OUT) ::  Rb                  !< boundary layer resistance
-
-    !
-    !-- const
+!
+!-- const
 
     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) 
-    !-- @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/diffc*kappa_stab*MAX(0.01_wp,ustar))
-    !      Rb = 1.0_wp / (kappa_stab*MAX(0.1_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.1_wp,ustar))
-    !ELSE
+!
+!-- 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/diffc*kappa_stab*MAX(0.01_wp,ustar))
+!--        Rb = 1.0_wp / (kappa_stab*MAX(0.1_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.1_wp,ustar))
+!--  ELSE
     Rb = 5.0_wp / MAX( 0.01_wp, ustar ) * ( thk / diffc )**0.67_wp
-    !END IF
+!--  END IF
 
  END SUBROUTINE get_rb_cell
@@ -5835,6 +5520,6 @@
     IMPLICIT NONE
 
-    !
-    !-- in/out 
+!
+!-- in/out 
 
     REAL(wp), INTENT(IN) ::  q                      !< specific humidity [(kg water)/(kg air)]
@@ -5842,16 +5527,13 @@
 
     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
-
 
 
@@ -5880,21 +5562,18 @@
     IMPLICIT NONE
 
-    !
-    !-- in/out
+!
+!-- 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:
+!
+!-- saturation vapor pressure:
     e_sat_w = p0 * exp( 17.67_wp * ( t - 273.16_wp ) / ( t - 29.66_wp ) )    !< [Pa]
 
  END FUNCTION saturationvaporpressure
-
 
 
@@ -5912,6 +5591,6 @@
     IMPLICIT NONE
 
-    !
-    !-- in/out 
+!
+!-- in/out 
 
     REAL(wp), INTENT(IN) ::  q    !< specific humidity [(kg water)/(kg air)]
@@ -5920,11 +5599,9 @@
 
     REAL(wp) ::  rh               !< relative humidity [%]
-
-    !
-    !-- relative humidity:
+!
+!-- relative humidity:
     rh = watervaporpartialpressure( q, p ) / saturationvaporpressure( t ) * 100.0_wp
 
  END FUNCTION relativehumidity_from_specifichumidity
-
 
      
