tau_v_s.f

Go to the documentation of this file.

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

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

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

! Dummy arguments
!---------------------------------------------------------------------//
! TAU_V_s
      DOUBLE PRECISION, INTENT(OUT) :: lTAU_V_s(dimension_3, dimension_m)

! Local Variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: IJK, J, IJKN
! Additional indices needed for shear case
      INTEGER :: I, IJKE, IJKW, IJKNE, IJKNW
! Phase index
      INTEGER :: M, L
! Average volume fraction
      DOUBLE PRECISION :: EPSA, EPStmp
! Source terms (Surface)
      DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
! Shearing variables
      DOUBLE PRECISION :: Source_diff, Diffco_e, Diffco_w
!---------------------------------------------------------------------//

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

!!$omp  parallel do &
!!$omp  private(J, IJK, IJKN, I, IJKE, IJKW, IJKNE, IJKNW,            &
!!$omp&         EPSA, EPStmp, SSX, SSY, SSZ, SBV,                     &
!!$omp&         Source_diff, Diffco_e, Diffco_w)
!!$omp&  schedule(static)
         DO ijk = ijkstart3, ijkend3

! Skip walls where some values are undefined.
            IF(wall_at(ijk)) cycle
            j = j_of(ijk)
            ijkn = north_of(ijk)

            IF (kt_type_enum == ghd_2007) THEN
               epstmp = zero
               DO l = 1, smax
                  epstmp = epstmp + avg_y(ep_s(ijk,l),ep_s(ijkn,l),j)
               ENDDO
               epsa = epstmp
            ELSE
               epsa = avg_y(ep_s(ijk,m),ep_s(ijkn,m),j)
            ENDIF

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

               IF((.NOT.cartesian_grid).OR.(cg_safe_mode(4)==1)) THEN
                  CALL calc_reg_tau_v_s(ijk, m, ssx, ssy, ssz, sbv)
               ELSE
                  CALL calc_cg_tau_v_s(ijk, m, ssx, ssy, ssz, sbv)
               ENDIF


! Source terms from shear stress
               IF (shear) THEN
                  i = i_of(ijk)
                  ijke = east_of(ijk)
                  ijkn = north_of(ijk)
                  ijkne = east_of(ijkn)
                  ijkw = west_of(ijk)
                  ijknw = north_of(ijkw)
                  diffco_e = avg_y_h((avg_x_h(epmu_s(ijk,m),epmu_s(ijke,m),i)),&
                                     (avg_x_h(epmu_s(ijkn,m),epmu_s(ijkne,m),i)),j)*ayz_v(ijk)
                  diffco_w=avg_y_h((avg_x_h(epmu_s(ijk,m),epmu_s(ijkw,m),i)),&
                                   (avg_x_h(epmu_s(ijkn,m),epmu_s(ijknw,m),i)),j)*ayz_v(ijkw)
                  source_diff=(2d0*v_sh/xlength)*(diffco_e-diffco_w)
               ELSE
                  source_diff=0d0
               ENDIF

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

      ENDDO   ! end do m=1,mmax

      call send_recv(ltau_v_s,2)
      RETURN
      END SUBROUTINE calc_tau_v_s


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

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

      USE fun_avg
      USE compar
      USE geometry
      USE indices
      USE functions
      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

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: I, J, K, IM, KM, JP
      INTEGER :: IJKE, IJKW, IJKN, IJKT, IJKB
      INTEGER :: IJKNE, IJKNW, IJKTN, IJKBN
      INTEGER :: IJPK, IJMK, IMJK, IJKM
      INTEGER :: IMJPK, IJPKM
!---------------------------------------------------------------------//

      j = j_of(ijk)
      jp = jp1(j)
      i = i_of(ijk)
      im = im1(i)
      k = k_of(ijk)
      km = km1(k)

      ijpk = jp_of(ijk)
      ijmk = jm_of(ijk)
      ijkm = km_of(ijk)
      imjk = im_of(ijk)
      imjpk = im_of(ijpk)
      ijpkm = jp_of(ijkm)

      ijke = east_of(ijk)
      ijkn = north_of(ijk)
      ijkne = east_of(ijkn)
      ijkw = west_of(ijk)
      ijknw = north_of(ijkw)
      ijkt = top_of(ijk)
      ijktn = north_of(ijkt)
      ijkb = bottom_of(ijk)
      ijkbn = north_of(ijkb)

! Surface forces

! bulk viscosity term
! part of d/dy (tau_yy) xdxdydz =>
!         d/dy (lambda.trcD) xdxdydz =>
! delta (lambda.trcD)Ap|N-S  : at (i, j+1 - j-1, k)
      sbv = (eplambda_s(ijkn,m)*trd_s(ijkn,m)-&
             eplambda_s(ijk,m)*trd_s(ijk,m))*axz(ijk)

! shear stress terms at i, j+1/2, k
! part of 1/x d/dx(x.tau_xy) xdxdydz =>
!         1/x d/dx (x.mu.du/dy) xdxdydz =>
! delta (x.mu.du/dy)Ayz |E-W : at (i+1/2 - i-1/2, j+1/2, k)
      ssx = avg_y_h(avg_x_h(epmu_s(ijk,m),epmu_s(ijke,m),i),&
                    avg_x_h(epmu_s(ijkn,m),epmu_s(ijkne,m),i),j)*&
               (u_s(ijpk,m)-u_s(ijk,m))*ody_n(j)*ayz_v(ijk) - &
            avg_y_h(avg_x_h(epmu_s(ijkw,m),epmu_s(ijk,m),im),&
                    avg_x_h(epmu_s(ijknw,m),epmu_s(ijkn,m),im),j)*&
               (u_s(imjpk,m)-u_s(imjk,m))*ody_n(j)*ayz_v(imjk)

! part of d/dy (tau_xy) xdxdydz =>
!         d/dy (mu.dv/dy) xdxdydz =>
! delta (mu.dv/dx)Axz |N-S : at (i, j+1 - j-1, k)
      ssy = epmu_s(ijkn,m)*(v_s(ijpk,m)-v_s(ijk,m))*ody(jp)*axz_v(ijk) - &
            epmu_s(ijk,m)*(v_s(ijk,m)-v_s(ijmk,m))*ody(j)*axz_v(ijmk)

! part of 1/x d/dz (tau_xz) xdxdydz =>
!         1/x d/dz (mu.dw/dy) xdxdydz =>
! delta (mu.dw/dx)Axy |T-B : at (i, j+1/2, k+1/2 - k-1/2)
      ssz = avg_y_h(avg_z_h(epmu_s(ijk,m),epmu_s(ijkt,m),k),&
                    avg_z_h(epmu_s(ijkn,m),epmu_s(ijktn,m),k),j)*&
               (w_s(ijpk,m)-w_s(ijk,m))*ody_n(j)*axy_v(ijk) - &
            avg_y_h(avg_z_h(epmu_s(ijkb,m),epmu_s(ijk,m),km),&
                    avg_z_h(epmu_s(ijkbn,m),epmu_s(ijkn,m),km),j)*&
               (w_s(ijpkm,m)-w_s(ijkm,m))*ody_n(j)*axy_v(ijkm)

      RETURN
      END SUBROUTINE calc_reg_tau_v_s


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: calc_cg_tau_v_s                                         C
!  Purpose: Cross terms in the gradient of stress in V_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_v_s(IJK, M, SSX, SSY, SSZ, SBV)

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

      USE fun_avg
      USE compar
      USE geometry
      USE indices
      USE functions

! for cutcell:
! wall velocity for slip conditions
      USE bc, only: bc_hw_s, bc_uw_s, bc_vw_s, bc_ww_s
      USE bc
      USE quadric
      USE cutcell
      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

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: I, J, K, IM, KM, JP
      INTEGER :: IJKE, IJKW, IJKN, IJKT, IJKB
      INTEGER :: IJKNE, IJKNW, IJKTN, IJKBN
      INTEGER :: IJPK, IJMK, IMJK, IJKM
      INTEGER :: IMJPK, IJPKM

! for cartesian grid implementation:
      DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
      LOGICAL :: U_NODE_AT_NE,U_NODE_AT_NW,U_NODE_AT_SE,U_NODE_AT_SW
      LOGICAL :: W_NODE_AT_TN,W_NODE_AT_TS,W_NODE_AT_BN,W_NODE_AT_BS
      DOUBLE PRECISION :: dudy_at_E,dudy_at_W
      DOUBLE PRECISION :: dwdy_at_T,dwdy_at_B
      DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Wi,Sx,Sy,Sz
      DOUBLE PRECISION :: MU_S_CUT,SSX_CUT,SSZ_CUT
      DOUBLE PRECISION :: UW_s,VW_s,WW_s
      INTEGER :: BCV
      INTEGER :: BCT
!---------------------------------------------------------------------//

      j = j_of(ijk)
      jp = jp1(j)
      i = i_of(ijk)
      im = im1(i)
      k = k_of(ijk)
      km = km1(k)

      ijpk = jp_of(ijk)
      ijmk = jm_of(ijk)
      ijkm = km_of(ijk)
      imjk = im_of(ijk)
      imjpk = im_of(ijpk)
      ijpkm = jp_of(ijkm)

      ijke = east_of(ijk)
      ijkn = north_of(ijk)
      ijkne = east_of(ijkn)
      ijkw = west_of(ijk)
      ijknw = north_of(ijkw)
      ijkt = top_of(ijk)
      ijktn = north_of(ijkt)
      ijkb = bottom_of(ijk)
      ijkbn = north_of(ijkb)

! Surface forces
! bulk viscosity term
      sbv = (eplambda_s(ijkn,m)*trd_s(ijkn,m)) * axz_v(ijk) - &
            (eplambda_s(ijk,m) *trd_s(ijk,m) ) * axz_v(ijmk)

! shear stress terms
      IF(.NOT.cut_v_cell_at(ijk)) THEN
         ssx = avg_y_h(avg_x_h(epmu_s(ijk,m),epmu_s(ijke,m),i),&
                       avg_x_h(epmu_s(ijkn,m),epmu_s(ijkne,m),i),j)*&
                  (u_s(ijpk,m)-u_s(ijk,m))*oneody_n_u(ijk)*ayz_v(ijk) - &
               avg_y_h(avg_x_h(epmu_s(ijkw,m),epmu_s(ijk,m),im),&
                       avg_x_h(epmu_s(ijknw,m),epmu_s(ijkn,m),im),j)*&
                  (u_s(imjpk,m)-u_s(imjk,m))*oneody_n_u(imjk)*ayz_v(imjk)
         ssy = epmu_s(ijkn,m)*(v_s(ijpk,m)-v_s(ijk,m))*&
                  oneody_n_v(ijk)*axz_v(ijk) - &
               epmu_s(ijk,m)*(v_s(ijk,m)-v_s(ijmk,m))*&
                  oneody_n_v(ijmk)*axz_v(ijmk)
         IF(do_k) THEN
           ssz = avg_y_h(avg_z_h(epmu_s(ijk,m),epmu_s(ijkt,m),k),&
                         avg_z_h(epmu_s(ijkn,m),epmu_s(ijktn,m),k),j)*&
                    (w_s(ijpk,m)-w_s(ijk,m))*oneody_n_w(ijk)*axy_v(ijk) - &
                 avg_y_h(avg_z_h(epmu_s(ijkb,m),epmu_s(ijk,m),km),&
                         avg_z_h(epmu_s(ijkbn,m),epmu_s(ijkn,m),km),j)*&
                    (w_s(ijpkm,m)-w_s(ijkm,m))*oneody_n_w(ijkm)*axy_v(ijkm)
         ELSE
           ssz = zero
         ENDIF

! cut cell modifications
!---------------------------------------------------------------------//
      ELSE
         bcv = bc_v_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_vs = .true.
               noc_vs     = .true.
               uw_s = zero
               vw_s = zero
               ww_s = zero
            CASE (cg_fsw)
               cut_tau_vs = .false.
               noc_vs     = .false.
               uw_s = zero
               vw_s = zero
               ww_s = zero
            CASE(cg_psw)
               IF(bc_hw_s(bc_v_id(ijk),m)==undefined) THEN   ! same as NSW
                  cut_tau_vs = .true.
                  noc_vs     = .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_v_id(ijk),m)==zero) THEN   ! same as FSW
                  cut_tau_vs = .false.
                  noc_vs     = .false.
                  uw_s = zero
                  vw_s = zero
                  ww_s = zero
               ELSE                              ! partial slip
                  cut_tau_vs = .false.
                  noc_vs     = .false.
               ENDIF
            CASE (none)
! equivalent of setting tau_v_s to zero at this ijk cell
               ssx = zero
               ssy = zero
               ssz = zero
               sbv = zero
               RETURN
         END SELECT

         IF(cut_tau_vs) THEN
            mu_s_cut = (vol(ijk)*epmu_s(ijk,m) + &
                        vol(ijpk)*epmu_s(ijkn,m))/&
                       (vol(ijk) + vol(ijpk))
         ELSE
            mu_s_cut = zero
         ENDIF

! SSX:
         u_node_at_ne = ((.NOT.blocked_u_cell_at(ijpk)).AND.&
                         (.NOT.wall_u_at(ijpk)))
         u_node_at_se = ((.NOT.blocked_u_cell_at(ijk)).AND.&
                         (.NOT.wall_u_at(ijk)))
         u_node_at_nw = ((.NOT.blocked_u_cell_at(imjpk)).AND.&
                         (.NOT.wall_u_at(imjpk)))
         u_node_at_sw = ((.NOT.blocked_u_cell_at(imjk)).AND.&
                         (.NOT.wall_u_at(imjk)))

         IF(u_node_at_ne.AND.u_node_at_se) THEN
            ui = half * (u_s(ijpk,m) + u_s(ijk,m))
            xi = half * (x_u(ijpk) + x_u(ijk))
            yi = half * (y_u(ijpk) + y_u(ijk))
            zi = half * (z_u(ijpk) + z_u(ijk))
            sx = x_u(ijpk) - x_u(ijk)
            sy = y_u(ijpk) - y_u(ijk)
            sz = z_u(ijpk) - z_u(ijk)
            CALL get_del_h(ijk, 'V_MOMENTUM', xi, yi, zi, del_h, &
               nx, ny, nz)
            dudy_at_e = (u_s(ijpk,m) - u_s(ijk,m)) * oneody_n_u(ijk)
            IF(noc_vs) dudy_at_e = dudy_at_e - ((ui-uw_s) * &
                          oneody_n_u(ijk)/del_h*(sx*nx+sz*nz))
         ELSE
            dudy_at_e =  zero
         ENDIF

         IF(u_node_at_nw.AND.u_node_at_sw) THEN
            ui = half * (u_s(imjpk,m) + u_s(imjk,m))
            xi = half * (x_u(imjpk) + x_u(imjk))
            yi = half * (y_u(imjpk) + y_u(imjk))
            zi = half * (z_u(imjpk) + z_u(imjk))
            sx = x_u(imjpk) - x_u(imjk)
            sy = y_u(imjpk) - y_u(imjk)
            sz = z_u(imjpk) - z_u(imjk)
            CALL get_del_h(ijk, 'V_MOMENTUM', xi, yi, zi, del_h, &
               nx, ny, nz)
            dudy_at_w =  (u_s(imjpk,m)-u_s(imjk,m))*oneody_n_u(imjk)
            IF(noc_vs) dudy_at_w = dudy_at_w - ((ui-uw_s) * &
                          oneody_n_u(imjk)/del_h*(sx*nx+sz*nz))
         ELSE
            dudy_at_w =  zero
         ENDIF

         IF(u_node_at_se) THEN
            CALL get_del_h(ijk, 'V_MOMENTUM', x_u(ijk), y_u(ijk), &
               z_u(ijk), del_h, nx, ny, nz)
            ssx_cut = -mu_s_cut * (u_s(ijk,m) - uw_s) / &
                       del_h * (ny*nx) * area_v_cut(ijk)
         ELSE
            ssx_cut =  zero
         ENDIF

         ssx = avg_y_h(avg_x_h(epmu_s(ijk,m),epmu_s(ijke,m),i),&
                       avg_x_h(epmu_s(ijkn,m),epmu_s(ijkne,m),i),j)*&
                  dudy_at_e*ayz_v(ijk) - &
               avg_y_h(avg_x_h(epmu_s(ijkw,m),epmu_s(ijk,m),im),&
                       avg_x_h(epmu_s(ijknw,m),epmu_s(ijkn,m),im),j)*&
                  dudy_at_w*ayz_v(imjk) + ssx_cut

! SSY:
         CALL get_del_h(ijk,'V_MOMENTUM', x_v(ijk), y_v(ijk),&
            z_v(ijk), del_h, nx, ny, nz)
         ssy = epmu_s(ijkn,m)*(v_s(ijpk,m)-v_s(ijk,m))*&
                  oneody_n_v(ijk)*axz_v(ijk) - &
               epmu_s(ijk,m)*(v_s(ijk,m)-v_s(ijmk,m))*&
                  oneody_n_v(ijmk)*axz_v(ijmk) -&
               mu_s_cut * (v_s(ijk,m) - vw_s)/del_h * &
                  (ny**2) * area_v_cut(ijk)

! SSZ:
         IF(do_k) THEN
            w_node_at_tn = ((.NOT.blocked_w_cell_at(ijpk)).AND.&
                            (.NOT.wall_w_at(ijpk)))
            w_node_at_ts = ((.NOT.blocked_w_cell_at(ijk)).AND.&
                            (.NOT.wall_w_at(ijk)))
            w_node_at_bn = ((.NOT.blocked_w_cell_at(ijpkm)).AND.&
                            (.NOT.wall_w_at(ijpkm)))
            w_node_at_bs = ((.NOT.blocked_w_cell_at(ijkm)).AND.&
                            (.NOT.wall_w_at(ijkm)))

            IF(w_node_at_tn.AND.w_node_at_ts) THEN
               wi = half * (w_s(ijpk,m) + w_s(ijk,m))
               xi = half * (x_w(ijpk) + x_w(ijk))
               yi = half * (y_w(ijpk) + y_w(ijk))
               zi = half * (z_w(ijpk) + z_w(ijk))
               sx = x_w(ijpk) - x_w(ijk)
               sy = y_w(ijpk) - y_w(ijk)
               sz = z_w(ijpk) - z_w(ijk)
               CALL get_del_h(ijk, 'V_MOMENTUM', xi, yi, zi, del_h, &
                  nx, ny, nz)
               dwdy_at_t = (w_s(ijpk,m) - w_s(ijk,m)) * oneody_n_w(ijk)
               IF(noc_vs) dwdy_at_t = dwdy_at_t - ((wi-ww_s) * &
                             oneody_n_w(ijk)/del_h*(sx*nx+sz*nz))
            ELSE
               dwdy_at_t =  zero
            ENDIF

            IF(w_node_at_bn.AND.w_node_at_bs) THEN
               wi = half * (w_s(ijpkm,m) + w_s(ijkm,m))
               xi = half * (x_w(ijpkm) + x_w(ijkm))
               yi = half * (y_w(ijpkm) + y_w(ijkm))
               zi = half * (z_w(ijpkm) + z_w(ijkm))
               sx = x_w(ijpkm) - x_w(ijkm)
               sy = y_w(ijpkm) - y_w(ijkm)
               sz = z_w(ijpkm) - z_w(ijkm)
               CALL get_del_h(ijk, 'V_MOMENTUM', xi, yi, zi, del_h, &
                  nx, ny, nz)
               dwdy_at_b = (w_s(ijpkm,m)-w_s(ijkm,m)) * oneody_n_w(ijkm)
               IF(noc_vs) dwdy_at_b = dwdy_at_b - ((wi-ww_s) * &
                             oneody_n_w(ijkm)/del_h*(sx*nx+sz*nz))
            ELSE
               dwdy_at_b =  zero
            ENDIF

            IF(w_node_at_ts) THEN
               CALL get_del_h(ijk, 'V_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 * (ny*nz) * area_v_cut(ijk)
            ELSE
               ssz_cut =  zero
            ENDIF

            ssz = avg_y_h(avg_z_h(epmu_s(ijk,m),epmu_s(ijkt,m),k),&
                          avg_z_h(epmu_s(ijkn,m),epmu_s(ijktn,m),k),j)*&
                     dwdy_at_t*axy_v(ijk) - &
                  avg_y_h(avg_z_h(epmu_s(ijkb,m),epmu_s(ijk,m),km),&
                          avg_z_h(epmu_s(ijkbn,m),epmu_s(ijkn,m),km),j)*&
                     dwdy_at_b*axy_v(ijkm) + ssz_cut
         ELSE
            ssz = zero
         ENDIF  ! end if do_k

      ENDIF  ! end if/else cut_cell

      RETURN
      END SUBROUTINE calc_cg_tau_v_s