Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 184)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 187)
@@ -42,4 +42,8 @@
 Changed:
 -------
+Modification of the integrated version of the profile function for momentum 
+for unstable stratification; more consistent flux definitions. (wall_fluxes, 
+production_e)
+
 Strict grid matching along z is not needed for mg-solver. (check_parameters)
 
@@ -79,9 +83,14 @@
 (init_1d_model)
 
-advec_s_ups, advec_u_ups, advec_v_ups, advec_w_ups, calc_spectra, check_open, check_parameters, cpu_statistics, init_1d_model, init_3d_model, modules, palm, parin, poisfft, read_var_list, read_3d_binary, transpose, write_var_list, write_3d_binary
+advec_s_ups, advec_u_ups, advec_v_ups, advec_w_ups, calc_spectra, check_open, check_parameters, cpu_statistics, init_1d_model, init_3d_model, modules, palm, parin, poisfft, production_e, read_var_list, read_3d_binary, transpose, wall_fluxes, write_var_list, write_3d_binary
 
 
 Errors:
 ------
+Bugfix: change definition of us_wall from 1D to 2D. Tests showed that this decreases 
+u* by some 10% and increases TKE and momentum fluxes by some 10% because friction 
+was underestimated and momentum fluxes were wrongly calculated due to the bug. 
+(prandtl_fluxes, wall_fluxes)
+
 Bugfix: calculation of horizontal fluxes at vertical walls (diffusion_s)
 
@@ -106,20 +115,27 @@
 
 Introduce prefix_chr to ensure unique dvrp_file path.
+
 small bugfixes for user_interface sample code (comments):
 - initialize ustvst with 0.0 as it is now computed only until nxr and nyn
 - two ALLOCATE statements moved from user_read_restart_data back to user_init
 - remove 'READ (13) u2_av' statement in user_read_restart_data
+
 Bugfix: remove IF statement in plant_canopy_model_ij (plant_canopy_model)
+
 Bugfix: divide sums(k,8) (e) and sums(k,34) (e*) by ngp_2dh_s_inner(k,sr) 
 (like other scalars) (flow_statistics)
+
 Bugfix: dopr_time_count was written on the binary file, which caused that
 NetCDF files newly created by restart files (no append of existing files!)
 contained uneccessary time levels. (read_3d_binary, write_3d_binary)
+
 Bugfix: extra '*' removed in user_statistics sample code (user_interface)
+
 Bugfix: a stop command was missing in some cases of the parallel branch (local_stop)
+
 Bugfix in volume flow control for non-cyclic boundary conditions (pres)
+
 Bugfix: misplaced #endif directives (combine_plot_fields)
 
+check_parameters, diffusion_s, flow_statistics, init_dvrp, init_3d_model, local_stop, plant_canopy_model, poismg, prandtl_fluxes, pres, read_3d_binary, user_interface, wall_fluxes, write_3d_binary
 
-check_parameters, diffusion_s, flow_statistics, init_dvrp, init_3d_model, local_stop, plant_canopy_model, poismg, pres, read_3d_binary, user_interface, write_3d_binary
-
Index: palm/trunk/SOURCE/prandtl_fluxes.f90
===================================================================
--- palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 184)
+++ palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 187)
@@ -4,5 +4,20 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Bugfix: change definition of us_wall from 1D to 2D:
+!         Modification of the evaluation of the vertical turbulent momentum 
+!         fluxes u'w' and v'w'; the first usws that is computed corresponds 
+!         to -u'w'/u* and not as priorily assumed to (-u'w')**0.5, the first 
+!         vsws that is computed corresponds to -v'w'/u* and not as priorily
+!         assumed to (-v'w')**0.5. Therefore, the intermediate result for 
+!         usws has to be multiplied by -u* instead by itself in order to 
+!         get u'w'. Accordingly, the intermediate result for vsws has to be 
+!         multiplied by -u* instead by itself in order to get v'w'. As u* 
+!         is calculated for the position of a scalar an additional 
+!         interpolation of u* to the position of u and v, respectively, 
+!         is necessary. As u* is not determined for the ghost points on 
+!         each PE, an additional exchange of information from neighbouring 
+!         PEs is necessary.!  
+! Change: Modification of the integrated version of the profile function for 
+!         momentum for unstable stratification
 !
 ! Former revisions:
@@ -81,19 +96,9 @@
 !--             Unstable stratification
                 a = SQRT( 1.0 - 16.0 * rif(j,i) )
-                b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * 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 ( a == 1.0  .OR.  b == 1.0 )  THEN
-                   ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
-                                     LOG( z_p / z0(j,i) ) +                   &
-                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
-                                                                   )
-                ELSE
-                   ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
-                                 LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
-                                                                   )
-                ENDIF
+                b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p )
+
+                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (         &
+                          LOG( z_p / z0(j,i) ) -                           &
+                          2.0 * LOG( ( 1.0 + a ) / ( 1.0 + b ) ) )
              ENDIF
 
@@ -168,25 +173,21 @@
 !
 !--          Unstable stratification
-             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
-             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * 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 ( a == 1.0  .OR.  b == 1.0 )  THEN
-                us(j,i) = kappa * uv_total / (                                &
-                                     LOG( z_p / z0(j,i) ) +                   &
-                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
-                                             )
-             ELSE
-                us(j,i) = kappa * uv_total / (                               &
-                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
-                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
-                                             )
-             ENDIF
+             a = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
+             b = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) )
+
+             us(j,i) = kappa * uv_total / (                                  &
+                       LOG( z_p / z0(j,i) ) -                                &
+                       LOG( ( 1.0 + a )**2 * ( 1.0 + a**2 ) / (              &
+                            ( 1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +          &
+                            2.0 * ( ATAN( a ) - ATAN( b ) )                  &
+                                           )
           ENDIF
        ENDDO
     ENDDO
 
+!
+!-- Values of us at ghost point locations are needed for the evaluation of usws 
+!-- and vsws.
+    CALL exchange_horiz_2d( us )
 !
 !-- Compute u'w' for the total model domain.
@@ -212,23 +213,15 @@
 !
 !--          Unstable stratification
-             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
-             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * 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 ( a == 1.0  .OR.  B == 1.0 )  THEN
-                usws(j,i) = kappa * u(k+1,j,i) / (                           &
-                                        LOG( z_p / z0(j,i) ) +               &
-                                        5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
+             a = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
+             b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
+
+             usws(j,i) = kappa * u(k+1,j,i) / (                           &
+                         LOG( z_p / z0(j,i) ) -                           &
+                         LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
+                              (1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +      &
+                              2.0 * ( ATAN( a ) - ATAN( b ) )             &
                                                  )
-             ELSE
-                usws(j,i) = kappa * u(k+1,j,i) / (                           &
-                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
-                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
-                                                 )
-             ENDIF
           ENDIF
-          usws(j,i) = -usws(j,i) * ABS( usws(j,i) )
+          usws(j,i) = -usws(j,i) * 0.5 * ( us(j,i-1) + us(j,i) )
        ENDDO
     ENDDO
@@ -257,23 +250,15 @@
 !
 !--          Unstable stratification
-             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
-             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * 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 ( a == 1.0  .OR.  b == 1.0 )  THEN
-                vsws(j,i) = kappa * v(k+1,j,i) / (                           &
-                                        LOG( z_p / z0(j,i) ) +               &
-                                        5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
+             a = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
+             b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
+
+             vsws(j,i) = kappa * v(k+1,j,i) / (                           &
+                         LOG( z_p / z0(j,i) ) -                           &
+                         LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
+                              (1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +      &
+                              2.0 * ( ATAN( a ) - ATAN( b ) )             &
                                                  )
-             ELSE
-                vsws(j,i) = kappa * v(k+1,j,i) / (                           &
-                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
-                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
-                                                 )
-             ENDIF
           ENDIF
-          vsws(j,i) = -vsws(j,i) * ABS( vsws(j,i) )
+          vsws(j,i) = -vsws(j,i) * 0.5 * ( us(j-1,i) + us(j,i) )
        ENDDO
     ENDDO
@@ -317,20 +302,10 @@
 !
 !--                Unstable stratification
-                   a = SQRT( 1.0 - 16.0 * rif(j,i) )
-                   b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * 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 ( a == 1.0  .OR.  b == 1.0 )  THEN
-                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (         &
-                                     LOG( z_p / z0(j,i) ) +                   &
-                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
-                                                                    )
-                   ELSE
-                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (          &
-                                  LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
-                                                                    )
-                   ENDIF
+                   a = SQRT( 1.0 - 16.0 * rif(j,i) ) 
+                   b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p ) 
+ 
+                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (   &
+                             LOG( z_p / z0(j,i) ) -                    &
+                              2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
                 ENDIF
 
@@ -341,7 +316,6 @@
 
 !
-!-- Exchange the boundaries for u* and the momentum fluxes (fluxes only for
-!-- completeness's sake).
-    CALL exchange_horiz_2d( us )
+!-- Exchange the boundaries for the momentum fluxes (only for sake of 
+!-- completeness)
     CALL exchange_horiz_2d( usws )
     CALL exchange_horiz_2d( vsws )
Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 184)
+++ palm/trunk/SOURCE/production_e.f90	(revision 187)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Change: add 'minus' sign to fluxes obtained from subroutine wall_fluxes_e for
+!         consistency with subroutine wall_fluxes
 !
 ! Former revisions:
@@ -164,10 +165,10 @@
                       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)
-!                      dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
+                      dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
+!                      dudy = - wall_e_y(j,i) * usvs(k,j,i) / 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)
-!                      dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
+                      dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
+!                      dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
                    ELSE
                       dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
@@ -180,10 +181,10 @@
                       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)
-!                      dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
+                      dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
+!                      dvdx = - wall_e_x(j,i) * vsus(k,j,i) / 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)
-!                      dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
+                      dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
+!                      dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
                    ELSE
                       dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
@@ -220,8 +221,8 @@
 
                       IF ( wall_e_y(j,i) /= 0.0 )  THEN
-                         dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
-!                         dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
-                         dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
-!                         dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
+                         dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
+!                         dudy = - wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
+                         dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
+!                         dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
                       ELSE
                          dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
@@ -232,8 +233,8 @@
 
                       IF ( wall_e_x(j,i) /= 0.0 )  THEN
-                         dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
-!                         dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
-                         dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
-!                         dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
+                         dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
+!                         dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
+                         dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
+!                         dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
                       ELSE
                          dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
@@ -629,8 +630,8 @@
                 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)
+                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)
+                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) - &
@@ -643,8 +644,8 @@
                 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)
+                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)
+                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) - &
@@ -679,6 +680,6 @@
 
                 IF ( wall_e_y(j,i) /= 0.0 )  THEN
-                   dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
-                   dwdy = wall_e_y(j,i) * wsvs(k) / 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) - &
@@ -689,6 +690,6 @@
 
                 IF ( wall_e_x(j,i) /= 0.0 )  THEN
-                   dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
-                   dwdx = wall_e_x(j,i) * wsus(k) / 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) - &
Index: palm/trunk/SOURCE/wall_fluxes.f90
===================================================================
--- palm/trunk/SOURCE/wall_fluxes.f90	(revision 184)
+++ palm/trunk/SOURCE/wall_fluxes.f90	(revision 187)
@@ -3,5 +3,20 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Bugfix: change definition of us_wall from 1D to 2D:
+!         Modification of the evaluation of the vertical turbulent momentum 
+!         fluxes u'w' and v'w'; the first usws that is computed corresponds 
+!         to -u'w'/u* and not as priorily assumed to (-u'w')**0.5, the first 
+!         vsws that is computed corresponds to -v'w'/u* and not as priorily
+!         assumed to (-v'w')**0.5. Therefore, the intermediate result for 
+!         usws has to be multiplied by -u* instead by itself in order to 
+!         get u'w'. Accordingly, the intermediate result for vsws has to be 
+!         multiplied by -u* instead by itself in order to get v'w'.
+!         This requires the calculation of us_wall (and vel_total, u_i, v_i, ws)
+!         also in wall_fluxes_e.
+! Bugfix: storage of rifs to rifs_wall in wall_fluxes_e removed
+! Change: add 'minus' sign to fluxes produced by subroutine wall_fluxes_e for 
+!         consistency with subroutine wall_fluxes
+! Change: Modification of the integrated version of the profile function for 
+!         momentum for unstable stratification
 !
 ! Former revisions:
@@ -69,5 +84,5 @@
 !
 !--             All subsequent variables are computed for the respective
-!--             location where the relevant variable is defined
+!--             location where the respective flux is defined.
                 DO  k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
 
@@ -109,22 +124,13 @@
 !
 !--                   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
-
+                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+                      us_wall = kappa * vel_total / (                          &
+                           LOG( zp / z0(j,i) ) -                               &
+                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                    )
                    ENDIF
 
@@ -160,27 +166,16 @@
 !
 !--                   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
-
+                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+                      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) ) -                               &
+                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                                            )
                    ENDIF
-                   wall_flux(k,j,i) = -wall_flux(k,j,i) * ABS(wall_flux(k,j,i))
+                   wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
 
 !
@@ -226,5 +221,5 @@
 !
 !--    All subsequent variables are computed for the respective location where 
-!--    the relevant variable is defined
+!--    the respective flux is defined.
        DO  k = nzb_w, nzt_w
 
@@ -266,22 +261,13 @@
 !
 !--          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
-
+             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+             us_wall = kappa * vel_total / (                          &
+                  LOG( zp / z0(j,i) ) -                               &
+                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                           )
           ENDIF
 
@@ -315,27 +301,16 @@
 !
 !--          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
-
+             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+             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) ) -                               &
+                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                                   )
           ENDIF
-          wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) )
+          wall_flux(k) = -wall_flux(k) * us_wall
 
 !
@@ -371,5 +346,6 @@
 
        INTEGER ::  i, j, k, kk, wall_index
-       REAL    ::  a, b, c1, c2, h1, h2, vel_zp, zp
+       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
+                   ws, zp
 
        REAL ::  rifs
@@ -388,10 +364,8 @@
              IF ( wall(j,i) /= 0.0 )  THEN
 !
-!--             All subsequent variables are computed for the respective
-!--             location where the relevant variable is defined
+!--             All subsequent variables are computed for scalar locations.
                 DO  k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
-
-!
-!--                (1) Compute rifs
+!
+!--                (1) Compute rifs, u_i, v_i, and ws
                    IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
                       kk = nzb_diff_s_inner(j,i)-1
@@ -404,14 +378,41 @@
                                  )
 
-!
-!--                Skip (2) to (4) of wall_fluxes, because here rifs is
-!--                already available from (1)
-
+                   u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
+                   v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
+                   ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
+!
+!--                (2) Compute wall-parallel absolute velocity vel_total and
+!--                interpolate appropriate velocity component vel_zp.
+                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
+                   vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
+!
+!--                (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 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+                      us_wall = kappa * vel_total / (                          &
+                           LOG( zp / z0(j,i) ) -                               &
+                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                    )
+                   ENDIF
+
+!
+!--                Skip step (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
@@ -425,29 +426,15 @@
 !
 !--                   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
-
+                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+                      wall_flux(k,j,i) = kappa * vel_zp / (                    &
+                           LOG( zp / z0(j,i) ) -                               &
+                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                          )
                    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
+                   wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
 
                 ENDDO
@@ -476,5 +463,6 @@
 
        INTEGER ::  i, j, k, kk, nzb_w, nzt_w, wall_index
-       REAL    ::  a, b, c1, c2, h1, h2, vel_zp, zp
+       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
+                   ws, zp
 
        REAL ::  rifs
@@ -488,10 +476,9 @@
 
 !
-!--    All subsequent variables are computed for the respective location where 
-!--    the relevant variable is defined
+!--    All subsequent variables are computed for scalar locations.
        DO  k = nzb_w, nzt_w
 
 !
-!--       (1) Compute rifs
+!--       (1) Compute rifs, u_i, v_i, and ws
           IF ( k == nzb_w )  THEN
              kk = nzb_w
@@ -504,10 +491,44 @@
                         )
 
-!
-!--       Skip (2) to (4) of wall_fluxes, because here rifs is already available
-!--       from (1)
-
+          u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
+          v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
+          ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
+!
+!--       (2) Compute wall-parallel absolute velocity vel_total and
+!--       interpolate appropriate velocity component vel_zp.
+          vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
+          vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
+!
+!--       (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 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+             us_wall = kappa * vel_total / (                          &
+                  LOG( zp / z0(j,i) ) -                               &
+                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                           )
+          ENDIF
+
+!
+!--       Skip step (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')
+!--       First interpolate the velocity (this is different from 
+!--       subroutine wall_fluxes because fluxes in subroutine 
+!--       wall_fluxes_e are defined at scalar locations).
           vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
                                  b * ( v(k,j,i) + v(k,j+1,i) ) +  &
@@ -525,29 +546,15 @@
 !
 !--          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
-
+             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+             wall_flux(k) = kappa * vel_zp / (                        &
+                  LOG( zp / z0(j,i) ) -                               &
+                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                 )
           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
+          wall_flux(k) = - wall_flux(k) * us_wall
 
        ENDDO
