tau_u_s.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_Tau_U_s                                            C
!  Purpose: Cross terms in the gradient of stress in U_s momentum      C
!                                                                      C
!  Comments: This routine calculates the components of the solids      C
!  phase viscous stress tensor of the u-momentum equation that cannot  C
!  be cast in the form: mu.grad(u). These components are stored in the C
!  passed variable, which is then used as a source of the u-momentum   C
!  equation.                                                           C
!  For greater details see calc_tau_u_g.                               C
!                                                                      C
!                                                                      C
!  Author: M. Syamlal                                 Date: 19-DEC-96  C
!  Reviewer:                                          Date:            C
!                                                                      C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_tau_u_s(lTAU_U_S)

! Modules
!---------------------------------------------------------------------//
      USE param, only: dimension_3, dimension_m
      USE param1, only: zero
      USE physprop, only: smax, mmax
      USE fldvar, only: v_s
      USE fldvar, only: ep_s
      USE toleranc, only: dil_ep_s
      USE run, only: kt_type_enum, ghd_2007
      USE run, only: shear
      USE vshear, only: vsh
      USE cutcell, only: cartesian_grid, cg_safe_mode

      USE compar
      USE geometry
      USE indices
      USE functions
      USE sendrecv
      USE fun_avg
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! TAU_U_s
      DOUBLE PRECISION, INTENT(OUT) :: lTAU_U_s(dimension_3, dimension_m)

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: IJK, I, IJKE
! Phase index
      INTEGER :: M, L
! Average volume fraction
      DOUBLE PRECISION :: EPSA, EPStmp
! Source terms (Surface)
      DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
! Cylindrical source terms
      DOUBLE PRECISION :: Vtzb
!---------------------------------------------------------------------//


      DO m = 1, mmax
         IF(kt_type_enum == ghd_2007 .AND. m /= mmax) cycle

         IF (shear) THEN
!!$omp  parallel do private(IJK)
            DO ijk = ijkstart3, ijkend3
               IF (fluid_at(ijk)) THEN
                  v_s(ijk,m)=v_s(ijk,m)+vsh(ijk)
               ENDIF
            ENDDO
         ENDIF

!!$omp  parallel do &
!!$omp  private(I, IJK, IJKE,                                         &
!!$omp&         EPSA, EPStmp, SSX, SSY, SSZ, SBV, VTZB)               &
!!$omp&  schedule(static)
         DO ijk = ijkstart3, ijkend3

! Skip walls where some values are undefined.
            IF(wall_at(ijk)) cycle
            i = i_of(ijk)
            ijke = east_of(ijk)

            IF (kt_type_enum == ghd_2007) THEN
               epstmp = zero
               DO l = 1, smax
                 epstmp = epstmp + avg_x(ep_s(ijk,l),ep_s(ijke,l),i)
               ENDDO
               epsa = epstmp
            ELSE
               epsa = avg_x(ep_s(ijk,m),ep_s(ijke,m),i)
            ENDIF

            IF ( .NOT.sip_at_e(ijk) .AND. epsa>dil_ep_s) THEN

               IF ((.NOT.cartesian_grid) .OR. cg_safe_mode(3) ==1) THEN
                  CALL calc_reg_tau_u_s(ijk, m, ssx, ssy, ssz, sbv, vtzb)
               ELSE
                  vtzb = zero   ! no cylindrical terms in CG
                  CALL calc_cg_tau_u_s(ijk, m, ssx, ssy, ssz, sbv)
               ENDIF

! Add the terms
               ltau_u_s(ijk,m) = sbv + ssx + ssy + ssz + vtzb*vol_u(ijk)
            ELSE
               ltau_u_s(ijk,m) = zero
            ENDIF   ! end if/else .not.sip_at_e .and. epsa>dil_ep_s
         ENDDO   ! end do ijk
!!$omp end parallel do

         IF (shear) THEN
!!$omp  parallel do private(IJK)
            DO ijk = ijkstart3, ijkend3
               IF (fluid_at(ijk)) THEN
                  v_s(ijk,m)=v_s(ijk,m)-vsh(ijk)
               ENDIF
            ENDDO
         ENDIF

      ENDDO   ! end do m=1,mmax

      call send_recv(ltau_u_s,2)
      RETURN
      END SUBROUTINE calc_tau_u_s


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: calc_reg_tau_u_s                                        C
!  Purpose: Cross terms in the gradient of stress in U_s momentum      C
!  based on NON cartesian grid cut cell.                               C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_reg_tau_u_s(IJK, M, SSX, SSY, SSZ, SBV, VTZB)

! Modules
!---------------------------------------------------------------------//
      USE param1, only: zero, half
      USE fldvar, only: u_s, v_s, w_s
      USE visc_s, only: epmu_s, eplambda_s
      USE visc_s, only: trd_s

      USE compar
      USE geometry
      USE indices
      USE functions
      USE fun_avg

      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! current ijk index
      INTEGER, INTENT(IN) :: IJK
! solids phase index
      INTEGER, INTENT(IN) :: M
! Source terms (surface)
      DOUBLE PRECISION, INTENT(OUT) :: SSX
      DOUBLE PRECISION, INTENT(OUT) :: SSY
      DOUBLE PRECISION, INTENT(OUT) :: SSZ
      DOUBLE PRECISION, INTENT(OUT) :: SBV
! Cylindrical source terms
      DOUBLE PRECISION, INTENT(OUT) :: VTZB

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: I, J, K, JM, KM, IP
      INTEGER :: IJKE, IJKN, IJKS, IJKT, IJKB
      INTEGER :: IJKNE, IJKSE, IJKTE, IJKBE
      INTEGER :: IPJK, IMJK, IJKM, IJMK
      INTEGER :: IPJMK, IPJKM
! Average viscosity
      DOUBLE PRECISION :: MU_ste, MU_sbe, MUSA
! Average dW/Xdz
      DOUBLE PRECISION :: dWoXdz
!---------------------------------------------------------------------//

      i = i_of(ijk)
      ip = ip1(i)
      j = j_of(ijk)
      jm = jm1(j)
      k = k_of(ijk)
      km = km1(k)

      ipjk = ip_of(ijk)
      imjk = im_of(ijk)
      ijmk = jm_of(ijk)
      ijkm = km_of(ijk)
      ipjkm = ip_of(ijkm)
      ipjmk = jm_of(ipjk)

      ijke = east_of(ijk)
      ijkn = north_of(ijk)
      ijkne = east_of(ijkn)
      ijks = south_of(ijk)
      ijkse = east_of(ijks)
      ijkt = top_of(ijk)
      ijkte = east_of(ijkt)
      ijkb = bottom_of(ijk)
      ijkbe = east_of(ijkb)

! Surface forces at i+1/2, j, k
! bulk viscosity term
! combines part of 1/x d/dx (x.tau_xx) xdxdydz and -tau_zz/x xdxdydz =>
! combines 1/x d/dx (x.lambda.trcD) xdxdydz - (lambda/x.trcD) xdxdydz =>
!              d/dx (lambda.trcD) xdxdydz
! delta (lambda.trcD)Ap |E-W : at (i+1 - i-1), j, k
      sbv = (eplambda_s(ijke,m)*trd_s(ijke,m)-&
             eplambda_s(ijk,m)*trd_s(ijk,m))*ayz(ijk)

! shear stress terms at i+1/2, j, k
! part of 1/x d/dx(x.tau_xx) xdxdydz =>
!         1/x d/dx (x.mu.du/dx) xdxdydz =>
! delta (mu du/dx)Ayz |E-W : at (i+1 - i-1), j, k
      ssx = epmu_s(ijke,m)*(u_s(ipjk,m)-u_s(ijk,m))*&
               odx(ip)*ayz_u(ijk) - &
            epmu_s(ijk,m)*(u_s(ijk,m)-u_s(imjk,m))*&
               odx(i)*ayz_u(imjk)

! part of d/dy (tau_xy) xdxdydz =>
!         d/dy (mu.dv/dx) xdxdydz =>
! delta (mu.dv/dx)Axz |N-S : at i+1/2, (j+1/2 - j-1/2), k
      ssy = avg_x_h(avg_y_h(epmu_s(ijk,m),epmu_s(ijkn,m),j),&
                    avg_y_h(epmu_s(ijke,m),epmu_s(ijkne,m),j),i)*&
               (v_s(ipjk,m)-v_s(ijk,m))*odx_e(i)*axz_u(ijk) - &
            avg_x_h(avg_y_h(epmu_s(ijks,m),epmu_s(ijk,m),jm),&
                    avg_y_h(epmu_s(ijkse,m),epmu_s(ijke,m),jm),i)*&
               (v_s(ipjmk,m)-v_s(ijmk,m))*odx_e(i)*axz_u(ijmk)

! part of 1/x d/dz (tau_xz) xdxdydz =>
!         1/x d/dz (mu.dw/dx) xdxdydz =>
! delta (mu.dw/dx)Axy |T-B : at i+1/2, j, (k+1/2 - k-1/2)
      mu_ste = avg_x_h(avg_z_h(epmu_s(ijk,m),epmu_s(ijkt,m),k),&
                       avg_z_h(epmu_s(ijke,m),epmu_s(ijkte,m),k),i)
      mu_sbe = avg_x_h(avg_z_h(epmu_s(ijkb,m),epmu_s(ijk,m),km),&
                       avg_z_h(epmu_s(ijkbe,m),epmu_s(ijke,m),km),i)
      ssz = mu_ste*(w_s(ipjk,m)-w_s(ijk,m))*odx_e(i)*&
                        axy_u(ijk) - &
            mu_sbe*(w_s(ipjkm,m)-w_s(ijkm,m))*odx_e(i)*&
                        axy_u(ijkm)


! Special terms for cylindrical coordinates
      IF (cylindrical) THEN
! part of 1/x d/dz (tau_xz) xdxdydz =>
!         1/x d/dz (mu.(-w/x)) xdxdydz =>
! delta (mu(-w/x))Axy |T-B : at i+1/2, j, (k+1/2- k-1/2)
         ssz = ssz - (&
               mu_ste*(half*(w_s(ipjk,m)+w_s(ijk,m))*&
                  ox_e(i))*axy_u(ijk)-&
               mu_sbe*(half*(w_s(ipjkm,m)+w_s(ijkm,m))*&
                  ox_e(i))*axy_u(ijkm))

! part of -tau_zz/x xdxdydz =>
!         -(2.mu/x).(1/x).dw/dz xdxdydz
! delta (2.mu/x.1/x.dw/dz)Vp |p : at i+1/2, j, k
         musa = avg_x(epmu_s(ijk,m),epmu_s(ijke,m),i)
         dwoxdz = half*&
                  ((w_s(ijk,m)-w_s(ijkm,m))*ox(i)*odz(k)+&
                   (w_s(ipjk,m)-w_s(ipjkm,m))*ox(ip)*odz(k))
         vtzb = -2.d0*musa*ox_e(i)*dwoxdz
      ELSE
         vtzb = zero
      ENDIF

      RETURN
      END SUBROUTINE calc_reg_tau_u_s


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: calc_cg_tau_u_s                                         C
!  Purpose: Cross terms in the gradient of stress in U_s momentum      C
!  based on cartesian grid cut cell.                                   C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 01-Jul-09  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_cg_tau_u_s(IJK, M, SSX, SSY, SSZ, SBV)

! Modules
!---------------------------------------------------------------------//

      USE compar
      USE fldvar, only: u_s, v_s, w_s
      USE fun_avg
      USE functions
      USE geometry
      USE indices
      USE param1, only: zero, half, undefined
      USE visc_s, only: epmu_s, eplambda_s
      USE visc_s, only: trd_s

! for cutcell:
! wall velocity for slip conditions
      USE bc, only: bc_hw_s, bc_uw_s, bc_vw_s, bc_ww_s, cg_nsw, cg_fsw, cg_psw, cg_mi, none, bc_type_enum
      USE quadric
      USE cutcell

! Dummy arguments
!---------------------------------------------------------------------//
! current ijk index
      INTEGER, INTENT(IN) :: IJK
! solids phase index
      INTEGER, INTENT(IN) :: M
! Source terms (surface)
      DOUBLE PRECISION, INTENT(OUT) :: SSX
      DOUBLE PRECISION, INTENT(OUT) :: SSY
      DOUBLE PRECISION, INTENT(OUT) :: SSZ
      DOUBLE PRECISION, INTENT(OUT) :: SBV

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: I, J, K, JM, KM, IP
      INTEGER :: IJKE, IJKN, IJKS, IJKT, IJKB
      INTEGER :: IJKNE, IJKSE, IJKTE, IJKBE
      INTEGER :: IPJK, IMJK, IJKM, IJMK
      INTEGER :: IPJMK, IPJKM
! Average viscosity
      DOUBLE PRECISION :: MU_ste, MU_sbe

! for cartesian grid implementation:
      DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
      LOGICAL :: V_NODE_AT_NE, V_NODE_AT_NW, V_NODE_AT_SE, V_NODE_AT_SW
      LOGICAL :: W_NODE_AT_TE, W_NODE_AT_TW, W_NODE_AT_BE, W_NODE_AT_BW
      DOUBLE PRECISION :: dvdx_at_N, dvdx_at_S
      DOUBLE PRECISION :: dwdx_at_T, dwdx_at_B
      DOUBLE PRECISION :: Xi, Yi, Zi, Vi, Wi, Sx, Sy, Sz
      DOUBLE PRECISION :: MU_S_CUT, SSY_CUT, SSZ_CUT
      DOUBLE PRECISION :: UW_s, VW_s, WW_s
      INTEGER :: BCV
      INTEGER :: BCT
!---------------------------------------------------------------------//

      i = i_of(ijk)
      ip = ip1(i)
      j = j_of(ijk)
      jm = jm1(j)
      k = k_of(ijk)
      km = km1(k)

      ipjk = ip_of(ijk)
      imjk = im_of(ijk)
      ijmk = jm_of(ijk)
      ijkm = km_of(ijk)
      ipjkm = ip_of(ijkm)
      ipjmk = jm_of(ipjk)

      ijke = east_of(ijk)
      ijkn = north_of(ijk)
      ijkne = east_of(ijkn)
      ijks = south_of(ijk)
      ijkse = east_of(ijks)
      ijkt = top_of(ijk)
      ijkte = east_of(ijkt)
      ijkb = bottom_of(ijk)
      ijkbe = east_of(ijkb)


! bulk viscosity term
      sbv = (eplambda_s(ijke,m)*trd_s(ijke,m)) * ayz_u(ijk) - &
            (eplambda_s(ijk,m) *trd_s(ijk,m) ) * ayz_u(imjk)

! shear stress terms
      IF(.NOT.cut_u_cell_at(ijk))   THEN
         ssx = epmu_s(ijke,m)*(u_s(ipjk,m)-u_s(ijk,m))*&
                  oneodx_e_u(ijk)*ayz_u(ijk) - &
               epmu_s(ijk,m)*(u_s(ijk,m)-u_s(imjk,m))*&
                  oneodx_e_u(imjk)*ayz_u(imjk)
         ssy = avg_x_h(avg_y_h(epmu_s(ijk,m),epmu_s(ijkn,m),j),&
                       avg_y_h(epmu_s(ijke,m),epmu_s(ijkne,m),j),i)*&
                  (v_s(ipjk,m)-v_s(ijk,m))*&
                  oneodx_e_v(ijk)*axz_u(ijk) - &
               avg_x_h(avg_y_h(epmu_s(ijks,m),epmu_s(ijk,m),jm),&
                       avg_y_h(epmu_s(ijkse,m),epmu_s(ijke,m),jm),i)*&
                  (v_s(ipjmk,m)-v_s(ijmk,m))*&
                  oneodx_e_v(ijmk)*axz_u(ijmk)
         IF(do_k) THEN
            mu_ste = avg_x_h(avg_z_h(epmu_s(ijk,m),epmu_s(ijkt,m),k),&
                             avg_z_h(epmu_s(ijke,m),epmu_s(ijkte,m),k),i)
            mu_sbe = avg_x_h(avg_z_h(epmu_s(ijkb,m),epmu_s(ijk,m),km),&
                             avg_z_h(epmu_s(ijkbe,m),epmu_s(ijke,m),km),i)
            ssz = mu_ste*(w_s(ipjk,m)-w_s(ijk,m))*&
                     oneodx_e_w(ijk)*axy_u(ijk) - &
                  mu_sbe*(w_s(ipjkm,m)-w_s(ijkm,m))*&
                     oneodx_e_w(ijkm)*axy_u(ijkm)
         ELSE
            ssz = zero
         ENDIF

! cut cell modifications
!---------------------------------------------------------------------//
      ELSE

         bcv = bc_u_id(ijk)
         IF(bcv > 0 ) THEN
            bct = bc_type_enum(bcv)
         ELSE
            bct = none
         ENDIF

         SELECT CASE (bct)
            CASE (cg_nsw,cg_mi)
               cut_tau_us = .true.
               noc_us     = .true.
               uw_s = zero
               vw_s = zero
               ww_s = zero
            CASE (cg_fsw)
               cut_tau_us = .false.
               noc_us     = .false.
               uw_s = zero
               vw_s = zero
               ww_s = zero
            CASE(cg_psw)
               IF(bc_hw_s(bc_u_id(ijk),m)==undefined) THEN   ! same as NSW
                  cut_tau_us = .true.
                  noc_us     = .true.
                  uw_s = bc_uw_s(bcv,m)
                  vw_s = bc_vw_s(bcv,m)
                  ww_s = bc_ww_s(bcv,m)
               ELSEIF(bc_hw_s(bc_u_id(ijk),m)==zero) THEN   ! same as FSW
                  cut_tau_us = .false.
                  noc_us     = .false.
                  uw_s = zero
                  vw_s = zero
                  ww_s = zero
               ELSE                              ! partial slip
                  cut_tau_us = .false.
                  noc_us     = .false.
               ENDIF
            CASE (none)
! equivalent of setting tau_u_s to zero at this ijk cell
               ssx = zero
               ssy = zero
               ssz = zero
               sbv = zero
               RETURN
         END SELECT

         IF(cut_tau_us) THEN
            mu_s_cut = (vol(ijk)*epmu_s(ijk,m) + &
                        vol(ipjk)*epmu_s(ijke,m))/&
                       (vol(ijk) + vol(ipjk))
         ELSE
            mu_s_cut = zero
         ENDIF


! SSX:
         CALL get_del_h(ijk, 'U_MOMENTUM', x_u(ijk), &
            y_u(ijk), z_u(ijk), del_h, nx, ny, nz)
         ssx = epmu_s(ijke,m)*(u_s(ipjk,m)-u_s(ijk,m))*&
                  oneodx_e_u(ijk)*ayz_u(ijk) - &
               epmu_s(ijk,m)*(u_s(ijk,m)-u_s(imjk,m))*&
                  oneodx_e_u(imjk)*ayz_u(imjk) - &
               mu_s_cut * (u_s(ijk,m) - uw_s) / del_h * &
                  (nx**2) * area_u_cut(ijk)

! SSY:
         v_node_at_ne = ((.NOT.blocked_v_cell_at(ipjk)).AND.&
                         (.NOT.wall_v_at(ipjk)))
         v_node_at_nw = ((.NOT.blocked_v_cell_at(ijk)).AND.&
                         (.NOT.wall_v_at(ijk)))
         v_node_at_se = ((.NOT.blocked_v_cell_at(ipjmk)).AND.&
                         (.NOT.wall_v_at(ipjmk)))
         v_node_at_sw = ((.NOT.blocked_v_cell_at(ijmk)).AND.&
                         (.NOT.wall_v_at(ijmk)))

         IF(v_node_at_ne.AND.v_node_at_nw) THEN
            vi = half * (v_s(ipjk,m) + v_s(ijk,m))
            xi = half * (x_v(ipjk) + x_v(ijk))
            yi = half * (y_v(ipjk) + y_v(ijk))
            zi = half * (z_v(ipjk) + z_v(ijk))
            sx = x_v(ipjk) - x_v(ijk)
            sy = y_v(ipjk) - y_v(ijk)
            sz = z_v(ipjk) - z_v(ijk)
            CALL get_del_h(ijk, 'U_MOMENTUM', xi, yi, zi, &
               del_h, nx, ny, nz)
            dvdx_at_n = (v_s(ipjk,m) - v_s(ijk,m)) * &
                         oneodx_e_v(ijk)
            IF(noc_us) dvdx_at_n = dvdx_at_n - ((vi-vw_s) *&
                          oneodx_e_v(ijk) /del_h*(sy*ny+sz*nz))
         ELSE
            dvdx_at_n =  zero
         ENDIF

         IF(v_node_at_se.AND.v_node_at_sw) THEN
            vi = half * (v_s(ipjmk,m) + v_s(ijmk,m))
            xi = half * (x_v(ipjmk) + x_v(ijmk))
            yi = half * (y_v(ipjmk) + y_v(ijmk))
            zi = half * (z_v(ipjmk) + z_v(ijmk))
            sx = x_v(ipjmk) - x_v(ijmk)
            sy = y_v(ipjmk) - y_v(ijmk)
            sz = z_v(ipjmk) - z_v(ijmk)
            CALL get_del_h(ijk, 'U_MOMENTUM', xi, yi, zi, &
               del_h, nx, ny, nz)
            dvdx_at_s = (v_s(ipjmk,m)-v_s(ijmk,m))*oneodx_e_v(ijmk)
            IF(noc_us) dvdx_at_s = dvdx_at_s - ((vi-vw_s) * &
                          oneodx_e_v(ijmk)/del_h*(sy*ny+sz*nz))
         ELSE
            dvdx_at_s =  zero
         ENDIF

         IF(v_node_at_nw) THEN
            CALL get_del_h(ijk, 'U_MOMENTUM', x_v(ijk), &
               y_v(ijk), z_v(ijk), del_h, nx, ny, nz)
            ssy_cut = -mu_s_cut * (v_s(ijk,m) - vw_s)/del_h*&
                       (nx*ny) * area_u_cut(ijk)
         ELSE
            ssy_cut =  zero
         ENDIF

         ssy = avg_x_h(avg_y_h(epmu_s(ijk,m),epmu_s(ijkn,m),j),&
                       avg_y_h(epmu_s(ijke,m),epmu_s(ijkne,m),j),i)*&
                  dvdx_at_n*axz_u(ijk) - &
               avg_x_h(avg_y_h(epmu_s(ijks,m),epmu_s(ijk,m),jm),&
                       avg_y_h(epmu_s(ijkse,m),epmu_s(ijke,m),jm),i)*&
                  dvdx_at_s*axz_u(ijmk) + ssy_cut

! SSZ:
         IF(do_k) THEN
            w_node_at_te = ((.NOT.blocked_w_cell_at(ipjk)).AND.&
                            (.NOT.wall_w_at(ipjk)))
            w_node_at_tw = ((.NOT.blocked_w_cell_at(ijk)).AND.&
                            (.NOT.wall_w_at(ijk)))
            w_node_at_be = ((.NOT.blocked_w_cell_at(ipjkm)).AND.&
                            (.NOT.wall_w_at(ipjkm)))
            w_node_at_bw = ((.NOT.blocked_w_cell_at(ijkm)).AND.&
                            (.NOT.wall_w_at(ijkm)))

            IF(w_node_at_te.AND.w_node_at_tw) THEN
               wi = half * (w_s(ipjk,m) + w_s(ijk,m))
               xi = half * (x_w(ipjk) + x_w(ijk))
               yi = half * (y_w(ipjk) + y_w(ijk))
               zi = half * (z_w(ipjk) + z_w(ijk))
               sx = x_w(ipjk) - x_w(ijk)
               sy = y_w(ipjk) - y_w(ijk)
               sz = z_w(ipjk) - z_w(ijk)
               CALL get_del_h(ijk, 'U_MOMENTUM', xi, yi, zi, &
                  del_h, nx, ny, nz)
               dwdx_at_t = (w_s(ipjk,m) - w_s(ijk,m)) * &
                           oneodx_e_w(ijk)
               IF(noc_us) dwdx_at_t = dwdx_at_t - ((wi-ww_s)*&
                             oneodx_e_w(ijk)/del_h*(sy*ny+sz*nz))
            ELSE
               dwdx_at_t =  zero
            ENDIF

            IF(w_node_at_be.AND.w_node_at_bw) THEN
               wi = half * (w_s(ipjkm,m) + w_s(ijkm,m))
               xi = half * (x_w(ipjkm) + x_w(ijkm))
               yi = half * (y_w(ipjkm) + y_w(ijkm))
               zi = half * (z_w(ipjkm) + z_w(ijkm))
               sx = x_w(ipjkm) - x_w(ijkm)
               sy = y_w(ipjkm) - y_w(ijkm)
               sz = z_w(ipjkm) - z_w(ijkm)
               CALL get_del_h(ijk, 'U_MOMENTUM', xi, yi, &
                  zi, del_h, nx, ny, nz)
               dwdx_at_b = (w_s(ipjkm,m) - w_s(ijkm,m)) * &
                           oneodx_e_w(ijkm)
               IF(noc_us) dwdx_at_b = dwdx_at_b - ((wi-ww_s)*&
                             oneodx_e_w(ijkm)/del_h*(sy*ny+sz*nz))
            ELSE
               dwdx_at_b =  zero
            ENDIF

            IF(w_node_at_tw) THEN
               CALL get_del_h(ijk, 'U_MOMENTUM', x_w(ijk), &
                  y_w(ijk), z_w(ijk), del_h, nx, ny, nz)
               ssz_cut = -mu_s_cut * (w_s(ijk,m) - ww_s) / &
                         del_h * (nx*nz) * area_u_cut(ijk)
            ELSE
               ssz_cut =  zero
            ENDIF

            mu_ste = avg_x_h(avg_z_h(epmu_s(ijk,m),epmu_s(ijkt,m),k),&
                             avg_z_h(epmu_s(ijke,m),epmu_s(ijkte,m),k),i)
            mu_sbe = avg_x_h(avg_z_h(epmu_s(ijkb,m),epmu_s(ijk,m),km),&
                             avg_z_h(epmu_s(ijkbe,m),epmu_s(ijke,m),km),i)
            ssz = mu_ste*dwdx_at_t*axy_u(ijk) -  &
                  mu_sbe*dwdx_at_b*axy_u(ijkm) + ssz_cut
         ELSE
            ssz = zero
         ENDIF  ! end if do_k

      ENDIF  ! end if/else cut_cell

      RETURN
      END SUBROUTINE calc_cg_tau_u_s