source: palm/trunk/SOURCE/wall_fluxes.f90 @ 54

 SUBROUTINE wall_fluxes( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
!------------------------------------------------------------------------------!
! Actual revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $Id: wall_fluxes.f90 54 2007-03-08 00:00:02Z raasch $
! 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, 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
             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) = 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



 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