get_cut_cell_flags.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: SET_3D_CUT_CELL_FLAGS                                  C
!  Purpose: Set flags for scalar cut cells, based on intersection      C
!  of the grid with the quadric(s)                                     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 set_3d_cut_cell_flags

      USE compar, ONLY: mype, pe_io, ijkstart3, ijkend3, istart1, iend1, jstart1, jend1, kstart1, kend1
      USE cutcell
      USE functions, ONLY: funijk, ip_of, jp_of, kp_of, bottom_of, south_of, west_of, fluid_at, is_on_mype_wobnd
      USE geometry, ONLY: do_i, do_j, do_k, imin1, imax3, jmin1, jmax3, kmin1, kmax3, no_k, vol, axy, axz, ayz, dx, dy, dz, flag
      USE indices, ONLY: i_of, j_of, k_of
      USE param, only: dimension_3
      USE param1, only: zero, half, one, undefined
      USE polygon, ONLY: n_polygon
      USE quadric, ONLY: tol_f
      USE sendrecv
      USE cut_cell_preproc, ONLY: cad_intersect, clean_intersect, clean_intersect_scalar, eval_f, intersect
      USE vtk, ONLY: global_var_allocated, grid_info_printed_on_screen

      IMPLICIT NONE
      INTEGER :: IJK,I,J,K
      INTEGER :: TOTAL_NUMBER_OF_INTERSECTIONS
      INTEGER :: NODE,N_N1,N_N2,Q_ID
      INTEGER :: MIN_INTERSECTIONS,MAX_INTERSECTIONS
      LOGICAL :: CLIP_FLAG
      INTEGER :: BCID

      INTEGER, DIMENSION(DIMENSION_3,15) :: OLD_CONNECTIVITY
      DOUBLE PRECISION, DIMENSION(:), allocatable :: X_OLD_POINT,Y_OLD_POINT,Z_OLD_POINT

      INTEGER, DIMENSION(6) :: NB

      INTEGER :: IJK_NB,NODE_NB,NN,N_NB,NC

      DOUBLE PRECISION :: X1,Y1,Z1,X2,Y2,Z2,D,TOT_VOL_NODE,TOT_VOL_CELLS

      DOUBLE PRECISION, DIMENSION(:,:), allocatable :: SCALAR_NODE_XYZ_TEMP

      DOUBLE PRECISION :: DIST, NORM1, NORM2, NORM3,Diagonal
      INTEGER :: IJK2, I1, I2, J1, J2, K1, K2, II, JJ, KK
      LOGICAL :: COND_1, COND_2

      allocate(x_old_point(dimension_max_cut_cell))
      allocate(y_old_point(dimension_max_cut_cell))
      allocate(z_old_point(dimension_max_cut_cell))
      allocate(scalar_node_xyz_temp(dimension_3, 3))

      IF(mype == pe_io) THEN
         WRITE(*,10)'INTERSECTING GEOMETRY WITH SCALAR CELLS...'
      ENDIF
10    FORMAT(1x,a)

      global_var_allocated = .false.
      grid_info_printed_on_screen = .false.

!  Set arrays for computing indices
!

!      CALL SET_INCREMENTS

      dx(imax3+1) = dx(imax3)
      dy(jmax3+1) = dy(jmax3)
      dz(kmax3+1) = dz(kmax3)

!     x-Location of U-momentum cells for original (uncut grid)
      IF (do_i) THEN
         xg_e(1) = zero
         DO i = imin1, imax3
            xg_e(i) = xg_e(i-1) + dx(i)
         END DO
      ENDIF


!     y-Location of V-momentum cells for original (uncut grid)
      IF (do_j) THEN
         yg_n(1) = zero
         DO j = jmin1, jmax3
            yg_n(j) = yg_n(j-1) + dy(j)
         END DO
      ENDIF


!     z-Location of W-momentum cells for original (uncut grid)
      IF (do_k) THEN
         zg_t(1) = zero
         DO k = kmin1, kmax3
            zg_t(k) = zg_t(k-1) + dz(k)
         END DO
      ENDIF

      partition = dble(mype)       ! ASSIGN processor ID (for vizualisation)

      CALL set_geometry1   ! INITIALIZE ALL VOLUMES AND AREAS


!======================================================================
!  Grid coordinates:
!======================================================================

      small_cell_at = .false.
      small_cell_flag = 0
      number_of_small_cells = 0

!======================================================================
!  Intersection between quadric Grid
!======================================================================

      intersect_x = .false.
      intersect_y = .false.
      intersect_z = .false.
      snap = .false.

      DO ijk = ijkstart3, ijkend3

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

         IF(no_k) THEN   ! 2D case

            interior_cell_at(ijk) = (     (i >= istart1 ).AND.(i <= iend1 )  &
                                     .AND.(j >= jstart1 ).AND.(j <= jend1 ) )

         ELSE            ! 3D case

            interior_cell_at(ijk) = (     (i >= istart1 ).AND.(i <= iend1 )  &
                                     .AND.(j >= jstart1 ).AND.(j <= jend1 )  &
                                     .AND.(k >= kstart1 ).AND.(k <= kend1 ) )

         ENDIF

      END DO

       number_of_nodes = 0

      CALL get_potential_cut_cells

!       POTENTIAL_CUT_CELL_AT=.TRUE.

      IF(.NOT.(use_stl.OR.use_msh)) THEN
         DO ijk = ijkstart3, ijkend3
            IF(potential_cut_cell_at(ijk))  CALL intersect(ijk,'SCALAR',xn_int(ijk),ye_int(ijk),zt_int(ijk))
!            CALL INTERSECT(IJK,'SCALAR',Xn_int(IJK),Ye_int(IJK),Zt_int(IJK))
         END DO

!======================================================================
!  Clean-up intersection flags by snapping intersection points to close
!  corner (within tolerance TOL_SNAP)
!======================================================================
         CALL clean_intersect_scalar

      ELSE
         CALL cad_intersect('SCALAR',xn_int,ye_int,zt_int)
      ENDIF
!        NUMBER_OF_NODES = 0

      DO ijk = ijkstart3, ijkend3

          IF(potential_cut_cell_at(ijk))  THEN

         CALL write_progress_bar(ijk,ijkend3 - ijkstart3 + 1,'C')

!======================================================================
!  Get coordinates of eight nodes
!======================================================================

         CALL get_cell_node_coordinates(ijk,'SCALAR')

         scalar_node_xyz_temp(ijk,1) = x_node(8)
         scalar_node_xyz_temp(ijk,2) = y_node(8)
         scalar_node_xyz_temp(ijk,3) = z_node(8)
!======================================================================
!  Initialize location of velocity nodes
!======================================================================

         x_u(ijk) = x_node(8)
         y_u(ijk) = half * (y_node(6) + y_node(8))
         z_u(ijk) = half * (z_node(4) + z_node(8))

         x_v(ijk) = half * (x_node(7) + x_node(8))
         y_v(ijk) = y_node(8)
         z_v(ijk) = half * (z_node(4) + z_node(8))

         x_w(ijk) = half * (x_node(7) + x_node(8))
         y_w(ijk) = half * (y_node(6) + y_node(8))
         z_w(ijk) = z_node(8)

         IF(interior_cell_at(ijk)) THEN
!======================================================================
!  Get Connectivity
!======================================================================

            IF(no_k) THEN
               min_intersections = 2
               max_intersections = 2
               n_n1 = 5
               n_n2 = 8

            ELSE
               min_intersections = 3
               max_intersections = 6
               n_n1 = 1
               n_n2 = 8
            ENDIF

            CALL get_connectivity(ijk,'SCALAR',number_of_new_points,number_of_nodes(ijk),connectivity,&
            x_new_point,y_new_point,z_new_point,total_number_of_intersections,xn_int,ye_int,zt_int)

            IF(total_number_of_intersections < min_intersections ) THEN   ! Not a cut cell

               q_id = 1
               CALL eval_f('QUADRIC',x_node(0),y_node(0),z_node(0),q_id,f_node(0),clip_flag)

               CALL eval_f('POLYGON',x_node(0),y_node(0),z_node(0),n_polygon,f_node(0),clip_flag)

               CALL eval_f('USR_DEF',x_node(0),y_node(0),z_node(0),n_usr_def,f_node(0),clip_flag)

!               CALL EVAL_STL_FCT_AT('SCALAR',IJK,0,F_NODE(0),CLIP_FLAG,BCID)

               IF(use_msh.OR.use_stl) THEN
!                  CALL EVAL_STL_FCT_AT('SCALAR',IJK,1,F_NODE(1),CLIP_FLAG,BCID)
!                  CALL EVAL_STL_FCT_AT('SCALAR',IJK,8,F_NODE(8),CLIP_FLAG,BCID)
!                  F_NODE(0) = HALF*(F_NODE(1)+F_NODE(8))
! Pick the first defined value of F_AT at one of the cell corners
! This should avoid wrongly setting a cell next to the global ghost layer as blocked cell
! when the average between undefined and negative is a positive number
! This occurs for external flows only due to F_AT being initialized as UNDEFINED
! which is a positive number (considered outside the fluid region).
                  DO node=1,8
                     CALL eval_stl_fct_at('SCALAR',ijk,node,f_node(node),clip_flag,bcid)
                     IF(f_node(node)/=undefined.AND.f_node(node)/=zero) THEN
                        f_node(0) = f_node(node)
                        EXIT
                     ENDIF

                  ENDDO
               ENDIF


!  Uncomment the line below to force the cell to be a fluid cell
!  Currently, F_NODE(0) > 0 ==> Blocked cell
!             F_NODE(0) < 0 ==> Fluid cell

!               IF(TOTAL_NUMBER_OF_INTERSECTIONS==-1) F_NODE(0) = -ONE

               bc_id(ijk) = 0

!               IF(F_NODE(0) < TOL_F) THEN
               IF(f_node(0) < zero) THEN
                  IF((flag(ijk)>=100).AND.(flag(ijk)<=102)) THEN
                     blocked_cell_at(ijk) = .true.
                  ELSE
                     blocked_cell_at(ijk) = .false.
                     standard_cell_at(ijk) = .true.           ! Regular fluid cell
                  ENDIF
               ELSE
                  flag(ijk) = 100
                  blocked_cell_at(ijk) = .true.               ! Blocked fluid cell
                  standard_cell_at(ijk) = .false.
                  axy(ijk) = zero
                  axz(ijk) = zero
                  ayz(ijk) = zero
                  vol(ijk) = zero

                  axy(bottom_of(ijk)) = zero
                  axz(south_of(ijk)) = zero
                  ayz(west_of(ijk)) = zero
               ENDIF

               IF(no_k) THEN
                  number_of_nodes(ijk) = 4
                  connectivity(ijk,1) = ijk_of_node(5)
                  connectivity(ijk,2) = ijk_of_node(6)
                  connectivity(ijk,3) = ijk_of_node(8)
                  connectivity(ijk,4) = ijk_of_node(7)
               ELSE
                  number_of_nodes(ijk) = 8
                  DO node = n_n1,n_n2
                     connectivity(ijk,node) = ijk_of_node(node)
                  END DO
               ENDIF


            ELSE IF(total_number_of_intersections > max_intersections ) THEN

               IF(print_warnings) THEN

                  WRITE(*,*) 'TOO MANY INTERSECTIONS FOUND IN CELL IJK = ',ijk
                  WRITE(*,*) 'MAXIMUM NUMBER OF INTERSECTIONS = ',max_intersections
                  WRITE(*,*) 'TOTAL NUMBER OF INTERSECTIONS = ',total_number_of_intersections
                  WRITE(*,*) 'REMOVING SCALAR CELL FROM COMPUTATION.'
!                  WRITE(*,*) 'MFIX WILL EXIT NOW.'
!                  CALL MFIX_EXIT(MYPE)

               ENDIF

               blocked_cell_at(ijk) = .true.               ! Blocked fluid cell

            ELSE                                         ! Cut cell

               cut_cell_at(ijk) = .true.
               blocked_cell_at(ijk) = .false.
               standard_cell_at(ijk) = .false.

               DO node = n_n1,n_n2
                  IF(f_node(node) < - tol_f) THEN
                     IF(.NOT.snap(ijk_of_node(node))) THEN
                        number_of_nodes(ijk) = number_of_nodes(ijk) + 1
                        connectivity(ijk,number_of_nodes(ijk)) = ijk_of_node(node)
                     ENDIF
                  ENDIF
               END DO

               CALL get_cut_cell_volume_and_areas(ijk,'SCALAR',number_of_nodes(ijk),connectivity,&
               x_new_point,y_new_point,z_new_point)

            ENDIF

         ENDIF        ! Interior cell

      ENDIF
      END DO          ! IJK Loop


      call send_receive_1d_logical(small_cell_at,2)
      call send_receive_1d_logical(standard_cell_at,2)
      call send_receive_1d_logical(blocked_cell_at,2)
      call send_receive_1d_logical(cut_cell_at,2)

! Consolidate blocked cells
      DO ijk = ijkstart3, ijkend3
         IF(blocked_cell_at(ijk)) THEN
            standard_cell_at(ijk) = .false.
            cut_cell_at(ijk)      = .false.

            flag(ijk)             = 100

            vol(ijk)              = zero
            axy(ijk)              = zero
            axz(ijk)              = zero
            ayz(ijk)              = zero

            axy(bottom_of(ijk))   = zero
            axz(south_of(ijk))    = zero
            ayz(west_of(ijk))     = zero
         ENDIF
      ENDDO

      call send_receive_1d_logical(small_cell_at,2)
      call send_receive_1d_logical(standard_cell_at,2)
      call send_receive_1d_logical(blocked_cell_at,2)
      call send_receive_1d_logical(cut_cell_at,2)

      call send_recv(flag,2)
      call send_recv(vol,2)
      call send_recv(axy,2)
      call send_recv(axz,2)
      call send_recv(ayz,2)

!      print*,'Before removing duplicate nodes:'
      DO ijk = ijkstart3, ijkend3,-1
         IF(cut_cell_at(ijk)) THEN
            print*,'==========================================================='
            print*,'IJK,  I,J=',ijk,i_of(ijk),j_of(ijk)
            print*,'NUMBER_OF_NODES=',number_of_nodes(ijk)
            DO node = 1,number_of_nodes(ijk)
               IF(connectivity(ijk,node)>ijkend3) THEN
                  print*,'CNCT=',node,connectivity(ijk,node), &
                       x_new_point(connectivity(ijk,node)-ijkend3),y_new_point(connectivity(ijk,node)-ijkend3)
               ELSE
                  print*,'CNCT=',node,connectivity(ijk,node)
               ENDIF
            ENDDO
            print*,''
         ENDIF
      ENDDO



      old_connectivity = connectivity
!      X_OLD_POINT      = X_NEW_POINT
!      Y_OLD_POINT      = Y_NEW_POINT
!      Z_OLD_POINT      = Z_NEW_POINT


!======================================================================
!  Removing duplicate new points
!  Up to here, a cut face node that is shared among neighbor cells had
!  a different index, and appeared as a duplicate node.
!  The list will be updated and duplicate nodes will be removed from
!  the connectivity list
!  This is done to ensure that we can assign a volume surrounding each
!  node, that will be used to compute solids volume fraction for MPPIC
!======================================================================

!      print*,'Removing duplicate nodes...'

      DO ijk = ijkstart3, ijkend3
         IF(cut_cell_at(ijk)) THEN          ! for each cut cell, identify neibhor cells that are also cut cells
                                             ! Look east and north in 2D, and also Top in 3D

         diagonal = dsqrt(dx(i_of(ijk))**2 + dy(j_of(ijk))**2 + dz(k_of(ijk))**2)

            nb(1) = ip_of(ijk)
            nb(2) = jp_of(ijk)

            IF(no_k) THEN
               n_nb = 2
            ELSE
               n_nb = 3
               nb(3) = kp_of(ijk)
            ENDIF


            DO nn = 1,n_nb
               IF(cut_cell_at(nb(nn))) THEN   ! For two neighbor cut cells, compare each cut-face node to remove duplicates

                  ijk_nb = nb(nn)

!                  print*,'comparing:',IJK,' and',IJK_NB

                  DO node = 1,number_of_nodes(ijk)
                     IF(old_connectivity(ijk,node)>ijkend3) THEN  ! node belongs to the cut-face

                        x1 = x_new_point(old_connectivity(ijk,node)-ijkend3)
                        y1 = y_new_point(old_connectivity(ijk,node)-ijkend3)
                        z1 = z_new_point(old_connectivity(ijk,node)-ijkend3)


                        DO node_nb = 1,number_of_nodes(ijk_nb)
                           IF(old_connectivity(ijk_nb,node_nb)>ijkend3) THEN  ! node belongs to the cut-face

                              x2 = x_new_point(old_connectivity(ijk_nb,node_nb)-ijkend3)
                              y2 = y_new_point(old_connectivity(ijk_nb,node_nb)-ijkend3)
                              z2 = z_new_point(old_connectivity(ijk_nb,node_nb)-ijkend3)

                              ! compare coordinates of cut-face nodes
                              d = (x2-x1)**2 + (y2-y1)**2 + (z2-z1)**2

                              ! Duplicate nodes have identical coordinates (within tolerance TOL_MERGE times diagonal)
                              IF(d<tol_merge*diagonal) THEN ! keep the smallest node ID
                                 !  print*,'DULICATE NODES:', &
                                 ! NODE,NODE_NB,OLD_CONNECTIVITY(IJK,NODE),OLD_CONNECTIVITY(IJK_NB,NODE_NB)
                                 nc = min(old_connectivity(ijk,node),old_connectivity(ijk_nb,node_nb))
                                 connectivity(ijk   ,node   ) = nc
                                 connectivity(ijk_nb,node_nb) = nc

                              ENDIF
                           ENDIF
                        ENDDO

                     ENDIF
                  ENDDO

               ENDIF
            ENDDO

         ENDIF
      ENDDO

      ALLOCATE(scalar_node_xyz(dimension_3 + number_of_new_points,3))
      ALLOCATE(ovol_around_node(dimension_3 + number_of_new_points))
      ALLOCATE(scalar_node_atwall(dimension_3 + number_of_new_points))

      !first fill with standard nodes
      DO ijk = ijkstart3, ijkend3
         scalar_node_xyz(ijk,1:3)  = scalar_node_xyz_temp(ijk,1:3)
      ENDDO

      !now fill with cut-face nodes
      DO ijk = 1,number_of_new_points
         scalar_node_xyz(ijkend3+ijk,1) = x_new_point(ijk)
         scalar_node_xyz(ijkend3+ijk,2) = y_new_point(ijk)
         scalar_node_xyz(ijkend3+ijk,3) = z_new_point(ijk)
      ENDDO



      scalar_node_atwall(:)  = .true.
      !Rahul:
      !One could either set all scalar_node_atwall to false and
      !then set to true the nodes that are found as on the wall or outside
      !the doamin
      !In some situations, a node can be deemed to be both inside
      !and outside the domain depending upon which cell you look from.
      !Think of a small cell. According to the small cells, all the nodes
      !are outside the domain, but from the perspective or surrounding
      !cut-cells, some of the nodes of this small cell are within the domain.

      !so I'm setting all the points as being outside the domain.
      !if a point is ever found to be in the domain, then it will stay away
      !and further tests will not be able to revert it.
      ijkloop: DO ijk = ijkstart3, ijkend3
         i = i_of(ijk)
         j = j_of(ijk)
         k = k_of(ijk)
         IF(.not. is_on_mype_wobnd(i,j,k)) cycle

         i1 = i-1
         i2 = i
         j1 = j-1
         j2 = j

         IF(no_k) THEN
            k1 = k
            k2 = k
         ELSE
            k1 = k-1
            k2 = k
         ENDIF
         !Convention used to number node numbers is described below

         ! i=1, j=2           i=2, j=2
         !   _____________________
         !   |                   |
         !   |  I = 2, J = 2     |
         !   |___________________|
         ! i=1, j=1           i=2, j=1

         !Let's say the scalar cell with I = 2 and J = 2, i.e., the
         !first scalar cell in either direction.
         !then this scalar cell's node numbering in x- direction
         !will go from 1 to 2 and likewise in other directions.
         DO kk = k1, k2
            DO jj = j1, j2
               DO ii = i1, i2
                  ijk2 = funijk(ii, jj, kk)
                  cond_1 = .false.
                  cond_2 = .false.
                  !if it was already found to be inside the domain, then don't bother
                  IF(.not.scalar_node_atwall(ijk2))  cycle

                  IF(.not.fluid_at(ijk)) THEN
                     !do nothing
                  ELSE
                     !IJK is a fluid scalar cell. Check if it is
                     !a cut-cell or not
                     IF(cut_cell_at(ijk)) THEN
                        CALL get_del_h_des(ijk,'SCALAR', &
                        & scalar_node_xyz(ijk2,1), scalar_node_xyz(ijk2,2), &
                        & scalar_node_xyz(ijk2,3), &
                        & dist, norm1, norm2, norm3, .true.)
                        IF(dist.GE.zero) THEN
                           scalar_node_atwall(ijk2)  = .false.
                           cond_1 =  .true.
                        ENDIF
                     ELSE
                        scalar_node_atwall(ijk2)  = .false.
                        cond_2 = .true.
                     ENDIF
                  ENDIF
                  !if(II == 1) then
                  ! write(*,'(10x, A, 4(2x,i10),3(2x,L2))') &
                  ! 'I1,J1, I, J, ATWALL, COND1, COND2 =  ', II, JJ, I, J,  SCALAR_NODE_ATWALL(IJK2), COND_1, COND_2
                  !endif
               ENDDO
            ENDDO
         ENDDO

      ENDDO ijkloop

!      print*,'After removing duplicate nodes:'
      DO ijk = ijkstart3, ijkend3,-1
         print*,'==========================================================='
         print*,'IJK,  I,J=',ijk,i_of(ijk),j_of(ijk)
         print*,'NUMBER_OF_NODES=',number_of_nodes(ijk)
         DO node = 1,number_of_nodes(ijk)
            print*,'CNCT=',node,connectivity(ijk,node), &
                 scalar_node_xyz(connectivity(ijk,node),1),scalar_node_xyz(connectivity(ijk,node),2)
         ENDDO
         print*,''
      ENDDO



      ovol_around_node = zero !UNDEFINED

      tot_vol_cells = zero

      DO ijk = ijkstart3, ijkend3

         DO node = 1,number_of_nodes(ijk)
            nc = connectivity(ijk,node)
            ovol_around_node(nc) = ovol_around_node(nc) + vol(ijk)/number_of_nodes(ijk)
         ENDDO

!         print*,'IJK,VOL=',IJK,VOL(IJK)


         IF(interior_cell_at(ijk)) tot_vol_cells = tot_vol_cells + vol(ijk)


      ENDDO


      tot_vol_node= zero

      DO ijk = ijkstart3, ijkend3 + number_of_new_points    ! Loop over all nodes

         IF(ovol_around_node(ijk)>zero) THEN
!             print*,'NODE,VOL=',IJK,Ovol_around_node(IJK)
             tot_vol_node = tot_vol_node + ovol_around_node(ijk)
             ovol_around_node(ijk) = one / ovol_around_node(ijk)    ! Store One/volume
         ENDIF


      ENDDO

      IF(tot_vol_cells>zero) THEN
         IF(dabs(one-dabs(tot_vol_node/tot_vol_cells))>1.0d-6) THEN
            WRITE(*,*)'==========================================================='
            WRITE(*,*)'Total volume around nodes=',tot_vol_node
            WRITE(*,*)'Total volume =',tot_vol_cells
            WRITE(*,*)'The two volumes above should be the same.'
            WRITE(*,*)'==========================================================='
            CALL mfix_exit(mype)
         ENDIF
      ENDIF

      deallocate(x_old_point)
      deallocate(y_old_point)
      deallocate(z_old_point)
      deallocate(scalar_node_xyz_temp)

      RETURN
      END SUBROUTINE set_3d_cut_cell_flags

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: SET_3D_CUT_U_CELL_FLAGS                                C
!  Purpose: Set flags for U-Momentum cut cells, based on intersection  C
!  of the grid with the quadric(s)                                     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 set_3d_cut_u_cell_flags

      USE compar, ONLY: mype, pe_io, ijkstart3, ijkend3
      USE cut_cell_preproc, ONLY: cad_intersect, clean_intersect, clean_intersect_scalar, eval_f, intersect
      USE cutcell
      USE geometry, ONLY: no_k, axy_u, ayz_u, vol_u, axz_u
      USE param1, ONLY: half, one, zero
      USE polygon, ONLY: n_polygon
      USE quadric, ONLY: tol_f

      IMPLICIT NONE
      INTEGER :: IJK
      INTEGER :: TOTAL_NUMBER_OF_INTERSECTIONS
      INTEGER :: NODE,N_N1,N_N2,Q_ID
      INTEGER :: MIN_INTERSECTIONS,MAX_INTERSECTIONS
      LOGICAL :: CLIP_FLAG
      INTEGER :: BCID

      IF(mype == pe_io) THEN
         WRITE(*,10)'INTERSECTING GEOMETRY WITH U-MOMENTUM CELLS...'
      ENDIF
10    FORMAT(1x,a)
!======================================================================
!  Intersection between quadric Grid
!======================================================================

      intersect_x = .false.
      intersect_y = .false.
      intersect_z = .false.
      snap = .false.
      tol_snap = zero

      IF(.NOT.(use_stl.OR.use_msh)) THEN
         DO ijk = ijkstart3, ijkend3
!             IF(POTENTIAL_CUT_CELL_AT(IJK))  CALL INTERSECT(IJK,'U_MOMENTUM',Xn_U_int(IJK),Ye_U_int(IJK),Zt_U_int(IJK))
            CALL intersect(ijk,'U_MOMENTUM',xn_u_int(ijk),ye_u_int(ijk),zt_u_int(ijk))
         END DO

!======================================================================
!  Clean-up intersection flags in preparaton of small cells removal
!======================================================================
         DO ijk = ijkstart3, ijkend3
            IF(interior_cell_at(ijk)) THEN
!               IF(POTENTIAL_CUT_CELL_AT(IJK))  CALL CLEAN_INTERSECT(IJK,'U_MOMENTUM',Xn_U_int(IJK),Ye_U_int(IJK),Zt_U_int(IJK))
               CALL clean_intersect(ijk,'U_MOMENTUM',xn_u_int(ijk),ye_u_int(ijk),zt_u_int(ijk))
            ENDIF
         END DO

      ELSE
         CALL cad_intersect('U_MOMENTUM',xn_u_int,ye_u_int,zt_u_int)
      ENDIF


      number_of_new_u_points = 0

      DO ijk = ijkstart3, ijkend3

         CALL write_progress_bar(ijk,ijkend3 - ijkstart3 + 1,'C')

         IF(interior_cell_at(ijk)) THEN

!======================================================================
!  Get coordinates of eight nodes
!======================================================================

            CALL get_cell_node_coordinates(ijk,'U_MOMENTUM')

!======================================================================
!  Initialize location of velocity nodes at center of E,W,T faces
!======================================================================

                  x_u_ec(ijk) = x_node(8)
                  y_u_ec(ijk) = half * (y_node(6) + y_node(8))
                  z_u_ec(ijk) = half * (z_node(4) + z_node(8))

                  x_u_nc(ijk) = half * (x_node(7) + x_node(8))
                  y_u_nc(ijk) = y_node(8)
                  z_u_nc(ijk) = half * (z_node(4) + z_node(8))

                  x_u_tc(ijk) = half * (x_node(7) + x_node(8))
                  y_u_tc(ijk) = half * (y_node(6) + y_node(8))
                  z_u_tc(ijk) = z_node(8)

!======================================================================
!  Get Connectivity
!======================================================================

            IF(no_k) THEN
               min_intersections = 2
               max_intersections = 2
               n_n1 = 5
               n_n2 = 8

            ELSE
               min_intersections = 3
               max_intersections = 6
               n_n1 = 1
               n_n2 = 8
            ENDIF

            CALL get_connectivity(ijk,'U_MOMENTUM',number_of_new_u_points,number_of_u_nodes(ijk),connectivity_u,&
            x_new_u_point,y_new_u_point,z_new_u_point,total_number_of_intersections,xn_u_int,ye_u_int,zt_u_int)


            IF(total_number_of_intersections < min_intersections ) THEN   ! Not a cut cell

               q_id = 1
               CALL eval_f('QUADRIC',x_node(0),y_node(0),z_node(0),q_id,f_node(0),clip_flag)

               CALL eval_f('POLYGON',x_node(0),y_node(0),z_node(0),n_polygon,f_node(0),clip_flag)

               CALL eval_f('USR_DEF',x_node(0),y_node(0),z_node(0),n_usr_def,f_node(0),clip_flag)

               CALL eval_stl_fct_at('U_MOMENTUM',ijk,0,f_node(0),clip_flag,bcid)

               IF(total_number_of_intersections==-1) f_node(0) = -one
               bc_u_id(ijk) = 0

               IF(f_node(0) <= zero) THEN
                  blocked_u_cell_at(ijk) = .false.
                  standard_u_cell_at(ijk) = .true.          ! Regular fluid cell

                  x_u_ec(ijk) = x_node(8)
                  y_u_ec(ijk) = half * (y_node(6) + y_node(8))
                  z_u_ec(ijk) = half * (z_node(4) + z_node(8))

                  x_u_nc(ijk) = half * (x_node(7) + x_node(8))
                  y_u_nc(ijk) = y_node(8)
                  z_u_nc(ijk) = half * (z_node(4) + z_node(8))

                  x_u_tc(ijk) = half * (x_node(7) + x_node(8))
                  y_u_tc(ijk) = half * (y_node(6) + y_node(8))
                  z_u_tc(ijk) = z_node(8)

               ELSE
                  blocked_u_cell_at(ijk) = .true.           ! Blocked fluid cell
                  standard_u_cell_at(ijk) = .false.
                  axy_u(ijk) = zero
                  axz_u(ijk) = zero
                  ayz_u(ijk) = zero
                  vol_u(ijk) = zero
               ENDIF

               IF(no_k) THEN
                  number_of_u_nodes(ijk) = 4
                  connectivity_u(ijk,1) = ijk_of_node(5)
                  connectivity_u(ijk,2) = ijk_of_node(6)
                  connectivity_u(ijk,3) = ijk_of_node(8)
                  connectivity_u(ijk,4) = ijk_of_node(7)
               ELSE
                  number_of_u_nodes(ijk) = 8
                  DO node = n_n1,n_n2
                     connectivity_u(ijk,node) = ijk_of_node(node)
                  END DO
               ENDIF

            ELSE IF(total_number_of_intersections > max_intersections) THEN

               IF(print_warnings) THEN
                  WRITE(*,*) 'TOO MANY INTERSECTIONS FOUND IN U-CELL IJK = ',ijk
                  WRITE(*,*) 'MAXIMUM NUMBER OF INTERSECTIONS = ',max_intersections
                  WRITE(*,*) 'TOTAL NUMBER OF INTERSECTIONS = ',total_number_of_intersections
                  WRITE(*,*) 'REMOVING U-CELL FROM COMPUTATION.'
!                  WRITE(*,*) 'MFIX WILL EXIT NOW.'
!                  CALL MFIX_EXIT(MYPE)

               ENDIF

               blocked_u_cell_at(ijk) = .true.           ! Blocked fluid cell
               standard_u_cell_at(ijk) = .false.
               cut_u_cell_at(ijk) = .false.
               axy_u(ijk) = zero
               axz_u(ijk) = zero
               ayz_u(ijk) = zero
               vol_u(ijk) = zero

            ELSE                                         ! Cut cell

               cut_u_cell_at(ijk) = .true.
               blocked_u_cell_at(ijk) = .false.
               standard_u_cell_at(ijk) = .false.

               DO node = n_n1,n_n2
                  IF(f_node(node) < - tol_f) THEN
                     number_of_u_nodes(ijk) = number_of_u_nodes(ijk) + 1
                     connectivity_u(ijk,number_of_u_nodes(ijk)) = ijk_of_node(node)
                  ENDIF
               END DO

               CALL get_cut_cell_volume_and_areas(ijk,'U_MOMENTUM',number_of_u_nodes(ijk),connectivity_u,&
               x_new_u_point,y_new_u_point,z_new_u_point)

            ENDIF

         ENDIF

      END DO

      DO ijk = ijkstart3, ijkend3

         CALL get_cell_node_coordinates(ijk,'U_MOMENTUM')

         delx_ue(ijk) = x_node(8) - x_u(ijk)
         delx_uw(ijk) = x_u(ijk) - x_node(1)

         dely_un(ijk) = y_node(8) - y_u(ijk)
         dely_us(ijk) = y_u(ijk) - y_node(1)

         delz_ut(ijk) = z_node(8) - z_u(ijk)
         delz_ub(ijk) = z_u(ijk) - z_node(1)

      ENDDO

      RETURN

      END SUBROUTINE set_3d_cut_u_cell_flags

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: SET_3D_CUT_V_CELL_FLAGS                                C
!  Purpose: Set flags for V-Momentum cut cells, based on intersection  C
!  of the grid with the quadric(s)                                     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 set_3d_cut_v_cell_flags

      USE compar, ONLY: mype, pe_io, ijkstart3, ijkend3
      USE cut_cell_preproc, ONLY: cad_intersect, clean_intersect, clean_intersect_scalar, eval_f, intersect
      USE cutcell
      USE geometry, ONLY: no_k, axy_v, axz_v, ayz_v, vol_v
      USE param1, ONLY: half, one, zero
      USE polygon, ONLY: n_polygon
      USE quadric, ONLY: tol_f

      IMPLICIT NONE
      INTEGER :: IJK
      INTEGER :: TOTAL_NUMBER_OF_INTERSECTIONS
      INTEGER :: NODE,N_N1,N_N2,Q_ID
      INTEGER :: MIN_INTERSECTIONS,MAX_INTERSECTIONS
      LOGICAL :: CLIP_FLAG
      INTEGER :: BCID

      IF(mype == pe_io) THEN
         WRITE(*,*)'INTERSECTING GEOMETRY WITH V-MOMENTUM CELLS...'
      ENDIF
!======================================================================
!  Intersection between quadric Grid
!======================================================================

      intersect_x = .false.
      intersect_y = .false.
      intersect_z = .false.
      snap = .false.

      IF(.NOT.(use_stl.OR.use_msh)) THEN
         DO ijk = ijkstart3, ijkend3
!             IF(POTENTIAL_CUT_CELL_AT(IJK))  CALL INTERSECT(IJK,'V_MOMENTUM',Xn_V_int(IJK),Ye_V_int(IJK),Zt_V_int(IJK))
            CALL intersect(ijk,'V_MOMENTUM',xn_v_int(ijk),ye_v_int(ijk),zt_v_int(ijk))
         END DO

!======================================================================
!  Clean-up intersection flags in preparaton of small cells removal
!======================================================================
         DO ijk = ijkstart3, ijkend3
            IF(interior_cell_at(ijk)) THEN
!                IF(POTENTIAL_CUT_CELL_AT(IJK)) CALL CLEAN_INTERSECT(IJK,'V_MOMENTUM',Xn_V_int(IJK),Ye_V_int(IJK),Zt_V_int(IJK))
               CALL clean_intersect(ijk,'V_MOMENTUM',xn_v_int(ijk),ye_v_int(ijk),zt_v_int(ijk))
            ENDIF
         END DO

      ELSE
         CALL cad_intersect('V_MOMENTUM',xn_v_int,ye_v_int,zt_v_int)
      ENDIF


      number_of_new_v_points = 0

      DO ijk = ijkstart3, ijkend3
         CALL write_progress_bar(ijk,ijkend3 - ijkstart3 + 1,'C')

         IF(interior_cell_at(ijk)) THEN

!======================================================================
!  Get coordinates of eight nodes
!======================================================================

            CALL get_cell_node_coordinates(ijk,'V_MOMENTUM')

!======================================================================
!  Initialize location of velocity nodes at center of E,W,T faces
!======================================================================

            x_v_ec(ijk) = x_node(8)
            y_v_ec(ijk) = half * (y_node(6) + y_node(8))
            z_v_ec(ijk) = half * (z_node(4) + z_node(8))

            x_v_nc(ijk) = half * (x_node(7) + x_node(8))
            y_v_nc(ijk) = y_node(8)
            z_v_nc(ijk) = half * (z_node(4) + z_node(8))

            x_v_tc(ijk) = half * (x_node(7) + x_node(8))
            y_v_tc(ijk) = half * (y_node(6) + y_node(8))
            z_v_tc(ijk) = z_node(8)

!======================================================================
!  Get Connectivity
!======================================================================

            IF(no_k) THEN
               min_intersections = 2
               max_intersections = 2
               n_n1 = 5
               n_n2 = 8

            ELSE
               min_intersections = 3
               max_intersections = 6
               n_n1 = 1
               n_n2 = 8
            ENDIF


            CALL get_connectivity(ijk,'V_MOMENTUM',number_of_new_v_points,number_of_v_nodes(ijk),connectivity_v,&
            x_new_v_point,y_new_v_point,z_new_v_point,total_number_of_intersections,xn_v_int,ye_v_int,zt_v_int)


            IF(total_number_of_intersections < min_intersections ) THEN   ! Not a cut cell

               q_id = 1
               CALL eval_f('QUADRIC',x_node(0),y_node(0),z_node(0),q_id,f_node(0),clip_flag)

               CALL eval_f('POLYGON',x_node(0),y_node(0),z_node(0),n_polygon,f_node(0),clip_flag)

               CALL eval_f('USR_DEF',x_node(0),y_node(0),z_node(0),n_usr_def,f_node(0),clip_flag)

               CALL eval_stl_fct_at('V_MOMENTUM',ijk,0,f_node(0),clip_flag,bcid)

               IF(total_number_of_intersections==-1) f_node(0) = -one
               bc_v_id(ijk) = 0

               IF(f_node(0) <= zero) THEN
                  blocked_v_cell_at(ijk) = .false.
                  standard_v_cell_at(ijk) = .true.          ! Regular fluid cell

                  x_v_ec(ijk) = x_node(8)
                  y_v_ec(ijk) = half * (y_node(6) + y_node(8))
                  z_v_ec(ijk) = half * (z_node(4) + z_node(8))

                  x_v_nc(ijk) = half * (x_node(7) + x_node(8))
                  y_v_nc(ijk) = y_node(8)
                  z_v_nc(ijk) = half * (z_node(4) + z_node(8))

                  x_v_tc(ijk) = half * (x_node(7) + x_node(8))
                  y_v_tc(ijk) = half * (y_node(6) + y_node(8))
                  z_v_tc(ijk) = z_node(8)

               ELSE
                  blocked_v_cell_at(ijk) = .true.           ! Blocked fluid cell
                  standard_v_cell_at(ijk) = .false.
                  axy_v(ijk) = zero
                  axz_v(ijk) = zero
                  ayz_v(ijk) = zero
                  vol_v(ijk) = zero
               ENDIF

               IF(no_k) THEN
                  number_of_v_nodes(ijk) = 4
                  connectivity_v(ijk,1) = ijk_of_node(5)
                  connectivity_v(ijk,2) = ijk_of_node(6)
                  connectivity_v(ijk,3) = ijk_of_node(8)
                  connectivity_v(ijk,4) = ijk_of_node(7)
               ELSE
                  number_of_v_nodes(ijk) = 8
                  DO node = n_n1,n_n2
                     connectivity_v(ijk,node) = ijk_of_node(node)
                  END DO
               ENDIF

            ELSE IF(total_number_of_intersections > max_intersections ) THEN

               IF(print_warnings) THEN
                  WRITE(*,*) 'TOO MANY INTERSECTIONS FOUND IN V-CELL IJK = ',ijk
                  WRITE(*,*) 'MAXIMUM NUMBER OF INTERSECTIONS = ',max_intersections
                  WRITE(*,*) 'TOTAL NUMBER OF INTERSECTIONS = ',total_number_of_intersections
                  WRITE(*,*) 'REMOVING V-CELL FROM COMPUTATION.'
!                  WRITE(*,*) 'MFIX WILL EXIT NOW.'
!                  CALL MFIX_EXIT(MYPE)

               ENDIF

               blocked_v_cell_at(ijk) = .true.           ! Blocked fluid cell
               standard_v_cell_at(ijk) = .false.
               cut_v_cell_at(ijk) = .false.
               axy_v(ijk) = zero
               axz_v(ijk) = zero
               ayz_v(ijk) = zero
               vol_v(ijk) = zero

            ELSE                                         ! Cut cell

               cut_v_cell_at(ijk) = .true.
               blocked_v_cell_at(ijk) = .false.
               standard_v_cell_at(ijk) = .false.

               DO node = n_n1,n_n2
                  IF(f_node(node) < - tol_f) THEN
                     number_of_v_nodes(ijk) = number_of_v_nodes(ijk) + 1
                     connectivity_v(ijk,number_of_v_nodes(ijk)) = ijk_of_node(node)
                  ENDIF
               END DO

               CALL get_cut_cell_volume_and_areas(ijk,'V_MOMENTUM',number_of_v_nodes(ijk),connectivity_v,&
               x_new_v_point,y_new_v_point,z_new_v_point)

            ENDIF

         ENDIF
      END DO

      DO ijk = ijkstart3, ijkend3

         CALL get_cell_node_coordinates(ijk,'V_MOMENTUM')

         delx_ve(ijk) = x_node(8) - x_v(ijk)
         delx_vw(ijk) = x_v(ijk) - x_node(1)

         dely_vn(ijk) = y_node(8) - y_v(ijk)
         dely_vs(ijk) = y_v(ijk) - y_node(1)

         delz_vt(ijk) = z_node(8) - z_v(ijk)
         delz_vb(ijk) = z_v(ijk) - z_node(1)

      ENDDO

      RETURN

      END SUBROUTINE set_3d_cut_v_cell_flags

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: SET_3D_CUT_W_CELL_FLAGS                                C
!  Purpose: Set flags for W-Momentum cut cells, based on intersection  C
!  of the grid with the quadric(s)                                     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 set_3d_cut_w_cell_flags

      USE compar, ONLY: mype, pe_io, ijkstart3, ijkend3
      USE cut_cell_preproc, ONLY: cad_intersect, clean_intersect, clean_intersect_scalar, eval_f, intersect
      USE cutcell
      USE geometry, ONLY: axy_w, axz_w, ayz_w, vol_w
      USE param1, ONLY: half, one, zero
      USE polygon, ONLY: n_polygon
      USE quadric, ONLY: tol_f

      IMPLICIT NONE
      INTEGER :: IJK
      INTEGER :: TOTAL_NUMBER_OF_INTERSECTIONS
      INTEGER :: NODE,Q_ID
      LOGICAL :: CLIP_FLAG
      INTEGER :: BCID

      IF(mype == pe_io) THEN
         WRITE(*,10)'INTERSECTING GEOMETRY WITH W-MOMENTUM CELLS...'
      ENDIF
10    FORMAT(1x,a)
!======================================================================
!  Intersection between quadric Grid
!======================================================================

      intersect_x = .false.
      intersect_y = .false.
      intersect_z = .false.
      snap = .false.

      IF(.NOT.(use_stl.OR.use_msh)) THEN
         DO ijk = ijkstart3, ijkend3
!             IF(POTENTIAL_CUT_CELL_AT(IJK)) CALL INTERSECT(IJK,'W_MOMENTUM',Xn_W_int(IJK),Ye_W_int(IJK),Zt_W_int(IJK))
            CALL intersect(ijk,'W_MOMENTUM',xn_w_int(ijk),ye_w_int(ijk),zt_w_int(ijk))
         END DO

!======================================================================
!  Clean-up intersection flags in preparaton of small cells removal
!======================================================================
         DO ijk = ijkstart3, ijkend3
            IF(interior_cell_at(ijk)) THEN
!               IF(POTENTIAL_CUT_CELL_AT(IJK)) CALL CLEAN_INTERSECT(IJK,'W_MOMENTUM',Xn_W_int(IJK),Ye_W_int(IJK),Zt_W_int(IJK))
               CALL clean_intersect(ijk,'W_MOMENTUM',xn_w_int(ijk),ye_w_int(ijk),zt_w_int(ijk))
            ENDIF
         END DO

      ELSE
         CALL cad_intersect('W_MOMENTUM',xn_w_int,ye_w_int,zt_w_int)
      ENDIF


      number_of_new_w_points = 0

      DO ijk = ijkstart3, ijkend3

         CALL write_progress_bar(ijk,ijkend3 - ijkstart3 + 1,'C')

         IF(interior_cell_at(ijk)) THEN

!======================================================================
!  Get coordinates of eight nodes
!======================================================================

            CALL get_cell_node_coordinates(ijk,'W_MOMENTUM')

!======================================================================
!  Initialize location of velocity nodes at center of E,W,T faces
!======================================================================

            x_w_ec(ijk) = x_node(8)
            y_w_ec(ijk) = half * (y_node(6) + y_node(8))
            z_w_ec(ijk) = half * (z_node(4) + z_node(8))

            x_w_nc(ijk) = half * (x_node(7) + x_node(8))
            y_w_nc(ijk) = y_node(8)
            z_w_nc(ijk) = half * (z_node(4) + z_node(8))

            x_w_tc(ijk) = half * (x_node(7) + x_node(8))
            y_w_tc(ijk) = half * (y_node(6) + y_node(8))
            z_w_tc(ijk) = z_node(8)

!======================================================================
!  Get Connectivity
!======================================================================

            CALL get_connectivity(ijk,'W_MOMENTUM',number_of_new_w_points,number_of_w_nodes(ijk),connectivity_w,&
            x_new_w_point,y_new_w_point,z_new_w_point,total_number_of_intersections,xn_w_int,ye_w_int,zt_w_int)


            IF(total_number_of_intersections < 3 ) THEN   ! Not a cut cell

               q_id = 1
               CALL eval_f('QUADRIC',x_node(0),y_node(0),z_node(0),q_id,f_node(0),clip_flag)

               CALL eval_f('POLYGON',x_node(0),y_node(0),z_node(0),n_polygon,f_node(0),clip_flag)

               CALL eval_f('USR_DEF',x_node(0),y_node(0),z_node(0),n_usr_def,f_node(0),clip_flag)

               CALL eval_stl_fct_at('W_MOMENTUM',ijk,0,f_node(0),clip_flag,bcid)

               IF(total_number_of_intersections==-1) f_node(0) = -one
               bc_w_id(ijk) = 0

               IF(f_node(0) <= zero) THEN
                  blocked_w_cell_at(ijk) = .false.
                  standard_w_cell_at(ijk) = .true.          ! Regular fluid cell

                  x_w_ec(ijk) = x_node(8)
                  y_w_ec(ijk) = half * (y_node(6) + y_node(8))
                  z_w_ec(ijk) = half * (z_node(4) + z_node(8))

                  x_w_nc(ijk) = half * (x_node(7) + x_node(8))
                  y_w_nc(ijk) = y_node(8)
                  z_w_nc(ijk) = half * (z_node(4) + z_node(8))

                  x_w_tc(ijk) = half * (x_node(7) + x_node(8))
                  y_w_tc(ijk) = half * (y_node(6) + y_node(8))
                  z_w_tc(ijk) = z_node(8)

               ELSE
                  blocked_w_cell_at(ijk) = .true.           ! Blocked fluid cell
                  standard_w_cell_at(ijk) = .false.
                  axy_w(ijk) = zero
                  axz_w(ijk) = zero
                  ayz_w(ijk) = zero
                  vol_w(ijk) = zero
               ENDIF

               number_of_w_nodes(ijk) = 8
               DO node = 1,8
                  connectivity_w(ijk,node) = ijk_of_node(node)
               END DO

            ELSE IF(total_number_of_intersections > 6 ) THEN

               IF(print_warnings) THEN
                  WRITE(*,*) 'TOO MANY INTERSECTIONS FOUND IN W-CELL IJK = ',ijk
                  WRITE(*,*) 'MAXIMUM NUMBER OF INTERSECTIONS = ',6
                  WRITE(*,*) 'TOTAL NUMBER OF INTERSECTIONS = ',total_number_of_intersections
                  WRITE(*,*) 'REMOVING W-CELL FROM COMPUTATION.'
!                  WRITE(*,*) 'MFIX WILL EXIT NOW.'
!                  CALL MFIX_EXIT(MYPE)

               ENDIF

               blocked_w_cell_at(ijk) = .true.           ! Blocked fluid cell
               standard_w_cell_at(ijk) = .false.
               cut_w_cell_at(ijk) = .false.
               axy_w(ijk) = zero
               axz_w(ijk) = zero
               ayz_w(ijk) = zero
               vol_w(ijk) = zero

            ELSE                                         ! Cut cell

               cut_w_cell_at(ijk) = .true.

               DO node = 1,8
                  IF(f_node(node) < - tol_f) THEN
                     number_of_w_nodes(ijk) = number_of_w_nodes(ijk) + 1
                     connectivity_w(ijk,number_of_w_nodes(ijk)) = ijk_of_node(node)
                  ENDIF
               END DO


               CALL get_cut_cell_volume_and_areas(ijk,'W_MOMENTUM',number_of_w_nodes(ijk),connectivity_w,&
               x_new_w_point,y_new_w_point,z_new_w_point)

            ENDIF

         ENDIF
      END DO

      DO ijk = ijkstart3, ijkend3

         CALL get_cell_node_coordinates(ijk,'W_MOMENTUM')

         delx_we(ijk) = x_node(8) - x_w(ijk)
         delx_ww(ijk) = x_w(ijk) - x_node(1)

         dely_wn(ijk) = y_node(8) - y_w(ijk)
         dely_ws(ijk) = y_w(ijk) - y_node(1)

         delz_wt(ijk) = z_node(8) - z_w(ijk)
         delz_wb(ijk) = z_w(ijk) - z_node(1)

      ENDDO

      RETURN

      END SUBROUTINE set_3d_cut_w_cell_flags

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: SET_3D_CUT_CELL_TREATMENT_FLAGS                        C
!  Purpose: Set flags for scalar cut cells, based on intersection      C
!  of the grid with the quadric(s)                                     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 set_3d_cut_cell_treatment_flags

      USE compar, only: mype, pe_io, ijkstart3, ijkend3
      USE cutcell
      USE functions, ONLY: funijk
      USE geometry, ONLY: do_k
      USE indices, ONLY: i_of, j_of, k_of

      IMPLICIT NONE
      INTEGER :: IJK,I,J,K,IM,IP,JM,JP,KM,KP
      INTEGER :: IMJK,IPJK,IJMK,IJPK,IJKM,IJKP

      IF(mype == pe_io) THEN
         WRITE(*,*)'SETTING CUT CELL TREATMENT FLAGS...'
      ENDIF
!======================================================================
!  Set flags identifying cells requiring cut cell treatment:
!  These are the cut cells and their neighbours
!======================================================================

      cut_treatment_at   = .false.
      cut_u_treatment_at = .false.
      cut_v_treatment_at = .false.
      cut_w_treatment_at = .false.

      call send_receive_1d_logical(cut_cell_at,2)
      call send_receive_1d_logical(cut_u_cell_at,2)
      call send_receive_1d_logical(cut_v_cell_at,2)
      call send_receive_1d_logical(cut_w_cell_at,2)

      DO ijk = ijkstart3, ijkend3

         CALL write_progress_bar(ijk,ijkend3 - ijkstart3 + 1,'C')

         IF(interior_cell_at(ijk)) THEN

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

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

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

            imjk   = funijk(im,j,k)
            ijmk   = funijk(i,jm,k)
            ijkm   = funijk(i,j,km)

            ipjk   = funijk(ip,j,k)
            ijpk   = funijk(i,jp,k)
            ijkp   = funijk(i,j,kp)

            cut_treatment_at(ijk) = (cut_cell_at(ijk ).OR.   &
                                     cut_cell_at(imjk).OR.   &
                                     cut_cell_at(ipjk).OR.   &
                                     cut_cell_at(ijmk).OR.   &
                                     cut_cell_at(ijpk))

            cut_u_treatment_at(ijk) = (cut_u_cell_at(ijk ).OR.   &
                                       cut_u_cell_at(imjk).OR.   &
                                       cut_u_cell_at(ipjk).OR.   &
                                       cut_u_cell_at(ijmk).OR.   &
                                       cut_u_cell_at(ijpk))

            cut_v_treatment_at(ijk) = (cut_v_cell_at(ijk ).OR.   &
                                       cut_v_cell_at(imjk).OR.   &
                                       cut_v_cell_at(ipjk).OR.   &
                                       cut_v_cell_at(ijmk).OR.   &
                                       cut_v_cell_at(ijpk))

            IF(do_k) THEN

               cut_treatment_at(ijk) = (cut_treatment_at(ijk).OR.   &
                                            cut_cell_at(ijkm).OR.   &
                                            cut_cell_at(ijkp))

               cut_u_treatment_at(ijk) = (cut_u_treatment_at(ijk).OR.   &
                                              cut_u_cell_at(ijkm).OR.   &
                                              cut_u_cell_at(ijkp))

               cut_v_treatment_at(ijk) = (cut_v_treatment_at(ijk).OR.   &
                                              cut_v_cell_at(ijkm).OR.   &
                                              cut_v_cell_at(ijkp))

               cut_w_treatment_at(ijk) = (cut_w_cell_at(ijk ).OR.   &
                                          cut_w_cell_at(imjk).OR.   &
                                          cut_w_cell_at(ipjk).OR.   &
                                          cut_w_cell_at(ijmk).OR.   &
                                          cut_w_cell_at(ijpk).OR.   &
                                          cut_w_cell_at(ijkm).OR.   &
                                          cut_w_cell_at(ijkp))

            ENDIF

         ENDIF

         IF(.NOT.re_indexing) THEN
            IF(blocked_cell_at(ijk)) cut_treatment_at(ijk) = .false.
            IF(blocked_u_cell_at(ijk)) cut_u_treatment_at(ijk) = .false.
            IF(blocked_v_cell_at(ijk)) cut_v_treatment_at(ijk) = .false.
            IF(blocked_w_cell_at(ijk)) cut_w_treatment_at(ijk) = .false.
         ENDIF

      END DO

      IF(cg_safe_mode(1)==1) cut_treatment_at   = .false.
      IF(cg_safe_mode(3)==1) cut_u_treatment_at = .false.
      IF(cg_safe_mode(4)==1) cut_v_treatment_at = .false.
      IF(cg_safe_mode(5)==1) cut_w_treatment_at = .false.

      RETURN

      END SUBROUTINE set_3d_cut_cell_treatment_flags

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: SET_GHOST_CELL_FLAGS                                   C
!  Purpose: Set flags for ghost cell flags                             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 set_ghost_cell_flags

      USE compar, ONLY: ijkstart3, ijkend3
      USE compar, ONLY: iend1, iend2, iend3, jstart1, jstart2, jstart3, kstart1, kstart2, kstart3
      USE compar, only: istart1, istart2, istart3, jend1, jend2, jend3, kend1, kend2, kend3
      USE cutcell
      USE functions, ONLY: funijk
      USE functions, ONLY: west_of, east_of, south_of, north_of, bottom_of, top_of
      USE functions, ONLY: im_of, ip_of, jm_of, jp_of, km_of, kp_of
      USE geometry, ONLY: imax1, imin1, jmax1, jmin1, kmax1, kmin1, vol, vol_u, vol_v, vol_w, axy, axz, ayz, ayz_u, ayz_v, ayz_w, axy_u, axy_v, axy_w, axz_u, axz_v, axz_w, flag, do_k, flag_e, flag_n, flag_t

      USE bc
      USE sendrecv

      IMPLICIT NONE
      INTEGER :: IJK,I,J,K,I23,J23,K23
      INTEGER :: BCID
      INTEGER :: IPJK,IMJK,IJPK,IJMK,IJKP,IJKM

!     EAST BOUNDARY
      i = iend1

      IF(i==imax1) THEN

         DO i23 = iend2,iend3
            DO k = kstart3,kend3
               DO j = jstart3, jend3

                  ijk  = funijk(i,j,k)
                  ipjk = funijk(i23,j,k)

                  blocked_cell_at(ipjk) = blocked_cell_at(ijk)
                  blocked_u_cell_at(ipjk) = blocked_u_cell_at(ijk)
                  blocked_v_cell_at(ipjk) = blocked_v_cell_at(ijk)
                  blocked_w_cell_at(ipjk) = blocked_w_cell_at(ijk)

                  IF(blocked_cell_at(ipjk)) flag(ipjk) = 100

                  vol(ipjk)   = vol(ijk)
                  vol_u(ipjk) = vol_u(ijk)
                  vol_v(ipjk) = vol_v(ijk)
                  vol_w(ipjk) = vol_w(ijk)

                  ayz(ipjk)   = ayz(ijk)
                  ayz_u(ipjk) = ayz_u(ijk)
                  ayz_v(ipjk) = ayz_v(ijk)
                  ayz_w(ipjk) = ayz_w(ijk)

                  axy(ipjk)   = axy(ijk)
                  axy_u(ipjk) = axy_u(ijk)
                  axy_v(ipjk) = axy_v(ijk)
                  axy_w(ipjk) = axy_w(ijk)

                  axz(ipjk)   = axz(ijk)
                  axz_u(ipjk) = axz_u(ijk)
                  axz_v(ipjk) = axz_v(ijk)
                  axz_w(ipjk) = axz_w(ijk)

                  bc_id(ipjk)   = bc_id(ijk)
                  bc_u_id(ipjk) = bc_u_id(ijk)
                  bc_v_id(ipjk) = bc_v_id(ijk)
                  bc_w_id(ipjk) = bc_w_id(ijk)

               END DO
            END DO
         END DO

      ENDIF

!     WEST BOUNDARY
      i = istart1

      IF(i==imin1) THEN

         DO i23 = istart3,istart2
            DO k = kstart3,kend3
               DO j = jstart3, jend3

                  ijk  = funijk(i,j,k)
                  imjk = funijk(i23,j,k)

                  blocked_cell_at(imjk) = blocked_cell_at(ijk)
                  blocked_u_cell_at(imjk) = blocked_u_cell_at(ijk)
                  blocked_v_cell_at(imjk) = blocked_v_cell_at(ijk)
                  blocked_w_cell_at(imjk) = blocked_w_cell_at(ijk)

                  IF(blocked_cell_at(imjk)) flag(imjk) = 100

                  vol(imjk)   = vol(ijk)
                  vol_u(imjk) = vol_u(ijk)
                  vol_v(imjk) = vol_v(ijk)
                  vol_w(imjk) = vol_w(ijk)

                  ayz(imjk)   = ayz(ijk)
                  ayz_u(imjk) = ayz_u(ijk)
                  ayz_v(imjk) = ayz_v(ijk)
                  ayz_w(imjk) = ayz_w(ijk)

                  axy(imjk)   = axy(ijk)
                  axy_u(imjk) = axy_u(ijk)
                  axy_v(imjk) = axy_v(ijk)
                  axy_w(imjk) = axy_w(ijk)

                  axz(imjk)   = axz(ijk)
                  axz_u(imjk) = axz_u(ijk)
                  axz_v(imjk) = axz_v(ijk)
                  axz_w(imjk) = axz_w(ijk)

                  bc_id(imjk)   = bc_id(ijk)
                  bc_u_id(imjk) = bc_u_id(ijk)
                  bc_v_id(imjk) = bc_v_id(ijk)
                  bc_w_id(imjk) = bc_w_id(ijk)

               END DO
            END DO
         END DO

      ENDIF

!     NORTH BOUNDARY
      j = jend1

      IF(j==jmax1) THEN

         DO j23 = jend2,jend3
            DO k = kstart3,kend3
               DO i = istart3, iend3

                  ijk  = funijk(i,j,k)
                  ijpk = funijk(i,j23,k)

                  blocked_cell_at(ijpk) = blocked_cell_at(ijk)
                  blocked_u_cell_at(ijpk) = blocked_u_cell_at(ijk)
                  blocked_v_cell_at(ijpk) = blocked_v_cell_at(ijk)
                  blocked_w_cell_at(ijpk) = blocked_w_cell_at(ijk)

                  IF(blocked_cell_at(ijpk)) flag(ijpk) = 100

                  vol(ijpk)   = vol(ijk)
                  vol_u(ijpk) = vol_u(ijk)
                  vol_v(ijpk) = vol_v(ijk)
                  vol_w(ijpk) = vol_w(ijk)

                  ayz(ijpk)   = ayz(ijk)
                  ayz_u(ijpk) = ayz_u(ijk)
                  ayz_v(ijpk) = ayz_v(ijk)
                  ayz_w(ijpk) = ayz_w(ijk)

                  axy(ijpk)   = axy(ijk)
                  axy_u(ijpk) = axy_u(ijk)
                  axy_v(ijpk) = axy_v(ijk)
                  axy_w(ijpk) = axy_w(ijk)

                  axz(ijpk)   = axz(ijk)
                  axz_u(ijpk) = axz_u(ijk)
                  axz_v(ijpk) = axz_v(ijk)
                  axz_w(ijpk) = axz_w(ijk)

                  bc_id(ijpk)   = bc_id(ijk)
                  bc_u_id(ijpk) = bc_u_id(ijk)
                  bc_v_id(ijpk) = bc_v_id(ijk)
                  bc_w_id(ijpk) = bc_w_id(ijk)

               END DO
            END DO
         END DO

      ENDIF
!     SOUTH BOUNDARY
      j = jstart1

      IF(j==jmin1) THEN

         DO j23 = jstart3,jstart2
            DO k = kstart3,kend3
               DO i = istart3, iend3

                  ijk  = funijk(i,j,k)
                  ijmk = funijk(i,j23,k)

                  blocked_cell_at(ijmk) = blocked_cell_at(ijk)
                  blocked_u_cell_at(ijmk) = blocked_u_cell_at(ijk)
                  blocked_v_cell_at(ijmk) = blocked_v_cell_at(ijk)
                  blocked_w_cell_at(ijmk) = blocked_w_cell_at(ijk)

                  IF(blocked_cell_at(ijmk)) flag(ijmk) = 100

                  vol(ijmk)   = vol(ijk)
                  vol_u(ijmk) = vol_u(ijk)
                  vol_v(ijmk) = vol_v(ijk)
                  vol_w(ijmk) = vol_w(ijk)

                  ayz(ijmk)   = ayz(ijk)
                  ayz_u(ijmk) = ayz_u(ijk)
                  ayz_v(ijmk) = ayz_v(ijk)
                  ayz_w(ijmk) = ayz_w(ijk)

                  axy(ijmk)   = axy(ijk)
                  axy_u(ijmk) = axy_u(ijk)
                  axy_v(ijmk) = axy_v(ijk)
                  axy_w(ijmk) = axy_w(ijk)

                  axz(ijmk)   = axz(ijk)
                  axz_u(ijmk) = axz_u(ijk)
                  axz_v(ijmk) = axz_v(ijk)
                  axz_w(ijmk) = axz_w(ijk)

                  bc_id(ijmk)   = bc_id(ijk)
                  bc_u_id(ijmk) = bc_u_id(ijk)
                  bc_v_id(ijmk) = bc_v_id(ijk)
                  bc_w_id(ijmk) = bc_w_id(ijk)

               END DO
            END DO
         END DO

      ENDIF

      IF(do_k) THEN

!        TOP BOUNDARY
         k = kend1

         IF(k==kmax1) THEN

            DO k23=kend2,kend3
               DO j = jstart3,jend3
                  DO i = istart3, iend3

                     ijk  = funijk(i,j,k)
                     ijkp = funijk(i,j,k23)

                     blocked_cell_at(ijkp) = blocked_cell_at(ijk)
                     blocked_u_cell_at(ijkp) = blocked_u_cell_at(ijk)
                     blocked_v_cell_at(ijkp) = blocked_v_cell_at(ijk)
                     blocked_w_cell_at(ijkp) = blocked_w_cell_at(ijk)

                     IF(blocked_cell_at(ijkp)) flag(ijkp) = 100

                     vol(ijkp)   = vol(ijk)
                     vol_u(ijkp) = vol_u(ijk)
                     vol_v(ijkp) = vol_v(ijk)
                     vol_w(ijkp) = vol_w(ijk)

                     ayz(ijkp)   = ayz(ijk)
                     ayz_u(ijkp) = ayz_u(ijk)
                     ayz_v(ijkp) = ayz_v(ijk)
                     ayz_w(ijkp) = ayz_w(ijk)

                     axy(ijkp)   = axy(ijk)
                     axy_u(ijkp) = axy_u(ijk)
                     axy_v(ijkp) = axy_v(ijk)
                     axy_w(ijkp) = axy_w(ijk)

                     axz(ijkp)   = axz(ijk)
                     axz_u(ijkp) = axz_u(ijk)
                     axz_v(ijkp) = axz_v(ijk)
                     axz_w(ijkp) = axz_w(ijk)

                     bc_id(ijkp)   = bc_id(ijk)
                     bc_u_id(ijkp) = bc_u_id(ijk)
                     bc_v_id(ijkp) = bc_v_id(ijk)
                     bc_w_id(ijkp) = bc_w_id(ijk)

                  END DO
               END DO
            END DO

         ENDIF

!        BOTTOM BOUNDARY
         k = kstart1

         IF(k==kmin1) THEN

            DO k23 = kstart3,kstart2
               DO j = jstart3,jend3
                  DO i = istart3, iend3

                     ijk  = funijk(i,j,k)
                     ijkm = funijk(i,j,k23)

                     blocked_cell_at(ijkm) = blocked_cell_at(ijk)
                     blocked_u_cell_at(ijkm) = blocked_u_cell_at(ijk)
                     blocked_v_cell_at(ijkm) = blocked_v_cell_at(ijk)
                     blocked_w_cell_at(ijkm) = blocked_w_cell_at(ijk)

                     IF(blocked_cell_at(ijkm)) flag(ijkm) = 100

                     vol(ijkm)   = vol(ijk)
                     vol_u(ijkm) = vol_u(ijk)
                     vol_v(ijkm) = vol_v(ijk)
                     vol_w(ijkm) = vol_w(ijk)

                     ayz(ijkm)   = ayz(ijk)
                     ayz_u(ijkm) = ayz_u(ijk)
                     ayz_v(ijkm) = ayz_v(ijk)
                     ayz_w(ijkm) = ayz_w(ijk)

                     axy(ijkm)   = axy(ijk)
                     axy_u(ijkm) = axy_u(ijk)
                     axy_v(ijkm) = axy_v(ijk)
                     axy_w(ijkm) = axy_w(ijk)

                     axz(ijkm)   = axz(ijk)
                     axz_u(ijkm) = axz_u(ijk)
                     axz_v(ijkm) = axz_v(ijk)
                     axz_w(ijkm) = axz_w(ijk)

                     bc_id(ijkm)   = bc_id(ijk)
                     bc_u_id(ijkm) = bc_u_id(ijk)
                     bc_v_id(ijkm) = bc_v_id(ijk)
                     bc_w_id(ijkm) = bc_w_id(ijk)

                  END DO
               END DO
            END DO

         ENDIF

      ENDIF

      DO ijk = ijkstart3, ijkend3

         IF(blocked_u_cell_at(ijk)) flag_e(ijk)=100
         IF(blocked_v_cell_at(ijk)) flag_n(ijk)=100
         IF(blocked_w_cell_at(ijk)) flag_t(ijk)=100

      ENDDO

!     BLOCKED CELLS WERE ASSIGNED THE FLAG 100 DURING PRE_PROCESSING
!     THIS IS INCORRECT FOR FREEE-SLIP AND PARTIAL-SLIP
!     THE FLAG IS OVERWRITTEN HERE TO ACCOUNT FOR NSW,FSW AND PSW BCs

      DO ijk = ijkstart3, ijkend3

         IF(flag(ijk)==100) THEN

            imjk   = im_of(ijk)
            ijmk   = jm_of(ijk)
            ijkm   = km_of(ijk)

            ipjk   = ip_of(ijk)
            ijpk   = jp_of(ijk)
            ijkp   = kp_of(ijk)

            IF(cut_cell_at(imjk)) THEN
               bcid=bc_id(imjk)
               IF(bcid>0) THEN
                  IF(is_nsw(bc_type_enum(bcid))) flag(ijk) = 100
                  IF(is_fsw(bc_type_enum(bcid))) flag(ijk) = 101
                  IF(is_psw(bc_type_enum(bcid))) flag(ijk) = 102
               ENDIF

            ELSEIF(cut_cell_at(ipjk)) THEN
               bcid=bc_id(ipjk)
               IF(bcid>0) THEN
                  IF(is_nsw(bc_type_enum(bcid))) flag(ijk) = 100
                  IF(is_fsw(bc_type_enum(bcid))) flag(ijk) = 101
                  IF(is_psw(bc_type_enum(bcid))) flag(ijk) = 102
               ENDIF

            ELSEIF(cut_cell_at(ijmk)) THEN
               bcid=bc_id(ijmk)
               IF(bcid>0) THEN
                  IF(is_nsw(bc_type_enum(bcid))) flag(ijk) = 100
                  IF(is_fsw(bc_type_enum(bcid))) flag(ijk) = 101
                  IF(is_psw(bc_type_enum(bcid))) flag(ijk) = 102
               ENDIF

            ELSEIF(cut_cell_at(ijpk)) THEN
               bcid=bc_id(ijpk)
               IF(bcid>0) THEN
                  IF(is_nsw(bc_type_enum(bcid))) flag(ijk) = 100
                  IF(is_fsw(bc_type_enum(bcid))) flag(ijk) = 101
                  IF(is_psw(bc_type_enum(bcid))) flag(ijk) = 102
               ENDIF

            ELSEIF(cut_cell_at(ijkm)) THEN
               bcid=bc_id(ijkm)
               IF(bcid>0) THEN
                  IF(is_nsw(bc_type_enum(bcid))) flag(ijk) = 100
                  IF(is_fsw(bc_type_enum(bcid))) flag(ijk) = 101
                  IF(is_psw(bc_type_enum(bcid))) flag(ijk) = 102
               ENDIF

            ELSEIF(cut_cell_at(ijkp)) THEN
               bcid=bc_id(ijkp)
               IF(bcid>0) THEN
                  IF(is_nsw(bc_type_enum(bcid))) flag(ijk) = 100
                  IF(is_fsw(bc_type_enum(bcid))) flag(ijk) = 101
                  IF(is_psw(bc_type_enum(bcid))) flag(ijk) = 102
               ENDIF

            ENDIF
         ENDIF
      ENDDO

      call send_recv(flag,2)
      RETURN

      END SUBROUTINE set_ghost_cell_flags

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: GET_POTENTIAL_CUT_CELLS                                  C
!  Purpose: Set flags for scalar cut cells, based on intersection      C
!  of the grid with the quadric(s)                                     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_potential_cut_cells

      USE compar, only: ijkstart3, ijkend3, mype, pe_io
      USE cutcell
      USE functions, only: funijk, bottom_of, south_of, west_of
      USE geometry, ONLY: dx, dy, dz, do_k, imin3, imax3, jmin3, jmax3, kmin3, kmax3, flag, axy, axz, ayz, vol, no_k
      USE indices, only: i_of, j_of, k_of
      USE polygon, ONLY: n_polygon
      USE quadric, ONLY: tol_f
      USE param, ONLY: dimension_3
      USE param1, ONLY: half, zero

      IMPLICIT NONE
      INTEGER :: IJK,I,J,K,II,JJ,KK
      INTEGER :: I1,I2,J1,J2,K1,K2
      INTEGER :: NODE,N_N1,N_N2
      LOGICAL :: ALL_NEGATIVE,ALL_POSITIVE,CLIP_FLAG
      INTEGER :: Q_ID,BCID

      INTEGER :: IJK_NB,NUMBER_OF_POTENTIAL_CUT_CELLS

      DOUBLE PRECISION :: xc,yc,zc,fc

      LOGICAL, DIMENSION(DIMENSION_3) ::POSITIVE_F_AT

      potential_cut_cell_at=.true.

      RETURN  ! This subroutine is currently disabled

      IF(mype == pe_io) THEN
         WRITE(*,10)'ESTIMATING POTENTIAL SCALAR CUT CELLS...'
      ENDIF
10    FORMAT(1x,a)

!======================================================================
!  Evaluate f at cell center and store where f>0
!======================================================================

      DO ijk = ijkstart3, ijkend3

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

         xc = xg_e(i) - half*dx(i)
         yc = yg_n(j) - half*dy(j)

         IF(do_k) THEN
            zc = zg_t(k) - half*dz(k)
         ELSE
            zc = zero
         ENDIF

          q_id = 1
          CALL eval_f('QUADRIC',xc,yc,zc,q_id,fc,clip_flag)

          CALL eval_f('POLYGON',xc,yc,zc,n_polygon,fc,clip_flag)

          CALL eval_f('USR_DEF',xc,yc,zc,n_usr_def,fc,clip_flag)

          x_node(15) = xc
          y_node(15) = yc
          z_node(15) = zc
          CALL eval_stl_fct_at('SCALAR',ijk,15,fc,clip_flag,bcid)

          IF(fc>tol_f) THEN
             positive_f_at(ijk)=.true.
          ELSE
             positive_f_at(ijk)=.false.
          ENDIF

       ENDDO

      DO ijk = ijkstart3, ijkend3

         IF(interior_cell_at(ijk)) THEN

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

            i1 = max(i - 2,imin3)
            i2 = min(i + 2,imax3)
            j1 = max(j - 2,jmin3)
            j2 = min(j + 2,jmax3)

            IF(do_k) THEN
               k1 = max(k - 2,kmin3)
               k2 = min(k + 2,kmax3)
            ELSE
               k1=k
               k2=k
            ENDIF

            IF(positive_f_at(ijk)) THEN
               all_positive=.true.
               DO kk=k1,k2
                  DO jj=j1,j2
                     DO ii=i1,i2
                        ijk_nb = funijk(ii,jj,kk)
                        IF(.NOT.positive_f_at(ijk_nb)) THEN
                           all_positive=.false.
                        ENDIF
                     ENDDO
                  ENDDO
                ENDDO

                IF(all_positive) THEN
                   potential_cut_cell_at(ijk)=.false.
                   flag(ijk) = 100
                   blocked_cell_at(ijk) = .true.               ! Blocked fluid cell
                   standard_cell_at(ijk) = .false.
                   axy(ijk) = zero
                   axz(ijk) = zero
                   ayz(ijk) = zero
                   vol(ijk) = zero

                   axy(bottom_of(ijk)) = zero
                   axz(south_of(ijk)) = zero
                   ayz(west_of(ijk)) = zero
                ENDIF

            ELSE
               all_negative=.true.
               DO kk=k1,k2
                  DO jj=j1,j2
                     DO ii=i1,i2
                        ijk_nb = funijk(ii,jj,kk)
                        IF(positive_f_at(ijk_nb)) THEN
                           all_negative=.false.
                        ENDIF
                     ENDDO
                  ENDDO
                ENDDO

                IF(all_negative) THEN
                   potential_cut_cell_at(ijk)=.false.
                   blocked_cell_at(ijk) = .false.
                   standard_cell_at(ijk) = .true.           ! Regular fluid cell
                ENDIF

            ENDIF

            IF((flag(ijk)>=100).AND.(flag(ijk)<=102)) THEN
               potential_cut_cell_at(ijk)=.false.
               blocked_cell_at(ijk) = .true.
               standard_cell_at(ijk) = .false.          ! Blocked cell = wall cell
            ENDIF

         ENDIF

      END DO

      number_of_potential_cut_cells = 0

      IF(no_k) THEN
         n_n1 = 5
         n_n2 = 8
      ELSE
         n_n1 = 1
         n_n2 = 8
      ENDIF

      DO ijk=ijkstart3,ijkend3
         IF(potential_cut_cell_at(ijk)) THEN
            number_of_potential_cut_cells = number_of_potential_cut_cells + 1
         ELSE
            CALL get_cell_node_coordinates(ijk,'SCALAR')

            IF(no_k) THEN
               number_of_nodes(ijk) = 4
               connectivity(ijk,1) = ijk_of_node(5)
               connectivity(ijk,2) = ijk_of_node(6)
               connectivity(ijk,3) = ijk_of_node(8)
               connectivity(ijk,4) = ijk_of_node(7)
            ELSE
               number_of_nodes(ijk) = 8
               DO node = n_n1,n_n2
                  connectivity(ijk,node) = ijk_of_node(node)
               END DO
            ENDIF

            x_u(ijk) = x_node(8)
            y_u(ijk) = half * (y_node(6) + y_node(8))
            z_u(ijk) = half * (z_node(4) + z_node(8))

            x_v(ijk) = half * (x_node(7) + x_node(8))
            y_v(ijk) = y_node(8)
            z_v(ijk) = half * (z_node(4) + z_node(8))

            x_w(ijk) = half * (x_node(7) + x_node(8))
            y_w(ijk) = half * (y_node(6) + y_node(8))
            z_w(ijk) = z_node(8)

         ENDIF
      ENDDO

!      call SEND_RECEIVE_1D_LOGICAL(SNAP,2)
      IF(mype == pe_io) THEN
         WRITE(*,*)'DONE ESTIMATING POTENTIAL SCALAR CUT CELLS.',number_of_potential_cut_cells,ijkend3
      ENDIF

      RETURN
      END SUBROUTINE get_potential_cut_cells