get_delh.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name:  GET_DEL_H                                             C
!  Purpose: Finds the normal distance and unit vector from a cut face  C
!  to any point (x0,y0,z0)                                             C
!  The unit normal vector points away from the boundary,               C
!  towards the fluid.                                                  C
!  This subroutine must be called from a cut-cell                      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_del_h(IJK,TYPE_OF_CELL,X0,Y0,Z0,Del_H,Nx,Ny,Nz)

      USE param
      USE param1
      USE parallel
      USE constant
      USE run
      USE toleranc
      USE geometry
      USE indices
      USE compar
      USE sendrecv
      USE quadric
      USE cutcell
      USE functions

      IMPLICIT NONE
      CHARACTER (LEN=*) :: TYPE_OF_CELL
      DOUBLE PRECISION:: X0,Y0,Z0,XREF,YREF,ZREF
      INTEGER :: IJK,I,J,K
      DOUBLE PRECISION :: Del_H,Diagonal
      DOUBLE PRECISION :: Nx,Ny,Nz

      SELECT CASE (type_of_cell)
         CASE('SCALAR')

            IF(.NOT.cut_cell_at(ijk)) THEN
               WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H:'
               WRITE(*,*)' SCALAR CELL',ijk,' IS NOT A CUT CELL'
               WRITE(*,*)' MFiX will exit now.'
               CALL mfix_exit(mype)
            ENDIF

            nx = normal_s(ijk,1)
            ny = normal_s(ijk,2)
            nz = normal_s(ijk,3)

            xref = refp_s(ijk,1)
            yref = refp_s(ijk,2)
            zref = refp_s(ijk,3)

         CASE('U_MOMENTUM')

            IF(.NOT.cut_u_cell_at(ijk)) THEN
               WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H:'
               WRITE(*,*)' U-MOMENTUM CELL',ijk,' IS NOT A CUT CELL'
               WRITE(*,*)' MFiX will exit now.'
               CALL mfix_exit(mype)
            ENDIF

            nx = normal_u(ijk,1)
            ny = normal_u(ijk,2)
            nz = normal_u(ijk,3)

            xref = refp_u(ijk,1)
            yref = refp_u(ijk,2)
            zref = refp_u(ijk,3)

            IF(wall_u_at(ijk)) THEN
               nx = zero
               ny = zero
               nz = zero
               del_h = undefined
               RETURN
            ENDIF

         CASE('V_MOMENTUM')

            IF(.NOT.cut_v_cell_at(ijk)) THEN
               WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H:'
               WRITE(*,*)' V-MOMENTUM CELL',ijk,' IS NOT A CUT CELL'
               WRITE(*,*)' MFiX will exit now.'
               CALL mfix_exit(mype)
            ENDIF

            nx = normal_v(ijk,1)
            ny = normal_v(ijk,2)
            nz = normal_v(ijk,3)

            xref = refp_v(ijk,1)
            yref = refp_v(ijk,2)
            zref = refp_v(ijk,3)

            IF(wall_v_at(ijk)) THEN
               nx = zero
               ny = zero
               nz = zero
               del_h = undefined
               RETURN
            ENDIF

         CASE('W_MOMENTUM')

            IF(.NOT.cut_w_cell_at(ijk)) THEN
               WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H:'
               WRITE(*,*)' W-MOMENTUM CELL',ijk,' IS NOT A CUT CELL'
               WRITE(*,*)' MFiX will exit now.'
               CALL mfix_exit(mype)
            ENDIF

            nx = normal_w(ijk,1)
            ny = normal_w(ijk,2)
            nz = normal_w(ijk,3)

            xref = refp_w(ijk,1)
            yref = refp_w(ijk,2)
            zref = refp_w(ijk,3)

            IF(wall_w_at(ijk)) THEN
               nx = zero
               ny = zero
               nz = zero
               del_h = undefined
               RETURN
            ENDIF

         CASE DEFAULT
            WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H:'
            WRITE(*,*)'UNKNOWN TYPE OF CELL:',type_of_cell
            WRITE(*,*)'ACCEPTABLE TYPES ARE:'
            WRITE(*,*)'SCALAR'
            WRITE(*,*)'U_MOMENTUM'
            WRITE(*,*)'V_MOMENTUM'
            WRITE(*,*)'W_MOMENTUM'
            CALL mfix_exit(mype)
      END SELECT


      del_h = nx * (x0 - xref) + ny * (y0 - yref) + nz * (z0 - zref)

!======================================================================
! Negative values of DEL_H are not accepted
!======================================================================

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

       IF(no_k) THEN
          diagonal = sqrt(dx(i)**2 + dy(j)**2 )
       ELSE
          diagonal = sqrt(dx(i)**2 + dy(j)**2 + dz(k)**2)
       ENDIF

      IF (del_h <= tol_delh * diagonal) THEN

         del_h = undefined
         nx = zero
         ny = zero
         nz = zero

      ENDIF

      RETURN


      END SUBROUTINE get_del_h

  SUBROUTINE get_del_h_des(IJK,TYPE_OF_CELL,X0,Y0,Z0,Del_H,Nx,Ny,Nz, allow_neg_dist)

      USE param
      USE param1
      USE parallel
      USE constant
      USE run
      USE toleranc
      USE geometry
      USE indices
      USE compar
      USE sendrecv
      USE quadric
      USE cutcell
      USE functions

      IMPLICIT NONE
      CHARACTER (LEN=*) :: TYPE_OF_CELL
      DOUBLE PRECISION:: X0,Y0,Z0,XREF,YREF,ZREF
      INTEGER :: IJK,I,J,K
      DOUBLE PRECISION :: Del_H,Diagonal
      DOUBLE PRECISION :: Nx,Ny,Nz

      LOGICAL :: ALLOW_NEG_DIST

      SELECT CASE (type_of_cell)
         CASE('SCALAR')

            IF(.NOT.cut_cell_at(ijk)) THEN
               WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H_DES:'
               WRITE(*,*)' SCALAR CELL',ijk,' IS NOT A CUT CELL'
               WRITE(*,*)' I, J, K =',i_of(ijk), j_of(ijk), k_of(ijk)
               WRITE(*,*)' MFiX will exit now.'
               CALL mfix_exit(mype)
            ENDIF

            nx = normal_s(ijk,1)
            ny = normal_s(ijk,2)
            nz = normal_s(ijk,3)

            xref = refp_s(ijk,1)
            yref = refp_s(ijk,2)
            zref = refp_s(ijk,3)

         CASE('U_MOMENTUM')

            IF(.NOT.cut_u_cell_at(ijk)) THEN
               WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H_DES:'
               WRITE(*,*)' U-MOMENTUM CELL',ijk,' IS NOT A CUT CELL'
               WRITE(*,*)' MFiX will exit now.'
               CALL mfix_exit(mype)
            ENDIF

            nx = normal_u(ijk,1)
            ny = normal_u(ijk,2)
            nz = normal_u(ijk,3)

            xref = refp_u(ijk,1)
            yref = refp_u(ijk,2)
            zref = refp_u(ijk,3)

            IF(wall_u_at(ijk)) THEN
               nx = zero
               ny = zero
               nz = zero
               del_h = undefined
               RETURN
            ENDIF

         CASE('V_MOMENTUM')

            IF(.NOT.cut_v_cell_at(ijk)) THEN
               WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H_DES:'
               WRITE(*,*)' V-MOMENTUM CELL',ijk,' IS NOT A CUT CELL'
               WRITE(*,*)' MFiX will exit now.'
               CALL mfix_exit(mype)
            ENDIF

            nx = normal_v(ijk,1)
            ny = normal_v(ijk,2)
            nz = normal_v(ijk,3)

            xref = refp_v(ijk,1)
            yref = refp_v(ijk,2)
            zref = refp_v(ijk,3)

            IF(wall_v_at(ijk)) THEN
               nx = zero
               ny = zero
               nz = zero
               del_h = undefined
               RETURN
            ENDIF

         CASE('W_MOMENTUM')

            IF(.NOT.cut_w_cell_at(ijk)) THEN
               WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H_DES:'
               WRITE(*,*)' W-MOMENTUM CELL',ijk,' IS NOT A CUT CELL'
               WRITE(*,*)' MFiX will exit now.'
               CALL mfix_exit(mype)
            ENDIF

            nx = normal_w(ijk,1)
            ny = normal_w(ijk,2)
            nz = normal_w(ijk,3)

            xref = refp_w(ijk,1)
            yref = refp_w(ijk,2)
            zref = refp_w(ijk,3)

            IF(wall_w_at(ijk)) THEN
               nx = zero
               ny = zero
               nz = zero
               del_h = undefined
               RETURN
            ENDIF

         CASE DEFAULT
            WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H_DES:'
            WRITE(*,*)'UNKNOWN TYPE OF CELL:',type_of_cell
            WRITE(*,*)'ACCEPTABLE TYPES ARE:'
            WRITE(*,*)'SCALAR'
            WRITE(*,*)'U_MOMENTUM'
            WRITE(*,*)'V_MOMENTUM'
            WRITE(*,*)'W_MOMENTUM'
            CALL mfix_exit(mype)
      END SELECT

      del_h = nx * (x0 - xref) + ny * (y0 - yref) + nz * (z0 - zref)

!======================================================================
! Negative values of DEL_H are not accepted
!======================================================================

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

       IF(.NOT.allow_neg_dist) THEN
          diagonal = sqrt(dx(i)**2 + dy(j)**2 + dz(k)**2)

          IF (del_h <= tol_delh * diagonal) THEN
             del_h = undefined
             nx = zero
             ny = zero
             nz = zero

          ENDIF
       ENDIF

      RETURN

      END SUBROUTINE get_del_h_des
!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: STORE_CUT_FACE_INFO                                    C
!  Purpose: Compute and store unit normal vector and reference point   C
!           Defining a cut face                                        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 store_cut_face_info(IJK,TYPE_OF_CELL,N_CUT_FACE_NODES,COORD_CUT_FACE_NODES,X_MEAN,Y_MEAN,Z_MEAN)

      USE param
      USE param1
      USE parallel
      USE constant
      USE run
      USE toleranc
      USE geometry
      USE indices
      USE compar
      USE sendrecv
      USE quadric
      USE cutcell
      USE functions

      IMPLICIT NONE
      CHARACTER (LEN=*) :: TYPE_OF_CELL
      INTEGER :: IJK
      INTEGER :: N_CUT_FACE_NODES
      DOUBLE PRECISION, DIMENSION(3,15) :: COORD_CUT_FACE_NODES
      DOUBLE PRECISION :: X_MEAN,Y_MEAN,Z_MEAN
      DOUBLE PRECISION, DIMENSION(3) :: NN,V,N1,N2
      DOUBLE PRECISION :: NORM,N1_dot_N2
      DOUBLE PRECISION, DIMENSION(3) :: VECTEMP,VECTEMP2

     IF(no_k) THEN  ! 2D case

        nn(1) = coord_cut_face_nodes(2,1) - coord_cut_face_nodes(2,2) ! y1-y2
        nn(2) = coord_cut_face_nodes(1,2) - coord_cut_face_nodes(1,1) ! x2-x1
        nn(3) = zero

     ELSE

!======================================================================
! Make sure there are a least three points along the plane
!======================================================================

        IF(n_cut_face_nodes < 3) THEN
           WRITE(*,*)' ERROR IN SUBROUTINE STORE_CUT_FACE_INFO:'
           WRITE(*,*)' CUT FACE HAS LESS THAN 3 NODES.'
           WRITE(*,*)' MFIX WILL EXIT NOW.'
           CALL mfix_exit(mype)
        END IF

!======================================================================
!  Find tentative unit normal vector
!  and reverse sign if necessary
!  (unit vector must be pointing towards the fluid)
!======================================================================

        vectemp  = coord_cut_face_nodes(:,2)-coord_cut_face_nodes(:,1)
        vectemp2 = coord_cut_face_nodes(:,3)-coord_cut_face_nodes(:,1)
        nn = cross_product(vectemp,vectemp2)

     ENDIF

      norm = sqrt(dot_product(nn(:),nn(:)))
      nn = nn / norm

      v(1) = x_mean - coord_cut_face_nodes(1,1)
      v(2) = y_mean - coord_cut_face_nodes(2,1)
      v(3) = z_mean - coord_cut_face_nodes(3,1)

      IF (dot_product(nn,v) < zero) nn = - nn

      IF(n_cut_face_nodes > 3) THEN     ! FOR 3D geometry, check normal of plane defined by nodes 1,2, and 4

         n1 = nn  ! Keep copy of previous N (nodes 1,2,3)

        vectemp  = coord_cut_face_nodes(:,2)-coord_cut_face_nodes(:,1)
        vectemp2 = coord_cut_face_nodes(:,4)-coord_cut_face_nodes(:,1)

         n2 = cross_product(vectemp,vectemp2)

         norm = sqrt(dot_product(n2(:),n2(:)))
         n2 = n2 / norm

         v(1) = x_mean - coord_cut_face_nodes(1,1)
         v(2) = y_mean - coord_cut_face_nodes(2,1)
         v(3) = z_mean - coord_cut_face_nodes(3,1)

         IF (dot_product(n2,v) < zero) n2 = - n2

      ENDIF

      n1_dot_n2 = dot_product(n1,n2)
      debug_cg(ijk,1)=n1_dot_n2

      IF(n1_dot_n2<0.99) THEN

!         What should be done when the unit vectors are different ?

      ENDIF

!======================================================================
! Store unit normal vector and reference point
!======================================================================

      SELECT CASE (type_of_cell)
         CASE('SCALAR')

            normal_s(ijk,:) = nn
            refp_s(ijk,:)   = coord_cut_face_nodes(:,1)

            IF(do_k) CALL test_del_h(ijk,'SCALAR') ! test for negative del_H

         CASE('U_MOMENTUM')

            normal_u(ijk,:) = nn
            refp_u(ijk,:)   = coord_cut_face_nodes(:,1)

         CASE('V_MOMENTUM')

            normal_v(ijk,:) = nn
            refp_v(ijk,:)   = coord_cut_face_nodes(:,1)

         CASE('W_MOMENTUM')

            normal_w(ijk,:) = nn
            refp_w(ijk,:)   = coord_cut_face_nodes(:,1)

         CASE DEFAULT
            WRITE(*,*)'SUBROUTINE: STORE_CUT_FACE_INFO'
            WRITE(*,*)'UNKNOWN TYPE OF CELL:',type_of_cell
            WRITE(*,*)'ACCEPTABLE TYPES ARE:'
            WRITE(*,*)'SCALAR'
            WRITE(*,*)'U_MOMENTUM'
            WRITE(*,*)'V_MOMENTUM'
            WRITE(*,*)'W_MOMENTUM'
            CALL mfix_exit(mype)
      END SELECT

      RETURN

      END SUBROUTINE store_cut_face_info

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: TEST_DEL_H                                             C
!  Purpose: tests the computation of wall distance                     C
!           If a negative distance is detected, the normal vector      C
!           is inverted                                                C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 22-Feb-12  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Revision Number #                                  Date: ##-###-##  C
!  Author: #                                                           C
!  Purpose: #                                                          C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
  SUBROUTINE test_del_h(IJK,TYPE_OF_CELL)

      USE param
      USE param1
      USE parallel
      USE constant
      USE run
      USE toleranc
      USE geometry
      USE indices
      USE compar
      USE sendrecv
      USE quadric
      USE cutcell

      IMPLICIT NONE
      CHARACTER (LEN=*) :: TYPE_OF_CELL
      DOUBLE PRECISION:: X_COPY,Y_COPY,Z_COPY
      INTEGER :: IJK
      INTEGER :: NODE,N_N1,N_N2,NN
      DOUBLE PRECISION :: Del_H
      DOUBLE PRECISION :: Nx,Ny,Nz

      LOGICAL :: ALLOW_NEG_DIST = .true.  ! forces GET_DEL_H_DES to output negative delh
                                           ! i.e. do not let the subroutine overwrite negative values

! This subroutine tests values of del_H for nodes defining a cut cell.
! Only nodes that are in the fluid region are tested.
! Nodes belonging to the cut face should return zero (or near zero) values and are not tested.
! If a negative del_H is detected, the unit normal vector is reversed.

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

      CALL get_cell_node_coordinates(ijk,type_of_cell)

      DO node = 1,number_of_nodes(ijk)
         IF(connectivity(ijk,node)<=ijkend3) THEN  ! node does not belong to the cut-face
                                                    ! i.e. is in the fluid region

            DO nn = n_n1,n_n2                     ! get node coordinate
               IF(connectivity(ijk,node) == ijk_of_node(nn)) THEN
                  x_copy = x_node(nn)
                  y_copy = y_node(nn)
                  z_copy = z_node(nn)
                  EXIT
               ENDIF
            ENDDO

!           Compute del_H
            CALL get_del_h_des(ijk,type_of_cell,x_copy,y_copy,z_copy,del_h,nx,ny,nz, allow_neg_dist)


            IF(del_h<zero) THEN


               IF(print_warnings.AND.mype==pe_io) THEN
                  WRITE(*,*) ' Warning: Negative delh detected in scalar cell :',ijk
                  WRITE(*,*) ' Location (X,Y,Z) = ',x_copy,y_copy,z_copy
                  WRITE(*,*) ' Reverting unit normal vector.'
               ENDIF

               normal_s(ijk,1) = -normal_s(ijk,1)
               normal_s(ijk,2) = -normal_s(ijk,2)
               normal_s(ijk,3) = -normal_s(ijk,3)

            ENDIF

         ENDIF
      ENDDO

      RETURN

      END SUBROUTINE test_del_h

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name:  GET_DISTANCE_TO_WALL                                  C
!  Purpose: Finds the distance fraom any scalar cell to the closest    C
!  wall                                                                C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 16-May-13  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Revision Number #                                  Date: ##-###-##  C
!  Author: #                                                           C
!  Purpose: #                                                          C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
  SUBROUTINE get_distance_to_wall

      USE param
      USE param1
      USE parallel
      USE constant
      USE run
      USE toleranc
      USE geometry
      USE indices
      USE compar
      USE sendrecv
      USE quadric
      USE cutcell
      USE functions

      USE mpi_utility

      IMPLICIT NONE
!      CHARACTER (LEN=*) :: TYPE_OF_CELL
      DOUBLE PRECISION:: X0,Y0,Z0,XREF,YREF,ZREF
      INTEGER :: IJK,NN

      DOUBLE PRECISION :: D_TO_CUT, D_TO_PE_REF

      INTEGER :: N_CUT_CELLS
      INTEGER :: LIST_OF_CUT_CELLS(dimension_3)

      INTEGER :: iproc,IERR,IJK_OFFSET,nb,n1,n2
      INTEGER :: GLOBAL_N_CUT_CELLS
      INTEGER, DIMENSION(0:numPEs-1) :: disp,rcount
      LOGICAL, DIMENSION(0:numPEs-1) :: ALREADY_VISITED
      DOUBLE PRECISION, DIMENSION(0:numPEs-1,3) :: PE_REFP,ALL_PE_REFP
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: LOCAL_REFP_S,GLOBAL_REFP_S

      IF(mype==pe_io) WRITE(*,*)'COMPUTING WALL DISTANCE...'

! Count the number of cut cells
      n_cut_cells = 0
      DO ijk = ijkstart3, ijkend3
         IF(interior_cell_at(ijk).AND.cut_cell_at(ijk)) THEN
            n_cut_cells = n_cut_cells + 1
            list_of_cut_cells(n_cut_cells) = ijk
         ENDIF
      ENDDO

! Store cut cell reference points (centroid of cut cell face)
! and save a reference point for the entire processor
      ALLOCATE (local_refp_s(n_cut_cells,3))

      n_cut_cells = 0
      pe_refp(mype,1) = zero
      pe_refp(mype,2) = zero
      pe_refp(mype,3) = zero

      DO ijk = ijkstart3, ijkend3
         IF(interior_cell_at(ijk).AND.cut_cell_at(ijk)) THEN
            n_cut_cells = n_cut_cells + 1
            local_refp_s(n_cut_cells,1) = refp_s(ijk,1)
            local_refp_s(n_cut_cells,2) = refp_s(ijk,2)
            local_refp_s(n_cut_cells,3) = refp_s(ijk,3)

            pe_refp(mype,1) = pe_refp(mype,1) + refp_s(ijk,1)
            pe_refp(mype,2) = pe_refp(mype,2) + refp_s(ijk,2)
            pe_refp(mype,3) = pe_refp(mype,3) + refp_s(ijk,3)
         ENDIF
      ENDDO

      IF(n_cut_cells>0) THEN
         pe_refp(mype,1) = pe_refp(mype,1) / n_cut_cells
         pe_refp(mype,2) = pe_refp(mype,2) / n_cut_cells
         pe_refp(mype,3) = pe_refp(mype,3) / n_cut_cells
      ELSE
         pe_refp(mype,1) = undefined
         pe_refp(mype,2) = undefined
         pe_refp(mype,3) = undefined
      ENDIF

!======================================================================
! Now, gather the local reference points to head node
! to get a global list, and broadcast it to each processor.
!======================================================================

!======================================================================
! First get the offset and build the rcount and disp arrays.
! rcount is the number of elements to be gathered.
! disp is the displacement of the variable size gather,
! i.e. the cumulative sum at a given procesor.
!======================================================================

      CALL allgather_1i (n_cut_cells,rcount,ierr)

      IF (mype == 0) THEN
         ijk_offset = 0
      ELSE
         ijk_offset = 0
         DO iproc=0,mype-1
            ijk_offset = ijk_offset + rcount(iproc)
         ENDDO
      ENDIF

      CALL allgather_1i (ijk_offset,disp,ierr)

      CALL allgather_1d (pe_refp(mype,1),all_pe_refp(:,1),ierr)
      CALL allgather_1d (pe_refp(mype,2),all_pe_refp(:,2),ierr)
      CALL allgather_1d (pe_refp(mype,3),all_pe_refp(:,3),ierr)

!======================================================================
! Get the global number of cut cells and allocate the
! global reference point array. Each processor gets its own
! copy of all cut cell reference points !!
!======================================================================
      CALL global_all_sum(n_cut_cells, global_n_cut_cells,  pe_io, ierr )
      ALLOCATE (global_refp_s(global_n_cut_cells,3))

!======================================================================
! For a serial run, the global and local arrays are the same.
! for a parallel run, first gather on head node, then broadcast to all.
!======================================================================

      IF(numpes==1) THEN  ! Serial run
         global_refp_s =  local_refp_s
      ELSE !Parallel run
         call gatherv_1d( local_refp_s(:,1), n_cut_cells, global_refp_s(:,1), rcount, disp, pe_io, ierr )
         call gatherv_1d( local_refp_s(:,2), n_cut_cells, global_refp_s(:,2), rcount, disp, pe_io, ierr )
         call gatherv_1d( local_refp_s(:,3), n_cut_cells, global_refp_s(:,3), rcount, disp, pe_io, ierr )

         call bcast(global_refp_s(:,1))
         call bcast(global_refp_s(:,2))
         call bcast(global_refp_s(:,3))
      ENDIF


      already_visited(:) = .false.
!======================================================================
!  Loop through all scalar cells, grouped by processor.
!  Compute the distance to each cut cell and take the minimum
!  This is doe in three passes, in order of likelyhood to find a wall:
!  1) Within local processor.
!  2) Within neighboring processors (from send2 schedule)
!  3) Within all other processors
!
!  The idea is that it is very likely that cut cells located
!  on current or neighbor processors will contribute to the
!  wall distance calculation, whereas cut cells located
!  on the other processors will not contribute
!  During the third pass, the entire processor is skipped
!  if a reference point (average of all cut faces reference
!  points) is farther than the current wall distance.
!  This should speed up the brute calculation of going
!  explicitely through all cut cells.
!  However, it is not guaranteed that dwall will be
!  independent of the grid partitionning. This is considered
!  very unlikely at the moment
!  Change DWALL_BRUTE_FORCE to .TRUE. to force brute-force
!  calculation.
!======================================================================

      DO ijk = ijkstart3, ijkend3

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

         IF(interior_cell_at(ijk)) THEN
            dwall(ijk) = undefined

            IF(.NOT.blocked_cell_at(ijk)) THEN

! Get coordinates of cell center
               CALL get_cell_node_coordinates(ijk,'SCALAR')

               x0 = x_node(0)
               y0 = y_node(0)
               z0 = z_node(0)

!======================================================================
! First pass: Loop through local cut cells
!======================================================================

               already_visited(mype) = .true.

               DO nn = 1,n_cut_cells

                  xref = local_refp_s(nn,1)
                  yref = local_refp_s(nn,2)
                  zref = local_refp_s(nn,3)

                  d_to_cut = sqrt((x0 - xref)**2 + (y0 - yref)**2 + (z0 - zref)**2)

                  dwall(ijk) = min(dwall(ijk),d_to_cut)

               ENDDO


!======================================================================
! Second pass: Loop through neighbor processors (use send2 schedule)
!======================================================================

               DO nb=1,nsend2
                  iproc = sendproc2(nb)
                  already_visited(iproc) = .true.
                  n1 = disp(iproc)+1
                  n2 = n1 + rcount(iproc) - 1

                  DO nn = n1,n2

                     xref = global_refp_s(nn,1)
                     yref = global_refp_s(nn,2)
                     zref = global_refp_s(nn,3)

                     d_to_cut = sqrt((x0 - xref)**2 + (y0 - yref)**2 + (z0 - zref)**2)

                     dwall(ijk) = min(dwall(ijk),d_to_cut)

                  ENDDO

               ENDDO

!======================================================================
! Third pass: Loop through all other processors
! skip already visited processors (myPE and its neigbhors)
! First test if the average reference point is farther than dwall
! if this is true, then skip the entire processor cut cells.
!======================================================================


               DO iproc=0,numpes-1
                  IF(already_visited(iproc)) cycle

                  xref = all_pe_refp(iproc,1)
                  yref = all_pe_refp(iproc,2)
                  zref = all_pe_refp(iproc,3)

                  d_to_pe_ref = sqrt((x0 - xref)**2 + (y0 - yref)**2 + (z0 - zref)**2)

                  IF((dwall(ijk) < d_to_pe_ref).AND. &
                     (.NOT.dwall_brute_force  )) THEN
                     cycle
                  ELSE
                     n1 = disp(iproc)+1
                     n2 = n1 + rcount(iproc) - 1

                     DO nn = n1,n2

                        xref = global_refp_s(nn,1)
                        yref = global_refp_s(nn,2)
                        zref = global_refp_s(nn,3)

                        d_to_cut = sqrt((x0 - xref)**2 + (y0 - yref)**2 + (z0 - zref)**2)

                        dwall(ijk) = min(dwall(ijk),d_to_cut)

                     ENDDO
                  ENDIF
               ENDDO

            ENDIF

         ENDIF

      ENDDO

      call send_recv(dwall,2)
!======================================================================
!  Deallocate arrays before leaving
!======================================================================
      DEALLOCATE (local_refp_s)
      DEALLOCATE (global_refp_s)
      RETURN

      END SUBROUTINE get_distance_to_wall