Index: palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_mod.f90	(revision 2688)
+++ palm/trunk/SOURCE/surface_mod.f90	(revision 2696)
@@ -1,5 +1,5 @@
 !> @file surface_mod.f90
 !------------------------------------------------------------------------------!
-! This file is part of PALM.
+! This file is part of the PALM model system.
 !
 ! PALM is free software: you can redistribute it and/or modify it under the
@@ -26,4 +26,16 @@
 ! -----------------
 ! $Id$
+! - Missing code for css added to surf_*, handling of surface_csflux updated (FK)
+! - Bugfixes in reading/writing restart data in case several surface types are 
+!   present at the same time (MS)
+! - Implementation of chemistry module (FK)
+! - Allocation of pt1 and qv1 now done for all surface types (MS)
+! - Revised classification of surface types
+! - Introduce offset values to simplify index determination of surface elements
+! - Dimensions of function get_topo_top_index (MS)
+! - added variables for USM 
+! - adapted to changes in USM (RvT)
+! 
+! 2688 2017-12-12 17:27:04Z Giersch
 ! Allocation and initialization of the latent heat flux (qsws) at the top of 
 ! the ocean domain in case of coupled runs. In addtion, a double if-query has 
@@ -84,4 +96,8 @@
 !> In addition, a further derived data structure is defined, in order to set
 !> boundary conditions at surfaces.  
+!> @todo For the moment, downward-facing surfaces are only classified as  
+!>        default type 
+!> @todo Clean up urban-surface variables (some of them are not used any more) 
+!> @todo Revise chemistry surface flux part (reduce loops?!) 
 !------------------------------------------------------------------------------!
  MODULE surface_mod
@@ -91,4 +107,8 @@
                momentumflux_input_conversion
 
+#if defined( __chem )
+    USE chem_modules
+#endif
+
     USE control_parameters               
 
@@ -126,5 +146,9 @@
     TYPE surf_type
 
+       INTEGER(iwp) :: ioff                                !< offset value in x-direction, used to determine index of surface element
+       INTEGER(iwp) :: joff                                !< offset value in y-direction, used to determine index of surface element
+       INTEGER(iwp) :: koff                                !< offset value in z-direction, used to determine index of surface element
        INTEGER(iwp) :: ns                                  !< number of surface elements on the PE
+
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< x-index linking to the PALM 3D-grid  
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< y-index linking to the PALM 3D-grid    
@@ -154,6 +178,7 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  z0q       !< roughness length for humidity
 
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt1       !< Specific humidity at first grid level (required for cloud_physics = .T. or cloud_droplets = .T.)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qv1       !< Potential temperature at first grid level (required for cloud_physics = .T. or cloud_droplets = .T.)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt1       !< Potential temperature at first grid level
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qv1       !< Specific humidity at first grid level
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  css     !< scaling parameter chemical species
 !
 !--    Define arrays for surface fluxes
@@ -169,4 +194,6 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  nrsws     !< surface flux nr
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  sasws     !< surface flux salinity
+       
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  cssws   !< surface flux chemical species
 !
 !--    Required for horizontal walls in production_e
@@ -179,25 +206,52 @@
 !
 !--    Variables required for LSM as well as for USM
-       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  pt_surface        !< skin-surface temperature
-       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net_l         !< net radiation 
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h          !< heat conductivity of soil (W/m/K) 
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_def      !< default heat conductivity of soil (W/m/K)   
-       
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  nzt_pavement  !< top index for pavement in soil
+       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  albedo_type   !< albedo type, for each fraction (wall,green,window or vegetation,pavement water)
+
        LOGICAL, DIMENSION(:), ALLOCATABLE  ::  building_surface    !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet)
+       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  building_covered    !< flag indicating that buildings are on top of orography, only used for vertical surfaces in LSM
        LOGICAL, DIMENSION(:), ALLOCATABLE  ::  pavement_surface    !< flag parameter for pavements
        LOGICAL, DIMENSION(:), ALLOCATABLE  ::  water_surface       !< flag parameter for water surfaces
-       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  vegetation_surface     !< flag parameter for natural land surfaces
+       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  vegetation_surface  !< flag parameter for natural land surfaces
+
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  albedo            !< broadband albedo for each surface fraction (LSM: vegetation, water, pavement; USM: wall, green, window)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  emissivity        !< emissivity of the surface, for each fraction  (LSM: vegetation, water, pavement; USM: wall, green, window)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  frac              !< relative surface fraction (LSM: vegetation, water, pavement; USM: wall, green, window)
+
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  aldif             !< albedo for longwave diffusive radiation, solar angle of 60Â°
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  aldir             !< albedo for longwave direct radiation, solar angle of 60Â°
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  asdif             !< albedo for shortwave diffusive radiation, solar angle of 60Â°
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  asdir             !< albedo for shortwave direct radiation, solar angle of 60Â°
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rrtm_aldif        !< albedo for longwave diffusive radiation, solar angle of 60Â°
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rrtm_aldir        !< albedo for longwave direct radiation, solar angle of 60Â°
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rrtm_asdif        !< albedo for shortwave diffusive radiation, solar angle of 60Â°
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rrtm_asdir        !< albedo for shortwave direct radiation, solar angle of 60Â°
+
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  pt_surface        !< skin-surface temperature
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net           !< net radiation 
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net_l         !< net radiation, used in USM
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h          !< heat conductivity of soil/ wall (W/m/K) 
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_green    !< heat conductivity of green soil (W/m/K) 
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_window   !< heat conductivity of windows (W/m/K) 
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_def      !< default heat conductivity of soil (W/m/K)   
+
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_lw_in           !< incoming longwave radiation
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_lw_out          !< emitted longwave radiation
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_sw_in           !< incoming shortwave radiation
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_sw_out          !< emitted shortwave radiation
+       
+
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_liq               !< liquid water coverage (of vegetated area)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_veg               !< vegetation coverage
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  f_sw_in             !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ghf              !< ground heat flux
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ghf                 !< ground heat flux
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  g_d                 !< coefficient for dependence of r_canopy on water vapour pressure deficit
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  lai                 !< leaf area index
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surface_u    !< coupling between surface and soil (depends on vegetation type) (W/m2/K)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surface_s    !< coupling between surface and soil (depends on vegetation type) (W/m2/K)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_liq         !< surface flux of latent heat (liquid water portion)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_soil        !< surface flux of latent heat (soil portion)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_veg         !< surface flux of latent heat (vegetation portion)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_liq            !< surface flux of latent heat (liquid water portion)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_soil           !< surface flux of latent heat (soil portion)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_veg            !< surface flux of latent heat (vegetation portion)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a                 !< aerodynamic resistance 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_canopy            !< canopy resistance
@@ -224,12 +278,20 @@
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  surface_types   !< array of types of wall parameters
 
-       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  isroof_surf         !< flag indication roof surfaces
-
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  albedo_surf         !< albedo of the surface
+       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  isroof_surf          !< flag indicating roof surfaces
+       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  ground_level         !< flag indicating ground floor level surfaces 
+
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  target_temp_summer  !< indoor target temperature summer
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  target_temp_winter  !< indoor target temperature summer        
+
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_surface           !< heat capacity of the wall surface skin (J/m2/K)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  emiss_surf          !< emissivity of the wall surface
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_surface_green     !< heat capacity of the green surface skin (J/m2/K)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_surface_window    !< heat capacity of the window surface skin (J/m2/K)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surf         !< heat conductivity between air and surface (W/m2/K)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  roughness_wall      !< roughness relative to concrete
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surf_green   !< heat conductivity between air and green surface (W/m2/K)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surf_window  !< heat conductivity between air and window surface (W/m2/K)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  thickness_wall      !< thickness of the wall, roof and soil layers
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  thickness_green     !< thickness of the green wall, roof and soil layers
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  thickness_window    !< thickness of the window wall, roof and soil layers
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  transmissivity      !< transmissivity of windows
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutsl           !< reflected shortwave radiation for local surface in i-th reflection
@@ -238,14 +300,16 @@
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  tt_surface_m        !< surface temperature tendency (K)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  tt_surface_window_m !< window surface temperature tendency (K)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  tt_surface_green_m  !< green surface temperature tendency (K)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  wshf                !< kinematic wall heat flux of sensible heat (actually no longer needed)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  wshf_eb             !< wall heat flux of sensible heat in wall normal direction
 
-
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb             !< wall ground heat flux
-
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_in_sw           !< incoming shortwave radiation
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_out_sw          !< emitted shortwave radiation
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_in_lw           !< incoming longwave radiation
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_out_lw          !< emitted longwave radiation
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_window      !< window ground heat flux
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_green       !< green ground heat flux
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  iwghf_eb            !< indoor wall ground heat flux
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  iwghf_eb_window     !< indoor window ground heat flux
+
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_lw_out_change_0
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinsw            !< shortwave radiation falling to local surface including radiation from reflections
@@ -262,4 +326,18 @@
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_wall_m         !< t_wall prognostic array 
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw                !< wall layer depths (m)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_window      !< volumetric heat capacity of the window material ( J m-3 K-1 ) (= 2.19E6)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_window         !< window grid spacing (center-center)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_window        !< 1/dz_window
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_window_stag    !< window grid spacing (edge-edge)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_window_stag   !< 1/dz_window_stag
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_window_m       !< t_window prognostic array 
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw_window         !< window layer depths (m)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_green       !< volumetric heat capacity of the green material ( J m-3 K-1 ) (= 2.19E6)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_green          !< green grid spacing (center-center)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_green         !< 1/dz_green
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_green_stag     !< green grid spacing (edge-edge)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_green_stag    !< 1/dz_green_stag
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_green_m        !< t_green prognostic array 
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw_green          !< green layer depths (m)
 
 
@@ -279,8 +357,18 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface  
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_av       !< average of wghf_eb
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_window_av  !< average of wghf_eb window
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_green_av   !< average of wghf_eb window
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  iwghf_eb_av        !< indoor average of wghf_eb
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  iwghf_eb_window_av !< indoor average of wghf_eb window
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  wshf_eb_av       !< average of wshf_eb
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_surf_av        !< average of wall surface temperature (K)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_surf_window_av !< average of window surface temperature (K)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_surf_green_av  !< average of green wall surface temperature (K)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_surf_whole_av  !< average of whole wall surface temperature (K)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_surf_10cm_av   !< average of the near surface temperature (K)
 
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  t_wall_av      !< Average of t_wall
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  t_window_av    !< Average of t_window
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  t_green_av     !< Average of t_green
 
     END TYPE surf_type
@@ -305,4 +393,5 @@
     INTERFACE get_topography_top_index
        MODULE PROCEDURE get_topography_top_index
+       MODULE PROCEDURE get_topography_top_index_ji
     END  INTERFACE get_topography_top_index
 
@@ -331,4 +420,9 @@
     END INTERFACE surface_last_actions
 
+    INTERFACE surface_restore_elements
+       MODULE PROCEDURE surface_restore_elements_1d
+       MODULE PROCEDURE surface_restore_elements_2d
+    END INTERFACE surface_restore_elements
+
 !
 !-- Public variables
@@ -338,6 +432,7 @@
 !-- Public subroutines and functions
     PUBLIC get_topography_top_index, init_bc, init_surfaces,                   &
-           init_surface_arrays, surface_read_restart_data,                     &
-           surface_write_restart_data, surface_last_actions
+           init_surface_arrays, surface_read_restart_data,                     & 
+           surface_restore_elements, surface_write_restart_data,               &
+           surface_last_actions
 
 
@@ -473,4 +568,6 @@
        INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v !< number of vertically-aligned urban surfaces 
 
+       LOGICAL ::  building                       !< flag indicating building grid point
+       LOGICAL ::  terrain                        !< flag indicating natural terrain grid point
 
        num_def_h = 0
@@ -481,4 +578,10 @@
        num_usm_v = 0
 !
+!--    Surfaces are classified according to the input data read from static
+!--    input file. If no input file is present, all surfaces are classified
+!--    either as natural, urban, or default, depending on the setting of 
+!--    land_surface and urban_surface. To control this, use the control
+!--    flag topo_no_distinct
+!
 !--    Count number of horizontal surfaces on local domain 
        DO  i = nxl, nxr
@@ -491,17 +594,34 @@
 !--                Check if grid point adjoins to any upward-facing horizontal
 !--                surface, e.g. the Earth surface, plane roofs, or ceilings.
+
                    IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )  THEN
 !
+!--                   Determine flags indicating terrain or building. 
+                      terrain  = BTEST( wall_flags_0(k-1,j,i), 5 )  .OR.       &
+                                 topo_no_distinct
+                      building = BTEST( wall_flags_0(k-1,j,i), 6 )  .OR.       &
+                                 topo_no_distinct
+!
 !--                   Land-surface type
-                      IF ( land_surface )  THEN
+                      IF ( land_surface  .AND.  terrain )  THEN
                          num_lsm_h    = num_lsm_h    + 1 
 !
 !--                   Urban surface tpye
-                      ELSEIF ( urban_surface )  THEN 
+                      ELSEIF ( urban_surface  .AND.  building )  THEN
                          num_usm_h    = num_usm_h    + 1 
 !
 !--                   Default-surface type
-                      ELSE
-                         num_def_h(0) = num_def_h(0) + 1 
+                      ELSEIF ( .NOT. land_surface    .AND.                     &
+                               .NOT. urban_surface )  THEN
+                               
+                         num_def_h(0) = num_def_h(0) + 1
+!
+!--                   Unclassifified surface-grid point. Give error message.
+                      ELSE 
+                         WRITE( message_string, * )                           &
+                                          'Unclassified upward-facing ' //    &
+                                          'surface element at '//             &
+                                          'grid point (k,j,i) = ', k, j, i
+                         CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 )
                       ENDIF
 
@@ -531,10 +651,29 @@
 !--                Northward-facing
                    IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) )  THEN
-                      IF ( urban_surface )  THEN
+!
+!--                   Determine flags indicating terrain or building
+
+                      terrain  = BTEST( wall_flags_0(k,j-1,i), 5 )  .OR.       &
+                                 topo_no_distinct
+                      building = BTEST( wall_flags_0(k,j-1,i), 6 )   .OR.       &
+                                 topo_no_distinct
+                      IF (  land_surface  .AND.  terrain )  THEN
+                         num_lsm_v(0) = num_lsm_v(0) + 1 
+                      ELSEIF ( urban_surface  .AND.  building )  THEN
                          num_usm_v(0) = num_usm_v(0) + 1 
-                      ELSEIF ( land_surface )  THEN
-                         num_lsm_v(0) = num_lsm_v(0) + 1 
-                      ELSE
+!
+!--                   Default-surface type
+                      ELSEIF ( .NOT. land_surface    .AND.                     &
+                               .NOT. urban_surface )  THEN
                          num_def_v(0) = num_def_v(0) + 1 
+!
+!--                   Unclassifified surface-grid point. Give error message.
+                      ELSE 
+                         WRITE( message_string, * )                           &
+                                          'Unclassified northward-facing ' // &
+                                          'surface element at '//             &
+                                          'grid point (k,j,i) = ', k, j, i
+                         CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 ) 
+
                       ENDIF
                    ENDIF
@@ -542,10 +681,28 @@
 !--                Southward-facing
                    IF ( .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) )  THEN
-                      IF ( urban_surface )  THEN
+!
+!--                   Determine flags indicating terrain or building
+                      terrain  = BTEST( wall_flags_0(k,j+1,i), 5 )  .OR.       &
+                                 topo_no_distinct
+                      building = BTEST( wall_flags_0(k,j+1,i), 6 )  .OR.       &
+                                 topo_no_distinct
+                      IF (  land_surface  .AND.  terrain )  THEN
+                         num_lsm_v(1) = num_lsm_v(1) + 1 
+                      ELSEIF ( urban_surface  .AND.  building )  THEN
                          num_usm_v(1) = num_usm_v(1) + 1 
-                      ELSEIF ( land_surface )  THEN
-                         num_lsm_v(1) = num_lsm_v(1) + 1 
-                      ELSE
+!
+!--                   Default-surface type
+                      ELSEIF ( .NOT. land_surface    .AND.                     &
+                               .NOT. urban_surface )  THEN
                          num_def_v(1) = num_def_v(1) + 1 
+!
+!--                   Unclassifified surface-grid point. Give error message.
+                      ELSE 
+                         WRITE( message_string, * )                           &
+                                          'Unclassified southward-facing ' // &
+                                          'surface element at '//             &
+                                          'grid point (k,j,i) = ', k, j, i
+                         CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 ) 
+
                       ENDIF
                    ENDIF
@@ -553,10 +710,28 @@
 !--                Eastward-facing
                    IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) )  THEN
-                      IF ( urban_surface )  THEN
+!
+!--                   Determine flags indicating terrain or building
+                      terrain  = BTEST( wall_flags_0(k,j,i-1), 5 )  .OR.       &
+                                 topo_no_distinct
+                      building = BTEST( wall_flags_0(k,j,i-1), 6 )  .OR.       &
+                                 topo_no_distinct
+                      IF (  land_surface  .AND.  terrain )  THEN
+                         num_lsm_v(2) = num_lsm_v(2) + 1 
+                      ELSEIF ( urban_surface  .AND.  building )  THEN
                          num_usm_v(2) = num_usm_v(2) + 1 
-                      ELSEIF ( land_surface )  THEN
-                         num_lsm_v(2) = num_lsm_v(2) + 1 
-                      ELSE
+!
+!--                   Default-surface type
+                      ELSEIF ( .NOT. land_surface    .AND.                     &
+                               .NOT. urban_surface )  THEN
                          num_def_v(2) = num_def_v(2) + 1 
+!
+!--                   Unclassifified surface-grid point. Give error message.
+                      ELSE 
+                         WRITE( message_string, * )                           &
+                                          'Unclassified eastward-facing ' //  &
+                                          'surface element at '//             &
+                                          'grid point (k,j,i) = ', k, j, i
+                         CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 ) 
+
                       ENDIF
                    ENDIF
@@ -564,10 +739,28 @@
 !--                Westward-facing
                    IF ( .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) )  THEN
-                      IF ( urban_surface )  THEN
-                         num_usm_v(3) = num_usm_v(3) + 1
-                      ELSEIF ( land_surface )  THEN
-                         num_lsm_v(3) = num_lsm_v(3) + 1  
-                      ELSE
+!
+!--                   Determine flags indicating terrain or building
+                      terrain  = BTEST( wall_flags_0(k,j,i+1), 5 )  .OR.       &
+                                 topo_no_distinct
+                      building = BTEST( wall_flags_0(k,j,i+1), 6 )  .OR.       &
+                                 topo_no_distinct
+                      IF (  land_surface  .AND.  terrain )  THEN
+                         num_lsm_v(3) = num_lsm_v(3) + 1 
+                      ELSEIF ( urban_surface  .AND.  building )  THEN
+                         num_usm_v(3) = num_usm_v(3) + 1 
+!
+!--                   Default-surface type
+                      ELSEIF ( .NOT. land_surface    .AND.                     &
+                               .NOT. urban_surface )  THEN
                          num_def_v(3) = num_def_v(3) + 1 
+!
+!--                   Unclassifified surface-grid point. Give error message.
+                      ELSE 
+                         WRITE( message_string, * )                           &
+                                          'Unclassified westward-facing ' //  &
+                                          'surface element at '//             &
+                                          'grid point (k,j,i) = ', k, j, i
+                         CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 ) 
+
                       ENDIF
                    ENDIF
@@ -651,5 +844,5 @@
 !--    Default type.
        DO  l = 0, 3
-          CALL allocate_surface_attributes_v ( surf_def_v(l), .FALSE.,         &
+          CALL allocate_surface_attributes_v ( surf_def_v(l),                  &
                                                nys, nyn, nxl, nxr )
        ENDDO
@@ -657,5 +850,5 @@
 !--    Natural type
        DO  l = 0, 3
-          CALL allocate_surface_attributes_v ( surf_lsm_v(l), .TRUE.,          &
+          CALL allocate_surface_attributes_v ( surf_lsm_v(l),                  &
                                                nys, nyn, nxl, nxr )
        ENDDO
@@ -663,5 +856,5 @@
 !--    Urban type
        DO  l = 0, 3
-          CALL allocate_surface_attributes_v ( surf_usm_v(l), .FALSE.,         &
+          CALL allocate_surface_attributes_v ( surf_usm_v(l),                  &
                                                nys, nyn, nxl, nxr )
        ENDDO
@@ -716,5 +909,4 @@
 !--    Surface-parallel wind velocity
        ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
-!      ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) )
 !
 !--    Roughness
@@ -746,4 +938,7 @@
        ALLOCATE ( surfaces%shf(1:surfaces%ns) )    
 !
+!--    surface temperature
+       ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) ) 
+!
 !--    Characteristic humidity and surface flux of latent heat
        IF ( humidity )  THEN
@@ -758,10 +953,16 @@
        ENDIF 
 !
-!--    When cloud physics is used, arrays for storing potential temperature and
-!--    specific humidity at first grid level are required. 
-       IF ( cloud_physics  .OR.  cloud_droplets )  THEN
-          ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
-          ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
-       ENDIF
+!--    Scaling parameter (cs*) and surface flux of chemical species
+#if defined ( __chem )
+       IF ( air_chemistry )  THEN
+          ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns)   )   
+          ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) ) 
+       ENDIF 
+#endif 
+!
+!--    Arrays for storing potential temperature and
+!--    specific humidity at first grid level
+       ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
+       ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
 !
 !--       
@@ -836,4 +1037,11 @@
        ENDIF 
 !
+!--    Chemical species flux
+#if defined( __chem )
+       IF ( air_chemistry )  THEN
+          ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) ) 
+       ENDIF 
+#endif 
+!
 !--       
        IF ( cloud_physics .AND. microphysics_morrison)  THEN
@@ -858,5 +1066,5 @@
 !> Allocating memory for vertical surface types. 
 !------------------------------------------------------------------------------!
-    SUBROUTINE allocate_surface_attributes_v( surfaces, lsm,                   &
+    SUBROUTINE allocate_surface_attributes_v( surfaces,                        &
                                               nys_l, nyn_l, nxl_l, nxr_l )
 
@@ -867,6 +1075,4 @@
        INTEGER(iwp) ::  nxl_l  !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
        INTEGER(iwp) ::  nxr_l  !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
-
-       LOGICAL         ::  lsm      !< flag indicating data type of natural land surface
 
        TYPE(surf_type) ::  surfaces !< respective surface type
@@ -909,10 +1115,9 @@
        ALLOCATE ( surfaces%us(1:surfaces%ns) )
 !
-!--    Allocate Obukhov length and bulk Richardson number. Only required 
+!--    Allocate Obukhov length and bulk Richardson number. Actually, at
+!--    vertical surfaces these are only required for natural surfaces.  
 !--    for natural land surfaces
-       IF ( lsm )  THEN
-          ALLOCATE( surfaces%ol(1:surfaces%ns)  ) 
-          ALLOCATE( surfaces%rib(1:surfaces%ns) ) 
-       ENDIF
+       ALLOCATE( surfaces%ol(1:surfaces%ns)  ) 
+       ALLOCATE( surfaces%rib(1:surfaces%ns) ) 
 !
 !--    Allocate arrays for surface momentum fluxes for u and v. For u at north- 
@@ -930,5 +1135,8 @@
 !--    Characteristic temperature and surface flux of sensible heat
        ALLOCATE ( surfaces%ts(1:surfaces%ns)  )    
-       ALLOCATE ( surfaces%shf(1:surfaces%ns) )    
+       ALLOCATE ( surfaces%shf(1:surfaces%ns) )
+!
+!--    surface temperature
+       ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) ) 
 !
 !--    Characteristic humidity and surface flux of latent heat
@@ -942,5 +1150,18 @@
           ALLOCATE ( surfaces%ss(1:surfaces%ns)   )   
           ALLOCATE ( surfaces%ssws(1:surfaces%ns) ) 
+       ENDIF
+!
+!--    Scaling parameter (cs*) and surface flux of chemical species
+#if defined( __chem )
+       IF ( air_chemistry )  THEN
+             ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns)   )   
+             ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) ) 
        ENDIF 
+#endif
+!
+!--    Arrays for storing potential temperature and
+!--    specific humidity at first grid level
+       ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
+       ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
 
        IF ( cloud_physics .AND. microphysics_seifert)  THEN
@@ -966,5 +1187,10 @@
 ! Description:
 ! ------------
-!> Initialize surface elements. 
+!> Initialize surface elements, i.e. set initial values for surface fluxes, 
+!> friction velocity, calcuation of start/end indices, etc. . 
+!> Please note, further initialization concerning 
+!> special surface characteristics, e.g. soil- and vegatation type, 
+!> building type, etc., is done in the land-surface and urban-surface module, 
+!> respectively.  
 !------------------------------------------------------------------------------!
     SUBROUTINE init_surfaces
@@ -1001,4 +1227,63 @@
        INTEGER(iwp), DIMENSION(0:3) ::  start_index_usm_v !< dummy to determing local start index in surface type for given (j,i), for vertical urban surfaces
 
+       LOGICAL ::  building     !< flag indicating building grid point
+       LOGICAL ::  terrain      !< flag indicating natural terrain grid point
+
+!
+!--    Set offset indices, i.e. index difference between surface element and 
+!--    surface-bounded grid point.
+!--    Upward facing - no horizontal offsets
+       surf_def_h(0:2)%ioff = 0
+       surf_def_h(0:2)%joff = 0
+
+       surf_lsm_h%ioff = 0
+       surf_lsm_h%joff = 0
+
+       surf_usm_h%ioff = 0
+       surf_usm_h%joff = 0
+!
+!--    Upward facing vertical offsets
+       surf_def_h(0)%koff   = -1
+       surf_lsm_h%koff      = -1
+       surf_usm_h%koff      = -1
+!
+!--    Downward facing vertical offset
+       surf_def_h(1:2)%koff = 1
+!
+!--    Vertical surfaces - no vertical offset
+       surf_def_v(0:3)%koff = 0
+       surf_lsm_v(0:3)%koff = 0
+       surf_usm_v(0:3)%koff = 0
+!
+!--    North- and southward facing - no offset in x
+       surf_def_v(0:1)%ioff = 0
+       surf_lsm_v(0:1)%ioff = 0
+       surf_usm_v(0:1)%ioff = 0
+!
+!--    Northward facing offset in y
+       surf_def_v(0)%joff = -1
+       surf_lsm_v(0)%joff = -1
+       surf_usm_v(0)%joff = -1
+!
+!--    Southward facing offset in y
+       surf_def_v(1)%joff = 1
+       surf_lsm_v(1)%joff = 1
+       surf_usm_v(1)%joff = 1
+
+!
+!--    East- and westward facing - no offset in y
+       surf_def_v(2:3)%joff = 0
+       surf_lsm_v(2:3)%joff = 0
+       surf_usm_v(2:3)%joff = 0
+!
+!--    Eastward facing offset in x
+       surf_def_v(2)%ioff = -1
+       surf_lsm_v(2)%ioff = -1
+       surf_usm_v(2)%ioff = -1
+!
+!--    Westward facing offset in y
+       surf_def_v(3)%ioff = 1
+       surf_lsm_v(3)%ioff = 1
+       surf_usm_v(3)%ioff = 1
 
 !
@@ -1039,10 +1324,15 @@
 !--                Upward-facing surface. Distinguish between differet surface types.
 !--                To do, think about method to flag natural and non-natural 
-!--                surfaces. Only to ask for land_surface or urban surface 
-!--                is just a work-around. 
+!--                surfaces. 
                    IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )  THEN 
 !
+!--                   Determine flags indicating terrain or building
+                      terrain  = BTEST( wall_flags_0(k-1,j,i), 5 )  .OR.       &
+                                 topo_no_distinct
+                      building = BTEST( wall_flags_0(k-1,j,i), 6 )  .OR.       &
+                                 topo_no_distinct
+!
 !--                   Natural surface type          
-                      IF ( land_surface )  THEN
+                      IF ( land_surface  .AND.  terrain )  THEN
                          CALL initialize_horizontal_surfaces( k, j, i,         &
                                                               surf_lsm_h,      &
@@ -1052,5 +1342,5 @@
 !
 !--                   Urban surface tpye
-                      ELSEIF ( urban_surface )  THEN
+                      ELSEIF ( urban_surface  .AND.  building )  THEN
                          CALL initialize_horizontal_surfaces( k, j, i,         &
                                                               surf_usm_h,      &
@@ -1069,5 +1359,7 @@
                    ENDIF  
 !
-!--                downward-facing surface, first, model top
+!--                downward-facing surface, first, model top. Please note, 
+!--                for the moment, downward-facing surfaces are always of 
+!--                default type
                    IF ( k == nzt  .AND.  use_top_fluxes )  THEN
                       CALL initialize_top( k, j, i, surf_def_h(2),             &
@@ -1087,5 +1379,11 @@
 !                  Start with northward-facing surface.
                    IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) )  THEN
-                      IF ( urban_surface )  THEN
+!
+!--                   Determine flags indicating terrain or building
+                      terrain  = BTEST( wall_flags_0(k,j-1,i), 5 )  .OR.       &
+                                 topo_no_distinct
+                      building = BTEST( wall_flags_0(k,j-1,i), 6 )  .OR.       &
+                                 topo_no_distinct
+                      IF ( urban_surface  .AND.  building )  THEN
                          CALL initialize_vertical_surfaces( 0, k, j, i,        &
                                                             surf_usm_v(0),     &
@@ -1094,5 +1392,5 @@
                                                             .FALSE., .FALSE.,  &             
                                                             .FALSE., .TRUE. ) 
-                      ELSEIF ( land_surface )  THEN
+                      ELSEIF ( land_surface  .AND.  terrain )  THEN
                          CALL initialize_vertical_surfaces( 0, k, j, i,        &
                                                             surf_lsm_v(0),     &
@@ -1113,5 +1411,11 @@
 !--                southward-facing surface
                    IF ( .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) )  THEN
-                      IF ( urban_surface )  THEN
+!
+!--                   Determine flags indicating terrain or building
+                      terrain  = BTEST( wall_flags_0(k,j+1,i), 5 )  .OR.       &
+                                 topo_no_distinct
+                      building = BTEST( wall_flags_0(k,j+1,i), 6 )  .OR.       &
+                                 topo_no_distinct
+                      IF ( urban_surface  .AND.  building )  THEN
                          CALL initialize_vertical_surfaces( 1, k, j, i,        &
                                                             surf_usm_v(1),     &
@@ -1120,5 +1424,5 @@
                                                             .FALSE., .FALSE.,  &
                                                             .TRUE., .FALSE. )
-                      ELSEIF ( land_surface )  THEN
+                      ELSEIF ( land_surface  .AND.  terrain )  THEN
                          CALL initialize_vertical_surfaces( 1, k, j, i,        &
                                                             surf_lsm_v(1),     &
@@ -1139,5 +1443,11 @@
 !--                eastward-facing surface
                    IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) )  THEN
-                      IF ( urban_surface )  THEN
+!
+!--                   Determine flags indicating terrain or building
+                      terrain  = BTEST( wall_flags_0(k,j,i-1), 5 )  .OR.       &
+                                 topo_no_distinct
+                      building = BTEST( wall_flags_0(k,j,i-1), 6 )  .OR.       &
+                                 topo_no_distinct
+                      IF ( urban_surface  .AND.  building )  THEN
                          CALL initialize_vertical_surfaces( 2, k, j, i,        &
                                                             surf_usm_v(2),     &
@@ -1146,5 +1456,5 @@
                                                             .TRUE., .FALSE.,   &
                                                             .FALSE., .FALSE. ) 
-                      ELSEIF ( land_surface )  THEN
+                      ELSEIF ( land_surface  .AND.  terrain )  THEN
                          CALL initialize_vertical_surfaces( 2, k, j, i,        &
                                                             surf_lsm_v(2),     &
@@ -1165,5 +1475,11 @@
 !--                westward-facing surface
                    IF ( .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) )  THEN
-                      IF ( urban_surface )  THEN
+!
+!--                   Determine flags indicating terrain or building
+                      terrain  = BTEST( wall_flags_0(k,j,i+1), 5 )  .OR.       &
+                                 topo_no_distinct
+                      building = BTEST( wall_flags_0(k,j,i+1), 6 )  .OR.       &
+                                 topo_no_distinct
+                      IF ( urban_surface  .AND.  building )  THEN
                          CALL initialize_vertical_surfaces( 3, k, j, i,        &
                                                             surf_usm_v(3),     &
@@ -1172,5 +1488,5 @@
                                                            .FALSE., .TRUE.,    &
                                                            .FALSE., .FALSE. ) 
-                      ELSEIF ( land_surface )  THEN
+                      ELSEIF ( land_surface  .AND.  terrain )  THEN
                          CALL initialize_vertical_surfaces( 3, k, j, i,        &
                                                             surf_lsm_v(3),     &
@@ -1285,4 +1601,5 @@
        ENDDO
 
+
        CONTAINS
 
@@ -1305,4 +1622,6 @@
              INTEGER(iwp)  ::  num_h            !< current number of surface element
              INTEGER(iwp)  ::  num_h_kji        !< dummy increment
+             INTEGER(iwp)  ::  lsp              !< running index chemical species
+             INTEGER(iwp)  ::  lsp_pr           !< running index chemical species??
 
              LOGICAL       ::  upward_facing    !< flag indicating upward-facing surface
@@ -1356,6 +1675,6 @@
 !
 !--          Initialization in case of constant profiles
-             ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )&
-             THEN
+             ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0&
+                .OR.  INDEX(initializing_actions, 'inifor' ) /= 0 )  THEN
 
                 surf%ol(num_h)   = surf%z_mo(num_h) / zeta_min
@@ -1403,4 +1722,13 @@
 
              IF ( passive_scalar )  surf%ss(num_h) = 0.0_wp
+
+#if defined( __chem )
+             DO  lsp = 1, nvar
+                IF ( air_chemistry )  surf%css(lsp,num_h) = 0.0_wp
+             ENDDO
+#endif
+!
+!--          Set initial value for surface temperature
+             surf%pt_surface(num_h) = pt_surface
 !
 !--          Inititalize surface fluxes of sensible and latent heat, as well as
@@ -1414,5 +1742,5 @@
 !--                   if random_heatflux is set. This case, shf is initialized later.
                       IF ( .NOT. random_heatflux )  THEN
-                         surf%shf(num_h) = surface_heatflux *               &
+                         surf%shf(num_h) = surface_heatflux *                  &
                                                  heatflux_input_conversion(nzb)
 !
@@ -1420,5 +1748,5 @@
 !--                      prescribed wall heatflux
                          IF ( k-1 /= 0 )  THEN
-                            surf%shf(num_h) = wall_heatflux(0) *            &
+                            surf%shf(num_h) = wall_heatflux(0) *               &
                                                  heatflux_input_conversion(k-1)
                          ENDIF
@@ -1476,4 +1804,31 @@
                 ENDIF
 
+#if defined( __chem )
+                IF ( air_chemistry )  THEN
+                   lsp_pr = 1
+                   DO  WHILE ( TRIM( surface_csflux_name( lsp_pr ) ) /= 'novalue' )   !<'novalue' is the default
+                      DO  lsp = 1, nvar
+!
+!--                      Assign surface flux for each variable species
+                         IF ( TRIM( spc_names(lsp) ) == TRIM( surface_csflux_name(lsp_pr) ) )  THEN   
+                            IF ( upward_facing )  THEN
+                               IF ( constant_csflux(lsp_pr) )  THEN
+                                  surf%cssws(lsp,num_h) = surface_csflux(lsp_pr)
+
+                                  IF ( k-1 /= 0 )                                    &
+                                     surf%cssws(lsp,num_h) = wall_csflux(lsp,0)
+                               ELSE
+                                  surf%cssws(lsp,num_h) = 0.0_wp
+                               ENDIF
+                            ELSE
+                               surf%cssws(lsp,num_h) = wall_csflux(lsp,5)
+                            ENDIF
+                         ENDIF
+                      ENDDO
+                      lsp_pr = lsp_pr + 1
+                   ENDDO
+                ENDIF
+#endif
+
                 IF ( ocean )  THEN
                    IF ( upward_facing )  THEN 
@@ -1508,4 +1863,5 @@
              INTEGER(iwp)  ::  num_h            !< current number of surface element
              INTEGER(iwp)  ::  num_h_kji        !< dummy increment
+             INTEGER(iwp)  ::  lsp              !< running index for chemical species
 
              TYPE( surf_type ) :: surf          !< respective surface type
@@ -1543,4 +1899,13 @@
                 surf%ssws(num_h) = top_scalarflux
 !
+!--          Prescribe top chemical species' flux
+#if defined( __chem )
+             DO  lsp = 1, nvar
+                IF ( air_chemistry  .AND.  constant_top_csflux(lsp) )  THEN 
+                   surf%cssws(lsp,num_h) = top_csflux(lsp)
+                ENDIF
+             ENDDO
+#endif
+!
 !--          Prescribe top salinity flux
              IF ( ocean .AND. constant_top_salinityflux)                       &
@@ -1575,5 +1940,5 @@
              IMPLICIT NONE 
 
-             INTEGER(iwp)  ::  component !<
+             INTEGER(iwp)  ::  component       !< index of wall_fluxes_ array for respective orientation 
              INTEGER(iwp)  ::  i               !< running index x-direction
              INTEGER(iwp)  ::  j               !< running index x-direction
@@ -1582,4 +1947,5 @@
              INTEGER(iwp)  ::  num_v           !< current number of surface element
              INTEGER(iwp)  ::  num_v_kji       !< current number of surface element at (j,i)
+             INTEGER(iwp)  ::  lsp             !< running index for chemical species
 
 
@@ -1609,20 +1975,20 @@
 !--          etc., on surface type (further below)
              IF ( north_facing )  THEN
-                surf%facing(num_v) = IBSET( surf%facing(num_v), 0 )  
+                surf%facing(num_v) = 5 !IBSET( surf%facing(num_v), 0 )  
                 component          = 4
              ENDIF
 
              IF ( south_facing )  THEN
-                surf%facing(num_v) = IBSET( surf%facing(num_v), 1 ) 
+                surf%facing(num_v) = 6 !IBSET( surf%facing(num_v), 1 ) 
                 component          = 3
              ENDIF
 
              IF ( east_facing )  THEN
-                surf%facing(num_v) = IBSET( surf%facing(num_v), 2 )
+                surf%facing(num_v) = 7 !IBSET( surf%facing(num_v), 2 )
                 component          = 2
              ENDIF
 
              IF ( west_facing )  THEN
-                surf%facing(num_v) = IBSET( surf%facing(num_v), 3 ) 
+                surf%facing(num_v) = 8 !IBSET( surf%facing(num_v), 3 ) 
                 component          = 1
              ENDIF
@@ -1647,4 +2013,7 @@
              surf%ts(num_v)    = 0.0_wp
              surf%shf(num_v)   = wall_heatflux(component)
+!
+!--          Set initial value for surface temperature
+             surf%pt_surface(num_v) = pt_surface
 
              IF ( humidity )  THEN
@@ -1673,4 +2042,14 @@
                 surf%ssws(num_v) = wall_scalarflux(component)
              ENDIF
+
+#if defined( __chem )
+             IF ( air_chemistry )  THEN        
+                DO  lsp = 1, nvar
+                   surf%css(lsp,num_v)   = 0.0_wp
+                   surf%cssws(lsp,num_v) = wall_csflux(lsp,component)
+                ENDDO
+             ENDIF
+#endif
+
 !
 !--          So far, salinityflux at vertical surfaces is simply zero 
@@ -1692,15 +2071,13 @@
 !> Determines topography-top index at given (j,i)-position.  
 !------------------------------------------------------------------------------!
-    FUNCTION get_topography_top_index( j, i, grid )
-
-       USE kinds
+    FUNCTION get_topography_top_index_ji( j, i, grid )
 
        IMPLICIT NONE 
 
-       CHARACTER(LEN=*) ::  grid                      !< flag to distinquish between staggered grids
-       INTEGER(iwp)     ::  i                         !< grid index in x-dimension
-       INTEGER(iwp)     ::  ibit                      !< bit position where topography information is stored on respective grid
-       INTEGER(iwp)     ::  j                         !< grid index in y-dimension
-       INTEGER(iwp)     ::  get_topography_top_index  !< topography top index
+       CHARACTER(LEN=*) ::  grid                         !< flag to distinquish between staggered grids
+       INTEGER(iwp)     ::  i                            !< grid index in x-dimension
+       INTEGER(iwp)     ::  ibit                         !< bit position where topography information is stored on respective grid
+       INTEGER(iwp)     ::  j                            !< grid index in y-dimension
+       INTEGER(iwp)     ::  get_topography_top_index_ji  !< topography top index
 
        SELECT CASE ( TRIM( grid ) )
@@ -1729,9 +2106,57 @@
        END SELECT
 
-       get_topography_top_index = MAXLOC(                                      &
+       get_topography_top_index_ji = MAXLOC(                                   &
                                      MERGE( 1, 0,                              &
                                             BTEST( wall_flags_0(:,j,i), ibit ) &
                                           ), DIM = 1                           &
-                                        ) - 1
+                                           ) - 1
+
+       RETURN
+
+    END FUNCTION get_topography_top_index_ji
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Determines topography-top index at each (j,i)-position.  
+!------------------------------------------------------------------------------!
+    FUNCTION get_topography_top_index( grid )
+
+       IMPLICIT NONE 
+
+       CHARACTER(LEN=*) ::  grid                      !< flag to distinquish between staggered grids
+       INTEGER(iwp)     ::  ibit                      !< bit position where topography information is stored on respective grid
+       INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) ::  get_topography_top_index  !< topography top index
+
+       SELECT CASE ( TRIM( grid ) )
+
+          CASE ( 's'     )
+             ibit = 12
+          CASE ( 'u'     )
+             ibit = 14
+          CASE ( 'v'     )
+             ibit = 16
+          CASE ( 'w'     )
+             ibit = 18
+          CASE ( 's_out' )
+             ibit = 24
+          CASE ( 'u_out' )
+             ibit = 26
+          CASE ( 'v_out' )
+             ibit = 27
+          CASE ( 'w_out' )
+             ibit = 28
+          CASE DEFAULT
+!
+!--          Set default to scalar grid
+             ibit = 12 
+
+       END SELECT
+
+       get_topography_top_index(nys:nyn,nxl:nxr) = MAXLOC(                     &
+                         MERGE( 1, 0,                                          &
+                                 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), ibit )&
+                              ), DIM = 1                                       &
+                                                         ) - 1
 
        RETURN
@@ -1754,4 +2179,5 @@
        INTEGER(iwp)                 ::  j    !< running index y-direction
        INTEGER(iwp)                 ::  l    !< index surface type orientation
+       INTEGER(iwp)                 ::  lsp  !< running index chemical species
        INTEGER(iwp)                 ::  m    !< running index for surface elements on individual surface array
        INTEGER(iwp), DIMENSION(0:3) ::  mm   !< running index for surface elements on gathered surface array
@@ -1781,5 +2207,5 @@
        DO  l = 0, 3
           surf_v(l)%ns = ns_v_on_file(l)
-          CALL allocate_surface_attributes_v( surf_v(l), .FALSE.,              &
+          CALL allocate_surface_attributes_v( surf_v(l),                       &
                                               nys, nyn, nxl, nxr )
        ENDDO
@@ -1823,4 +2249,16 @@
                    IF ( ALLOCATED( surf_def_h(l)%ssws ) )                      &
                       surf_h(l)%ssws(mm(l))    = surf_def_h(l)%ssws(m)
+#if defined( __chem )
+                   IF ( ALLOCATED( surf_def_h(l)%css ) )  THEN 
+                      DO  lsp = 1,nvar
+                         surf_h(l)%css(lsp,mm(l)) = surf_def_h(l)%css(lsp,m)
+                      ENDDO
+                   ENDIF
+                   IF ( ALLOCATED( surf_def_h(l)%cssws ) )  THEN 
+                      DO  lsp = 1,nvar
+                         surf_h(l)%cssws(lsp,mm(l)) = surf_def_h(l)%cssws(lsp,m)
+                      ENDDO
+                   ENDIF
+#endif
                    IF ( ALLOCATED( surf_def_h(l)%ncsws ) )                     &
                       surf_h(l)%ncsws(mm(l))   = surf_def_h(l)%ncsws(m)
@@ -1866,4 +2304,16 @@
                       IF ( ALLOCATED( surf_lsm_h%ssws ) )                      &
                          surf_h(0)%ssws(mm(0))    = surf_lsm_h%ssws(m)
+#if defined( __chem )
+                      IF ( ALLOCATED( surf_lsm_h%css ) )  THEN                  
+                         DO  lsp = 1, nvar
+                            surf_h(0)%css(lsp,mm(0)) = surf_lsm_h%css(lsp,m)
+                         ENDDO
+                      ENDIF
+                      IF ( ALLOCATED( surf_lsm_h%cssws ) )  THEN
+                         DO  lsp = 1, nvar
+                            surf_h(0)%cssws(lsp,mm(0)) = surf_lsm_h%cssws(lsp,m)
+                         ENDDO 
+                      ENDIF
+#endif
                       IF ( ALLOCATED( surf_lsm_h%ncsws ) )                     &
                          surf_h(0)%ncsws(mm(0))   = surf_lsm_h%ncsws(m)
@@ -1909,4 +2359,16 @@
                       IF ( ALLOCATED( surf_usm_h%ssws ) )                      &
                          surf_h(0)%ssws(mm(0))    = surf_usm_h%ssws(m)
+#if defined( __chem )
+                      IF ( ALLOCATED( surf_usm_h%css ) )  THEN             
+                         DO lsp = 1, nvar
+                            surf_h(0)%css(lsp,mm(0)) = surf_usm_h%css(lsp,m)
+                         ENDDO
+                      ENDIF
+                      IF ( ALLOCATED( surf_usm_h%cssws ) )  THEN             
+                         DO lsp = 1, nvar
+                            surf_h(0)%cssws(lsp,mm(0)) = surf_usm_h%cssws(lsp,m)
+                         ENDDO
+                      ENDIF
+#endif
                       IF ( ALLOCATED( surf_usm_h%ncsws ) )                     &
                          surf_h(0)%ncsws(mm(0))   = surf_usm_h%ncsws(m)
@@ -1926,15 +2388,27 @@
 
           ENDDO
-          IF ( l == 0 )  THEN
-             surf_h(l)%start_index = MAX( surf_def_h(l)%start_index,           &
-                                          surf_lsm_h%start_index,              &
-                                          surf_usm_h%start_index )
-             surf_h(l)%end_index   = MAX( surf_def_h(l)%end_index,             &
-                                          surf_lsm_h%end_index,                &
-                                          surf_usm_h%end_index )
-          ELSE
-             surf_h(l)%start_index = surf_def_h(l)%start_index
-             surf_h(l)%end_index   = surf_def_h(l)%end_index
-          ENDIF
+!
+!--       Gather start- and end indices                                        
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                IF ( surf_def_h(l)%start_index(j,i) <=                         &
+                     surf_def_h(l)%end_index(j,i) )  THEN
+                   surf_h(l)%start_index(j,i) = surf_def_h(l)%start_index(j,i)
+                   surf_h(l)%end_index(j,i)   = surf_def_h(l)%end_index(j,i)
+                ENDIF
+                IF ( l == 0 )  THEN
+                   IF ( surf_lsm_h%start_index(j,i) <=                         &
+                        surf_lsm_h%end_index(j,i) )  THEN
+                      surf_h(l)%start_index(j,i) = surf_lsm_h%start_index(j,i)
+                      surf_h(l)%end_index(j,i)   = surf_lsm_h%end_index(j,i)
+                   ENDIF
+                   IF ( surf_usm_h%start_index(j,i) <=                         &
+                        surf_usm_h%end_index(j,i) )  THEN
+                      surf_h(l)%start_index(j,i) = surf_usm_h%start_index(j,i)
+                      surf_h(l)%end_index(j,i)   = surf_usm_h%end_index(j,i)
+                   ENDIF
+                ENDIF
+             ENDDO
+          ENDDO
        ENDDO
 
@@ -1972,4 +2446,16 @@
                    IF ( ALLOCATED( surf_def_v(l)%ssws ) )                      &
                       surf_v(l)%ssws(mm(l))    = surf_def_v(l)%ssws(m)
+#if defined( __chem )
+                   IF ( ALLOCATED( surf_def_v(l)%css ) )  THEN               
+                      DO  lsp = 1, nvar
+                         surf_v(l)%css(lsp,mm(l)) = surf_def_v(l)%css(lsp,m)
+                      ENDDO
+                   ENDIF
+                   IF ( ALLOCATED( surf_def_v(l)%cssws ) )  THEN               
+                      DO  lsp = 1, nvar
+                         surf_v(l)%cssws(lsp,mm(l)) = surf_def_v(l)%cssws(lsp,m)
+                      ENDDO
+                   ENDIF
+#endif
                    IF ( ALLOCATED( surf_def_v(l)%ncsws ) )                     &
                       surf_v(l)%ncsws(mm(l))   = surf_def_v(l)%ncsws(m)
@@ -2020,4 +2506,16 @@
                    IF ( ALLOCATED( surf_lsm_v(l)%ssws ) )                      &
                       surf_v(l)%ssws(mm(l))    = surf_lsm_v(l)%ssws(m)
+#if defined( __chem )
+                   IF ( ALLOCATED( surf_lsm_v(l)%css ) )  THEN              
+                      DO  lsp = 1, nvar
+                         surf_v(l)%css(lsp,mm(l)) = surf_lsm_v(l)%css(lsp,m)
+                      ENDDO
+                   ENDIF
+                   IF ( ALLOCATED( surf_lsm_v(l)%cssws ) )  THEN              
+                      DO  lsp = 1, nvar
+                         surf_v(l)%cssws(lsp,mm(l)) = surf_lsm_v(l)%cssws(lsp,m)
+                      ENDDO
+                   ENDIF
+#endif
                    IF ( ALLOCATED( surf_lsm_v(l)%ncsws ) )                     &
                       surf_v(l)%ncsws(mm(l))   = surf_lsm_v(l)%ncsws(m)
@@ -2068,4 +2566,16 @@
                    IF ( ALLOCATED( surf_usm_v(l)%ssws ) )                      &
                       surf_v(l)%ssws(mm(l))    = surf_usm_v(l)%ssws(m)
+#if defined( __chem )
+                   IF ( ALLOCATED( surf_usm_v(l)%css ) )  THEN              
+                      DO  lsp = 1, nvar
+                         surf_v(l)%css(lsp,mm(l)) = surf_usm_v(l)%css(lsp,m)
+                      ENDDO
+                   ENDIF
+                   IF ( ALLOCATED( surf_usm_v(l)%cssws ) )  THEN              
+                      DO  lsp = 1, nvar
+                         surf_v(l)%cssws(lsp,mm(l)) = surf_usm_v(l)%cssws(lsp,m)
+                      ENDDO
+                   ENDIF
+#endif
                    IF ( ALLOCATED( surf_usm_v(l)%ncsws ) )                     &
                       surf_v(l)%ncsws(mm(l))   = surf_usm_v(l)%ncsws(m)
@@ -2171,4 +2681,14 @@
              WRITE ( 14 )  surf_h(l)%ssws
           ENDIF
+#if defined ( __chem )
+          WRITE ( 14 )  'surf_h(' // dum // ')%css                 '
+          IF ( ALLOCATED ( surf_h(l)%css ) )  THEN 
+             WRITE ( 14 )  surf_h(l)%css
+          ENDIF
+          WRITE ( 14 )  'surf_h(' // dum // ')%cssws               '
+          IF ( ALLOCATED ( surf_h(l)%cssws ) )  THEN 
+             WRITE ( 14 )  surf_h(l)%cssws
+          ENDIF
+#endif
           WRITE ( 14 )  'surf_h(' // dum // ')%qcsws               '  
           IF ( ALLOCATED ( surf_h(l)%qcsws ) )  THEN
@@ -2254,4 +2774,14 @@
              WRITE ( 14 )  surf_v(l)%ssws
           ENDIF
+#if defined( __chem )
+          WRITE ( 14 )  'surf_v(' // dum // ')%css                 '  
+          IF ( ALLOCATED ( surf_v(l)%css ) )  THEN
+             WRITE ( 14 )  surf_v(l)%css
+          ENDIF
+          WRITE ( 14 )  'surf_v(' // dum // ')%cssws               '  
+          IF ( ALLOCATED ( surf_v(l)%cssws ) )  THEN
+             WRITE ( 14 )  surf_v(l)%cssws
+          ENDIF
+#endif
           WRITE ( 14 )  'surf_v(' // dum // ')%qcsws               '  
           IF ( ALLOCATED ( surf_v(l)%qcsws ) )  THEN
@@ -2399,5 +2929,5 @@
        DO  l = 0, 3
           surf_v(l)%ns = ns_v_on_file(l)
-          CALL allocate_surface_attributes_v( surf_v(l), .FALSE.,              &
+          CALL allocate_surface_attributes_v( surf_v(l),                       &
                                               nys_on_file, nyn_on_file,        &
                                               nxl_on_file, nxr_on_file )
@@ -2486,4 +3016,12 @@
                       IF ( ALLOCATED( surf_h(0)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_h(0)%ssws
+#if defined( __chem )
+                   CASE ( 'surf_h(0)%css' )
+                      IF ( ALLOCATED( surf_h(0)%css )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_h(0)%css
+                   CASE ( 'surf_h(0)%cssws' )         
+                      IF ( ALLOCATED( surf_h(0)%cssws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_h(0)%cssws
+#endif
                    CASE ( 'surf_h(0)%qcsws' )         
                       IF ( ALLOCATED( surf_h(0)%qcsws )  .AND.  kk == 1 )      &
@@ -2554,4 +3092,12 @@
                       IF ( ALLOCATED( surf_h(1)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_h(1)%ssws
+#if defined( __chem )
+                   CASE ( 'surf_h(1)%css' )
+                      IF ( ALLOCATED( surf_h(1)%css )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_h(1)%css
+                   CASE ( 'surf_h(1)%cssws' )         
+                      IF ( ALLOCATED( surf_h(1)%cssws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_h(1)%cssws
+#endif
                    CASE ( 'surf_h(1)%qcsws' )         
                       IF ( ALLOCATED( surf_h(1)%qcsws )  .AND.  kk == 1 )      &
@@ -2622,4 +3168,12 @@
                       IF ( ALLOCATED( surf_h(2)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_h(2)%ssws
+#if defined( __chem )
+                   CASE ( 'surf_h(2)%css' )
+                      IF ( ALLOCATED( surf_h(2)%css )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_h(2)%css
+                   CASE ( 'surf_h(2)%cssws' )         
+                      IF ( ALLOCATED( surf_h(2)%cssws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_h(2)%cssws
+#endif
                    CASE ( 'surf_h(2)%qcsws' )         
                       IF ( ALLOCATED( surf_h(2)%qcsws )  .AND.  kk == 1 )      &
@@ -2686,4 +3240,12 @@
                       IF ( ALLOCATED( surf_v(0)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_v(0)%ssws
+#if defined( __chem )
+                   CASE ( 'surf_v(0)%css' ) 
+                      IF ( ALLOCATED( surf_v(0)%css )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_v(0)%css
+                   CASE ( 'surf_v(0)%cssws' )         
+                      IF ( ALLOCATED( surf_v(0)%cssws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_v(0)%cssws
+#endif
                    CASE ( 'surf_v(0)%qcsws' )         
                       IF ( ALLOCATED( surf_v(0)%qcsws )  .AND.  kk == 1 )      &
@@ -2757,4 +3319,12 @@
                       IF ( ALLOCATED( surf_v(1)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_v(1)%ssws
+#if defined( __chem )
+                   CASE ( 'surf_v(1)%css' ) 
+                      IF ( ALLOCATED( surf_v(1)%css )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_v(1)%css
+                   CASE ( 'surf_v(1)%cssws' )         
+                      IF ( ALLOCATED( surf_v(1)%cssws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_v(1)%cssws
+#endif
                    CASE ( 'surf_v(1)%qcsws' )         
                       IF ( ALLOCATED( surf_v(1)%qcsws )  .AND.  kk == 1 )      &
@@ -2828,4 +3398,12 @@
                       IF ( ALLOCATED( surf_v(2)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_v(2)%ssws
+#if defined( __chem )
+                   CASE ( 'surf_v(2)%css' ) 
+                      IF ( ALLOCATED( surf_v(2)%css )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_v(2)%css
+                   CASE ( 'surf_v(2)%cssws' )         
+                      IF ( ALLOCATED( surf_v(2)%cssws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_v(2)%cssws
+#endif
                    CASE ( 'surf_v(2)%qcsws' )         
                       IF ( ALLOCATED( surf_v(2)%qcsws )  .AND.  kk == 1 )      &
@@ -2899,4 +3477,12 @@
                       IF ( ALLOCATED( surf_v(3)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_v(3)%ssws
+#if defined( __chem )
+                   CASE ( 'surf_v(3)%css' ) 
+                      IF ( ALLOCATED( surf_v(3)%css )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_v(3)%css
+                   CASE ( 'surf_v(3)%cssws' )         
+                      IF ( ALLOCATED( surf_v(3)%cssws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_v(3)%cssws
+#endif
                    CASE ( 'surf_v(3)%qcsws' )         
                       IF ( ALLOCATED( surf_v(3)%qcsws )  .AND.  kk == 1 )      &
@@ -3047,9 +3633,10 @@
              INTEGER(iwp)      ::  m_file      !< respective surface-element index of current surface array 
              INTEGER(iwp)      ::  m_target    !< respecitve surface-element index of surface array on file
+             INTEGER(iwp)      ::  lsp         !< running index chemical species
 
              TYPE( surf_type ) ::  surf_target !< target surface type
              TYPE( surf_type ) ::  surf_file   !< surface type on file
 
-             IF ( SCAN( TRIM( field_chr ), '%us' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%us' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%us )  .AND.                        &
                      ALLOCATED( surf_file%us   ) )                             & 
@@ -3057,5 +3644,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%ol' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%ol' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%ol )  .AND.                        &
                      ALLOCATED( surf_file%ol   ) )                             & 
@@ -3063,5 +3650,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%usws' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%usws' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%usws )  .AND.                      &
                      ALLOCATED( surf_file%usws   ) )                           & 
@@ -3069,5 +3656,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%vsws' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%vsws' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%vsws )  .AND.                      &
                      ALLOCATED( surf_file%vsws   ) )                           & 
@@ -3075,5 +3662,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%ts' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%ts' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%ts )  .AND.                        &
                      ALLOCATED( surf_file%ts   ) )                             & 
@@ -3081,5 +3668,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%shf' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%shf' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%shf )  .AND.                       &
                      ALLOCATED( surf_file%shf   ) )                            & 
@@ -3087,5 +3674,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%qs' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%qs' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%qs )  .AND.                        &
                      ALLOCATED( surf_file%qs   ) )                             & 
@@ -3093,5 +3680,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%qsws' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%qsws' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%qsws )  .AND.                      &
                      ALLOCATED( surf_file%qsws   ) )                           & 
@@ -3099,5 +3686,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%ss' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%ss' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%ss )  .AND.                        &
                      ALLOCATED( surf_file%ss   ) )                             & 
@@ -3105,5 +3692,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%ssws' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%ssws' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%ssws )  .AND.                      &
                      ALLOCATED( surf_file%ssws   ) )                           & 
@@ -3111,5 +3698,24 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%qcs' ) /= 0 )  THEN 
+#if defined( __chem )
+             IF ( INDEX( TRIM( field_chr ), '%css' ) /= 0 )  THEN 
+                IF ( ALLOCATED( surf_target%css )  .AND.                     &
+                     ALLOCATED( surf_file%css   ) )  THEN                  
+                   DO  lsp = 1, nvar
+                      surf_target%css(lsp,m_target) = surf_file%css(lsp,m_file)
+                   ENDDO
+                ENDIF
+             ENDIF
+             IF ( INDEX( TRIM( field_chr ), '%cssws' ) /= 0 )  THEN 
+                IF ( ALLOCATED( surf_target%cssws )  .AND.                     &
+                     ALLOCATED( surf_file%cssws   ) )  THEN                  
+                   DO  lsp = 1, nvar
+                      surf_target%cssws(lsp,m_target) = surf_file%cssws(lsp,m_file)
+                   ENDDO
+                ENDIF
+             ENDIF
+#endif
+
+             IF ( INDEX( TRIM( field_chr ), '%qcs' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%qcs )  .AND.                       &
                      ALLOCATED( surf_file%qcs   ) )                            & 
@@ -3117,5 +3723,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%qcsws' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%qcsws' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%qcsws )  .AND.                     &
                      ALLOCATED( surf_file%qcsws   ) )                          & 
@@ -3123,5 +3729,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%ncs' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%ncs' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%ncs )  .AND.                       &
                      ALLOCATED( surf_file%ncs   ) )                            & 
@@ -3129,5 +3735,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%ncsws' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%ncsws' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%ncsws )  .AND.                     &
                      ALLOCATED( surf_file%ncsws   ) )                          & 
@@ -3135,5 +3741,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%qrs' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%qrs' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%qrs )  .AND.                       &
                      ALLOCATED( surf_file%qrs   ) )                            & 
@@ -3141,5 +3747,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%qrsws' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%qrsws' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%qrsws )  .AND.                     &
                      ALLOCATED( surf_file%qrsws   ) )                          & 
@@ -3147,5 +3753,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%nrs' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%nrs' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%nrs )  .AND.                       &
                      ALLOCATED( surf_file%nrs   ) )                            & 
@@ -3153,5 +3759,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%nrsws' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%nrsws' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%nrsws )  .AND.                     &
                      ALLOCATED( surf_file%nrsws   ) )                          & 
@@ -3159,5 +3765,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%sasws' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%sasws' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%sasws )  .AND.                     &
                      ALLOCATED( surf_file%sasws   ) )                          & 
@@ -3165,5 +3771,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%mom_uv' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%mom_uv' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%mom_flux_uv )  .AND.               &
                      ALLOCATED( surf_file%mom_flux_uv   ) )                    & 
@@ -3172,5 +3778,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%mom_w' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%mom_w' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%mom_flux_w )  .AND.                &
                      ALLOCATED( surf_file%mom_flux_w   ) )                     & 
@@ -3179,5 +3785,5 @@
              ENDIF
 
-             IF ( SCAN( TRIM( field_chr ), '%mom_tke' ) /= 0 )  THEN 
+             IF ( INDEX( TRIM( field_chr ), '%mom_tke' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%mom_flux_tke )  .AND.              &
                      ALLOCATED( surf_file%mom_flux_tke   ) )                   & 
@@ -3214,4 +3820,124 @@
     END SUBROUTINE surface_last_actions
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Routine maps surface data read from file after restart - 1D arrays
+!------------------------------------------------------------------------------!
+    SUBROUTINE surface_restore_elements_1d( surf_target, surf_file,            &
+                                            start_index_c,                     &
+                                            start_index_on_file,               &
+                                            end_index_on_file,                 &
+                                            nxlc, nysc, nxlf, nxrf, nysf, nynf,&
+                                            nys_on_file, nyn_on_file,          &
+                                            nxl_on_file,nxr_on_file )
+
+       IMPLICIT NONE
+    
+       INTEGER(iwp) ::  i         !< running index along x-direction, refers to former domain size
+       INTEGER(iwp) ::  ic        !< running index along x-direction, refers to current domain size
+       INTEGER(iwp) ::  j         !< running index along y-direction, refers to former domain size
+       INTEGER(iwp) ::  jc        !< running index along y-direction, refers to former domain size        
+       INTEGER(iwp) ::  m         !< surface-element index on file
+       INTEGER(iwp) ::  mm        !< surface-element index on current subdomain
+       INTEGER(iwp) ::  nxlc      !< index of left boundary on current subdomain
+       INTEGER(iwp) ::  nxlf      !< index of left boundary on former subdomain 
+       INTEGER(iwp) ::  nxrf      !< index of right boundary on former subdomain
+       INTEGER(iwp) ::  nysc      !< index of north boundary on current subdomain
+       INTEGER(iwp) ::  nynf      !< index of north boundary on former subdomain
+       INTEGER(iwp) ::  nysf      !< index of south boundary on former subdomain
+
+       INTEGER(iwp) ::  nxl_on_file !< leftmost index on file
+       INTEGER(iwp) ::  nxr_on_file !< rightmost index on file
+       INTEGER(iwp) ::  nyn_on_file !< northmost index on file
+       INTEGER(iwp) ::  nys_on_file !< southmost index on file
+
+       INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) ::  start_index_c             
+       INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: & 
+                            start_index_on_file   !< start index of surface elements on file
+       INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: & 
+                            end_index_on_file     !< end index of surface elements on file
+       
+       REAL(wp), DIMENSION(:) ::  surf_target !< target surface type
+       REAL(wp), DIMENSION(:) ::  surf_file   !< surface type on file
+       
+       ic = nxlc
+       DO  i = nxlf, nxrf
+          jc = nysc
+          DO  j = nysf, nynf
+
+             mm = start_index_c(jc,ic)
+             DO  m = start_index_on_file(j,i), end_index_on_file(j,i)
+                surf_target(mm) = surf_file(m)
+                mm = mm + 1
+             ENDDO
+
+             jc = jc + 1
+          ENDDO
+          ic = ic + 1
+       ENDDO
+
+
+    END SUBROUTINE surface_restore_elements_1d
+   
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Routine maps surface data read from file after restart - 2D arrays
+!------------------------------------------------------------------------------! 
+    SUBROUTINE surface_restore_elements_2d( surf_target, surf_file,            &
+                                            start_index_c,                     &
+                                            start_index_on_file,               &
+                                            end_index_on_file,                 &
+                                            nxlc, nysc, nxlf, nxrf, nysf, nynf,&
+                                            nys_on_file, nyn_on_file,          &
+                                            nxl_on_file,nxr_on_file )
+
+       IMPLICIT NONE
+    
+       INTEGER(iwp) ::  i         !< running index along x-direction, refers to former domain size
+       INTEGER(iwp) ::  ic        !< running index along x-direction, refers to current domain size
+       INTEGER(iwp) ::  j         !< running index along y-direction, refers to former domain size
+       INTEGER(iwp) ::  jc        !< running index along y-direction, refers to former domain size        
+       INTEGER(iwp) ::  m         !< surface-element index on file
+       INTEGER(iwp) ::  mm        !< surface-element index on current subdomain
+       INTEGER(iwp) ::  nxlc      !< index of left boundary on current subdomain
+       INTEGER(iwp) ::  nxlf      !< index of left boundary on former subdomain 
+       INTEGER(iwp) ::  nxrf      !< index of right boundary on former subdomain
+       INTEGER(iwp) ::  nysc      !< index of north boundary on current subdomain
+       INTEGER(iwp) ::  nynf      !< index of north boundary on former subdomain
+       INTEGER(iwp) ::  nysf      !< index of south boundary on former subdomain
+
+       INTEGER(iwp) ::  nxl_on_file !< leftmost index on file
+       INTEGER(iwp) ::  nxr_on_file !< rightmost index on file
+       INTEGER(iwp) ::  nyn_on_file !< northmost index on file
+       INTEGER(iwp) ::  nys_on_file !< southmost index on file
+
+       INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) ::  start_index_c !< start index of surface type
+       INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: & 
+                            start_index_on_file   !< start index of surface elements on file
+       INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: & 
+                            end_index_on_file     !< end index of surface elements on file
+       
+       REAL(wp), DIMENSION(:,:) ::  surf_target !< target surface type
+       REAL(wp), DIMENSION(:,:) ::  surf_file   !< surface type on file
+       
+       ic = nxlc
+       DO  i = nxlf, nxrf
+          jc = nysc
+          DO  j = nysf, nynf
+             mm = start_index_c(jc,ic)
+             DO  m = start_index_on_file(j,i), end_index_on_file(j,i)
+                surf_target(:,mm) = surf_file(:,m)
+                mm = mm + 1
+             ENDDO
+
+             jc = jc + 1
+          ENDDO
+          ic = ic + 1
+       ENDDO
+
+    END SUBROUTINE surface_restore_elements_2d
+
 
  END MODULE surface_mod
