Index: palm/trunk/SOURCE/wall_fluxes.f90
===================================================================
--- palm/trunk/SOURCE/wall_fluxes.f90	(revision 52)
+++ palm/trunk/SOURCE/wall_fluxes.f90	(revision 52)
@@ -0,0 +1,143 @@
+ SUBROUTINE wall_fluxes( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
+!------------------------------------------------------------------------------!
+! Actual revisions:
+! -----------------
+! 
+!
+! Former revisions:
+! -----------------
+! $Id$
+! Initial version (2007/03/07)
+!
+! Description:
+! ------------
+! Calculates momentum fluxes at vertical walls 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, nzb_w, nzt_w
+    REAL    ::  a, b, c1, c2, h1, h2, delta_p
+
+    REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
+
+    REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
+
+
+    delta_p   = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
+    wall_flux = 0.0
+
+!
+!-- 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,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
+!
+!--    (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  &
+                                        )
+       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) ) )
+!
+!--       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( delta_p / z0(j,i) ) +                         &
+                       5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p &
+                                                   )
+          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 delta_p/L (corresponds to neutral Richardson flux number 
+!--        rifs)
+       rifs = -1.0 * delta_p * 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( delta_p / z0(j,i) ) +                                &
+               5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p        &
+                                                                         )
+       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) ) )
+!
+!--       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( 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( (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,NINT(a+2*b+3*c1+4*c2)) = rifs
+
+    ENDDO
+ END SUBROUTINE wall_fluxes
