Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 52)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 53)
@@ -2,5 +2,5 @@
 ---
 
-Wall functions for vertical walls now include diabatic conditions. New subroutine wall_fluxes. New 4D-array rif_wall.
+Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.
 
 new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files
Index: palm/trunk/SOURCE/diffusion_u.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_u.f90	(revision 52)
+++ palm/trunk/SOURCE/diffusion_u.f90	(revision 53)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! wall functions now include diabatic conditions, call of routine wall_fluxes
+! Wall functions now include diabatic conditions, call of routine wall_fluxes
 ! 
 ! Former revisions:
Index: palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_v.f90	(revision 52)
+++ palm/trunk/SOURCE/diffusion_v.f90	(revision 53)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! wall functions now include diabatic conditions, call of routine wall_fluxes
+! Wall functions now include diabatic conditions, call of routine wall_fluxes
 !
 ! Former revisions:
Index: palm/trunk/SOURCE/diffusion_w.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_w.f90	(revision 52)
+++ palm/trunk/SOURCE/diffusion_w.f90	(revision 53)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! wall functions now include diabatic conditions, call of routine wall_fluxes
+! Wall functions now include diabatic conditions, call of routine wall_fluxes
 !
 ! Former revisions:
Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 52)
+++ palm/trunk/SOURCE/production_e.f90	(revision 53)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Wall functions now include diabatic conditions, call of routine wall_fluxes
 !
 ! Former revisions:
@@ -67,5 +67,7 @@
 
        REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
-                   k1, k2, theta, temp, usvs, vsus, wsus, wsvs
+                   k1, k2, theta, temp
+
+       REAL, DIMENSION(nzb:nzt+1) ::  usvs, vsus, wsus, wsvs
 
 
@@ -119,48 +121,39 @@
                    k = nzb_diff_s_inner(j,i) - 1
                    dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
+                   dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                  u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+                   dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
+                   dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                  v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+                   dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
+
                    IF ( wall_e_y(j,i) /= 0.0 )  THEN
-                      usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) &
-                             / LOG( 0.5 * dy / z0(j,i) )
-                      usvs = usvs * ABS( usvs )
-                      dudy = wall_e_y(j,i) * usvs / km(k,j,i)
+                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
+                                          usvs, 1.0, 0.0, 0.0, 0.0 )
+
+                      dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
+                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
+                                          wsvs, 0.0, 0.0, 1.0, 0.0 )
+                      dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
                    ELSE
                       dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                      dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                      w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
                    ENDIF
-                   dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                  u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
 
                    IF ( wall_e_x(j,i) /= 0.0 )  THEN
-                      vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) &
-                             / LOG( 0.5 * dx / z0(j,i))
-                      vsus = vsus * ABS( vsus )
-                      dvdx = wall_e_x(j,i) * vsus / km(k,j,i)
+                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
+                                          vsus, 0.0, 1.0, 0.0, 0.0 )
+                      dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
+                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
+                                          wsus, 0.0, 0.0, 0.0, 1.0 )
+                      dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
                    ELSE
                       dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                   ENDIF
-                   dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
-                   dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                  v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
-
-                   IF ( wall_e_x(j,i) /= 0.0 )  THEN
-                      wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &
-                             / LOG( 0.5 * dx / z0(j,i))
-                      wsus = wsus * ABS( wsus )
-                      dwdx = wall_e_x(j,i) * wsus / km(k,j,i)
-                   ELSE
                       dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
                    ENDIF
-                   IF ( wall_e_y(j,i) /= 0.0 ) THEN
-                      wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &
-                             / LOG( 0.5 * dy / z0(j,i))
-                      wsvs = wsvs * ABS( wsvs )
-                      dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)
-                   ELSE
-                      dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                      w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                   ENDIF
-                   dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
                    def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
@@ -172,64 +165,40 @@
                    tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
 
-                ENDIF
-
-             ENDDO
-
-!
-!--          3 - Wird nur ausgefuehrt, wenn mindestens ein Niveau
-!--          zwischen 2 und 4 liegt, d.h. ab einer Gebaeudemindest-
-!--          hoehe von 2 dz. 'Nur Wand: Wall functions'
-             DO  j = nys, nyn
-
-                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
-                THEN
+
+!
+!--                3 - Wird nur ausgefuehrt, wenn mindestens ein Niveau
+!--                zwischen 2 und 4 liegt, d.h. ab einer Gebaeudemindest-
+!--                hoehe von 2 dz. Die Wall fluxes sind in diesem Fall
+!--                schon unter (2) berechnet worden. 'Nur Wand: Wall functions'
 
                    DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
 
                       dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
+                      dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+                      dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                      dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+                      dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
+
                       IF ( wall_e_y(j,i) /= 0.0 )  THEN
-                         usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) &
-                                / LOG( 0.5 * dy / z0(j,i) )
-                         usvs = usvs * ABS( usvs )
-                         dudy = wall_e_y(j,i) * usvs / km(k,j,i)
+                         dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
+                         dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
                       ELSE
                          dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
                                          u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                         dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
                       ENDIF
-                      dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
 
                       IF ( wall_e_x(j,i) /= 0.0 )  THEN
-                         vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) &
-                                / LOG( 0.5 * dx / z0(j,i))
-                         vsus = vsus * ABS( vsus )
-                         dvdx = wall_e_x(j,i) * vsus / km(k,j,i)
+                         dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
+                         dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
                       ELSE
                          dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
                                          v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                      ENDIF
-                      dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                      dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                      IF ( wall_e_x(j,i) /= 0.0 )  THEN
-                         wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &
-                                / LOG( 0.5 * dx / z0(j,i))
-                         wsus = wsus * ABS( wsus )
-                         dwdx = wall_e_x(j,i) * wsus / km(k,j,i)
-                      ELSE
                          dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
                                          w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
                       ENDIF
-                      IF ( wall_e_y(j,i) /= 0.0 )  THEN
-                         wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &
-                                / LOG( 0.5 * dy / z0(j,i))
-                         wsvs = wsvs * ABS( wsvs )
-                         dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)
-                      ELSE
-                         dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                      ENDIF
-                      dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
                       def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
@@ -423,5 +392,5 @@
                 DO  j = nys, nyn
 
-                   k = nzb_diff_s_inner(j,i)
+                   k = nzb_diff_s_inner(j,i)-1
 
                    IF ( .NOT. cloud_physics )  THEN
@@ -505,5 +474,7 @@
 
        REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
-                   k1, k2, theta, temp, usvs, vsus, wsus,wsvs
+                   k1, k2, theta, temp
+
+       REAL, DIMENSION(nzb:nzt+1) ::  usvs, vsus, wsus,wsvs
 
 !
@@ -549,48 +520,38 @@
 
              dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
+             dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                            u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+             dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+             dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                            v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+             dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
+
              IF ( wall_e_y(j,i) /= 0.0 )  THEN
-                usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) &
-                       / LOG( 0.5 * dy / z0(j,i) )
-                usvs = usvs * ABS( usvs )
-                dudy = wall_e_y(j,i) * usvs / km(k,j,i)
+                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
+                                    usvs, 1.0, 0.0, 0.0, 0.0 )
+                dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
+                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
+                                    wsvs, 0.0, 0.0, 1.0, 0.0 )
+                dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
              ELSE
                 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
              ENDIF
-             dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                            u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
 
              IF ( wall_e_x(j,i) /= 0.0 )  THEN
-                vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) &
-                       / LOG( 0.5 * dx / z0(j,i))
-                vsus = vsus * ABS( vsus )
-                dvdx = wall_e_x(j,i) * vsus / km(k,j,i)
+                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
+                                    vsus, 0.0, 1.0, 0.0, 0.0 )
+                dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
+                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
+                                    wsus, 0.0, 0.0, 0.0, 1.0 )
+                dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
              ELSE
                 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-             ENDIF
-             dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-             dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                            v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
-
-             IF ( wall_e_x(j,i) /= 0.0 )  THEN
-                wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &
-                       / LOG( 0.5 * dx / z0(j,i))
-                wsus = wsus * ABS( wsus )
-                dwdx = wall_e_x(j,i) * wsus / km(k,j,i)
-             ELSE
                 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
              ENDIF
-             IF ( wall_e_y(j,i) /= 0.0 )  THEN
-                wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &
-                       / LOG( 0.5 * dy / z0(j,i))
-                wsvs = wsvs * ABS( wsvs )
-                dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)
-             ELSE
-                dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-             ENDIF
-             dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
              def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
@@ -605,52 +566,35 @@
 !--          3 - Wird nur ausgefuehrt, wenn mindestens ein Niveau
 !--          zwischen 2 und 4 liegt, d.h. ab einer Gebaeudemindest-
-!--          hoehe von 2 dz. 'Nur Wand: Wall functions'
+!--          hoehe von 2 dz. Di Wall fluxes sind in diesem Fall bereits vorab
+!--          unter (2) berechnet worden. 'Nur Wand: Wall functions'
              DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
 
                 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
+                dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                               u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+                dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                               v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+                dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
+
                 IF ( wall_e_y(j,i) /= 0.0 )  THEN
-                   usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) &
-                          / LOG( 0.5 * dy / z0(j,i) )
-                   usvs = usvs * ABS( usvs )
-                   dudy = wall_e_y(j,i) * usvs / km(k,j,i)
+                   dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
+                   dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
                 ELSE
                    dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                ENDIF
-                dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                               u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+                   dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                ENDIF
 
                 IF ( wall_e_x(j,i) /= 0.0 )  THEN
-                   vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) &
-                          / LOG( 0.5 * dx / z0(j,i))
-                   vsus = vsus * ABS( vsus )
-                   dvdx = wall_e_x(j,i) * vsus / km(k,j,i)
+                   dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
+                   dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
                 ELSE
                    dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                ENDIF
-                dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                               v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                IF ( wall_e_x(j,i) /= 0.0 )  THEN
-                   wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &
-                          / LOG( 0.5 * dx / z0(j,i))
-                   wsus = wsus * ABS( wsus )
-                   dwdx = wall_e_x(j,i) * wsus / km(k,j,i)
-                ELSE
                    dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
                 ENDIF
-                IF ( wall_e_y(j,i) /= 0.0 )  THEN
-                   wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &
-                          / LOG( 0.5 * dy / z0(j,i))
-                   wsvs = wsvs * ABS( wsvs )
-                   dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)
-                ELSE
-                   dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                ENDIF
-                dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
                 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
Index: palm/trunk/SOURCE/wall_fluxes.f90
===================================================================
--- palm/trunk/SOURCE/wall_fluxes.f90	(revision 52)
+++ palm/trunk/SOURCE/wall_fluxes.f90	(revision 53)
@@ -26,6 +26,6 @@
     IMPLICIT NONE
 
-    INTEGER ::  i, j, k, nzb_w, nzt_w
-    REAL    ::  a, b, c1, c2, h1, h2, delta_p
+    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
@@ -34,6 +34,7 @@
 
 
-    delta_p   = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
-    wall_flux = 0.0
+    zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
+    wall_flux  = 0.0
+    wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
 
 !
@@ -41,35 +42,45 @@
 !-- 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,NINT(a+2*b+3*c1+4*c2))
-       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) +  &
-            a * 0.25 * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) + &
-            b * 0.25 * ( 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
+       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( delta_p / z0(j,i) ) +                    &
-                    5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p  &
+          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 / delta_p * z0(j,i) ) )
+          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
+
 !
 !--       If a borderline case occurs, the formula for stable stratification
@@ -77,20 +88,21 @@
 !--       argument of the logarithm.
           IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
-             us_wall = kappa * vel_total / (                           &
-                       LOG( delta_p / z0(j,i) ) +                         &
-                       5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p &
-                                                   )
+             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
+                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
+                                           )
           ELSE
-             us_wall = kappa * vel_total / (                           &
+             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 delta_p/L (corresponds to neutral Richardson flux number 
+
+       ENDIF
+
+!
+!--    (4) Compute zp/L (corresponds to neutral Richardson flux number 
 !--        rifs)
-       rifs = -1.0 * delta_p * kappa * g * wspts / &
-                       ( pt_i * ( us_wall**3 + 1E-30 ) )
+       rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * ( us_wall**3 + 1E-30 ) )
+
 !
 !--    Limit the value range of the Richardson numbers.
@@ -101,19 +113,23 @@
        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( delta_p / z0(j,i) ) +                                &
-               5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p        &
-                                                                         )
-       ELSE
+          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 / delta_p * z0(j,i) ) )
+          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
+
 !
 !--       If a borderline case occurs, the formula for stable stratification
@@ -121,16 +137,17 @@
 !--       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( delta_p / z0(j,i) ) +                                &
-                  5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p        &
-                                                                            )
+             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 ) )                            &
-                                                                            )
+             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
        wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) )
@@ -138,6 +155,115 @@
 !
 !--    store rifs for next time step
-       rif_wall(k,j,i,NINT(a+2*b+3*c1+4*c2)) = rifs
+       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)
+!
+! Description:
+! ------------
+! Calculates momentum fluxes at vertical walls for routine production_e
+! assuming Monin-Obukhov 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, kk, nzb_w, nzt_w, wall_index
+    REAL    ::  a, b, c1, c2, h1, h2, 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')
+       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
+             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
+       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
