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