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