get_alpha.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: GET_3D_ALPHA_U_CUT_CELL                                C
!  Purpose: Calculate the correction term alpha_U                      C
!           for U-momentum cut cells,                                  C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 21-Feb-08  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Revision Number #                                  Date: ##-###-##  C
!  Author: #                                                           C
!  Purpose: #                                                          C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
  SUBROUTINE get_3d_alpha_u_cut_cell

      USE bc
      USE compar, ONLY: iend1, jend1, kend1, ijkstart3, ijkend3, mype, pe_io
      USE cutcell
      USE exit, only: mfix_exit
      USE functions, ONLY: funijk
      USE geometry, ONLY: do_k, no_k, ayz, ayz_u, flag_e
      USE indices, ONLY: i_of, j_of, k_of
      USE param1, ONLY: half, one, undefined, zero

      IMPLICIT NONE
      DOUBLE PRECISION:: Xe,Ye,Ze,Xn,Yn,Zn,Xt,Yt,Zt
      INTEGER :: I,J,K,IM,JM,KM,IP,JP,KP,IJK,IMJK,IJMK,IPJK,IJPK,IJKM,IJKP
      DOUBLE PRECISION :: Xi
      DOUBLE PRECISION :: Sx,Sy,Sz
      DOUBLE PRECISION :: Nx,Ny,Nz
      DOUBLE PRECISION :: DELH_ec , DELH_e
      DOUBLE PRECISION :: DELH_nc , DELH_n
      DOUBLE PRECISION :: DELH_tc , DELH_t
      LOGICAL :: V_NODE_AT_NE,V_NODE_AT_NW
      LOGICAL :: W_NODE_AT_TE,W_NODE_AT_TW
      INTEGER :: BCV
      INTEGER ::BCT

      IF(mype == pe_io) THEN
         WRITE(*,10)'COMPUTING INTERPOLATION FACTORS IN U-MOMENTUM CELLS...'
      ENDIF
10    FORMAT(1x,a)

      delh_u = undefined

      theta_ue  = half
      theta_ue_bar = half
      theta_u_ne  = half
      theta_u_nw  = half
      theta_u_te  = half
      theta_u_tw  = half
      alpha_ue_c  = one
      noc_u_e  = zero
      theta_un  = half
      theta_un_bar = half
      alpha_un_c  = one
      noc_u_n  = zero
      theta_ut  = half
      theta_ut_bar = half
      alpha_ut_c  = one
      noc_u_t  = zero

      a_upg_e = ayz_u
      a_upg_w = ayz_u


      DO ijk = ijkstart3, ijkend3
         IF(interior_cell_at(ijk)) THEN

            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)

            im = i - 1
            jm = j - 1
            km = k - 1

            ip = i + 1
            jp = j + 1
            kp = k + 1

            imjk = funijk(im,j,k)
            ipjk = funijk(ip,j,k)
            ijmk = funijk(i,jm,k)
            ijpk = funijk(i,jp,k)
            ijkm = funijk(i,j,km)
            ijkp = funijk(i,j,kp)

            CALL get_cell_node_coordinates(ijk,'U_MOMENTUM')

!======================================================================
!  Get Interpolation factors at East face
!======================================================================

            theta_ue(ijk) = delx_ue(ijk) / (delx_ue(ijk) + delx_uw(ipjk))
            theta_ue_bar(ijk) = one - theta_ue(ijk)

!======================================================================
!  Get Interpolation factors at North face
!======================================================================

            v_node_at_nw = ((.NOT.blocked_v_cell_at(ijk)).AND.(.NOT.wall_v_at(ijk)))
            v_node_at_ne = ((.NOT.blocked_v_cell_at(ipjk)).AND.(.NOT.wall_v_at(ipjk)))

            IF(v_node_at_nw.AND.v_node_at_ne) THEN
                theta_u_ne(ijk) = (x_u_nc(ijk) - x_v(ijk) ) / (x_v(ipjk) - x_v(ijk))
                theta_u_nw(ijk) = one - theta_u_ne(ijk)
            ELSE IF (v_node_at_ne.AND.(.NOT.v_node_at_nw)) THEN
               IF(no_k) THEN
                  xi = xn_u_int(ijk)
               ELSE
                  xi = half * (xn_u_int(ijk) + xn_u_int(ijkm))
               ENDIF
               theta_u_ne(ijk) = (x_u_nc(ijk) - xi) / (x_v(ipjk) - xi)
               theta_u_nw(ijk) = one - theta_u_ne(ijk)
            ELSE IF ((.NOT.v_node_at_ne).AND.v_node_at_nw) THEN
               IF(no_k) THEN
                  xi = xn_u_int(ijk)
               ELSE
                  xi = half * (xn_u_int(ijk) + xn_u_int(ijkm))
               ENDIF
               theta_u_ne(ijk) = (x_u_nc(ijk) - x_v(ijk) ) / (xi - x_v(ijk))
               theta_u_nw(ijk) = one - theta_u_ne(ijk)
            ELSE
               theta_u_ne(ijk) = zero
               theta_u_nw(ijk) = zero
            ENDIF


            IF ((theta_u_ne(ijk)>=one).OR.(theta_u_ne(ijk)<=zero)) THEN
               theta_u_ne(ijk) = half
               theta_u_nw(ijk) = half
            ENDIF

!======================================================================
!  Get Interpolation factors at Top face
!======================================================================

            IF(do_k) THEN

               w_node_at_tw = ((.NOT.blocked_w_cell_at(ijk)).AND.(.NOT.wall_w_at(ijk)))
               w_node_at_te = ((.NOT.blocked_w_cell_at(ipjk)).AND.(.NOT.wall_w_at(ipjk)))

               IF(w_node_at_tw.AND.w_node_at_te) THEN
                  theta_u_te(ijk) = (x_u_tc(ijk) - x_w(ijk) ) / (x_w(ipjk) - x_w(ijk))
                  theta_u_tw(ijk) = one - theta_u_te(ijk)

               ELSE IF (w_node_at_te.AND.(.NOT.w_node_at_tw)) THEN
                  xi = half * (xn_u_int(ijk) + xn_u_int(ijmk))
                  theta_u_te(ijk) = (x_u_tc(ijk) - xi) / (x_w(ipjk) - xi)
                  theta_u_tw(ijk) = one - theta_u_te(ijk)
               ELSE IF ((.NOT.w_node_at_te).AND.w_node_at_tw) THEN
                  xi = half * (xn_u_int(ijk) + xn_u_int(ijmk))
                  theta_u_te(ijk) = (x_u_tc(ijk) - x_w(ijk) ) / (xi - x_w(ijk))
                  theta_u_tw(ijk) = one - theta_u_te(ijk)
               ELSE
                  theta_u_te(ijk) = zero
                  theta_u_tw(ijk) = zero
               ENDIF

               IF ((theta_u_te(ijk)>=one).OR.(theta_u_te(ijk)<=zero)) THEN
                  theta_u_te(ijk) = half
                  theta_u_tw(ijk) = half
               ENDIF
            ENDIF

            bcv = bc_u_id(ijk)

            IF(bcv>0) THEN
               bct = bc_type_enum(bcv)
            ELSE
               bct = blank
            ENDIF

            IF(bct ==cg_nsw.OR.bct ==cg_psw) THEN

               CALL get_del_h(ijk,'U_MOMENTUM',x_u(ijk),y_u(ijk),z_u(ijk),delh_u(ijk),nx,ny,nz)

               IF(delh_u(ijk)<zero.AND.(.NOT.wall_u_at(ijk))) THEN
                  WRITE(*,*) 'NEGATIVE DELH-U AT XYZ=', x_u(ijk),y_u(ijk),z_u(ijk)
                  WRITE(*,*) 'DELH_U=', delh_u(ijk)
                  WRITE(*,*) 'AYZ=', ayz(ijk)
                  WRITE(*,*) 'MFIX WILL EXIT NOW.'
                  CALL mfix_exit(mype)
               ENDIF

!======================================================================
!  Get Alpha Correction factors at East face
!======================================================================

! Location of the interpolated velocity along east face
! Xe has not moved based on definition of Theta_Ue

               xe = x_node(8)
               ye = theta_ue_bar(ijk) * y_u(ijk) + theta_ue(ijk) * y_u(ipjk)
               ze = theta_ue_bar(ijk) * z_u(ijk) + theta_ue(ijk) * z_u(ipjk)

               CALL get_del_h(ijk,'U_MOMENTUM',x_u_ec(ijk),y_u_ec(ijk),z_u_ec(ijk),delh_ec,nx,ny,nz)
               CALL get_del_h(ijk,'U_MOMENTUM',xe,ye,ze,delh_e,nx,ny,nz)

               IF((delh_ec == undefined).OR.(delh_e == undefined)) THEN
                  alpha_ue_c(ijk) = one
               ELSE
                  alpha_ue_c(ijk) = dmin1(alpha_max , delh_ec / delh_e )
               ENDIF


               IF(blocked_u_cell_at(ipjk)) alpha_ue_c(ijk) = zero
               IF(wall_u_at(ijk).AND.wall_u_at(ipjk)) alpha_ue_c(ijk) = zero
               IF(i == iend1) alpha_ue_c(ijk) = one

               IF(alpha_ue_c(ijk)<zero) THEN
                 WRITE(*,*) 'NEGATIVE ALPHA_Ue_c at IJK=',ijk
                 WRITE(*,*) 'MFIX WILL EXIT NOW.'
                 CALL mfix_exit(mype)
               ENDIF

!======================================================================
!  Get Non-Ortogonality correction factor at East face
!======================================================================

               sx = x_u(ipjk) - x_u(ijk)
               sy = y_u(ipjk) - y_u(ijk)
               sz = z_u(ipjk) - z_u(ijk)

               noc_u_e(ijk) = (sy * ny + sz * nz)/(sx * delh_e)

               IF(blocked_u_cell_at(ipjk)) noc_u_e(ijk) = zero
               IF(wall_u_at(ijk).AND.wall_u_at(ipjk)) noc_u_e(ijk) = zero
               IF(i == iend1) noc_u_e(ijk) = zero

!======================================================================
!  Get Alpha Correction factors at North face
!======================================================================

               theta_un(ijk) = dely_un(ijk) / (dely_un(ijk) + dely_us(ijpk))
               theta_un_bar(ijk) = one - theta_un(ijk)

               IF ((theta_un(ijk)>=one).OR.(theta_un(ijk)<=zero)) THEN
                  theta_un(ijk) = half
                  theta_un_bar(ijk) = half
               ENDIF


! Location of the interpolated velocity along North face
! Yn has not moved based on definition of Theta_Un


               xn = theta_un_bar(ijk) * x_u(ijk) + theta_un(ijk) * x_u(ijpk)
               yn = y_node(8)
               zn = theta_un_bar(ijk) * z_u(ijk) + theta_un(ijk) * z_u(ijpk)

               CALL get_del_h(ijk,'U_MOMENTUM',x_u_nc(ijk),y_u_nc(ijk),z_u_nc(ijk),delh_nc,nx,ny,nz)
               CALL get_del_h(ijk,'U_MOMENTUM',xn,yn,zn,delh_n,nx,ny,nz)

               IF((delh_nc == undefined).OR.(delh_n == undefined)) THEN
                  alpha_un_c(ijk) = one
               ELSE
                  alpha_un_c(ijk) = dmin1(alpha_max , delh_nc / delh_n )
               ENDIF


               IF(blocked_u_cell_at(ijpk)) alpha_un_c(ijk) = zero
               IF(wall_u_at(ijk).AND.wall_u_at(ijpk)) alpha_un_c(ijk) = zero
               IF(j == jend1) alpha_un_c(ijk) = one

               IF(alpha_un_c(ijk)<zero) alpha_un_c(ijk) = one

!======================================================================
!  Get Non-Ortogonality correction factor at North face
!======================================================================

               sx = x_u(ijpk) - x_u(ijk)
               sy = y_u(ijpk) - y_u(ijk)
               sz = z_u(ijpk) - z_u(ijk)

               noc_u_n(ijk) = (sx * nx + sz * nz)/(sy * delh_n)

               IF(blocked_u_cell_at(ijpk)) noc_u_n(ijk) = zero
               IF(wall_u_at(ijk).AND.wall_u_at(ijpk)) noc_u_n(ijk) = zero
               IF(j == jend1) noc_u_n(ijk) = zero


               IF(do_k) THEN
!======================================================================
!  Get Alpha Correction factors at Top face
!======================================================================

                  theta_ut(ijk) = delz_ut(ijk) / (delz_ut(ijk) + delz_ub(ijkp))
                  theta_ut_bar(ijk) = one - theta_ut(ijk)

                  IF ((theta_ut(ijk)>=one).OR.(theta_ut(ijk)<=zero)) THEN
                     theta_ut(ijk) = half
                     theta_ut_bar(ijk) = half
                  ENDIF


! Location of the interpolated velocity along Top face
! Zt has not moved based on definition of Theta_Ut

                  xt = theta_ut_bar(ijk) * x_u(ijk) + theta_ut(ijk) * x_u(ijkp)
                  yt = theta_ut_bar(ijk) * y_u(ijk) + theta_ut(ijk) * y_u(ijkp)
                  zt = z_node(8)

                  CALL get_del_h(ijk,'U_MOMENTUM',x_u_tc(ijk),y_u_tc(ijk),z_u_tc(ijk),delh_tc,nx,ny,nz)
                  CALL get_del_h(ijk,'U_MOMENTUM',xt,yt,zt,delh_t,nx,ny,nz)

                  IF((delh_tc == undefined).OR.(delh_t == undefined)) THEN
                     alpha_ut_c(ijk) = one
                  ELSE
                     alpha_ut_c(ijk) = dmin1(alpha_max , delh_tc / delh_t )
                  ENDIF

                  IF(blocked_u_cell_at(ijkp)) alpha_ut_c(ijk) = zero
                  IF(wall_u_at(ijk).AND.wall_u_at(ijkp)) alpha_ut_c(ijk) = zero
                  IF(k == kend1) alpha_ut_c(ijk) = one

                  IF(alpha_ut_c(ijk)<zero) alpha_ut_c(ijk) = one

!======================================================================
!  Get Non-Ortogonality correction factor at Top face
!======================================================================

                  sx = x_u(ijkp) - x_u(ijk)
                  sy = y_u(ijkp) - y_u(ijk)
                  sz = z_u(ijkp) - z_u(ijk)

                  noc_u_t(ijk) = (sx * nx + sy * ny)/(sz * delh_t)


                  IF(blocked_u_cell_at(ijkp)) noc_u_t(ijk) = zero
                  IF(wall_u_at(ijk).AND.wall_u_at(ijkp)) noc_u_t(ijk) = zero
                  IF(k == kend1) noc_u_t(ijk) = zero

               ENDIF

            ENDIF



!======================================================================
!  Get Surfaces used to compute pressure gradient
!======================================================================

            IF (  cut_u_cell_at(ijk)  )  THEN

               SELECT CASE (pg_option)

                  CASE(0)                               ! do nothing
                                                        ! will be treated below

                  CASE(1)

                     a_upg_e(ijk) = dmax1(ayz_u(ijk),ayz_u(imjk))
                     a_upg_w(ijk) = a_upg_e(ijk)

                  CASE(2)

                     a_upg_e(ijk) = ayz_u(ijk)
                     a_upg_w(ijk) = ayz_u(imjk)

                  CASE DEFAULT

                     WRITE(*,*)'INVALID PRESSURE GRADIENT OPTION:',pg_option
                     WRITE(*,*)'PG_OPTION SHOULD BE SET EQUAL TO 0, 1 OR 2.'
                     WRITE(*,*)'PLEASE VERIFY MFIX.DAT FILE AND TRY AGAIN.'
                     WRITE(*,*) 'MFIX WILL EXIT NOW.'
                     CALL mfix_exit(mype)

               END SELECT


               IF (blocked_cell_at(ijk).OR.blocked_cell_at(ipjk)) THEN
                  IF(.NOT.wall_u_at(ijk)) THEN
                     wall_u_at(ijk) = .true.
                     flag_e(ijk) = 0
                     IF(print_warnings) THEN
                        WRITE(*,*) 'WARNING: ONLY ONE PRESSURE NODE DETECTED IN U-MOMENTUM CELL IJK =',ijk
                        WRITE(*,*) 'RESETTING U-MOMENTUM CELL AS U_WALL CELL.'
                     ENDIF
!                     write(*,*) 'MFiX will exit now.'
!                     CALL MFIX_EXIT(myPE)
                  ENDIF
               ENDIF

            ELSE   ! Regular cell

               a_upg_e(ijk) = ayz_u(ijk)
               a_upg_w(ijk) = ayz_u(ijk)

            ENDIF

         ENDIF
      END DO

      IF(pg_option==0) THEN
         a_upg_e = ayz
         a_upg_w = ayz
      ENDIF

      RETURN

      END SUBROUTINE get_3d_alpha_u_cut_cell

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: GET_3D_ALPHA_V_CUT_CELL                                C
!  Purpose: Calculate the correction term alpha_V                      C
!           for V-momentum cut cells,                                  C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 21-Feb-08  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Revision Number #                                  Date: ##-###-##  C
!  Author: #                                                           C
!  Purpose: #                                                          C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
  SUBROUTINE get_3d_alpha_v_cut_cell

      USE bc
      USE compar, ONLY: iend1, jend1, kend1, ijkstart3, ijkend3, mype, pe_io
      USE cutcell
      USE functions, ONLY: funijk
      USE exit, only: mfix_exit
      USE geometry, ONLY: do_k, no_k, axz, axz_v, flag_n
      USE indices, ONLY: i_of, j_of, k_of
      USE param1, ONLY: half, one, undefined, zero

      IMPLICIT NONE
      DOUBLE PRECISION:: Xe,Ye,Ze,Xn,Yn,Zn,Xt,Yt,Zt
      INTEGER :: I,J,K,IM,JM,KM,IP,JP,KP,IJK,IMJK,IJMK,IPJK,IJPK,IJKM,IJKP
      DOUBLE PRECISION :: Yi
      DOUBLE PRECISION :: Sx,Sy,Sz
      DOUBLE PRECISION :: Nx,Ny,Nz
      DOUBLE PRECISION :: DELH_ec , DELH_e
      DOUBLE PRECISION :: DELH_nc , DELH_n
      DOUBLE PRECISION :: DELH_tc , DELH_t
      LOGICAL :: U_NODE_AT_NE, U_NODE_AT_SE
      LOGICAL :: W_NODE_AT_NT, W_NODE_AT_ST
      INTEGER :: BCV
      INTEGER ::BCT

      IF(mype == pe_io) THEN
         WRITE(*,10)'COMPUTING INTERPOLATION FACTORS IN V-MOMENTUM CELLS...'
      ENDIF
10    FORMAT(1x,a)

      delh_v = undefined

      theta_v_ne = half
      theta_v_se  = half
      theta_vn  = half
      theta_vn_bar = half
      theta_v_nt  = half
      theta_v_st = half
      theta_ve  = half
      theta_ve_bar = half
      alpha_ve_c  = one
      noc_v_e  = zero
      alpha_vn_c  = one
      noc_v_n  = zero
      theta_vt  = half
      theta_vt_bar = half
      alpha_vt_c  = one
      noc_v_t  = zero

      a_vpg_n = axz_v
      a_vpg_s = axz_v


      DO ijk = ijkstart3, ijkend3
         IF(interior_cell_at(ijk)) THEN

            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)

            im = i - 1
            jm = j - 1
            km = k - 1

            ip = i + 1
            jp = j + 1
            kp = k + 1

            imjk = funijk(im,j,k)
            ipjk = funijk(ip,j,k)
            ijmk = funijk(i,jm,k)
            ijpk = funijk(i,jp,k)
            ijkm = funijk(i,j,km)
            ijkp = funijk(i,j,kp)

            CALL get_cell_node_coordinates(ijk,'V_MOMENTUM')

!======================================================================
!  Get Interpolation factors at East face
!======================================================================

            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)))

            IF(u_node_at_se.AND.u_node_at_ne) THEN
               theta_v_ne(ijk) = (y_v_ec(ijk) - y_u(ijk)    ) / (y_u(ijpk) - y_u(ijk))
               theta_v_se(ijk) = one - theta_v_ne(ijk)
            ELSE IF (u_node_at_se.AND.(.NOT.u_node_at_ne)) THEN
               IF(no_k) THEN
                  yi = ye_v_int(ijk)
               ELSE
                  yi = half * (ye_v_int(ijk) + ye_v_int(ijkm))
               ENDIF
               theta_v_ne(ijk) = (y_v_ec(ijk) - y_u(ijk)    ) / (yi - y_u(ijk))
               theta_v_se(ijk) = one - theta_v_ne(ijk)
            ELSE IF ((.NOT.u_node_at_se).AND.u_node_at_ne) THEN
               IF(no_k) THEN
                  yi = ye_v_int(ijk)
               ELSE
                  yi = half * (ye_v_int(ijk) + ye_v_int(ijkm))
               ENDIF
               theta_v_ne(ijk) = (y_v_ec(ijk) - yi    ) / (y_u(ijpk) - yi)
               theta_v_se(ijk) = one - theta_v_ne(ijk)
            ELSE
               theta_v_ne(ijk) = zero
               theta_v_se(ijk) = zero
            ENDIF

            IF ((theta_v_ne(ijk)>=one).OR.(theta_v_ne(ijk)<=zero)) THEN
               theta_v_ne(ijk) = half
               theta_v_se(ijk) = half
            ENDIF

!======================================================================
!  Get Interpolation factors at North face
!======================================================================

            theta_vn(ijk) = dely_vn(ijk) / (dely_vn(ijk) + dely_vs(ijpk))
            theta_vn_bar(ijk) = one - theta_vn(ijk)

            IF ((theta_vn(ijk)>=one).OR.(theta_vn(ijk)<=zero)) THEN
               theta_vn(ijk) = half
               theta_vn_bar(ijk) = half
            ENDIF

            IF(do_k) THEN
!======================================================================
!  Get Interpolation factors at Top face
!======================================================================

               w_node_at_nt = ((.NOT.blocked_w_cell_at(ijpk)).AND.(.NOT.wall_w_at(ijpk)))
               w_node_at_st = ((.NOT.blocked_w_cell_at(ijk)).AND.(.NOT.wall_w_at(ijk)))

               IF(w_node_at_st.AND.w_node_at_nt) THEN
                  theta_v_nt(ijk) = (y_v_tc(ijk) - y_w(ijk)    ) / (y_w(ijpk) - y_w(ijk))
                  theta_v_st(ijk) = one - theta_v_nt(ijk)
               ELSE IF (w_node_at_st.AND.(.NOT.w_node_at_nt)) THEN
                  yi = half * (ye_v_int(ijk) + ye_v_int(imjk))
                  theta_v_nt(ijk) = (y_v_tc(ijk) - y_w(ijk)    ) / (yi - y_w(ijk))
                  theta_v_st(ijk) = one - theta_v_nt(ijk)
               ELSE IF ((.NOT.w_node_at_st).AND.w_node_at_nt) THEN
                  yi = half * (ye_v_int(ijk) + ye_v_int(imjk))
                  theta_v_nt(ijk) = (y_v_tc(ijk) - yi ) / (y_w(ijpk) - yi)
                  theta_v_st(ijk) = one - theta_v_nt(ijk)
               ELSE
                  theta_v_nt(ijk) = zero
                  theta_v_st(ijk) = zero
               ENDIF

               IF ((theta_v_nt(ijk)>=one).OR.(theta_v_nt(ijk)<=zero)) THEN
                  theta_v_nt(ijk) = half
                  theta_v_st(ijk) = half
               ENDIF

            ENDIF

            bcv = bc_v_id(ijk)
            IF(bcv>0) THEN
               bct = bc_type_enum(bcv)
            ELSE
               bct = blank
            ENDIF

            IF(bct ==cg_nsw.OR.bct ==cg_psw) THEN

               CALL get_del_h(ijk,'V_MOMENTUM',x_v(ijk),y_v(ijk),z_v(ijk),delh_v(ijk),nx,ny,nz)

               IF(delh_v(ijk)<zero.AND.(.NOT.wall_v_at(ijk))) THEN
                 WRITE(*,*) 'NEGATIVE DELH-V AT XYZ=', x_v(ijk),y_v(ijk),z_v(ijk)
                 WRITE(*,*) 'DELH_V=', delh_v(ijk)
                 WRITE(*,*) 'AYZ=', axz(ijk)
                 WRITE(*,*) 'MFIX WILL EXIT NOW.'
                 CALL mfix_exit(mype)
              ENDIF

!======================================================================
!  Get Alpha Correction factors at East face
!======================================================================

               theta_ve(ijk) = delx_ve(ijk) / (delx_ve(ijk) + delx_vw(ipjk))
               theta_ve_bar(ijk) = one - theta_ve(ijk)

            IF ((theta_ve(ijk)>=one).OR.(theta_ve(ijk)<=zero)) THEN
               theta_ve(ijk) = half
               theta_ve_bar(ijk) = half
            ENDIF


! Location of the interpolated velocity along East face
! Xe has not moved based on definition of Theta_Ve
               xe = x_node(8)
               ye = theta_ve_bar(ijk) * y_v(ijk) + theta_ve(ijk) * y_v(ipjk)
               ze = theta_ve_bar(ijk) * z_v(ijk) + theta_ve(ijk) * z_v(ipjk)


               CALL get_del_h(ijk,'V_MOMENTUM',x_v_ec(ijk),y_v_ec(ijk),z_v_ec(ijk),delh_ec,nx,ny,nz)

               CALL get_del_h(ijk,'V_MOMENTUM',xe,ye,ze,delh_e,nx,ny,nz)

               IF((delh_ec == undefined).OR.(delh_e == undefined)) THEN
                  alpha_ve_c(ijk) = one
               ELSE
                  alpha_ve_c(ijk) = dmin1(alpha_max , delh_ec / delh_e )
               ENDIF


               IF(blocked_v_cell_at(ipjk)) alpha_ve_c(ijk) = zero
               IF(wall_v_at(ijk).AND.wall_v_at(ipjk)) alpha_ve_c(ijk) = zero
               IF(i == iend1) alpha_ve_c(ijk) = one

               IF(alpha_ve_c(ijk)<zero) alpha_ve_c(ijk) = one

!======================================================================
!  Get Non-Ortogonality correction factor at East face
!======================================================================

               sx = x_v(ipjk) - x_v(ijk)
               sy = y_v(ipjk) - y_v(ijk)
               sz = z_v(ipjk) - z_v(ijk)

               noc_v_e(ijk) = (sy * ny + sz * nz)/(sx * delh_e)

               IF(blocked_v_cell_at(ipjk)) noc_v_e(ijk) = zero
               IF(wall_v_at(ijk).AND.wall_v_at(ipjk)) noc_v_e(ijk) = zero
               IF(i == iend1) noc_v_e(ijk) = zero

!======================================================================
!  Get Alpha Correction factors at North face
!======================================================================

! Location of the interpolated velocity along north face
! Yn has not moved based on definition of Theta_V

               xn = theta_vn_bar(ijk) * x_v(ijk) + theta_vn(ijk) * x_v(ijpk)
               yn = y_node(8)
               zn = theta_vn_bar(ijk) * z_v(ijk) + theta_vn(ijk) * z_v(ijpk)

               CALL get_del_h(ijk,'V_MOMENTUM',x_v_nc(ijk),y_v_nc(ijk),z_v_nc(ijk),delh_nc,nx,ny,nz)

               CALL get_del_h(ijk,'V_MOMENTUM',xn,yn,zn,delh_n,nx,ny,nz)

               IF((delh_nc == undefined).OR.(delh_n == undefined)) THEN
                  alpha_vn_c(ijk) = one
               ELSE
                  alpha_vn_c(ijk) = dmin1(alpha_max , delh_nc / delh_n )
               ENDIF

               IF(blocked_v_cell_at(ijpk)) alpha_vn_c(ijk) = zero
               IF(wall_v_at(ijk).AND.wall_v_at(ijpk)) alpha_vn_c(ijk) = zero
               IF(j == jend1) alpha_vn_c(ijk) = one

               IF(alpha_vn_c(ijk)<zero) THEN
                 WRITE(*,*) 'NEGATIVE ALPHA_Vn_c at IJK=',ijk
                 WRITE(*,*) 'MFIX WILL EXIT NOW.'
                 CALL mfix_exit(mype)
               ENDIF

!======================================================================
!  Get Non-Ortogonality correction factor at North face
!======================================================================

               sx = x_v(ijpk) - x_v(ijk)
               sy = y_v(ijpk) - y_v(ijk)
               sz = z_v(ijpk) - z_v(ijk)

               noc_v_n(ijk) = (sx * nx + sz * nz)/(sy * delh_n)

               IF(blocked_v_cell_at(ijpk)) noc_v_n(ijk) = zero
               IF(wall_v_at(ijk).AND.wall_v_at(ijpk)) noc_v_n(ijk) = zero
               IF(j == jend1) noc_v_n(ijk) = zero

               IF(do_k) THEN

!======================================================================
!  Get Alpha Correction factors at Top face
!======================================================================

                  theta_vt(ijk) = delz_vt(ijk) / (delz_vt(ijk) + delz_vb(ijkp))
                  theta_vt_bar(ijk) = one - theta_vt(ijk)

                  IF ((theta_vt(ijk)>=one).OR.(theta_vt(ijk)<=zero)) THEN
                     theta_vt(ijk) = half
                     theta_vt_bar(ijk) = half
                  ENDIF

! Location of the interpolated velocity along Top face
! Zt has not moved based on definition of Theta_Ut

                  xt = theta_vt_bar(ijk) * x_v(ijk) + theta_vt(ijk) * x_v(ijkp)
                  yt = theta_vt_bar(ijk) * y_v(ijk) + theta_vt(ijk) * y_v(ijkp)
                  zt = z_node(8)

                  CALL get_del_h(ijk,'V_MOMENTUM',x_v_tc(ijk),y_v_tc(ijk),z_v_tc(ijk),delh_tc,nx,ny,nz)

                  CALL get_del_h(ijk,'V_MOMENTUM',xt,yt,zt,delh_t,nx,ny,nz)

                  IF((delh_tc == undefined).OR.(delh_t == undefined)) THEN
                     alpha_vt_c(ijk) = one
                  ELSE
                     alpha_vt_c(ijk) = dmin1(alpha_max , delh_tc / delh_t )
                  ENDIF

                  IF(blocked_v_cell_at(ijkp)) alpha_vt_c(ijk) = zero
                  IF(wall_v_at(ijk).AND.wall_v_at(ijkp)) alpha_vt_c(ijk) = zero
                  IF(k == kend1) alpha_vt_c(ijk) = one

                  IF(alpha_vt_c(ijk)<zero) alpha_vt_c(ijk) = one

!======================================================================
!  Get Non-Ortogonality correction factor at Top face
!======================================================================

                  sx = x_v(ijkp) - x_v(ijk)
                  sy = y_v(ijkp) - y_v(ijk)
                  sz = z_v(ijkp) - z_v(ijk)

                  noc_v_t(ijk) = (sx * nx + sy * ny)/(sz * delh_t)


                  IF(blocked_v_cell_at(ijkp)) noc_v_t(ijk) = zero
                  IF(wall_v_at(ijk).AND.wall_v_at(ijkp)) noc_v_t(ijk) = zero
                  IF(k == kend1) noc_v_t(ijk) = zero

               ENDIF

            ENDIF

!======================================================================
!  Get Surfaces used to compute pressure gradient
!======================================================================

            IF (  cut_v_cell_at(ijk)  )  THEN

               SELECT CASE (pg_option)

                  CASE(0)                               ! do nothing
                                                     ! will be treated below

                  CASE(1)

                     a_vpg_n(ijk) = dmax1(axz_v(ijk),axz_v(ijmk))
                     a_vpg_s(ijk) = a_vpg_n(ijk)

                  CASE(2)

                     a_vpg_n(ijk) = axz_v(ijk)
                     a_vpg_s(ijk) = axz_v(ijmk)


                  CASE DEFAULT

                     WRITE(*,*)'INVALID PRESSURE GRADIENT OPTION:',pg_option
                     WRITE(*,*)'PG_OPTION SHOULD BE SET EQUAL TO 0, 1 OR 2.'
                     WRITE(*,*)'PLEASE VERIFY MFIX.DAT FILE AND TRY AGAIN.'
                     WRITE(*,*) 'MFIX WILL EXIT NOW.'
                     CALL mfix_exit(mype)

               END SELECT


               IF (blocked_cell_at(ijk).OR.blocked_cell_at(ijpk)) THEN
                  IF(.NOT.wall_v_at(ijk)) THEN
                     wall_v_at(ijk) = .true.
                     flag_n(ijk) = 0
                     IF(print_warnings) THEN
                        WRITE(*,*) 'WARNING: ONLY ONE PRESSURE NODE DETECTED IN V-MOMENTUM CELL IJK =',ijk
                        WRITE(*,*) 'RESETTING U-MOMENTUM CELL AS V_WALL CELL.'
                     ENDIF
!                     write(*,*) 'MFiX will exit now.'
!                     CALL MFIX_EXIT(myPE)
                  ENDIF
               ENDIF

            ELSE

               a_vpg_n(ijk) = axz_v(ijk)
               a_vpg_s(ijk) = axz_v(ijk)

            ENDIF

         ENDIF
      END DO

      IF(pg_option==0) THEN
         a_vpg_n = axz
         a_vpg_s = axz
      ENDIF

      RETURN

      END SUBROUTINE get_3d_alpha_v_cut_cell

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: GET_3D_ALPHA_W_CUT_CELL                                C
!  Purpose: Calculate the correction term alpha_W                      C
!           for W-momentum cut cells,                                  C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 21-Feb-08  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Revision Number #                                  Date: ##-###-##  C
!  Author: #                                                           C
!  Purpose: #                                                          C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
  SUBROUTINE get_3d_alpha_w_cut_cell

      USE bc
      USE compar, ONLY: iend1, jend1, kend1, ijkstart3, ijkend3, mype, pe_io
      USE cutcell
      USE exit, only: mfix_exit
      USE functions, ONLY: funijk
      USE geometry, ONLY: axy, axy_w, flag_t
      USE indices, ONLY: i_of, j_of, k_of
      USE param1, ONLY: half, one, undefined, zero

      IMPLICIT NONE
      DOUBLE PRECISION:: Xe,Ye,Ze,Xn,Yn,Zn,Xt,Yt,Zt
      INTEGER :: I,J,K,IM,JM,KM,IP,JP,KP,IJK,IMJK,IJMK,IPJK,IJPK,IJKM,IJKP
      DOUBLE PRECISION :: Zi
      DOUBLE PRECISION :: Sx,Sy,Sz
      DOUBLE PRECISION :: Nx,Ny,Nz
      DOUBLE PRECISION :: DELH_ec , DELH_e
      DOUBLE PRECISION :: DELH_nc , DELH_n
      DOUBLE PRECISION :: DELH_tc , DELH_t
      LOGICAL :: U_NODE_AT_TE, U_NODE_AT_BE
      LOGICAL :: V_NODE_AT_TN, V_NODE_AT_BN
      INTEGER :: BCV
      INTEGER ::BCT

      IF(mype == pe_io) THEN
         WRITE(*,10)'COMPUTING INTERPOLATION FACTORS IN W-MOMENTUM CELLS...'
      ENDIF
10    FORMAT(1x,a)

      delh_w = undefined

      theta_w_te = half
      theta_w_be = half
      theta_w_tn = half
      theta_w_bn = half
      theta_wt = half
      theta_wt_bar = half
      theta_we = half
      theta_we_bar = half
      alpha_we_c = one
      noc_w_e = zero
      theta_wn = half
      theta_wn_bar = half
      alpha_wn_c = one
      noc_w_n = zero
      alpha_wt_c = one
      noc_w_t = zero
      a_wpg_t = axy_w
      a_wpg_b = axy_w

      DO ijk = ijkstart3, ijkend3
         IF(interior_cell_at(ijk)) THEN

            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)

            im = i - 1
            jm = j - 1
            km = k - 1

            ip = i + 1
            jp = j + 1
            kp = k + 1

            imjk = funijk(im,j,k)
            ipjk = funijk(ip,j,k)
            ijmk = funijk(i,jm,k)
            ijpk = funijk(i,jp,k)
            ijkm = funijk(i,j,km)
            ijkp = funijk(i,j,kp)

            CALL get_cell_node_coordinates(ijk,'W_MOMENTUM')

!======================================================================
!  Get Interpolation factors at East face
!======================================================================

            u_node_at_te = ((.NOT.blocked_u_cell_at(ijkp)).AND.(.NOT.wall_u_at(ijkp)))
            u_node_at_be = ((.NOT.blocked_u_cell_at(ijk)).AND.(.NOT.wall_u_at(ijk)))


            IF(u_node_at_te.AND.u_node_at_be) THEN
               theta_w_te(ijk) = (z_w_ec(ijk) - z_u(ijk)    ) / (z_u(ijkp) - z_u(ijk))
               theta_w_be(ijk) = one - theta_w_te(ijk)
            ELSE IF (u_node_at_be.AND.(.NOT.u_node_at_te)) THEN
               zi = half * (zt_w_int(ijk) + zt_w_int(ijmk))
               theta_w_te(ijk) = (z_w_ec(ijk) - z_u(ijk)    ) / (zi - z_u(ijk))
               theta_w_be(ijk) = one - theta_w_te(ijk)
            ELSE IF ((.NOT.u_node_at_be).AND.u_node_at_te) THEN
               zi = half * (zt_w_int(ijk) + zt_w_int(ijmk))
               theta_w_te(ijk) = (z_w_ec(ijk) - zi) / (z_u(ijkp) - zi)
               theta_w_be(ijk) = one - theta_w_te(ijk)
            ELSE
               theta_w_te(ijk) = zero
               theta_w_be(ijk) = zero
            ENDIF

            IF ((theta_w_te(ijk)>=one).OR.(theta_w_te(ijk)<=zero)) THEN
               theta_w_te(ijk) = half
               theta_w_be(ijk) = half
            ENDIF

!======================================================================
!  Get Interpolation factors at North face
!======================================================================

            v_node_at_tn = ((.NOT.blocked_v_cell_at(ijkp)).AND.(.NOT.wall_v_at(ijkp)))
            v_node_at_bn = ((.NOT.blocked_v_cell_at(ijk)).AND.(.NOT.wall_v_at(ijk)))


            IF(v_node_at_tn.AND.v_node_at_bn) THEN
               theta_w_tn(ijk) = (z_w_nc(ijk) - z_v(ijk)    ) / (z_v(ijkp) - z_v(ijk))
               theta_w_bn(ijk) = one - theta_w_tn(ijk)
            ELSE IF (v_node_at_bn.AND.(.NOT.v_node_at_tn)) THEN
               zi = half * (zt_w_int(ijk) + zt_w_int(imjk))
               theta_w_tn(ijk) = (z_w_nc(ijk) - z_v(ijk)    ) / (zi - z_v(ijk))
               theta_w_bn(ijk) = one - theta_w_bn(ijk)
            ELSE IF ((.NOT.v_node_at_bn).AND.v_node_at_tn) THEN
               zi = half * (zt_w_int(ijk) + zt_w_int(imjk))
               theta_w_tn(ijk) = (z_w_nc(ijk) - zi) / (z_v(ijkp) - zi)
               theta_w_bn(ijk) = one - theta_w_bn(ijk)
            ELSE
               theta_w_tn(ijk) = zero
               theta_w_bn(ijk) = zero
            ENDIF

            IF ((theta_w_tn(ijk)>=one).OR.(theta_w_tn(ijk)<=zero)) THEN
               theta_w_tn(ijk) = half
               theta_w_bn(ijk) = half
            ENDIF

!======================================================================
!  Get Interpolation factors at Top face
!======================================================================

            theta_wt(ijk) = delz_wt(ijk) / (delz_wt(ijk) + delz_wb(ijkp))
            theta_wt_bar(ijk) = one - theta_wt(ijk)

            IF ((theta_wt(ijk)>=one).OR.(theta_wt(ijk)<=zero)) THEN
               theta_wt(ijk) = half
               theta_wt_bar(ijk) = half
            ENDIF

            bcv = bc_w_id(ijk)
            IF(bcv>0) THEN
               bct = bc_type_enum(bcv)
            ELSE
               bct = blank
            ENDIF

            IF(bct ==cg_nsw.OR.bct ==cg_psw) THEN

               CALL get_del_h(ijk,'W_MOMENTUM',x_w(ijk),y_w(ijk),z_w(ijk),delh_w(ijk),nx,ny,nz)

               IF(delh_w(ijk)<zero.AND.(.NOT.wall_w_at(ijk))) THEN
                  WRITE(*,*) 'NEGATIVE DELH-W AT XYZ=', x_w(ijk),y_w(ijk),z_w(ijk)
                  WRITE(*,*) 'DELH_W=', delh_w(ijk)
                  WRITE(*,*) 'AXY=', axy(ijk)
                  WRITE(*,*) 'MFIX WILL EXIT NOW.'
                  CALL mfix_exit(mype)
               ENDIF

!======================================================================
!  Get Alpha Correction factors at East face
!======================================================================

               theta_we(ijk) = delx_we(ijk) / (delx_we(ijk) + delx_ww(ipjk))
               theta_we_bar(ijk) = one - theta_we(ijk)

               IF ((theta_we(ijk)>=one).OR.(theta_we(ijk)<=zero)) THEN
                  theta_we(ijk) = half
                  theta_we_bar(ijk) = half
               ENDIF

! Location of the interpolated velocity along East face
! Xe has not moved based on definition of Theta_Ve
               xe = x_node(8)
               ye = theta_we_bar(ijk) * y_w(ijk) +theta_we(ijk) * y_w(ipjk)
               ze = theta_we_bar(ijk) * z_w(ijk) +theta_we(ijk) * z_w(ipjk)

               CALL get_del_h(ijk,'W_MOMENTUM',x_w_ec(ijk),y_w_ec(ijk),z_w_ec(ijk),delh_ec,nx,ny,nz)

               CALL get_del_h(ijk,'W_MOMENTUM',xe,ye,ze,delh_e,nx,ny,nz)


               IF((delh_ec == undefined).OR.(delh_e == undefined)) THEN
                  alpha_we_c(ijk) = one
               ELSE
                  alpha_we_c(ijk) = dmin1(alpha_max , delh_ec / delh_e )
               ENDIF

               IF(blocked_w_cell_at(ipjk)) alpha_we_c(ijk) = zero
               IF(wall_w_at(ijk).AND.wall_w_at(ipjk)) alpha_we_c(ijk) = zero
               IF(i == iend1) alpha_we_c(ijk) = one

               IF(alpha_we_c(ijk)<zero) alpha_we_c(ijk) = one

!======================================================================
!  Get Non-Ortogonality correction factor at East face
!======================================================================

               sx = x_w(ipjk) - x_w(ijk)
               sy = y_w(ipjk) - y_w(ijk)
               sz = z_w(ipjk) - z_w(ijk)

               noc_w_e(ijk) = (sy * ny + sz * nz)/(sx * delh_e)


               IF(blocked_w_cell_at(ipjk)) noc_w_e(ijk) = zero
               IF(wall_w_at(ijk).AND.wall_w_at(ipjk)) noc_w_e(ijk) = zero
               IF(i == iend1) noc_w_e(ijk) = zero

!======================================================================
!  Get Alpha Correction factors at North face
!======================================================================

               theta_wn(ijk) = dely_wn(ijk) / (dely_wn(ijk) + dely_ws(ijpk))
               theta_wn_bar(ijk) = one - theta_wn(ijk)

               IF ((theta_wn(ijk)>=one).OR.(theta_wn(ijk)<=zero)) THEN
                  theta_wn(ijk) = half
                  theta_wn_bar(ijk) = half
               ENDIF


! Location of the interpolated velocity along north face
! Yn has not moved based on definition of Theta_V

               xn = theta_wn_bar(ijk) * x_w(ijk) + theta_wn(ijk) * x_w(ijpk)
               yn = y_node(8)
               zn = theta_wn_bar(ijk) * z_w(ijk) + theta_wn(ijk) * z_w(ijpk)

               CALL get_del_h(ijk,'W_MOMENTUM',x_w_nc(ijk),y_w_nc(ijk),z_w_nc(ijk),delh_nc,nx,ny,nz)

               CALL get_del_h(ijk,'W_MOMENTUM',xn,yn,zn,delh_n,nx,ny,nz)

               IF((delh_nc == undefined).OR.(delh_n == undefined)) THEN
                  alpha_wn_c(ijk) = one
               ELSE
                  alpha_wn_c(ijk) = dmin1(alpha_max , delh_nc / delh_n )
               ENDIF

               IF(blocked_w_cell_at(ijpk)) alpha_wn_c(ijk) = zero
               IF(wall_w_at(ijk).AND.wall_w_at(ijpk)) alpha_wn_c(ijk) = zero
               IF(j == jend1) alpha_wn_c(ijk) = one

               IF(alpha_wn_c(ijk)<zero) alpha_wn_c(ijk) = one

!======================================================================
!  Get Non-Ortogonality correction factor at North face
!======================================================================

               sx = x_w(ijpk) - x_w(ijk)
               sy = y_w(ijpk) - y_w(ijk)
               sz = z_w(ijpk) - z_w(ijk)

               noc_w_n(ijk) = (sx * nx + sz * nz)/(sy * delh_n)

               IF(blocked_w_cell_at(ijpk)) noc_w_n(ijk) = zero
               IF(wall_w_at(ijk).AND.wall_w_at(ijpk)) noc_w_n(ijk) = zero
               IF(j == jend1) noc_w_n(ijk) = zero

!======================================================================
!  Get Alpha Correction factors at Top face
!======================================================================

! Location of the interpolated velocity along Top face
! Zt has not moved based on definition of Theta_Ut

               xt = theta_wt_bar(ijk) * x_w(ijk) + theta_wt(ijk) * x_w(ijkp)
               yt = theta_wt_bar(ijk) * y_w(ijk) + theta_wt(ijk) * y_w(ijkp)
               zt = z_node(8)

               CALL get_del_h(ijk,'W_MOMENTUM',x_w_tc(ijk),y_w_tc(ijk),z_w_tc(ijk),delh_tc,nx,ny,nz)

               CALL get_del_h(ijk,'W_MOMENTUM',xt,yt,zt,delh_t,nx,ny,nz)

               IF((delh_tc == undefined).OR.(delh_t == undefined)) THEN
                  alpha_wt_c(ijk) = one
               ELSE
                  alpha_wt_c(ijk) = dmin1(alpha_max , delh_tc / delh_t )
               ENDIF

               IF(blocked_w_cell_at(ijkp)) alpha_wt_c(ijk) = zero
               IF(wall_w_at(ijk).AND.wall_w_at(ijkp)) alpha_wt_c(ijk) = zero
               IF(k == kend1) alpha_wt_c(ijk) = one

               IF(alpha_wt_c(ijk)<zero) THEN
                  WRITE(*,*) 'NEGATIVE ALPHA_Wt_c at IJK=',ijk
                  WRITE(*,*) 'MFIX WILL EXIT NOW.'
                  CALL mfix_exit(mype)
               ENDIF

!======================================================================
!  Get Non-Ortogonality correction factor at Top face
!======================================================================

               sx = x_w(ijkp) - x_w(ijk)
               sy = y_w(ijkp) - y_w(ijk)
               sz = z_w(ijkp) - z_w(ijk)

               noc_w_t(ijk) = (sx * nx + sy * ny)/(sz * delh_t)

               IF(blocked_w_cell_at(ijkp)) noc_w_t(ijk) = zero
               IF(wall_w_at(ijk).AND.wall_w_at(ijkp)) noc_w_t(ijk) = zero
               IF(k == kend1) noc_w_t(ijk) = zero

            ENDIF

!======================================================================
!  Get Surfaces used to compute pressure gradient
!======================================================================

            IF (  cut_w_cell_at(ijk)  )  THEN

               SELECT CASE (pg_option)

                  CASE(0)                               ! do nothing
                                                     ! will be treated below

                  CASE(1)

                     a_wpg_t(ijk) = dmax1(axy_w(ijk),axy_w(ijkm))
                     a_wpg_b(ijk) = a_wpg_t(ijk)

                  CASE(2)

                     a_wpg_t(ijk) = axy_w(ijk)
                     a_wpg_b(ijk) = axy_w(ijkm)

                  CASE DEFAULT

                     WRITE(*,*)'INVALID PRESSURE GRADIENT OPTION:',pg_option
                     WRITE(*,*)'PG_OPTION SHOULD BE SET EQUAL TO 0, 1 OR 2.'
                     WRITE(*,*)'PLEASE VERIFY MFIX.DAT FILE AND TRY AGAIN.'
                     WRITE(*,*) 'MFIX WILL EXIT NOW.'
                     CALL mfix_exit(mype)

               END SELECT

               IF (blocked_cell_at(ijk).OR.blocked_cell_at(ijkp)) THEN
                  IF(.NOT.wall_w_at(ijk)) THEN
                     wall_w_at(ijk) = .true.
                     flag_t(ijk) = 0
                     IF(print_warnings) THEN
                        WRITE(*,*) 'WARNING: ONLY ONE PRESSURE NODE DETECTED IN W-MOMENTUM CELL IJK =',ijk
                        WRITE(*,*) 'RESETTING U-MOMENTUM CELL AS W_WALL CELL.'
                     ENDIF
!                     write(*,*) 'MFiX will exit now.'
!                     CALL MFIX_EXIT(myPE)
                  ENDIF
               ENDIF

            ELSE

               a_wpg_t(ijk) = axy_w(ijk)
               a_wpg_b(ijk) = axy_w(ijk)

            ENDIF

         ENDIF
      END DO

      IF(pg_option==0) THEN
         a_wpg_t = axy
         a_wpg_b = axy
      ENDIF

      RETURN

      END SUBROUTINE get_3d_alpha_w_cut_cell