Index: palm/trunk/SOURCE/wall_fluxes.f90
===================================================================
--- palm/trunk/SOURCE/wall_fluxes.f90	(revision 55)
+++ palm/trunk/SOURCE/wall_fluxes.f90	(revision 56)
@@ -1,3 +1,3 @@
- SUBROUTINE wall_fluxes( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
+ MODULE wall_fluxes_mod
 !------------------------------------------------------------------------------!
 ! Actual revisions:
@@ -15,162 +15,346 @@
 ! similarity.
 ! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
-!------------------------------------------------------------------------------!
-
-    USE arrays_3d
-    USE control_parameters
-    USE grid_variables
-    USE indices
-    USE statistics
-    USE user
-
-    IMPLICIT NONE
-
-    INTEGER ::  i, j, k, nzb_w, nzt_w, wall_index
-    REAL    ::  a, b, c1, c2, h1, h2, zp
-
-    REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
-
-    REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
-
-
-    zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
-    wall_flux  = 0.0
-    wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
-
-!
-!-- All subsequent variables are computed for the respective location where 
-!-- the relevant variable is defined
-    DO  k = nzb_w, nzt_w
-
-!
-!--    (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
-       rifs  = rif_wall(k,j,i,wall_index)
-
-       u_i   = a * u(k,j,i) + c1 * 0.25 * &
-               ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
-
-       v_i   = b * v(k,j,i) + c2 * 0.25 * &
-               ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
-
-       ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                             &
-                   a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
-                 + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
-                                               )
-       pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + b * pt(k,j-1,i)  &
-                       + ( c1 + c2 ) * pt(k+1,j,i) )
-
-       pts   = pt_i - hom(k,1,4,0)
-       wspts = ws * pts
-
-!
-!--    (2) Compute wall-parallel absolute velocity vel_total
-       vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 )
-
-!
-!--    (3) Compute wall friction velocity us_wall
-       IF ( rifs >= 0.0 )  THEN
-
-!
-!--       Stable stratification (and neutral)
-          us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +               &
-                                          5.0 * rifs * ( zp - z0(j,i) ) / zp  &
-                                        )
-       ELSE
-
-!
-!--       Unstable stratification
-          h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
-          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
-
-!
-!--       If a borderline case occurs, the formula for stable stratification
-!--       must be used anyway, or else a zero division would occur in the 
-!--       argument of the logarithm.
-          IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
+! The all-gridpoint version of wall_fluxes_e is not used so far, because
+! it gives slightly different results from the ij-version for some unknown
+! reason.
+!------------------------------------------------------------------------------!
+    PRIVATE
+    PUBLIC wall_fluxes, wall_fluxes_e
+    
+    INTERFACE wall_fluxes
+       MODULE PROCEDURE wall_fluxes
+       MODULE PROCEDURE wall_fluxes_ij
+    END INTERFACE wall_fluxes
+    
+    INTERFACE wall_fluxes_e
+       MODULE PROCEDURE wall_fluxes_e
+       MODULE PROCEDURE wall_fluxes_e_ij
+    END INTERFACE wall_fluxes_e
+ 
+ CONTAINS
+
+!------------------------------------------------------------------------------!
+! Call for all grid points
+!------------------------------------------------------------------------------!
+    SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, ixp, jyp, nzb_uvw_inner, &
+                            nzb_uvw_outer, wall )
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE statistics
+       USE user
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, ixp, j, jyp, k, wall_index
+
+       INTEGER, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1) ::  nzb_uvw_inner, &
+                                                       nzb_uvw_outer
+       REAL ::  a, b, c1, c2, h1, h2, zp
+       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
+
+       REAL, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1)           ::  wall
+       REAL, DIMENSION(nzb:nzt+1,nys:nyn+jyp,nxl:nxr+ixp) ::  wall_flux
+
+
+       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
+       wall_flux  = 0.0
+       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
+
+       DO  i = nxl, nxr+ixp
+          DO  j = nys, nyn+jyp
+
+             IF ( wall(j,i) /= 0.0 )  THEN
+!
+!--             All subsequent variables are computed for the respective
+!--             location where the relevant variable is defined
+                DO  k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
+
+!
+!--                (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
+                   rifs  = rif_wall(k,j,i,wall_index)
+
+                   u_i   = a * u(k,j,i) + c1 * 0.25 * &
+                           ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
+
+                   v_i   = b * v(k,j,i) + c2 * 0.25 * &
+                           ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
+
+                   ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                   &
+                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
+                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
+                                                           )
+                   pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + &
+                                   b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
+
+                   pts   = pt_i - hom(k,1,4,0)
+                   wspts = ws * pts
+
+!
+!--                (2) Compute wall-parallel absolute velocity vel_total
+                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
+
+!
+!--                (3) Compute wall friction velocity us_wall
+                   IF ( rifs >= 0.0 )  THEN
+
+!
+!--                   Stable stratification (and neutral)
+                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
+                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
+                                                    )
+                   ELSE
+
+!
+!--                   Unstable stratification
+                      h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ))
+
+!
+!--                   If a borderline case occurs, the formula for stable
+!--                   stratification must be used anyway, or else a zero
+!--                   division would occur in the argument of the logarithm.
+                      IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
+                         us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
+                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
+                                                       )
+                      ELSE
+                         us_wall = kappa * vel_total / (                       &
+                            LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + &
+                            2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                  &
+                                                       )
+                      ENDIF
+
+                   ENDIF
+
+!
+!--                (4) Compute zp/L (corresponds to neutral Richardson flux
+!--                    number rifs)
+                   rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
+                                                        ( us_wall**3 + 1E-30 ) )
+
+!
+!--                Limit the value range of the Richardson numbers.
+!--                This is necessary for very small velocities (u,w --> 0),
+!--                because the absolute value of rif can then become very
+!--                large, which in consequence would result in very large
+!--                shear stresses and very small momentum fluxes (both are
+!--                generally unrealistic).
+                   IF ( rifs < rif_min )  rifs = rif_min
+                   IF ( rifs > rif_max )  rifs = rif_max
+
+!
+!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
+                   IF ( rifs >= 0.0 )  THEN
+
+!
+!--                   Stable stratification (and neutral)
+                      wall_flux(k,j,i) = kappa *                               &
+                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
+                              (  LOG( zp / z0(j,i) ) +                         &
+                                 5.0 * rifs * ( zp - z0(j,i) ) / zp            &
+                              )
+                   ELSE
+
+!
+!--                   Unstable stratification
+                      h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ))
+
+!
+!--                   If a borderline case occurs, the formula for stable
+!--                   stratification must be used anyway, or else a zero
+!--                   division would occur in the argument of the logarithm. 
+                      IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
+                         wall_flux(k,j,i) = kappa *                            &
+                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
+                              ( LOG( zp / z0(j,i) ) +                          &
+                                5.0 * rifs * ( zp - z0(j,i) ) / zp             &
+                              )
+                      ELSE
+                         wall_flux(k,j,i) = kappa *                            &
+                             ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /  &
+                             ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) )&
+                               + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )             &
+                             )
+                      ENDIF
+
+                   ENDIF
+                   wall_flux(k,j,i) = -wall_flux(k,j,i) * ABS(wall_flux(k,j,i))
+
+!
+!--                store rifs for next time step
+                   rif_wall(k,j,i,wall_index) = rifs
+
+                ENDDO
+
+             ENDIF
+
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE wall_fluxes
+
+
+
+!------------------------------------------------------------------------------!
+! Call for all grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE statistics
+       USE user
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k, nzb_w, nzt_w, wall_index
+       REAL    ::  a, b, c1, c2, h1, h2, zp
+
+       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
+
+       REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
+
+
+       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
+       wall_flux  = 0.0
+       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
+
+!
+!--    All subsequent variables are computed for the respective location where 
+!--    the relevant variable is defined
+       DO  k = nzb_w, nzt_w
+
+!
+!--       (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
+          rifs  = rif_wall(k,j,i,wall_index)
+
+          u_i   = a * u(k,j,i) + c1 * 0.25 * &
+                  ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
+
+          v_i   = b * v(k,j,i) + c2 * 0.25 * &
+                  ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
+
+          ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                            &
+                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
+                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
+                                                  )
+          pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + b * pt(k,j-1,i)  &
+                          + ( c1 + c2 ) * pt(k+1,j,i) )
+
+          pts   = pt_i - hom(k,1,4,0)
+          wspts = ws * pts
+
+!
+!--       (2) Compute wall-parallel absolute velocity vel_total
+          vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 )
+
+!
+!--       (3) Compute wall friction velocity us_wall
+          IF ( rifs >= 0.0 )  THEN
+
+!
+!--          Stable stratification (and neutral)
              us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
                                             5.0 * rifs * ( zp - z0(j,i) ) / zp &
                                            )
           ELSE
-             us_wall = kappa * vel_total / (                                   &
+
+!
+!--          Unstable stratification
+             h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+             h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
+
+!
+!--          If a borderline case occurs, the formula for stable stratification
+!--          must be used anyway, or else a zero division would occur in the 
+!--          argument of the logarithm.
+             IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
+                us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +          &
+                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
+                                              )
+             ELSE
+                us_wall = kappa * vel_total / (                                &
                             LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + &
                             2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                  &
-                                           )
+                                              )
+             ENDIF
+
           ENDIF
 
-       ENDIF
-
-!
-!--    (4) Compute zp/L (corresponds to neutral Richardson flux number 
-!--        rifs)
-       rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * ( us_wall**3 + 1E-30 ) )
-
-!
-!--    Limit the value range of the Richardson numbers.
-!--    This is necessary for very small velocities (u,w --> 0), because
-!--    the absolute value of rif can then become very large, which in
-!--    consequence would result in very large shear stresses and very
-!--    small momentum fluxes (both are generally unrealistic).
-       IF ( rifs < rif_min )  rifs = rif_min
-       IF ( rifs > rif_max )  rifs = rif_max
-
-!
-!--    (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
-       IF ( rifs >= 0.0 )  THEN
-
-!
-!--       Stable stratification (and neutral)
-          wall_flux(k) = kappa *                                          &
-                         ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
-                         (  LOG( zp / z0(j,i) ) +                         &
-                            5.0 * rifs * ( zp - z0(j,i) ) / zp            &
-                         )
-       ELSE
-
-!
-!--       Unstable stratification
-          h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
-          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
-
-!
-!--       If a borderline case occurs, the formula for stable stratification
-!--       must be used anyway, or else a zero division would occur in the 
-!--       argument of the logarithm. 
-          IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
+!
+!--       (4) Compute zp/L (corresponds to neutral Richardson flux number 
+!--           rifs)
+          rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * (us_wall**3 + 1E-30) )
+
+!
+!--       Limit the value range of the Richardson numbers.
+!--       This is necessary for very small velocities (u,w --> 0), because
+!--       the absolute value of rif can then become very large, which in
+!--       consequence would result in very large shear stresses and very
+!--       small momentum fluxes (both are generally unrealistic).
+          IF ( rifs < rif_min )  rifs = rif_min
+          IF ( rifs > rif_max )  rifs = rif_max
+
+!
+!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
+          IF ( rifs >= 0.0 )  THEN
+
+!
+!--          Stable stratification (and neutral)
              wall_flux(k) = kappa *                                          &
                             ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
-                            ( LOG( zp / z0(j,i) ) +                          &
-                              5.0 * rifs * ( zp - z0(j,i) ) / zp             &
+                            (  LOG( zp / z0(j,i) ) +                         &
+                               5.0 * rifs * ( zp - z0(j,i) ) / zp            &
                             )
           ELSE
-             wall_flux(k) = kappa *                                            &
-                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /   &
-                            ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) &
-                              + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )              &
-                            )
+
+!
+!--          Unstable stratification
+             h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+             h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
+
+!
+!--          If a borderline case occurs, the formula for stable stratification
+!--          must be used anyway, or else a zero division would occur in the 
+!--          argument of the logarithm. 
+             IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
+                wall_flux(k) = kappa *                                         &
+                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
+                              ( LOG( zp / z0(j,i) ) +                          &
+                                5.0 * rifs * ( zp - z0(j,i) ) / zp             &
+                              )
+             ELSE
+                wall_flux(k) = kappa *                                         &
+                             ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /  &
+                             ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) )&
+                               + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )             &
+                             )
+             ENDIF
+
           ENDIF
-
-       ENDIF
-       wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) )
-
-!
-!--    store rifs for next time step
-       rif_wall(k,j,i,wall_index) = rifs
-
-    ENDDO
-
- END SUBROUTINE wall_fluxes
-
-
-
- SUBROUTINE wall_fluxes_e( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
-!------------------------------------------------------------------------------!
-! Actual revisions:
-! -----------------
-! 
-!
-! Former revisions:
-! -----------------
-! Initial version (2007/03/07)
-!
+          wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) )
+
+!
+!--       store rifs for next time step
+          rif_wall(k,j,i,wall_index) = rifs
+
+       ENDDO
+
+    END SUBROUTINE wall_fluxes_ij
+
+
+
+!------------------------------------------------------------------------------!
+! Call for all grid points
+!------------------------------------------------------------------------------!
+    SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
+
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
@@ -180,90 +364,198 @@
 !------------------------------------------------------------------------------!
 
-    USE arrays_3d
-    USE control_parameters
-    USE grid_variables
-    USE indices
-    USE statistics
-    USE user
-
-    IMPLICIT NONE
-
-    INTEGER ::  i, j, k, kk, nzb_w, nzt_w, wall_index
-    REAL    ::  a, b, c1, c2, h1, h2, vel_zp, zp
-
-    REAL ::  rifs
-
-    REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
-
-
-    zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
-    wall_flux  = 0.0
-    wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
-
-!
-!-- All subsequent variables are computed for the respective location where 
-!-- the relevant variable is defined
-    DO  k = nzb_w, nzt_w
-
-!
-!--    (1) Compute rifs
-       IF ( k == nzb_w )  THEN
-          kk = nzb_w
-       ELSE
-          kk = k-1
-       ENDIF
-       rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                         &
-                       a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
-                       c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
-                     )
-
-!
-!--    Skip (2) to (4) of wall_fluxes, because here rifs is already available
-!--    from (1)
-
-!
-!--    (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
-       vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
-                              b * ( v(k,j,i) + v(k,j+1,i) ) +  &
-                        (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )    &
-                      )
-
-       IF ( rifs >= 0.0 )  THEN
-
-!
-!--       Stable stratification (and neutral)
-          wall_flux(k) = kappa *  vel_zp / &
-                         ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
-       ELSE
-
-!
-!--       Unstable stratification
-          h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
-          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
-
-!
-!--       If a borderline case occurs, the formula for stable stratification
-!--       must be used anyway, or else a zero division would occur in the 
-!--       argument of the logarithm. 
-          IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
-             wall_flux(k) = kappa * vel_zp /                                 &
-                            ( LOG( zp / z0(j,i) ) +                          &
-                              5.0 * rifs * ( zp - z0(j,i) ) / zp             &
-                            )
-          ELSE
-             wall_flux(k) = kappa * vel_zp /                                   &
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE statistics
+       USE user
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k, kk, wall_index
+       REAL    ::  a, b, c1, c2, h1, h2, vel_zp, zp
+
+       REAL ::  rifs
+
+       REAL, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1)   ::  wall
+       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
+
+
+       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
+       wall_flux  = 0.0
+       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+
+             IF ( wall(j,i) /= 0.0 )  THEN
+!
+!--             All subsequent variables are computed for the respective
+!--             location where the relevant variable is defined
+                DO  k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
+
+!
+!--                (1) Compute rifs
+                   IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
+                      kk = nzb_diff_s_inner(j,i)-1
+                   ELSE
+                      kk = k-1
+                   ENDIF
+                   rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                &
+                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
+                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
+                                 )
+
+!
+!--                Skip (2) to (4) of wall_fluxes, because here rifs is
+!--                already available from (1)
+
+!
+!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
+                   vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
+                                          b * ( v(k,j,i) + v(k,j+1,i) ) +  &
+                                    (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )    &
+                                  )
+
+                   IF ( rifs >= 0.0 )  THEN
+
+!
+!--                   Stable stratification (and neutral)
+                      wall_flux(k,j,i) = kappa *  vel_zp / &
+                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
+                   ELSE
+
+!
+!--                   Unstable stratification
+                      h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ))
+
+!
+!--                   If a borderline case occurs, the formula for stable
+!--                   stratification must be used anyway, or else a zero
+!--                   division would occur in the argument of the logarithm. 
+                      IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
+                         wall_flux(k,j,i) = kappa * vel_zp /                 &
+                                        ( LOG( zp / z0(j,i) ) +              &
+                                          5.0 * rifs * ( zp - z0(j,i) ) / zp &
+                                        )
+                      ELSE
+                         wall_flux(k,j,i) = kappa * vel_zp /                   &
                             ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) &
                               + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )              &
                             )
+                      ENDIF
+
+                   ENDIF
+                   wall_flux(k,j,i) = wall_flux(k,j,i) * ABS( wall_flux(k,j,i) )
+
+!
+!--                Store rifs for next time step
+                   rif_wall(k,j,i,wall_index) = rifs
+
+                ENDDO
+
+             ENDIF
+
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE wall_fluxes_e
+
+
+
+!------------------------------------------------------------------------------!
+! Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE statistics
+       USE user
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k, kk, nzb_w, nzt_w, wall_index
+       REAL    ::  a, b, c1, c2, h1, h2, vel_zp, zp
+
+       REAL ::  rifs
+
+       REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
+
+
+       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
+       wall_flux  = 0.0
+       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
+
+!
+!--    All subsequent variables are computed for the respective location where 
+!--    the relevant variable is defined
+       DO  k = nzb_w, nzt_w
+
+!
+!--       (1) Compute rifs
+          IF ( k == nzb_w )  THEN
+             kk = nzb_w
+          ELSE
+             kk = k-1
           ENDIF
-
-       ENDIF
-       wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) )
-
-!
-!--    store rifs for next time step
-       rif_wall(k,j,i,wall_index) = rifs
-
-    ENDDO
-
- END SUBROUTINE wall_fluxes_e
+          rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                         &
+                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
+                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
+                        )
+
+!
+!--       Skip (2) to (4) of wall_fluxes, because here rifs is already available
+!--       from (1)
+
+!
+!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
+          vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
+                                 b * ( v(k,j,i) + v(k,j+1,i) ) +  &
+                           (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )    &
+                         )
+
+          IF ( rifs >= 0.0 )  THEN
+
+!
+!--          Stable stratification (and neutral)
+             wall_flux(k) = kappa *  vel_zp / &
+                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
+          ELSE
+
+!
+!--          Unstable stratification
+             h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+             h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
+
+!
+!--          If a borderline case occurs, the formula for stable stratification
+!--          must be used anyway, or else a zero division would occur in the 
+!--          argument of the logarithm. 
+             IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
+                wall_flux(k) = kappa * vel_zp /                     &
+                               ( LOG( zp / z0(j,i) ) +              &
+                                 5.0 * rifs * ( zp - z0(j,i) ) / zp &
+                               )
+             ELSE
+                wall_flux(k) = kappa * vel_zp /                                &
+                            ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) &
+                              + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )              &
+                            )
+             ENDIF
+
+          ENDIF
+          wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) )
+
+!
+!--       Store rifs for next time step
+          rif_wall(k,j,i,wall_index) = rifs
+
+       ENDDO
+
+    END SUBROUTINE wall_fluxes_e_ij
+
+ END MODULE wall_fluxes_mod
