File: RELATIVE:/../../../mfix.git/model/cartesian_grid/get_connectivity.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: GET_CONNECTIVITY                                       C
4     !  Purpose: Set flags for saclar cut cells, based on intersection      C
5     !  of the grid with the quadric(s)                                     C
6     !                                                                      C
7     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  Revision Number #                                  Date: ##-###-##  C
11     !  Author: #                                                           C
12     !  Purpose: #                                                          C
13     !                                                                      C
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
15       SUBROUTINE GET_CONNECTIVITY(IJK,TYPE_OF_CELL,N_NEW_POINTS,N_NODES,CONNECT,X_NP,Y_NP,Z_NP,TOTAL_NUMBER_OF_INTERSECTIONS,&
16                  X_intersect,Y_intersect,Z_intersect)
17     
18           USE compar, ONLY: ijkend3
19           USE cutcell
20           USE functions, ONLY: FUNIJK
21           USE geometry, ONLY: DO_K, NO_K
22           USE indices, ONLY: I_OF, J_OF, K_OF
23           USE polygon, ONLY: n_polygon
24           USE quadric, ONLY: tol_f
25     
26           IMPLICIT NONE
27           INTEGER :: I,J,K,IM,JM,KM
28           INTEGER :: IJK,IMJK,IJMK,IJKM,IMJMK,IMJKM,IJMKM,IMJMKM
29           LOGICAL :: CLIP_FLAG
30           LOGICAL, DIMENSION(8) :: CORNER_INTERSECTION
31           INTEGER :: TOTAL_NUMBER_OF_INTERSECTIONS,NUMBER_OF_EDGE_INTERSECTIONS
32           INTEGER :: NUMBER_OF_CORNER_INTERSECTIONS,MAX_CORNER_INTERSECTIONS
33           INTEGER :: N_NODES,N_NEW_POINTS
34           INTEGER :: NODE,N_N1,N_N2,Q_ID,BCID
35           INTEGER, DIMENSION(DIMENSION_3,15) :: CONNECT
36           DOUBLE PRECISION, DIMENSION(DIMENSION_MAX_CUT_CELL) :: X_NP,Y_NP,Z_NP
37           DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: X_intersect,Y_intersect,Z_intersect
38           CHARACTER (LEN=*) :: TYPE_OF_CELL
39     
40     !======================================================================
41     !  Get coordinates of eight nodes
42     !======================================================================
43     
44           CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
45     
46           I = I_OF(IJK)
47           J = J_OF(IJK)
48           K = K_OF(IJK)
49     
50           IM = I - 1
51           JM = J - 1
52           KM = K - 1
53     
54           IMJK   = FUNIJK(IM,J,K)
55           IJMK   = FUNIJK(I,JM,K)
56           IJKM   = FUNIJK(I,J,KM)
57     
58           IMJMK  = FUNIJK(IM,JM,K)
59           IMJKM  = FUNIJK(IM,J,KM)
60           IJMKM  = FUNIJK(I,JM,KM)
61     
62           IMJMKM = FUNIJK(IM,JM,KM)
63     
64     
65     !======================================================================
66     !  Evaluate Quadric at all corners
67     !======================================================================
68     
69              N_NODES = 0
70     
71              CORNER_INTERSECTION = .FALSE.
72              NUMBER_OF_CORNER_INTERSECTIONS = 0
73              NUMBER_OF_EDGE_INTERSECTIONS = 0
74     
75              IF(NO_K) THEN
76                 N_N1 = 5
77                 N_N2 = 8
78              ELSE
79                 N_N1 = 1
80                 N_N2 = 8
81              ENDIF
82     
83              DO NODE = N_N1,N_N2
84     
85                 Q_ID = 1
86                 CALL EVAL_F('QUADRIC',X_NODE(NODE),Y_NODE(NODE),Z_NODE(NODE),Q_ID,F_NODE(NODE),CLIP_FLAG)
87     
88                 CALL EVAL_F('POLYGON',X_NODE(NODE),Y_NODE(NODE),Z_NODE(NODE),N_POLYGON,F_NODE(NODE),CLIP_FLAG)
89     
90                 CALL EVAL_F('USR_DEF',X_NODE(NODE),Y_NODE(NODE),Z_NODE(NODE),N_USR_DEF,F_NODE(NODE),CLIP_FLAG)
91     
92                 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,NODE,F_NODE(NODE),CLIP_FLAG,BCID)
93     
94                 IF (ABS(F_NODE(NODE)) < TOL_F ) THEN
95                    CORNER_INTERSECTION(NODE) = .TRUE.
96                    NUMBER_OF_CORNER_INTERSECTIONS = NUMBER_OF_CORNER_INTERSECTIONS + 1
97                    N_NODES = N_NODES + 1
98                    CONNECT(IJK,N_NODES) = IJK_OF_NODE(NODE)
99                 ENDIF
100     
101     
102                 IF(SNAP(IJK_OF_NODE(NODE))) THEN
103                    CORNER_INTERSECTION(NODE) = .TRUE.
104                    NUMBER_OF_CORNER_INTERSECTIONS = NUMBER_OF_CORNER_INTERSECTIONS + 1
105                    N_NODES = N_NODES + 1
106                    CONNECT(IJK,N_NODES) = IJK_OF_NODE(NODE)
107                 ENDIF
108     
109     
110              END DO
111     
112     !======================================================================
113     !  Count the number of edge intersections (excluding corner intersections)
114     !  For each new edge intersection found:
115     !     - Increment the total number of points
116     !     - Store the location of the additional point
117     !======================================================================
118     
119                 IF(DO_K) THEN
120     
121                    IF(INTERSECT_X(IJMKM)) THEN  ! Edge 1  = Nodes 1-2
122                       IF((.NOT.CORNER_INTERSECTION(1)).AND.(.NOT.CORNER_INTERSECTION(2))) THEN
123                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
124                          N_NEW_POINTS = N_NEW_POINTS + 1
125                          X_NP(N_NEW_POINTS) = X_intersect(IJMKM)
126                          Y_NP(N_NEW_POINTS) = Y_NODE(1)
127                          Z_NP(N_NEW_POINTS) = Z_NODE(1)
128                          N_NODES = N_NODES + 1
129                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
130                       ENDIF
131                    ENDIF
132     
133                    IF(INTERSECT_Y(IJKM)) THEN  ! Edge 2  = Nodes 2-4
134                       IF((.NOT.CORNER_INTERSECTION(2)).AND.(.NOT.CORNER_INTERSECTION(4))) THEN
135                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
136                          N_NEW_POINTS = N_NEW_POINTS + 1
137                          X_NP(N_NEW_POINTS) = X_NODE(2)
138                          Y_NP(N_NEW_POINTS) = Y_intersect(IJKM)
139                          Z_NP(N_NEW_POINTS) = Z_NODE(2)
140                          N_NODES = N_NODES + 1
141                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
142                       ENDIF
143                    ENDIF
144     
145                    IF(INTERSECT_X(IJKM)) THEN  ! Edge 3  = Nodes 3-4
146                       IF((.NOT.CORNER_INTERSECTION(3)).AND.(.NOT.CORNER_INTERSECTION(4))) THEN
147                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
148                          N_NEW_POINTS = N_NEW_POINTS + 1
149                          X_NP(N_NEW_POINTS) = X_intersect(IJKM)
150                          Y_NP(N_NEW_POINTS) = Y_NODE(3)
151                          Z_NP(N_NEW_POINTS) = Z_NODE(3)
152                          N_NODES = N_NODES + 1
153                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
154                       ENDIF
155                    ENDIF
156     
157                    IF(INTERSECT_Y(IMJKM)) THEN  ! Edge 4  = Nodes 1-3
158                       IF((.NOT.CORNER_INTERSECTION(1)).AND.(.NOT.CORNER_INTERSECTION(3))) THEN
159                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
160                          N_NEW_POINTS = N_NEW_POINTS + 1
161                          X_NP(N_NEW_POINTS) = X_NODE(1)
162                          Y_NP(N_NEW_POINTS) = Y_intersect(IMJKM)
163                          Z_NP(N_NEW_POINTS) = Z_NODE(1)
164                          N_NODES = N_NODES + 1
165                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
166                       ENDIF
167                    ENDIF
168     
169                 ENDIF
170     
171     
172                 IF(INTERSECT_X(IJMK)) THEN  ! Edge 5  = Nodes 5-6
173                    IF((.NOT.CORNER_INTERSECTION(5)).AND.(.NOT.CORNER_INTERSECTION(6))) THEN
174                       NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
175                       N_NEW_POINTS = N_NEW_POINTS + 1
176                       X_NP(N_NEW_POINTS) = X_intersect(IJMK)
177                       Y_NP(N_NEW_POINTS) = Y_NODE(5)
178                       Z_NP(N_NEW_POINTS) = Z_NODE(5)
179                       N_NODES = N_NODES + 1
180                       CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
181                    ENDIF
182                 ENDIF
183     
184                 IF(INTERSECT_Y(IJK)) THEN  ! Edge 6  = Nodes 6-8
185                    IF((.NOT.CORNER_INTERSECTION(6)).AND.(.NOT.CORNER_INTERSECTION(8))) THEN
186                       NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
187                       N_NEW_POINTS = N_NEW_POINTS + 1
188                       X_NP(N_NEW_POINTS) = X_NODE(6)
189                       Y_NP(N_NEW_POINTS) = Y_intersect(IJK)
190                       Z_NP(N_NEW_POINTS) = Z_NODE(6)
191                       N_NODES = N_NODES + 1
192                       CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
193                    ENDIF
194                 ENDIF
195     
196                 IF(INTERSECT_X(IJK)) THEN  ! Edge 7  = Nodes 7-8
197                    IF((.NOT.CORNER_INTERSECTION(7)).AND.(.NOT.CORNER_INTERSECTION(8))) THEN
198                       NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
199                       N_NEW_POINTS = N_NEW_POINTS + 1
200                       X_NP(N_NEW_POINTS) = X_intersect(IJK)
201                       Y_NP(N_NEW_POINTS) = Y_NODE(7)
202                       Z_NP(N_NEW_POINTS) = Z_NODE(7)
203                       N_NODES = N_NODES + 1
204                       CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
205                    ENDIF
206                 ENDIF
207     
208                 IF(INTERSECT_Y(IMJK)) THEN  ! Edge 8  = Nodes 5-7
209                    IF((.NOT.CORNER_INTERSECTION(5)).AND.(.NOT.CORNER_INTERSECTION(7))) THEN
210                       NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
211                       N_NEW_POINTS = N_NEW_POINTS + 1
212                       X_NP(N_NEW_POINTS) = X_NODE(5)
213                       Y_NP(N_NEW_POINTS) = Y_intersect(IMJK)
214                       Z_NP(N_NEW_POINTS) = Z_NODE(5)
215                       N_NODES = N_NODES + 1
216                       CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
217                    ENDIF
218                 ENDIF
219     
220     
221                 IF(DO_K) THEN
222     
223                    IF(INTERSECT_Z(IMJMK)) THEN  ! Edge 9  = Nodes 1-5
224                       IF((.NOT.CORNER_INTERSECTION(1)).AND.(.NOT.CORNER_INTERSECTION(5))) THEN
225                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
226                          N_NEW_POINTS = N_NEW_POINTS + 1
227                          X_NP(N_NEW_POINTS) = X_NODE(1)
228                          Y_NP(N_NEW_POINTS) = Y_NODE(1)
229                          Z_NP(N_NEW_POINTS) = Z_intersect(IMJMK)
230                          N_NODES = N_NODES + 1
231                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
232                       ENDIF
233                    ENDIF
234     
235                    IF(INTERSECT_Z(IJMK))  THEN  ! Edge 10 = Nodes 2-6
236                       IF((.NOT.CORNER_INTERSECTION(2)).AND.(.NOT.CORNER_INTERSECTION(6))) THEN
237                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
238                          N_NEW_POINTS = N_NEW_POINTS + 1
239                          X_NP(N_NEW_POINTS) = X_NODE(2)
240                          Y_NP(N_NEW_POINTS) = Y_NODE(2)
241                          Z_NP(N_NEW_POINTS) = Z_intersect(IJMK)
242                          N_NODES = N_NODES + 1
243                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
244                       ENDIF
245                    ENDIF
246     
247                    IF(INTERSECT_Z(IJK))  THEN  ! Edge 11 = Nodes 4-8
248                       IF((.NOT.CORNER_INTERSECTION(4)).AND.(.NOT.CORNER_INTERSECTION(8))) THEN
249                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
250                          N_NEW_POINTS = N_NEW_POINTS + 1
251                          X_NP(N_NEW_POINTS) = X_NODE(4)
252                          Y_NP(N_NEW_POINTS) = Y_NODE(4)
253                          Z_NP(N_NEW_POINTS) = Z_intersect(IJK)
254                          N_NODES = N_NODES + 1
255                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
256                       ENDIF
257                    ENDIF
258     
259                    IF(INTERSECT_Z(IMJK)) THEN  ! Edge 12 = Nodes 3-7
260                       IF((.NOT.CORNER_INTERSECTION(3)).AND.(.NOT.CORNER_INTERSECTION(7))) THEN
261                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
262                          N_NEW_POINTS = N_NEW_POINTS + 1
263                          X_NP(N_NEW_POINTS) = X_NODE(3)
264                          Y_NP(N_NEW_POINTS) = Y_NODE(3)
265                          Z_NP(N_NEW_POINTS) = Z_intersect(IMJK)
266                          N_NODES = N_NODES + 1
267                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
268                       ENDIF
269                    ENDIF
270     
271                 ENDIF
272     
273     !======================================================================
274     !  Count the total number of intersections (corner and edge intersections)
275     !======================================================================
276     
277                 IF(NO_K) THEN
278                    MAX_CORNER_INTERSECTIONS = 2
279                 ELSE
280                    MAX_CORNER_INTERSECTIONS = 4
281                 ENDIF
282     
283     
284                IF(NUMBER_OF_CORNER_INTERSECTIONS > MAX_CORNER_INTERSECTIONS) THEN
285                   IF(PRINT_WARNINGS) THEN
286                      WRITE(*,*)'WARNING:',NUMBER_OF_CORNER_INTERSECTIONS,&
287                             ' CORNER INTERSECTIONS DETECTED IN CELL IJK=',IJK
288                      WRITE(*,*)'THIS USUALLY INDICATE A FALSE CUT-CELL (CORNER CELL)'
289     !                 WRITE(*,*)'RESETTING NUMBER_OF_CORNER_INTERSECTIONS TO', MAX_CORNER_INTERSECTIONS
290                      WRITE(*,*)'REMOVING CUT CELL'
291                   ENDIF
292     !              NUMBER_OF_CORNER_INTERSECTIONS = MAX_CORNER_INTERSECTIONS
293     
294                   ! Force the total number of intersections to be zero, and therefore, the cell will be considered as a non-cut cell
295     !              NUMBER_OF_CORNER_INTERSECTIONS = -NUMBER_OF_EDGE_INTERSECTIONS
296     
297                   ! Force the total number of intersections to be -one, and therefore, the cell will be considered as a non-cut cell
298                   NUMBER_OF_CORNER_INTERSECTIONS = -NUMBER_OF_EDGE_INTERSECTIONS -1
299     
300     
301                ENDIF
302     
303                 TOTAL_NUMBER_OF_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + NUMBER_OF_CORNER_INTERSECTIONS
304     
305           RETURN
306     
307           END SUBROUTINE GET_CONNECTIVITY
308     
309     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
310     !                                                                      C
311     !  Module name: GET_CELL_NODE_COORDINATES                              C
312     !  Purpose: Get the cell corners (x,y,z) coordinates                   C
313     !                                                                      C
314     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
315     !  Reviewer:                                          Date:            C
316     !                                                                      C
317     !  Revision Number #                                  Date: ##-###-##  C
318     !  Author: #                                                           C
319     !  Purpose: #                                                          C
320     !                                                                      C
321     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
322       SUBROUTINE GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
323     
324           USE compar, ONLY: mype
325           USE cutcell
326           USE functions, ONLY: FUNIJK
327           USE geometry, ONLY: NO_K, dx, dy, dz
328           USE indices, ONLY: I_OF, J_OF, K_OF
329     
330           IMPLICIT NONE
331           CHARACTER (LEN=*) :: TYPE_OF_CELL
332           DOUBLE PRECISION:: Xw,Xe,Yn,Ys,Zb,Zt
333           INTEGER :: I,J,K,IP,JP,KP,IM,JM,KM
334           INTEGER :: IJK,IMJK,IJMK,IJKM,IMJMK,IMJKM,IJMKM,IMJMKM
335     
336           I = I_OF(IJK)
337           J = J_OF(IJK)
338           K = K_OF(IJK)
339     
340           IP = I + 1
341           JP = J + 1
342           KP = K + 1
343     
344           IM = I - 1
345           JM = J - 1
346           KM = K - 1
347     
348           IMJK   = FUNIJK(IM,J,K)
349           IJMK   = FUNIJK(I,JM,K)
350           IJKM   = FUNIJK(I,J,KM)
351     
352           IMJMK  = FUNIJK(IM,JM,K)
353           IMJKM  = FUNIJK(IM,J,KM)
354           IJMKM  = FUNIJK(I,JM,KM)
355     
356           IMJMKM = FUNIJK(IM,JM,KM)
357     
358           IJK_OF_NODE(0) = IJK
359           IJK_OF_NODE(1) = IMJMKM
360           IJK_OF_NODE(2) = IJMKM
361           IJK_OF_NODE(3) = IMJKM
362           IJK_OF_NODE(4) = IJKM
363           IJK_OF_NODE(5) = IMJMK
364           IJK_OF_NODE(6) = IJMK
365           IJK_OF_NODE(7) = IMJK
366           IJK_OF_NODE(8) = IJK
367           IJK_OF_NODE(15) = IJK
368     
369     
370           SELECT CASE (TYPE_OF_CELL)
371              CASE('SCALAR')
372                 Xw = XG_E(I) - DX(I)          ! west face location
373                 Xe = XG_E(I)                  ! east face location
374     
375                 Ys = YG_N(J) - DY(J)          ! south face location
376                 Yn = YG_N(J)                  ! north face location
377     
378                 IF(NO_K) THEN
379                    Zb = ZERO                  ! bottom face location
380                    Zt = ZERO                  ! top face location
381                 ELSE
382                    Zb = ZG_T(K) - DZ(K)       ! bottom face location
383                    Zt = ZG_T(K)               ! top face location
384                 ENDIF
385     
386              CASE('U_MOMENTUM')
387                 Xw = XG_E(I) - HALF * DX(I)   ! west face location
388                 Xe = XG_E(I) + HALF * DX(IP)  ! east face location
389     
390                 Ys = YG_N(J) - DY(J)          ! south face location
391                 Yn = YG_N(J)                  ! north face location
392     
393                 IF(NO_K) THEN
394                    Zb = ZERO                  ! bottom face location
395                    Zt = ZERO                  ! top face location
396                 ELSE
397                    Zb = ZG_T(K) - DZ(K)       ! bottom face location
398                    Zt = ZG_T(K)               ! top face location
399                 ENDIF
400     
401              CASE('V_MOMENTUM')
402                 Xw = XG_E(I) - DX(I)          ! west face location
403                 Xe = XG_E(I)                  ! east face location
404     
405                 Ys = YG_N(J) - HALF * DY(J)   ! south face location
406                 Yn = YG_N(J) + HALF * DY(JP)  ! north face location
407     
408                 IF(NO_K) THEN
409                    Zb = ZERO                  ! bottom face location
410                    Zt = ZERO                  ! top face location
411                 ELSE
412                    Zb = ZG_T(K) - DZ(K)       ! bottom face location
413                    Zt = ZG_T(K)               ! top face location
414                 ENDIF
415     
416              CASE('W_MOMENTUM')
417                 Xw = XG_E(I) - DX(I)          ! west face location
418                 Xe = XG_E(I)                  ! east face location
419     
420                 Ys = YG_N(J) - DY(J)          ! south face location
421                 Yn = YG_N(J)                  ! north face location
422     
423                 Zb = ZG_T(K) - HALF * DZ(K)   ! bottom face location
424                 Zt = ZG_T(K) + HALF * DZ(KP)  ! top face location
425     
426     
427              CASE DEFAULT
428                 WRITE(*,*)'SUBROUTINE: GET_CELL_NODE_COORDINATES'
429                 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
430                 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
431                 WRITE(*,*)'SCALAR'
432                 WRITE(*,*)'U_MOMENTUM'
433                 WRITE(*,*)'V_MOMENTUM'
434                 WRITE(*,*)'W_MOMENTUM'
435                 CALL MFIX_EXIT(myPE)
436           END SELECT
437     
438     !     CELL CENTER :
439           X_NODE(0) = HALF * ( xw + xe )
440           Y_NODE(0) = HALF * ( ys + yn )
441           Z_NODE(0) = HALF * ( zb + zt )
442     
443     !     NODE 1 : IM,JM,KM
444           X_NODE(1) = xw
445           Y_NODE(1) = ys
446           Z_NODE(1) = zb
447     
448     !     NODE 2 : I,JM,KM
449           X_NODE(2) = xe
450           Y_NODE(2) = ys
451           Z_NODE(2) = zb
452     
453     !     NODE 3 : IM,J,KM
454           X_NODE(3) = xw
455           Y_NODE(3) = yn
456           Z_NODE(3) = zb
457     
458     !     NODE 4 : I,J,KM
459           X_NODE(4) = xe
460           Y_NODE(4) = yn
461           Z_NODE(4) = zb
462     
463     !     NODE 5 : IM,JM,K
464           X_NODE(5) = xw
465           Y_NODE(5) = ys
466           Z_NODE(5) = zt
467     
468     !     NODE 6 : I,JM,K
469           X_NODE(6) = xe
470           Y_NODE(6) = ys
471           Z_NODE(6) = zt
472     
473     !     NODE 7 : IM,J,K
474           X_NODE(7) = xw
475           Y_NODE(7) = yn
476           Z_NODE(7) = zt
477     
478     !     NODE 8 : I,J,K
479           X_NODE(8) = xe
480           Y_NODE(8) = yn
481           Z_NODE(8) = zt
482     
483           RETURN
484     
485           END SUBROUTINE GET_CELL_NODE_COORDINATES
486     
487     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
488     !                                                                      C
489     !  Module name: GET_GLOBAL_CELL_NODE_COORDINATES                       C
490     !  Purpose: Get the cell corners (x,y,z) coordinates                   C
491     !                                                                      C
492     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
493     !  Reviewer:                                          Date:            C
494     !                                                                      C
495     !  Revision Number #                                  Date: ##-###-##  C
496     !  Author: #                                                           C
497     !  Purpose: #                                                          C
498     !                                                                      C
499     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
500       SUBROUTINE GET_GLOBAL_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
501     
502           USE compar, ONLY: MYPE
503           USE cutcell
504           USE functions, ONLY: FUNIJK_GL
505           USE geometry, ONLY: NO_K, dx, dy, dz
506           USE vtk, ONLY: GLOBAL_I_OF, GLOBAL_J_OF, GLOBAL_K_OF
507     
508           IMPLICIT NONE
509           CHARACTER (LEN=*) :: TYPE_OF_CELL
510           DOUBLE PRECISION:: Xw,Xe,Yn,Ys,Zb,Zt
511           INTEGER :: I,J,K,IP,JP,KP,IM,JM,KM
512           INTEGER :: IJK,IMJK,IJMK,IJKM,IMJMK,IMJKM,IJMKM,IMJMKM
513     
514           I = GLOBAL_I_OF(IJK)
515           J = GLOBAL_J_OF(IJK)
516           K = GLOBAL_K_OF(IJK)
517     
518           IP = I + 1
519           JP = J + 1
520           KP = K + 1
521     
522           IM = I - 1
523           JM = J - 1
524           KM = K - 1
525     
526           IMJK   = FUNIJK_GL(IM,J,K)
527           IJMK   = FUNIJK_GL(I,JM,K)
528           IJKM   = FUNIJK_GL(I,J,KM)
529     
530           IMJMK  = FUNIJK_GL(IM,JM,K)
531           IMJKM  = FUNIJK_GL(IM,J,KM)
532           IJMKM  = FUNIJK_GL(I,JM,KM)
533     
534           IMJMKM = FUNIJK_GL(IM,JM,KM)
535     
536           IJK_OF_NODE(1) = IMJMKM
537           IJK_OF_NODE(2) = IJMKM
538           IJK_OF_NODE(3) = IMJKM
539           IJK_OF_NODE(4) = IJKM
540           IJK_OF_NODE(5) = IMJMK
541           IJK_OF_NODE(6) = IJMK
542           IJK_OF_NODE(7) = IMJK
543           IJK_OF_NODE(8) = IJK
544     
545           SELECT CASE (TYPE_OF_CELL)
546              CASE('SCALAR')
547                 Xw = XG_E(I) - DX(I)          ! west face location
548                 Xe = XG_E(I)                  ! east face location
549     
550                 Ys = YG_N(J) - DY(J)          ! south face location
551                 Yn = YG_N(J)                  ! north face location
552     
553                 IF(NO_K) THEN
554                    Zb = ZERO                  ! bottom face location
555                    Zt = ZERO                  ! top face location
556                 ELSE
557                    Zb = ZG_T(K) - DZ(K)       ! bottom face location
558                    Zt = ZG_T(K)               ! top face location
559                 ENDIF
560     
561              CASE('U_MOMENTUM')
562                 Xw = XG_E(I) - HALF * DX(I)   ! west face location
563                 Xe = XG_E(I) + HALF * DX(IP)  ! east face location
564     
565                 Ys = YG_N(J) - DY(J)          ! south face location
566                 Yn = YG_N(J)                  ! north face location
567     
568                 IF(NO_K) THEN
569                    Zb = ZERO                  ! bottom face location
570                    Zt = ZERO                  ! top face location
571                 ELSE
572                    Zb = ZG_T(K) - DZ(K)       ! bottom face location
573                    Zt = ZG_T(K)               ! top face location
574                 ENDIF
575     
576              CASE('V_MOMENTUM')
577                 Xw = XG_E(I) - DX(I)          ! west face location
578                 Xe = XG_E(I)                  ! east face location
579     
580                 Ys = YG_N(J) - HALF * DY(J)   ! south face location
581                 Yn = YG_N(J) + HALF * DY(JP)  ! north face location
582     
583                 IF(NO_K) THEN
584                    Zb = ZERO                  ! bottom face location
585                    Zt = ZERO                  ! top face location
586                 ELSE
587                    Zb = ZG_T(K) - DZ(K)       ! bottom face location
588                    Zt = ZG_T(K)               ! top face location
589                 ENDIF
590     
591              CASE('W_MOMENTUM')
592                 Xw = XG_E(I) - DX(I)          ! west face location
593                 Xe = XG_E(I)                  ! east face location
594     
595                 Ys = YG_N(J) - DY(J)          ! south face location
596                 Yn = YG_N(J)                  ! north face location
597     
598                 Zb = ZG_T(K) - HALF * DZ(K)   ! bottom face location
599                 Zt = ZG_T(K) + HALF * DZ(KP)  ! top face location
600     
601     
602              CASE DEFAULT
603                 WRITE(*,*)'SUBROUTINE: GET_CELL_NODE_COORDINATES'
604                 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
605                 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
606                 WRITE(*,*)'SCALAR'
607                 WRITE(*,*)'U_MOMENTUM'
608                 WRITE(*,*)'V_MOMENTUM'
609                 WRITE(*,*)'W_MOMENTUM'
610                 CALL MFIX_EXIT(myPE)
611           END SELECT
612     
613     !     CELL CENTER :
614           X_NODE(0) = HALF * ( xw + xe )
615           Y_NODE(0) = HALF * ( ys + yn )
616           Z_NODE(0) = HALF * ( zb + zt )
617     
618     !     NODE 1 : IM,JM,KM
619           X_NODE(1) = xw
620           Y_NODE(1) = ys
621           Z_NODE(1) = zb
622     
623     !     NODE 2 : I,JM,KM
624           X_NODE(2) = xe
625           Y_NODE(2) = ys
626           Z_NODE(2) = zb
627     
628     !     NODE 3 : IM,J,KM
629           X_NODE(3) = xw
630           Y_NODE(3) = yn
631           Z_NODE(3) = zb
632     
633     !     NODE 4 : I,J,KM
634           X_NODE(4) = xe
635           Y_NODE(4) = yn
636           Z_NODE(4) = zb
637     
638     !     NODE 5 : IM,JM,K
639           X_NODE(5) = xw
640           Y_NODE(5) = ys
641           Z_NODE(5) = zt
642     
643     !     NODE 6 : I,JM,K
644           X_NODE(6) = xe
645           Y_NODE(6) = ys
646           Z_NODE(6) = zt
647     
648     !     NODE 7 : IM,J,K
649           X_NODE(7) = xw
650           Y_NODE(7) = yn
651           Z_NODE(7) = zt
652     
653     !     NODE 8 : I,J,K
654           X_NODE(8) = xe
655           Y_NODE(8) = yn
656           Z_NODE(8) = zt
657     
658           RETURN
659     
660           END SUBROUTINE GET_GLOBAL_CELL_NODE_COORDINATES
661