Index: palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_mod.f90	(revision 4598)
+++ palm/trunk/SOURCE/surface_mod.f90	(revision 4671)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Implementation of downward facing USM and LSM surfaces
+!
+! 4598 2020-07-10 10:13:23Z suehring
 ! Revise surface-element mapping in mpi-io restart branch
 !
@@ -543,7 +546,7 @@
     TYPE (surf_type), DIMENSION(0:2), TARGET ::  surf_def_h  !< horizontal default surfaces (Up, Down, and Top)
     TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_def_v  !< vertical default surfaces (North, South, East, West)
-    TYPE (surf_type)                , TARGET ::  surf_lsm_h  !< horizontal natural land surfaces, so far only upward-facing
+    TYPE (surf_type), DIMENSION(0:1), TARGET ::  surf_lsm_h  !< horizontal natural land surfaces (Up, Down)
     TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_lsm_v  !< vertical land surfaces (North, South, East, West)
-    TYPE (surf_type)                , TARGET ::  surf_usm_h  !< horizontal urban surfaces, so far only upward-facing
+    TYPE (surf_type), DIMENSION(0:1), TARGET ::  surf_usm_h  !< horizontal urban surfaces (Up, Down)
     TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_usm_v  !< vertical urban surfaces (North, South, East, West)
 
@@ -855,8 +858,10 @@
     INTEGER(iwp) ::  k         !< running index z-direction
     INTEGER(iwp) ::  l         !< index variable for surface facing
-    INTEGER(iwp) ::  num_lsm_h !< number of horizontally-aligned natural surfaces
-    INTEGER(iwp) ::  num_usm_h !< number of horizontally-aligned urban surfaces
+    INTEGER(iwp) ::  kk        !< auxiliary index z-direction
+    INTEGER(iwp) ::  kd        !< direction index
 
     INTEGER(iwp), DIMENSION(0:2) ::  num_def_h !< number of horizontally-aligned default surfaces
+    INTEGER(iwp), DIMENSION(0:1) ::  num_lsm_h !< number of horizontally-aligned natural surfaces
+    INTEGER(iwp), DIMENSION(0:1) ::  num_usm_h !< number of horizontally-aligned urban surfaces
     INTEGER(iwp), DIMENSION(0:3) ::  num_def_v !< number of vertically-aligned default surfaces
     INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v !< number of vertically-aligned natural surfaces
@@ -891,49 +896,51 @@
              IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
 !
-!--             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_total_0(k-1,j,i), 0 ) )  THEN
-!
-!--                Determine flags indicating a terrain surface, a building surface,
-                   terrain  = BTEST( wall_flags_total_0(k-1,j,i), 5 )  .OR.  topo_no_distinct
-                   building = BTEST( wall_flags_total_0(k-1,j,i), 6 )  .OR.  topo_no_distinct
-!
-!--                Unresolved_building indicates a surface with equal height as terrain but with a
-!--                non-grid resolved building on top. These surfaces will be flagged as urban
-!--                surfaces.
-                   unresolved_building = BTEST( wall_flags_total_0(k-1,j,i), 5 )  .AND.            &
-                                         BTEST( wall_flags_total_0(k-1,j,i), 6 )
-!
-!--                Land-surface type
-                   IF ( land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
-                      num_lsm_h    = num_lsm_h    + 1
-!
-!--                Urban surface tpye
-                   ELSEIF ( urban_surface  .AND.  building )  THEN
-                      num_usm_h = num_usm_h + 1
-!
-!--                Default-surface type
-                   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.
+!--             Check if grid point adjoins to any upward- and downward-facing horizontal surface,
+!--             e.g. the Earth surface, plane roofs, or ceilings.
+                DO kk = k-1, k+1, 2
+!
+!--                Check for top-fluxes
+                   IF ( kk == nzt+1  .AND.  use_top_fluxes )  THEN
+                      num_def_h(2) = num_def_h(2) + 1
                    ELSE
-                      WRITE( message_string, * ) 'Unclassified upward-facing surface element '//   &
-                                                 'at grid point (k,j,i) = ', k, j, i
-                      CALL message( 'surface_mod', 'PA0698', 1, 2, myid, 6, 0 )
+!
+!--                   set direction index of the potential surface
+                      kd = MERGE( 0, 1, kk == k-1 )
+!
+!--                   test the adjacent grid cell
+                      IF ( .NOT. BTEST( wall_flags_total_0(kk,j,i), 0 ) )  THEN
+!
+!--                      Determine flags indicating a terrain surface, a building surface,
+                         terrain  = BTEST( wall_flags_total_0(kk,j,i), 5 )  .OR.  topo_no_distinct
+                         building = BTEST( wall_flags_total_0(kk,j,i), 6 )  .OR.  topo_no_distinct
+!
+!--                      Unresolved_building indicates a surface with equal height as terrain but with a
+!--                      non-grid resolved building on top. These surfaces will be flagged as urban
+!--                      surfaces.
+                         unresolved_building = BTEST( wall_flags_total_0(kk,j,i), 5 )  .AND.            &
+                                               BTEST( wall_flags_total_0(kk,j,i), 6 )
+!
+!--                      Land-surface type
+                         IF ( land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
+                            num_lsm_h(kd)    = num_lsm_h(kd)    + 1
+!
+!--                      Urban surface tpye
+                         ELSEIF ( urban_surface  .AND.  building )  THEN
+                            num_usm_h(kd) = num_usm_h(kd) + 1
+!
+!--                      Default-surface type
+                         ELSEIF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
+                            num_def_h(kd) = num_def_h(kd) + 1
+!
+!--                      Unclassifified surface-grid point. Give error message.
+                         ELSE
+                            WRITE( message_string, * ) 'Unclassified ',                            &
+                                  TRIM(MERGE('  upward','downward',kk==0)),                        &
+                                  '-facing surface element at grid point (k,j,i) = ', k, j, i
+                            CALL message( 'surface_mod', 'PA0698', 1, 2, myid, 6, 0 )
+                         ENDIF
+                      ENDIF
                    ENDIF
-
-                ENDIF
-!
-!--             Check for top-fluxes
-                IF ( k == nzt  .AND.  use_top_fluxes )  THEN
-                   num_def_h(2) = num_def_h(2) + 1
-!
-!--             Check for any other downward-facing surface. So far only for default surface type.
-                ELSEIF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i), 0 ) )  THEN
-                   num_def_h(1) = num_def_h(1) + 1
-                ENDIF
-
+                ENDDO
              ENDIF
           ENDDO
@@ -1075,9 +1082,15 @@
     surf_def_h(2)%ns = num_def_h(2)
 !
-!-- Horizontal surface, natural type, so far only upward-facing
-    surf_lsm_h%ns    = num_lsm_h
-!
-!-- Horizontal surface, urban type, so far only upward-facing
-    surf_usm_h%ns    = num_usm_h
+!-- Horizontal surface, natural type, upward facing
+    surf_lsm_h(0)%ns    = num_lsm_h(0)
+!
+!-- Horizontal surface, natural type, downward facing
+    surf_lsm_h(1)%ns    = num_lsm_h(1)
+!
+!-- Horizontal surface, urban type, upward facing
+    surf_usm_h(0)%ns    = num_usm_h(0)
+!
+!-- Horizontal surface, urban type, downward facing
+    surf_usm_h(1)%ns    = num_usm_h(1)
 !
 !-- Vertical surface, default type, northward facing
@@ -1127,8 +1140,12 @@
 !
 !-- Allocate required attributes for horizontal surfaces - natural type.
-    CALL allocate_surface_attributes_h ( surf_lsm_h, nys, nyn, nxl, nxr )
+    DO  l = 0, 1
+       CALL allocate_surface_attributes_h ( surf_lsm_h(l), nys, nyn, nxl, nxr )
+    ENDDO
 !
 !-- Allocate required attributes for horizontal surfaces - urban type.
-    CALL allocate_surface_attributes_h ( surf_usm_h, nys, nyn, nxl, nxr )
+    DO  l = 0, 1
+       CALL allocate_surface_attributes_h ( surf_usm_h(l), nys, nyn, nxl, nxr )
+    ENDDO
 
 !
@@ -1182,7 +1199,7 @@
     !$ACC COPYIN(surf_def_h(0:2)) &
     !$ACC COPYIN(surf_def_v(0:3)) &
-    !$ACC COPYIN(surf_lsm_h) &
+    !$ACC COPYIN(surf_lsm_h(0:1)) &
     !$ACC COPYIN(surf_lsm_v(0:3)) &
-    !$ACC COPYIN(surf_usm_h) &
+    !$ACC COPYIN(surf_usm_h(0:1)) &
     !$ACC COPYIN(surf_usm_v(0:3))
 !
@@ -1199,5 +1216,7 @@
 !
 !-- Copy data in surf_lsm_h
-    CALL enter_surface_attributes_h( surf_lsm_h )
+    DO  l = 0, 1
+       CALL enter_surface_attributes_h( surf_lsm_h(l) )
+    ENDDO
 !
 !-- Copy data in surf_lsm_v(0:3)
@@ -1207,5 +1226,7 @@
 !
 !-- Copy data in surf_usm_h
-    CALL enter_surface_attributes_h( surf_usm_h )
+    DO  l = 0, 1
+       CALL enter_surface_attributes_h( surf_usm_h(l) )
+    ENDDO
 !
 !-- Copy data in surf_usm_v(0:3)
@@ -1241,5 +1262,7 @@
 !
 !-- Delete data in surf_lsm_h
-    CALL exit_surface_attributes_h( surf_lsm_h )
+    DO  l = 0, 1
+       CALL exit_surface_attributes_h( surf_lsm_h(l) )
+    ENDDO
 !
 !-- Delete data in surf_lsm_v(0:3)
@@ -1249,5 +1272,7 @@
 !
 !-- Delete data in surf_usm_h
-    CALL exit_surface_attributes_h( surf_usm_h )
+    DO  l = 0, 1
+       CALL exit_surface_attributes_h( surf_usm_h(l) )
+    ENDDO
 !
 !-- Delete data in surf_usm_v(0:3)
@@ -1259,7 +1284,7 @@
     !$ACC DELETE(surf_def_h(0:2)) &
     !$ACC DELETE(surf_def_v(0:3)) &
-    !$ACC DELETE(surf_lsm_h) &
+    !$ACC DELETE(surf_lsm_h(0:1)) &
     !$ACC DELETE(surf_lsm_v(0:3)) &
-    !$ACC DELETE(surf_usm_h) &
+    !$ACC DELETE(surf_usm_h(0:1)) &
     !$ACC DELETE(surf_usm_v(0:3))
 
@@ -2178,16 +2203,18 @@
     INTEGER(iwp) ::  j  !< running index y-direction
     INTEGER(iwp) ::  k  !< running index z-direction
-
-    INTEGER(iwp)  ::  start_index_lsm_h  !< dummy to determing local start index in surface type for given (j,i), for horizontal
-                                         !< natural surfaces
-    INTEGER(iwp)  ::  start_index_usm_h  !< dummy to determing local start index in surface type for given (j,i), for horizontal
-                                         !< urban surfaces
-
-    INTEGER(iwp)  ::  num_lsm_h      !< current number of horizontal surface element, natural type
-    INTEGER(iwp)  ::  num_lsm_h_kji  !< dummy to determing local end index in surface type for given (j,i), for for horizonal
-                                     !< natural surfaces
-    INTEGER(iwp)  ::  num_usm_h      !< current number of horizontal surface element, urban type
-    INTEGER(iwp)  ::  num_usm_h_kji  !< dummy to determing local end index in surface type for given (j,i), for for horizonal urban
-                                     !< surfaces
+    INTEGER(iwp) ::  kk !< auxiliary index z-direction
+    INTEGER(iwp) ::  l  !< direction index
+
+
+    INTEGER(iwp), DIMENSION(0:1)  ::  start_index_lsm_h !< dummy to determing local start index in surface type for given (j,i), for horizontal
+                                                        !< natural surfaces
+    INTEGER(iwp), DIMENSION(0:1)  ::  start_index_usm_h !< dummy to determing local start index in surface type for given (j,i), for horizontal
+                                                        !< urban surfaces
+    INTEGER(iwp), DIMENSION(0:1)  ::  num_lsm_h         !< current number of horizontal surface element, natural type
+    INTEGER(iwp), DIMENSION(0:1)  ::  num_lsm_h_kji     !< dummy to determing local end index in surface type for given (j,i), for for horizonal
+                                                        !< natural surfaces
+    INTEGER(iwp), DIMENSION(0:1)  ::  num_usm_h         !< current number of horizontal surface element, urban type
+    INTEGER(iwp), DIMENSION(0:1)  ::  num_usm_h_kji     !< dummy to determing local end index in surface type for given (j,i), for for horizonal urban
+                                                        !< surfaces
 
     INTEGER(iwp), DIMENSION(0:2) ::  num_def_h          !< current number of horizontal surface element, default type
@@ -2197,13 +2224,13 @@
                                                         !< for horizontal default surfaces
 
-    INTEGER(iwp), DIMENSION(0:3) ::  num_def_v      !< current number of vertical surface element, default type
-    INTEGER(iwp), DIMENSION(0:3) ::  num_def_v_kji  !< dummy to determing local end index in surface type for given (j,i),
-                                                    !< for vertical default surfaces
-    INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v      !< current number of vertical surface element, natural type
-    INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v_kji  !< dummy to determing local end index in surface type for given (j,i),
-                                                    !< for vertical natural surfaces
-    INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v      !< current number of vertical surface element, urban type
-    INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v_kji  !< dummy to determing local end index in surface type for given (j,i),
-                                                    !< for vertical urban surfaces
+    INTEGER(iwp), DIMENSION(0:3) ::  num_def_v          !< current number of vertical surface element, default type
+    INTEGER(iwp), DIMENSION(0:3) ::  num_def_v_kji      !< dummy to determing local end index in surface type for given (j,i),
+                                                        !< for vertical default surfaces
+    INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v          !< current number of vertical surface element, natural type
+    INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v_kji      !< dummy to determing local end index in surface type for given (j,i),
+                                                        !< for vertical natural surfaces
+    INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v          !< current number of vertical surface element, urban type
+    INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v_kji      !< dummy to determing local end index in surface type for given (j,i),
+                                                        !< for vertical urban surfaces
 
     INTEGER(iwp), DIMENSION(0:3) ::  start_index_def_v  !< dummy to determing local start index in surface type for given (j,i),
@@ -2214,8 +2241,8 @@
                                                         !< for vertical urban surfaces
 
-    LOGICAL ::  building             !< flag indicating building grid point
-    LOGICAL ::  terrain              !< flag indicating natural terrain grid point
-    LOGICAL ::  unresolved_building  !< flag indicating a grid point where actually a building is defined but not resolved by the
-                                     !< vertical grid
+    LOGICAL ::  building                                !< flag indicating building grid point
+    LOGICAL ::  terrain                                 !< flag indicating natural terrain grid point
+    LOGICAL ::  unresolved_building                     !< flag indicating a grid point where actually a building is defined but not resolved by the
+                                                        !< vertical grid
 !
 !-- Set offset indices, i.e. index difference between surface element and surface-bounded grid point.
@@ -2224,17 +2251,19 @@
     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
+    surf_lsm_h(0:1)%ioff = 0
+    surf_lsm_h(0:1)%joff = 0
+
+    surf_usm_h(0:1)%ioff = 0
+    surf_usm_h(0:1)%joff = 0
 !
 !-- Upward facing vertical offsets
     surf_def_h(0)%koff = -1
-    surf_lsm_h%koff    = -1
-    surf_usm_h%koff    = -1
+    surf_lsm_h(0)%koff = -1
+    surf_usm_h(0)%koff = -1
 !
 !-- Downward facing vertical offset
     surf_def_h(1:2)%koff = 1
+    surf_lsm_h(1)%koff   = 1
+    surf_usm_h(1)%koff   = 1
 !
 !-- Vertical surfaces - no vertical offset
@@ -2279,8 +2308,8 @@
     num_def_v(0:3) = 1
 
-    num_lsm_h      = 1
+    num_lsm_h(0:1) = 1
     num_lsm_v(0:3) = 1
 
-    num_usm_h      = 1
+    num_usm_h(0:1) = 1
     num_usm_v(0:3) = 1
 
@@ -2288,8 +2317,8 @@
     start_index_def_v(0:3) = 1
 
-    start_index_lsm_h      = 1
+    start_index_lsm_h(0:1) = 1
     start_index_lsm_v(0:3) = 1
 
-    start_index_usm_h      = 1
+    start_index_usm_h(0:1) = 1
     start_index_usm_v(0:3) = 1
 
@@ -2308,46 +2337,48 @@
              IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
 !
-!--             Upward-facing surface. Distinguish between differet surface types.
-!--             To do, think about method to flag natural and non-natural surfaces.
-                IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i), 0 ) )  THEN
-!
-!--                Determine flags indicating terrain or building
-                   terrain  = BTEST( wall_flags_total_0(k-1,j,i), 5 )  .OR.  topo_no_distinct
-                   building = BTEST( wall_flags_total_0(k-1,j,i), 6 )  .OR.  topo_no_distinct
-
-!
-!--                Unresolved_building indicates a surface with equal height as terrain but with a
-!--                non-grid resolved building on top. These surfaces will be flagged as urban
-!--                surfaces.
-                   unresolved_building = BTEST( wall_flags_total_0(k-1,j,i), 5 ) .AND.             &
-                                         BTEST( wall_flags_total_0(k-1,j,i), 6 )
-!
-!--                Natural surface type
-                   IF ( land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
-                      CALL initialize_horizontal_surfaces( k, j, i, surf_lsm_h, num_lsm_h,         &
-                                                           num_lsm_h_kji, .TRUE., .FALSE. )
-!
-!--                Urban surface tpye
-                   ELSEIF ( urban_surface  .AND.  building )  THEN
-                      CALL initialize_horizontal_surfaces( k, j, i, surf_usm_h, num_usm_h,         &
-                                                           num_usm_h_kji, .TRUE., .FALSE. )
-!
-!--                Default surface type
+!--             Check if grid point adjoins to any upward- and downward-facing horizontal surface,
+!--             e.g. the Earth surface, plane roofs, or ceilings.
+                DO kk = k-1, k+1, 2
+!
+!--                Check for top-fluxes first
+                   IF ( kk == nzt+1  .AND.  use_top_fluxes )  THEN
+                      CALL initialize_top( k, j, i, surf_def_h(2), num_def_h(2), num_def_h_kji(2) )
                    ELSE
-                      CALL initialize_horizontal_surfaces( k, j, i, surf_def_h(0), num_def_h(0),   &
-                                                           num_def_h_kji(0), .TRUE., .FALSE. )
+!
+!--                   set direction index of the potential surface
+                      l = MERGE( 0, 1, kk == k-1 )
+!
+!--                   Upward- or donward facing surface. Distinguish between differet surface types.
+!--                   To do, think about method to flag natural and non-natural surfaces.
+                      IF ( .NOT. BTEST( wall_flags_total_0(kk,j,i), 0 ) )  THEN
+!
+!--                      Determine flags indicating terrain or building
+                         terrain  = BTEST( wall_flags_total_0(kk,j,i), 5 )  .OR.  topo_no_distinct
+                         building = BTEST( wall_flags_total_0(kk,j,i), 6 )  .OR.  topo_no_distinct
+!
+!--                      Unresolved_building indicates a surface with equal height as terrain but with a
+!--                      non-grid resolved building on top. These surfaces will be flagged as urban
+!--                      surfaces.
+                         unresolved_building = BTEST( wall_flags_total_0(kk,j,i), 5 ) .AND.       &
+                                               BTEST( wall_flags_total_0(kk,j,i), 6 )
+!
+!--                      Natural surface type
+                         IF ( land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
+                            CALL initialize_horizontal_surfaces( k, j, i, surf_lsm_h(l), num_lsm_h(l), &
+                                                                 num_lsm_h_kji(l), .TRUE., .FALSE. )
+!
+!--                      Urban surface tpye
+                         ELSEIF ( urban_surface  .AND.  building )  THEN
+                            CALL initialize_horizontal_surfaces( k, j, i, surf_usm_h(l), num_usm_h(l), &
+                                                                 num_usm_h_kji(l), .TRUE., .FALSE. )
+!
+!--                      Default surface type
+                         ELSE
+                            CALL initialize_horizontal_surfaces( k, j, i, surf_def_h(l), num_def_h(l), &
+                                                                 num_def_h_kji(l), .TRUE., .FALSE. )
+                         ENDIF
+                      ENDIF
                    ENDIF
-                ENDIF
-!
-!--             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), num_def_h(2), num_def_h_kji(2) )
-!
-!--             Check for any other downward-facing surface. So far only for default surface type.
-                ELSEIF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i), 0 ) )  THEN
-                   CALL initialize_horizontal_surfaces( k, j, i, surf_def_h(1), num_def_h(1),      &
-                                                        num_def_h_kji(1), .FALSE., .TRUE. )
-                ENDIF
+                ENDDO
 !
 !--             Check for vertical walls and, if required, initialize it.
@@ -2485,13 +2516,17 @@
           start_index_def_h(2)           = surf_def_h(2)%end_index(j,i) + 1
 !
-!--       Horizontal natural land surfaces
-          surf_lsm_h%start_index(j,i)    = start_index_lsm_h
-          surf_lsm_h%end_index(j,i)      = surf_lsm_h%start_index(j,i) + num_lsm_h_kji - 1
-          start_index_lsm_h              = surf_lsm_h%end_index(j,i) + 1
-!
-!--       Horizontal urban surfaces
-          surf_usm_h%start_index(j,i)    = start_index_usm_h
-          surf_usm_h%end_index(j,i)      = surf_usm_h%start_index(j,i) + num_usm_h_kji - 1
-          start_index_usm_h              = surf_usm_h%end_index(j,i) + 1
+!--       Upward- and downward-facing horizontal land and urban surfaces
+          DO l = 0, 1
+!
+!--          Horizontal natural land surfaces
+             surf_lsm_h(l)%start_index(j,i) = start_index_lsm_h(l)
+             surf_lsm_h(l)%end_index(j,i)   = surf_lsm_h(l)%start_index(j,i) + num_lsm_h_kji(l) - 1
+             start_index_lsm_h(l)           = surf_lsm_h(l)%end_index(j,i) + 1
+!
+!--          Horizontal urban surfaces
+             surf_usm_h(l)%start_index(j,i) = start_index_usm_h(l)
+             surf_usm_h(l)%end_index(j,i)   = surf_usm_h(l)%start_index(j,i) + num_usm_h_kji(l) - 1
+             start_index_usm_h(l)           = surf_usm_h(l)%end_index(j,i) + 1
+          ENDDO
 
 !
@@ -3138,125 +3173,125 @@
              ENDDO
 
-             IF ( l == 0 )  THEN
-                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
-                   IF ( ALLOCATED( surf_lsm_h%us ) )  surf_h(0)%us(mm(0)) = surf_lsm_h%us(m)
-                   IF ( ALLOCATED( surf_lsm_h%ts ) )  surf_h(0)%ts(mm(0)) = surf_lsm_h%ts(m)
-                   IF ( ALLOCATED( surf_lsm_h%qs ) )  surf_h(0)%qs(mm(0)) = surf_lsm_h%qs(m)
-                   IF ( ALLOCATED( surf_lsm_h%ss ) )  surf_h(0)%ss(mm(0)) = surf_lsm_h%ss(m)
-                   IF ( ALLOCATED( surf_lsm_h%qcs ) )  surf_h(0)%qcs(mm(0)) = surf_lsm_h%qcs(m)
-                   IF ( ALLOCATED( surf_lsm_h%ncs ) )  surf_h(0)%ncs(mm(0)) = surf_lsm_h%ncs(m)
-                   IF ( ALLOCATED( surf_lsm_h%qis ) )  surf_h(0)%qis(mm(0)) = surf_lsm_h%qis(m)
-                   IF ( ALLOCATED( surf_lsm_h%nis ) )  surf_h(0)%nis(mm(0)) = surf_lsm_h%nis(m)
-                   IF ( ALLOCATED( surf_lsm_h%qrs ) )  surf_h(0)%qrs(mm(0)) = surf_lsm_h%qrs(m)
-                   IF ( ALLOCATED( surf_lsm_h%nrs ) )  surf_h(0)%nrs(mm(0)) = surf_lsm_h%nrs(m)
-                   IF ( ALLOCATED( surf_lsm_h%ol ) )  surf_h(0)%ol(mm(0)) = surf_lsm_h%ol(m)
-                   IF ( ALLOCATED( surf_lsm_h%rib ) )  surf_h(0)%rib(mm(0)) = surf_lsm_h%rib(m)
-                   IF ( ALLOCATED( surf_lsm_h%pt_surface ) )                                       &
-                      surf_h(l)%pt_surface(mm(l)) = surf_lsm_h%pt_surface(m)
-                   IF ( ALLOCATED( surf_def_h(l)%q_surface ) )                                     &
-                      surf_h(l)%q_surface(mm(l)) = surf_lsm_h%q_surface(m)
-                   IF ( ALLOCATED( surf_def_h(l)%vpt_surface ) )                                   &
-                      surf_h(l)%vpt_surface(mm(l)) = surf_lsm_h%vpt_surface(m)
-                   IF ( ALLOCATED( surf_lsm_h%usws ) )  surf_h(0)%usws(mm(0)) = surf_lsm_h%usws(m)
-                   IF ( ALLOCATED( surf_lsm_h%vsws ) )  surf_h(0)%vsws(mm(0)) = surf_lsm_h%vsws(m)
-                   IF ( ALLOCATED( surf_lsm_h%shf ) )  surf_h(0)%shf(mm(0)) = surf_lsm_h%shf(m)
-                   IF ( ALLOCATED( surf_lsm_h%qsws ) )  surf_h(0)%qsws(mm(0)) = surf_lsm_h%qsws(m)
-                   IF ( ALLOCATED( surf_lsm_h%ssws ) )  surf_h(0)%ssws(mm(0)) = surf_lsm_h%ssws(m)
-                   IF ( ALLOCATED( surf_lsm_h%css ) )  THEN
+             IF ( l < 2 )  THEN
+                DO  m = surf_lsm_h(l)%start_index(j,i), surf_lsm_h(l)%end_index(j,i)
+                   IF ( ALLOCATED( surf_lsm_h(l)%us ) )  surf_h(0)%us(mm(0)) = surf_lsm_h(l)%us(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%ts ) )  surf_h(0)%ts(mm(0)) = surf_lsm_h(l)%ts(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%qs ) )  surf_h(0)%qs(mm(0)) = surf_lsm_h(l)%qs(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%ss ) )  surf_h(0)%ss(mm(0)) = surf_lsm_h(l)%ss(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%qcs ) )  surf_h(0)%qcs(mm(0)) = surf_lsm_h(l)%qcs(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%ncs ) )  surf_h(0)%ncs(mm(0)) = surf_lsm_h(l)%ncs(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%qis ) )  surf_h(0)%qis(mm(0)) = surf_lsm_h(l)%qis(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%nis ) )  surf_h(0)%nis(mm(0)) = surf_lsm_h(l)%nis(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%qrs ) )  surf_h(0)%qrs(mm(0)) = surf_lsm_h(l)%qrs(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%nrs ) )  surf_h(0)%nrs(mm(0)) = surf_lsm_h(l)%nrs(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%ol ) )  surf_h(0)%ol(mm(0)) = surf_lsm_h(l)%ol(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%rib ) )  surf_h(0)%rib(mm(0)) = surf_lsm_h(l)%rib(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%pt_surface ) )                                       &
+                      surf_h(l)%pt_surface(mm(l)) = surf_lsm_h(l)%pt_surface(m)
+                   IF ( ALLOCATED( surf_def_h(l)%q_surface ) )                                        &
+                      surf_h(l)%q_surface(mm(l)) = surf_lsm_h(l)%q_surface(m)
+                   IF ( ALLOCATED( surf_def_h(l)%vpt_surface ) )                                      &
+                      surf_h(l)%vpt_surface(mm(l)) = surf_lsm_h(l)%vpt_surface(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%usws ) )  surf_h(0)%usws(mm(0)) = surf_lsm_h(l)%usws(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%vsws ) )  surf_h(0)%vsws(mm(0)) = surf_lsm_h(l)%vsws(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%shf ) )  surf_h(0)%shf(mm(0)) = surf_lsm_h(l)%shf(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%qsws ) )  surf_h(0)%qsws(mm(0)) = surf_lsm_h(l)%qsws(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%ssws ) )  surf_h(0)%ssws(mm(0)) = surf_lsm_h(l)%ssws(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%css ) )  THEN
                       DO  lsp = 1, nvar
-                         surf_h(0)%css(lsp,mm(0)) = surf_lsm_h%css(lsp,m)
+                         surf_h(0)%css(lsp,mm(0)) = surf_lsm_h(l)%css(lsp,m)
                       ENDDO
                    ENDIF
-                   IF ( ALLOCATED( surf_lsm_h%cssws ) )  THEN
+                   IF ( ALLOCATED( surf_lsm_h(l)%cssws ) )  THEN
                       DO  lsp = 1, nvar
-                         surf_h(0)%cssws(lsp,mm(0)) = surf_lsm_h%cssws(lsp,m)
+                         surf_h(0)%cssws(lsp,mm(0)) = surf_lsm_h(l)%cssws(lsp,m)
                       ENDDO
                    ENDIF
-                   IF ( ALLOCATED( surf_lsm_h%qcsws ) )                                            &
-                      surf_h(0)%qcsws(mm(0)) = surf_lsm_h%qcsws(m)
-                   IF ( ALLOCATED( surf_lsm_h%qisws ) )                                            &
-                      surf_h(0)%qisws(mm(0)) = surf_lsm_h%qisws(m)
-                   IF ( ALLOCATED( surf_lsm_h%qrsws ) )                                            &
-                      surf_h(0)%qrsws(mm(0)) = surf_lsm_h%qrsws(m)
-                   IF ( ALLOCATED( surf_lsm_h%ncsws ) )                                            &
-                      surf_h(0)%ncsws(mm(0)) = surf_lsm_h%ncsws(m)
-                   IF ( ALLOCATED( surf_lsm_h%nisws ) )                                            &
-                      surf_h(0)%nisws(mm(0)) = surf_lsm_h%nisws(m)
-                   IF ( ALLOCATED( surf_lsm_h%nrsws ) )                                            &
-                      surf_h(0)%nrsws(mm(0)) = surf_lsm_h%nrsws(m)
-                   IF ( ALLOCATED( surf_lsm_h%sasws ) )                                            &
-                     surf_h(0)%sasws(mm(0)) = surf_lsm_h%sasws(m)
-
-                   mm(0) = mm(0) + 1
+                   IF ( ALLOCATED( surf_lsm_h(l)%qcsws ) )                                            &
+                      surf_h(0)%qcsws(mm(0)) = surf_lsm_h(l)%qcsws(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%qisws ) )                                            &
+                      surf_h(0)%qisws(mm(0)) = surf_lsm_h(l)%qisws(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%qrsws ) )                                            &
+                      surf_h(0)%qrsws(mm(0)) = surf_lsm_h(l)%qrsws(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%ncsws ) )                                            &
+                      surf_h(0)%ncsws(mm(0)) = surf_lsm_h(l)%ncsws(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%nisws ) )                                            &
+                      surf_h(0)%nisws(mm(0)) = surf_lsm_h(l)%nisws(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%nrsws ) )                                            &
+                      surf_h(0)%nrsws(mm(0)) = surf_lsm_h(l)%nrsws(m)
+                   IF ( ALLOCATED( surf_lsm_h(l)%sasws ) )                                            &
+                     surf_h(0)%sasws(mm(0)) = surf_lsm_h(l)%sasws(m)
+
+                   mm(l) = mm(l) + 1
 
                 ENDDO
 
-                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
-                   IF ( ALLOCATED( surf_usm_h%us ) )                                               &
-                      surf_h(0)%us(mm(0)) = surf_usm_h%us(m)
-                   IF ( ALLOCATED( surf_usm_h%ts ) )                                               &
-                      surf_h(0)%ts(mm(0)) = surf_usm_h%ts(m)
-                   IF ( ALLOCATED( surf_usm_h%qs ) )                                               &
-                      surf_h(0)%qs(mm(0)) = surf_usm_h%qs(m)
-                   IF ( ALLOCATED( surf_usm_h%ss ) )                                               &
-                      surf_h(0)%ss(mm(0)) = surf_usm_h%ss(m)
-                   IF ( ALLOCATED( surf_usm_h%qcs ) )                                              &
-                      surf_h(0)%qcs(mm(0)) = surf_usm_h%qcs(m)
-                   IF ( ALLOCATED( surf_usm_h%ncs ) )                                              &
-                      surf_h(0)%ncs(mm(0)) = surf_usm_h%ncs(m)
-                   IF ( ALLOCATED( surf_usm_h%qis ) )                                              &
-                      surf_h(0)%qis(mm(0)) = surf_usm_h%qis(m)
-                   IF ( ALLOCATED( surf_usm_h%nis ) )                                              &
-                      surf_h(0)%nis(mm(0)) = surf_usm_h%nis(m)
-                   IF ( ALLOCATED( surf_usm_h%qrs ) )                                              &
-                      surf_h(0)%qrs(mm(0)) = surf_usm_h%qrs(m)
-                   IF ( ALLOCATED( surf_usm_h%nrs ) )                                              &
-                      surf_h(0)%nrs(mm(0)) = surf_usm_h%nrs(m)
-                   IF ( ALLOCATED( surf_usm_h%ol ) )                                               &
-                      surf_h(0)%ol(mm(0)) = surf_usm_h%ol(m)
-                   IF ( ALLOCATED( surf_usm_h%rib ) )                                              &
-                      surf_h(0)%rib(mm(0)) = surf_usm_h%rib(m)
-                   IF ( ALLOCATED( surf_usm_h%pt_surface ) )                                       &
-                      surf_h(l)%pt_surface(mm(l)) = surf_usm_h%pt_surface(m)
-                    IF ( ALLOCATED( surf_usm_h%q_surface ) )                                       &
-                      surf_h(l)%q_surface(mm(l)) = surf_usm_h%q_surface(m)
-                   IF ( ALLOCATED( surf_usm_h%vpt_surface ) )                                      &
-                      surf_h(l)%vpt_surface(mm(l)) = surf_usm_h%vpt_surface(m)
-                   IF ( ALLOCATED( surf_usm_h%usws ) )                                             &
-                      surf_h(0)%usws(mm(0)) = surf_usm_h%usws(m)
-                   IF ( ALLOCATED( surf_usm_h%vsws ) )                                             &
-                      surf_h(0)%vsws(mm(0)) = surf_usm_h%vsws(m)
-                   IF ( ALLOCATED( surf_usm_h%shf ) )                                              &
-                      surf_h(0)%shf(mm(0)) = surf_usm_h%shf(m)
-                   IF ( ALLOCATED( surf_usm_h%qsws ) )                                             &
-                      surf_h(0)%qsws(mm(0)) = surf_usm_h%qsws(m)
-                   IF ( ALLOCATED( surf_usm_h%ssws ) )                                             &
-                      surf_h(0)%ssws(mm(0)) = surf_usm_h%ssws(m)
-                   IF ( ALLOCATED( surf_usm_h%css ) )  THEN
+                DO  m = surf_usm_h(l)%start_index(j,i), surf_usm_h(l)%end_index(j,i)
+                   IF ( ALLOCATED( surf_usm_h(l)%us ) )                                               &
+                      surf_h(0)%us(mm(0)) = surf_usm_h(l)%us(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%ts ) )                                               &
+                      surf_h(0)%ts(mm(0)) = surf_usm_h(l)%ts(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%qs ) )                                               &
+                      surf_h(0)%qs(mm(0)) = surf_usm_h(l)%qs(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%ss ) )                                               &
+                      surf_h(0)%ss(mm(0)) = surf_usm_h(l)%ss(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%qcs ) )                                              &
+                      surf_h(0)%qcs(mm(0)) = surf_usm_h(l)%qcs(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%ncs ) )                                              &
+                      surf_h(0)%ncs(mm(0)) = surf_usm_h(l)%ncs(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%qis ) )                                              &
+                      surf_h(0)%qis(mm(0)) = surf_usm_h(l)%qis(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%nis ) )                                              &
+                      surf_h(0)%nis(mm(0)) = surf_usm_h(l)%nis(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%qrs ) )                                              &
+                      surf_h(0)%qrs(mm(0)) = surf_usm_h(l)%qrs(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%nrs ) )                                              &
+                      surf_h(0)%nrs(mm(0)) = surf_usm_h(l)%nrs(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%ol ) )                                               &
+                      surf_h(0)%ol(mm(0)) = surf_usm_h(l)%ol(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%rib ) )                                              &
+                      surf_h(0)%rib(mm(0)) = surf_usm_h(l)%rib(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%pt_surface ) )                                       &
+                      surf_h(l)%pt_surface(mm(l)) = surf_usm_h(l)%pt_surface(m)
+                    IF ( ALLOCATED( surf_usm_h(l)%q_surface ) )                                       &
+                      surf_h(l)%q_surface(mm(l)) = surf_usm_h(l)%q_surface(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%vpt_surface ) )                                      &
+                      surf_h(l)%vpt_surface(mm(l)) = surf_usm_h(l)%vpt_surface(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%usws ) )                                             &
+                      surf_h(0)%usws(mm(0)) = surf_usm_h(l)%usws(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%vsws ) )                                             &
+                      surf_h(0)%vsws(mm(0)) = surf_usm_h(l)%vsws(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%shf ) )                                              &
+                      surf_h(0)%shf(mm(0)) = surf_usm_h(l)%shf(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%qsws ) )                                             &
+                      surf_h(0)%qsws(mm(0)) = surf_usm_h(l)%qsws(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%ssws ) )                                             &
+                      surf_h(0)%ssws(mm(0)) = surf_usm_h(l)%ssws(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%css ) )  THEN
                       DO lsp = 1, nvar
-                         surf_h(0)%css(lsp,mm(0)) = surf_usm_h%css(lsp,m)
+                         surf_h(0)%css(lsp,mm(0)) = surf_usm_h(l)%css(lsp,m)
                       ENDDO
                    ENDIF
-                   IF ( ALLOCATED( surf_usm_h%cssws ) )  THEN
+                   IF ( ALLOCATED( surf_usm_h(l)%cssws ) )  THEN
                       DO lsp = 1, nvar
-                         surf_h(0)%cssws(lsp,mm(0)) = surf_usm_h%cssws(lsp,m)
+                         surf_h(0)%cssws(lsp,mm(0)) = surf_usm_h(l)%cssws(lsp,m)
                       ENDDO
                    ENDIF
-                   IF ( ALLOCATED( surf_usm_h%qcsws ) )                                            &
-                      surf_h(0)%qcsws(mm(0)) = surf_usm_h%qcsws(m)
-                   IF ( ALLOCATED( surf_usm_h%qisws ) )                                            &
-                      surf_h(0)%qisws(mm(0)) = surf_usm_h%qisws(m)
-                   IF ( ALLOCATED( surf_usm_h%qrsws ) )                                            &
-                      surf_h(0)%qrsws(mm(0)) = surf_usm_h%qrsws(m)
-                   IF ( ALLOCATED( surf_usm_h%ncsws ) )                                            &
-                      surf_h(0)%ncsws(mm(0))   = surf_usm_h%ncsws(m)
-                   IF ( ALLOCATED( surf_usm_h%nrsws ) )                                            &
-                      surf_h(0)%nrsws(mm(0)) = surf_usm_h%nrsws(m)
-                   IF ( ALLOCATED( surf_usm_h%nisws ) )                                            &
-                      surf_h(0)%nisws(mm(0)) = surf_usm_h%nisws(m)
-                   IF ( ALLOCATED( surf_usm_h%sasws ) )                                            &
-                     surf_h(0)%sasws(mm(0)) = surf_usm_h%sasws(m)
-
-                   mm(0) = mm(0) + 1
+                   IF ( ALLOCATED( surf_usm_h(l)%qcsws ) )                                            &
+                      surf_h(0)%qcsws(mm(0)) = surf_usm_h(l)%qcsws(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%qisws ) )                                            &
+                      surf_h(0)%qisws(mm(0)) = surf_usm_h(l)%qisws(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%qrsws ) )                                            &
+                      surf_h(0)%qrsws(mm(0)) = surf_usm_h(l)%qrsws(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%ncsws ) )                                            &
+                      surf_h(0)%ncsws(mm(0))   = surf_usm_h(l)%ncsws(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%nrsws ) )                                            &
+                      surf_h(0)%nrsws(mm(0)) = surf_usm_h(l)%nrsws(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%nisws ) )                                            &
+                      surf_h(0)%nisws(mm(0)) = surf_usm_h(l)%nisws(m)
+                   IF ( ALLOCATED( surf_usm_h(l)%sasws ) )                                            &
+                     surf_h(0)%sasws(mm(0)) = surf_usm_h(l)%sasws(m)
+
+                   mm(l) = mm(l) + 1
 
                 ENDDO
@@ -3280,9 +3315,9 @@
                 surf_h(l)%end_index(j,i) = surf_h(l)%end_index(j,i) + 1
              ENDDO
-             IF ( l == 0 )  THEN
-                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
+             IF ( l < 2 )  THEN
+                DO  m = surf_lsm_h(l)%start_index(j,i), surf_lsm_h(l)%end_index(j,i)
                    surf_h(l)%end_index(j,i) = surf_h(l)%end_index(j,i) + 1
                 ENDDO
-                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
+                DO  m = surf_usm_h(l)%start_index(j,i), surf_usm_h(l)%end_index(j,i)
                    surf_h(l)%end_index(j,i) = surf_h(l)%end_index(j,i) + 1
                 ENDDO
@@ -4977,8 +5012,11 @@
 !--          Determine type of surface element, i.e. default, natural, urban, at current grid point.
              surf_match_def  = surf_def_h(l)%end_index(jc,ic) >= surf_def_h(l)%start_index(jc,ic)
-             surf_match_lsm  = ( surf_lsm_h%end_index(jc,ic)  >= surf_lsm_h%start_index(jc,ic) )   &
-                                .AND.  l == 0
-             surf_match_usm  = ( surf_usm_h%end_index(jc,ic)  >= surf_usm_h%start_index(jc,ic) )   &
-                                 .AND.  l == 0
+             IF ( l < 2 ) THEN
+                surf_match_lsm  = surf_lsm_h(l)%end_index(jc,ic) >= surf_lsm_h(l)%start_index(jc,ic)
+                surf_match_usm  = surf_usm_h(l)%end_index(jc,ic) >= surf_usm_h(l)%start_index(jc,ic)
+             ELSE
+                surf_match_lsm  = .FALSE.
+                surf_match_usm  = .FALSE.
+             ENDIF
 !
 !--          Write restart data onto default-type surfaces if required.
@@ -5002,8 +5040,8 @@
 !--          surfaces at (j,i). An example would be bridges.
              IF ( surf_match_lsm )  THEN
-                mm = surf_lsm_h%start_index(jc,ic)
+                mm = surf_lsm_h(l)%start_index(jc,ic)
                 DO  m = surf_h(l)%start_index(j,i), surf_h(l)%end_index(j,i)
-                   IF ( surf_lsm_h%end_index(jc,ic) >= mm )                                        &
-                      CALL restore_surface_elements( surf_lsm_h, mm, surf_h(l), m )
+                   IF ( surf_lsm_h(l)%end_index(jc,ic) >= mm )                                        &
+                      CALL restore_surface_elements( surf_lsm_h(l), mm, surf_h(l), m )
                    mm = mm + 1
                 ENDDO
@@ -5012,8 +5050,8 @@
 !--          Same for urban-type surfaces
              IF ( surf_match_usm )  THEN
-                mm = surf_usm_h%start_index(jc,ic)
+                mm = surf_usm_h(l)%start_index(jc,ic)
                 DO  m = surf_h(l)%start_index(j,i), surf_h(l)%end_index(j,i)
-                   IF ( surf_usm_h%end_index(jc,ic) >= mm )                                        &
-                      CALL restore_surface_elements( surf_usm_h, mm, surf_h(l), m )
+                   IF ( surf_usm_h(l)%end_index(jc,ic) >= mm )                                        &
+                      CALL restore_surface_elements( surf_usm_h(l), mm, surf_h(l), m )
                    mm = mm + 1
                 ENDDO
@@ -5606,12 +5644,12 @@
     DO  i = nxl, nxr
        DO  j = nys, nyn
-          surf_match_lsm  = surf_lsm_h%end_index(j,i) >= surf_lsm_h%start_index(j,i)
-          surf_match_usm  = surf_usm_h%end_index(j,i) >= surf_usm_h%start_index(j,i)
+          surf_match_lsm  = surf_lsm_h(l)%end_index(j,i) >= surf_lsm_h(l)%start_index(j,i)
+          surf_match_usm  = surf_usm_h(l)%end_index(j,i) >= surf_usm_h(l)%start_index(j,i)
 
           IF ( surf_match_lsm )  THEN
-             mm = surf_lsm_h%start_index(j,i)
+             mm = surf_lsm_h(l)%start_index(j,i)
              DO  m = surf_h(0)%start_index(j,i), surf_h(0)%end_index(j,i)
-                IF ( surf_lsm_h%end_index(j,i) >= mm )                                             &
-                   CALL restore_surface_elements( surf_lsm_h, mm, surf_h(0), m )
+                IF ( surf_lsm_h(l)%end_index(j,i) >= mm )                                             &
+                   CALL restore_surface_elements( surf_lsm_h(l), mm, surf_h(0), m )
                 mm = mm + 1
              ENDDO
@@ -5619,8 +5657,8 @@
 
           IF ( surf_match_usm )  THEN
-             mm = surf_usm_h%start_index(j,i)
+             mm = surf_usm_h(l)%start_index(j,i)
              DO  m = surf_h(0)%start_index(j,i), surf_h(0)%end_index(j,i)
-                IF ( surf_usm_h%end_index(j,i) >= mm )                                             &
-                   CALL restore_surface_elements( surf_usm_h, mm, surf_h(0), m )
+                IF ( surf_usm_h(l)%end_index(j,i) >= mm )                                             &
+                   CALL restore_surface_elements( surf_usm_h(l), mm, surf_h(0), m )
                 mm = mm + 1
              ENDDO
@@ -5802,6 +5840,6 @@
 !
 !-- Horizontal surfaces
-    ns_h_on_file(0) = surf_def_h(0)%ns + surf_lsm_h%ns + surf_usm_h%ns
-    ns_h_on_file(1) = surf_def_h(1)%ns
+    ns_h_on_file(0) = surf_def_h(0)%ns + surf_lsm_h(0)%ns + surf_usm_h(0)%ns
+    ns_h_on_file(1) = surf_def_h(1)%ns + surf_lsm_h(1)%ns + surf_usm_h(1)%ns
     ns_h_on_file(2) = surf_def_h(2)%ns
 !
