Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 1323)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 1353)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL constants provided with KIND-attribute 
 !
 ! Former revisions:
@@ -105,5 +105,5 @@
                f2_mg, f3_mg, l_grid, l_wall, zu, zw
         
-    USE control_parameters,                                                             &
+    USE control_parameters,                                                    &
         ONLY:  bc_lr, bc_ns, building_height, building_length_x,               &
                building_length_y, building_wall_left, building_wall_south,     &
@@ -201,13 +201,13 @@
 !
 !-- Allocate grid arrays
-    ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1), &
+    ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1),        &
               dzw(1:nzt+1), l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1) )
 
 !
 !-- Compute height of u-levels from constant grid length and dz stretch factors
-    IF ( dz == -1.0 )  THEN
+    IF ( dz == -1.0_wp )  THEN
        message_string = 'missing dz'
        CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 ) 
-    ELSEIF ( dz <= 0.0 )  THEN
+    ELSEIF ( dz <= 0.0_wp )  THEN
        WRITE( message_string, * ) 'dz=',dz,' <= 0.0'
        CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 )
@@ -223,10 +223,10 @@
 
        IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
-          zu(0) = 0.0
-      !    zu(0) = - dz * 0.5
+          zu(0) = 0.0_wp
+      !    zu(0) = - dz * 0.5_wp
        ELSE
-          zu(0) = - dz * 0.5
+          zu(0) = - dz * 0.5_wp
        ENDIF
-       zu(1) =   dz * 0.5
+       zu(1) =   dz * 0.5_wp
 
        dz_stretch_level_index = nzt+1
@@ -246,9 +246,9 @@
 !--    ground the first u- and w-level (k=0) are defined at same height (z=0). 
 !--    The top w-level is extrapolated linearly.
-       zw(0) = 0.0
+       zw(0) = 0.0_wp
        DO  k = 1, nzt
-          zw(k) = ( zu(k) + zu(k+1) ) * 0.5
-       ENDDO
-       zw(nzt+1) = zw(nzt) + 2.0 * ( zu(nzt+1) - zw(nzt) )
+          zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
+       ENDDO
+       zw(nzt+1) = zw(nzt) + 2.0_wp * ( zu(nzt+1) - zw(nzt) )
 
     ELSE
@@ -259,6 +259,6 @@
 !--    w-level are defined at same height, but staggered from the second level.
 !--    The second u-level (k=1) corresponds to the top of the Prandtl-layer.
-       zu(nzt+1) =   dz * 0.5
-       zu(nzt)   = - dz * 0.5
+       zu(nzt+1) =   dz * 0.5_wp
+       zu(nzt)   = - dz * 0.5_wp
 
        dz_stretch_level_index = 0
@@ -281,7 +281,7 @@
 !--    consistency, since w and all scalar variables are defined up tp nzt+1.
        zw(nzt+1) = dz
-       zw(nzt)   = 0.0
+       zw(nzt)   = 0.0_wp
        DO  k = 0, nzt
-          zw(k) = ( zu(k) + zu(k+1) ) * 0.5
+          zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
        ENDDO
 
@@ -299,11 +299,11 @@
     DO  k = 1, nzt+1
        dzu(k)  = zu(k) - zu(k-1)
-       ddzu(k) = 1.0 / dzu(k)
+       ddzu(k) = 1.0_wp / dzu(k)
        dzw(k)  = zw(k) - zw(k-1)
-       ddzw(k) = 1.0 / dzw(k)
+       ddzw(k) = 1.0_wp / dzw(k)
     ENDDO
 
     DO  k = 1, nzt
-       dd2zu(k) = 1.0 / ( dzu(k) + dzu(k+1) )
+       dd2zu(k) = 1.0_wp / ( dzu(k) + dzu(k+1) )
     ENDDO
     
@@ -340,6 +340,6 @@
        nzt_l = nzt
        DO  l = maximum_grid_level-1, 1, -1
-           dzu_mg(nzb+1,l) = 2.0 * dzu_mg(nzb+1,l+1)
-           dzw_mg(nzb+1,l) = 2.0 * dzw_mg(nzb+1,l+1)
+           dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1)
+           dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1)
            nzt_l = nzt_l / 2
            DO  k = 2, nzt_l+1
@@ -353,15 +353,15 @@
        dy_l  = dy
        DO  l = maximum_grid_level, 1, -1
-          ddx2_mg(l) = 1.0 / dx_l**2
-          ddy2_mg(l) = 1.0 / dy_l**2
+          ddx2_mg(l) = 1.0_wp / dx_l**2
+          ddy2_mg(l) = 1.0_wp / dy_l**2
           DO  k = nzb+1, nzt_l
-             f2_mg(k,l) = 1.0 / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
-             f3_mg(k,l) = 1.0 / ( dzu_mg(k,l)   * dzw_mg(k,l) )
-             f1_mg(k,l) = 2.0 * ( ddx2_mg(l) + ddy2_mg(l) ) + &
+             f2_mg(k,l) = 1.0_wp / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
+             f3_mg(k,l) = 1.0_wp / ( dzu_mg(k,l)   * dzw_mg(k,l) )
+             f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) + &
                           f2_mg(k,l) + f3_mg(k,l)
           ENDDO
           nzt_l = nzt_l / 2
-          dx_l  = dx_l * 2.0
-          dy_l  = dy_l * 2.0
+          dx_l  = dx_l * 2.0_wp
+          dy_l  = dy_l * 2.0_wp
        ENDDO
 
@@ -370,10 +370,10 @@
 !
 !-- Compute the reciprocal values of the horizontal grid lengths.
-    ddx = 1.0 / dx
-    ddy = 1.0 / dy
+    ddx = 1.0_wp / dx
+    ddy = 1.0_wp / dy
     dx2 = dx * dx
     dy2 = dy * dy
-    ddx2 = 1.0 / dx2
-    ddy2 = 1.0 / dy2
+    ddx2 = 1.0_wp / dx2
+    ddy2 = 1.0_wp / dy2
 
 !
@@ -433,6 +433,6 @@
     nzb_w_inner = nzb;  nzb_w_outer = nzb
 
-    rflags_s_inner = 1.0
-    rflags_invers  = 1.0
+    rflags_s_inner = 1.0_wp
+    rflags_invers  = 1.0_wp
 
 !
@@ -453,8 +453,8 @@
     nzb_diff_u = nzb_diff;  nzb_diff_v = nzb_diff
 
-    wall_e_x = 0.0;  wall_e_y = 0.0;  wall_u = 0.0;  wall_v = 0.0
-    wall_w_x = 0.0;  wall_w_y = 0.0
-    fwxp = 1.0;  fwxm = 1.0;  fwyp = 1.0;  fwym = 1.0
-    fxp  = 1.0;  fxm  = 1.0;  fyp  = 1.0;  fym  = 1.0
+    wall_e_x = 0.0_wp;  wall_e_y = 0.0_wp;  wall_u = 0.0_wp;  wall_v = 0.0_wp
+    wall_w_x = 0.0_wp;  wall_w_y = 0.0_wp
+    fwxp = 1.0_wp;  fwxm = 1.0_wp;  fwyp = 1.0_wp;  fwym = 1.0_wp
+    fxp  = 1.0_wp;  fxm  = 1.0_wp;  fyp  = 1.0_wp;  fym  = 1.0_wp
 
 !
@@ -471,10 +471,10 @@
     DO  k = 1, nzt
        vertical_influence(k) = MIN ( INT( l_grid(k) / &
-                     ( wall_adjustment_factor * dzw(k) ) + 0.5 ), nzt - k )
+                     ( wall_adjustment_factor * dzw(k) ) + 0.5_wp ), nzt - k )
     ENDDO
 
     DO  k = 1, MAXVAL( nzb_s_inner )
-       IF ( l_grid(k) > 1.5 * dx * wall_adjustment_factor .OR.  &
-            l_grid(k) > 1.5 * dy * wall_adjustment_factor )  THEN
+       IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.  &
+            l_grid(k) > 1.5_wp * dy * wall_adjustment_factor )  THEN
           WRITE( message_string, * ) 'grid anisotropy exceeds ', &
                                      'threshold given by only local', &
@@ -582,26 +582,26 @@
 !--       Street canyon size has to meet some requirements
           IF ( canyon_width_x /= 9999999.9_wp )  THEN
-             IF ( ( cxl < 1 ) .OR. ( cxr > nx-1 ) .OR. ( cwx < 3 ) .OR.  &
+             IF ( ( cxl < 1 ) .OR. ( cxr > nx-1 ) .OR. ( cwx < 3 ) .OR.        &
                ( ch < 3 ) )  THEN
-                WRITE( message_string, * ) 'inconsistent canyon parameters:', &
-                                           '&cxl=', cxl, 'cxr=', cxr,         &
-                                           'cwx=', cwx,                       &
+                WRITE( message_string, * ) 'inconsistent canyon parameters:',  &
+                                           '&cxl=', cxl, 'cxr=', cxr,          &
+                                           'cwx=', cwx,                        &
                                            'ch=', ch, 'nx=', nx, 'ny=', ny
                 CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 ) 
              ENDIF
           ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
-             IF ( ( cys < 1 ) .OR. ( cyn > ny-1 ) .OR. ( cwy < 3 ) .OR.  &
+             IF ( ( cys < 1 ) .OR. ( cyn > ny-1 ) .OR. ( cwy < 3 ) .OR.        &
                ( ch < 3 ) )  THEN
-                WRITE( message_string, * ) 'inconsistent canyon parameters:', &
-                                           '&cys=', cys, 'cyn=', cyn,         &
-                                           'cwy=', cwy,                       &
+                WRITE( message_string, * ) 'inconsistent canyon parameters:',  &
+                                           '&cys=', cys, 'cyn=', cyn,          &
+                                           'cwy=', cwy,                        &
                                            'ch=', ch, 'nx=', nx, 'ny=', ny
                 CALL message( 'init_grid', 'PA0206', 1, 2, 0, 6, 0 ) 
              ENDIF
           ENDIF
-          IF ( canyon_width_x /= 9999999.9_wp .AND. canyon_width_y /= 9999999.9_wp ) &
-               THEN
-             message_string = 'inconsistent canyon parameters:' //     &  
-                              '&street canyon can only be oriented' // &
+          IF ( canyon_width_x /= 9999999.9_wp .AND.                            &                 
+               canyon_width_y /= 9999999.9_wp )  THEN
+             message_string = 'inconsistent canyon parameters:' //             &   
+                              '&street canyon can only be oriented' //         &
                               '&either in x- or in y-direction'
              CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 )
@@ -679,5 +679,5 @@
 !-- scheme have to be reduced up to nzt (required at the lateral boundaries).
     nzb_max = MAXVAL( nzb_local )
-    IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR.    &
+    IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR.             &
          inflow_n .OR. outflow_n .OR. inflow_s .OR. outflow_s )  THEN
          nzb_max = nzt
@@ -693,7 +693,7 @@
 !--    Consistency checks
        IF ( MINVAL( nzb_local ) < 0  .OR.  MAXVAL( nzb_local ) > nz + 1 )  THEN
-          WRITE( message_string, * ) 'nzb_local values are outside the',      &
-                                'model domain',                               &
-                                '&MINVAL( nzb_local ) = ', MINVAL(nzb_local), &
+          WRITE( message_string, * ) 'nzb_local values are outside the',       &
+                                'model domain',                                &
+                                '&MINVAL( nzb_local ) = ', MINVAL(nzb_local),  &
                                 '&MAXVAL( nzb_local ) = ', MAXVAL(nzb_local)
           CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 )
@@ -701,7 +701,7 @@
 
        IF ( bc_lr == 'cyclic' )  THEN
-          IF ( ANY( nzb_local(:,-1) /= nzb_local(:,nx)   )  .OR. &
+          IF ( ANY( nzb_local(:,-1) /= nzb_local(:,nx)   )  .OR.               &
                ANY( nzb_local(:,0)  /= nzb_local(:,nx+1) ) )  THEN
-             message_string = 'nzb_local does not fulfill cyclic' // &
+             message_string = 'nzb_local does not fulfill cyclic' //           &
                               ' boundary condition in x-direction'
              CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 )
@@ -709,7 +709,7 @@
        ENDIF
        IF ( bc_ns == 'cyclic' )  THEN
-          IF ( ANY( nzb_local(-1,:) /= nzb_local(ny,:)   )  .OR. &
+          IF ( ANY( nzb_local(-1,:) /= nzb_local(ny,:)   )  .OR.               &
                ANY( nzb_local(0,:)  /= nzb_local(ny+1,:) ) )  THEN
-             message_string = 'nzb_local does not fulfill cyclic' // &
+             message_string = 'nzb_local does not fulfill cyclic' //           &
                               ' boundary condition in y-direction'
              CALL message( 'init_grid', 'PA0212', 1, 2, 0, 6, 0 )
@@ -770,5 +770,5 @@
        DO  j = -1, ny + 1
           DO  i = 0, nx
-             nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i), &
+             nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i),             &
                                  nzb_local(j,i+1) )
           ENDDO
@@ -776,5 +776,5 @@
        DO  i = nxl, nxr
           DO  j = nys, nyn
-             nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &
+             nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i),             &
                                      nzb_tmp(j+1,i) )
           ENDDO
@@ -812,5 +812,5 @@
        DO  i = nxl, nxr
           DO  j = nys, nyn
-             nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &
+             nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i),             &
                                      nzb_tmp(j+1,i) )
           ENDDO
@@ -844,5 +844,5 @@
        DO  j = nys, nyn
           DO  i = nxl, nxr
-             nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i), &
+             nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i),             &
                                      nzb_tmp(j,i+1) )
           ENDDO
@@ -891,6 +891,6 @@
           DO  j = nysg, nyng
              DO  k = nzb, nzt+1
-                IF ( k <= nzb_s_inner(j,i) )  rflags_s_inner(k,j,i) = 0.0
-                IF ( k <= nzb_s_inner(j,i) )  rflags_invers(j,i,k)  = 0.0
+                IF ( k <= nzb_s_inner(j,i) )  rflags_s_inner(k,j,i) = 0.0_wp
+                IF ( k <= nzb_s_inner(j,i) )  rflags_invers(j,i,k)  = 0.0_wp
              ENDDO
           ENDDO
@@ -938,22 +938,22 @@
 !--       u-component
           IF ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) )  THEN
-             wall_u(j,i) = 1.0   ! north wall (location of adjacent fluid)
-             fym(j,i)    = 0.0
-             fyp(j,i)    = 1.0
+             wall_u(j,i) = 1.0_wp   ! north wall (location of adjacent fluid)
+             fym(j,i)    = 0.0_wp
+             fyp(j,i)    = 1.0_wp
           ELSEIF ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) )  THEN
-             wall_u(j,i) = 1.0   ! south wall (location of adjacent fluid)
-             fym(j,i)    = 1.0
-             fyp(j,i)    = 0.0
+             wall_u(j,i) = 1.0_wp   ! south wall (location of adjacent fluid)
+             fym(j,i)    = 1.0_wp
+             fyp(j,i)    = 0.0_wp
           ENDIF
 !
 !--       v-component
           IF ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) )  THEN
-             wall_v(j,i) = 1.0   ! rigth wall (location of adjacent fluid)
-             fxm(j,i)    = 0.0
-             fxp(j,i)    = 1.0
+             wall_v(j,i) = 1.0_wp   ! rigth wall (location of adjacent fluid)
+             fxm(j,i)    = 0.0_wp
+             fxp(j,i)    = 1.0_wp
           ELSEIF ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) )  THEN
-             wall_v(j,i) = 1.0   ! left wall (location of adjacent fluid)
-             fxm(j,i)    = 1.0
-             fxp(j,i)    = 0.0
+             wall_v(j,i) = 1.0_wp   ! left wall (location of adjacent fluid)
+             fxm(j,i)    = 1.0_wp
+             fxp(j,i)    = 0.0_wp
           ENDIF
 !
@@ -961,24 +961,24 @@
 !--       production of tke
           IF ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) )  THEN
-             wall_e_y(j,i) =  1.0   ! north wall (location of adjacent fluid)
-             wall_w_y(j,i) =  1.0
-             fwym(j,i)     =  0.0
-             fwyp(j,i)     =  1.0
+             wall_e_y(j,i) =  1.0_wp   ! north wall (location of adjacent fluid)
+             wall_w_y(j,i) =  1.0_wp
+             fwym(j,i)     =  0.0_wp
+             fwyp(j,i)     =  1.0_wp
           ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) )  THEN
-             wall_e_y(j,i) = -1.0   ! south wall (location of adjacent fluid)
-             wall_w_y(j,i) =  1.0
-             fwym(j,i)     =  1.0
-             fwyp(j,i)     =  0.0
+             wall_e_y(j,i) = -1.0_wp   ! south wall (location of adjacent fluid)
+             wall_w_y(j,i) =  1.0_wp
+             fwym(j,i)     =  1.0_wp
+             fwyp(j,i)     =  0.0_wp
           ENDIF
           IF ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) )  THEN
-             wall_e_x(j,i) =  1.0   ! right wall (location of adjacent fluid)
-             wall_w_x(j,i) =  1.0
-             fwxm(j,i)     =  0.0
-             fwxp(j,i)     =  1.0
+             wall_e_x(j,i) =  1.0_wp   ! right wall (location of adjacent fluid)
+             wall_w_x(j,i) =  1.0_wp
+             fwxm(j,i)     =  0.0_wp
+             fwxp(j,i)     =  1.0_wp
           ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) )  THEN
-             wall_e_x(j,i) = -1.0   ! left wall (location of adjacent fluid)
-             wall_w_x(j,i) =  1.0
-             fwxm(j,i)     =  1.0
-             fwxp(j,i)     =  0.0
+             wall_e_x(j,i) = -1.0_wp   ! left wall (location of adjacent fluid)
+             wall_w_x(j,i) =  1.0_wp
+             fwxm(j,i)     =  1.0_wp
+             fwxp(j,i)     =  0.0_wp
           ENDIF
 !
@@ -1399,11 +1399,11 @@
 !--             North wall (y distance)
                 DO  k = wall_n(j,i), nzb_si
-                   l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), 0.5 * dy )
+                   l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), 0.5_wp * dy )
                 ENDDO
 !
 !--             Above North wall (yz distance)
                 DO  k = nzb_si + 1, nzb_si + vi
-                   l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i),     &
-                                          SQRT( 0.25 * dy**2 + &
+                   l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i),                     &
+                                          SQRT( 0.25_wp * dy**2 +              &
                                           ( zu(k) - zw(nzb_si) )**2 ) )
                 ENDDO
@@ -1413,12 +1413,12 @@
                    DO  k = corner_nl(j,i), nzb_si
                       l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1), &
-                                               0.5 * SQRT( dx**2 + dy**2 ) )
+                                               0.5_wp * SQRT( dx**2 + dy**2 ) )
                    ENDDO
 !
 !--                Above Northleft corner (xyz distance)
                    DO  k = nzb_si + 1, nzb_si + vi
-                      l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1),             &
-                                               SQRT( 0.25 * (dx**2 + dy**2) + &
-                                               ( zu(k) - zw(nzb_si) )**2 ) )
+                      l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1),              &
+                                            SQRT( 0.25_wp * (dx**2 + dy**2) +  &
+                                            ( zu(k) - zw(nzb_si) )**2 ) )
                    ENDDO
                 ENDIF
@@ -1427,13 +1427,13 @@
                 IF ( corner_nr(j,i) > 0 )  THEN
                    DO  k = corner_nr(j,i), nzb_si
-                       l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &
-                                                0.5 * SQRT( dx**2 + dy**2 ) )
+                       l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1),             &
+                                                0.5_wp * SQRT( dx**2 + dy**2 ) )
                    ENDDO
 !
 !--                Above northright corner (xyz distance)
                    DO  k = nzb_si + 1, nzb_si + vi
-                      l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &
-                                               SQRT( 0.25 * (dx**2 + dy**2) + &
-                                               ( zu(k) - zw(nzb_si) )**2 ) )
+                      l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1),              &
+                                            SQRT( 0.25_wp * (dx**2 + dy**2) +  &
+                                            ( zu(k) - zw(nzb_si) )**2 ) )
                    ENDDO
                 ENDIF
@@ -1444,12 +1444,11 @@
 !--             South wall (y distance)
                 DO  k = wall_s(j,i), nzb_si
-                   l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i), 0.5 * dy )
+                   l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i), 0.5_wp * dy )
                 ENDDO
 !
 !--             Above south wall (yz distance)
-                DO  k = nzb_si + 1, &
-                        nzb_si + vi
-                   l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i),     &
-                                          SQRT( 0.25 * dy**2 + &
+                DO  k = nzb_si + 1, nzb_si + vi
+                   l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i),                     &
+                                          SQRT( 0.25_wp * dy**2 +              &
                                           ( zu(k) - zw(nzb_si) )**2 ) )
                 ENDDO
@@ -1458,13 +1457,13 @@
                 IF ( corner_sl(j,i) > 0 )  THEN
                    DO  k = corner_sl(j,i), nzb_si
-                      l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1), &
-                                               0.5 * SQRT( dx**2 + dy**2 ) )
+                      l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1),              &
+                                               0.5_wp * SQRT( dx**2 + dy**2 ) )
                    ENDDO
 !
 !--                Above southleft corner (xyz distance)
                    DO  k = nzb_si + 1, nzb_si + vi
-                      l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1),             &
-                                               SQRT( 0.25 * (dx**2 + dy**2) + &
-                                               ( zu(k) - zw(nzb_si) )**2 ) )
+                      l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1),              &
+                                            SQRT( 0.25_wp * (dx**2 + dy**2) +  &
+                                            ( zu(k) - zw(nzb_si) )**2 ) )
                    ENDDO
                 ENDIF
@@ -1473,13 +1472,13 @@
                 IF ( corner_sr(j,i) > 0 )  THEN
                    DO  k = corner_sr(j,i), nzb_si
-                      l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1), &
-                                               0.5 * SQRT( dx**2 + dy**2 ) )
+                      l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1),              &
+                                               0.5_wp * SQRT( dx**2 + dy**2 ) )
                    ENDDO
 !
 !--                Above southright corner (xyz distance)
                    DO  k = nzb_si + 1, nzb_si + vi
-                      l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1),             &
-                                               SQRT( 0.25 * (dx**2 + dy**2) + &
-                                               ( zu(k) - zw(nzb_si) )**2 ) )
+                      l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1),              &
+                                            SQRT( 0.25_wp * (dx**2 + dy**2) +  &
+                                            ( zu(k) - zw(nzb_si) )**2 ) )
                    ENDDO
                 ENDIF
@@ -1491,28 +1490,28 @@
 !--             Left wall (x distance)
                 DO  k = wall_l(j,i), nzb_si
-                   l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1), 0.5 * dx )
+                   l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1), 0.5_wp * dx )
                 ENDDO
 !
 !--             Above left wall (xz distance)
                 DO  k = nzb_si + 1, nzb_si + vi
-                   l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1),     &
-                                          SQRT( 0.25 * dx**2 + &
+                   l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1),                     &
+                                       SQRT( 0.25_wp * dx**2 +                 &
+                                       ( zu(k) - zw(nzb_si) )**2 ) )
+                ENDDO
+             ENDIF
+
+             IF ( wall_r(j,i) > 0 )  THEN
+!
+!--             Right wall (x distance)
+                DO  k = wall_r(j,i), nzb_si
+                   l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), 0.5_wp * dx )
+                ENDDO
+!
+!--             Above right wall (xz distance)
+                DO  k = nzb_si + 1, nzb_si + vi
+                   l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1),                     &
+                                          SQRT( 0.25_wp * dx**2 +              &
                                           ( zu(k) - zw(nzb_si) )**2 ) )
                 ENDDO
-             ENDIF
-
-             IF ( wall_r(j,i) > 0 )  THEN
-!
-!--             Right wall (x distance)
-                DO  k = wall_r(j,i), nzb_si
-                   l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), 0.5 * dx )
-                ENDDO
-!
-!--             Above right wall (xz distance)
-                DO  k = nzb_si + 1, nzb_si + vi
-                   l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1),     &
-                                          SQRT( 0.25 * dx**2 + &
-                                          ( zu(k) - zw(nzb_si) )**2 ) )
-                ENDDO
 
              ENDIF
