tau_w_s.f

Go to the documentation of this file.

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

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

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

! Dummy arguments
!---------------------------------------------------------------------//
! TAU_W_s
      DOUBLE PRECISION, INTENT(OUT) :: lTAU_W_s(dimension_3, dimension_m)

! Local Variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: IJK, K, IJKT
! Phase index
      INTEGER :: M, L
! Average volume fraction
      DOUBLE PRECISION :: EPSA, EPStmp
! Source terms (Surface)
      DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
! Cylindrical source terms (Volumetric)
      DOUBLE PRECISION :: Vxz
!---------------------------------------------------------------------//

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

!!$omp  parallel do &
!!$omp  private(IJK, K, IJKT,                                         &
!!$omp&         EPSA, EPStmp, SSX, SSY, SSZ, SBV, VXZ)                &
!!$omp& schedule(static)

         DO ijk = ijkstart3, ijkend3

! Skip walls where some values are undefined.
            IF(wall_at(ijk)) cycle

            k = k_of(ijk)
            ijkt = top_of(ijk)

            IF (kt_type_enum == ghd_2007) THEN
               epstmp = zero
               DO l = 1, smax
                  epstmp = epstmp + avg_z(ep_s(ijk,l),ep_s(ijkt,l),k)
               ENDDO
               epsa = epstmp
            ELSE
               epsa = avg_z(ep_s(ijk,m),ep_s(ijkt,m),k)
            ENDIF

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

               IF((.NOT.cartesian_grid).OR.(cg_safe_mode(5)==1)) THEN
                   CALL calc_reg_tau_w_s(ijk, m, ssx, ssy, ssz, sbv, vxz)
               ELSE
                   vxz = zero
                   CALL calc_cg_tau_w_s(ijk, m, ssx, ssy, ssz, sbv)
               ENDIF

! Add the terms
               ltau_w_s(ijk,m) = sbv + ssx + ssy + ssz + vxz*vol_w(ijk)
            ELSE
               ltau_w_s(ijk,m) = zero
            ENDIF   ! end if/else .not.sip_at_t .and. epsa>dil_ep_s

         ENDDO   ! end do ijk
!!$omp end parallel do

      ENDDO   ! end do m=1,mmax

      call send_recv(ltau_w_s,2)
      RETURN
      END SUBROUTINE calc_tau_w_s


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

! Modules
!---------------------------------------------------------------------//
      USE param1, only: zero, half, one, 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 functions
      USE fun_avg
      USE compar
      USE geometry
      USE indices
      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) :: VXZ

! Local Variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: I, J, K, IM, JM, KP
      INTEGER :: IJKP, IMJK, IJMK, IJKM
      INTEGER :: IJKE, IJKW, IJKN, IJKS, IJKT
      INTEGER :: IJKNT, IJKST, IJKTE, IJKTW
      INTEGER :: IJMKP, IMJKP
! Average velocity gradients
      DOUBLE PRECISION :: duodz
!---------------------------------------------------------------------//

      k = k_of(ijk)
      j = j_of(ijk)
      i = i_of(ijk)
      im = im1(i)
      jm = jm1(j)
      kp = kp1(k)

      ijkp = kp_of(ijk)
      imjk = im_of(ijk)
      ijmk = jm_of(ijk)
      ijkm = km_of(ijk)
      ijmkp = jm_of(ijkp)
      imjkp = kp_of(imjk)

      ijkt = top_of(ijk)
      ijkn = north_of(ijk)
      ijknt = top_of(ijkn)
      ijks = south_of(ijk)
      ijkst = top_of(ijks)
      ijke = east_of(ijk)
      ijkte = east_of(ijkt)
      ijkw = west_of(ijk)
      ijktw = west_of(ijkt)

! Surface forces

! bulk viscosity term
! part of 1/x d/dz (tau_zz) xdxdydz =>
!         1/x d/dz (lambda.trcD) xdxdydz=>
! delta (lambda.trcD)Ap |T-B : at (i, j, k+1 - k-1)
      sbv = (eplambda_s(ijkt,m)*trd_s(ijkt,m)-&
             eplambda_s(ijk,m)*trd_s(ijk,m))*axy(ijk)

! shear stress terms
! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
! part of (tau_xz/x + 1/x d/dx (x tau_xz) ) xdxdydz =>
!         1/x d/dx(mu.du/dz) xdxdydz =>
! delta (mu/x du/dz)Ayz |E-W : at (i+1/2-i-1/2, j, k+1/2)
      ssx = avg_z_h(avg_x_h(epmu_s(ijk,m),epmu_s(ijke,m),i),&
                    avg_x_h(epmu_s(ijkt,m),epmu_s(ijkte,m),i),k)*&
               (u_s(ijkp,m)-u_s(ijk,m))*ox_e(i)*odz_t(k)*ayz_w(ijk) - &
            avg_z_h(avg_x_h(epmu_s(ijkw,m),epmu_s(ijk,m),im),&
                    avg_x_h(epmu_s(ijktw,m),epmu_s(ijkt,m),im),k)*&
               (u_s(imjkp,m)-u_s(imjk,m))*odz_t(k)*dy(j)*(half*(dz(k)+dz(kp)))
! Same as oX_E(IM)*AYZ_W(IMJK), but avoids singularity

! part of d/dy (tau_zy) xdxdydz =>
!         d/dy (mu/x dv/dz) xdxdydz =>
! delta (mu/x dv/dz)Axz |N-S : at (i, j+1/2 - j-1/2, k+1/2)
      ssy = avg_z_h(avg_y_h(epmu_s(ijk,m),epmu_s(ijkn,m),j),&
                    avg_y_h(epmu_s(ijkt,m),epmu_s(ijknt,m),j),k)*&
               (v_s(ijkp,m)-v_s(ijk,m))*ox(i)*odz_t(k)*axz_w(ijk) - &
            avg_z_h(avg_y_h(epmu_s(ijks,m),epmu_s(ijk,m),jm),&
                    avg_y_h(epmu_s(ijkst,m),epmu_s(ijkt,m),jm),k)*&
               (v_s(ijmkp,m)-v_s(ijmk,m))*ox(i)*odz_t(k)*axz_w(ijmk)

! part of 1/x d/dz (tau_zz) xdxdydz =>
!         1/x d/dz (mu/x dw/dz) xdxdydz =>
! delta (mu/x dw/dz)Axy |T-B : at (i, j, k+1 - k-1)
      ssz = epmu_s(ijkt,m)*(w_s(ijkp,m)-w_s(ijk,m))*&
               ox(i)*odz(kp)*axy_w(ijk) - &
            epmu_s(ijk,m)*(w_s(ijk,m)-w_s(ijkm,m))*&
               ox(i)*odz(k)* axy_w(ijkm)

! Special terms for cylindrical coordinates
      vxz = zero
      IF (cylindrical) THEN

! part of 1/x d/dz (tau_zz) xdxdydz =>
!         1/x d/dz (2.mu/x u) xdxdydz =>
! delta (2.mu/x u)Axy |T-B : at (i, j, k+1 - k-1)
         ssz = ssz + &
            epmu_s(ijkt,m)*(u_s(ijkp,m)+u_s(imjkp,m))*ox(i)*axy_w(ijk) - &
            epmu_s(ijk,m)*(u_s(ijk,m)+u_s(imjk,m))*ox(i)*axy_w(ijkm)

! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
!         1/x (mu/x du/dz) xdxdydz =>
! delta (1/x mu/x du/dz)Vp : at (i, j, k+1/2)
         IF (ox_e(im) /= undefined) THEN
            duodz = (u_s(imjkp,m)-u_s(imjk,m))*ox_e(im)*odz_t(k)
         ELSE
            duodz = zero
         ENDIF
         vxz = avg_z(epmu_s(ijk,m),epmu_s(ijkt,m),k)*ox(i)*half*&
               ((u_s(ijkp,m)-u_s(ijk,m))*ox_e(i)*odz_t(k)+duodz)
      ENDIF

      RETURN
      END SUBROUTINE calc_reg_tau_w_s


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

! Modules
!---------------------------------------------------------------------//
      USE param1, only: zero, half, one, 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 functions
      USE fun_avg
      USE compar
      USE geometry
      USE indices
      USE sendrecv

! 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, JM, KP
      INTEGER :: IJKP, IMJK, IJMK, IJKM
      INTEGER :: IJKE, IJKW, IJKN, IJKS, IJKT
      INTEGER :: IJKNT, IJKST, IJKTE, IJKTW
      INTEGER :: IJMKP, IMJKP

! for cartesian grid implementation:
      DOUBLE PRECISION :: DEL_H, Nx, Ny, Nz
      LOGICAL :: U_NODE_AT_ET, U_NODE_AT_EB, U_NODE_AT_WT, U_NODE_AT_WB
      LOGICAL :: V_NODE_AT_NT, V_NODE_AT_NB, V_NODE_AT_ST, V_NODE_AT_SB
      DOUBLE PRECISION :: dudz_at_E, dudz_at_W
      DOUBLE PRECISION :: dvdz_at_N, dvdz_at_S
      DOUBLE PRECISION :: MU_S_CUT, SSX_CUT, SSY_CUT
      DOUBLE PRECISION :: Xi, Yi, Zi, Ui, Vi, Sx, Sy, Sz
      DOUBLE PRECISION :: UW_s, VW_s, WW_s
      INTEGER :: BCV
      INTEGER :: BCT
!---------------------------------------------------------------------//

      k = k_of(ijk)
      j = j_of(ijk)
      i = i_of(ijk)
      im = im1(i)
      jm = jm1(j)
      kp = kp1(k)

      ijkp = kp_of(ijk)
      imjk = im_of(ijk)
      ijmk = jm_of(ijk)
      ijkm = km_of(ijk)
      ijmkp = jm_of(ijkp)
      imjkp = kp_of(imjk)

      ijkt = top_of(ijk)
      ijkn = north_of(ijk)
      ijknt = top_of(ijkn)
      ijks = south_of(ijk)
      ijkst = top_of(ijks)
      ijke = east_of(ijk)
      ijkte = east_of(ijkt)
      ijkw = west_of(ijk)
      ijktw = west_of(ijkt)

! Surface forces
! bulk viscosity term
      sbv =  (eplambda_s(ijkt,m)*trd_s(ijkt,m)) * axy_w(ijk) - &
             (eplambda_s(ijk,m) *trd_s(ijk,m) ) * axy_w(ijkm)

! shear stress terms
      IF(.NOT.cut_w_cell_at(ijk)) THEN
         ssx = avg_z_h(avg_x_h(epmu_s(ijk,m),epmu_s(ijke,m),i),&
                       avg_x_h(epmu_s(ijkt,m),epmu_s(ijkte,m),i),k)*&
                  (u_s(ijkp,m)-u_s(ijk,m))*oneodz_t_u(ijk)*ayz_w(ijk) - &
               avg_z_h(avg_x_h(epmu_s(ijkw,m),epmu_s(ijk,m),im),&
                       avg_x_h(epmu_s(ijktw,m),epmu_s(ijkt,m),im),k)*&
                  (u_s(imjkp,m)-u_s(imjk,m))*oneodz_t_u(imjk)*ayz_w(imjk)
         ssy = avg_z_h(avg_y_h(epmu_s(ijk,m),epmu_s(ijkn,m),j),&
                       avg_y_h(epmu_s(ijkt,m),epmu_s(ijknt,m),j),k)*&
                  (v_s(ijkp,m)-v_s(ijk,m))*oneodz_t_v(ijk)*axz_w(ijk) - &
               avg_z_h(avg_y_h(epmu_s(ijks,m),epmu_s(ijk,m),jm),&
                       avg_y_h(epmu_s(ijkst,m),epmu_s(ijkt,m),jm),k)*&
                  (v_s(ijmkp,m)-v_s(ijmk,m))*oneodz_t_v(ijmk)*axz_w(ijmk)
         ssz = epmu_s(ijkt,m)*(w_s(ijkp,m)-w_s(ijk,m))*&
                  oneodz_t_w(ijk)*axy_w(ijk) - &
               epmu_s(ijk,m)*(w_s(ijk,m)-w_s(ijkm,m))*&
                  oneodz_t_w(ijkm)*axy_w(ijkm)

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

         bcv = bc_w_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_ws = .true.
               noc_ws     = .true.
               uw_s = zero
               vw_s = zero
               ww_s = zero
            CASE (cg_fsw)
               cut_tau_ws = .false.
               noc_ws     = .false.
               uw_s = zero
               vw_s = zero
               ww_s = zero
            CASE(cg_psw)
               IF(bc_hw_s(bc_w_id(ijk),m)==undefined) THEN   ! same as NSW
                  cut_tau_ws = .true.
                  noc_ws     = .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_w_id(ijk),m)==zero) THEN   ! same as FSW
                  cut_tau_ws = .false.
                  noc_ws     = .false.
                  uw_s = zero
                  vw_s = zero
                  ww_s = zero
               ELSE                              ! partial slip
                  cut_tau_ws = .false.
                  noc_ws     = .false.
               ENDIF
            CASE (none)
! equivalent of setting tau_w_s to zero at this ijk cell
               ssx = zero
               ssy = zero
               ssz = zero
               sbv = zero
               RETURN
         END SELECT

         IF(cut_tau_ws) THEN
            mu_s_cut = (vol(ijk)*epmu_s(ijk,m) + &
                        vol(ijkp)*epmu_s(ijkt,m))/(vol(ijk) + vol(ijkp))
         ELSE
            mu_s_cut = zero
         ENDIF

! SSX:
         u_node_at_et = ((.NOT.blocked_u_cell_at(ijkp)).AND.&
                         (.NOT.wall_u_at(ijkp)))
         u_node_at_eb = ((.NOT.blocked_u_cell_at(ijk)).AND.&
                         (.NOT.wall_u_at(ijk)))
         u_node_at_wt = ((.NOT.blocked_u_cell_at(imjkp)).AND.&
                         (.NOT.wall_u_at(imjkp)))
         u_node_at_wb = ((.NOT.blocked_u_cell_at(imjk)).AND.&
                         (.NOT.wall_u_at(imjk)))

         IF(u_node_at_et.AND.u_node_at_eb) THEN
            ui = half * (u_s(ijkp,m) + u_s(ijk,m))
            xi = half * (x_u(ijkp) + x_u(ijk))
            yi = half * (y_u(ijkp) + y_u(ijk))
            zi = half * (z_u(ijkp) + z_u(ijk))
            sx = x_u(ijkp) - x_u(ijk)
            sy = y_u(ijkp) - y_u(ijk)
            sz = z_u(ijkp) - z_u(ijk)
            CALL get_del_h(ijk, 'W_MOMENTUM', xi, yi, zi, del_h, &
               nx, ny, nz)
            dudz_at_e = (u_s(ijkp,m) - u_s(ijk,m)) * oneodz_t_u(ijk)
            IF(noc_ws) dudz_at_e = dudz_at_e - ((ui-uw_s) * &
                          oneodz_t_u(ijk)/del_h*(sx*nx+sy*ny))
         ELSE
            dudz_at_e =  zero
         ENDIF

         IF(u_node_at_wt.AND.u_node_at_wb) THEN
            ui = half * (u_s(imjkp,m) + u_s(imjk,m))
            xi = half * (x_u(imjkp) + x_u(imjk))
            yi = half * (y_u(imjkp) + y_u(imjk))
            zi = half * (z_u(imjkp) + z_u(imjk))
            sx = x_u(imjkp) - x_u(imjk)
            sy = y_u(imjkp) - y_u(imjk)
            sz = z_u(imjkp) - z_u(imjk)
            CALL get_del_h(ijk, 'W_MOMENTUM', xi, yi, zi, del_h, &
               nx, ny, nz)
            dudz_at_w = (u_s(imjkp,m) - u_s(imjk,m)) * oneodz_t_u(imjk)
            IF(noc_ws) dudz_at_w = dudz_at_w - ((ui-uw_s) * &
                          oneodz_t_u(imjk)/del_h*(sx*nx+sy*ny))
         ELSE
            dudz_at_w =  zero
         ENDIF

         IF(u_node_at_eb) THEN
            CALL get_del_h(ijk, 'W_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 * &
                      (nz*nx) * area_w_cut(ijk)
         ELSE
            ssx_cut =  zero
         ENDIF

         ssx = avg_z_h(avg_x_h(epmu_s(ijk,m),epmu_s(ijke,m),i),&
                       avg_x_h(epmu_s(ijkt,m),epmu_s(ijkte,m),i),k)*&
                  dudz_at_e*ayz_w(ijk) - &
               avg_z_h(avg_x_h(epmu_s(ijkw,m),epmu_s(ijk,m),im),&
                       avg_x_h(epmu_s(ijktw,m),epmu_s(ijkt,m),im),k)*&
                  dudz_at_w*ayz_w(imjk) + ssx_cut

! SSY:
         v_node_at_nt = ((.NOT.blocked_v_cell_at(ijkp)).AND.&
                         (.NOT.wall_v_at(ijkp)))
         v_node_at_nb = ((.NOT.blocked_v_cell_at(ijk)).AND.&
                         (.NOT.wall_v_at(ijk)))
         v_node_at_st = ((.NOT.blocked_v_cell_at(ijmkp)).AND.&
                         (.NOT.wall_v_at(ijmkp)))
         v_node_at_sb = ((.NOT.blocked_v_cell_at(ijmk)).AND.&
                         (.NOT.wall_v_at(ijmk)))

         IF(v_node_at_nt.AND.v_node_at_nb) THEN
            vi = half * (v_s(ijkp,m) + v_s(ijk,m))
            xi = half * (x_v(ijkp) + x_v(ijk))
            yi = half * (y_v(ijkp) + y_v(ijk))
            zi = half * (z_v(ijkp) + z_v(ijk))
            sx = x_v(ijkp) - x_v(ijk)
            sy = y_v(ijkp) - y_v(ijk)
            sz = z_v(ijkp) - z_v(ijk)
            CALL get_del_h(ijk, 'W_MOMENTUM', xi, yi, zi, del_h, &
               nx, ny, nz)
            dvdz_at_n = (v_s(ijkp,m) - v_s(ijk,m)) * oneodz_t_v(ijk)
            IF(noc_ws) dvdz_at_n = dvdz_at_n - ((vi-vw_s) * &
                          oneodz_t_v(ijk)/del_h*(sx*nx+sy*ny))
         ELSE
            dvdz_at_n =  zero
         ENDIF

         IF(v_node_at_st.AND.v_node_at_sb) THEN
            vi = half * (v_s(ijmkp,m) + v_s(ijmk,m))
            xi = half * (x_v(ijmkp) + x_v(ijmk))
            yi = half * (y_v(ijmkp) + y_v(ijmk))
            zi = half * (z_v(ijmkp) + z_v(ijmk))
            sx = x_v(ijmkp) - x_v(ijmk)
            sy = y_v(ijmkp) - y_v(ijmk)
            sz = z_v(ijmkp) - z_v(ijmk)
            CALL get_del_h(ijk, 'W_MOMENTUM', xi, yi, zi, del_h, &
               nx, ny, nz)
            dvdz_at_s = (v_s(ijmkp,m) - v_s(ijmk,m)) * oneodz_t_v(ijmk)
            IF(noc_ws) dvdz_at_s = dvdz_at_s - ((vi-vw_s) * &
                          oneodz_t_v(ijmk)/del_h*(sx*nx+sy*ny))
         ELSE
            dvdz_at_s =  zero
         ENDIF

         IF(v_node_at_nb) THEN
            CALL get_del_h(ijk, 'W_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 * &
                      (nz*ny) * area_w_cut(ijk)
         ELSE
            ssy_cut =  zero
         ENDIF

         ssy = avg_z_h(avg_y_h(epmu_s(ijk,m),epmu_s(ijkn,m),j),&
                       avg_y_h(epmu_s(ijkt,m),epmu_s(ijknt,m),j),k)*&
                  dvdz_at_n*axz_w(ijk) - &
               avg_z_h(avg_y_h(epmu_s(ijks,m),epmu_s(ijk,m),jm),&
                       avg_y_h(epmu_s(ijkst,m),epmu_s(ijkt,m),jm),k)*&
                  dvdz_at_s*axz_w(ijmk) + ssy_cut

! SSZ:
         CALL get_del_h(ijk, 'W_MOMENTUM', x_w(ijk), y_w(ijk), &
            z_w(ijk), del_h, nx, ny, nz)
         ssz = epmu_s(ijkt,m)*(w_s(ijkp,m)-w_s(ijk,m))*&
                  oneodz_t_w(ijk)*axy_w(ijk) - &
               epmu_s(ijk,m)*(w_s(ijk,m)-w_s(ijkm,m))*&
                  ox(i)*oneodz_t_w(ijkm)*axy_w(ijkm) &
               -mu_s_cut * (w_s(ijk,m) - ww_s) / del_h * &
                  (nz**2) * area_w_cut(ijk)

      ENDIF  ! end if/else cut_cell

      RETURN
      END SUBROUTINE calc_cg_tau_w_s