Index: /palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_3d.f90	(revision 3813)
+++ /palm/trunk/SOURCE/data_output_3d.f90	(revision 3814)
@@ -298,5 +298,5 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  nzub, nzut
+        ONLY:  nz_urban_b, nz_urban_t
 
     USE turbulence_closure_mod,                                                &
@@ -723,6 +723,6 @@
 !
 !--             For urban model quantities, it is required to re-allocate local_pf
-                nzb_do = nzub
-                nzt_do = nzut
+                nzb_do = nz_urban_b
+                nzt_do = nz_urban_t
 
                 DEALLOCATE ( local_pf )
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3813)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3814)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! Rename exported nzu, nzp and related variables due to name conflict
+! 
+! 3771 2019-02-28 12:19:33Z raasch
 ! rrtmg preprocessor for directives moved/added, save attribute added to temporary
 ! pointers to avoid compiler warnings about outlived pointer targets,
@@ -584,5 +587,5 @@
 
     USE basic_constants_and_equations_mod,                                     &
-        ONLY:  c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant,      &
+        ONLY:  c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, sigma_sb, &
                barometric_formula
 
@@ -708,8 +711,6 @@
                                                          /)
 
-    REAL(wp), PARAMETER :: sigma_sb       = 5.67037321E-8_wp !< Stefan-Boltzmann constant
-
-    INTEGER(iwp) :: albedo_type  = 9999999, & !< Albedo surface type
-                    dots_rad     = 0          !< starting index for timeseries output 
+    INTEGER(iwp) :: albedo_type  = 9999999_iwp, &     !< Albedo surface type
+                    dots_rad     = 0_iwp              !< starting index for timeseries output 
 
     LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
@@ -912,8 +913,9 @@
 !
 !-- Parameters of urban and land surface models
-    INTEGER(iwp)                                   ::  nzu                                !< number of layers of urban surface (will be calculated)
-    INTEGER(iwp)                                   ::  nzp                                !< number of layers of plant canopy (will be calculated)
-    INTEGER(iwp)                                   ::  nzub,nzut                          !< bottom and top layer of urban surface (will be calculated)
-    INTEGER(iwp)                                   ::  nzpt                               !< top layer of plant canopy (will be calculated)
+    INTEGER(iwp)                                   ::  nz_urban                           !< number of layers of urban surface (will be calculated)
+    INTEGER(iwp)                                   ::  nz_plant                           !< number of layers of plant canopy (will be calculated)
+    INTEGER(iwp)                                   ::  nz_urban_b                         !< bottom layer of urban surface (will be calculated)
+    INTEGER(iwp)                                   ::  nz_urban_t                         !< top layer of urban surface (will be calculated)
+    INTEGER(iwp)                                   ::  nz_plant_t                         !< top layer of plant canopy (will be calculated)
 !-- parameters of urban and land surface models
     INTEGER(iwp), PARAMETER                        ::  nzut_free = 3                      !< number of free layers above top of of topography
@@ -1104,5 +1106,5 @@
     INTEGER(iwp)                                   ::  nmrtfa           !< dimmension of array allocated for storage of mrt
     INTEGER(iwp)                                   ::  msvf, mcsf, mmrtf!< mod for swapping the growing array
-    INTEGER(iwp), PARAMETER                        ::  gasize = 100000  !< initial size of growing arrays
+    INTEGER(iwp), PARAMETER                        ::  gasize = 100000_iwp  !< initial size of growing arrays
     REAL(wp), PARAMETER                            ::  grow_factor = 1.4_wp !< growth factor of growing arrays
     INTEGER(iwp)                                   ::  nsvfl            !< number of svf for local processor
@@ -1298,5 +1300,5 @@
            iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
            iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
-           nsurf_type, nzub, nzut, nzu, pch, nsurf,                            &
+           nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf,                 &
            idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
            radiation_interactions, startwall, startland, endland, endwall,     &
@@ -2843,6 +2845,6 @@
           IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1   = 0.0_wp
 
-          pt1_l = SUM( pt(nzut,nys:nyn,nxl:nxr) )
-          IF ( bulk_cloud_model  .OR.  cloud_droplets  )  ql1_l = SUM( ql(nzut,nys:nyn,nxl:nxr) )
+          pt1_l = SUM( pt(nz_urban_t,nys:nyn,nxl:nxr) )
+          IF ( bulk_cloud_model  .OR.  cloud_droplets  )  ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) )
 
 #if defined( __parallel )      
@@ -2866,5 +2868,5 @@
 #endif
 
-          IF ( bulk_cloud_model  .OR.  cloud_droplets  )  pt1 = pt1 + lv_d_cp / exner(nzut) * ql1
+          IF ( bulk_cloud_model  .OR.  cloud_droplets  )  pt1 = pt1 + lv_d_cp / exner(nz_urban_t) * ql1
 !
 !--       Finally, divide by number of grid points
@@ -2905,5 +2907,5 @@
              IF ( average_radiation ) THEN       
 
-                k = nzut
+                k = nz_urban_t
 
                 surf%rad_sw_in  = solar_constant * sky_trans * zenith(0)
@@ -3023,6 +3025,6 @@
           IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1   = 0.0_wp
 
-          pt1_l = SUM( pt(nzut,nys:nyn,nxl:nxr) )
-          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1_l = SUM( ql(nzut,nys:nyn,nxl:nxr) )
+          pt1_l = SUM( pt(nz_urban_t,nys:nyn,nxl:nxr) )
+          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) )
 
 #if defined( __parallel )      
@@ -3044,5 +3046,5 @@
           IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1 = ql1_l
 #endif
-          IF ( bulk_cloud_model  .OR.  cloud_droplets )  pt1 = pt1 + lv_d_cp / exner(nzut+1) * ql1
+          IF ( bulk_cloud_model  .OR.  cloud_droplets )  pt1 = pt1 + lv_d_cp / exner(nz_urban_t+1) * ql1
 !
 !--       Finally, divide by number of grid points
@@ -3086,5 +3088,5 @@
                 surf%rad_net = net_radiation
 
-                surf%rad_lw_in  = 0.8_wp * sigma_sb * (pt1 * exner(nzut+1))**4
+                surf%rad_lw_in  = 0.8_wp * sigma_sb * (pt1 * exner(nz_urban_t+1))**4
 
                 surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4   &
@@ -4979,5 +4981,5 @@
      REAL(wp), DIMENSION(3)            :: sunorig_grid       !< grid squashed solar direction unit vector (zyx)
      REAL(wp), DIMENSION(0:nsurf_type) :: costheta           !< direct irradiance factor of solar angle
-     REAL(wp), DIMENSION(nzub:nzut)    :: pchf_prep          !< precalculated factor for canopy temperature tendency
+     REAL(wp), DIMENSION(nz_urban_b:nz_urban_t)    :: pchf_prep          !< precalculated factor for canopy temperature tendency
      REAL(wp), PARAMETER               :: alpha = 0._wp      !< grid rotation (TODO: add to namelist or remove)
      REAL(wp)                          :: pc_box_area, pc_abs_frac, pc_abs_eff
@@ -5002,6 +5004,6 @@
 
      IF ( plant_canopy )  THEN
-         pchf_prep(:) = r_d * exner(nzub:nzut)                                 &
-                     / (c_p * hyp(nzub:nzut) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T)
+         pchf_prep(:) = r_d * exner(nz_urban_b:nz_urban_t)                                 &
+                     / (c_p * hyp(nz_urban_b:nz_urban_t) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T)
      ENDIF
 
@@ -5923,5 +5925,5 @@
         ENDDO
 !
-!--    Find nzub, nzut, nzu via wall_flag_0 array (nzb_s_inner will be
+!--    Find nz_urban_b, nz_urban_t, nz_urban via wall_flag_0 array (nzb_s_inner will be
 !--    removed later). The following contruct finds the lowest / largest index
 !--    for any upward-facing wall (see bit 12).
@@ -5976,23 +5978,23 @@
 
 #if defined( __parallel )
-       CALL MPI_AllReduce(nzubl, nzub, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )
+       CALL MPI_AllReduce(nzubl, nz_urban_b, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )
        IF ( ierr /= 0 ) THEN
-           WRITE(9,*) 'Error MPI_AllReduce11:', ierr, nzubl, nzub
+           WRITE(9,*) 'Error MPI_AllReduce11:', ierr, nzubl, nz_urban_b
            FLUSH(9)
        ENDIF
-       CALL MPI_AllReduce(nzutl, nzut, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
+       CALL MPI_AllReduce(nzutl, nz_urban_t, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
        IF ( ierr /= 0 ) THEN
-           WRITE(9,*) 'Error MPI_AllReduce12:', ierr, nzutl, nzut
+           WRITE(9,*) 'Error MPI_AllReduce12:', ierr, nzutl, nz_urban_t
            FLUSH(9)
        ENDIF
-       CALL MPI_AllReduce(nzptl, nzpt, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
+       CALL MPI_AllReduce(nzptl, nz_plant_t, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
        IF ( ierr /= 0 ) THEN
-           WRITE(9,*) 'Error MPI_AllReduce13:', ierr, nzptl, nzpt
+           WRITE(9,*) 'Error MPI_AllReduce13:', ierr, nzptl, nz_plant_t
            FLUSH(9)
        ENDIF
 #else
-       nzub = nzubl
-       nzut = nzutl
-       nzpt = nzptl
+       nz_urban_b = nzubl
+       nz_urban_t = nzutl
+       nz_plant_t = nzptl
 #endif
 !
@@ -6002,17 +6004,17 @@
 !--    are active. ABS (...) is required because the default value of 
 !--    dz_stretch_level_start is -9999999.9_wp (negative). 
-       IF ( ABS( dz_stretch_level_start(1) ) <= zw(nzut) ) THEN
+       IF ( ABS( dz_stretch_level_start(1) ) <= zw(nz_urban_t) ) THEN
           WRITE( message_string, * ) 'The lowest level where vertical ',       &
                                      'stretching is applied have to be ',      &
-                                     'greater than ', zw(nzut)
+                                     'greater than ', zw(nz_urban_t)
           CALL message( 'radiation_interaction_init', 'PA0496', 1, 2, 0, 6, 0 )
        ENDIF 
 !
 !--    global number of urban and plant layers
-       nzu = nzut - nzub + 1
-       nzp = nzpt - nzub + 1
+       nz_urban = nz_urban_t - nz_urban_b + 1
+       nz_plant = nz_plant_t - nz_urban_b + 1
 !
 !--    check max_raytracing_dist relative to urban surface layer height
-       mrl = 2.0_wp * nzu * dz(1)
+       mrl = 2.0_wp * nz_urban * dz(1)
 !--    set max_raytracing_dist to double the urban surface layer height, if not set
        IF ( max_raytracing_dist == -999.0_wp ) THEN
@@ -6063,5 +6065,5 @@
        IF ( npcbl > 0 )  THEN
            ALLOCATE( pcbl(iz:ix, 1:npcbl) )
-           ALLOCATE( gridpcbl(nzub:nzpt,nys:nyn,nxl:nxr) )
+           ALLOCATE( gridpcbl(nz_urban_b:nz_plant_t,nys:nyn,nxl:nxr) )
            pcbl = -1
            gridpcbl(:,:,:) = 0
@@ -6125,10 +6127,10 @@
           ENDIF
 
-          CALL MPI_Win_allocate(INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nzu*nny*nnx, &
+          CALL MPI_Win_allocate(INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nz_urban*nny*nnx, &
                                     kind=MPI_ADDRESS_KIND), STORAGE_SIZE(1_iwp)/8,  &
                                 minfo, comm2d, gridsurf_rma_p, win_gridsurf, ierr)
           IF ( ierr /= 0 ) THEN
               WRITE(9,*) 'Error MPI_Win_allocate1:', ierr, &
-                 INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nzu*nny*nnx,kind=MPI_ADDRESS_KIND), &
+                 INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nz_urban*nny*nnx,kind=MPI_ADDRESS_KIND), &
                  STORAGE_SIZE(1_iwp)/8, win_gridsurf
               FLUSH(9)
@@ -6147,9 +6149,9 @@
 !--       pointer and then redirecting a multidimensional pointer to it works
 !--       fine.
-          CALL c_f_pointer(gridsurf_rma_p, gridsurf_rma, (/ nsurf_type_u*nzu*nny*nnx /))
-          gridsurf(0:nsurf_type_u-1, nzub:nzut, nys:nyn, nxl:nxr) =>                &
-                     gridsurf_rma(1:nsurf_type_u*nzu*nny*nnx)
+          CALL c_f_pointer(gridsurf_rma_p, gridsurf_rma, (/ nsurf_type_u*nz_urban*nny*nnx /))
+          gridsurf(0:nsurf_type_u-1, nz_urban_b:nz_urban_t, nys:nyn, nxl:nxr) =>                &
+                     gridsurf_rma(1:nsurf_type_u*nz_urban*nny*nnx)
 #else
-          ALLOCATE(gridsurf(0:nsurf_type_u-1,nzub:nzut,nys:nyn,nxl:nxr) )
+          ALLOCATE(gridsurf(0:nsurf_type_u-1,nz_urban_b:nz_urban_t,nys:nyn,nxl:nxr) )
 #endif
           gridsurf(:,:,:,:) = -999
@@ -6511,5 +6513,5 @@
         IF ( plant_canopy )  THEN
             ALLOCATE( plantt(0:(nx+1)*(ny+1)-1) )
-            maxboxesg = nx + ny + nzp + 1
+            maxboxesg = nx + ny + nz_plant + 1
             max_track_len = nx + ny + 1
 !--         temporary arrays storing values for csf calculation during raytracing
@@ -6564,5 +6566,5 @@
 !--             optimization of memory should be done
 !--             Argument X of function STORAGE_SIZE(X) needs arbitrary REAL(wp) value, set to 1.0_wp for now
-                size_lad_rma = STORAGE_SIZE(1.0_wp)/8*nnx*nny*nzp
+                size_lad_rma = STORAGE_SIZE(1.0_wp)/8*nnx*nny*nz_plant
                 CALL MPI_Win_allocate(size_lad_rma, STORAGE_SIZE(1.0_wp)/8, minfo, comm2d, &
                                         lad_s_rma_p, win_lad, ierr)
@@ -6572,16 +6574,16 @@
                     FLUSH(9)
                 ENDIF
-                CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzp*nny*nnx /))
-                sub_lad(nzub:nzpt, nys:nyn, nxl:nxr) => lad_s_rma(1:nzp*nny*nnx)
+                CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nz_plant*nny*nnx /))
+                sub_lad(nz_urban_b:nz_plant_t, nys:nyn, nxl:nxr) => lad_s_rma(1:nz_plant*nny*nnx)
             ELSE
-                ALLOCATE(sub_lad(nzub:nzpt, nys:nyn, nxl:nxr))
+                ALLOCATE(sub_lad(nz_urban_b:nz_plant_t, nys:nyn, nxl:nxr))
             ENDIF
 #else
             plantt = RESHAPE( pct(nys:nyn,nxl:nxr), (/(nx+1)*(ny+1)/) )
-            ALLOCATE(sub_lad(nzub:nzpt, nys:nyn, nxl:nxr))
+            ALLOCATE(sub_lad(nz_urban_b:nz_plant_t, nys:nyn, nxl:nxr))
 #endif
             plantt_max = MAXVAL(plantt)
-            ALLOCATE( rt2_track(2, max_track_len), rt2_track_lad(nzub:plantt_max, max_track_len), &
-                      rt2_track_dist(0:max_track_len), rt2_dist(plantt_max-nzub+2) )
+            ALLOCATE( rt2_track(2, max_track_len), rt2_track_lad(nz_urban_b:plantt_max, max_track_len), &
+                      rt2_track_dist(0:max_track_len), rt2_dist(plantt_max-nz_urban_b+2) )
 
             sub_lad(:,:,:) = 0._wp
@@ -6590,5 +6592,5 @@
                     k = get_topography_top_index_ji( j, i, 's' )
 
-                    sub_lad(k:nzpt, j, i) = lad_s(0:nzpt-k, j, i)
+                    sub_lad(k:nz_plant_t, j, i) = lad_s(0:nz_plant_t-k, j, i)
                 ENDDO
             ENDDO
@@ -6608,10 +6610,10 @@
                 
             ELSE
-                ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nzp-1) )
-                CALL MPI_AllGather( sub_lad, nnx*nny*nzp, MPI_REAL, &
-                                    sub_lad_g, nnx*nny*nzp, MPI_REAL, comm2d, ierr )
+                ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nz_plant-1) )
+                CALL MPI_AllGather( sub_lad, nnx*nny*nz_plant, MPI_REAL, &
+                                    sub_lad_g, nnx*nny*nz_plant, MPI_REAL, comm2d, ierr )
                 IF ( ierr /= 0 ) THEN
                     WRITE(9,*) 'Error MPI_AllGather3:', ierr, SIZE(sub_lad), &
-                                nnx*nny*nzp, SIZE(sub_lad_g), nnx*nny*nzp
+                                nnx*nny*nz_plant, SIZE(sub_lad_g), nnx*nny*nz_plant
                     FLUSH(9)
                 ENDIF
@@ -7551,5 +7553,5 @@
 #if defined( __parallel )
                         lad_ip(ncsb) = ip
-                        lad_disp(ncsb) = (box(3)-px*nnx)*(nny*nzp) + (box(2)-py*nny)*nzp + box(1)-nzub
+                        lad_disp(ncsb) = (box(3)-px*nnx)*(nny*nz_plant) + (box(2)-py*nny)*nz_plant + box(1)-nz_urban_b
 #endif
                     ENDIF
@@ -7595,5 +7597,5 @@
                     lad_s_target = lad_s_ray(i)
                 ELSE
-                    lad_s_target = sub_lad_g(lad_ip(i)*nnx*nny*nzp + lad_disp(i))
+                    lad_s_target = sub_lad_g(lad_ip(i)*nnx*nny*nz_plant + lad_disp(i))
                 ENDIF
 #else
@@ -7736,5 +7738,5 @@
          rt2_track_dist(0) = 0._wp
          rt2_track_lad(:,:) = 0._wp
-         nly = plantt_max - nzub + 1
+         nly = plantt_max - nz_urban_b + 1
       ENDIF
 
@@ -7871,5 +7873,5 @@
                IF ( plantt(ig) < lowest_lad )  CYCLE
 
-               wdisp = (rt2_track(2,i)-px*nnx)*(nny*nzp) + (rt2_track(1,i)-py*nny)*nzp + lowest_lad-nzub
+               wdisp = (rt2_track(2,i)-px*nnx)*(nny*nz_plant) + (rt2_track(1,i)-py*nny)*nz_plant + lowest_lad-nz_urban_b
                wcount = plantt(ig)-lowest_lad+1
                ! TODO send request ASAP - even during raytracing
@@ -7897,11 +7899,11 @@
                py = rt2_track(1,i)/nny
                ip = px*pdims(2)+py
-               ig = ip*nnx*nny*nzp + (rt2_track(2,i)-px*nnx)*(nny*nzp) + (rt2_track(1,i)-py*nny)*nzp
-               rt2_track_lad(nzub:plantt_max, i) = sub_lad_g(ig:ig+nly-1)
+               ig = ip*nnx*nny*nz_plant + (rt2_track(2,i)-px*nnx)*(nny*nz_plant) + (rt2_track(1,i)-py*nny)*nz_plant
+               rt2_track_lad(nz_urban_b:plantt_max, i) = sub_lad_g(ig:ig+nly-1)
             ENDDO
          ENDIF
 #else
          DO  i = 1, ntrack
-            rt2_track_lad(nzub:plantt_max, i) = sub_lad(rt2_track(1,i), rt2_track(2,i), nzub:plantt_max)
+            rt2_track_lad(nz_urban_b:plantt_max, i) = sub_lad(rt2_track(1,i), rt2_track(2,i), nz_urban_b:plantt_max)
          ENDDO
 #endif
@@ -7942,6 +7944,6 @@
 !--      Assert that we have space allocated for CSFs
 !--      
-         maxboxes = (ntrack + MAX(CEILING(origin(1)-.5_wp) - nzub,          &
-                                  nzut - CEILING(origin(1)-.5_wp))) * nrays
+         maxboxes = (ntrack + MAX(CEILING(origin(1)-.5_wp) - nz_urban_b,          &
+                                  nz_urban_t - CEILING(origin(1)-.5_wp))) * nrays
          IF ( ncsfl + maxboxes > ncsfla )  THEN
 !--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
@@ -7955,5 +7957,5 @@
 !--      Calculate transparencies and store new CSFs
 !--      
-         zbottom = REAL(nzub, wp) - .5_wp
+         zbottom = REAL(nz_urban_b, wp) - .5_wp
          ztop = REAL(plantt_max, wp) + .5_wp
 
@@ -8116,6 +8118,6 @@
          py = y / nny
          iproc = px * pdims(2) + py
-         target_displ = ((x-px*nnx) * nny + y - py*nny ) * nzu * nsurf_type_u +&
-                        ( z-nzub ) * nsurf_type_u + d
+         target_displ = ((x-px*nnx) * nny + y - py*nny ) * nz_urban * nsurf_type_u +&
+                        ( z-nz_urban_b ) * nsurf_type_u + d
 !
 !--      Send MPI_Get request to obtain index target_surfl(i)
Index: /palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3813)
+++ /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3814)
@@ -480,5 +480,5 @@
                iz, iy, ix,  nsurf, idsvf, ndsvf,                               &
                idcsf, ndcsf, kdcsf, pct,                                       &
-               nzub, nzut, unscheduled_radiation_calls
+               nz_urban_b, nz_urban_t, unscheduled_radiation_calls
 
     USE statistics,                                                            &
@@ -8585,5 +8585,5 @@
             DO i = nxl, nxr
                DO j = nys, nyn
-                  DO k = nzub, min(nzut,naheatlayers)
+                  DO k = nz_urban_b, min(nz_urban_t,naheatlayers)
                      IF ( k > get_topography_top_index_ji( j, i, 's' ) ) THEN
 !
