File: N:\mfix\model\cartesian_grid\get_cut_cell_flags.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: SET_3D_CUT_CELL_FLAGS                                  C
4     !  Purpose: Set flags for scalar cut cells, based on intersection      C
5     !  of the grid with the quadric(s)                                     C
6     !                                                                      C
7     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  Revision Number #                                  Date: ##-###-##  C
11     !  Author: #                                                           C
12     !  Purpose: #                                                          C
13     !                                                                      C
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
15       SUBROUTINE SET_3D_CUT_CELL_FLAGS
16     
17           USE compar, ONLY: mype, pe_io, ijkstart3, ijkend3, istart1, iend1, jstart1, jend1, kstart1, kend1
18           USE cutcell
19           USE functions, ONLY: funijk, ip_of, jp_of, kp_of, bottom_of, south_of, west_of, fluid_at, is_on_mype_wobnd
20           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
21           USE indices, ONLY: i_of, j_of, k_of
22           USE param, only: dimension_3
23           USE param1, only: zero, half, one, undefined
24           USE polygon, ONLY: n_polygon
25           USE quadric, ONLY: tol_f
26           USE sendrecv
27           USE cut_cell_preproc, ONLY: cad_intersect, clean_intersect, clean_intersect_scalar, eval_f, intersect
28           USE vtk, ONLY: GLOBAL_VAR_ALLOCATED, GRID_INFO_PRINTED_ON_SCREEN
29     
30           IMPLICIT NONE
31           INTEGER :: IJK,I,J,K
32           INTEGER :: TOTAL_NUMBER_OF_INTERSECTIONS
33           INTEGER :: NODE,N_N1,N_N2,Q_ID
34           INTEGER :: MIN_INTERSECTIONS,MAX_INTERSECTIONS
35           LOGICAL :: CLIP_FLAG
36           INTEGER :: BCID
37     
38           INTEGER, DIMENSION(DIMENSION_3,15) :: OLD_CONNECTIVITY
39           DOUBLE PRECISION, DIMENSION(:), allocatable :: X_OLD_POINT,Y_OLD_POINT,Z_OLD_POINT
40     
41           INTEGER, DIMENSION(6) :: NB
42     
43           INTEGER :: IJK_NB,NODE_NB,NN,N_NB,NC
44     
45           DOUBLE PRECISION :: X1,Y1,Z1,X2,Y2,Z2,D,TOT_VOL_NODE,TOT_VOL_CELLS
46     
47           DOUBLE PRECISION, DIMENSION(:,:), allocatable :: SCALAR_NODE_XYZ_TEMP
48     
49           DOUBLE PRECISION :: DIST, NORM1, NORM2, NORM3,Diagonal
50           INTEGER :: IJK2, I1, I2, J1, J2, K1, K2, II, JJ, KK
51           LOGICAL :: COND_1, COND_2
52     
53           allocate(X_OLD_POINT(DIMENSION_MAX_CUT_CELL))
54           allocate(Y_OLD_POINT(DIMENSION_MAX_CUT_CELL))
55           allocate(Z_OLD_POINT(DIMENSION_MAX_CUT_CELL))
56           allocate(SCALAR_NODE_XYZ_TEMP(DIMENSION_3, 3))
57     
58           IF(MyPE == PE_IO) THEN
59              WRITE(*,10)'INTERSECTING GEOMETRY WITH SCALAR CELLS...'
60           ENDIF
61     10    FORMAT(1X,A)
62     
63           GLOBAL_VAR_ALLOCATED = .FALSE.
64           GRID_INFO_PRINTED_ON_SCREEN = .FALSE.
65     
66     !  Set arrays for computing indices
67     !
68     
69     !      CALL SET_INCREMENTS
70     
71           DX(IMAX3+1) = DX(IMAX3)
72           DY(JMAX3+1) = DY(JMAX3)
73           DZ(KMAX3+1) = DZ(KMAX3)
74     
75     !     x-Location of U-momentum cells for original (uncut grid)
76           IF (DO_I) THEN
77              XG_E(1) = ZERO
78              DO I = IMIN1, IMAX3
79                 XG_E(I) = XG_E(I-1) + DX(I)
80              END DO
81           ENDIF
82     
83     
84     !     y-Location of V-momentum cells for original (uncut grid)
85           IF (DO_J) THEN
86              YG_N(1) = ZERO
87              DO J = JMIN1, JMAX3
88                 YG_N(J) = YG_N(J-1) + DY(J)
89              END DO
90           ENDIF
91     
92     
93     !     z-Location of W-momentum cells for original (uncut grid)
94           IF (DO_K) THEN
95              ZG_T(1) = ZERO
96              DO K = KMIN1, KMAX3
97                 ZG_T(K) = ZG_T(K-1) + DZ(K)
98              END DO
99           ENDIF
100     
101           PARTITION = DBLE(myPE)       ! ASSIGN processor ID (for vizualisation)
102     
103           CALL SET_GEOMETRY1   ! INITIALIZE ALL VOLUMES AND AREAS
104     
105     
106     !======================================================================
107     !  Grid coordinates:
108     !======================================================================
109     
110           SMALL_CELL_AT = .FALSE.
111           SMALL_CELL_FLAG = 0
112           NUMBER_OF_SMALL_CELLS = 0
113     
114     !======================================================================
115     !  Intersection between quadric Grid
116     !======================================================================
117     
118           INTERSECT_X = .FALSE.
119           INTERSECT_Y = .FALSE.
120           INTERSECT_Z = .FALSE.
121           SNAP = .FALSE.
122     
123           DO IJK = IJKSTART3, IJKEND3
124     
125              I = I_OF(IJK)
126              J = J_OF(IJK)
127              K = K_OF(IJK)
128     
129              IF(NO_K) THEN   ! 2D case
130     
131                 INTERIOR_CELL_AT(IJK) = (     (I >= ISTART1 ).AND.(I <= IEND1 )  &
132                                          .AND.(J >= JSTART1 ).AND.(J <= JEND1 ) )
133     
134              ELSE            ! 3D case
135     
136                 INTERIOR_CELL_AT(IJK) = (     (I >= ISTART1 ).AND.(I <= IEND1 )  &
137                                          .AND.(J >= JSTART1 ).AND.(J <= JEND1 )  &
138                                          .AND.(K >= KSTART1 ).AND.(K <= KEND1 ) )
139     
140              ENDIF
141     
142           END DO
143     
144            NUMBER_OF_NODES = 0
145     
146           CALL GET_POTENTIAL_CUT_CELLS
147     
148     !       POTENTIAL_CUT_CELL_AT=.TRUE.
149     
150           IF(.NOT.(USE_STL.OR.USE_MSH)) THEN
151              DO IJK = IJKSTART3, IJKEND3
152                 IF(POTENTIAL_CUT_CELL_AT(IJK))  CALL INTERSECT(IJK,'SCALAR',Xn_int(IJK),Ye_int(IJK),Zt_int(IJK))
153     !            CALL INTERSECT(IJK,'SCALAR',Xn_int(IJK),Ye_int(IJK),Zt_int(IJK))
154              END DO
155     
156     !======================================================================
157     !  Clean-up intersection flags by snapping intersection points to close
158     !  corner (within tolerance TOL_SNAP)
159     !======================================================================
160              CALL CLEAN_INTERSECT_SCALAR
161     
162           ELSE
163              CALL CAD_INTERSECT('SCALAR',Xn_int,Ye_int,Zt_int)
164           ENDIF
165     !        NUMBER_OF_NODES = 0
166     
167           DO IJK = IJKSTART3, IJKEND3
168     
169               IF(POTENTIAL_CUT_CELL_AT(IJK))  THEN
170     
171              CALL WRITE_PROGRESS_BAR(IJK,IJKEND3 - IJKSTART3 + 1,'C')
172     
173     !======================================================================
174     !  Get coordinates of eight nodes
175     !======================================================================
176     
177              CALL GET_CELL_NODE_COORDINATES(IJK,'SCALAR')
178     
179              SCALAR_NODE_XYZ_TEMP(IJK,1) = X_NODE(8)
180              SCALAR_NODE_XYZ_TEMP(IJK,2) = Y_NODE(8)
181              SCALAR_NODE_XYZ_TEMP(IJK,3) = Z_NODE(8)
182     !======================================================================
183     !  Initialize location of velocity nodes
184     !======================================================================
185     
186              X_U(IJK) = X_NODE(8)
187              Y_U(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
188              Z_U(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
189     
190              X_V(IJK) = HALF * (X_NODE(7) + X_NODE(8))
191              Y_V(IJK) = Y_NODE(8)
192              Z_V(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
193     
194              X_W(IJK) = HALF * (X_NODE(7) + X_NODE(8))
195              Y_W(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
196              Z_W(IJK) = Z_NODE(8)
197     
198              IF(INTERIOR_CELL_AT(IJK)) THEN
199     !======================================================================
200     !  Get Connectivity
201     !======================================================================
202     
203                 IF(NO_K) THEN
204                    MIN_INTERSECTIONS = 2
205                    MAX_INTERSECTIONS = 2
206                    N_N1 = 5
207                    N_N2 = 8
208     
209                 ELSE
210                    MIN_INTERSECTIONS = 3
211                    MAX_INTERSECTIONS = 6
212                    N_N1 = 1
213                    N_N2 = 8
214                 ENDIF
215     
216                 CALL GET_CONNECTIVITY(IJK,'SCALAR',NUMBER_OF_NEW_POINTS,NUMBER_OF_NODES(IJK),CONNECTIVITY,&
217                 X_NEW_POINT,Y_NEW_POINT,Z_NEW_POINT,TOTAL_NUMBER_OF_INTERSECTIONS,Xn_int,Ye_int,Zt_int)
218     
219                 IF(TOTAL_NUMBER_OF_INTERSECTIONS < MIN_INTERSECTIONS ) THEN   ! Not a cut cell
220     
221                    Q_ID = 1
222                    CALL EVAL_F('QUADRIC',X_NODE(0),Y_NODE(0),Z_NODE(0),Q_ID,F_NODE(0),CLIP_FLAG)
223     
224                    CALL EVAL_F('POLYGON',X_NODE(0),Y_NODE(0),Z_NODE(0),N_POLYGON,F_NODE(0),CLIP_FLAG)
225     
226                    CALL EVAL_F('USR_DEF',X_NODE(0),Y_NODE(0),Z_NODE(0),N_USR_DEF,F_NODE(0),CLIP_FLAG)
227     
228     !               CALL EVAL_STL_FCT_AT('SCALAR',IJK,0,F_NODE(0),CLIP_FLAG,BCID)
229     
230                    IF(USE_MSH.OR.USE_STL) THEN
231     !                  CALL EVAL_STL_FCT_AT('SCALAR',IJK,1,F_NODE(1),CLIP_FLAG,BCID)
232     !                  CALL EVAL_STL_FCT_AT('SCALAR',IJK,8,F_NODE(8),CLIP_FLAG,BCID)
233     !                  F_NODE(0) = HALF*(F_NODE(1)+F_NODE(8))
234     ! Pick the first defined value of F_AT at one of the cell corners
235     ! This should avoid wrongly setting a cell next to the global ghost layer as blocked cell
236     ! when the average between undefined and negative is a positive number
237     ! This occurs for external flows only due to F_AT being initialized as UNDEFINED
238     ! which is a positive number (considered outside the fluid region).
239                       DO NODE=1,8
240                          CALL EVAL_STL_FCT_AT('SCALAR',IJK,NODE,F_NODE(NODE),CLIP_FLAG,BCID)
241                          IF(F_NODE(NODE)/=UNDEFINED.AND.F_NODE(NODE)/=ZERO) THEN
242                             F_NODE(0) = F_NODE(NODE)
243                             EXIT
244                          ENDIF
245     
246                       ENDDO
247                    ENDIF
248     
249     
250     !  Uncomment the line below to force the cell to be a fluid cell
251     !  Currently, F_NODE(0) > 0 ==> Blocked cell
252     !             F_NODE(0) < 0 ==> Fluid cell
253     
254     !               IF(TOTAL_NUMBER_OF_INTERSECTIONS==-1) F_NODE(0) = -ONE
255     
256                    BC_ID(IJK) = 0
257     
258     !               IF(F_NODE(0) < TOL_F) THEN
259                    IF(F_NODE(0) < ZERO) THEN
260                       IF((FLAG(IJK)>=100).AND.(FLAG(IJK)<=102)) THEN
261                          BLOCKED_CELL_AT(IJK) = .TRUE.
262                       ELSE
263                          BLOCKED_CELL_AT(IJK) = .FALSE.
264                          STANDARD_CELL_AT(IJK) = .TRUE.           ! Regular fluid cell
265                       ENDIF
266                    ELSE
267                       FLAG(IJK) = 100
268                       BLOCKED_CELL_AT(IJK) = .TRUE.               ! Blocked fluid cell
269                       STANDARD_CELL_AT(IJK) = .FALSE.
270                       AXY(IJK) = ZERO
271                       AXZ(IJK) = ZERO
272                       AYZ(IJK) = ZERO
273                       VOL(IJK) = ZERO
274     
275                       AXY(BOTTOM_OF(IJK)) = ZERO
276                       AXZ(SOUTH_OF(IJK)) = ZERO
277                       AYZ(WEST_OF(IJK)) = ZERO
278                    ENDIF
279     
280                    IF(NO_K) THEN
281                       NUMBER_OF_NODES(IJK) = 4
282                       CONNECTIVITY(IJK,1) = IJK_OF_NODE(5)
283                       CONNECTIVITY(IJK,2) = IJK_OF_NODE(6)
284                       CONNECTIVITY(IJK,3) = IJK_OF_NODE(8)
285                       CONNECTIVITY(IJK,4) = IJK_OF_NODE(7)
286                    ELSE
287                       NUMBER_OF_NODES(IJK) = 8
288                       DO NODE = N_N1,N_N2
289                          CONNECTIVITY(IJK,NODE) = IJK_OF_NODE(NODE)
290                       END DO
291                    ENDIF
292     
293     
294                 ELSE IF(TOTAL_NUMBER_OF_INTERSECTIONS > MAX_INTERSECTIONS ) THEN
295     
296                    IF(PRINT_WARNINGS) THEN
297     
298                       WRITE(*,*) 'TOO MANY INTERSECTIONS FOUND IN CELL IJK = ',IJK
299                       WRITE(*,*) 'MAXIMUM NUMBER OF INTERSECTIONS = ',MAX_INTERSECTIONS
300                       WRITE(*,*) 'TOTAL NUMBER OF INTERSECTIONS = ',TOTAL_NUMBER_OF_INTERSECTIONS
301                       WRITE(*,*) 'REMOVING SCALAR CELL FROM COMPUTATION.'
302     !                  WRITE(*,*) 'MFIX WILL EXIT NOW.'
303     !                  CALL MFIX_EXIT(MYPE)
304     
305                    ENDIF
306     
307                    BLOCKED_CELL_AT(IJK) = .TRUE.               ! Blocked fluid cell
308     
309                 ELSE                                         ! Cut cell
310     
311                    CUT_CELL_AT(IJK) = .TRUE.
312                    BLOCKED_CELL_AT(IJK) = .FALSE.
313                    STANDARD_CELL_AT(IJK) = .FALSE.
314     
315                    DO NODE = N_N1,N_N2
316                       IF(F_NODE(NODE) < - TOL_F) THEN
317                          IF(.NOT.SNAP(IJK_OF_NODE(NODE))) THEN
318                             NUMBER_OF_NODES(IJK) = NUMBER_OF_NODES(IJK) + 1
319                             CONNECTIVITY(IJK,NUMBER_OF_NODES(IJK)) = IJK_OF_NODE(NODE)
320                          ENDIF
321                       ENDIF
322                    END DO
323     
324                    CALL GET_CUT_CELL_VOLUME_AND_AREAS(IJK,'SCALAR',NUMBER_OF_NODES(IJK),CONNECTIVITY,&
325                    X_NEW_POINT,Y_NEW_POINT,Z_NEW_POINT)
326     
327                 ENDIF
328     
329              ENDIF        ! Interior cell
330     
331           ENDIF
332           END DO          ! IJK Loop
333     
334     
335           call SEND_RECEIVE_1D_LOGICAL(SMALL_CELL_AT,2)
336           call SEND_RECEIVE_1D_LOGICAL(STANDARD_CELL_AT,2)
337           call SEND_RECEIVE_1D_LOGICAL(BLOCKED_CELL_AT,2)
338           call SEND_RECEIVE_1D_LOGICAL(CUT_CELL_AT,2)
339     
340     ! Consolidate blocked cells
341           DO IJK = IJKSTART3, IJKEND3
342              IF(BLOCKED_CELL_AT(IJK)) THEN
343                 STANDARD_CELL_AT(IJK) = .FALSE.
344                 CUT_CELL_AT(IJK)      = .FALSE.
345     
346                 FLAG(IJK)             = 100
347     
348                 VOL(IJK)              = ZERO
349                 AXY(IJK)              = ZERO
350                 AXZ(IJK)              = ZERO
351                 AYZ(IJK)              = ZERO
352     
353                 AXY(BOTTOM_OF(IJK))   = ZERO
354                 AXZ(SOUTH_OF(IJK))    = ZERO
355                 AYZ(WEST_OF(IJK))     = ZERO
356              ENDIF
357           ENDDO
358     
359           call SEND_RECEIVE_1D_LOGICAL(SMALL_CELL_AT,2)
360           call SEND_RECEIVE_1D_LOGICAL(STANDARD_CELL_AT,2)
361           call SEND_RECEIVE_1D_LOGICAL(BLOCKED_CELL_AT,2)
362           call SEND_RECEIVE_1D_LOGICAL(CUT_CELL_AT,2)
363     
364           call send_recv(FLAG,2)
365           call send_recv(VOL,2)
366           call send_recv(AXY,2)
367           call send_recv(AXZ,2)
368           call send_recv(AYZ,2)
369     
370     !      print*,'Before removing duplicate nodes:'
371           DO IJK = IJKSTART3, IJKEND3,-1
372              IF(CUT_CELL_AT(IJK)) THEN
373                 print*,'==========================================================='
374                 print*,'IJK,  I,J=',IJK,I_OF(IJK),J_OF(IJK)
375                 print*,'NUMBER_OF_NODES=',NUMBER_OF_NODES(IJK)
376                 DO NODE = 1,NUMBER_OF_NODES(IJK)
377                    IF(CONNECTIVITY(IJK,NODE)>IJKEND3) THEN
378                       print*,'CNCT=',NODE,CONNECTIVITY(IJK,NODE), &
379                            X_NEW_POINT(CONNECTIVITY(IJK,NODE)-IJKEND3),Y_NEW_POINT(CONNECTIVITY(IJK,NODE)-IJKEND3)
380                    ELSE
381                       print*,'CNCT=',NODE,CONNECTIVITY(IJK,NODE)
382                    ENDIF
383                 ENDDO
384                 print*,''
385              ENDIF
386           ENDDO
387     
388     
389     
390           OLD_CONNECTIVITY = CONNECTIVITY
391     !      X_OLD_POINT      = X_NEW_POINT
392     !      Y_OLD_POINT      = Y_NEW_POINT
393     !      Z_OLD_POINT      = Z_NEW_POINT
394     
395     
396     !======================================================================
397     !  Removing duplicate new points
398     !  Up to here, a cut face node that is shared among neighbor cells had
399     !  a different index, and appeared as a duplicate node.
400     !  The list will be updated and duplicate nodes will be removed from
401     !  the connectivity list
402     !  This is done to ensure that we can assign a volume surrounding each
403     !  node, that will be used to compute solids volume fraction for MPPIC
404     !======================================================================
405     
406     !      print*,'Removing duplicate nodes...'
407     
408           DO IJK = IJKSTART3, IJKEND3
409              IF(CUT_CELL_AT(IJK)) THEN          ! for each cut cell, identify neibhor cells that are also cut cells
410                                                  ! Look east and north in 2D, and also Top in 3D
411     
412              Diagonal = DSQRT(DX(I_OF(IJK))**2 + DY(J_OF(IJK))**2 + DZ(K_OF(IJK))**2)
413     
414                 NB(1) = IP_OF(IJK)
415                 NB(2) = JP_OF(IJK)
416     
417                 IF(NO_K) THEN
418                    N_NB = 2
419                 ELSE
420                    N_NB = 3
421                    NB(3) = KP_OF(IJK)
422                 ENDIF
423     
424     
425                 DO NN = 1,N_NB
426                    IF(CUT_CELL_AT(NB(NN))) THEN   ! For two neighbor cut cells, compare each cut-face node to remove duplicates
427     
428                       IJK_NB = NB(NN)
429     
430     !                  print*,'comparing:',IJK,' and',IJK_NB
431     
432                       DO NODE = 1,NUMBER_OF_NODES(IJK)
433                          IF(OLD_CONNECTIVITY(IJK,NODE)>IJKEND3) THEN  ! node belongs to the cut-face
434     
435                             X1 = X_NEW_POINT(OLD_CONNECTIVITY(IJK,NODE)-IJKEND3)
436                             Y1 = Y_NEW_POINT(OLD_CONNECTIVITY(IJK,NODE)-IJKEND3)
437                             Z1 = Z_NEW_POINT(OLD_CONNECTIVITY(IJK,NODE)-IJKEND3)
438     
439     
440                             DO NODE_NB = 1,NUMBER_OF_NODES(IJK_NB)
441                                IF(OLD_CONNECTIVITY(IJK_NB,NODE_NB)>IJKEND3) THEN  ! node belongs to the cut-face
442     
443                                   X2 = X_NEW_POINT(OLD_CONNECTIVITY(IJK_NB,NODE_NB)-IJKEND3)
444                                   Y2 = Y_NEW_POINT(OLD_CONNECTIVITY(IJK_NB,NODE_NB)-IJKEND3)
445                                   Z2 = Z_NEW_POINT(OLD_CONNECTIVITY(IJK_NB,NODE_NB)-IJKEND3)
446     
447                                   ! compare coordinates of cut-face nodes
448                                   D = (X2-X1)**2 + (Y2-Y1)**2 + (Z2-Z1)**2
449     
450                                   ! Duplicate nodes have identical coordinates (within tolerance TOL_MERGE times diagonal)
451                                   IF(D<TOL_MERGE*Diagonal) THEN ! keep the smallest node ID
452                                      !  print*,'DULICATE NODES:', &
453                                      ! NODE,NODE_NB,OLD_CONNECTIVITY(IJK,NODE),OLD_CONNECTIVITY(IJK_NB,NODE_NB)
454                                      NC = MIN(OLD_CONNECTIVITY(IJK,NODE),OLD_CONNECTIVITY(IJK_NB,NODE_NB))
455                                      CONNECTIVITY(IJK   ,NODE   ) = NC
456                                      CONNECTIVITY(IJK_NB,NODE_NB) = NC
457     
458                                   ENDIF
459                                ENDIF
460                             ENDDO
461     
462                          ENDIF
463                       ENDDO
464     
465                    ENDIF
466                 ENDDO
467     
468              ENDIF
469           ENDDO
470     
471           ALLOCATE(SCALAR_NODE_XYZ(DIMENSION_3 + NUMBER_OF_NEW_POINTS,3))
472           ALLOCATE(Ovol_around_node(DIMENSION_3 + NUMBER_OF_NEW_POINTS))
473           ALLOCATE(SCALAR_NODE_ATWALL(DIMENSION_3 + NUMBER_OF_NEW_POINTS))
474     
475           !first fill with standard nodes
476           DO IJK = IJKSTART3, IJKEND3
477              SCALAR_NODE_XYZ(IJK,1:3)  = SCALAR_NODE_XYZ_TEMP(IJK,1:3)
478           ENDDO
479     
480           !now fill with cut-face nodes
481           DO IJK = 1,NUMBER_OF_NEW_POINTS
482              SCALAR_NODE_XYZ(IJKEND3+IJK,1) = X_NEW_POINT(IJK)
483              SCALAR_NODE_XYZ(IJKEND3+IJK,2) = Y_NEW_POINT(IJK)
484              SCALAR_NODE_XYZ(IJKEND3+IJK,3) = Z_NEW_POINT(IJK)
485           ENDDO
486     
487     
488     
489           SCALAR_NODE_ATWALL(:)  = .true.
490           !Rahul:
491           !One could either set all scalar_node_atwall to false and
492           !then set to true the nodes that are found as on the wall or outside
493           !the doamin
494           !In some situations, a node can be deemed to be both inside
495           !and outside the domain depending upon which cell you look from.
496           !Think of a small cell. According to the small cells, all the nodes
497           !are outside the domain, but from the perspective or surrounding
498           !cut-cells, some of the nodes of this small cell are within the domain.
499     
500           !so I'm setting all the points as being outside the domain.
501           !if a point is ever found to be in the domain, then it will stay away
502           !and further tests will not be able to revert it.
503           IJKLOOP: DO IJK = IJKSTART3, IJKEND3
504              I = I_OF(IJK)
505              J = J_OF(IJK)
506              K = K_OF(IJK)
507              IF(.not. IS_ON_myPE_wobnd(I,J,K)) cycle
508     
509              I1 = I-1
510              I2 = I
511              J1 = J-1
512              J2 = J
513     
514              IF(NO_K) THEN
515                 K1 = K
516                 K2 = K
517              ELSE
518                 K1 = K-1
519                 K2 = K
520              ENDIF
521              !Convention used to number node numbers is described below
522     
523              ! i=1, j=2           i=2, j=2
524              !   _____________________
525              !   |                   |
526              !   |  I = 2, J = 2     |
527              !   |___________________|
528              ! i=1, j=1           i=2, j=1
529     
530              !Let's say the scalar cell with I = 2 and J = 2, i.e., the
531              !first scalar cell in either direction.
532              !then this scalar cell's node numbering in x- direction
533              !will go from 1 to 2 and likewise in other directions.
534              DO KK = K1, K2
535                 DO JJ = J1, J2
536                    DO II = I1, I2
537                       IJK2 = funijk(II, JJ, KK)
538                       COND_1 = .false.
539                       COND_2 = .false.
540                       !if it was already found to be inside the domain, then don't bother
541                       IF(.not.SCALAR_NODE_ATWALL(IJK2))  CYCLE
542     
543                       IF(.not.FLUID_AT(IJK)) THEN
544                          !do nothing
545                       ELSE
546                          !IJK is a fluid scalar cell. Check if it is
547                          !a cut-cell or not
548                          IF(CUT_CELL_AT(IJK)) THEN
549                             CALL GET_DEL_H_DES(IJK,'SCALAR', &
550                             & SCALAR_NODE_XYZ(IJK2,1), SCALAR_NODE_XYZ(IJK2,2), &
551                             & SCALAR_NODE_XYZ(IJK2,3), &
552                             & DIST, NORM1, NORM2, NORM3, .true.)
553                             IF(DIST.GE.ZERO) THEN
554                                SCALAR_NODE_ATWALL(IJK2)  = .false.
555                                COND_1 =  .true.
556                             ENDIF
557                          ELSE
558                             SCALAR_NODE_ATWALL(IJK2)  = .false.
559                             COND_2 = .true.
560                          ENDIF
561                       ENDIF
562                       !if(II == 1) then
563                       ! write(*,'(10x, A, 4(2x,i10),3(2x,L2))') &
564                       ! 'I1,J1, I, J, ATWALL, COND1, COND2 =  ', II, JJ, I, J,  SCALAR_NODE_ATWALL(IJK2), COND_1, COND_2
565                       !endif
566                    ENDDO
567                 ENDDO
568              ENDDO
569     
570           ENDDO IJKLOOP
571     
572     !      print*,'After removing duplicate nodes:'
573           DO IJK = IJKSTART3, IJKEND3,-1
574              print*,'==========================================================='
575              print*,'IJK,  I,J=',IJK,I_OF(IJK),J_OF(IJK)
576              print*,'NUMBER_OF_NODES=',NUMBER_OF_NODES(IJK)
577              DO NODE = 1,NUMBER_OF_NODES(IJK)
578                 print*,'CNCT=',NODE,CONNECTIVITY(IJK,NODE), &
579                      SCALAR_NODE_XYZ(CONNECTIVITY(IJK,NODE),1),SCALAR_NODE_XYZ(CONNECTIVITY(IJK,NODE),2)
580              ENDDO
581              print*,''
582           ENDDO
583     
584     
585     
586           Ovol_around_node = ZERO !UNDEFINED
587     
588           TOT_VOL_CELLS = ZERO
589     
590           DO IJK = IJKSTART3, IJKEND3
591     
592              DO NODE = 1,NUMBER_OF_NODES(IJK)
593                 NC = CONNECTIVITY(IJK,NODE)
594                 Ovol_around_node(NC) = Ovol_around_node(NC) + VOL(IJK)/NUMBER_OF_NODES(IJK)
595              ENDDO
596     
597     !         print*,'IJK,VOL=',IJK,VOL(IJK)
598     
599     
600              IF(INTERIOR_CELL_AT(IJK)) TOT_VOL_CELLS = TOT_VOL_CELLS + VOL(IJK)
601     
602     
603           ENDDO
604     
605     
606           TOT_VOL_NODE= ZERO
607     
608           DO IJK = IJKSTART3, IJKEND3 + NUMBER_OF_NEW_POINTS    ! Loop over all nodes
609     
610              IF(Ovol_around_node(IJK)>ZERO) THEN
611     !             print*,'NODE,VOL=',IJK,Ovol_around_node(IJK)
612                  TOT_VOL_NODE = TOT_VOL_NODE + Ovol_around_node(IJK)
613                  Ovol_around_node(IJK) = ONE / Ovol_around_node(IJK)    ! Store One/volume
614              ENDIF
615     
616     
617           ENDDO
618     
619           IF(TOT_VOL_CELLS>ZERO) THEN
620              IF(DABS(ONE-DABS(TOT_VOL_NODE/TOT_VOL_CELLS))>1.0D-6) THEN
621                 WRITE(*,*)'==========================================================='
622                 WRITE(*,*)'Total volume around nodes=',TOT_VOL_NODE
623                 WRITE(*,*)'Total volume =',TOT_VOL_CELLS
624                 WRITE(*,*)'The two volumes above should be the same.'
625                 WRITE(*,*)'==========================================================='
626                 CALL MFIX_EXIT(MYPE)
627              ENDIF
628           ENDIF
629     
630           deallocate(X_OLD_POINT)
631           deallocate(Y_OLD_POINT)
632           deallocate(Z_OLD_POINT)
633           deallocate(SCALAR_NODE_XYZ_TEMP)
634     
635           RETURN
636           END SUBROUTINE SET_3D_CUT_CELL_FLAGS
637     
638     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
639     !                                                                      C
640     !  Module name: SET_3D_CUT_U_CELL_FLAGS                                C
641     !  Purpose: Set flags for U-Momentum cut cells, based on intersection  C
642     !  of the grid with the quadric(s)                                     C
643     !                                                                      C
644     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
645     !  Reviewer:                                          Date:            C
646     !                                                                      C
647     !  Revision Number #                                  Date: ##-###-##  C
648     !  Author: #                                                           C
649     !  Purpose: #                                                          C
650     !                                                                      C
651     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
652       SUBROUTINE SET_3D_CUT_U_CELL_FLAGS
653     
654           USE compar, ONLY: mype, pe_io, ijkstart3, ijkend3
655           USE cut_cell_preproc, ONLY: cad_intersect, clean_intersect, clean_intersect_scalar, eval_f, intersect
656           USE cutcell
657           USE geometry, ONLY: no_k, axy_u, ayz_u, vol_u, axz_u
658           USE param1, ONLY: half, one, zero
659           USE polygon, ONLY: n_polygon
660           USE quadric, ONLY: tol_f
661     
662           IMPLICIT NONE
663           INTEGER :: IJK
664           INTEGER :: TOTAL_NUMBER_OF_INTERSECTIONS
665           INTEGER :: NODE,N_N1,N_N2,Q_ID
666           INTEGER :: MIN_INTERSECTIONS,MAX_INTERSECTIONS
667           LOGICAL :: CLIP_FLAG
668           INTEGER :: BCID
669     
670           IF(MyPE == PE_IO) THEN
671              WRITE(*,10)'INTERSECTING GEOMETRY WITH U-MOMENTUM CELLS...'
672           ENDIF
673     10    FORMAT(1X,A)
674     !======================================================================
675     !  Intersection between quadric Grid
676     !======================================================================
677     
678           INTERSECT_X = .FALSE.
679           INTERSECT_Y = .FALSE.
680           INTERSECT_Z = .FALSE.
681           SNAP = .FALSE.
682           TOL_SNAP = ZERO
683     
684           IF(.NOT.(USE_STL.OR.USE_MSH)) THEN
685              DO IJK = IJKSTART3, IJKEND3
686     !             IF(POTENTIAL_CUT_CELL_AT(IJK))  CALL INTERSECT(IJK,'U_MOMENTUM',Xn_U_int(IJK),Ye_U_int(IJK),Zt_U_int(IJK))
687                 CALL INTERSECT(IJK,'U_MOMENTUM',Xn_U_int(IJK),Ye_U_int(IJK),Zt_U_int(IJK))
688              END DO
689     
690     !======================================================================
691     !  Clean-up intersection flags in preparaton of small cells removal
692     !======================================================================
693              DO IJK = IJKSTART3, IJKEND3
694                 IF(INTERIOR_CELL_AT(IJK)) THEN
695     !               IF(POTENTIAL_CUT_CELL_AT(IJK))  CALL CLEAN_INTERSECT(IJK,'U_MOMENTUM',Xn_U_int(IJK),Ye_U_int(IJK),Zt_U_int(IJK))
696                    CALL CLEAN_INTERSECT(IJK,'U_MOMENTUM',Xn_U_int(IJK),Ye_U_int(IJK),Zt_U_int(IJK))
697                 ENDIF
698              END DO
699     
700           ELSE
701              CALL CAD_INTERSECT('U_MOMENTUM',Xn_U_int,Ye_U_int,Zt_U_int)
702           ENDIF
703     
704     
705           NUMBER_OF_NEW_U_POINTS = 0
706     
707           DO IJK = IJKSTART3, IJKEND3
708     
709              CALL WRITE_PROGRESS_BAR(IJK,IJKEND3 - IJKSTART3 + 1,'C')
710     
711              IF(INTERIOR_CELL_AT(IJK)) THEN
712     
713     !======================================================================
714     !  Get coordinates of eight nodes
715     !======================================================================
716     
717                 CALL GET_CELL_NODE_COORDINATES(IJK,'U_MOMENTUM')
718     
719     !======================================================================
720     !  Initialize location of velocity nodes at center of E,W,T faces
721     !======================================================================
722     
723                       X_U_ec(IJK) = X_NODE(8)
724                       Y_U_ec(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
725                       Z_U_ec(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
726     
727                       X_U_nc(IJK) = HALF * (X_NODE(7) + X_NODE(8))
728                       Y_U_nc(IJK) = Y_NODE(8)
729                       Z_U_nc(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
730     
731                       X_U_tc(IJK) = HALF * (X_NODE(7) + X_NODE(8))
732                       Y_U_tc(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
733                       Z_U_tc(IJK) = Z_NODE(8)
734     
735     !======================================================================
736     !  Get Connectivity
737     !======================================================================
738     
739                 IF(NO_K) THEN
740                    MIN_INTERSECTIONS = 2
741                    MAX_INTERSECTIONS = 2
742                    N_N1 = 5
743                    N_N2 = 8
744     
745                 ELSE
746                    MIN_INTERSECTIONS = 3
747                    MAX_INTERSECTIONS = 6
748                    N_N1 = 1
749                    N_N2 = 8
750                 ENDIF
751     
752                 CALL GET_CONNECTIVITY(IJK,'U_MOMENTUM',NUMBER_OF_NEW_U_POINTS,NUMBER_OF_U_NODES(IJK),CONNECTIVITY_U,&
753                 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)
754     
755     
756                 IF(TOTAL_NUMBER_OF_INTERSECTIONS < MIN_INTERSECTIONS ) THEN   ! Not a cut cell
757     
758                    Q_ID = 1
759                    CALL EVAL_F('QUADRIC',X_NODE(0),Y_NODE(0),Z_NODE(0),Q_ID,F_NODE(0),CLIP_FLAG)
760     
761                    CALL EVAL_F('POLYGON',X_NODE(0),Y_NODE(0),Z_NODE(0),N_POLYGON,F_NODE(0),CLIP_FLAG)
762     
763                    CALL EVAL_F('USR_DEF',X_NODE(0),Y_NODE(0),Z_NODE(0),N_USR_DEF,F_NODE(0),CLIP_FLAG)
764     
765                    CALL EVAL_STL_FCT_AT('U_MOMENTUM',IJK,0,F_NODE(0),CLIP_FLAG,BCID)
766     
767                    IF(TOTAL_NUMBER_OF_INTERSECTIONS==-1) F_NODE(0) = -ONE
768                    BC_U_ID(IJK) = 0
769     
770                    IF(F_NODE(0) <= ZERO) THEN
771                       BLOCKED_U_CELL_AT(IJK) = .FALSE.
772                       STANDARD_U_CELL_AT(IJK) = .TRUE.          ! Regular fluid cell
773     
774                       X_U_ec(IJK) = X_NODE(8)
775                       Y_U_ec(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
776                       Z_U_ec(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
777     
778                       X_U_nc(IJK) = HALF * (X_NODE(7) + X_NODE(8))
779                       Y_U_nc(IJK) = Y_NODE(8)
780                       Z_U_nc(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
781     
782                       X_U_tc(IJK) = HALF * (X_NODE(7) + X_NODE(8))
783                       Y_U_tc(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
784                       Z_U_tc(IJK) = Z_NODE(8)
785     
786                    ELSE
787                       BLOCKED_U_CELL_AT(IJK) = .TRUE.           ! Blocked fluid cell
788                       STANDARD_U_CELL_AT(IJK) = .FALSE.
789                       AXY_U(IJK) = ZERO
790                       AXZ_U(IJK) = ZERO
791                       AYZ_U(IJK) = ZERO
792                       VOL_U(IJK) = ZERO
793                    ENDIF
794     
795                    IF(NO_K) THEN
796                       NUMBER_OF_U_NODES(IJK) = 4
797                       CONNECTIVITY_U(IJK,1) = IJK_OF_NODE(5)
798                       CONNECTIVITY_U(IJK,2) = IJK_OF_NODE(6)
799                       CONNECTIVITY_U(IJK,3) = IJK_OF_NODE(8)
800                       CONNECTIVITY_U(IJK,4) = IJK_OF_NODE(7)
801                    ELSE
802                       NUMBER_OF_U_NODES(IJK) = 8
803                       DO NODE = N_N1,N_N2
804                          CONNECTIVITY_U(IJK,NODE) = IJK_OF_NODE(NODE)
805                       END DO
806                    ENDIF
807     
808                 ELSE IF(TOTAL_NUMBER_OF_INTERSECTIONS > MAX_INTERSECTIONS) THEN
809     
810                    IF(PRINT_WARNINGS) THEN
811                       WRITE(*,*) 'TOO MANY INTERSECTIONS FOUND IN U-CELL IJK = ',IJK
812                       WRITE(*,*) 'MAXIMUM NUMBER OF INTERSECTIONS = ',MAX_INTERSECTIONS
813                       WRITE(*,*) 'TOTAL NUMBER OF INTERSECTIONS = ',TOTAL_NUMBER_OF_INTERSECTIONS
814                       WRITE(*,*) 'REMOVING U-CELL FROM COMPUTATION.'
815     !                  WRITE(*,*) 'MFIX WILL EXIT NOW.'
816     !                  CALL MFIX_EXIT(MYPE)
817     
818                    ENDIF
819     
820                    BLOCKED_U_CELL_AT(IJK) = .TRUE.           ! Blocked fluid cell
821                    STANDARD_U_CELL_AT(IJK) = .FALSE.
822                    CUT_U_CELL_AT(IJK) = .FALSE.
823                    AXY_U(IJK) = ZERO
824                    AXZ_U(IJK) = ZERO
825                    AYZ_U(IJK) = ZERO
826                    VOL_U(IJK) = ZERO
827     
828                 ELSE                                         ! Cut cell
829     
830                    CUT_U_CELL_AT(IJK) = .TRUE.
831                    BLOCKED_U_CELL_AT(IJK) = .FALSE.
832                    STANDARD_U_CELL_AT(IJK) = .FALSE.
833     
834                    DO NODE = N_N1,N_N2
835                       IF(F_NODE(NODE) < - TOL_F) THEN
836                          NUMBER_OF_U_NODES(IJK) = NUMBER_OF_U_NODES(IJK) + 1
837                          CONNECTIVITY_U(IJK,NUMBER_OF_U_NODES(IJK)) = IJK_OF_NODE(NODE)
838                       ENDIF
839                    END DO
840     
841                    CALL GET_CUT_CELL_VOLUME_AND_AREAS(IJK,'U_MOMENTUM',NUMBER_OF_U_NODES(IJK),CONNECTIVITY_U,&
842                    X_NEW_U_POINT,Y_NEW_U_POINT,Z_NEW_U_POINT)
843     
844                 ENDIF
845     
846              ENDIF
847     
848           END DO
849     
850           DO IJK = IJKSTART3, IJKEND3
851     
852              CALL GET_CELL_NODE_COORDINATES(IJK,'U_MOMENTUM')
853     
854              DELX_Ue(IJK) = X_NODE(8) - X_U(IJK)
855              DELX_Uw(IJK) = X_U(IJK) - X_NODE(1)
856     
857              DELY_Un(IJK) = Y_NODE(8) - Y_U(IJK)
858              DELY_Us(IJK) = Y_U(IJK) - Y_NODE(1)
859     
860              DELZ_Ut(IJK) = Z_NODE(8) - Z_U(IJK)
861              DELZ_Ub(IJK) = Z_U(IJK) - Z_NODE(1)
862     
863           ENDDO
864     
865           RETURN
866     
867           END SUBROUTINE SET_3D_CUT_U_CELL_FLAGS
868     
869     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
870     !                                                                      C
871     !  Module name: SET_3D_CUT_V_CELL_FLAGS                                C
872     !  Purpose: Set flags for V-Momentum cut cells, based on intersection  C
873     !  of the grid with the quadric(s)                                     C
874     !                                                                      C
875     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
876     !  Reviewer:                                          Date:            C
877     !                                                                      C
878     !  Revision Number #                                  Date: ##-###-##  C
879     !  Author: #                                                           C
880     !  Purpose: #                                                          C
881     !                                                                      C
882     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
883       SUBROUTINE SET_3D_CUT_V_CELL_FLAGS
884     
885           USE compar, ONLY: mype, pe_io, ijkstart3, ijkend3
886           USE cut_cell_preproc, ONLY: cad_intersect, clean_intersect, clean_intersect_scalar, eval_f, intersect
887           USE cutcell
888           USE geometry, ONLY: no_k, axy_v, axz_v, ayz_v, vol_v
889           USE param1, ONLY: half, one, zero
890           USE polygon, ONLY: n_polygon
891           USE quadric, ONLY: tol_f
892     
893           IMPLICIT NONE
894           INTEGER :: IJK
895           INTEGER :: TOTAL_NUMBER_OF_INTERSECTIONS
896           INTEGER :: NODE,N_N1,N_N2,Q_ID
897           INTEGER :: MIN_INTERSECTIONS,MAX_INTERSECTIONS
898           LOGICAL :: CLIP_FLAG
899           INTEGER :: BCID
900     
901           IF(MyPE == PE_IO) THEN
902              WRITE(*,*)'INTERSECTING GEOMETRY WITH V-MOMENTUM CELLS...'
903           ENDIF
904     !======================================================================
905     !  Intersection between quadric Grid
906     !======================================================================
907     
908           INTERSECT_X = .FALSE.
909           INTERSECT_Y = .FALSE.
910           INTERSECT_Z = .FALSE.
911           SNAP = .FALSE.
912     
913           IF(.NOT.(USE_STL.OR.USE_MSH)) THEN
914              DO IJK = IJKSTART3, IJKEND3
915     !             IF(POTENTIAL_CUT_CELL_AT(IJK))  CALL INTERSECT(IJK,'V_MOMENTUM',Xn_V_int(IJK),Ye_V_int(IJK),Zt_V_int(IJK))
916                 CALL INTERSECT(IJK,'V_MOMENTUM',Xn_V_int(IJK),Ye_V_int(IJK),Zt_V_int(IJK))
917              END DO
918     
919     !======================================================================
920     !  Clean-up intersection flags in preparaton of small cells removal
921     !======================================================================
922              DO IJK = IJKSTART3, IJKEND3
923                 IF(INTERIOR_CELL_AT(IJK)) THEN
924     !                IF(POTENTIAL_CUT_CELL_AT(IJK)) CALL CLEAN_INTERSECT(IJK,'V_MOMENTUM',Xn_V_int(IJK),Ye_V_int(IJK),Zt_V_int(IJK))
925                    CALL CLEAN_INTERSECT(IJK,'V_MOMENTUM',Xn_V_int(IJK),Ye_V_int(IJK),Zt_V_int(IJK))
926                 ENDIF
927              END DO
928     
929           ELSE
930              CALL CAD_INTERSECT('V_MOMENTUM',Xn_V_int,Ye_V_int,Zt_V_int)
931           ENDIF
932     
933     
934           NUMBER_OF_NEW_V_POINTS = 0
935     
936           DO IJK = IJKSTART3, IJKEND3
937              CALL WRITE_PROGRESS_BAR(IJK,IJKEND3 - IJKSTART3 + 1,'C')
938     
939              IF(INTERIOR_CELL_AT(IJK)) THEN
940     
941     !======================================================================
942     !  Get coordinates of eight nodes
943     !======================================================================
944     
945                 CALL GET_CELL_NODE_COORDINATES(IJK,'V_MOMENTUM')
946     
947     !======================================================================
948     !  Initialize location of velocity nodes at center of E,W,T faces
949     !======================================================================
950     
951                 X_V_ec(IJK) = X_NODE(8)
952                 Y_V_ec(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
953                 Z_V_ec(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
954     
955                 X_V_nc(IJK) = HALF * (X_NODE(7) + X_NODE(8))
956                 Y_V_nc(IJK) = Y_NODE(8)
957                 Z_V_nc(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
958     
959                 X_V_tc(IJK) = HALF * (X_NODE(7) + X_NODE(8))
960                 Y_V_tc(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
961                 Z_V_tc(IJK) = Z_NODE(8)
962     
963     !======================================================================
964     !  Get Connectivity
965     !======================================================================
966     
967                 IF(NO_K) THEN
968                    MIN_INTERSECTIONS = 2
969                    MAX_INTERSECTIONS = 2
970                    N_N1 = 5
971                    N_N2 = 8
972     
973                 ELSE
974                    MIN_INTERSECTIONS = 3
975                    MAX_INTERSECTIONS = 6
976                    N_N1 = 1
977                    N_N2 = 8
978                 ENDIF
979     
980     
981                 CALL GET_CONNECTIVITY(IJK,'V_MOMENTUM',NUMBER_OF_NEW_V_POINTS,NUMBER_OF_V_NODES(IJK),CONNECTIVITY_V,&
982                 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)
983     
984     
985                 IF(TOTAL_NUMBER_OF_INTERSECTIONS < MIN_INTERSECTIONS ) THEN   ! Not a cut cell
986     
987                    Q_ID = 1
988                    CALL EVAL_F('QUADRIC',X_NODE(0),Y_NODE(0),Z_NODE(0),Q_ID,F_NODE(0),CLIP_FLAG)
989     
990                    CALL EVAL_F('POLYGON',X_NODE(0),Y_NODE(0),Z_NODE(0),N_POLYGON,F_NODE(0),CLIP_FLAG)
991     
992                    CALL EVAL_F('USR_DEF',X_NODE(0),Y_NODE(0),Z_NODE(0),N_USR_DEF,F_NODE(0),CLIP_FLAG)
993     
994                    CALL EVAL_STL_FCT_AT('V_MOMENTUM',IJK,0,F_NODE(0),CLIP_FLAG,BCID)
995     
996                    IF(TOTAL_NUMBER_OF_INTERSECTIONS==-1) F_NODE(0) = -ONE
997                    BC_V_ID(IJK) = 0
998     
999                    IF(F_NODE(0) <= ZERO) THEN
1000                       BLOCKED_V_CELL_AT(IJK) = .FALSE.
1001                       STANDARD_V_CELL_AT(IJK) = .TRUE.          ! Regular fluid cell
1002     
1003                       X_V_ec(IJK) = X_NODE(8)
1004                       Y_V_ec(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
1005                       Z_V_ec(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
1006     
1007                       X_V_nc(IJK) = HALF * (X_NODE(7) + X_NODE(8))
1008                       Y_V_nc(IJK) = Y_NODE(8)
1009                       Z_V_nc(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
1010     
1011                       X_V_tc(IJK) = HALF * (X_NODE(7) + X_NODE(8))
1012                       Y_V_tc(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
1013                       Z_V_tc(IJK) = Z_NODE(8)
1014     
1015                    ELSE
1016                       BLOCKED_V_CELL_AT(IJK) = .TRUE.           ! Blocked fluid cell
1017                       STANDARD_V_CELL_AT(IJK) = .FALSE.
1018                       AXY_V(IJK) = ZERO
1019                       AXZ_V(IJK) = ZERO
1020                       AYZ_V(IJK) = ZERO
1021                       VOL_V(IJK) = ZERO
1022                    ENDIF
1023     
1024                    IF(NO_K) THEN
1025                       NUMBER_OF_V_NODES(IJK) = 4
1026                       CONNECTIVITY_V(IJK,1) = IJK_OF_NODE(5)
1027                       CONNECTIVITY_V(IJK,2) = IJK_OF_NODE(6)
1028                       CONNECTIVITY_V(IJK,3) = IJK_OF_NODE(8)
1029                       CONNECTIVITY_V(IJK,4) = IJK_OF_NODE(7)
1030                    ELSE
1031                       NUMBER_OF_V_NODES(IJK) = 8
1032                       DO NODE = N_N1,N_N2
1033                          CONNECTIVITY_V(IJK,NODE) = IJK_OF_NODE(NODE)
1034                       END DO
1035                    ENDIF
1036     
1037                 ELSE IF(TOTAL_NUMBER_OF_INTERSECTIONS > MAX_INTERSECTIONS ) THEN
1038     
1039                    IF(PRINT_WARNINGS) THEN
1040                       WRITE(*,*) 'TOO MANY INTERSECTIONS FOUND IN V-CELL IJK = ',IJK
1041                       WRITE(*,*) 'MAXIMUM NUMBER OF INTERSECTIONS = ',MAX_INTERSECTIONS
1042                       WRITE(*,*) 'TOTAL NUMBER OF INTERSECTIONS = ',TOTAL_NUMBER_OF_INTERSECTIONS
1043                       WRITE(*,*) 'REMOVING V-CELL FROM COMPUTATION.'
1044     !                  WRITE(*,*) 'MFIX WILL EXIT NOW.'
1045     !                  CALL MFIX_EXIT(MYPE)
1046     
1047                    ENDIF
1048     
1049                    BLOCKED_V_CELL_AT(IJK) = .TRUE.           ! Blocked fluid cell
1050                    STANDARD_V_CELL_AT(IJK) = .FALSE.
1051                    CUT_V_CELL_AT(IJK) = .FALSE.
1052                    AXY_V(IJK) = ZERO
1053                    AXZ_V(IJK) = ZERO
1054                    AYZ_V(IJK) = ZERO
1055                    VOL_V(IJK) = ZERO
1056     
1057                 ELSE                                         ! Cut cell
1058     
1059                    CUT_V_CELL_AT(IJK) = .TRUE.
1060                    BLOCKED_V_CELL_AT(IJK) = .FALSE.
1061                    STANDARD_V_CELL_AT(IJK) = .FALSE.
1062     
1063                    DO NODE = N_N1,N_N2
1064                       IF(F_NODE(NODE) < - TOL_F) THEN
1065                          NUMBER_OF_V_NODES(IJK) = NUMBER_OF_V_NODES(IJK) + 1
1066                          CONNECTIVITY_V(IJK,NUMBER_OF_V_NODES(IJK)) = IJK_OF_NODE(NODE)
1067                       ENDIF
1068                    END DO
1069     
1070                    CALL GET_CUT_CELL_VOLUME_AND_AREAS(IJK,'V_MOMENTUM',NUMBER_OF_V_NODES(IJK),CONNECTIVITY_V,&
1071                    X_NEW_V_POINT,Y_NEW_V_POINT,Z_NEW_V_POINT)
1072     
1073                 ENDIF
1074     
1075              ENDIF
1076           END DO
1077     
1078           DO IJK = IJKSTART3, IJKEND3
1079     
1080              CALL GET_CELL_NODE_COORDINATES(IJK,'V_MOMENTUM')
1081     
1082              DELX_Ve(IJK) = X_NODE(8) - X_V(IJK)
1083              DELX_Vw(IJK) = X_V(IJK) - X_NODE(1)
1084     
1085              DELY_Vn(IJK) = Y_NODE(8) - Y_V(IJK)
1086              DELY_Vs(IJK) = Y_V(IJK) - Y_NODE(1)
1087     
1088              DELZ_Vt(IJK) = Z_NODE(8) - Z_V(IJK)
1089              DELZ_Vb(IJK) = Z_V(IJK) - Z_NODE(1)
1090     
1091           ENDDO
1092     
1093           RETURN
1094     
1095           END SUBROUTINE SET_3D_CUT_V_CELL_FLAGS
1096     
1097     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1098     !                                                                      C
1099     !  Module name: SET_3D_CUT_W_CELL_FLAGS                                C
1100     !  Purpose: Set flags for W-Momentum cut cells, based on intersection  C
1101     !  of the grid with the quadric(s)                                     C
1102     !                                                                      C
1103     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1104     !  Reviewer:                                          Date:            C
1105     !                                                                      C
1106     !  Revision Number #                                  Date: ##-###-##  C
1107     !  Author: #                                                           C
1108     !  Purpose: #                                                          C
1109     !                                                                      C
1110     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1111       SUBROUTINE SET_3D_CUT_W_CELL_FLAGS
1112     
1113           USE compar, ONLY: mype, pe_io, ijkstart3, ijkend3
1114           USE cut_cell_preproc, ONLY: cad_intersect, clean_intersect, clean_intersect_scalar, eval_f, intersect
1115           USE cutcell
1116           USE geometry, ONLY: axy_w, axz_w, ayz_w, vol_w
1117           USE param1, ONLY: half, one, zero
1118           USE polygon, ONLY: n_polygon
1119           USE quadric, ONLY: tol_f
1120     
1121           IMPLICIT NONE
1122           INTEGER :: IJK
1123           INTEGER :: TOTAL_NUMBER_OF_INTERSECTIONS
1124           INTEGER :: NODE,Q_ID
1125           LOGICAL :: CLIP_FLAG
1126           INTEGER :: BCID
1127     
1128           IF(MyPE == PE_IO) THEN
1129              WRITE(*,10)'INTERSECTING GEOMETRY WITH W-MOMENTUM CELLS...'
1130           ENDIF
1131     10    FORMAT(1X,A)
1132     !======================================================================
1133     !  Intersection between quadric Grid
1134     !======================================================================
1135     
1136           INTERSECT_X = .FALSE.
1137           INTERSECT_Y = .FALSE.
1138           INTERSECT_Z = .FALSE.
1139           SNAP = .FALSE.
1140     
1141           IF(.NOT.(USE_STL.OR.USE_MSH)) THEN
1142              DO IJK = IJKSTART3, IJKEND3
1143     !             IF(POTENTIAL_CUT_CELL_AT(IJK)) CALL INTERSECT(IJK,'W_MOMENTUM',Xn_W_int(IJK),Ye_W_int(IJK),Zt_W_int(IJK))
1144                 CALL INTERSECT(IJK,'W_MOMENTUM',Xn_W_int(IJK),Ye_W_int(IJK),Zt_W_int(IJK))
1145              END DO
1146     
1147     !======================================================================
1148     !  Clean-up intersection flags in preparaton of small cells removal
1149     !======================================================================
1150              DO IJK = IJKSTART3, IJKEND3
1151                 IF(INTERIOR_CELL_AT(IJK)) THEN
1152     !               IF(POTENTIAL_CUT_CELL_AT(IJK)) CALL CLEAN_INTERSECT(IJK,'W_MOMENTUM',Xn_W_int(IJK),Ye_W_int(IJK),Zt_W_int(IJK))
1153                    CALL CLEAN_INTERSECT(IJK,'W_MOMENTUM',Xn_W_int(IJK),Ye_W_int(IJK),Zt_W_int(IJK))
1154                 ENDIF
1155              END DO
1156     
1157           ELSE
1158              CALL CAD_INTERSECT('W_MOMENTUM',Xn_W_int,Ye_W_int,Zt_W_int)
1159           ENDIF
1160     
1161     
1162           NUMBER_OF_NEW_W_POINTS = 0
1163     
1164           DO IJK = IJKSTART3, IJKEND3
1165     
1166              CALL WRITE_PROGRESS_BAR(IJK,IJKEND3 - IJKSTART3 + 1,'C')
1167     
1168              IF(INTERIOR_CELL_AT(IJK)) THEN
1169     
1170     !======================================================================
1171     !  Get coordinates of eight nodes
1172     !======================================================================
1173     
1174                 CALL GET_CELL_NODE_COORDINATES(IJK,'W_MOMENTUM')
1175     
1176     !======================================================================
1177     !  Initialize location of velocity nodes at center of E,W,T faces
1178     !======================================================================
1179     
1180                 X_W_ec(IJK) = X_NODE(8)
1181                 Y_W_ec(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
1182                 Z_W_ec(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
1183     
1184                 X_W_nc(IJK) = HALF * (X_NODE(7) + X_NODE(8))
1185                 Y_W_nc(IJK) = Y_NODE(8)
1186                 Z_W_nc(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
1187     
1188                 X_W_tc(IJK) = HALF * (X_NODE(7) + X_NODE(8))
1189                 Y_W_tc(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
1190                 Z_W_tc(IJK) = Z_NODE(8)
1191     
1192     !======================================================================
1193     !  Get Connectivity
1194     !======================================================================
1195     
1196                 CALL GET_CONNECTIVITY(IJK,'W_MOMENTUM',NUMBER_OF_NEW_W_POINTS,NUMBER_OF_W_NODES(IJK),CONNECTIVITY_W,&
1197                 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)
1198     
1199     
1200                 IF(TOTAL_NUMBER_OF_INTERSECTIONS < 3 ) THEN   ! Not a cut cell
1201     
1202                    Q_ID = 1
1203                    CALL EVAL_F('QUADRIC',X_NODE(0),Y_NODE(0),Z_NODE(0),Q_ID,F_NODE(0),CLIP_FLAG)
1204     
1205                    CALL EVAL_F('POLYGON',X_NODE(0),Y_NODE(0),Z_NODE(0),N_POLYGON,F_NODE(0),CLIP_FLAG)
1206     
1207                    CALL EVAL_F('USR_DEF',X_NODE(0),Y_NODE(0),Z_NODE(0),N_USR_DEF,F_NODE(0),CLIP_FLAG)
1208     
1209                    CALL EVAL_STL_FCT_AT('W_MOMENTUM',IJK,0,F_NODE(0),CLIP_FLAG,BCID)
1210     
1211                    IF(TOTAL_NUMBER_OF_INTERSECTIONS==-1) F_NODE(0) = -ONE
1212                    BC_W_ID(IJK) = 0
1213     
1214                    IF(F_NODE(0) <= ZERO) THEN
1215                       BLOCKED_W_CELL_AT(IJK) = .FALSE.
1216                       STANDARD_W_CELL_AT(IJK) = .TRUE.          ! Regular fluid cell
1217     
1218                       X_W_ec(IJK) = X_NODE(8)
1219                       Y_W_ec(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
1220                       Z_W_ec(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
1221     
1222                       X_W_nc(IJK) = HALF * (X_NODE(7) + X_NODE(8))
1223                       Y_W_nc(IJK) = Y_NODE(8)
1224                       Z_W_nc(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
1225     
1226                       X_W_tc(IJK) = HALF * (X_NODE(7) + X_NODE(8))
1227                       Y_W_tc(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
1228                       Z_W_tc(IJK) = Z_NODE(8)
1229     
1230                    ELSE
1231                       BLOCKED_W_CELL_AT(IJK) = .TRUE.           ! Blocked fluid cell
1232                       STANDARD_W_CELL_AT(IJK) = .FALSE.
1233                       AXY_W(IJK) = ZERO
1234                       AXZ_W(IJK) = ZERO
1235                       AYZ_W(IJK) = ZERO
1236                       VOL_W(IJK) = ZERO
1237                    ENDIF
1238     
1239                    NUMBER_OF_W_NODES(IJK) = 8
1240                    DO NODE = 1,8
1241                       CONNECTIVITY_W(IJK,NODE) = IJK_OF_NODE(NODE)
1242                    END DO
1243     
1244                 ELSE IF(TOTAL_NUMBER_OF_INTERSECTIONS > 6 ) THEN
1245     
1246                    IF(PRINT_WARNINGS) THEN
1247                       WRITE(*,*) 'TOO MANY INTERSECTIONS FOUND IN W-CELL IJK = ',IJK
1248                       WRITE(*,*) 'MAXIMUM NUMBER OF INTERSECTIONS = ',6
1249                       WRITE(*,*) 'TOTAL NUMBER OF INTERSECTIONS = ',TOTAL_NUMBER_OF_INTERSECTIONS
1250                       WRITE(*,*) 'REMOVING W-CELL FROM COMPUTATION.'
1251     !                  WRITE(*,*) 'MFIX WILL EXIT NOW.'
1252     !                  CALL MFIX_EXIT(MYPE)
1253     
1254                    ENDIF
1255     
1256                    BLOCKED_W_CELL_AT(IJK) = .TRUE.           ! Blocked fluid cell
1257                    STANDARD_W_CELL_AT(IJK) = .FALSE.
1258                    CUT_W_CELL_AT(IJK) = .FALSE.
1259                    AXY_W(IJK) = ZERO
1260                    AXZ_W(IJK) = ZERO
1261                    AYZ_W(IJK) = ZERO
1262                    VOL_W(IJK) = ZERO
1263     
1264                 ELSE                                         ! Cut cell
1265     
1266                    CUT_W_CELL_AT(IJK) = .TRUE.
1267     
1268                    DO NODE = 1,8
1269                       IF(F_NODE(NODE) < - TOL_F) THEN
1270                          NUMBER_OF_W_NODES(IJK) = NUMBER_OF_W_NODES(IJK) + 1
1271                          CONNECTIVITY_W(IJK,NUMBER_OF_W_NODES(IJK)) = IJK_OF_NODE(NODE)
1272                       ENDIF
1273                    END DO
1274     
1275     
1276                    CALL GET_CUT_CELL_VOLUME_AND_AREAS(IJK,'W_MOMENTUM',NUMBER_OF_W_NODES(IJK),CONNECTIVITY_W,&
1277                    X_NEW_W_POINT,Y_NEW_W_POINT,Z_NEW_W_POINT)
1278     
1279                 ENDIF
1280     
1281              ENDIF
1282           END DO
1283     
1284           DO IJK = IJKSTART3, IJKEND3
1285     
1286              CALL GET_CELL_NODE_COORDINATES(IJK,'W_MOMENTUM')
1287     
1288              DELX_We(IJK) = X_NODE(8) - X_W(IJK)
1289              DELX_Ww(IJK) = X_W(IJK) - X_NODE(1)
1290     
1291              DELY_Wn(IJK) = Y_NODE(8) - Y_W(IJK)
1292              DELY_Ws(IJK) = Y_W(IJK) - Y_NODE(1)
1293     
1294              DELZ_Wt(IJK) = Z_NODE(8) - Z_W(IJK)
1295              DELZ_Wb(IJK) = Z_W(IJK) - Z_NODE(1)
1296     
1297           ENDDO
1298     
1299           RETURN
1300     
1301           END SUBROUTINE SET_3D_CUT_W_CELL_FLAGS
1302     
1303     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1304     !                                                                      C
1305     !  Module name: SET_3D_CUT_CELL_TREATMENT_FLAGS                        C
1306     !  Purpose: Set flags for scalar cut cells, based on intersection      C
1307     !  of the grid with the quadric(s)                                     C
1308     !                                                                      C
1309     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1310     !  Reviewer:                                          Date:            C
1311     !                                                                      C
1312     !  Revision Number #                                  Date: ##-###-##  C
1313     !  Author: #                                                           C
1314     !  Purpose: #                                                          C
1315     !                                                                      C
1316     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1317       SUBROUTINE SET_3D_CUT_CELL_TREATMENT_FLAGS
1318     
1319           USE compar, only: mype, pe_io, ijkstart3, ijkend3
1320           USE cutcell
1321           USE functions, ONLY: funijk
1322           USE geometry, ONLY: do_k
1323           USE indices, ONLY: i_of, j_of, k_of
1324     
1325           IMPLICIT NONE
1326           INTEGER :: IJK,I,J,K,IM,IP,JM,JP,KM,KP
1327           INTEGER :: IMJK,IPJK,IJMK,IJPK,IJKM,IJKP
1328     
1329           IF(MyPE == PE_IO) THEN
1330              WRITE(*,*)'SETTING CUT CELL TREATMENT FLAGS...'
1331           ENDIF
1332     !======================================================================
1333     !  Set flags identifying cells requiring cut cell treatment:
1334     !  These are the cut cells and their neighbours
1335     !======================================================================
1336     
1337           CUT_TREATMENT_AT   = .FALSE.
1338           CUT_U_TREATMENT_AT = .FALSE.
1339           CUT_V_TREATMENT_AT = .FALSE.
1340           CUT_W_TREATMENT_AT = .FALSE.
1341     
1342           call SEND_RECEIVE_1D_LOGICAL(CUT_CELL_AT,2)
1343           call SEND_RECEIVE_1D_LOGICAL(CUT_U_CELL_AT,2)
1344           call SEND_RECEIVE_1D_LOGICAL(CUT_V_CELL_AT,2)
1345           call SEND_RECEIVE_1D_LOGICAL(CUT_W_CELL_AT,2)
1346     
1347           DO IJK = IJKSTART3, IJKEND3
1348     
1349              CALL WRITE_PROGRESS_BAR(IJK,IJKEND3 - IJKSTART3 + 1,'C')
1350     
1351              IF(INTERIOR_CELL_AT(IJK)) THEN
1352     
1353                 I = I_OF(IJK)
1354                 J = J_OF(IJK)
1355                 K = K_OF(IJK)
1356     
1357                 IP = I + 1
1358                 JP = J + 1
1359                 KP = K + 1
1360     
1361                 IM = I - 1
1362                 JM = J - 1
1363                 KM = K - 1
1364     
1365                 IMJK   = FUNIJK(IM,J,K)
1366                 IJMK   = FUNIJK(I,JM,K)
1367                 IJKM   = FUNIJK(I,J,KM)
1368     
1369                 IPJK   = FUNIJK(IP,J,K)
1370                 IJPK   = FUNIJK(I,JP,K)
1371                 IJKP   = FUNIJK(I,J,KP)
1372     
1373                 CUT_TREATMENT_AT(IJK) = (CUT_CELL_AT(IJK ).OR.   &
1374                                          CUT_CELL_AT(IMJK).OR.   &
1375                                          CUT_CELL_AT(IPJK).OR.   &
1376                                          CUT_CELL_AT(IJMK).OR.   &
1377                                          CUT_CELL_AT(IJPK))
1378     
1379                 CUT_U_TREATMENT_AT(IJK) = (CUT_U_CELL_AT(IJK ).OR.   &
1380                                            CUT_U_CELL_AT(IMJK).OR.   &
1381                                            CUT_U_CELL_AT(IPJK).OR.   &
1382                                            CUT_U_CELL_AT(IJMK).OR.   &
1383                                            CUT_U_CELL_AT(IJPK))
1384     
1385                 CUT_V_TREATMENT_AT(IJK) = (CUT_V_CELL_AT(IJK ).OR.   &
1386                                            CUT_V_CELL_AT(IMJK).OR.   &
1387                                            CUT_V_CELL_AT(IPJK).OR.   &
1388                                            CUT_V_CELL_AT(IJMK).OR.   &
1389                                            CUT_V_CELL_AT(IJPK))
1390     
1391                 IF(DO_K) THEN
1392     
1393                    CUT_TREATMENT_AT(IJK) = (CUT_TREATMENT_AT(IJK).OR.   &
1394                                                 CUT_CELL_AT(IJKM).OR.   &
1395                                                 CUT_CELL_AT(IJKP))
1396     
1397                    CUT_U_TREATMENT_AT(IJK) = (CUT_U_TREATMENT_AT(IJK).OR.   &
1398                                                   CUT_U_CELL_AT(IJKM).OR.   &
1399                                                   CUT_U_CELL_AT(IJKP))
1400     
1401                    CUT_V_TREATMENT_AT(IJK) = (CUT_V_TREATMENT_AT(IJK).OR.   &
1402                                                   CUT_V_CELL_AT(IJKM).OR.   &
1403                                                   CUT_V_CELL_AT(IJKP))
1404     
1405                    CUT_W_TREATMENT_AT(IJK) = (CUT_W_CELL_AT(IJK ).OR.   &
1406                                               CUT_W_CELL_AT(IMJK).OR.   &
1407                                               CUT_W_CELL_AT(IPJK).OR.   &
1408                                               CUT_W_CELL_AT(IJMK).OR.   &
1409                                               CUT_W_CELL_AT(IJPK).OR.   &
1410                                               CUT_W_CELL_AT(IJKM).OR.   &
1411                                               CUT_W_CELL_AT(IJKP))
1412     
1413                 ENDIF
1414     
1415              ENDIF
1416     
1417              IF(.NOT.RE_INDEXING) THEN
1418                 IF(BLOCKED_CELL_AT(IJK)) CUT_TREATMENT_AT(IJK) = .FALSE.
1419                 IF(BLOCKED_U_CELL_AT(IJK)) CUT_U_TREATMENT_AT(IJK) = .FALSE.
1420                 IF(BLOCKED_V_CELL_AT(IJK)) CUT_V_TREATMENT_AT(IJK) = .FALSE.
1421                 IF(BLOCKED_W_CELL_AT(IJK)) CUT_W_TREATMENT_AT(IJK) = .FALSE.
1422              ENDIF
1423     
1424           END DO
1425     
1426           IF(CG_SAFE_MODE(1)==1) CUT_TREATMENT_AT   = .FALSE.
1427           IF(CG_SAFE_MODE(3)==1) CUT_U_TREATMENT_AT = .FALSE.
1428           IF(CG_SAFE_MODE(4)==1) CUT_V_TREATMENT_AT = .FALSE.
1429           IF(CG_SAFE_MODE(5)==1) CUT_W_TREATMENT_AT = .FALSE.
1430     
1431           RETURN
1432     
1433           END SUBROUTINE SET_3D_CUT_CELL_TREATMENT_FLAGS
1434     
1435     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1436     !                                                                      C
1437     !  Module name: SET_GHOST_CELL_FLAGS                                   C
1438     !  Purpose: Set flags for ghost cell flags                             C
1439     !                                                                      C
1440     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1441     !  Reviewer:                                          Date:            C
1442     !                                                                      C
1443     !  Revision Number #                                  Date: ##-###-##  C
1444     !  Author: #                                                           C
1445     !  Purpose: #                                                          C
1446     !                                                                      C
1447     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1448       SUBROUTINE SET_GHOST_CELL_FLAGS
1449     
1450           USE compar, ONLY: ijkstart3, ijkend3
1451           USE compar, ONLY: iend1, iend2, iend3, jstart1, jstart2, jstart3, kstart1, kstart2, kstart3
1452           USE compar, only: istart1, istart2, istart3, jend1, jend2, jend3, kend1, kend2, kend3
1453           USE cutcell
1454           USE functions, ONLY: funijk
1455           USE functions, ONLY: WEST_OF, EAST_OF, SOUTH_OF, NORTH_OF, BOTTOM_OF, TOP_OF
1456           USE functions, ONLY: IM_OF, IP_OF, JM_OF, JP_OF, KM_OF, KP_OF
1457           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
1458     
1459           USE bc
1460           USE sendrecv
1461     
1462           IMPLICIT NONE
1463           INTEGER :: IJK,I,J,K,I23,J23,K23
1464           INTEGER :: BCID
1465           INTEGER :: IPJK,IMJK,IJPK,IJMK,IJKP,IJKM
1466     
1467     !     EAST BOUNDARY
1468           I = IEND1
1469     
1470           IF(I==IMAX1) THEN
1471     
1472              DO I23 = IEND2,IEND3
1473                 DO K = KSTART3,KEND3
1474                    DO J = JSTART3, JEND3
1475     
1476                       IJK  = FUNIJK(I,J,K)
1477                       IPJK = FUNIJK(I23,J,K)
1478     
1479                       BLOCKED_CELL_AT(IPJK) = BLOCKED_CELL_AT(IJK)
1480                       BLOCKED_U_CELL_AT(IPJK) = BLOCKED_U_CELL_AT(IJK)
1481                       BLOCKED_V_CELL_AT(IPJK) = BLOCKED_V_CELL_AT(IJK)
1482                       BLOCKED_W_CELL_AT(IPJK) = BLOCKED_W_CELL_AT(IJK)
1483     
1484                       IF(BLOCKED_CELL_AT(IPJK)) FLAG(IPJK) = 100
1485     
1486                       VOL(IPJK)   = VOL(IJK)
1487                       VOL_U(IPJK) = VOL_U(IJK)
1488                       VOL_V(IPJK) = VOL_V(IJK)
1489                       VOL_W(IPJK) = VOL_W(IJK)
1490     
1491                       AYZ(IPJK)   = AYZ(IJK)
1492                       AYZ_U(IPJK) = AYZ_U(IJK)
1493                       AYZ_V(IPJK) = AYZ_V(IJK)
1494                       AYZ_W(IPJK) = AYZ_W(IJK)
1495     
1496                       AXY(IPJK)   = AXY(IJK)
1497                       AXY_U(IPJK) = AXY_U(IJK)
1498                       AXY_V(IPJK) = AXY_V(IJK)
1499                       AXY_W(IPJK) = AXY_W(IJK)
1500     
1501                       AXZ(IPJK)   = AXZ(IJK)
1502                       AXZ_U(IPJK) = AXZ_U(IJK)
1503                       AXZ_V(IPJK) = AXZ_V(IJK)
1504                       AXZ_W(IPJK) = AXZ_W(IJK)
1505     
1506                       BC_ID(IPJK)   = BC_ID(IJK)
1507                       BC_U_ID(IPJK) = BC_U_ID(IJK)
1508                       BC_V_ID(IPJK) = BC_V_ID(IJK)
1509                       BC_W_ID(IPJK) = BC_W_ID(IJK)
1510     
1511                    END DO
1512                 END DO
1513              END DO
1514     
1515           ENDIF
1516     
1517     !     WEST BOUNDARY
1518           I = ISTART1
1519     
1520           IF(I==IMIN1) THEN
1521     
1522              DO I23 = ISTART3,ISTART2
1523                 DO K = KSTART3,KEND3
1524                    DO J = JSTART3, JEND3
1525     
1526                       IJK  = FUNIJK(I,J,K)
1527                       IMJK = FUNIJK(I23,J,K)
1528     
1529                       BLOCKED_CELL_AT(IMJK) = BLOCKED_CELL_AT(IJK)
1530                       BLOCKED_U_CELL_AT(IMJK) = BLOCKED_U_CELL_AT(IJK)
1531                       BLOCKED_V_CELL_AT(IMJK) = BLOCKED_V_CELL_AT(IJK)
1532                       BLOCKED_W_CELL_AT(IMJK) = BLOCKED_W_CELL_AT(IJK)
1533     
1534                       IF(BLOCKED_CELL_AT(IMJK)) FLAG(IMJK) = 100
1535     
1536                       VOL(IMJK)   = VOL(IJK)
1537                       VOL_U(IMJK) = VOL_U(IJK)
1538                       VOL_V(IMJK) = VOL_V(IJK)
1539                       VOL_W(IMJK) = VOL_W(IJK)
1540     
1541                       AYZ(IMJK)   = AYZ(IJK)
1542                       AYZ_U(IMJK) = AYZ_U(IJK)
1543                       AYZ_V(IMJK) = AYZ_V(IJK)
1544                       AYZ_W(IMJK) = AYZ_W(IJK)
1545     
1546                       AXY(IMJK)   = AXY(IJK)
1547                       AXY_U(IMJK) = AXY_U(IJK)
1548                       AXY_V(IMJK) = AXY_V(IJK)
1549                       AXY_W(IMJK) = AXY_W(IJK)
1550     
1551                       AXZ(IMJK)   = AXZ(IJK)
1552                       AXZ_U(IMJK) = AXZ_U(IJK)
1553                       AXZ_V(IMJK) = AXZ_V(IJK)
1554                       AXZ_W(IMJK) = AXZ_W(IJK)
1555     
1556                       BC_ID(IMJK)   = BC_ID(IJK)
1557                       BC_U_ID(IMJK) = BC_U_ID(IJK)
1558                       BC_V_ID(IMJK) = BC_V_ID(IJK)
1559                       BC_W_ID(IMJK) = BC_W_ID(IJK)
1560     
1561                    END DO
1562                 END DO
1563              END DO
1564     
1565           ENDIF
1566     
1567     !     NORTH BOUNDARY
1568           J = JEND1
1569     
1570           IF(J==JMAX1) THEN
1571     
1572              DO J23 = JEND2,JEND3
1573                 DO K = KSTART3,KEND3
1574                    DO I = ISTART3, IEND3
1575     
1576                       IJK  = FUNIJK(I,J,K)
1577                       IJPK = FUNIJK(I,J23,K)
1578     
1579                       BLOCKED_CELL_AT(IJPK) = BLOCKED_CELL_AT(IJK)
1580                       BLOCKED_U_CELL_AT(IJPK) = BLOCKED_U_CELL_AT(IJK)
1581                       BLOCKED_V_CELL_AT(IJPK) = BLOCKED_V_CELL_AT(IJK)
1582                       BLOCKED_W_CELL_AT(IJPK) = BLOCKED_W_CELL_AT(IJK)
1583     
1584                       IF(BLOCKED_CELL_AT(IJPK)) FLAG(IJPK) = 100
1585     
1586                       VOL(IJPK)   = VOL(IJK)
1587                       VOL_U(IJPK) = VOL_U(IJK)
1588                       VOL_V(IJPK) = VOL_V(IJK)
1589                       VOL_W(IJPK) = VOL_W(IJK)
1590     
1591                       AYZ(IJPK)   = AYZ(IJK)
1592                       AYZ_U(IJPK) = AYZ_U(IJK)
1593                       AYZ_V(IJPK) = AYZ_V(IJK)
1594                       AYZ_W(IJPK) = AYZ_W(IJK)
1595     
1596                       AXY(IJPK)   = AXY(IJK)
1597                       AXY_U(IJPK) = AXY_U(IJK)
1598                       AXY_V(IJPK) = AXY_V(IJK)
1599                       AXY_W(IJPK) = AXY_W(IJK)
1600     
1601                       AXZ(IJPK)   = AXZ(IJK)
1602                       AXZ_U(IJPK) = AXZ_U(IJK)
1603                       AXZ_V(IJPK) = AXZ_V(IJK)
1604                       AXZ_W(IJPK) = AXZ_W(IJK)
1605     
1606                       BC_ID(IJPK)   = BC_ID(IJK)
1607                       BC_U_ID(IJPK) = BC_U_ID(IJK)
1608                       BC_V_ID(IJPK) = BC_V_ID(IJK)
1609                       BC_W_ID(IJPK) = BC_W_ID(IJK)
1610     
1611                    END DO
1612                 END DO
1613              END DO
1614     
1615           ENDIF
1616     !     SOUTH BOUNDARY
1617           J = JSTART1
1618     
1619           IF(J==JMIN1) THEN
1620     
1621              DO J23 = JSTART3,JSTART2
1622                 DO K = KSTART3,KEND3
1623                    DO I = ISTART3, IEND3
1624     
1625                       IJK  = FUNIJK(I,J,K)
1626                       IJMK = FUNIJK(I,J23,K)
1627     
1628                       BLOCKED_CELL_AT(IJMK) = BLOCKED_CELL_AT(IJK)
1629                       BLOCKED_U_CELL_AT(IJMK) = BLOCKED_U_CELL_AT(IJK)
1630                       BLOCKED_V_CELL_AT(IJMK) = BLOCKED_V_CELL_AT(IJK)
1631                       BLOCKED_W_CELL_AT(IJMK) = BLOCKED_W_CELL_AT(IJK)
1632     
1633                       IF(BLOCKED_CELL_AT(IJMK)) FLAG(IJMK) = 100
1634     
1635                       VOL(IJMK)   = VOL(IJK)
1636                       VOL_U(IJMK) = VOL_U(IJK)
1637                       VOL_V(IJMK) = VOL_V(IJK)
1638                       VOL_W(IJMK) = VOL_W(IJK)
1639     
1640                       AYZ(IJMK)   = AYZ(IJK)
1641                       AYZ_U(IJMK) = AYZ_U(IJK)
1642                       AYZ_V(IJMK) = AYZ_V(IJK)
1643                       AYZ_W(IJMK) = AYZ_W(IJK)
1644     
1645                       AXY(IJMK)   = AXY(IJK)
1646                       AXY_U(IJMK) = AXY_U(IJK)
1647                       AXY_V(IJMK) = AXY_V(IJK)
1648                       AXY_W(IJMK) = AXY_W(IJK)
1649     
1650                       AXZ(IJMK)   = AXZ(IJK)
1651                       AXZ_U(IJMK) = AXZ_U(IJK)
1652                       AXZ_V(IJMK) = AXZ_V(IJK)
1653                       AXZ_W(IJMK) = AXZ_W(IJK)
1654     
1655                       BC_ID(IJMK)   = BC_ID(IJK)
1656                       BC_U_ID(IJMK) = BC_U_ID(IJK)
1657                       BC_V_ID(IJMK) = BC_V_ID(IJK)
1658                       BC_W_ID(IJMK) = BC_W_ID(IJK)
1659     
1660                    END DO
1661                 END DO
1662              END DO
1663     
1664           ENDIF
1665     
1666           IF(DO_K) THEN
1667     
1668     !        TOP BOUNDARY
1669              K = KEND1
1670     
1671              IF(K==KMAX1) THEN
1672     
1673                 DO K23=KEND2,KEND3
1674                    DO J = JSTART3,JEND3
1675                       DO I = ISTART3, IEND3
1676     
1677                          IJK  = FUNIJK(I,J,K)
1678                          IJKP = FUNIJK(I,J,K23)
1679     
1680                          BLOCKED_CELL_AT(IJKP) = BLOCKED_CELL_AT(IJK)
1681                          BLOCKED_U_CELL_AT(IJKP) = BLOCKED_U_CELL_AT(IJK)
1682                          BLOCKED_V_CELL_AT(IJKP) = BLOCKED_V_CELL_AT(IJK)
1683                          BLOCKED_W_CELL_AT(IJKP) = BLOCKED_W_CELL_AT(IJK)
1684     
1685                          IF(BLOCKED_CELL_AT(IJKP)) FLAG(IJKP) = 100
1686     
1687                          VOL(IJKP)   = VOL(IJK)
1688                          VOL_U(IJKP) = VOL_U(IJK)
1689                          VOL_V(IJKP) = VOL_V(IJK)
1690                          VOL_W(IJKP) = VOL_W(IJK)
1691     
1692                          AYZ(IJKP)   = AYZ(IJK)
1693                          AYZ_U(IJKP) = AYZ_U(IJK)
1694                          AYZ_V(IJKP) = AYZ_V(IJK)
1695                          AYZ_W(IJKP) = AYZ_W(IJK)
1696     
1697                          AXY(IJKP)   = AXY(IJK)
1698                          AXY_U(IJKP) = AXY_U(IJK)
1699                          AXY_V(IJKP) = AXY_V(IJK)
1700                          AXY_W(IJKP) = AXY_W(IJK)
1701     
1702                          AXZ(IJKP)   = AXZ(IJK)
1703                          AXZ_U(IJKP) = AXZ_U(IJK)
1704                          AXZ_V(IJKP) = AXZ_V(IJK)
1705                          AXZ_W(IJKP) = AXZ_W(IJK)
1706     
1707                          BC_ID(IJKP)   = BC_ID(IJK)
1708                          BC_U_ID(IJKP) = BC_U_ID(IJK)
1709                          BC_V_ID(IJKP) = BC_V_ID(IJK)
1710                          BC_W_ID(IJKP) = BC_W_ID(IJK)
1711     
1712                       END DO
1713                    END DO
1714                 END DO
1715     
1716              ENDIF
1717     
1718     !        BOTTOM BOUNDARY
1719              K = KSTART1
1720     
1721              IF(K==KMIN1) THEN
1722     
1723                 DO K23 = KSTART3,KSTART2
1724                    DO J = JSTART3,JEND3
1725                       DO I = ISTART3, IEND3
1726     
1727                          IJK  = FUNIJK(I,J,K)
1728                          IJKM = FUNIJK(I,J,K23)
1729     
1730                          BLOCKED_CELL_AT(IJKM) = BLOCKED_CELL_AT(IJK)
1731                          BLOCKED_U_CELL_AT(IJKM) = BLOCKED_U_CELL_AT(IJK)
1732                          BLOCKED_V_CELL_AT(IJKM) = BLOCKED_V_CELL_AT(IJK)
1733                          BLOCKED_W_CELL_AT(IJKM) = BLOCKED_W_CELL_AT(IJK)
1734     
1735                          IF(BLOCKED_CELL_AT(IJKM)) FLAG(IJKM) = 100
1736     
1737                          VOL(IJKM)   = VOL(IJK)
1738                          VOL_U(IJKM) = VOL_U(IJK)
1739                          VOL_V(IJKM) = VOL_V(IJK)
1740                          VOL_W(IJKM) = VOL_W(IJK)
1741     
1742                          AYZ(IJKM)   = AYZ(IJK)
1743                          AYZ_U(IJKM) = AYZ_U(IJK)
1744                          AYZ_V(IJKM) = AYZ_V(IJK)
1745                          AYZ_W(IJKM) = AYZ_W(IJK)
1746     
1747                          AXY(IJKM)   = AXY(IJK)
1748                          AXY_U(IJKM) = AXY_U(IJK)
1749                          AXY_V(IJKM) = AXY_V(IJK)
1750                          AXY_W(IJKM) = AXY_W(IJK)
1751     
1752                          AXZ(IJKM)   = AXZ(IJK)
1753                          AXZ_U(IJKM) = AXZ_U(IJK)
1754                          AXZ_V(IJKM) = AXZ_V(IJK)
1755                          AXZ_W(IJKM) = AXZ_W(IJK)
1756     
1757                          BC_ID(IJKM)   = BC_ID(IJK)
1758                          BC_U_ID(IJKM) = BC_U_ID(IJK)
1759                          BC_V_ID(IJKM) = BC_V_ID(IJK)
1760                          BC_W_ID(IJKM) = BC_W_ID(IJK)
1761     
1762                       END DO
1763                    END DO
1764                 END DO
1765     
1766              ENDIF
1767     
1768           ENDIF
1769     
1770           DO IJK = IJKSTART3, IJKEND3
1771     
1772              IF(BLOCKED_U_CELL_AT(IJK)) FLAG_E(IJK)=100
1773              IF(BLOCKED_V_CELL_AT(IJK)) FLAG_N(IJK)=100
1774              IF(BLOCKED_W_CELL_AT(IJK)) FLAG_T(IJK)=100
1775     
1776           ENDDO
1777     
1778     !     BLOCKED CELLS WERE ASSIGNED THE FLAG 100 DURING PRE_PROCESSING
1779     !     THIS IS INCORRECT FOR FREEE-SLIP AND PARTIAL-SLIP
1780     !     THE FLAG IS OVERWRITTEN HERE TO ACCOUNT FOR NSW,FSW AND PSW BCs
1781     
1782           DO IJK = IJKSTART3, IJKEND3
1783     
1784              IF(FLAG(IJK)==100) THEN
1785     
1786                 IMJK   = IM_OF(IJK)
1787                 IJMK   = JM_OF(IJK)
1788                 IJKM   = KM_OF(IJK)
1789     
1790                 IPJK   = IP_OF(IJK)
1791                 IJPK   = JP_OF(IJK)
1792                 IJKP   = KP_OF(IJK)
1793     
1794                 IF(CUT_CELL_AT(IMJK)) THEN
1795                    BCID=BC_ID(IMJK)
1796                    IF(BCID>0) THEN
1797                       IF(IS_NSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 100
1798                       IF(IS_FSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 101
1799                       IF(IS_PSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 102
1800                    ENDIF
1801     
1802                 ELSEIF(CUT_CELL_AT(IPJK)) THEN
1803                    BCID=BC_ID(IPJK)
1804                    IF(BCID>0) THEN
1805                       IF(IS_NSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 100
1806                       IF(IS_FSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 101
1807                       IF(IS_PSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 102
1808                    ENDIF
1809     
1810                 ELSEIF(CUT_CELL_AT(IJMK)) THEN
1811                    BCID=BC_ID(IJMK)
1812                    IF(BCID>0) THEN
1813                       IF(IS_NSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 100
1814                       IF(IS_FSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 101
1815                       IF(IS_PSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 102
1816                    ENDIF
1817     
1818                 ELSEIF(CUT_CELL_AT(IJPK)) THEN
1819                    BCID=BC_ID(IJPK)
1820                    IF(BCID>0) THEN
1821                       IF(IS_NSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 100
1822                       IF(IS_FSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 101
1823                       IF(IS_PSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 102
1824                    ENDIF
1825     
1826                 ELSEIF(CUT_CELL_AT(IJKM)) THEN
1827                    BCID=BC_ID(IJKM)
1828                    IF(BCID>0) THEN
1829                       IF(IS_NSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 100
1830                       IF(IS_FSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 101
1831                       IF(IS_PSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 102
1832                    ENDIF
1833     
1834                 ELSEIF(CUT_CELL_AT(IJKP)) THEN
1835                    BCID=BC_ID(IJKP)
1836                    IF(BCID>0) THEN
1837                       IF(IS_NSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 100
1838                       IF(IS_FSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 101
1839                       IF(IS_PSW(BC_TYPE_ENUM(BCID))) FLAG(IJK) = 102
1840                    ENDIF
1841     
1842                 ENDIF
1843              ENDIF
1844           ENDDO
1845     
1846           call send_recv(FLAG,2)
1847           RETURN
1848     
1849           END SUBROUTINE SET_GHOST_CELL_FLAGS
1850     
1851     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1852     !                                                                      C
1853     !  Module name: GET_POTENTIAL_CUT_CELLS                                  C
1854     !  Purpose: Set flags for scalar cut cells, based on intersection      C
1855     !  of the grid with the quadric(s)                                     C
1856     !                                                                      C
1857     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1858     !  Reviewer:                                          Date:            C
1859     !                                                                      C
1860     !  Revision Number #                                  Date: ##-###-##  C
1861     !  Author: #                                                           C
1862     !  Purpose: #                                                          C
1863     !                                                                      C
1864     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1865       SUBROUTINE GET_POTENTIAL_CUT_CELLS
1866     
1867           USE compar, only: ijkstart3, ijkend3, mype, pe_io
1868           USE cutcell
1869           USE functions, only: funijk, bottom_of, south_of, west_of
1870           USE geometry, ONLY: dx, dy, dz, do_k, imin3, imax3, jmin3, jmax3, kmin3, kmax3, flag, axy, axz, ayz, vol, no_k
1871           USE indices, only: i_of, j_of, k_of
1872           USE polygon, ONLY: n_polygon
1873           USE quadric, ONLY: tol_f
1874           USE param, ONLY: DIMENSION_3
1875           USE param1, ONLY: HALF, ZERO
1876     
1877           IMPLICIT NONE
1878           INTEGER :: IJK,I,J,K,II,JJ,KK
1879           INTEGER :: I1,I2,J1,J2,K1,K2
1880           INTEGER :: NODE,N_N1,N_N2
1881           LOGICAL :: ALL_NEGATIVE,ALL_POSITIVE,CLIP_FLAG
1882           INTEGER :: Q_ID,BCID
1883     
1884           INTEGER :: IJK_NB,NUMBER_OF_POTENTIAL_CUT_CELLS
1885     
1886           DOUBLE PRECISION :: xc,yc,zc,fc
1887     
1888           LOGICAL, DIMENSION(DIMENSION_3) ::POSITIVE_F_AT
1889     
1890           POTENTIAL_CUT_CELL_AT=.TRUE.
1891     
1892           RETURN  ! This subroutine is currently disabled
1893     
1894           IF(MyPE == PE_IO) THEN
1895              WRITE(*,10)'ESTIMATING POTENTIAL SCALAR CUT CELLS...'
1896           ENDIF
1897     10    FORMAT(1X,A)
1898     
1899     !======================================================================
1900     !  Evaluate f at cell center and store where f>0
1901     !======================================================================
1902     
1903           DO IJK = IJKSTART3, IJKEND3
1904     
1905              I = I_OF(IJK)
1906              J = J_OF(IJK)
1907              K = K_OF(IJK)
1908     
1909              xc = XG_E(I) - HALF*DX(I)
1910              yc = YG_N(J) - HALF*DY(J)
1911     
1912              IF(DO_K) THEN
1913                 zc = ZG_T(K) - HALF*DZ(K)
1914              ELSE
1915                 zc = zero
1916              ENDIF
1917     
1918               Q_ID = 1
1919               CALL EVAL_F('QUADRIC',xc,yc,zc,Q_ID,fc,CLIP_FLAG)
1920     
1921               CALL EVAL_F('POLYGON',xc,yc,zc,N_POLYGON,fc,CLIP_FLAG)
1922     
1923               CALL EVAL_F('USR_DEF',xc,yc,zc,N_USR_DEF,Fc,CLIP_FLAG)
1924     
1925               X_NODE(15) = xc
1926               Y_NODE(15) = yc
1927               Z_NODE(15) = zc
1928               CALL EVAL_STL_FCT_AT('SCALAR',IJK,15,fc,CLIP_FLAG,BCID)
1929     
1930               IF(fc>TOL_F) THEN
1931                  POSITIVE_F_AT(IJK)=.TRUE.
1932               ELSE
1933                  POSITIVE_F_AT(IJK)=.FALSE.
1934               ENDIF
1935     
1936            ENDDO
1937     
1938           DO IJK = IJKSTART3, IJKEND3
1939     
1940              IF(INTERIOR_CELL_AT(IJK)) THEN
1941     
1942                 I = I_OF(IJK)
1943                 J = J_OF(IJK)
1944                 K = K_OF(IJK)
1945     
1946                 I1 = MAX(I - 2,IMIN3)
1947                 I2 = MIN(I + 2,IMAX3)
1948                 J1 = MAX(J - 2,JMIN3)
1949                 J2 = MIN(J + 2,JMAX3)
1950     
1951                 IF(DO_K) THEN
1952                    K1 = MAX(K - 2,KMIN3)
1953                    K2 = MIN(K + 2,KMAX3)
1954                 ELSE
1955                    K1=K
1956                    K2=K
1957                 ENDIF
1958     
1959                 IF(POSITIVE_F_AT(IJK)) THEN
1960                    ALL_POSITIVE=.TRUE.
1961                    DO KK=K1,K2
1962                       DO JJ=J1,J2
1963                          DO II=I1,I2
1964                             IJK_NB = FUNIJK(II,JJ,KK)
1965                             IF(.NOT.POSITIVE_F_AT(IJK_NB)) THEN
1966                                ALL_POSITIVE=.FALSE.
1967                             ENDIF
1968                          ENDDO
1969                       ENDDO
1970                     ENDDO
1971     
1972                     IF(ALL_POSITIVE) THEN
1973                        POTENTIAL_CUT_CELL_AT(IJK)=.FALSE.
1974                        FLAG(IJK) = 100
1975                        BLOCKED_CELL_AT(IJK) = .TRUE.               ! Blocked fluid cell
1976                        STANDARD_CELL_AT(IJK) = .FALSE.
1977                        AXY(IJK) = ZERO
1978                        AXZ(IJK) = ZERO
1979                        AYZ(IJK) = ZERO
1980                        VOL(IJK) = ZERO
1981     
1982                        AXY(BOTTOM_OF(IJK)) = ZERO
1983                        AXZ(SOUTH_OF(IJK)) = ZERO
1984                        AYZ(WEST_OF(IJK)) = ZERO
1985                     ENDIF
1986     
1987                 ELSE
1988                    ALL_NEGATIVE=.TRUE.
1989                    DO KK=K1,K2
1990                       DO JJ=J1,J2
1991                          DO II=I1,I2
1992                             IJK_NB = FUNIJK(II,JJ,KK)
1993                             IF(POSITIVE_F_AT(IJK_NB)) THEN
1994                                ALL_NEGATIVE=.FALSE.
1995                             ENDIF
1996                          ENDDO
1997                       ENDDO
1998                     ENDDO
1999     
2000                     IF(ALL_NEGATIVE) THEN
2001                        POTENTIAL_CUT_CELL_AT(IJK)=.FALSE.
2002                        BLOCKED_CELL_AT(IJK) = .FALSE.
2003                        STANDARD_CELL_AT(IJK) = .TRUE.           ! Regular fluid cell
2004                     ENDIF
2005     
2006                 ENDIF
2007     
2008                 IF((FLAG(IJK)>=100).AND.(FLAG(IJK)<=102)) THEN
2009                    POTENTIAL_CUT_CELL_AT(IJK)=.FALSE.
2010                    BLOCKED_CELL_AT(IJK) = .TRUE.
2011                    STANDARD_CELL_AT(IJK) = .FALSE.          ! Blocked cell = wall cell
2012                 ENDIF
2013     
2014              ENDIF
2015     
2016           END DO
2017     
2018           NUMBER_OF_POTENTIAL_CUT_CELLS = 0
2019     
2020           IF(NO_K) THEN
2021              N_N1 = 5
2022              N_N2 = 8
2023           ELSE
2024              N_N1 = 1
2025              N_N2 = 8
2026           ENDIF
2027     
2028           DO IJK=IJKSTART3,IJKEND3
2029              IF(POTENTIAL_CUT_CELL_AT(IJK)) THEN
2030                 NUMBER_OF_POTENTIAL_CUT_CELLS = NUMBER_OF_POTENTIAL_CUT_CELLS + 1
2031              ELSE
2032                 CALL GET_CELL_NODE_COORDINATES(IJK,'SCALAR')
2033     
2034                 IF(NO_K) THEN
2035                    NUMBER_OF_NODES(IJK) = 4
2036                    CONNECTIVITY(IJK,1) = IJK_OF_NODE(5)
2037                    CONNECTIVITY(IJK,2) = IJK_OF_NODE(6)
2038                    CONNECTIVITY(IJK,3) = IJK_OF_NODE(8)
2039                    CONNECTIVITY(IJK,4) = IJK_OF_NODE(7)
2040                 ELSE
2041                    NUMBER_OF_NODES(IJK) = 8
2042                    DO NODE = N_N1,N_N2
2043                       CONNECTIVITY(IJK,NODE) = IJK_OF_NODE(NODE)
2044                    END DO
2045                 ENDIF
2046     
2047                 X_U(IJK) = X_NODE(8)
2048                 Y_U(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
2049                 Z_U(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
2050     
2051                 X_V(IJK) = HALF * (X_NODE(7) + X_NODE(8))
2052                 Y_V(IJK) = Y_NODE(8)
2053                 Z_V(IJK) = HALF * (Z_NODE(4) + Z_NODE(8))
2054     
2055                 X_W(IJK) = HALF * (X_NODE(7) + X_NODE(8))
2056                 Y_W(IJK) = HALF * (Y_NODE(6) + Y_NODE(8))
2057                 Z_W(IJK) = Z_NODE(8)
2058     
2059              ENDIF
2060           ENDDO
2061     
2062     !      call SEND_RECEIVE_1D_LOGICAL(SNAP,2)
2063           IF(MyPE == PE_IO) THEN
2064              WRITE(*,*)'DONE ESTIMATING POTENTIAL SCALAR CUT CELLS.',NUMBER_OF_POTENTIAL_CUT_CELLS,IJKEND3
2065           ENDIF
2066     
2067           RETURN
2068           END SUBROUTINE GET_POTENTIAL_CUT_CELLS
2069