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