File: /nfs/home/0/users/jenkins/mfix.git/model/cartesian_grid/cut_cell_preprocessing.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CUT_CELL_PREPROCESSING                                 C
4     !  Purpose: Perform the cut-cell preprocessing stage:                  C
5     !           Identify cut cells, define face areas, and volumes         C
6     !           Set flags                                                  C
7     !           Compute Interpolations factors                             C
8     !           Compute Non-orthogonality Corrections terms                C
9     !                                                                      C
10     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
11     !  Reviewer:                                          Date:            C
12     !                                                                      C
13     !  Revision Number #                                  Date: ##-###-##  C
14     !  Author: #                                                           C
15     !  Purpose: #                                                          C
16     !                                                                      C
17     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
18       SUBROUTINE CUT_CELL_PREPROCESSING
19     
20           USE param
21           USE param1
22           USE parallel
23           USE constant
24           USE run
25           USE toleranc
26           USE geometry
27           USE indices
28           USE compar
29           USE sendrecv
30           USE quadric
31           USE cutcell
32           USE vtk
33           use cdist
34           USE fldvar
35           USE functions
36     
37           IMPLICIT NONE
38     
39           INTEGER :: SAFE_MODE_COUNT
40           DOUBLE PRECISION :: CPU_PP_START,CPU_PP_END
41     
42           IF(.NOT.CG_HEADER_WAS_PRINTED) CALL PRINT_CG_HEADER
43     
44           CALL CPU_TIME (CPU_PP_START)
45     
46           CALL OPEN_CUT_CELL_FILES
47     
48           CALL ALLOCATE_CUT_CELL_ARRAYS
49     
50           CALL DEFINE_QUADRICS
51     
52           CALL SET_3D_CUT_CELL_FLAGS
53     
54           CALL GATHER_DATA
55     
56     !======================================================================
57     ! Gather Data  and writes surface(s) defined by all cut cells
58     !======================================================================
59     
60           IF(WRITE_VTK_FILES.AND.(.NOT.BDIST_IO)) THEN
61              CALL WRITE_CUT_SURFACE_VTK
62           ENDIF
63     
64     
65     
66           CALL SET_3D_CUT_U_CELL_FLAGS
67           CALL SET_3D_CUT_V_CELL_FLAGS
68           IF(DO_K) CALL SET_3D_CUT_W_CELL_FLAGS
69     
70           CALL SET_3D_CUT_CELL_TREATMENT_FLAGS
71     
72           CALL GET_3D_ALPHA_U_CUT_CELL
73           CALL GET_3D_ALPHA_V_CUT_CELL
74           IF(DO_K) CALL GET_3D_ALPHA_W_CUT_CELL
75     
76           CALL SET_GHOST_CELL_FLAGS
77     
78           CALL SET_ODXYZ_U_CUT_CELL
79           CALL SET_ODXYZ_V_CUT_CELL
80           IF(DO_K) CALL SET_ODXYZ_W_CUT_CELL
81     
82           CALL GET_U_MASTER_CELLS
83           CALL GET_V_MASTER_CELLS
84           IF(DO_K) CALL GET_W_MASTER_CELLS
85     
86           CALL SEND_RECEIVE_CUT_CELL_VARIABLES
87     
88           CALL GET_DISTANCE_TO_WALL
89     
90           CALL PRINT_GRID_STATISTICS
91     
92           CALL CG_GET_BC_AREA
93     
94           CALL FLOW_TO_VEL(.FALSE.)
95     
96           CALL CG_FLOW_TO_VEL
97     
98           CALL CONVERT_CG_MI_TO_PS
99     
100     
101           CALL CPU_TIME (CPU_PP_END)
102     
103     
104     
105           IF(myPE == PE_IO) THEN
106              WRITE(*,20)'CARTESIAN GRID PRE-PROCESSING COMPLETED IN ',CPU_PP_END - CPU_PP_START, ' SECONDS.'
107              WRITE(*,10)'============================================================================='
108           ENDIF
109     
110           IF(myPE == PE_IO) THEN
111     
112              SAFE_MODE_COUNT = SUM(CG_SAFE_MODE)
113     
114              IF(SAFE_MODE_COUNT>0) THEN
115     
116     
117                 WRITE(*,10)'######################################################################'
118                 WRITE(*,10)'######################################################################'
119                 WRITE(*,10)'##                                                                  ##'
120                 WRITE(*,10)'##                              ||                                  ##'
121                 WRITE(*,10)'##                              ||                                  ##'
122                 WRITE(*,10)'##                              \/                                  ##'
123                 WRITE(*,10)'##                                                                  ##'
124                 WRITE(*,10)'##  ===>   WARNING: RUNNING CARTESIAN GRID IN SAFE MODE !  <===     ##'
125                 WRITE(*,10)'##                                                                  ##'
126                 WRITE(*,10)'##  SAFE MODE ACTIVATED FOR :                                       ##'
127                 IF(CG_SAFE_MODE(1)==1) WRITE(*,10)'##                            - All scalar quantities               ##'
128                 IF(CG_SAFE_MODE(3)==1) WRITE(*,10)'##                            - X-Velocity (Gas and Solids)         ##'
129                 IF(CG_SAFE_MODE(4)==1) WRITE(*,10)'##                            - Y-Velocity (Gas and Solids)         ##'
130                 IF(CG_SAFE_MODE(5)==1) WRITE(*,10)'##                            - Z-Velocity (Gas and Solids)         ##'
131                 WRITE(*,10)'##                                                                  ##'
132                 WRITE(*,10)'##                              /\                                  ##'
133                 WRITE(*,10)'##                              ||                                  ##'
134                 WRITE(*,10)'##                              ||                                  ##'
135                 WRITE(*,10)'##                                                                  ##'
136                 WRITE(*,10)'######################################################################'
137                 WRITE(*,10)'######################################################################'
138     
139     
140     
141              ENDIF
142           ENDIF
143     
144     
145           RETURN
146     
147     10    FORMAT(1X,A)
148     20    FORMAT(1X,A,F8.2,A)
149     
150           END SUBROUTINE CUT_CELL_PREPROCESSING
151     
152     
153     
154     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
155     !                                                                      C
156     !  Module name: PRINT_CG_HEADER                                        C
157     !  Purpose: Display Cartesian-Grid Header on screen                    C
158     !                                                                      C
159     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
160     !  Reviewer:                                          Date:            C
161     !                                                                      C
162     !  Revision Number #                                  Date: ##-###-##  C
163     !  Author: #                                                           C
164     !  Purpose: #                                                          C
165     !                                                                      C
166     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
167       SUBROUTINE PRINT_CG_HEADER
168     
169           USE param
170           USE param1
171           USE parallel
172           USE compar
173           USE cutcell
174     
175           IMPLICIT NONE
176     
177     
178     
179     
180           IF(myPE == PE_IO) THEN
181     
182              WRITE(*,10)'============================================================================='
183              WRITE(*,10)'  ____          ___________    _   ____      ____    __________   __________  '
184              WRITE(*,10)' |    \        /           |  (_)  \   \    /   /   |          | |          | '
185              WRITE(*,10)' |     \      /      ______|  ___   \   \  /   /    |    ______| |    ______| '
186              WRITE(*,10)' |      \    /      |______  |   |   \   \/   /     |   |        |   |        '
187              WRITE(*,10)' |       \  /              | |   |    \      /  === |   |        |   |  ____  '
188              WRITE(*,10)' |   |\   \/   /|    ______| |   |    /      \      |   |        |   | |_   | '
189              WRITE(*,10)' |   | \      / |   |        |   |   /   /\   \     |   |______  |   |___|  | '
190              WRITE(*,10)' |   |  \    /  |   |        |   |  /   /  \   \    |          | |          | '
191              WRITE(*,10)' |___|   \__/   |___|        |___| /___/    \___\   |__________| |__________| '
192              WRITE(*,10)'                                                                              '
193              WRITE(*,10)'============================================================================='
194              WRITE(*,10)'MFIX WITH CARTESIAN GRID IMPLEMENTATION.'
195     
196     
197              IF(RE_INDEXING) THEN
198                 WRITE(*,10)'RE-INDEXING IS TURNED ON.'
199     !            IF(ADJUST_PROC_DOMAIN_SIZE) THEN
200     !               WRITE(*,10)'EACH PROCESSOR DOMAIN SIZE WILL BE ADJUSTED FOR BETTER LOAD BALANCING.'
201     !            ELSE
202     !               WRITE(*,10)'WARNING: PROCESSOR DOMAIN SIZE WILL BE UNIFORMLY DISTIBUTED.'
203     !               WRITE(*,10)'THIS COULD RESULT IN VERY POOR LOAD BALANCING.'
204     !            ENDIF
205              ELSE
206                 WRITE(*,10)'RE-INDEXING IS TURNED OFF.'
207              ENDIF
208     
209     
210     
211           ENDIF
212     
213           CG_HEADER_WAS_PRINTED = .TRUE.
214     
215     10    FORMAT(1X,A)
216     
217           RETURN
218     
219           END SUBROUTINE PRINT_CG_HEADER
220     
221     
222     
223     
224     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
225     !                                                                      C
226     !  Module name: EVAL_F                                                 C
227     !  Purpose: Evaluate the function f(x,y,z) defining the boundary       C
228     !                                                                      C
229     !  Author: Jeff Dietiker                               Date: 21-FEB-08 C
230     !  Reviewer:                                           Date:           C
231     !                                                                      C
232     !                                                                      C
233     !  Literature/Document References:                                     C
234     !                                                                      C
235     !  Variables referenced:                                               C
236     !  Variables modified:                                                 C
237     !                                                                      C
238     !  Local variables:                                                    C
239     !                                                                      C
240     !  Modified: ##                                        Date: ##-###-## C
241     !  Purpose: ##                                                         C
242     !                                                                      C
243     !  Literature/Document References: ##                                  C
244     !                                                                      C
245     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
246             SUBROUTINE EVAL_F(METHOD,x1,x2,x3,Q,f,CLIP_FLAG)
247     
248           USE parallel
249           USE compar
250           USE sendrecv
251           USE quadric
252           USE cutcell
253     
254           USE quadric
255     
256     
257           IMPLICIT NONE
258     
259           DOUBLE PRECISION x1,x2,x3
260           DOUBLE PRECISION f
261           INTEGER :: Q,Q_ID,BCID
262           LOGICAL :: CLIP_FLAG
263           CHARACTER (LEN = 7) :: METHOD
264           CHARACTER(LEN=9) :: GR
265     
266           INTEGER :: GROUP,GS,P
267     
268           DOUBLE PRECISION,DIMENSION(DIM_GROUP,0:DIM_QUADRIC) :: F_G
269     
270           SELECT CASE(METHOD)
271     
272              CASE('QUADRIC')
273     
274                 F_G = - UNDEFINED
275     
276                 DO GROUP = 1, N_GROUP
277                    GS = GROUP_SIZE(GROUP)
278                    GR = TRIM(GROUP_RELATION(GROUP))
279     
280                    DO P = 1 , GS
281                       Q_ID = GROUP_Q(GROUP,P)
282                       CALL GET_F_QUADRIC(x1,x2,x3,Q_ID,F_G(GROUP,P),CLIP_FLAG)
283                    ENDDO
284                    IF(GR == 'AND') THEN
285                       F_G(GROUP,0) = MAXVAL(F_G(GROUP,1:GS))
286                    ELSEIF(GR == 'OR') THEN
287                       F_G(GROUP,0) = MINVAL(F_G(GROUP,1:GS))
288                    ELSEIF(GR == 'PIECEWISE') THEN
289                       CALL REASSSIGN_QUADRIC(x1,x2,x3,GROUP,Q_ID)
290     !                  CLIP_FLAG=.FALSE.
291                       CALL GET_F_QUADRIC(x1,x2,x3,Q_ID,F_G(GROUP,0),CLIP_FLAG)
292     !                  CLIP_FLAG=.TRUE.
293                    ENDIF
294     
295                 ENDDO
296     
297                 f = F_G(1,0)
298     
299                 DO GROUP = 2, N_GROUP
300     
301                    GR = TRIM(RELATION_WITH_PREVIOUS(GROUP))
302     
303                    IF(GR =='AND') THEN
304                       f = DMAX1(f,F_G(GROUP,0))
305                    ELSEIF(GR =='OR') THEN
306                       f = DMIN1(f,F_G(GROUP,0))
307                    ENDIF
308     
309                 ENDDO
310     
311     
312              CASE('POLYGON')
313     
314                 CALL EVAL_POLY_FCT(x1,x2,x3,Q,f,CLIP_FLAG,BCID)
315     
316              CASE('USR_DEF')
317     
318                 CALL EVAL_USR_FCT(x1,x2,x3,Q,f,CLIP_FLAG)
319     
320     !         CASE('STL')
321     
322     !            CALL EVAL_STL_FCT(x1,x2,x3,Q,f,CLIP_FLAG,BCID)
323     
324              CASE DEFAULT
325     
326                 WRITE(*,*)'ERROR IN SUBROUTINE EVAL_F.'
327                 WRITE(*,*)'UNKNOWN METHOD:',METHOD
328                 WRITE(*,*)'ACCEPTABLE METHODS:'
329                 WRITE(*,*)'QUADRIC'
330                 WRITE(*,*)'POLYGON'
331                 WRITE(*,*)'USR_DEF'
332     !            WRITE(*,*)'STL'
333                 CALL MFIX_EXIT(myPE)
334     
335           END SELECT
336     
337     
338     
339           RETURN
340           END SUBROUTINE EVAL_F
341     
342     
343     
344     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
345     !                                                                      C
346     !  Module name: intersect_line                                         C
347     !  Purpose: Finds the intersection between the quadric surface ,       C
348     !           and the line (xa,ya,za) and (xb,yb,zb).                    C
349     !                                                                      C
350     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
351     !  Reviewer:                                          Date:            C
352     !                                                                      C
353     !  Revision Number #                                  Date: ##-###-##  C
354     !  Author: #                                                           C
355     !  Purpose: #                                                          C
356     !                                                                      C
357     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
358       SUBROUTINE INTERSECT_LINE(METHOD,xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_FLAG,xc,yc,zc)
359     
360           USE param
361           USE param1
362           USE parallel
363           USE constant
364           USE run
365           USE toleranc
366           USE geometry
367           USE indices
368           USE compar
369           USE sendrecv
370           USE quadric
371     
372           IMPLICIT NONE
373           DOUBLE PRECISION:: x1,y1,z1,x2,y2,z2,x3,y3,z3
374           DOUBLE PRECISION:: xa,ya,za,xb,yb,zb,xc,yc,zc
375           INTEGER :: Q_ID,niter
376           DOUBLE PRECISION :: x_intersection
377           DOUBLE PRECISION :: f1,f2,f3,fa,fb
378           DOUBLE PRECISION :: t1,t2,t3
379           LOGICAL :: CLIP_FLAG,CLIP_FLAG1,CLIP_FLAG2,CLIP_FLAG3,INTERSECT_FLAG
380           CHARACTER (LEN=7) ::METHOD
381     
382           x1 = xa    ! Initial guesses
383           y1 = ya
384           z1 = za
385           t1 = ZERO
386     
387           x2 = xb
388           y2 = yb
389           z2 = zb
390           t2 = ONE
391     
392           CALL EVAL_F(METHOD,x1,y1,z1,Q_ID,f1,CLIP_FLAG1)
393           CALL EVAL_F(METHOD,x2,y2,z2,Q_ID,f2,CLIP_FLAG2)
394     
395     !======================================================================
396     !  The line from (x1,y1,z1) and (x2,y2,z2) is parametrized
397     !  from t1 = ZERO to t2 = ONE
398     !======================================================================
399     
400           niter = 1
401     
402           CLIP_FLAG = (CLIP_FLAG1).AND.(CLIP_FLAG2)
403     
404           if(DABS(f1)<TOL_F) then  ! ignore intersection at corner
405             xc = UNDEFINED
406             yc = UNDEFINED
407             zc = UNDEFINED
408             INTERSECT_FLAG = .FALSE.
409           elseif(DABS(f2)<TOL_F) then ! ignore intersection at corner
410             xc = UNDEFINED
411             yc = UNDEFINED
412             zc = UNDEFINED
413             INTERSECT_FLAG = .FALSE.
414           elseif(f1*f2 < ZERO) then
415             niter = 0
416             f3 = 2.0d0*TOL_F
417             do while (   (abs(f3) > TOL_F)   .AND.   (niter<ITERMAX_INT)       )
418     
419               t3 = t1 - f1*(t2-t1)/(f2-f1)
420     
421               x3 = x1 + t3 * (x2 - x1)
422               y3 = y1 + t3 * (y2 - y1)
423               z3 = z1 + t3 * (z2 - z1)
424     
425               CALL EVAL_F(METHOD,x3,y3,z3,Q_ID,f3,CLIP_FLAG3)
426               if(f1*f3<0) then
427                 t2 = t3
428                 f2 = f3
429               else
430                 t1 = t3
431                 f1 = f3
432               endif
433               niter = niter + 1
434     
435             end do
436             if (niter < ITERMAX_INT) then
437                xc = x3
438                yc = y3
439                zc = z3
440               INTERSECT_FLAG = .TRUE.
441             else
442                WRITE(*,*)'   Subroutine intersect_line:'
443                WRITE(*,*)   'Unable to find the intersection of quadric:',Q_ID
444                WRITE(*,1000)'between (x1,y1,z1)= ', xa,ya,za
445                WRITE(*,1000)'   and  (x2,y2,z2)= ', xb,yb,zb
446                CALL EVAL_F(METHOD,xa,ya,za,Q_ID,fa,CLIP_FLAG1)
447                CALL EVAL_F(METHOD,xb,yb,zb,Q_ID,fb,CLIP_FLAG1)
448                WRITE(*,1000)'f(x1,y1,z1) = ', fa
449                WRITE(*,1000)'f(x2,y2,z2) = ', fb
450                WRITE(*,1000)'Current Location (x3,y3,z3)= ', x3,y3,z3
451                WRITE(*,1000)'Current value of abs(f) = ', DABS(f3)
452                WRITE(*,1000)'Tolerance = ', TOL_F
453                WRITE(*,*)'Maximum number of iterations = ', ITERMAX_INT
454                WRITE(*,*)   'Please increase the intersection tolerance, '
455                WRITE(*,*)   'or the maximum number of iterations, and try again.'
456                WRITE(*,*)   'MFiX will exit now.'
457                CALL MFIX_EXIT(myPE)
458                x_intersection = UNDEFINED
459                INTERSECT_FLAG = .FALSE.
460     
461             endif
462           else
463             xc = UNDEFINED
464             yc = UNDEFINED
465             zc = UNDEFINED
466             INTERSECT_FLAG = .FALSE.
467           endif
468     
469      1000 FORMAT(A,3(2X,G12.5))
470     
471     
472           RETURN
473     
474           END SUBROUTINE INTERSECT_LINE
475     
476     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
477     !                                                                      C
478     !  Module name: INTERSECT                                              C
479     !  Purpose: Intersects quadric with grid                               C
480     !                                                                      C
481     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
482     !  Reviewer:                                          Date:            C
483     !                                                                      C
484     !  Revision Number #                                  Date: ##-###-##  C
485     !  Author: #                                                           C
486     !  Purpose: #                                                          C
487     !                                                                      C
488     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
489       SUBROUTINE INTERSECT(IJK,TYPE_OF_CELL,Xi,Yi,Zi)
490     
491           USE param
492           USE param1
493           USE parallel
494           USE constant
495           USE run
496           USE toleranc
497           USE geometry
498           USE indices
499           USE compar
500           USE sendrecv
501           USE quadric
502           USE cutcell
503           USE polygon
504           USE functions
505     
506           IMPLICIT NONE
507           CHARACTER (LEN=*) :: TYPE_OF_CELL
508           INTEGER :: IJK,I,J,K,Q_ID,N_int_x,N_int_y,N_int_z,N_USR
509           DOUBLE PRECISION :: xa,ya,za,xb,yb,zb,xc,yc,zc
510           DOUBLE PRECISION :: Xi,Yi,Zi,Xc_backup,Yc_backup,Zc_backup
511           LOGICAL :: INTERSECT_FLAG
512     
513           Xi = UNDEFINED
514           Yi = UNDEFINED
515           Zi = UNDEFINED
516     
517     !======================================================================
518     !  Get coordinates of eight nodes
519     !======================================================================
520     
521           CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
522     
523     !======================================================================
524     !  Intersection with Edge 7 (node 7-8, Face North-Top):
525     !======================================================================
526           xa = X_NODE(7)
527           ya = Y_NODE(7)
528           za = Z_NODE(7)
529     
530           xb = X_NODE(8)
531           yb = Y_NODE(8)
532           zb = Z_NODE(8)
533     
534     
535           N_int_x = 0
536           INTERSECT_X(IJK) = .FALSE.
537           Q_ID = 1
538           CALL INTERSECT_LINE('QUADRIC',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_FLAG,xc,yc,zc)
539           IF ((INTERSECT_FLAG).AND.(xc/=Xi)) THEN
540              N_int_x = N_int_x + 1
541              INTERSECT_X(IJK) = .TRUE.
542              xc_backup = Xi
543              Xi = xc
544           ENDIF
545     
546           IF(N_int_x /= 1) THEN
547              Xi = UNDEFINED
548              INTERSECT_X(IJK) = .FALSE.
549           ENDIF
550     
551     
552           DO Q_ID = 1, N_POLYGON
553              CALL INTERSECT_LINE('POLYGON',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_X(IJK),xc,yc,zc)
554              IF(INTERSECT_X(IJK)) Xi = xc
555           ENDDO
556     
557           DO N_USR= 1, N_USR_DEF
558              CALL INTERSECT_LINE('USR_DEF',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_X(IJK),xc,yc,zc)
559              IF(INTERSECT_X(IJK)) Xi = xc
560           ENDDO
561     
562     !      IF(USE_STL) THEN
563     !         CALL INTERSECT_LINE_WITH_STL(xa,ya,za,xb,yb,zb,INTERSECT_X(IJK),xc,yc,zc)
564     !         IF(INTERSECT_X(IJK)) Xi = xc
565     !      ENDIF
566     
567           IF(TYPE_OF_CELL=='U_MOMENTUM') THEN
568              IF(SNAP(IJK)) THEN
569                 INTERSECT_X(IJK) = .TRUE.
570                 I = I_OF(IJK)
571                 Xi = XG_E(I)
572              ENDIF
573           ENDIF
574     
575     
576     
577     !======================================================================
578     !  Intersection with Edge 6 (node 6-8, Face East-Top):
579     !======================================================================
580           xa = X_NODE(6)
581           ya = Y_NODE(6)
582           za = Z_NODE(6)
583     
584           N_int_y = 0
585           INTERSECT_Y(IJK) = .FALSE.
586           Q_ID = 1
587              CALL INTERSECT_LINE('QUADRIC',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_FLAG,xc,yc,zc)
588              IF ((INTERSECT_FLAG).AND.(yc/=Yi)) THEN
589                 N_int_y = N_int_y + 1
590                 INTERSECT_Y(IJK) = .TRUE.
591                 yc_backup = Yi
592                 Yi = yc
593              ENDIF
594     
595           IF(N_int_y /= 1) THEN
596              Yi = UNDEFINED
597              INTERSECT_Y(IJK) = .FALSE.
598           ENDIF
599     
600           DO Q_ID = 1, N_POLYGON
601              CALL INTERSECT_LINE('POLYGON',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_Y(IJK),xc,yc,zc)
602              IF(INTERSECT_Y(IJK)) Yi = yc
603           ENDDO
604     
605           DO N_USR= 1, N_USR_DEF
606              CALL INTERSECT_LINE('USR_DEF',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_Y(IJK),xc,yc,zc)
607              IF(INTERSECT_Y(IJK)) Yi = yc
608           ENDDO
609     
610     !      IF(USE_STL) THEN
611     !         CALL INTERSECT_LINE_WITH_STL(xa,ya,za,xb,yb,zb,INTERSECT_Y(IJK),xc,yc,zc)
612     !         IF(INTERSECT_Y(IJK)) Yi = yc
613     !      ENDIF
614     
615           IF(TYPE_OF_CELL=='V_MOMENTUM') THEN
616              IF(SNAP(IJK)) THEN
617                 INTERSECT_Y(IJK) = .TRUE.
618                 J = J_OF(IJK)
619                 Yi = YG_N(J)
620              ENDIF
621           ENDIF
622     
623           IF(DO_K) THEN
624     !======================================================================
625     !  Intersection with Edge 11 (node 4-8, Face East-North):
626     !======================================================================
627              xa = X_NODE(4)
628              ya = Y_NODE(4)
629              za = Z_NODE(4)
630     
631              N_int_z = 0
632              INTERSECT_Z(IJK) = .FALSE.
633              Q_ID = 1
634                 CALL INTERSECT_LINE('QUADRIC',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_FLAG,xc,yc,zc)
635                 IF ((INTERSECT_FLAG).AND.(zc/=Zi)) THEN
636                    N_int_z = N_int_z + 1
637                    INTERSECT_Z(IJK) = .TRUE.
638                    zc_backup = Zi
639                    Zi = zc
640                 ENDIF
641     
642                 IF(N_int_z /= 1) THEN
643                    Zi = UNDEFINED
644                    INTERSECT_Z(IJK) = .FALSE.
645                 ENDIF
646     
647              DO Q_ID = 1, N_POLYGON
648                 CALL INTERSECT_LINE('POLYGON',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_Z(IJK),xc,yc,zc)
649                 IF(INTERSECT_Z(IJK)) Zi = zc
650              ENDDO
651     
652              DO N_USR= 1, N_USR_DEF
653                 CALL INTERSECT_LINE('USR_DEF',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_Z(IJK),xc,yc,zc)
654                 IF(INTERSECT_Z(IJK)) Zi = zc
655              ENDDO
656     
657        !      IF(USE_STL) THEN
658        !         CALL INTERSECT_LINE_WITH_STL(xa,ya,za,xb,yb,zb,INTERSECT_Z(IJK),xc,yc,zc)
659        !         IF(INTERSECT_Z(IJK)) Zi = zc
660        !      ENDIF
661     
662              IF(TYPE_OF_CELL=='W_MOMENTUM') THEN
663                 IF(SNAP(IJK)) THEN
664                    INTERSECT_Z(IJK) = .TRUE.
665                    K = K_OF(IJK)
666                    Zi = ZG_T(K)
667                 ENDIF
668              ENDIF
669     
670           ENDIF
671     
672     
673           IF(INTERSECT_X(IJK)) THEN
674              POTENTIAL_CUT_CELL_AT(IJK) = .TRUE.
675              POTENTIAL_CUT_CELL_AT(EAST_OF(IJK)) = .TRUE.
676           ENDIF
677     
678           IF(INTERSECT_Y(IJK)) THEN
679              POTENTIAL_CUT_CELL_AT(IJK) = .TRUE.
680              POTENTIAL_CUT_CELL_AT(NORTH_OF(IJK)) = .TRUE.
681           ENDIF
682     
683     
684           IF(INTERSECT_Z(IJK)) THEN
685              POTENTIAL_CUT_CELL_AT(IJK) = .TRUE.
686              POTENTIAL_CUT_CELL_AT(TOP_OF(IJK)) = .TRUE.
687           ENDIF
688     
689     
690           RETURN
691     
692           END SUBROUTINE INTERSECT
693     
694     
695     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
696     !                                                                      C
697     !  Module name: CLEAN_INTERSECT                                        C
698     !  Purpose: Remove Intersection flags in preparation of small cell     C
699     !           removal                                                    C
700     !                                                                      C
701     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
702     !  Reviewer:                                          Date:            C
703     !                                                                      C
704     !  Revision Number #                                  Date: ##-###-##  C
705     !  Author: #                                                           C
706     !  Purpose: #                                                          C
707     !                                                                      C
708     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
709       SUBROUTINE CLEAN_INTERSECT(IJK,TYPE_OF_CELL,Xi,Yi,Zi)
710     
711           USE param
712           USE param1
713           USE parallel
714           USE constant
715           USE run
716           USE toleranc
717           USE geometry
718           USE indices
719           USE compar
720           USE sendrecv
721           USE quadric
722           USE cutcell
723           USE polygon
724           USE STL
725           USE functions
726     
727           IMPLICIT NONE
728           CHARACTER (LEN=*) :: TYPE_OF_CELL
729           INTEGER :: IJK,I,J,K,IM,JM,KM,IP,JP,KP
730           INTEGER :: BCID
731           INTEGER :: IMJK,IPJK,IJMK,IJPK,IJKM,IJKP,IMJPK,IMJKP,IPJMK,IJMKP,IPJKM,IJPKM
732           DOUBLE PRECISION :: xa,ya,za,xb,yb,zb
733           DOUBLE PRECISION :: Xi,Yi,Zi
734           DOUBLE PRECISION :: DFC,DFC_MAX,F4,F6,F7,F8
735           LOGICAL :: CLIP_FLAG,CAD,F_TEST
736     
737     
738     ! When inputing geometry from CAD (STL or MSH file), the snapping procedure is
739     ! dependent on the value of F at the cell corners
740     ! For other gemoetry inputs (say quadrics), This is not needed, and the value
741     ! of F_TEST is set to .TRUE. here
742           CAD = USE_MSH.OR.USE_STL
743           F_TEST = .TRUE.
744     
745     !======================================================================
746     !  Get coordinates of eight nodes
747     !======================================================================
748     
749           CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
750     
751           I = I_OF(IJK)
752           J = J_OF(IJK)
753           K = K_OF(IJK)
754     
755           IM = I - 1
756           JM = J - 1
757           KM = K - 1
758     
759           IP = I + 1
760           JP = J + 1
761           KP = K + 1
762     
763           IMJK = FUNIJK(IM,J,K)
764           IPJK = FUNIJK(IP,J,K)
765           IJMK = FUNIJK(I,JM,K)
766           IJPK = FUNIJK(I,JP,K)
767           IJKM = FUNIJK(I,J,KM)
768           IJKP = FUNIJK(I,J,KP)
769     
770           IMJPK = FUNIJK(IM,JP,K)
771           IMJKP = FUNIJK(IM,J,KP)
772     
773           IPJMK = FUNIJK(IP,JM,K)
774           IJMKP = FUNIJK(I,JM,KP)
775     
776           IPJKM = FUNIJK(IP,J,KM)
777           IJPKM = FUNIJK(I,JP,KM)
778     
779     
780           IF(IMJK<1.OR.IMJK>DIMENSION_3) IMJK = IJK
781           IF(IPJK<1.OR.IPJK>DIMENSION_3) IPJK = IJK
782           IF(IJMK<1.OR.IJMK>DIMENSION_3) IJMK = IJK
783           IF(IJPK<1.OR.IJPK>DIMENSION_3) IJPK = IJK
784           IF(IJKM<1.OR.IJKM>DIMENSION_3) IJKM = IJK
785           IF(IJKP<1.OR.IJKP>DIMENSION_3) IJKP = IJK
786     
787           IF(IMJPK<1.OR.IMJPK>DIMENSION_3) IMJPK = IJK
788           IF(IMJKP<1.OR.IMJKP>DIMENSION_3) IMJKP = IJK
789     
790           IF(IPJMK<1.OR.IPJMK>DIMENSION_3) IPJMK = IJK
791           IF(IJMKP<1.OR.IJMKP>DIMENSION_3) IJMKP = IJK
792     
793           IF(IPJKM<1.OR.IPJKM>DIMENSION_3) IPJKM = IJK
794           IF(IJPKM<1.OR.IJPKM>DIMENSION_3) IJPKM = IJK
795     
796     !======================================================================
797     !  Clean Intersection with Edge 7 (node 7-8, Face North-Top):
798     !======================================================================
799     
800           xa = X_NODE(7)
801           ya = Y_NODE(7)
802           za = Z_NODE(7)
803     
804           xb = X_NODE(8)
805           yb = Y_NODE(8)
806           zb = Z_NODE(8)
807     
808           DFC_MAX = TOL_SNAP(1) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)  ! MAXIMUM DISTANCE FROM CORNER
809     
810           IF(INTERSECT_X(IJK)) THEN
811     
812              DFC = DABS(Xi-xa) ! DISTANCE FROM CORNER (NODE 7)
813     
814              IF(CAD) F_TEST = (F_AT(IMJK)/=ZERO)
815              IF(DFC < DFC_MAX.AND.F_TEST) THEN
816                 IF(PRINT_WARNINGS) THEN
817                    WRITE(*,*)'MERGING X-INTERSECTION ALONG EDGE 7 ONTO NODE 7'
818                    WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
819                 ENDIF
820     
821                 INTERSECT_X(IJK)  = .FALSE.
822                 IF(I>=IMIN1) THEN
823                    INTERSECT_X(IMJK) = .FALSE.
824                    INTERSECT_Y(IMJK)  = .FALSE.
825                    INTERSECT_Y(IMJPK) = .FALSE.
826                    IF(DO_K) INTERSECT_Z(IMJK)  = .FALSE.
827                    IF(DO_K) INTERSECT_Z(IMJKP) = .FALSE.
828     
829                    SNAP(IMJK) = .TRUE.
830                 ENDIF
831     
832              ENDIF
833     
834     
835              DFC = DABS(Xi-xb) ! DISTANCE FROM CORNER (NODE 8)
836     
837              IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
838              IF(DFC < DFC_MAX.AND.F_TEST) THEN
839                 IF(PRINT_WARNINGS) THEN
840                    WRITE(*,*)'MERGING X-INTERSECTION ALONG EDGE 7 ONTO NODE 8'
841                    WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
842                 ENDIF
843     
844     
845                 INTERSECT_X(IJK)  = .FALSE.
846                 INTERSECT_Y(IJK)  = .FALSE.
847                 IF(DO_K) INTERSECT_Z(IJK)  = .FALSE.
848                 SNAP(IJK) = .TRUE.
849     
850                 IF(I<=IMAX1) INTERSECT_X(IPJK) = .FALSE.
851                 IF(J<=JMAX1) INTERSECT_Y(IJPK) = .FALSE.
852                 IF(DO_K.AND.(K<=KMAX1)) INTERSECT_Z(IJKP) = .FALSE.
853     
854     
855              ENDIF
856     
857              IF(USE_STL.OR.USE_MSH) THEN
858                 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,7,F7,CLIP_FLAG,BCID)
859                 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
860                 IF(F7*F8>TOL_STL**2) INTERSECT_X(IJK)  = .FALSE.
861              ENDIF
862     
863           ENDIF
864     
865     
866     
867     
868     !======================================================================
869     !  Clean Intersection with Edge 6 (node 6-8, Face East-Top):
870     !======================================================================
871     
872           xa = X_NODE(6)
873           ya = Y_NODE(6)
874           za = Z_NODE(6)
875     
876           DFC_MAX = TOL_SNAP(2) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)  ! MAXIMUM DISTANCE FROM CORNER
877     
878     
879           IF(INTERSECT_Y(IJK)) THEN
880     
881              DFC = DABS(Yi-ya) ! DISTANCE FROM CORNER (NODE 6)
882     
883              IF(CAD) F_TEST = (F_AT(IJMK)/=ZERO)
884              IF(DFC < DFC_MAX.AND.F_TEST) THEN
885                 IF(PRINT_WARNINGS) THEN
886                    WRITE(*,*)'MERGING Y-INTERSECTION ALONG EDGE 6 ONTO NODE 6'
887                    WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
888                 ENDIF
889     
890                 INTERSECT_Y(IJK)  = .FALSE.
891                 IF(J>=JMIN1) THEN
892                    INTERSECT_X(IJMK)  = .FALSE.
893                    INTERSECT_X(IPJMK) = .FALSE.
894                    INTERSECT_Y(IJMK) = .FALSE.
895                    IF(DO_K) INTERSECT_Z(IJMK)  = .FALSE.
896                    IF(DO_K) INTERSECT_Z(IJMKP) = .FALSE.
897     
898                    SNAP(IJMK) = .TRUE.
899                 ENDIF
900     
901              ENDIF
902     
903     
904              DFC = DABS(Yi-yb) ! DISTANCE FROM CORNER (NODE 8)
905     
906              IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
907              IF(DFC < DFC_MAX.AND.F_TEST) THEN
908                 IF(PRINT_WARNINGS) THEN
909                    WRITE(*,*)'MERGING Y-INTERSECTION ALONG EDGE 6 ONTO NODE 8'
910                    WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
911                 ENDIF
912     
913                 INTERSECT_X(IJK)  = .FALSE.
914                 INTERSECT_Y(IJK)  = .FALSE.
915                 IF(DO_K) INTERSECT_Z(IJK)  = .FALSE.
916                 SNAP(IJK) = .TRUE.
917     
918                 IF(I<=IMAX1) INTERSECT_X(IPJK) = .FALSE.
919                 IF(J<=JMAX1) INTERSECT_Y(IJPK) = .FALSE.
920                 IF(DO_K.AND.(K<=KMAX1)) INTERSECT_Z(IJKP) = .FALSE.
921     
922     
923              ENDIF
924     
925     
926              IF(USE_STL.OR.USE_MSH) THEN
927                 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,6,F6,CLIP_FLAG,BCID)
928                 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
929                 IF(F6*F8>TOL_STL**2) INTERSECT_Y(IJK)  = .FALSE.
930              ENDIF
931     
932           ENDIF
933     
934     
935           IF(DO_K) THEN
936     !======================================================================
937     !  Intersection with Edge 11 (node 4-8, Face East-North):
938     !======================================================================
939     
940              xa = X_NODE(4)
941              ya = Y_NODE(4)
942              za = Z_NODE(4)
943     
944              DFC_MAX = TOL_SNAP(3) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)  ! MAXIMUM DISTANCE FROM CORNER
945     
946              IF(INTERSECT_Z(IJK)) THEN
947     
948                 DFC = DABS(Zi-Za) ! DISTANCE FROM CORNER (NODE 4)
949     
950                 IF(CAD) F_TEST = (F_AT(IJKM)/=ZERO)
951                 IF(DFC < DFC_MAX.AND.F_TEST) THEN
952                    IF(PRINT_WARNINGS) THEN
953                       WRITE(*,*)'MERGING Z-INTERSECTION ALONG EDGE 11 ONTO NODE 4'
954                       WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
955                    ENDIF
956     
957                    INTERSECT_Z(IJK)  = .FALSE.
958     
959                    IF(K>=KMIN1) THEN
960                       INTERSECT_X(IJKM)  = .FALSE.
961                       INTERSECT_X(IPJKM) = .FALSE.
962                       INTERSECT_Y(IJKM)  = .FALSE.
963                       INTERSECT_Y(IJPKM) = .FALSE.
964                       INTERSECT_Z(IJKM) = .FALSE.
965     
966                       SNAP(IJKM) = .TRUE.
967                    ENDIF
968     
969                 ENDIF
970     
971     
972                 DFC = DABS(Zi-Zb) ! DISTANCE FROM CORNER (NODE 8)
973     
974                 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
975                 IF(DFC < DFC_MAX.AND.F_TEST) THEN
976                    IF(PRINT_WARNINGS) THEN
977                       WRITE(*,*)'MERGING Z-INTERSECTION ALONG EDGE 11 ONTO NODE 8'
978                       WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
979                    ENDIF
980     
981                    INTERSECT_X(IJK)  = .FALSE.
982                    INTERSECT_Y(IJK)  = .FALSE.
983                    INTERSECT_Z(IJK)  = .FALSE.
984                    SNAP(IJK) = .TRUE.
985     !            F_AT(IJKM) = UNDEFINED
986     !            IF(F_AT(IJK)/=ZERO) SNAP(IJK) = .TRUE.
987     !               F_AT(IJK) = ZERO
988     
989                    IF(I<=IMAX1) INTERSECT_X(IPJK) = .FALSE.
990                    IF(J<=JMAX1) INTERSECT_Y(IJPK) = .FALSE.
991                    IF(K<=KMAX1) INTERSECT_Z(IJKP) = .FALSE.
992     
993     
994                 ENDIF
995     
996                 IF(USE_STL.OR.USE_MSH) THEN
997                    CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,4,F4,CLIP_FLAG,BCID)
998                    CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
999                    IF(F4*F8>TOL_STL**2) INTERSECT_Z(IJK)  = .FALSE.
1000                 ENDIF
1001     
1002              ENDIF
1003     
1004           ENDIF
1005     
1006     
1007           RETURN
1008     
1009           END SUBROUTINE CLEAN_INTERSECT
1010     
1011     
1012     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1013     !                                                                      C
1014     !  Module name: CLEAN_INTERSECT_SCALAR                                 C
1015     !  Purpose: Remove Intersection flags in preparation of small cell     C
1016     !           removal                                                    C
1017     !                                                                      C
1018     !  Author: Jeff Dietiker                              Date: 04-Dec-14  C
1019     !  Reviewer:                                          Date:            C
1020     !                                                                      C
1021     !  Revision Number #                                  Date: ##-###-##  C
1022     !  Author: #                                                           C
1023     !  Purpose: #                                                          C
1024     !                                                                      C
1025     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1026       SUBROUTINE CLEAN_INTERSECT_SCALAR
1027     
1028           USE param
1029           USE param1
1030           USE parallel
1031           USE constant
1032           USE run
1033           USE toleranc
1034           USE geometry
1035           USE indices
1036           USE compar
1037           USE sendrecv
1038           USE quadric
1039           USE cutcell
1040           USE polygon
1041           USE STL
1042           USE functions
1043     
1044           IMPLICIT NONE
1045           INTEGER :: IJK
1046     
1047     
1048     
1049           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
1050           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
1051           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
1052           call send_recv(Xn_int,2)
1053           call send_recv(Ye_int,2)
1054           call send_recv(Zt_int,2)
1055     
1056           DO IJK = IJKSTART3, IJKEND3
1057              CALL SET_SNAP_FLAG(IJK,'SCALAR',Xn_int(IJK),Ye_int(IJK),Zt_int(IJK))
1058           END DO
1059     
1060           call SEND_RECEIVE_1D_LOGICAL(SNAP,2)
1061           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
1062           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
1063           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
1064           call send_recv(Xn_int,2)
1065           call send_recv(Ye_int,2)
1066           call send_recv(Zt_int,2)
1067     
1068           DO IJK = IJKSTART3, IJKEND3
1069              CALL REMOVE_INTERSECT_FLAG(IJK)
1070           END DO
1071     
1072           call SEND_RECEIVE_1D_LOGICAL(SNAP,2)
1073           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
1074           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
1075           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
1076           call send_recv(Xn_int,2)
1077           call send_recv(Ye_int,2)
1078           call send_recv(Zt_int,2)
1079     
1080           RETURN
1081     
1082           END SUBROUTINE CLEAN_INTERSECT_SCALAR
1083     
1084     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1085     !                                                                      C
1086     !  Module name: SET_SNAP_FLAG                                          C
1087     !  Purpose: Set SNAP flag in preparation of intersection flag removal  C
1088     !           removal                                                    C
1089     !                                                                      C
1090     !  Author: Jeff Dietiker                              Date: 04-Dec-14  C
1091     !  Reviewer:                                          Date:            C
1092     !                                                                      C
1093     !  Revision Number #                                  Date: ##-###-##  C
1094     !  Author: #                                                           C
1095     !  Purpose: #                                                          C
1096     !                                                                      C
1097     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1098       SUBROUTINE SET_SNAP_FLAG(IJK,TYPE_OF_CELL,Xi,Yi,Zi)
1099     
1100           USE param
1101           USE param1
1102           USE parallel
1103           USE constant
1104           USE run
1105           USE toleranc
1106           USE geometry
1107           USE indices
1108           USE compar
1109           USE sendrecv
1110           USE quadric
1111           USE cutcell
1112           USE polygon
1113           USE STL
1114           USE functions
1115     
1116           IMPLICIT NONE
1117           CHARACTER (LEN=*) :: TYPE_OF_CELL
1118           INTEGER :: IJK,I,J,K,IM,JM,KM
1119           INTEGER :: BCID
1120           INTEGER :: IMJK,IJMK,IJKM
1121           DOUBLE PRECISION :: xa,ya,za,xb,yb,zb
1122           DOUBLE PRECISION :: Xi,Yi,Zi
1123           DOUBLE PRECISION :: DFC,DFC_MAX,F4,F6,F7,F8
1124           LOGICAL :: CLIP_FLAG,CAD,F_TEST
1125     
1126     
1127     ! When inputing geometry from CAD (STL or MSH file), the snapping procedure is
1128     ! dependent on the value of F at the cell corners
1129     ! For other gemoetry inputs (say quadrics), This is not needed, and the value
1130     ! of F_TEST is set to .TRUE. here
1131           CAD = USE_MSH.OR.USE_STL
1132           F_TEST = .TRUE.
1133     
1134     !======================================================================
1135     !  Get coordinates of eight nodes
1136     !======================================================================
1137     
1138           CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
1139     
1140           I = I_OF(IJK)
1141           J = J_OF(IJK)
1142           K = K_OF(IJK)
1143     
1144           IM = I - 1
1145           JM = J - 1
1146           KM = K - 1
1147     
1148           IMJK = FUNIJK(IM,J,K)
1149           IJMK = FUNIJK(I,JM,K)
1150           IJKM = FUNIJK(I,J,KM)
1151     
1152           IF(IMJK<1.OR.IMJK>DIMENSION_3) IMJK = IJK
1153           IF(IJMK<1.OR.IJMK>DIMENSION_3) IJMK = IJK
1154           IF(IJKM<1.OR.IJKM>DIMENSION_3) IJKM = IJK
1155     
1156     !======================================================================
1157     !  Clean Intersection with Edge 7 (node 7-8, Face North-Top):
1158     !======================================================================
1159     
1160           xa = X_NODE(7)
1161           ya = Y_NODE(7)
1162           za = Z_NODE(7)
1163     
1164           xb = X_NODE(8)
1165           yb = Y_NODE(8)
1166           zb = Z_NODE(8)
1167     
1168     
1169           DFC_MAX = TOL_SNAP(1) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)  ! MAXIMUM DISTANCE FROM CORNER
1170     
1171           IF(INTERSECT_X(IJK)) THEN
1172     
1173              DFC = DABS(Xi-xa) ! DISTANCE FROM CORNER (NODE 7)
1174     
1175              IF(CAD) F_TEST = (F_AT(IMJK)/=ZERO)
1176              IF(DFC < DFC_MAX.AND.F_TEST) THEN
1177                 IF(PRINT_WARNINGS) THEN
1178                    WRITE(*,*)'MERGING X-INTERSECTION ALONG EDGE 7 ONTO NODE 7'
1179                    WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1180                 ENDIF
1181     
1182                 IF(I>=IMIN1) SNAP(IMJK) = .TRUE.
1183     
1184              ENDIF
1185     
1186     
1187              DFC = DABS(Xi-xb) ! DISTANCE FROM CORNER (NODE 8)
1188     
1189              IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
1190              IF(DFC < DFC_MAX.AND.F_TEST) THEN
1191                 IF(PRINT_WARNINGS) THEN
1192                    WRITE(*,*)'MERGING X-INTERSECTION ALONG EDGE 7 ONTO NODE 8'
1193                    WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1194                 ENDIF
1195     
1196                 SNAP(IJK) = .TRUE.
1197     
1198              ENDIF
1199     
1200              IF(USE_STL.OR.USE_MSH) THEN
1201                 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,7,F7,CLIP_FLAG,BCID)
1202                 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
1203                 IF(F7*F8>TOL_STL**2) INTERSECT_X(IJK)  = .FALSE.
1204              ENDIF
1205     
1206           ENDIF
1207     
1208     
1209     !======================================================================
1210     !  Clean Intersection with Edge 6 (node 6-8, Face East-Top):
1211     !======================================================================
1212     
1213           xa = X_NODE(6)
1214           ya = Y_NODE(6)
1215           za = Z_NODE(6)
1216     
1217           DFC_MAX = TOL_SNAP(2) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)  ! MAXIMUM DISTANCE FROM CORNER
1218     
1219     
1220           IF(INTERSECT_Y(IJK)) THEN
1221     
1222              DFC = DABS(Yi-ya) ! DISTANCE FROM CORNER (NODE 6)
1223     
1224              IF(CAD) F_TEST = (F_AT(IJMK)/=ZERO)
1225              IF(DFC < DFC_MAX.AND.F_TEST) THEN
1226                 IF(PRINT_WARNINGS) THEN
1227                    WRITE(*,*)'MERGING Y-INTERSECTION ALONG EDGE 6 ONTO NODE 6'
1228                    WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1229                 ENDIF
1230     
1231                 IF(J>=JMIN1) SNAP(IJMK) = .TRUE.
1232     
1233              ENDIF
1234     
1235     
1236              DFC = DABS(Yi-yb) ! DISTANCE FROM CORNER (NODE 8)
1237     
1238              IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
1239              IF(DFC < DFC_MAX.AND.F_TEST) THEN
1240                 IF(PRINT_WARNINGS) THEN
1241                    WRITE(*,*)'MERGING Y-INTERSECTION ALONG EDGE 6 ONTO NODE 8'
1242                    WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1243                 ENDIF
1244     
1245                 SNAP(IJK) = .TRUE.
1246     
1247     
1248              ENDIF
1249     
1250     
1251              IF(USE_STL.OR.USE_MSH) THEN
1252                 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,6,F6,CLIP_FLAG,BCID)
1253                 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
1254                 IF(F6*F8>TOL_STL**2) INTERSECT_Y(IJK)  = .FALSE.
1255              ENDIF
1256     
1257           ENDIF
1258     
1259     
1260           IF(DO_K) THEN
1261     !======================================================================
1262     !  Intersection with Edge 11 (node 4-8, Face East-North):
1263     !======================================================================
1264     
1265              xa = X_NODE(4)
1266              ya = Y_NODE(4)
1267              za = Z_NODE(4)
1268     
1269              DFC_MAX = TOL_SNAP(3) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)  ! MAXIMUM DISTANCE FROM CORNER
1270     
1271              IF(INTERSECT_Z(IJK)) THEN
1272     
1273                 DFC = DABS(Zi-Za) ! DISTANCE FROM CORNER (NODE 4)
1274     
1275                 IF(CAD) F_TEST = (F_AT(IJKM)/=ZERO)
1276                 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1277                    IF(PRINT_WARNINGS) THEN
1278                       WRITE(*,*)'MERGING Z-INTERSECTION ALONG EDGE 11 ONTO NODE 4'
1279                       WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1280                    ENDIF
1281     
1282                    IF(K>=KMIN1) SNAP(IJKM) = .TRUE.
1283     
1284                 ENDIF
1285     
1286                 DFC = DABS(Zi-Zb) ! DISTANCE FROM CORNER (NODE 8)
1287     
1288                 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
1289                 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1290                    IF(PRINT_WARNINGS) THEN
1291                       WRITE(*,*)'MERGING Z-INTERSECTION ALONG EDGE 11 ONTO NODE 8'
1292                       WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1293                    ENDIF
1294     
1295                    SNAP(IJK) = .TRUE.
1296     
1297                 ENDIF
1298     
1299     
1300                 IF(USE_STL.OR.USE_MSH) THEN
1301                    CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,4,F4,CLIP_FLAG,BCID)
1302                    CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
1303                    IF(F4*F8>TOL_STL**2) INTERSECT_Z(IJK)  = .FALSE.
1304                 ENDIF
1305     
1306              ENDIF
1307     
1308           ENDIF
1309     
1310     
1311           RETURN
1312     
1313           END SUBROUTINE SET_SNAP_FLAG
1314     
1315           !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1316     !                                                                      C
1317     !  Module name: REMOVE_INTERSECT_FLAG                                  C
1318     !  Purpose: Remove Intersection flags                                  C
1319     !           removal                                                    C
1320     !                                                                      C
1321     !  Author: Jeff Dietiker                              Date: 04-Dec-14  C
1322     !  Reviewer:                                          Date:            C
1323     !                                                                      C
1324     !  Revision Number #                                  Date: ##-###-##  C
1325     !  Author: #                                                           C
1326     !  Purpose: #                                                          C
1327     !                                                                      C
1328     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1329       SUBROUTINE REMOVE_INTERSECT_FLAG(IJK)
1330     
1331           USE param
1332           USE param1
1333           USE parallel
1334           USE constant
1335           USE run
1336           USE toleranc
1337           USE geometry
1338           USE indices
1339           USE compar
1340           USE sendrecv
1341           USE quadric
1342           USE cutcell
1343           USE polygon
1344           USE STL
1345           USE functions
1346     
1347           IMPLICIT NONE
1348           INTEGER :: IJK,I,J,K,IP,JP,KP
1349           INTEGER :: IPJK,IJPK,IJKP
1350     
1351     
1352           I = I_OF(IJK)
1353           J = J_OF(IJK)
1354           K = K_OF(IJK)
1355     
1356           IP = I + 1
1357           JP = J + 1
1358           KP = K + 1
1359     
1360           IPJK = FUNIJK(IP,J,K)
1361           IJPK = FUNIJK(I,JP,K)
1362           IJKP = FUNIJK(I,J,KP)
1363     
1364           IF(IPJK<1.OR.IPJK>DIMENSION_3) IPJK = IJK
1365           IF(IJPK<1.OR.IJPK>DIMENSION_3) IJPK = IJK
1366           IF(IJKP<1.OR.IJKP>DIMENSION_3) IJKP = IJK
1367     
1368            IF(SNAP(IJK)) THEN
1369               INTERSECT_X(IJK) = .FALSE.
1370               INTERSECT_Y(IJK) = .FALSE.
1371               INTERSECT_Z(IJK) = .FALSE.
1372               IF(I<=IMAX1) INTERSECT_X(IPJK) = .FALSE.
1373               IF(J<=JMAX1) INTERSECT_Y(IJPK) = .FALSE.
1374               IF(DO_K.AND.(K<=KMAX1)) INTERSECT_Z(IJKP) = .FALSE.
1375            ENDIF
1376     
1377     
1378           RETURN
1379     
1380           END SUBROUTINE REMOVE_INTERSECT_FLAG
1381     
1382     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1383     !                                                                      C
1384     !  Module name: OPEN_CUT_CELL_FILES                                    C
1385     !  Purpose: Open CUT CELL related file                                 C
1386     !           and writes headers                                         C
1387     !                                                                      C
1388     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1389     !  Reviewer:                                          Date:            C
1390     !                                                                      C
1391     !  Revision Number #                                  Date: ##-###-##  C
1392     !  Author: #                                                           C
1393     !  Purpose: #                                                          C
1394     !                                                                      C
1395     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1396       SUBROUTINE OPEN_CUT_CELL_FILES
1397     
1398           USE cutcell
1399           USE compar
1400     
1401           IMPLICIT NONE
1402     
1403           IF(MyPE == PE_IO)  THEN
1404              OPEN(UNIT = UNIT_CUT_CELL_LOG, FILE = 'CUT_CELL.LOG')
1405           ENDIF
1406     
1407           RETURN
1408     
1409     
1410           END SUBROUTINE OPEN_CUT_CELL_FILES
1411     
1412     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1413     !                                                                      C
1414     !  Module name: CLOSE_CUT_CELL_FILES                                   C
1415     !  Purpose: Close CUT CELL related file                                C
1416     !                                                                      C
1417     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1418     !  Reviewer:                                          Date:            C
1419     !                                                                      C
1420     !  Revision Number #                                  Date: ##-###-##  C
1421     !  Author: #                                                           C
1422     !  Purpose: #                                                          C
1423     !                                                                      C
1424     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1425       SUBROUTINE CLOSE_CUT_CELL_FILES
1426     
1427           USE compar
1428           USE cutcell
1429     
1430           IMPLICIT NONE
1431     
1432           IF(MyPE == PE_IO) THEN
1433              CLOSE(UNIT_CUT_CELL_LOG)
1434           ENDIF
1435     
1436           RETURN
1437     
1438     
1439           END SUBROUTINE CLOSE_CUT_CELL_FILES
1440     
1441     
1442     
1443     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1444     !                                                                      C
1445     !  Module name: CAD_INTERSECT                                          C
1446     !  Purpose: Intersects CAD (STL file or MSH file) geometry with grid   C
1447     !                                                                      C
1448     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1449     !  Reviewer:                                          Date:            C
1450     !                                                                      C
1451     !  Revision Number #                                  Date: ##-###-##  C
1452     !  Author: #                                                           C
1453     !  Purpose: #                                                          C
1454     !                                                                      C
1455     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1456       SUBROUTINE CAD_INTERSECT(TYPE_OF_CELL,Xint,Yint,Zint)
1457     
1458           USE param
1459           USE param1
1460           USE parallel
1461           USE constant
1462           USE run
1463           USE toleranc
1464           USE geometry
1465           USE indices
1466           USE compar
1467           USE sendrecv
1468           USE quadric
1469           USE cutcell
1470           USE polygon
1471           USE stl
1472     
1473           USE mpi_utility
1474           USE functions
1475     
1476           IMPLICIT NONE
1477           CHARACTER (LEN=*) :: TYPE_OF_CELL
1478           INTEGER :: IJK,I,J,K
1479           INTEGER :: IM,IP,JM,JP,KM,KP,IMJK,IPJK,IJMK,IJPK,IJKM,IJKP
1480           INTEGER :: IJPKP,IPJKP,IPJPK
1481           DOUBLE PRECISION :: xa,ya,za,xb,yb,zb,xc,yc,zc
1482           LOGICAL :: INTERSECT_FLAG,INSIDE_FACET_a,INSIDE_FACET_b
1483     
1484           DOUBLE PRECISION :: X1,X2,Y1,Y2,Z1,Z2
1485     
1486           DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: Xint,Yint,Zint
1487     
1488           INTEGER :: N,I1,I2,J1,J2,K1,K2
1489     
1490           DOUBLE PRECISION :: X_OFFSET, Y_OFFSET, Z_OFFSET
1491     
1492           DOUBLE PRECISION, DIMENSION(3) :: N4,N6,N7,N8
1493           DOUBLE PRECISION :: CURRENT_F
1494     
1495           INTEGER :: N_UNDEFINED, NTOTAL_UNDEFINED,N_PROP
1496           INTEGER, PARAMETER :: N_PROPMAX=1000
1497           LOGICAL:: F_FOUND
1498     
1499     !      CHARACTER (LEN=3) :: CAD_PROPAGATE_ORDER
1500     
1501           INTERSECT_X = .FALSE.
1502           INTERSECT_Y = .FALSE.
1503           INTERSECT_Z = .FALSE.
1504     
1505           Xint = UNDEFINED
1506           Yint = UNDEFINED
1507           Zint = UNDEFINED
1508     
1509           F_AT = UNDEFINED
1510     
1511           SELECT CASE (TYPE_OF_CELL)
1512              CASE('SCALAR')
1513     
1514     
1515                 X_OFFSET = ZERO
1516                 Y_OFFSET = ZERO
1517                 Z_OFFSET = ZERO
1518     
1519              CASE('U_MOMENTUM')
1520     
1521                 X_OFFSET = HALF
1522                 Y_OFFSET = ZERO
1523                 Z_OFFSET = ZERO
1524     
1525              CASE('V_MOMENTUM')
1526     
1527                 X_OFFSET = ZERO
1528                 Y_OFFSET = HALF
1529                 Z_OFFSET = ZERO
1530     
1531              CASE('W_MOMENTUM')
1532     
1533                 X_OFFSET = ZERO
1534                 Y_OFFSET = ZERO
1535                 Z_OFFSET = HALF
1536     
1537     
1538              CASE DEFAULT
1539                 WRITE(*,*)'SUBROUTINE: GET_CELL_NODE_COORDINATES'
1540                 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
1541                 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
1542                 WRITE(*,*)'SCALAR'
1543                 WRITE(*,*)'U_MOMENTUM'
1544                 WRITE(*,*)'V_MOMENTUM'
1545                 WRITE(*,*)'W_MOMENTUM'
1546                 CALL MFIX_EXIT(myPE)
1547           END SELECT
1548     
1549     
1550           DO N = 1,N_FACETS
1551     
1552     
1553              X1 = MINVAL(VERTEX(1:3,1,N))
1554              X2 = MAXVAL(VERTEX(1:3,1,N))
1555              Y1 = MINVAL(VERTEX(1:3,2,N))
1556              Y2 = MAXVAL(VERTEX(1:3,2,N))
1557              Z1 = MINVAL(VERTEX(1:3,3,N))
1558              Z2 = MAXVAL(VERTEX(1:3,3,N))
1559     
1560     
1561              I1 = IEND3
1562              I2 = ISTART3
1563     
1564              IF(X2>=ZERO-DX(ISTART3).AND.X1<=XLENGTH+DX(IEND3)) THEN
1565                 DO I = ISTART3, IEND3
1566                    IP = I+1
1567                    IF(XG_E(I)+X_OFFSET*DX(IP)>=X1-TOL_STL) THEN
1568                       I1=I
1569                       EXIT
1570                    ENDIF
1571                 ENDDO
1572     
1573                 DO I = IEND3, ISTART3,-1
1574                    IP = I+1
1575                    IF(XG_E(I)-DX(I)+X_OFFSET*DX(IP)<=X2+TOL_STL) THEN
1576                       I2=I
1577                       EXIT
1578                    ENDIF
1579                 ENDDO
1580              ENDIF
1581     
1582     
1583              J1 = JEND3
1584              J2 = JSTART3
1585     
1586              IF(Y2>=ZERO-DY(JSTART3).AND.Y1<=YLENGTH+DY(JEND3)) THEN
1587                 DO J = JSTART3, JEND3
1588                    JP = J+1
1589                    IF(YG_N(J)+Y_OFFSET*DY(JP)>=Y1-TOL_STL) THEN
1590                       J1=J
1591                       EXIT
1592                    ENDIF
1593                 ENDDO
1594     
1595                 DO J = JEND3, JSTART3,-1
1596                    JP=J+1
1597                    IF(YG_N(J)-DY(J)+Y_OFFSET*DY(JP)<=Y2+TOL_STL) THEN
1598                       J2=J
1599                       EXIT
1600                    ENDIF
1601                 ENDDO
1602              ENDIF
1603     
1604              K1 = KEND3
1605              K2 = KSTART3
1606     
1607              IF(Z2>=ZERO-DZ(KSTART3).AND.Z1<=ZLENGTH+DZ(KEND3)) THEN
1608                 DO K = KSTART3, KEND3
1609                    KP=K+1
1610     
1611                    IF(ZG_T(K)+Z_OFFSET*DZ(KP)>=Z1-TOL_STL) THEN
1612                       K1=K
1613                       EXIT
1614                    ENDIF
1615                 ENDDO
1616     
1617                 DO K = KEND3, KSTART3,-1
1618                    KP = K+1
1619                    IF(ZG_T(K)-DZ(K)+Z_OFFSET*DZ(KP)<=Z2+TOL_STL) THEN
1620                       K2=K
1621                       EXIT
1622                    ENDIF
1623                 ENDDO
1624              ENDIF
1625     
1626     
1627     
1628     
1629              DO K=K1,K2
1630                 DO J=J1,J2
1631                    DO I=I1,I2
1632     
1633                       IJK = FUNIJK(I,J,K)
1634     
1635     
1636                       IM = MAX0(I - 1 ,ISTART3)
1637                       JM = MAX0(J - 1 ,JSTART3)
1638                       KM = MAX0(K - 1 ,KSTART3)
1639     
1640                       IP = MIN0(I + 1 ,IEND3)
1641                       JP = MIN0(J + 1 ,JEND3)
1642                       KP = MIN0(K + 1 ,KEND3)
1643     
1644     
1645                       IMJK = FUNIJK(IM,J,K)
1646                       IPJK = FUNIJK(IP,J,K)
1647                       IJMK = FUNIJK(I,JM,K)
1648                       IJPK = FUNIJK(I,JP,K)
1649                       IJKM = FUNIJK(I,J,KM)
1650                       IJKP = FUNIJK(I,J,KP)
1651     
1652                       IJPKP = FUNIJK(I,JP,KP)
1653                       IPJKP = FUNIJK(IP,J,KP)
1654                       IPJPK = FUNIJK(IP,JP,K)
1655     
1656     
1657     !======================================================================
1658     !  Get coordinates of eight nodes
1659     !======================================================================
1660     
1661                       CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
1662     
1663     !======================================================================
1664     !  Intersection with Edge 7 (node 7-8, Face North-Top):
1665     !======================================================================
1666                       xa = X_NODE(7)
1667                       ya = Y_NODE(7)
1668                       za = Z_NODE(7)
1669     
1670                       xb = X_NODE(8)
1671                       yb = Y_NODE(8)
1672                       zb = Z_NODE(8)
1673     
1674     ! Check if intersection occurs at corners
1675     
1676                       CALL IS_POINT_INSIDE_FACET(xa,ya,za,N,INSIDE_FACET_a)
1677     
1678                       IF(INSIDE_FACET_a) THEN   ! corner intersection at node 7
1679     
1680                          F_AT(IMJK) = ZERO
1681     
1682                          IF(TRIM(TYPE_OF_CELL).eq.'SCALAR')  CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1683     
1684                       ENDIF
1685     
1686     
1687                       CALL IS_POINT_INSIDE_FACET(xb,yb,zb,N,INSIDE_FACET_b)
1688     
1689                       IF(INSIDE_FACET_b) THEN   ! corner intersection at node 8
1690     
1691                          F_AT(IJK) = ZERO
1692     
1693                          IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1694     
1695                       ENDIF
1696     
1697     
1698     
1699     ! Check intersection within line 7-8, excluding corners
1700     
1701                       INTERSECT_FLAG = .FALSE.
1702     
1703                       IF(.NOT.(INTERSECT_X(IJK).OR.INSIDE_FACET_a.OR.INSIDE_FACET_b)) THEN
1704                          CALL INTERSECT_LINE_WITH_FACET(xa,ya,za,xb,yb,zb,N,INTERSECT_FLAG,xc,yc,zc)
1705                       ENDIF
1706     
1707                       IF(INTERSECT_FLAG) THEN
1708                          IF(INTERSECT_X(IJK)) THEN
1709                             IF(DABS(Xint(IJK)-xc)>TOL_STL) THEN
1710     
1711                                ! Ignore intersections when two intersections are detected on the same edge
1712                                INTERSECT_X(IJK) = .FALSE.
1713                             ENDIF
1714                          ELSE
1715                             INTERSECT_X(IJK) = .TRUE.
1716                             Xint(IJK) = xc
1717     
1718     ! Set values at corners if they are not zero
1719     
1720                             N7(1) = xa-xc
1721                             N7(2) = ya-yc
1722                             N7(3) = za-zc
1723     
1724                             IF(DABS(F_AT(IMJK))>TOL_F)   F_AT(IMJK) = -DOT_PRODUCT(N7,NORM_FACE(:,N))
1725     
1726                             N8(1) = xb-xc
1727                             N8(2) = yb-yc
1728                             N8(3) = zb-zc
1729     
1730                             IF(DABS(F_AT(IJK))>TOL_F)   F_AT(IJK) = -DOT_PRODUCT(N8,NORM_FACE(:,N))
1731     
1732     
1733                             IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1734                             IF(JP<=J2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJPK,N)
1735                             IF(KP<=K2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJKP,N)
1736                             IF(JP<=J2.AND.KP<=K2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJPKP,N)
1737                          ENDIF
1738                       ENDIF
1739     
1740                       IF(TYPE_OF_CELL=='U_MOMENTUM') THEN
1741                          IF(SNAP(IJK)) THEN
1742                             INTERSECT_X(IJK) = .TRUE.
1743                             Xn_int(IJK) = XG_E(I)
1744                          ENDIF
1745                       ENDIF
1746     
1747     !======================================================================
1748     !  Intersection with Edge 6 (node 6-8, Face East-Top):
1749     !======================================================================
1750                       xa = X_NODE(6)
1751                       ya = Y_NODE(6)
1752                       za = Z_NODE(6)
1753     
1754     ! Check if intersection occurs at corners
1755     
1756                       CALL IS_POINT_INSIDE_FACET(xa,ya,za,N,INSIDE_FACET_a)
1757     
1758                       IF(INSIDE_FACET_a) THEN   ! corner intersection at node 6
1759     
1760                          F_AT(IJMK) = ZERO
1761     
1762                          IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1763     
1764                       ENDIF
1765     
1766     
1767     
1768     ! Check intersection within line 6-8, excluding corners
1769     
1770     
1771     
1772                       INTERSECT_FLAG = .FALSE.
1773     
1774                       IF(.NOT.(INTERSECT_Y(IJK).OR.INSIDE_FACET_a.OR.INSIDE_FACET_b)) THEN
1775                          CALL INTERSECT_LINE_WITH_FACET(xa,ya,za,xb,yb,zb,N,INTERSECT_FLAG,xc,yc,zc)
1776                       ENDIF
1777     
1778     
1779                       IF(INTERSECT_FLAG) THEN
1780     
1781                          IF(INTERSECT_Y(IJK)) THEN
1782     
1783                             IF(DABS(Yint(IJK)-yc)>TOL_STL) THEN
1784     
1785                                INTERSECT_Y(IJK) = .FALSE. ! Ignore intersections when two intersections are detected on the same edge
1786     
1787                             ENDIF
1788     
1789                          ELSE
1790     
1791     
1792                             INTERSECT_Y(IJK) = .TRUE.
1793                             Yint(IJK) = yc
1794     
1795     ! Set values at corners if they are not zero
1796     
1797                             N6(1) = xa-xc
1798                             N6(2) = ya-yc
1799                             N6(3) = za-zc
1800     
1801                             IF(DABS(F_AT(IJMK))>TOL_F)   F_AT(IJMK) = -DOT_PRODUCT(N6,NORM_FACE(:,N))
1802     
1803                             N8(1) = xb-xc
1804                             N8(2) = yb-yc
1805                             N8(3) = zb-zc
1806     
1807                             IF(DABS(F_AT(IJK))>TOL_F)   F_AT(IJK) = -DOT_PRODUCT(N8,NORM_FACE(:,N))
1808     
1809     
1810                             IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1811                             IF(IP<=I2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IPJK,N)
1812                             IF(KP<=K2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJKP,N)
1813                             IF(IP<=I2.AND.KP<=K2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IPJKP,N)
1814     
1815                          ENDIF
1816     
1817                       ENDIF
1818     
1819                       IF(TYPE_OF_CELL=='V_MOMENTUM') THEN
1820                          IF(SNAP(IJK)) THEN
1821                             INTERSECT_Y(IJK) = .TRUE.
1822                             Ye_int(IJK) = YG_N(J)
1823                          ENDIF
1824                       ENDIF
1825     
1826                       IF(DO_K) THEN
1827     !======================================================================
1828     !  Intersection with Edge 11 (node 4-8, Face East-North):
1829     !======================================================================
1830                          xa = X_NODE(4)
1831                          ya = Y_NODE(4)
1832                          za = Z_NODE(4)
1833     
1834     ! Check if intersection occurs at corners
1835     
1836                          CALL IS_POINT_INSIDE_FACET(xa,ya,za,N,INSIDE_FACET_a)
1837     
1838                          IF(INSIDE_FACET_a) THEN   ! corner intersection at node 4
1839     
1840                             F_AT(IJKM) = ZERO
1841     
1842                             IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1843     
1844                          ENDIF
1845     
1846     
1847     
1848     
1849     ! Check intersection within line 4-8, excluding corners
1850     
1851     
1852                          INTERSECT_FLAG = .FALSE.
1853     
1854                          IF(.NOT.(INTERSECT_Z(IJK).OR.INSIDE_FACET_a.OR.INSIDE_FACET_b)) THEN
1855                             CALL INTERSECT_LINE_WITH_FACET(xa,ya,za,xb,yb,zb,N,INTERSECT_FLAG,xc,yc,zc)
1856                          ENDIF
1857     
1858                          IF(INTERSECT_FLAG) THEN
1859     
1860                             IF(INTERSECT_Z(IJK)) THEN
1861     
1862                                IF(DABS(Zint(IJK)-zc)>TOL_STL) THEN
1863     
1864                                   INTERSECT_Z(IJK) = .FALSE. ! Ignore intersections when two intersections are detected on the same edge
1865     
1866                                ENDIF
1867     
1868                             ELSE
1869     
1870     
1871                                INTERSECT_Z(IJK) = .TRUE.
1872                                Zint(IJK) = zc
1873     
1874     
1875     ! Set values at corners if they are not zero
1876     
1877                                N4(1) = xa-xc
1878                                N4(2) = ya-yc
1879                                N4(3) = za-zc
1880     
1881                                IF(DABS(F_AT(IJKM))>TOL_F)   F_AT(IJKM) = -DOT_PRODUCT(N4,NORM_FACE(:,N))
1882     
1883                                N8(1) = xb-xc
1884                                N8(2) = yb-yc
1885                                N8(3) = zb-zc
1886     
1887                                IF(DABS(F_AT(IJK))>TOL_F)   F_AT(IJK) = -DOT_PRODUCT(N8,NORM_FACE(:,N))
1888     
1889     
1890     
1891                                IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1892                                IF(IP<=I2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IPJK,N)
1893                                IF(JP<=J2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJPK,N)
1894                                IF(IP<=I2.AND.JP<=J2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IPJPK,N)
1895     
1896                             ENDIF
1897     
1898                          ENDIF
1899     
1900     
1901                          IF(TYPE_OF_CELL=='W_MOMENTUM') THEN
1902                             IF(SNAP(IJK)) THEN
1903                                INTERSECT_Z(IJK) = .TRUE.
1904                                Zt_int(IJK) = ZG_T(K)
1905                             ENDIF
1906                          ENDIF
1907     
1908                       ENDIF
1909     
1910     
1911                    ENDDO  ! I loop
1912                 ENDDO  ! J loop
1913              ENDDO  ! K loop
1914     
1915     
1916           ENDDO  ! Loop over facets
1917     
1918           CURRENT_F = UNDEFINED
1919     
1920     ! Overwrite small values to set them to zero
1921     
1922           DO IJK = IJKSTART3, IJKEND3
1923              IF(DABS(F_AT(IJK))<TOL_STL) THEN
1924                 F_AT(IJK)=ZERO
1925              ENDIF
1926           END DO
1927     
1928     ! Propagates node values to all interior cells
1929     ! in order defined by CAD_PROPAGATE_ORDER (outer loop)
1930     
1931     
1932           call send_recv(F_AT,2)
1933           call send_recv(Xint,2)
1934           call send_recv(Yint,2)
1935           call send_recv(Zint,2)
1936           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
1937           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
1938           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
1939     !======================================================================
1940     !  Clean-up intersection flags in preparaton of small cells removal
1941     !======================================================================
1942     
1943           DO IJK = IJKSTART3, IJKEND3
1944     
1945     !        IF(INTERIOR_CELL_AT(IJK)) THEN
1946     
1947     !            IF(POTENTIAL_CUT_CELL_AT(IJK))  CALL CLEAN_INTERSECT(IJK,'SCALAR',Xn_int(IJK),Ye_int(IJK),Zt_int(IJK))
1948     
1949                 CALL CLEAN_INTERSECT(IJK,TYPE_OF_CELL,Xint(IJK),Yint(IJK),Zint(IJK))
1950     
1951     !        ENDIF
1952     
1953           END DO
1954     
1955           call send_recv(F_AT,2)
1956           call SEND_RECEIVE_1D_LOGICAL(SNAP,2)
1957           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
1958           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
1959           call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
1960     
1961           SELECT CASE (CAD_PROPAGATE_ORDER)
1962     
1963           CASE ('   ')
1964     
1965     ! After intersecting the edges of the background mesh with the STL facets,
1966     ! the end points (i.e., cell corners) are assigned a value, called F_AT, where:
1967     ! F_AT = zero if the corner point is on a facet (within some tolerance TOL_STL),
1968     ! F_AT < zero if the corner point is inside  the fluid region,
1969     ! F_AT > zero if the corner point is outside the fluid region.
1970     ! At this point F_AT is only defined across edges that intersect the STL facets,
1971     ! and it must be propagated to all cell corners to determine if uncut cells
1972     ! are inside or outside the fluid region.
1973     
1974     ! Only F_AT values that are defined and not zero are propagated to their direct
1975     ! neighbors, if it is not already defined. The propagation is repeated
1976     ! at most N_PROPMAX. The loop is exited when all F_AT values are defined.
1977     ! N_PROPMAX could be increased for very large domains.
1978     ! The propagation of F_AT will stop anytime a boundary is encountered since F_AT
1979     ! changes sign across a boundary.
1980     !
1981              DO N_PROP=1,N_PROPMAX
1982     
1983                 DO IJK = IJKSTART3, IJKEND3
1984     
1985     ! Aaron
1986                    I = I_OF(IJK)
1987                    J = J_OF(IJK)
1988                    K = K_OF(IJK)
1989                    IF(.NOT.IS_ON_myPE_plus1layer(I,J,K))cycle
1990     ! End aaron
1991     
1992                    IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
1993     
1994                          IMJK = IM_OF(IJK)
1995                          IF(F_AT(IMJK)==UNDEFINED.AND.F_AT(IMJK)/=ZERO)  F_AT(IMJK)=F_AT(IJK)
1996     
1997                          IPJK = IP_OF(IJK)
1998                          IF(F_AT(IPJK)==UNDEFINED.AND.F_AT(IPJK)/=ZERO)  F_AT(IPJK)=F_AT(IJK)
1999     
2000                          IJMK = JM_OF(IJK)
2001                          IF(F_AT(IJMK)==UNDEFINED.AND.F_AT(IJMK)/=ZERO)  F_AT(IJMK)=F_AT(IJK)
2002     
2003                          IJPK = JP_OF(IJK)
2004                          IF(F_AT(IJPK)==UNDEFINED.AND.F_AT(IJPK)/=ZERO)  F_AT(IJPK)=F_AT(IJK)
2005     
2006                          IJKM = KM_OF(IJK)
2007                          IF(F_AT(IJKM)==UNDEFINED.AND.F_AT(IJKM)/=ZERO)  F_AT(IJKM)=F_AT(IJK)
2008     
2009                          IJKP = KP_OF(IJK)
2010                          IF(F_AT(IJKP)==UNDEFINED.AND.F_AT(IJKP)/=ZERO)  F_AT(IJKP)=F_AT(IJK)
2011     
2012                    ENDIF
2013     
2014                 ENDDO ! IJK Loop
2015     
2016     
2017     ! Communicate F_AT accross processors for DMP runs
2018                 call send_recv(F_AT,2)
2019     
2020     ! Count the number of undefined values of F_AT
2021     ! and exit loop if all values of F_AT have been propagated
2022                 N_UNDEFINED = 0
2023                 DO IJK = IJKSTART3, IJKEND3
2024                    IF(INTERIOR_CELL_AT(IJK).AND.F_AT(IJK)==UNDEFINED) N_UNDEFINED = N_UNDEFINED + 1
2025                 ENDDO
2026     
2027                 call global_all_sum( N_UNDEFINED, NTOTAL_UNDEFINED )
2028                 IF(NTOTAL_UNDEFINED==0) EXIT
2029     
2030              ENDDO ! N_PROP Loop
2031     
2032     
2033              call send_recv(F_AT,2)
2034     
2035     ! If a process still has undefined values of F_AT, this probably means
2036     ! that all cells belonging to that process are dead cells.
2037              IF(N_UNDEFINED>0) THEN
2038                 WRITE(*,*)'WARNING: UNABLE TO PROPAGATE F_AT ARRAY FROM myPE=.', MyPE
2039                 WRITE(*,*)'         THIS USUALLY INDICATE A RANK WITH NO ACTIVE CELL'
2040                 WRITE(*,*)'         YOU MAY NEED TO ADJUST THE GRID PARTITIONNING'
2041                 WRITE(*,*)'         TO GET BETTER LOAD_BALANCE.'
2042              ENDIF
2043     
2044     
2045     
2046           CASE ('IJK')
2047     
2048              DO I=ISTART3,IEND3
2049                 DO K=KSTART3,KEND3
2050                    F_FOUND = .FALSE.
2051                    DO J=JSTART3,JEND3
2052                       IJK=FUNIJK(I,J,K)
2053                       IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2054                          F_FOUND = .TRUE.
2055                          CURRENT_F = F_AT(IJK)
2056                          EXIT
2057                       ENDIF
2058                    ENDDO
2059                    IF(F_FOUND) THEN
2060                       DO J=JSTART3,JEND3
2061                          IJK=FUNIJK(I,J,K)
2062                          IF(F_AT(IJK)==UNDEFINED) THEN
2063                             F_AT(IJK)=CURRENT_F
2064                          ELSEIF(F_AT(IJK)/=ZERO) THEN
2065                             CURRENT_F = F_AT(IJK)
2066                          ENDIF
2067                       ENDDO
2068                    ENDIF
2069                 ENDDO
2070              ENDDO
2071     
2072              call send_recv(F_AT,2)
2073     
2074              DO J=JSTART3,JEND3
2075                 DO I=ISTART3,IEND3
2076                    F_FOUND = .FALSE.
2077                    DO K=KSTART3,KEND3
2078                       IJK=FUNIJK(I,J,K)
2079                       IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2080                          F_FOUND = .TRUE.
2081                          CURRENT_F = F_AT(IJK)
2082                          EXIT
2083                       ENDIF
2084                    ENDDO
2085                    IF(F_FOUND) THEN
2086                       DO K=KSTART3,KEND3
2087                          IJK=FUNIJK(I,J,K)
2088                          IF(F_AT(IJK)==UNDEFINED) THEN
2089                             F_AT(IJK)=CURRENT_F
2090                          ELSEIF(F_AT(IJK)/=ZERO) THEN
2091                             CURRENT_F = F_AT(IJK)
2092                          ENDIF
2093                       ENDDO
2094                    ENDIF
2095                 ENDDO
2096              ENDDO
2097     
2098              call send_recv(F_AT,2)
2099     
2100              DO K=KSTART3,KEND3
2101                 DO J=JSTART3,JEND3
2102                    F_FOUND = .FALSE.
2103                    DO I=ISTART3,IEND3
2104                       IJK=FUNIJK(I,J,K)
2105                       IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2106                          F_FOUND = .TRUE.
2107                          CURRENT_F = F_AT(IJK)
2108                          EXIT
2109                       ENDIF
2110                    ENDDO
2111                    IF(F_FOUND) THEN
2112                       DO I=ISTART3,IEND3
2113                          IJK=FUNIJK(I,J,K)
2114                          IF(F_AT(IJK)==UNDEFINED) THEN
2115                             F_AT(IJK)=CURRENT_F
2116                          ELSEIF(F_AT(IJK)/=ZERO) THEN
2117                             CURRENT_F = F_AT(IJK)
2118                          ENDIF
2119                       ENDDO
2120                    ENDIF
2121                 ENDDO
2122              ENDDO
2123     
2124              call send_recv(F_AT,2)
2125     
2126           CASE ('JKI')
2127     
2128              DO J=JSTART3,JEND3
2129                 DO I=ISTART3,IEND3
2130                    F_FOUND = .FALSE.
2131                    DO K=KSTART3,KEND3
2132                       IJK=FUNIJK(I,J,K)
2133                       IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2134                          F_FOUND = .TRUE.
2135                          CURRENT_F = F_AT(IJK)
2136                          EXIT
2137                       ENDIF
2138                    ENDDO
2139                    IF(F_FOUND) THEN
2140                       DO K=KSTART3,KEND3
2141                          IJK=FUNIJK(I,J,K)
2142                          IF(F_AT(IJK)==UNDEFINED) THEN
2143                             F_AT(IJK)=CURRENT_F
2144                          ELSEIF(F_AT(IJK)/=ZERO) THEN
2145                             CURRENT_F = F_AT(IJK)
2146                          ENDIF
2147                       ENDDO
2148                    ENDIF
2149                 ENDDO
2150              ENDDO
2151     
2152              call send_recv(F_AT,2)
2153     
2154              DO K=KSTART3,KEND3
2155                 DO J=JSTART3,JEND3
2156                    F_FOUND = .FALSE.
2157                    DO I=ISTART3,IEND3
2158                       IJK=FUNIJK(I,J,K)
2159                       IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2160                          F_FOUND = .TRUE.
2161                          CURRENT_F = F_AT(IJK)
2162                          EXIT
2163                       ENDIF
2164                    ENDDO
2165                    IF(F_FOUND) THEN
2166                       DO I=ISTART3,IEND3
2167                          IJK=FUNIJK(I,J,K)
2168                          IF(F_AT(IJK)==UNDEFINED) THEN
2169                             F_AT(IJK)=CURRENT_F
2170                          ELSEIF(F_AT(IJK)/=ZERO) THEN
2171                             CURRENT_F = F_AT(IJK)
2172                          ENDIF
2173                       ENDDO
2174                    ENDIF
2175                 ENDDO
2176              ENDDO
2177     
2178              call send_recv(F_AT,2)
2179     
2180              DO I=ISTART3,IEND3
2181                 DO K=KSTART3,KEND3
2182                    F_FOUND = .FALSE.
2183                    DO J=JSTART3,JEND3
2184                       IJK=FUNIJK(I,J,K)
2185                       IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2186                          F_FOUND = .TRUE.
2187                          CURRENT_F = F_AT(IJK)
2188                          EXIT
2189                       ENDIF
2190                    ENDDO
2191                    IF(F_FOUND) THEN
2192                       DO J=JSTART3,JEND3
2193                          IJK=FUNIJK(I,J,K)
2194                          IF(F_AT(IJK)==UNDEFINED) THEN
2195                             F_AT(IJK)=CURRENT_F
2196                          ELSEIF(F_AT(IJK)/=ZERO) THEN
2197                             CURRENT_F = F_AT(IJK)
2198                          ENDIF
2199                       ENDDO
2200                    ENDIF
2201                 ENDDO
2202              ENDDO
2203     
2204              call send_recv(F_AT,2)
2205     
2206     
2207           CASE ('KIJ')
2208     
2209     
2210     
2211              DO K=KSTART3,KEND3
2212                 DO J=JSTART3,JEND3
2213                    F_FOUND = .FALSE.
2214                    DO I=ISTART3,IEND3
2215                       IJK=FUNIJK(I,J,K)
2216                       IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2217                          F_FOUND = .TRUE.
2218                          CURRENT_F = F_AT(IJK)
2219                          EXIT
2220                       ENDIF
2221                    ENDDO
2222                    IF(F_FOUND) THEN
2223                       DO I=ISTART3,IEND3
2224                          IJK=FUNIJK(I,J,K)
2225                          IF(F_AT(IJK)==UNDEFINED) THEN
2226                             F_AT(IJK)=CURRENT_F
2227                          ELSEIF(F_AT(IJK)/=ZERO) THEN
2228                             CURRENT_F = F_AT(IJK)
2229                          ENDIF
2230                       ENDDO
2231                    ENDIF
2232                 ENDDO
2233              ENDDO
2234     
2235              call send_recv(F_AT,2)
2236     
2237              DO I=ISTART3,IEND3
2238                 DO K=KSTART3,KEND3
2239                    F_FOUND = .FALSE.
2240                    DO J=JSTART3,JEND3
2241                       IJK=FUNIJK(I,J,K)
2242                       IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2243                          F_FOUND = .TRUE.
2244                          CURRENT_F = F_AT(IJK)
2245                          EXIT
2246                       ENDIF
2247                    ENDDO
2248                    IF(F_FOUND) THEN
2249                       DO J=JSTART3,JEND3
2250                          IJK=FUNIJK(I,J,K)
2251                          IF(F_AT(IJK)==UNDEFINED) THEN
2252                             F_AT(IJK)=CURRENT_F
2253                          ELSEIF(F_AT(IJK)/=ZERO) THEN
2254                             CURRENT_F = F_AT(IJK)
2255                          ENDIF
2256                       ENDDO
2257                    ENDIF
2258                 ENDDO
2259              ENDDO
2260     
2261     
2262              call send_recv(F_AT,2)
2263     
2264              DO J=JSTART3,JEND3
2265                 DO I=ISTART3,IEND3
2266                    F_FOUND = .FALSE.
2267                    DO K=KSTART3,KEND3
2268                       IJK=FUNIJK(I,J,K)
2269                       IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2270                          F_FOUND = .TRUE.
2271                          CURRENT_F = F_AT(IJK)
2272                          EXIT
2273                       ENDIF
2274                    ENDDO
2275                    IF(F_FOUND) THEN
2276                       DO K=KSTART3,KEND3
2277                          IJK=FUNIJK(I,J,K)
2278                          IF(F_AT(IJK)==UNDEFINED) THEN
2279                             F_AT(IJK)=CURRENT_F
2280                          ELSEIF(F_AT(IJK)/=ZERO) THEN
2281                             CURRENT_F = F_AT(IJK)
2282                          ENDIF
2283                       ENDDO
2284                    ENDIF
2285                 ENDDO
2286              ENDDO
2287     
2288     
2289              call send_recv(F_AT,2)
2290     
2291              CASE DEFAULT
2292                 IF(myPE == PE_IO) THEN
2293                    WRITE(*,*)'CAD_INTERSECT.'
2294                    WRITE(*,*)'UNKNOWN CAD_PROPAGATE_ORDER:',CAD_PROPAGATE_ORDER
2295                    WRITE(*,*)'ACCEPTABLE VALUES:'
2296                    WRITE(*,*)'IJK'
2297                    WRITE(*,*)'JKI'
2298                    WRITE(*,*)'KIJ'
2299                 ENDIF
2300     !            CALL MFIX_EXIT(myPE)
2301     
2302           END SELECT
2303     
2304     
2305     
2306           call send_recv(F_AT,2)
2307     
2308           RETURN
2309     
2310           END SUBROUTINE CAD_INTERSECT
2311     
2312     
2313     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2314     !                                                                      C
2315     !  Module name: ADD_FACET                                              C
2316     !  Purpose: Add facet to list in IJK scalar cell                       C
2317     !                                                                      C
2318     !  Author: Jeff Dietiker                              Date: 15-Oct-13  C
2319     !  Reviewer:                                          Date:            C
2320     !                                                                      C
2321     !  Revision Number #                                  Date: ##-###-##  C
2322     !  Author: #                                                           C
2323     !  Purpose: #                                                          C
2324     !                                                                      C
2325     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
2326       SUBROUTINE ADD_FACET_AND_SET_BC_ID(IJK,N)
2327     
2328           USE param
2329           USE param1
2330           USE parallel
2331           USE constant
2332           USE run
2333           USE toleranc
2334           USE geometry
2335           USE indices
2336           USE compar
2337           USE sendrecv
2338           USE quadric
2339           USE cutcell
2340           USE polygon
2341           USE stl
2342           USE mpi_utility
2343     
2344           IMPLICIT NONE
2345           INTEGER :: IJK,N
2346     
2347     
2348           BC_ID(IJK) = BC_ID_STL_FACE(N)             ! Set tentative BC_ID
2349     
2350           IF(N_FACET_AT(IJK)<DIM_FACETS_PER_CELL) THEN
2351     
2352              N_FACET_AT(IJK) = N_FACET_AT(IJK) + 1
2353              LIST_FACET_AT(IJK,N_FACET_AT(IJK)) = N
2354     
2355           ELSE
2356     
2357              WRITE(*,*) ' FATAL ERROR: TOO MANY FACETS IN CELL: ', IJK
2358              CALL MFIX_EXIT(myPE)
2359     
2360           ENDIF
2361     
2362     
2363     
2364           END SUBROUTINE ADD_FACET_AND_SET_BC_ID
2365     
2366     
2367