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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: WRITE_DBG_VTU_AND_VTP_FILES                            C
4     !  Purpose: Writes the cell and particle data in VTK format            C
5     !           for debug regions only.                                    C
6     !                                                                      C
7     !  Author: Jeff Dietiker                              Date: 22-Jul-15  C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  Revision Number #                                  Date: ##-###-##  C
11     !  Author: #                                                           C
12     !  Purpose: #                                                          C
13     !                                                                      C
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
15           SUBROUTINE WRITE_DBG_VTU_AND_VTP_FILES
16     
17           use discretelement, only: DISCRETE_ELEMENT
18           use vtp, only: write_vtp_file
19           use vtk, only: DIMENSION_VTK
20     
21           IMPLICIT NONE
22           INTEGER :: LC
23     
24           DO LC = 1, DIMENSION_VTK
25              CALL WRITE_VTU_FILE(LC,1)
26              IF(DISCRETE_ELEMENT) CALL WRITE_VTP_FILE(LC,1)
27           ENDDO
28     
29           END SUBROUTINE WRITE_DBG_VTU_AND_VTP_FILES
30     
31     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
32     !                                                                      C
33     !  Module name: WRITE_VTU_FILE                                         C
34     !  Purpose: Writes the cell data grid in VTK format (Unstructured VTU) C
35     !           Binary format                                              C
36     !                                                                      C
37     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
38     !  Reviewer:                                          Date:            C
39     !                                                                      C
40     !  Revision Number #                                  Date: ##-###-##  C
41     !  Author: #                                                           C
42     !  Purpose: #                                                          C
43     !                                                                      C
44     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
45       SUBROUTINE WRITE_VTU_FILE(LCV,MODE)
46     
47           USE compar
48           USE constant
49           USE cutcell
50           USE discretelement, Only :  DISCRETE_ELEMENT
51           USE fldvar
52           USE functions
53           USE geometry
54           USE indices
55           USE mfix_pic
56           USE mpi_utility
57           USE output
58           USE parallel
59           USE parallel_mpi
60           USE param
61           USE param1
62           USE pgcor
63           USE pgcor
64           USE physprop
65           USE pscor
66           USE quadric
67           USE run
68           USE rxns
69           USE scalars
70           USE sendrecv
71           USE stl
72           USE toleranc
73           USE visc_s
74           USE vtk
75           USE vtp
76     
77           IMPLICIT NONE
78           INTEGER :: I,J,K,L,M,N,R,IJK,LCV
79     
80           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE ::  FACET_COUNT_DES, NEIGHBORING_FACET
81     
82           INTEGER :: SPECIES_COUNTER,LT
83     
84           CHARACTER (LEN=32) :: SUBM,SUBN,SUBR
85           CHARACTER (LEN=64) :: VAR_NAME
86     
87           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE ::  DP_BC_ID, IJK_ARRAY
88     
89           INTEGER :: PASS
90           INTEGER :: WRITE_HEADER = 1
91           INTEGER :: WRITE_DATA   = 2
92           INTEGER :: MODE   ! MODE = 0 : Write regular VTK region file
93                             ! MODE = 1 : Write debug   VTK region file (VTK_DBG_FILE = .TRUE.)
94     
95           ! There is nothing to write if we are not in adefined vtk region
96           VTK_REGION = LCV
97           IF(.NOT.VTK_DEFINED(VTK_REGION)) RETURN
98           IF(VTK_DATA(LCV)/='C') RETURN
99           IF(MODE==0.AND.(VTK_DBG_FILE(LCV))) RETURN
100           IF(MODE==1.AND.(.NOT.VTK_DBG_FILE(LCV))) RETURN
101     
102     !     Location of U-momentum cells for original (uncut grid)
103           IF (DO_I) THEN
104             XG_E(1) = ZERO
105             DO I = IMIN1, IMAX2
106                XG_E(I) = XG_E(I-1) + DX(I)
107             END DO
108           ENDIF
109     
110     !     Location of V-momentum cells for original (uncut grid)
111           IF (DO_J) THEN
112             YG_N(1) = ZERO
113             DO J = JMIN1, JMAX2
114                YG_N(J) = YG_N(J-1) + DY(J)
115             END DO
116           ENDIF
117     
118     !     Location of W-momentum cells for original (uncut grid)
119           IF (DO_K) THEN
120             ZG_T(1) = ZERO
121             DO K = KMIN1, KMAX2
122                ZG_T(K) = ZG_T(K-1) + DZ(K)
123             END DO
124           ELSE
125              ZG_T = ZERO
126           ENDIF
127     
128     
129           CALL SETUP_VTK_REGION
130     
131           CALL OPEN_VTU_FILE_BIN(MODE)
132     
133           IF(MODE==0) CALL OPEN_PVD_FILE
134     
135           CALL CLEAN_GEOMETRY
136     
137           DO PASS=WRITE_HEADER,WRITE_DATA
138     
139     
140              CALL WRITE_GEOMETRY_IN_VTU_BIN(PASS)
141     
142              IF(VTK_EP_g(VTK_REGION)) &
143                 CALL WRITE_SCALAR_IN_VTU_BIN('EP_G',EP_G,PASS)
144     
145              IF(VTK_P_g(VTK_REGION)) &
146                 CALL WRITE_SCALAR_IN_VTU_BIN('P_G',P_G,PASS)
147     
148              IF(VTK_P_star(VTK_REGION)) &
149                 CALL WRITE_SCALAR_IN_VTU_BIN('P_S',P_S,PASS)
150     
151              IF(VTK_VEL_g(VTK_REGION)) &
152                 CALL WRITE_VECTOR_IN_VTU_BIN('Gas_Velocity',U_G,V_G,W_G,PASS)
153     
154              IF(VTK_U_g(VTK_REGION)) &
155                 CALL WRITE_SCALAR_IN_VTU_BIN('U_G',P_G,PASS)
156     
157              IF(VTK_V_g(VTK_REGION)) &
158                 CALL WRITE_SCALAR_IN_VTU_BIN('V_G',P_G,PASS)
159     
160              IF(VTK_W_g(VTK_REGION)) &
161                 CALL WRITE_SCALAR_IN_VTU_BIN('W_G',P_G,PASS)
162     
163              DO M = 1,MMAX
164                 IF(VTK_VEL_s(VTK_REGION,M)) THEN
165                    WRITE(SUBM,*)M
166                    CALL WRITE_VECTOR_IN_VTU_BIN('Solids_Velocity_'//ADJUSTL(SUBM),U_S(:,M),V_S(:,M),W_S(:,M),PASS)
167                 ENDIF
168                 IF(VTK_U_s(VTK_REGION,M)) THEN
169                    WRITE(SUBM,*)M
170                    CALL WRITE_SCALAR_IN_VTU_BIN('U_s_'//ADJUSTL(SUBM),U_S(:,M),PASS)
171                 ENDIF
172                 IF(VTK_V_s(VTK_REGION,M)) THEN
173                    WRITE(SUBM,*)M
174                    CALL WRITE_SCALAR_IN_VTU_BIN('V_s_'//ADJUSTL(SUBM),V_S(:,M),PASS)
175                 ENDIF
176                 IF(VTK_W_s(VTK_REGION,M)) THEN
177                    WRITE(SUBM,*)M
178                    CALL WRITE_SCALAR_IN_VTU_BIN('W_s_'//ADJUSTL(SUBM),W_S(:,M),PASS)
179                 ENDIF
180              END DO
181     
182              DO M = 1,MMAX
183                 IF(VTK_ROP_s(VTK_REGION,M)) THEN
184                    WRITE(SUBM,*)M
185                    CALL WRITE_SCALAR_IN_VTU_BIN('Solids_density_'//ADJUSTL(SUBM),ROP_S(:,M),PASS)
186                 ENDIF
187              END DO
188     
189              IF(VTK_T_g(VTK_REGION)) &
190                 CALL WRITE_SCALAR_IN_VTU_BIN('Gas_temperature',T_g,PASS)
191     
192              DO M = 1,MMAX
193                 IF(VTK_T_s(VTK_REGION,M)) THEN
194                    WRITE(SUBM,*)M
195                    CALL WRITE_SCALAR_IN_VTU_BIN('Solids_temperature_'//ADJUSTL(SUBM),T_S(:,M),PASS)
196                 ENDIF
197              END DO
198     
199     
200              SPECIES_COUNTER = 0
201              DO N = 1,NMAX(0)
202                 IF(VTK_X_g(VTK_REGION,N)) THEN
203                    WRITE(SUBN,*)N
204                    IF(USE_RRATES) THEN
205                       SPECIES_COUNTER = SPECIES_COUNTER + 1
206                       VAR_NAME = ADJUSTL(SPECIES_NAME(SPECIES_COUNTER))
207                       LT = LEN_TRIM(ADJUSTL(SPECIES_NAME(SPECIES_COUNTER)))
208                    ELSE
209                       VAR_NAME = ADJUSTL(SPECIES_ALIAS_g(N))
210                       LT = LEN_TRIM(ADJUSTL(SPECIES_ALIAS_g(N)))
211                    ENDIF
212                    VAR_NAME = VAR_NAME(1:LT)//'_Gas_mass_fractions_'//ADJUSTL(SUBN)
213                    CALL WRITE_SCALAR_IN_VTU_BIN(VAR_NAME,X_g(:,N),PASS)
214                 ENDIF
215              END DO
216     
217             DO M = 1, MMAX
218                WRITE(SUBM,*)M
219                DO N = 1,NMAX(M)
220                   IF(VTK_X_s(VTK_REGION,M,N)) THEN
221                      WRITE(SUBN,*)N
222                      IF(USE_RRATES) THEN
223                         SPECIES_COUNTER = SPECIES_COUNTER + 1
224                         VAR_NAME = ADJUSTL(SPECIES_NAME(SPECIES_COUNTER))
225                         LT = LEN_TRIM(ADJUSTL(SPECIES_NAME(SPECIES_COUNTER)))
226                      ELSE
227                         VAR_NAME = ADJUSTL(SPECIES_ALIAS_s(M,N))
228                         LT = LEN_TRIM(ADJUSTL(SPECIES_ALIAS_s(M,N)))
229                      ENDIF
230                      VAR_NAME = VAR_NAME(1:LT)//'_Solids_mass_fractions_'//TRIM(ADJUSTL(SUBM))//'_'//ADJUSTL(SUBN)
231                      CALL WRITE_SCALAR_IN_VTU_BIN(VAR_NAME,X_s(:,M,N),PASS)
232                   ENDIF
233                END DO
234             END DO
235     
236             DO M = 1,MMAX
237                IF(VTK_Theta_m(VTK_REGION,M)) THEN
238                   WRITE(SUBM,*)M
239                   CALL WRITE_SCALAR_IN_VTU_BIN('Granular_temperature_'//ADJUSTL(SUBM),Theta_m(:,M),PASS)
240                ENDIF
241             END DO
242     
243             DO N = 1,NSCALAR
244                IF(VTK_Scalar(VTK_REGION,N)) THEN
245                   WRITE(SUBN,*)N
246                   VAR_NAME = 'Scalar_'//ADJUSTL(SUBN)
247                   CALL WRITE_SCALAR_IN_VTU_BIN(VAR_NAME,Scalar(:,N),PASS)
248                ENDIF
249             END DO
250     
251             DO R = 1,nRR
252                IF(VTK_RRate(VTK_REGION,R)) THEN
253                   WRITE(SUBR,*)R
254                   VAR_NAME = 'RRates_'//ADJUSTL(SUBR)
255                   CALL WRITE_SCALAR_IN_VTU_BIN(VAR_NAME,ReactionRates(:, R),PASS)
256                ENDIF
257            END DO
258     
259            IF(K_EPSILON) THEN
260               IF(VTK_K_Turb_G(VTK_REGION)) &
261                  CALL WRITE_SCALAR_IN_VTU_BIN('K_Turb_G',K_Turb_G,PASS)
262               IF(VTK_E_Turb_G(VTK_REGION)) &
263                  CALL WRITE_SCALAR_IN_VTU_BIN('E_Turb_G',E_Turb_G,PASS)
264            ENDIF
265     
266     
267            IF(VTK_VORTICITY(VTK_REGION).OR.VTK_LAMBDA_2(VTK_REGION)) THEN
268               CALL CALC_VORTICITY
269            ENDIF
270     
271            IF(VTK_VORTICITY(VTK_REGION)) &
272               CALL WRITE_SCALAR_IN_VTU_BIN('VORTICITY_MAG',VORTICITY,PASS)
273            IF(VTK_LAMBDA_2(VTK_REGION)) &
274               CALL WRITE_SCALAR_IN_VTU_BIN('LAMBDA_2',LAMBDA2,PASS)
275     
276     
277            IF(VTK_PARTITION(VTK_REGION)) &
278               CALL WRITE_SCALAR_IN_VTU_BIN('PARTITION',PARTITION,PASS)
279     
280     
281            IF(VTK_BC_ID(VTK_REGION)) THEN
282               Allocate(DP_BC_ID(DIMENSION_3))
283               DP_BC_ID = DBLE(BC_ID)
284               CALL WRITE_SCALAR_IN_VTU_BIN('BC_ID',DP_BC_ID,PASS)
285               DeAllocate(DP_BC_ID)
286            ENDIF
287     
288     
289            IF(VTK_DWALL(VTK_REGION)) &
290               CALL WRITE_SCALAR_IN_VTU_BIN('DISTANCE_TO_WALL',DWALL,PASS)
291     
292            IF(DISCRETE_ELEMENT.AND.USE_STL) THEN
293               IF(VTK_FACET_COUNT_DES(VTK_REGION)) THEN
294                  ALLOCATE(FACET_COUNT_DES(DIMENSION_3))
295     
296                  DO IJK = IJKSTART3, IJKEND3
297                     FACET_COUNT_DES(IJK) = LIST_FACET_AT_DES(IJK)%COUNT_FACETS
298                  ENDDO
299     
300                  CALL WRITE_SCALAR_IN_VTU_BIN('FACET_COUNT', FACET_COUNT_DES,PASS)
301                  DEALLOCATE(FACET_COUNT_DES)
302               ENDIF
303     
304               IF(VTK_NB_FACET_DES(VTK_REGION)) THEN
305                  ALLOCATE(NEIGHBORING_FACET(DIMENSION_3))
306     
307                  DO IJK = IJKSTART3, IJKEND3
308                     NEIGHBORING_FACET(IJK) = 1.0
309                     IF(NO_NEIGHBORING_FACET_DES(IJK)) NEIGHBORING_FACET(IJK) = 0.0
310                  ENDDO
311     
312                  CALL WRITE_SCALAR_IN_VTU_BIN('NEIGH_FACET', NEIGHBORING_FACET,PASS)
313                  DEALLOCATE(NEIGHBORING_FACET)
314               ENDIF
315            ENDIF
316     
317     
318            IF(VTK_IJK(VTK_REGION)) THEN
319              Allocate(IJK_ARRAY(DIMENSION_3))
320              DO IJK = IJKSTART3, IJKEND3
321                 IJK_ARRAY(IJK) = DBLE(IJK)
322              ENDDO
323              CALL WRITE_SCALAR_IN_VTU_BIN('IJK',IJK_ARRAY,PASS)
324              DeAllocate(IJK_ARRAY)
325           ENDIF
326     
327            IF(VTK_IJK(VTK_REGION)) &
328               CALL WRITE_VECTOR_IN_VTU_BIN('Scalar normal',NORMAL_S(:,1),NORMAL_S(:,2),NORMAL_S(:,3),PASS)
329     
330            DO N=1,15
331               IF(VTK_DEBUG(VTK_REGION,N)) THEN
332                  WRITE(SUBN,*)N
333                  VAR_NAME = 'DEBUG_'//ADJUSTL(SUBN)
334                  CALL WRITE_SCALAR_IN_VTU_BIN(VAR_NAME,DEBUG_CG(:,N),PASS)
335               ENDIF
336            ENDDO
337     
338     
339     
340           ENDDO ! PASS LOOP, EITHER HEADER OR DATA
341     
342     
343           CALL CLOSE_VTU_FILE_BIN(MODE)
344           IF(MODE==0) CALL UPDATE_AND_CLOSE_PVD_FILE
345     
346           call MPI_barrier(MPI_COMM_WORLD,mpierr)
347     
348     ! Update Frames
349           IF (myPE == PE_IO.AND.TIME_DEPENDENT_FILENAME) THEN
350              OPEN(CONVERT='BIG_ENDIAN',UNIT = VTU_FRAME_UNIT, FILE = TRIM(VTU_FRAME_FILENAME))
351              DO L = 1, DIMENSION_VTK
352                 IF(VTK_DEFINED(L)) WRITE(VTU_FRAME_UNIT,*) L,FRAME(L)
353              ENDDO
354              CLOSE(VTU_FRAME_UNIT)
355           ENDIF
356     
357          IF (FULL_LOG.AND.myPE == PE_IO) WRITE(*,20)' DONE.'
358     
359     20    FORMAT(A,1X/)
360           RETURN
361     
362           END SUBROUTINE WRITE_VTU_FILE
363     
364     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
365     !                                                                      C
366     !  Module name: OPEN_VTU_FILE                                          C
367     !  Purpose: Open a vtu file and writes the header                      C
368     !           Binary format                                              C
369     !                                                                      C
370     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
371     !  Reviewer:                                          Date:            C
372     !                                                                      C
373     !  Revision Number #                                  Date: ##-###-##  C
374     !  Author: #                                                           C
375     !  Purpose: #                                                          C
376     !                                                                      C
377     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
378       SUBROUTINE OPEN_VTU_FILE_BIN(MODE)
379     
380           USE compar
381           USE constant
382           USE cutcell
383           USE fldvar
384           USE functions
385           USE geometry
386           USE indices
387           USE mpi_utility
388           USE output
389           USE parallel
390           USE param
391           USE param1
392           USE quadric
393           USE run
394           USE sendrecv
395           USE toleranc
396           USE vtk
397           use cdist
398     
399           IMPLICIT NONE
400           LOGICAL :: VTU_FRAME_FILE_EXISTS
401           INTEGER :: ISTAT,BUFF1,BUFF2,L
402           INTEGER :: MODE   ! MODE = 0 : Write regular VTK region file
403                             ! MODE = 1 : Write debug   VTK region file (VTK_DBG_FILE = .TRUE.)
404     
405     
406           call MPI_barrier(MPI_COMM_WORLD,mpierr)
407     
408     ! Only open the file from head node when not using distributed I/O
409           IF (myPE /= PE_IO.AND.(.NOT.BDIST_IO)) RETURN
410     
411           IF(TIME_DEPENDENT_FILENAME) THEN
412              INQUIRE(FILE=VTU_FRAME_FILENAME,EXIST=VTU_FRAME_FILE_EXISTS)
413              IF(VTU_FRAME_FILE_EXISTS) THEN
414                 OPEN(CONVERT='BIG_ENDIAN',UNIT = VTU_FRAME_UNIT, FILE = TRIM(VTU_FRAME_FILENAME))
415                 DO L = 1, DIMENSION_VTK
416                    IF(VTK_DEFINED(L)) THEN
417                       READ(VTU_FRAME_UNIT,*)BUFF1,BUFF2
418                       FRAME(L)=BUFF2
419                    ENDIF
420                 ENDDO
421                 CLOSE(VTU_FRAME_UNIT)
422              ENDIF
423              IF(RESET_FRAME_AT_TIME_ZERO.AND.TIME==ZERO) THEN
424                 DO L = 1, DIMENSION_VTK
425                    IF(L==VTK_REGION) FRAME(L)=-1
426                 ENDDO
427              ENDIF
428              DO L = 1, DIMENSION_VTK
429                 IF(L==VTK_REGION) FRAME(L) = FRAME(L) + 1
430              ENDDO
431           ENDIF
432     
433           IF (BDIST_IO.AND.NUMBER_OF_VTK_CELLS>0) THEN
434     
435     
436     ! For distributed I/O, define the file name for each processor
437              IF(TIME_DEPENDENT_FILENAME.AND.MODE==0) THEN
438                 WRITE(VTU_FILENAME,20) TRIM(VTK_FILEBASE(VTK_REGION)),FRAME(VTK_REGION),MYPE
439              ELSE
440                 WRITE(VTU_FILENAME,25) TRIM(VTK_FILEBASE(VTK_REGION)),MYPE
441              ENDIF
442           ELSE
443              IF(MYPE.EQ.PE_IO) THEN
444                 IF(TIME_DEPENDENT_FILENAME.AND.MODE==0) THEN
445                    WRITE(VTU_FILENAME,30) TRIM(VTK_FILEBASE(VTK_REGION)),FRAME(VTK_REGION)
446                 ELSE
447                    WRITE(VTU_FILENAME,35) TRIM(VTK_FILEBASE(VTK_REGION))
448                 ENDIF
449              END IF
450           END IF
451     
452     ! Add the VTU directory path if necessary
453     
454           IF(TRIM(VTU_DIR)/='.') VTU_FILENAME='./'//TRIM(VTU_DIR)//'/'//VTU_FILENAME
455     
456     ! Echo
457           IF (FULL_LOG) THEN
458              IF (.NOT.BDIST_IO) THEN
459                 WRITE(*,10,ADVANCE='NO')' WRITING VTU FILE : ', TRIM(VTU_FILENAME),' .'
460              ELSE
461                 IF(myPE==PE_IO) WRITE(*,15,ADVANCE='NO')' EACH PROCESOR IS WRITING ITS OWN VTU FILE.'
462              ENDIF
463           ENDIF
464     
465     ! Open File
466     !      OPEN(CONVERT='BIG_ENDIAN',UNIT = VTU_UNIT, FILE = TRIM(VTU_FILENAME),FORM='BINARY',IOSTAT=ISTAT)
467     
468     
469           IF(NUMBER_OF_VTK_CELLS>0) THEN
470     
471              OPEN(CONVERT='BIG_ENDIAN',UNIT     = VTU_UNIT,           &
472                   FILE     = TRIM(VTU_FILENAME), &
473                   FORM     = 'UNFORMATTED',      &  ! works with gfortran 4.3.4 and ifort 10.1 but may not be supported by all compilers
474                                                     ! use 'BINARY' if 'UNFORMATTED' is not supported
475                   ACCESS   = 'STREAM',           &  ! works with gfortran 4.3.4 and ifort 10.1 but may not be supported by all compilers
476                                                     ! use 'SEQUENTIAL' if 'STREAM' is not supported
477                   ACTION   = 'WRITE',            &
478                   IOSTAT=ISTAT)
479     
480     
481              IF (ISTAT /= 0) THEN
482                 IF(DMP_LOG) WRITE(UNIT_LOG, 1001) VTU_FILENAME, VTU_UNIT,VTU_DIR
483                 IF(FULL_LOG.AND.myPE == PE_IO) WRITE(*, 1001) VTU_FILENAME, VTU_UNIT,VTU_DIR
484                 CALL MFIX_EXIT(myPE)
485              ENDIF
486     
487     
488         1001 FORMAT(/1X,70('*')//, ' From: OPEN_VTU_FILE',/,' Message: ',          &
489                 'Error opening vtu file. Terminating run.',/10X,          &
490                 'File name:  ',A,/10X,                                         &
491                 'DES_UNIT :  ',i4, /10X,                                       &
492                 'PLEASE VERIFY THAT VTU_DIR EXISTS: ', A, &
493                 /1X,70('*')/)
494     
495     
496        ! Write file Header
497              BUFFER='<?xml version="1.0"?>'
498              WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
499     
500     
501              WRITE(BUFFER,110)'<!-- Time =',TIME,' sec. -->'
502              WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
503     
504              BUFFER='<VTKFile type="UnstructuredGrid" version="0.1" byte_order="BigEndian">'
505              WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
506     
507              BUFFER='  <UnstructuredGrid>'
508              WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
509     
510           ENDIF
511     ! For distributed I/O, open .pvtu file that combines all *.vtu files for a given FRAME
512     ! this is a simple ASCII file
513     
514           IF (myPE == PE_IO.AND.BDIST_IO) THEN
515     
516              IF(TIME_DEPENDENT_FILENAME.AND.MODE==0) THEN
517                 WRITE(PVTU_FILENAME,40) TRIM(VTK_FILEBASE(VTK_REGION)),FRAME(VTK_REGION)
518              ELSE
519                 WRITE(PVTU_FILENAME,45) TRIM(VTK_FILEBASE(VTK_REGION))
520              ENDIF
521     
522              IF(TRIM(VTU_DIR)/='.') PVTU_FILENAME='./'//TRIM(VTU_DIR)//'/'//PVTU_FILENAME
523     
524              OPEN(CONVERT='BIG_ENDIAN',UNIT = PVTU_UNIT, FILE = TRIM(PVTU_FILENAME))
525     
526              WRITE(PVTU_UNIT,100) '<?xml version="1.0"?>'
527              WRITE(PVTU_UNIT,110) '<!-- Time =',TIME,' sec. -->'
528              WRITE(PVTU_UNIT,120) '<VTKFile type="PUnstructuredGrid"',&
529                       ' version="0.1" byte_order="BigEndian">'
530     
531              WRITE(PVTU_UNIT,100) '  <PUnstructuredGrid GhostLevel="0">'
532              WRITE(PVTU_UNIT,100) '      <PPoints>'
533              WRITE(PVTU_UNIT,100) '        <PDataArray type="Float32" Name="coordinates" NumberOfComponents="3" &
534                   &format="appended" offset=" 0" />'
535              WRITE(PVTU_UNIT,100) '      </PPoints>'
536              WRITE(PVTU_UNIT,100) ''
537              WRITE(PVTU_UNIT,100) '      <PCellData Scalars="scalars">'
538     
539           ENDIF
540     
541     100   FORMAT(A)
542     110   FORMAT(A,E14.7,A)
543     120   FORMAT(A,A)
544     10    FORMAT(/1X,3A)
545     15    FORMAT(/1X,A)
546     20    FORMAT(A,"_",I4.4,"_",I5.5,".vtu")
547     25    FORMAT(A,"_",I5.5,".vtu")
548     30    FORMAT(A,"_",I4.4,".vtu")
549     35    FORMAT(A,".vtu")
550     40    FORMAT(A,"_",I4.4,".pvtu")
551     45    FORMAT(A,".pvtu")
552     
553           RETURN
554     
555           END SUBROUTINE OPEN_VTU_FILE_BIN
556     
557     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
558     !                                                                      C
559     !  Module name: WRITE_GEOMETRY_IN_VTU_BIN                              C
560     !  Purpose: Write Geometry and connectivity in a vtu file              C
561     !           Binary format                                              C
562     !                                                                      C
563     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
564     !  Reviewer:                                          Date:            C
565     !                                                                      C
566     !  Revision Number #                                  Date: ##-###-##  C
567     !  Author: #                                                           C
568     !  Purpose: #                                                          C
569     !                                                                      C
570     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
571       SUBROUTINE WRITE_GEOMETRY_IN_VTU_BIN(PASS)
572     
573           USE, INTRINSIC :: iso_c_binding
574           USE param
575           USE param1
576           USE parallel
577           USE constant
578           USE run
579           USE toleranc
580           USE geometry
581           USE indices
582           USE compar
583           USE mpi_utility
584           USE sendrecv
585           USE quadric
586           USE cutcell
587           USE fldvar
588           USE vtk
589           USE cdist
590           USE functions
591     
592           IMPLICIT NONE
593     
594           INTEGER :: IJK,L
595           INTEGER :: OFFSET
596     
597           INTEGER :: CELL_TYPE
598     
599           REAL(c_float) :: float
600           INTEGER(c_int) :: int
601     
602           INTEGER ::     nbytes_xyz,nbytes_connectivity,nbytes_offset,nbytes_type
603           INTEGER ::     offset_xyz,offset_connectivity,offset_offset,offset_type
604     
605           INTEGER :: PASS
606           INTEGER :: WRITE_HEADER = 1
607           INTEGER :: WRITE_DATA   = 2
608     
609     
610     ! First a series of tags is written for the geometry (PASS=WRITE_HEADER)
611     !  - Coordinates
612     !  - Connectivity
613     !  - Connectivity offset
614     !  - cell types
615     !
616     
617     ! Since the data is appended (i.e., written after all tags), the
618     ! offset, in number of bytes must be specified.  The offset includes
619     ! the size of the data for each field, plus the size of the integer
620     ! that stores the number of bytes.  this is why the offset of a field
621     ! equals the offset of the previous field plus sizeof(int) plus the
622     ! number of bytes of the field.
623     
624     ! Next, the actual data is written for the geometry (PASS=WRITE_DATA)
625     ! The DATA is converted to single precision to save memory.
626     
627           IF (myPE == PE_IO.AND.(.NOT.BDIST_IO)) THEN
628     ! The number of points and number of VTK cells is now computed in
629     ! SETUP_VTK_REGION
630     
631     ! Number of bytes of each field
632              nbytes_xyz          = NUMBER_OF_POINTS * 3 * sizeof(float)
633     
634              nbytes_connectivity = 0
635              DO IJK = 1,IJKMAX3
636                 IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) THEN
637                       nbytes_connectivity = nbytes_connectivity + GLOBAL_NUMBER_OF_NODES(IJK)
638                 ENDIF
639              END DO
640              nbytes_connectivity = nbytes_connectivity * sizeof(int)
641     
642     
643              nbytes_offset       = NUMBER_OF_VTK_CELLS * sizeof(int)
644     
645              nbytes_type         = NUMBER_OF_VTK_CELLS * sizeof(int)
646     
647     
648     ! Offset of each field
649              offset_xyz = 0
650              offset_connectivity = offset_xyz          + sizeof(int) + nbytes_xyz
651              offset_offset       = offset_connectivity + sizeof(int) + nbytes_connectivity
652              offset_type         = offset_offset       + sizeof(int) + nbytes_offset
653     
654     
655              IF(PASS==WRITE_HEADER) THEN
656     
657                 WRITE(BUFFER,100)'    <Piece NumberOfPoints="',NUMBER_OF_POINTS,'" NumberOfCells="',NUMBER_OF_VTK_CELLS,'" >'
658                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
659     
660                 WRITE(BUFFER,110)'      <Points>'
661                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
662     
663                 WRITE(BUFFER,100)'        <DataArray type="Float32" Name="coordinates" NumberOfComponents="3" &
664                      &format="appended" offset="',offset_xyz,'" />'
665                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
666     
667                 WRITE(BUFFER,110)'      </Points>'
668                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
669     
670                 WRITE(BUFFER,110)'      <Cells>'
671                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
672     
673                 WRITE(BUFFER,100)'        <DataArray type="Int32" Name="connectivity" format="appended" offset="', &
674                      offset_connectivity,'" />'
675                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
676     
677                 WRITE(BUFFER,100)'        <DataArray type="Int32" Name="offsets" format="appended" offset="',offset_offset,'" />'
678                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
679     
680                 WRITE(BUFFER,100)'        <DataArray type="Int32" Name="types" format="appended" offset="',offset_type,'" />'
681                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
682     
683                 WRITE(BUFFER,110)'      </Cells>'
684                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
685     
686                 VTU_offset =  offset_type       + sizeof(int) + nbytes_type  ! Store offset for first variable to be written
687     
688                 WRITE(BUFFER,110)'      <CellData>'                          ! Preparing CellData
689                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
690     
691     
692     
693     
694              ELSEIF(PASS==WRITE_DATA) THEN
695     
696                 WRITE(BUFFER,110)'      </CellData>'
697                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
698     
699                 WRITE(BUFFER,110)'    </Piece>'
700                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
701     
702     
703                 WRITE(BUFFER,110)'  </UnstructuredGrid>'
704                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
705     
706     
707                 WRITE(BUFFER,110)'  <AppendedData encoding="raw">'
708                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
709     
710     
711     ! Starting raw binary data with an underscore
712     
713                 WRITE(BUFFER,110)'_'
714                 WRITE(VTU_UNIT)TRIM(BUFFER)
715     
716     ! X,Y,Z coordinates
717                 WRITE(VTU_UNIT) nbytes_xyz, (GLOBAL_COORDS_OF_POINTS(1:3,L), L = 1,NUMBER_OF_POINTS)
718     
719     ! Connectivity
720                 WRITE(VTU_UNIT) nbytes_connectivity
721     
722                 DO IJK = 1,IJKMAX3
723                    IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) THEN
724                       WRITE(VTU_UNIT) (GLOBAL_CLEANED_CONNECTIVITY(IJK,L)-1,L=1,GLOBAL_NUMBER_OF_NODES(IJK))
725                    ENDIF
726                 END DO
727     
728     ! Offsets
729                 WRITE(VTU_UNIT) nbytes_offset
730     
731                 OFFSET = 0
732                 DO IJK = 1,IJKMAX3
733                    IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) THEN
734                       OFFSET = OFFSET + GLOBAL_NUMBER_OF_NODES(IJK)
735                       WRITE(VTU_UNIT) OFFSET
736                    ENDIF
737                 END DO
738     
739     ! Types
740                 WRITE(VTU_UNIT)nbytes_type
741     
742                 IF(NO_K) THEN
743                    CELL_TYPE = 7
744                 ELSE
745                    CELL_TYPE = 41
746                 ENDIF
747     
748                 DO IJK = 1,IJKMAX3
749                    IF (BELONGS_TO_VTK_SUBDOMAIN(IJK))  WRITE(VTU_UNIT) CELL_TYPE
750                 END DO
751     
752     
753              ENDIF
754     
755     
756           ELSEIF(BDIST_IO) THEN
757     
758     ! For distributed I/O, it works the same as above, except, the data is local to each processor
759     ! First compute local number of cells and points
760     
761     ! The number of points and number of VTK cells is now computed in
762     ! SETUP_VTK_REGION
763     
764     ! Number of bytes of each field
765              nbytes_xyz          = NUMBER_OF_POINTS * 3 * sizeof(float)
766     
767              nbytes_connectivity = 0
768              DO IJK = 1,IJKEND3
769                 IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) THEN
770                       nbytes_connectivity = nbytes_connectivity + NUMBER_OF_NODES(IJK)
771                 ENDIF
772              END DO
773              nbytes_connectivity = nbytes_connectivity * sizeof(int)
774     
775     
776              nbytes_offset       = NUMBER_OF_VTK_CELLS * sizeof(int)
777     
778              nbytes_type         = NUMBER_OF_VTK_CELLS * sizeof(int)
779     
780     
781     ! Offset of each field
782              offset_xyz = 0
783              offset_connectivity = offset_xyz          + sizeof(int) + nbytes_xyz
784              offset_offset       = offset_connectivity + sizeof(int) + nbytes_connectivity
785              offset_type         = offset_offset       + sizeof(int) + nbytes_offset
786     
787     
788              IF(PASS==WRITE_HEADER) THEN
789     
790                 WRITE(BUFFER,100)'    <Piece NumberOfPoints="',NUMBER_OF_POINTS,'" NumberOfCells="',NUMBER_OF_VTK_CELLS,'" >'
791                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
792     
793                 WRITE(BUFFER,110)'      <Points>'
794                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
795     
796                 WRITE(BUFFER,100)'        <DataArray type="Float32" Name="coordinates" NumberOfComponents="3" &
797                      &format="appended" offset="',offset_xyz,'" />'
798                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
799     
800                 WRITE(BUFFER,110)'      </Points>'
801                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
802     
803                 WRITE(BUFFER,110)'      <Cells>'
804                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
805     
806                 WRITE(BUFFER,100)'        <DataArray type="Int32" Name="connectivity" format="appended" offset="', &
807                      offset_connectivity,'" />'
808                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
809     
810                 WRITE(BUFFER,100)'        <DataArray type="Int32" Name="offsets" format="appended" offset="',offset_offset,'" />'
811                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
812     
813                 WRITE(BUFFER,100)'        <DataArray type="Int32" Name="types" format="appended" offset="',offset_type,'" />'
814                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
815     
816                 WRITE(BUFFER,110)'      </Cells>'
817                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
818     
819                 VTU_offset =  offset_type       + sizeof(int) + nbytes_type  ! Store offset for first variable to be written
820     
821                 WRITE(BUFFER,110)'      <CellData>'                          ! Preparing CellData
822                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
823     
824     
825     
826     
827              ELSEIF(PASS==WRITE_DATA) THEN
828     
829                 WRITE(BUFFER,110)'      </CellData>'
830                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
831     
832                 WRITE(BUFFER,110)'    </Piece>'
833                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
834     
835     
836                 WRITE(BUFFER,110)'  </UnstructuredGrid>'
837                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
838     
839     
840                 WRITE(BUFFER,110)'  <AppendedData encoding="raw">'
841                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
842     
843     
844     ! Starting raw binary data with an underscore
845     
846                 WRITE(BUFFER,110)'_'
847                 WRITE(VTU_UNIT)TRIM(BUFFER)
848     
849     ! X,Y,Z coordinates
850                 WRITE(VTU_UNIT) nbytes_xyz, (COORDS_OF_POINTS(L,1:3), L = 1,NUMBER_OF_POINTS)
851     
852     ! Connectivity
853                 WRITE(VTU_UNIT) nbytes_connectivity
854     
855                 DO IJK = 1,IJKEND3
856                    IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) THEN
857                       WRITE(VTU_UNIT) (CLEANED_CONNECTIVITY(IJK,L)-1,L=1,NUMBER_OF_NODES(IJK))
858                    ENDIF
859                 END DO
860     
861     ! Offsets
862                 WRITE(VTU_UNIT) nbytes_offset
863     
864                 OFFSET = 0
865                 DO IJK = 1,IJKEND3
866                    IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) THEN
867                       OFFSET = OFFSET + NUMBER_OF_NODES(IJK)
868                       WRITE(VTU_UNIT) OFFSET
869                    ENDIF
870                 END DO
871     
872     ! Types
873                 WRITE(VTU_UNIT)nbytes_type
874     
875                 IF(NO_K) THEN
876                    CELL_TYPE = 7
877                 ELSE
878                    CELL_TYPE = 41
879                 ENDIF
880     
881                 DO IJK = 1,IJKEND3
882                    IF (BELONGS_TO_VTK_SUBDOMAIN(IJK))  WRITE(VTU_UNIT) CELL_TYPE
883                 END DO
884     
885     
886              ENDIF
887     
888     
889           ENDIF
890     
891     
892     100   FORMAT(A,I12,A,I12,A)
893     110   FORMAT(A)
894     
895           RETURN
896     
897           END SUBROUTINE WRITE_GEOMETRY_IN_VTU_BIN
898     
899     
900     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
901     !                                                                      C
902     !  Module name: WRITE_SCALAR_IN_VTU_BIN                                C
903     !  Purpose: Write Scalar variable in a vtu file                        C
904     !           Binary format                                              C
905     !                                                                      C
906     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
907     !  Reviewer:                                          Date:            C
908     !                                                                      C
909     !  Revision Number #                                  Date: ##-###-##  C
910     !  Author: #                                                           C
911     !  Purpose: #                                                          C
912     !                                                                      C
913     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
914       SUBROUTINE WRITE_SCALAR_IN_VTU_BIN(VAR_NAME,VAR,PASS)
915     
916           USE, INTRINSIC :: iso_c_binding
917           USE param
918           USE param1
919           USE parallel
920           USE constant
921           USE run
922           USE toleranc
923           USE geometry
924           USE indices
925           USE compar
926           USE mpi_utility
927           USE sendrecv
928           USE quadric
929           USE cutcell
930           USE fldvar
931           USE vtk
932           USE cdist
933           USE output
934           USE functions
935     
936           IMPLICIT NONE
937           INTEGER :: I,IJK
938     
939           CHARACTER (*) :: VAR_NAME
940           DOUBLE PRECISION, DIMENSION(DIMENSION_3) ::  VAR
941           DOUBLE PRECISION, ALLOCATABLE :: GLOBAL_VAR(:)
942           DOUBLE PRECISION, DIMENSION(DIMENSION_3) ::  TMP_VAR
943     
944           REAL(c_float) :: float
945     
946           INTEGER :: nbytes_scalar
947     
948           INTEGER :: PASS
949           INTEGER :: WRITE_HEADER = 1
950           INTEGER :: WRITE_DATA   = 2
951     
952           IF (.NOT.BDIST_IO) THEN
953     
954     ! For each scalar, write a tag, with corresponding offset
955     
956              nbytes_scalar = NUMBER_OF_VTK_CELLS * sizeof(float)
957     
958              IF(PASS==WRITE_HEADER) THEN
959     !           For each scalar, write a tag, with corresponding offset
960     
961                 DO I = 1,LEN_TRIM(VAR_NAME)
962                    IF(VAR_NAME(I:I) == ' ') VAR_NAME(I:I) = '_'
963                 ENDDO
964     
965                 WRITE(BUFFER,90)'        <DataArray type="Float32" Name="', &
966                      TRIM(VAR_NAME),'" format="appended" offset="',VTU_offset,'" />'
967                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
968     
969                 VTU_offset = VTU_offset + sizeof(float) + nbytes_scalar
970     
971     
972              ELSEIF(PASS==WRITE_DATA) THEN
973     !           and write the data, always preceded by its size in number of bytes
974     
975                 IF (myPE == PE_IO) THEN
976                    allocate (GLOBAL_VAR(ijkmax3))
977                 ELSE
978                    allocate (GLOBAL_VAR(1))
979                 ENDIF
980     
981                 IF(RE_INDEXING) THEN
982                    CALL UNSHIFT_DP_ARRAY(VAR,TMP_VAR)
983                    CALL gather (TMP_VAR,GLOBAL_VAR,root)
984                 ELSE
985                    CALL gather (VAR,GLOBAL_VAR,root)
986                 ENDIF
987     
988                 IF (myPE /= PE_IO) RETURN
989     
990     
991                 WRITE(VTU_UNIT) nbytes_scalar
992     
993                 DO IJK = 1,IJKMAX3
994                    IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) WRITE(VTU_UNIT) REAL(GLOBAL_VAR(IJK))
995                 ENDDO
996     
997     
998                 Deallocate (GLOBAL_VAR)
999     
1000              ENDIF
1001     
1002     
1003           ELSE ! BDIST_IO=.TRUE.
1004     
1005     
1006              nbytes_scalar = NUMBER_OF_VTK_CELLS * sizeof(float)
1007     
1008              IF(PASS==WRITE_HEADER) THEN
1009     !           For each scalar, write a tag, with corresponding offset
1010     
1011                 DO I = 1,LEN_TRIM(VAR_NAME)
1012                    IF(VAR_NAME(I:I) == ' ') VAR_NAME(I:I) = '_'
1013                 ENDDO
1014     
1015                 WRITE(BUFFER,90)'        <DataArray type="Float32" Name="', &
1016                      TRIM(VAR_NAME),'" format="appended" offset="',VTU_offset,'" />'
1017                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1018     
1019     
1020                 VTU_offset = VTU_offset + sizeof(float) + nbytes_scalar
1021     
1022     
1023              ELSEIF(PASS==WRITE_DATA) THEN
1024     !           and write the data, always preceded by its size in number of bytes
1025     
1026                 WRITE(VTU_UNIT) nbytes_scalar
1027     
1028                 DO IJK = 1,IJKEND3
1029                    IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) WRITE(VTU_UNIT) REAL(VAR(IJK))
1030                 ENDDO
1031     
1032              ENDIF
1033     
1034     
1035              IF (myPE == PE_IO) THEN       ! Update pvtu file with variable name
1036                 WRITE(PVTU_UNIT,90) '        <DataArray type="Float32" Name="', &
1037                      TRIM(VAR_NAME),'" format="appended" offset="',VTU_offset,'" />'
1038              ENDIF
1039     
1040     
1041           ENDIF
1042     
1043     
1044           IF (PASS==WRITE_DATA.AND.FULL_LOG.AND.myPE == PE_IO) WRITE(*,10,ADVANCE='NO')'.'
1045     
1046     10    FORMAT(A)
1047     90    FORMAT(A,A,A,I12,A)
1048     
1049           RETURN
1050     
1051           END SUBROUTINE WRITE_SCALAR_IN_VTU_BIN
1052     
1053     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1054     !                                                                      C
1055     !  Module name: WRITE_VECTOR_IN_VTU                                    C
1056     !  Purpose: Write Vector variable in a vtu file                        C
1057     !           Binary format                                              C
1058     !                                                                      C
1059     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1060     !  Reviewer:                                          Date:            C
1061     !                                                                      C
1062     !  Revision Number #                                  Date: ##-###-##  C
1063     !  Author: #                                                           C
1064     !  Purpose: #                                                          C
1065     !                                                                      C
1066     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1067       SUBROUTINE WRITE_VECTOR_IN_VTU_BIN(VAR_NAME,VARX,VARY,VARZ,PASS)
1068     
1069           USE, INTRINSIC :: iso_c_binding
1070           USE param
1071           USE param1
1072           USE parallel
1073           USE constant
1074           USE run
1075           USE toleranc
1076           USE geometry
1077           USE indices
1078           USE compar
1079           USE mpi_utility
1080           USE sendrecv
1081           USE quadric
1082           USE cutcell
1083           USE fldvar
1084           USE vtk
1085           USE cdist
1086           USE output
1087           USE functions
1088     
1089           IMPLICIT NONE
1090           INTEGER :: IJK
1091     
1092           CHARACTER (*) :: VAR_NAME
1093           DOUBLE PRECISION, DIMENSION(DIMENSION_3) ::  VARX,VARY,VARZ
1094           DOUBLE PRECISION, ALLOCATABLE :: GLOBAL_VARX(:),GLOBAL_VARY(:),GLOBAL_VARZ(:)
1095           DOUBLE PRECISION, DIMENSION(DIMENSION_3) ::  TMP_VAR
1096     
1097           REAL(c_float) :: float
1098     
1099           INTEGER :: nbytes_vector
1100     
1101           INTEGER :: PASS
1102           INTEGER :: WRITE_HEADER = 1
1103           INTEGER :: WRITE_DATA   = 2
1104     
1105           IF (.NOT.BDIST_IO) THEN
1106     
1107              nbytes_vector = NUMBER_OF_VTK_CELLS * 3 * sizeof(float)
1108     
1109              IF(PASS==WRITE_HEADER) THEN
1110     !           For each vector, write a tag, with corresponding offset
1111     
1112                 WRITE(BUFFER,90)'        <DataArray type="Float32" Name="', &
1113                      TRIM(VAR_NAME),'"  NumberOfComponents="3" format="appended" offset="',VTU_offset,'" />'
1114                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1115     
1116     
1117                 VTU_offset = VTU_offset + sizeof(float) + nbytes_vector
1118     
1119     
1120              ELSEIF(PASS==WRITE_DATA) THEN
1121     !           and write the data, always preceded by its size in number of bytes
1122     
1123                 IF (myPE == PE_IO) THEN
1124                    allocate (GLOBAL_VARX(ijkmax3))
1125                    allocate (GLOBAL_VARY(ijkmax3))
1126                    allocate (GLOBAL_VARZ(ijkmax3))
1127                 ELSE
1128                    allocate (GLOBAL_VARX(1))
1129                    allocate (GLOBAL_VARY(1))
1130                    allocate (GLOBAL_VARZ(1))
1131                 ENDIF
1132     
1133                 IF(RE_INDEXING) THEN
1134                    CALL UNSHIFT_DP_ARRAY(VARX,TMP_VAR)
1135                    call gather (TMP_VAR,GLOBAL_VARX,root)
1136     
1137                    CALL UNSHIFT_DP_ARRAY(VARY,TMP_VAR)
1138                    call gather (TMP_VAR,GLOBAL_VARY,root)
1139     
1140                    CALL UNSHIFT_DP_ARRAY(VARZ,TMP_VAR)
1141                    call gather (TMP_VAR,GLOBAL_VARZ,root)
1142     
1143                 ELSE
1144                    call gather (VARX,GLOBAL_VARX,root)
1145                    call gather (VARY,GLOBAL_VARY,root)
1146                    call gather (VARZ,GLOBAL_VARZ,root)
1147                 ENDIF
1148     
1149     
1150                 IF (myPE /= PE_IO) RETURN
1151     
1152     
1153                 WRITE(VTU_UNIT) nbytes_vector
1154     
1155                 DO IJK = 1,IJKMAX3
1156                    IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) THEN
1157                       WRITE(VTU_UNIT) REAL(GLOBAL_VARX(IJK)),REAL(GLOBAL_VARY(IJK)),REAL(GLOBAL_VARZ(IJK))
1158                    ENDIF
1159                 ENDDO
1160     
1161     
1162                 Deallocate (GLOBAL_VARX)
1163                 Deallocate (GLOBAL_VARY)
1164                 Deallocate (GLOBAL_VARZ)
1165     
1166              ENDIF
1167     
1168     
1169           ELSE ! BDIST_IO=.TRUE.
1170     
1171     
1172              nbytes_vector = NUMBER_OF_VTK_CELLS * 3 * sizeof(float)
1173     
1174              IF(PASS==WRITE_HEADER) THEN
1175     !           For each vector, write a tag, with corresponding offset
1176     
1177     
1178                 WRITE(BUFFER,90)'        <DataArray type="Float32" Name="', &
1179                      TRIM(VAR_NAME),'"  NumberOfComponents="3" format="appended" offset="',VTU_offset,'" />'
1180                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1181     
1182     
1183                 VTU_offset = VTU_offset + sizeof(float) + nbytes_vector
1184     
1185     
1186              ELSEIF(PASS==WRITE_DATA) THEN
1187     !           and write the data, always preceded by its size in number of bytes
1188     
1189                 WRITE(VTU_UNIT) nbytes_vector
1190     
1191                 DO IJK = 1,IJKEND3
1192                    IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) THEN
1193                       WRITE(VTU_UNIT) REAL(VARX(IJK)),REAL(VARY(IJK)),REAL(VARZ(IJK))
1194                    ENDIF
1195                 ENDDO
1196     
1197              ENDIF
1198     
1199     
1200              IF (myPE == PE_IO) THEN       ! Update pvtu file with variable name
1201                 WRITE(PVTU_UNIT,90)'        <DataArray type="Float32" Name="', &
1202                      TRIM(VAR_NAME),'"  NumberOfComponents="3" format="appended" offset="',VTU_offset,'" />'
1203              ENDIF
1204     
1205           ENDIF
1206     
1207     
1208           IF (PASS==WRITE_DATA.AND.FULL_LOG.AND.myPE == PE_IO) WRITE(*,10,ADVANCE='NO')'.'
1209     
1210     10    FORMAT(A)
1211     90    FORMAT(A,A,A,I12,A)
1212     
1213           RETURN
1214     
1215           END SUBROUTINE WRITE_VECTOR_IN_VTU_BIN
1216     
1217     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1218     !                                                                      C
1219     !  Module name: CLOSE_VTU_FILE_BIN                                     C
1220     !  Purpose: Close a vtu file                                           C
1221     !           Binary format                                              C
1222     !                                                                      C
1223     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1224     !  Reviewer:                                          Date:            C
1225     !                                                                      C
1226     !  Revision Number #                                  Date: ##-###-##  C
1227     !  Author: #                                                           C
1228     !  Purpose: #                                                          C
1229     !                                                                      C
1230     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1231       SUBROUTINE CLOSE_VTU_FILE_BIN(MODE)
1232     
1233           USE compar
1234           Use run
1235           USE vtk
1236           use cdist
1237           USE mpi_utility
1238     
1239           IMPLICIT NONE
1240     
1241           INTEGER:: N
1242           CHARACTER (LEN=32)  :: VTU_NAME
1243           INTEGER, DIMENSION(0:numPEs-1) :: ALL_VTK_CELL_COUNT
1244           INTEGER :: IERR
1245           INTEGER :: MODE   ! MODE = 0 : Write regular VTK region file
1246                             ! MODE = 1 : Write debug   VTK region file (VTK_DBG_FILE = .TRUE.)
1247     
1248     
1249           IF (myPE /= PE_IO.AND.(.NOT.BDIST_IO)) RETURN
1250     
1251           IF(NUMBER_OF_VTK_CELLS>0) THEN
1252     
1253     ! Write last tags and close the vtu file
1254               WRITE(BUFFER,110)'  </AppendedData>'
1255               WRITE(VTU_UNIT)END_REC//TRIM(BUFFER)//END_REC
1256     
1257               WRITE(BUFFER,110)'</VTKFile>'
1258               WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1259     
1260              CLOSE(VTU_UNIT)
1261     
1262           ENDIF
1263     
1264     ! Update pvtu file and close
1265     
1266           IF(BDIST_IO)  CALL allgather_1i (NUMBER_OF_VTK_CELLS,ALL_VTK_CELL_COUNT,IERR)
1267     
1268           IF (myPE == PE_IO.AND.BDIST_IO) THEN
1269              WRITE(PVTU_UNIT,100) '      </PCellData>'
1270     
1271              DO N = 0,NumPEs-1
1272                 IF(ALL_VTK_CELL_COUNT(N)>0) THEN
1273                    IF(TIME_DEPENDENT_FILENAME.AND.MODE==0) THEN
1274                       WRITE(VTU_NAME,20) TRIM(VTK_FILEBASE(VTK_REGION)),FRAME(VTK_REGION),N
1275                    ELSE
1276                       WRITE(VTU_NAME,25) TRIM(VTK_FILEBASE(VTK_REGION)),N
1277                    ENDIF
1278     
1279                    WRITE(PVTU_UNIT,110) '      <Piece Source="',TRIM(VTU_NAME),'"/>'
1280                 ENDIF
1281              ENDDO
1282     
1283     
1284              WRITE(PVTU_UNIT,100) '  </PUnstructuredGrid>'
1285              WRITE(PVTU_UNIT,100) '</VTKFile>'
1286              CLOSE(PVTU_UNIT)
1287           ENDIF
1288     
1289     
1290     20    FORMAT(A,"_",I4.4,"_",I5.5,".vtu")
1291     25    FORMAT(A,"_",I5.5,".vtu")
1292     
1293     100   FORMAT(A)
1294     110   FORMAT(A,A,A)
1295     
1296           RETURN
1297     
1298           END SUBROUTINE CLOSE_VTU_FILE_BIN
1299     
1300     
1301     
1302     
1303     
1304     
1305     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1306     !                                                                      C
1307     !  Module name: OPEN_PVD_FILE                                          C
1308     !  Purpose: Open a PVD file and writes the header                      C
1309     !                                                                      C
1310     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1311     !  Reviewer:                                          Date:            C
1312     !                                                                      C
1313     !  Revision Number #                                  Date: ##-###-##  C
1314     !  Author: #                                                           C
1315     !  Purpose: #                                                          C
1316     !                                                                      C
1317     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1318       SUBROUTINE OPEN_PVD_FILE
1319     
1320           USE param
1321           USE param1
1322           USE parallel
1323           USE constant
1324           USE run
1325           USE toleranc
1326           USE geometry
1327           USE indices
1328           USE compar
1329           USE sendrecv
1330           USE output
1331           USE quadric
1332           USE cutcell
1333           USE fldvar
1334           USE vtk
1335           USE functions
1336     
1337           IMPLICIT NONE
1338           LOGICAL :: PVD_EXISTS
1339     
1340           IF (myPE /= PE_IO) RETURN
1341     
1342           PVD_FILENAME = TRIM(VTK_FILEBASE(VTK_REGION)) // '.pvd'
1343     
1344     ! First, check if the file already exists.
1345     
1346           INQUIRE(FILE=PVD_FILENAME,EXIST=PVD_EXISTS)
1347     
1348     ! The first time this subroutine is executed, properly initialize the pvd file
1349     
1350           IF(.NOT.PVD_FILE_INITIALIZED(VTK_REGION)) THEN
1351     
1352              IF(RUN_TYPE == 'NEW'.OR.RUN_TYPE=='RESTART_2')THEN
1353                 ! For a new or RESTART_2 run, the pvd file should not exist, and is created with appropriate header
1354                 IF (.NOT.PVD_EXISTS) THEN
1355                    OPEN(CONVERT='BIG_ENDIAN',UNIT = PVD_UNIT, FILE = TRIM(PVD_FILENAME))
1356                    WRITE(PVD_UNIT,100) '<?xml version="1.0"?>'
1357                    WRITE(PVD_UNIT,100) '<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">'
1358                    WRITE(PVD_UNIT,100) '<Collection>'
1359     !               CALL UPDATE_AND_CLOSE_PVD_FILE
1360                    PVD_FILE_INITIALIZED(VTK_REGION)=.TRUE.
1361                 ELSE ! If the pvd file exists, print error message and exits
1362                    WRITE(*,1002) TRIM(PVD_FILENAME)
1363                    WRITE(UNIT_LOG, 1002) TRIM(PVD_FILENAME)
1364                    CALL MFIX_EXIT(myPE)
1365                 ENDIF
1366              ELSE
1367                 ! For a restart_1 run, the pvd file must exist
1368                 IF (.NOT.PVD_EXISTS) THEN
1369                    ! If the pvd file does not exist, print error message and exits
1370                    WRITE(*,1003) TRIM(PVD_FILENAME)
1371                    WRITE(UNIT_LOG, 1003) TRIM(PVD_FILENAME)
1372                    CALL MFIX_EXIT(myPE)
1373                 ELSE
1374                ! If it already exists, go to the bottom of the file and prepare to append data (remove last two lines)
1375                    OPEN(CONVERT='BIG_ENDIAN',UNIT=PVD_UNIT,FILE = TRIM(PVD_FILENAME),POSITION="APPEND",STATUS='OLD')
1376                    BACKSPACE(PVD_UNIT)
1377                    BACKSPACE(PVD_UNIT)
1378                    PVD_FILE_INITIALIZED(VTK_REGION)=.TRUE.
1379                 ENDIF
1380              ENDIF
1381           ELSE
1382              ! When properly initialized, open the file and go to the
1383              ! bottom of the file and prepare to append data (remove last two lines)
1384              OPEN(CONVERT='BIG_ENDIAN',UNIT=PVD_UNIT,FILE = TRIM(PVD_FILENAME),POSITION="APPEND",STATUS='OLD')
1385              BACKSPACE(PVD_UNIT)
1386              BACKSPACE(PVD_UNIT)
1387           ENDIF
1388     
1389     
1390     100   FORMAT(A)
1391     
1392     1002  FORMAT(/1X,70('*')/,' From: OPEN_PVD_FILE',/,' Message: ',       &
1393              A,' already exists in the run directory.',/10X,               &
1394                'This is not allowed for a new run.',/10X,                   &
1395              'Terminating run.',/1X,70('*')/)
1396     
1397     1003  FORMAT(/1X,70('*')/,' From: OPEN_PVD_FILE',/,' Message: ',       &
1398              A,' is missing from the  the run directory,',/10X,            &
1399                ' and must be present for a restart run.',/10X,              &
1400              'Terminating run.',/1X,70('*')/)
1401     
1402           RETURN
1403     
1404           END SUBROUTINE OPEN_PVD_FILE
1405     
1406     
1407     
1408     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1409     !                                                                      C
1410     !  Module name: UPDATE_AND_CLOSE_PVD_FILE                              C
1411     !  Purpose: Updates and close a pvd file                               C
1412     !                                                                      C
1413     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1414     !  Reviewer:                                          Date:            C
1415     !                                                                      C
1416     !  Revision Number #                                  Date: ##-###-##  C
1417     !  Author: #                                                           C
1418     !  Purpose: #                                                          C
1419     !                                                                      C
1420     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1421       SUBROUTINE UPDATE_AND_CLOSE_PVD_FILE
1422     
1423           USE compar
1424           use cdist
1425           USE run
1426           USE vtk
1427     
1428           IMPLICIT NONE
1429     
1430           CHARACTER (LEN=255)  :: FILENAME
1431           CHARACTER (LEN=5)   :: EXT
1432     
1433           IF (myPE /= PE_IO) RETURN
1434     
1435           IF(VTK_DATA(VTK_REGION)=='C') THEN
1436              EXT = '.pvtu'
1437           ELSEIF(VTK_DATA(VTK_REGION)=='P') THEN
1438              EXT = '.pvtp'
1439           ENDIF
1440     
1441     
1442           IF(.NOT.BDIST_IO) THEN
1443              FILENAME=VTU_FILENAME
1444           ELSE
1445              IF(TIME_DEPENDENT_FILENAME) THEN
1446                 WRITE(FILENAME,40) TRIM(VTK_FILEBASE(VTK_REGION)),FRAME(VTK_REGION),EXT
1447              ELSE
1448                 WRITE(FILENAME,45) TRIM(VTK_FILEBASE(VTK_REGION)),EXT
1449              ENDIF
1450              IF(TRIM(VTU_DIR)/='.') FILENAME='./'//TRIM(VTU_DIR)//'/'//FILENAME
1451           ENDIF
1452     
1453     
1454     ! Write the data to the file
1455              WRITE(PVD_UNIT,100)&
1456              '<DataSet timestep="',TIME,'" ',& ! simulation time
1457              'group="" part="0" ',& ! necessary file data
1458              'file="',TRIM(FILENAME),'"/>' ! file name of vtp
1459     
1460     ! Write the closing tags
1461              WRITE(PVD_UNIT,110)'</Collection>'
1462              WRITE(PVD_UNIT,110)'</VTKFile>'
1463     
1464              CLOSE(PVD_UNIT)
1465     
1466     
1467     ! 40    FORMAT(A,"_",I4.4,".pvtu")
1468     ! 45    FORMAT(A,".pvtu")
1469     40    FORMAT(A,"_",I4.4,A5)
1470     45    FORMAT(A,A5)
1471     100   FORMAT(6X,A,E14.7,5A)
1472     110   FORMAT(A)
1473     
1474     
1475           RETURN
1476     
1477           END SUBROUTINE UPDATE_AND_CLOSE_PVD_FILE
1478     
1479     
1480     
1481     
1482     
1483     
1484     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1485     !                                                                      C
1486     !  Module name: WRITE_CUT_SURFACE_VTK                                  C
1487     !  Purpose: Writes the cut cell surface in VTK format                  C
1488     !                                                                      C
1489     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1490     !  Reviewer:                                          Date:            C
1491     !                                                                      C
1492     !  Revision Number #                                  Date: ##-###-##  C
1493     !  Author: #                                                           C
1494     !  Purpose: #                                                          C
1495     !                                                                      C
1496     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1497       SUBROUTINE WRITE_CUT_SURFACE_VTK
1498     
1499           USE param
1500           USE param1
1501           USE parallel
1502           USE constant
1503           USE run
1504           USE toleranc
1505           USE geometry
1506           USE indices
1507           USE compar
1508           USE sendrecv
1509           USE quadric
1510           USE cutcell
1511           USE fldvar
1512           USE vtk
1513           USE polygon
1514           USE stl
1515           USE functions
1516     
1517           IMPLICIT NONE
1518     
1519           INTEGER :: L,IJK,NODE
1520           INTEGER :: POINT_ID,POLY_COUNT,FACE_ID,Q_ID
1521           INTEGER :: N_CUT_FACE_NODES
1522     
1523           INTEGER NUMBER_OF_FACES
1524           INTEGER NUMBER_OF_SURFACE_POINTS
1525     
1526           DOUBLE PRECISION, DIMENSION(3,15) :: COORD_CUT_FACE_NODES
1527           DOUBLE PRECISION, DIMENSION(3)    :: NORMAL
1528     
1529           INTEGER, DIMENSION(DIMENSION_MAX_CUT_CELL,6) ::  FACE_CONNECTIVITY
1530           INTEGER, DIMENSION(DIMENSION_MAX_CUT_CELL)   ::  NUMBER_OF_CUT_FACE_POINTS
1531     
1532           DOUBLE PRECISION, DIMENSION(DIMENSION_MAX_CUT_CELL) ::  X_FACE_POINT
1533           DOUBLE PRECISION, DIMENSION(DIMENSION_MAX_CUT_CELL) ::  Y_FACE_POINT
1534           DOUBLE PRECISION, DIMENSION(DIMENSION_MAX_CUT_CELL) ::  Z_FACE_POINT
1535     
1536           DOUBLE PRECISION :: X_COPY,Y_COPY,Z_COPY,F_COPY
1537     
1538           LOGICAL :: CLIP_FLAG
1539     
1540           CHARACTER (LEN=255) :: FILENAME
1541     
1542           LOGICAL :: CORNER_POINT
1543           INTEGER :: NODE_OF_CORNER, IERROR
1544     
1545           IF(myPE/=0) RETURN
1546     
1547     !======================================================================
1548     !  Set-up connectivity for each cell, i.e., regular cells and cut cells
1549     !======================================================================
1550     
1551           POLY_COUNT = 0
1552     
1553           NUMBER_OF_SURFACE_POINTS = 0
1554     
1555           NUMBER_OF_FACES = 0
1556     
1557           DO IJK = 1,IJKMAX3
1558     
1559              IF(GLOBAL_CUT_CELL_AT(IJK)) THEN
1560     
1561     !======================================================================
1562     !  Filter the connectivity to identify nodes belonging to cut face
1563     !======================================================================
1564     
1565     
1566                 NUMBER_OF_FACES = NUMBER_OF_FACES + 1
1567     
1568                 N_CUT_FACE_NODES = 0
1569     
1570                 CALL GET_GLOBAL_CELL_NODE_COORDINATES(IJK,'SCALAR')
1571     
1572                 DO L = 1, GLOBAL_NUMBER_OF_NODES(IJK)
1573                    IF(GLOBAL_CONNECTIVITY(IJK,L)>IJKMAX3) THEN   ! One of the new point
1574                       X_COPY = GLOBAL_X_NEW_POINT(GLOBAL_CONNECTIVITY(IJK,L)-IJKMAX3)
1575                       Y_COPY = GLOBAL_Y_NEW_POINT(GLOBAL_CONNECTIVITY(IJK,L)-IJKMAX3)
1576                       Z_COPY = GLOBAL_Z_NEW_POINT(GLOBAL_CONNECTIVITY(IJK,L)-IJKMAX3)
1577                       CORNER_POINT = .FALSE.
1578                    ELSE                                   ! An existing point
1579                       DO NODE = 1,8
1580                       CORNER_POINT = .TRUE.
1581                          IF(GLOBAL_CONNECTIVITY(IJK,L) == IJK_OF_NODE(NODE)) THEN
1582                             NODE_OF_CORNER = NODE
1583                             X_COPY = X_NODE(NODE)
1584                             Y_COPY = Y_NODE(NODE)
1585                             Z_COPY = Z_NODE(NODE)
1586     
1587                             IF (GLOBAL_SNAP(IJK_OF_NODE(NODE))) THEN ! One of the snapped corner point which now belongs to the cut face
1588                                N_CUT_FACE_NODES = N_CUT_FACE_NODES + 1
1589                                COORD_CUT_FACE_NODES(1,N_CUT_FACE_NODES) = X_COPY
1590                                COORD_CUT_FACE_NODES(2,N_CUT_FACE_NODES) = Y_COPY
1591                                COORD_CUT_FACE_NODES(3,N_CUT_FACE_NODES) = Z_COPY
1592                             ENDIF
1593                          ENDIF
1594                       END DO
1595     
1596                    ENDIF
1597     
1598                    IF(CORNER_POINT) THEN
1599                       Q_ID = 1
1600     
1601                       CALL EVAL_F('QUADRIC',X_COPY,Y_COPY,Z_COPY,Q_ID,F_COPY,CLIP_FLAG)
1602     
1603                       CALL EVAL_F('POLYGON',X_COPY,Y_COPY,Z_COPY,N_POLYGON,F_COPY,CLIP_FLAG)
1604     
1605                       CALL EVAL_F('USR_DEF',X_COPY,Y_COPY,Z_COPY,N_USR_DEF,F_COPY,CLIP_FLAG)
1606     
1607                       IF(USE_STL.OR.USE_MSH) F_COPY = GLOBAL_F_AT(IJK_OF_NODE(NODE_OF_CORNER))
1608     
1609     !                  CALL EVAL_STL_FCT_AT('SCALAR',IJK,NODE_OF_CORNER,F_COPY,CLIP_FLAG,BCID2)
1610                    ELSE
1611                       F_COPY = ZERO
1612                    ENDIF
1613     
1614                    IF (ABS(F_COPY) < TOL_F ) THEN ! belongs to cut face
1615                       N_CUT_FACE_NODES = N_CUT_FACE_NODES + 1
1616                       COORD_CUT_FACE_NODES(1,N_CUT_FACE_NODES) = X_COPY
1617                       COORD_CUT_FACE_NODES(2,N_CUT_FACE_NODES) = Y_COPY
1618                       COORD_CUT_FACE_NODES(3,N_CUT_FACE_NODES) = Z_COPY
1619                    ENDIF
1620     
1621                 END DO
1622     
1623                 CALL REORDER_POLYGON(N_CUT_FACE_NODES,COORD_CUT_FACE_NODES,NORMAL,IERROR)
1624     
1625                 NUMBER_OF_CUT_FACE_POINTS(NUMBER_OF_FACES) = N_CUT_FACE_NODES
1626                 POLY_COUNT = POLY_COUNT + N_CUT_FACE_NODES + 1
1627                 DO NODE = 1,N_CUT_FACE_NODES
1628                    NUMBER_OF_SURFACE_POINTS = NUMBER_OF_SURFACE_POINTS + 1
1629     
1630                    IF(NUMBER_OF_SURFACE_POINTS>=DIMENSION_MAX_CUT_CELL) THEN
1631                       WRITE(*,3000) 'ERROR IN SUBROUTINE WRITE_3DCUT_SURFACE_VTK:'
1632                       WRITE(*,3000) 'NUMBER_OF_SURFACE_POINTS>=DIMENSION_MAX_CUT_CELL:'
1633                       WRITE(*,3000) 'INCREASE VALUE OF FAC_DIM_MAX_CUT_CELL.'
1634                       WRITE(*,3010) 'CURRENT VALUE OF FAC_DIM_MAX_CUT_CELL =',FAC_DIM_MAX_CUT_CELL
1635                       WRITE(*,3020) 'CURRENT VALUE OF DIMENSION_MAX_CUT_CELL =',DIMENSION_MAX_CUT_CELL
1636                       WRITE(*,3000) 'MFiX will exit now.'
1637                       CALL MFIX_EXIT(myPE)
1638                    ENDIF
1639     
1640                    X_FACE_POINT(NUMBER_OF_SURFACE_POINTS) = COORD_CUT_FACE_NODES(1,NODE)
1641                    Y_FACE_POINT(NUMBER_OF_SURFACE_POINTS) = COORD_CUT_FACE_NODES(2,NODE)
1642                    Z_FACE_POINT(NUMBER_OF_SURFACE_POINTS) = COORD_CUT_FACE_NODES(3,NODE)
1643                    FACE_CONNECTIVITY(NUMBER_OF_FACES,NODE) = NUMBER_OF_SURFACE_POINTS
1644                 ENDDO
1645     
1646              ENDIF
1647     
1648           END DO
1649     
1650           FILENAME= TRIM(RUN_NAME) // '_boundary.vtk'
1651           FILENAME = TRIM(FILENAME)
1652           OPEN(CONVERT='BIG_ENDIAN',UNIT = 123, FILE = FILENAME)
1653           WRITE(123,1001)'# vtk DataFile Version 2.0'
1654           WRITE(123,1001)'3D CUT-CELL SURFACE'
1655           WRITE(123,1001)'ASCII'
1656     
1657           IF(NO_K) THEN   ! 2D GEOMETRY
1658              WRITE(123,1001)'DATASET UNSTRUCTURED_GRID'
1659           ELSE            ! 3D GEOMETRY
1660              WRITE(123,1001)'DATASET POLYDATA'
1661           ENDIF
1662     
1663           WRITE(123,1010)'POINTS ',NUMBER_OF_SURFACE_POINTS,' float'
1664     
1665           DO POINT_ID = 1,NUMBER_OF_SURFACE_POINTS
1666              WRITE(123,1020) X_FACE_POINT(POINT_ID),Y_FACE_POINT(POINT_ID),Z_FACE_POINT(POINT_ID)
1667           ENDDO
1668     
1669           IF(NO_K) THEN   ! 2D GEOMETRY
1670     
1671              WRITE(123,1030)'CELLS ',NUMBER_OF_FACES,POLY_COUNT
1672              DO FACE_ID = 1 , NUMBER_OF_FACES
1673                 WRITE(123,1040) NUMBER_OF_CUT_FACE_POINTS(FACE_ID),(FACE_CONNECTIVITY(FACE_ID,L)-1,&
1674                 L=1,NUMBER_OF_CUT_FACE_POINTS(FACE_ID))
1675              ENDDO
1676              WRITE(123,1030)'CELL_TYPES ',NUMBER_OF_FACES
1677              DO FACE_ID = 1 , NUMBER_OF_FACES
1678                 WRITE(123,1040) 3
1679              ENDDO
1680     
1681           ELSE            ! 3D GEOMETRY
1682     
1683              WRITE(123,1030)'POLYGONS ',NUMBER_OF_FACES,POLY_COUNT
1684              DO FACE_ID = 1 , NUMBER_OF_FACES
1685                 WRITE(123,1040) NUMBER_OF_CUT_FACE_POINTS(FACE_ID),(FACE_CONNECTIVITY(FACE_ID,L)-1,&
1686                 L=1,NUMBER_OF_CUT_FACE_POINTS(FACE_ID))
1687              ENDDO
1688     
1689           ENDIF
1690     
1691     1001  FORMAT(A)
1692     1010  FORMAT(A,I8,A)
1693     1020  FORMAT(3(E16.8,2X))
1694     1030  FORMAT(A,2(I8,2X))
1695     1040  FORMAT(20(I8,2X))
1696     3000  FORMAT(1X,A)
1697     3010  FORMAT(1X,A,F8.4)
1698     3020  FORMAT(1X,A,I8)
1699     3030  FORMAT(1X,A,A)
1700           CLOSE (123)
1701     
1702     
1703           WRITE(*,3030)'WROTE BOUNDARY IN VTK FILE : ',FILENAME
1704           RETURN
1705     
1706     
1707           END SUBROUTINE WRITE_CUT_SURFACE_VTK
1708     
1709     
1710     
1711     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1712     !                                                                      C
1713     !  Module name: GATHER_DATA                                            C
1714     !  Purpose: Gather data from all processes in preparation of           C
1715     !           Writing VTK files                                          C
1716     !                                                                      C
1717     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
1718     !  Reviewer:                                          Date:            C
1719     !                                                                      C
1720     !  Revision Number #                                  Date: ##-###-##  C
1721     !  Author: #                                                           C
1722     !  Purpose: #                                                          C
1723     !                                                                      C
1724     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1725       SUBROUTINE GATHER_DATA
1726     
1727           USE param
1728           USE param1
1729           USE parallel
1730           USE constant
1731           USE run
1732           USE toleranc
1733           USE geometry
1734           USE indices
1735           USE compar
1736           USE mpi_utility
1737           USE sendrecv
1738           USE quadric
1739           USE cutcell
1740           USE fldvar
1741           USE vtk
1742           USE functions
1743     
1744           IMPLICIT NONE
1745     
1746           INTEGER :: IJK,I,J,K,L
1747           INTEGER :: IJK_OFFSET
1748     
1749           INTEGER :: iproc,IERR
1750           INTEGER, DIMENSION(0:numPEs-1) :: disp,rcount
1751           INTEGER, DIMENSION(:,:), ALLOCATABLE :: SHIFTED_CONNECTIVITY
1752     
1753     !======================================================================
1754     !  parallel processing
1755     !======================================================================
1756     
1757           CALL allgather_1i (NUMBER_OF_NEW_POINTS,rcount,IERR)
1758     
1759           IF (myPE == 0) THEN
1760              IJK_OFFSET = 0
1761           ELSE
1762              IJK_OFFSET = 0
1763              DO iproc=0,myPE-1
1764                 IJK_OFFSET = IJK_OFFSET + rcount(iproc)
1765              ENDDO
1766           ENDIF
1767     
1768           CALL allgather_1i (IJK_OFFSET,disp,IERR)
1769     
1770           IF(.NOT.GLOBAL_VAR_ALLOCATED) THEN
1771     
1772              IF (myPE == PE_IO) THEN
1773                 allocate (GLOBAL_I_OF(ijkmax3))
1774                 allocate (GLOBAL_J_OF(ijkmax3))
1775                 allocate (GLOBAL_K_OF(ijkmax3))
1776                 allocate (GLOBAL_CONNECTIVITY(ijkmax3,15))
1777                 allocate (GLOBAL_NUMBER_OF_NODES(ijkmax3))
1778                 allocate (GLOBAL_INTERIOR_CELL_AT(ijkmax3))
1779                 allocate (GLOBAL_BLOCKED_CELL_AT(ijkmax3))
1780                 allocate (GLOBAL_STANDARD_CELL_AT(ijkmax3))
1781                 allocate (GLOBAL_CUT_CELL_AT(ijkmax3))
1782                 allocate (GLOBAL_SNAP(ijkmax3))
1783                 allocate (GLOBAL_X_NEW_POINT(ijkmax3))
1784                 allocate (GLOBAL_Y_NEW_POINT(ijkmax3))
1785                 allocate (GLOBAL_Z_NEW_POINT(ijkmax3))
1786                 allocate (GLOBAL_F_AT(ijkmax3))
1787     
1788              ELSE
1789                 allocate (GLOBAL_I_OF(1))
1790                 allocate (GLOBAL_J_OF(1))
1791                 allocate (GLOBAL_K_OF(1))
1792                 allocate (GLOBAL_CONNECTIVITY(1,15))
1793                 allocate (GLOBAL_NUMBER_OF_NODES(1))
1794                 allocate (GLOBAL_INTERIOR_CELL_AT(1))
1795                 allocate (GLOBAL_BLOCKED_CELL_AT(1))
1796                 allocate (GLOBAL_STANDARD_CELL_AT(1))
1797                 allocate (GLOBAL_CUT_CELL_AT(1))
1798                 allocate (GLOBAL_SNAP(1))
1799                 allocate (GLOBAL_X_NEW_POINT(1))
1800                 allocate (GLOBAL_Y_NEW_POINT(1))
1801                 allocate (GLOBAL_Z_NEW_POINT(1))
1802                 allocate (GLOBAL_F_AT(1))
1803              ENDIF
1804     
1805              GLOBAL_VAR_ALLOCATED = .TRUE.
1806     
1807           ENDIF
1808     
1809           IF(numPEs==1) THEN  ! Serial run
1810              GLOBAL_X_NEW_POINT(1:NUMBER_OF_NEW_POINTS) =  X_NEW_POINT(1:NUMBER_OF_NEW_POINTS)
1811              GLOBAL_Y_NEW_POINT(1:NUMBER_OF_NEW_POINTS) =  Y_NEW_POINT(1:NUMBER_OF_NEW_POINTS)
1812              IF(DO_K) GLOBAL_Z_NEW_POINT(1:NUMBER_OF_NEW_POINTS) =  Z_NEW_POINT(1:NUMBER_OF_NEW_POINTS)
1813           ELSE !Parallel run
1814              call gatherv_1d( X_NEW_POINT, NUMBER_OF_NEW_POINTS, GLOBAL_X_NEW_POINT, rcount, disp, PE_IO, ierr )
1815              call gatherv_1d( Y_NEW_POINT, NUMBER_OF_NEW_POINTS, GLOBAL_Y_NEW_POINT, rcount, disp, PE_IO, ierr )
1816              call gatherv_1d( Z_NEW_POINT, NUMBER_OF_NEW_POINTS, GLOBAL_Z_NEW_POINT, rcount, disp, PE_IO, ierr )
1817           ENDIF
1818     
1819           call global_sum(NUMBER_OF_NEW_POINTS, GLOBAL_NUMBER_OF_NEW_POINTS,  PE_IO, ierr )
1820     
1821           Allocate(  SHIFTED_CONNECTIVITY  (DIMENSION_3,15) )
1822     
1823           SHIFTED_CONNECTIVITY = CONNECTIVITY
1824     
1825           WHERE (SHIFTED_CONNECTIVITY > IJKEND3)
1826              SHIFTED_CONNECTIVITY = SHIFTED_CONNECTIVITY - IJKEND3 + IJKMAX3 + disp(myPE)
1827           END WHERE
1828     
1829           DO IJK = IJKSTART3,IJKEND3
1830              DO L=1,NUMBER_OF_NODES(IJK)
1831                 IF(CONNECTIVITY(IJK,L) <= IJKEND3) THEN
1832                    I = I_OF(CONNECTIVITY(IJK,L))
1833                    J = J_OF(CONNECTIVITY(IJK,L))
1834                    K = K_OF(CONNECTIVITY(IJK,L))
1835                    SHIFTED_CONNECTIVITY(IJK,L) = funijk_gl(I,J,K)
1836                 ENDIF
1837              ENDDO
1838           ENDDO
1839     
1840     
1841           GLOBAL_INTERIOR_CELL_AT = .FALSE.
1842           GLOBAL_BLOCKED_CELL_AT  = .FALSE.
1843           GLOBAL_CUT_CELL_AT      = .FALSE.
1844           call gather (I_OF,GLOBAL_I_OF,root)
1845           call gather (J_OF,GLOBAL_J_OF,root)
1846           call gather (K_OF,GLOBAL_K_OF,root)
1847           call gather (SHIFTED_CONNECTIVITY,GLOBAL_CONNECTIVITY,root)
1848           call gather (NUMBER_OF_NODES,GLOBAL_NUMBER_OF_NODES,root)
1849           call gather (INTERIOR_CELL_AT,GLOBAL_INTERIOR_CELL_AT,root)
1850           call gather (BLOCKED_CELL_AT,GLOBAL_BLOCKED_CELL_AT,root)
1851           call gather (STANDARD_CELL_AT,GLOBAL_STANDARD_CELL_AT,root)
1852           call gather (CUT_CELL_AT,GLOBAL_CUT_CELL_AT,root)
1853           call gather (SNAP,GLOBAL_SNAP,root)
1854           call gather (F_AT,GLOBAL_F_AT,root)
1855     
1856           Deallocate(  SHIFTED_CONNECTIVITY )
1857     
1858           IF (myPE == PE_IO) THEN
1859     
1860              POLY_COUNTER = 0
1861     
1862              NUMBER_OF_CELLS = 0
1863     
1864              NUMBER_OF_CUT_CELLS = 0
1865     
1866              NUMBER_OF_BLOCKED_CELLS = 0
1867     
1868              NUMBER_OF_STANDARD_CELLS = 0
1869     
1870              DO IJK = 1, IJKMAX3
1871     
1872                 IF (GLOBAL_INTERIOR_CELL_AT(IJK)) THEN
1873     
1874                    NUMBER_OF_CELLS = NUMBER_OF_CELLS + 1
1875     
1876                    IF (GLOBAL_BLOCKED_CELL_AT(IJK))  NUMBER_OF_BLOCKED_CELLS  = NUMBER_OF_BLOCKED_CELLS + 1
1877                    IF (GLOBAL_STANDARD_CELL_AT(IJK)) NUMBER_OF_STANDARD_CELLS = NUMBER_OF_STANDARD_CELLS + 1
1878                    IF (GLOBAL_CUT_CELL_AT(IJK))      NUMBER_OF_CUT_CELLS = NUMBER_OF_CUT_CELLS + 1
1879     
1880                    IF (.NOT.GLOBAL_BLOCKED_CELL_AT(IJK))   POLY_COUNTER = POLY_COUNTER + GLOBAL_NUMBER_OF_NODES(IJK) + 1
1881     
1882                 ENDIF
1883     
1884              END DO
1885     
1886     
1887              NUMBER_OF_POINTS = IJKMAX3 + GLOBAL_NUMBER_OF_NEW_POINTS
1888     
1889           ENDIF
1890     
1891           RETURN
1892     
1893     
1894           END SUBROUTINE GATHER_DATA
1895     
1896     
1897     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1898     !                                                                      C
1899     !  Module name: SETUP_VTK_NO_CUTCELL                                   C
1900     !  Purpose: Setup VTK data for the regular grid (no cut cells)         C
1901     !           This i scalled when CARTESIAN_GRID is .FALSE.              C
1902     !           and WRITE>VTK_FILES is .TRUE.              .               C
1903     !                                                                      C
1904     !  Author: Jeff Dietiker                              Date: 14-Jan-15  C
1905     !  Reviewer:                                          Date:            C
1906     !                                                                      C
1907     !  Revision Number #                                  Date: ##-###-##  C
1908     !  Author: #                                                           C
1909     !  Purpose: #                                                          C
1910     !                                                                      C
1911     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1912       SUBROUTINE SETUP_VTK_NO_CUTCELL
1913     
1914           USE param
1915           USE param1
1916           USE parallel
1917           USE constant
1918           USE run
1919           USE toleranc
1920           USE geometry
1921           USE indices
1922           USE compar
1923           USE mpi_utility
1924           USE sendrecv
1925           USE quadric
1926           USE cutcell
1927           USE fldvar
1928           USE vtk
1929           USE functions
1930     
1931           IMPLICIT NONE
1932     
1933           INTEGER :: IJK,I,J,K,L,NODE
1934     
1935           INTEGER, DIMENSION(:,:), ALLOCATABLE :: SHIFTED_CONNECTIVITY
1936     
1937     ! Only a few arrays need to be allocated here simce we do not need
1938     ! all Cartesian-grid arrays
1939           IF(.NOT.ALLOCATED(XG_E)) Allocate( XG_E(0:DIMENSION_I) )
1940           IF(.NOT.ALLOCATED(YG_N)) Allocate( YG_N(0:DIMENSION_J) )
1941           IF(.NOT.ALLOCATED(ZG_T)) Allocate( ZG_T(0:DIMENSION_K) )
1942     
1943           IF(.NOT.ALLOCATED(INTERIOR_CELL_AT)) THEN
1944              Allocate(  INTERIOR_CELL_AT  (DIMENSION_3) )
1945              INTERIOR_CELL_AT = .FALSE.
1946           ENDIF
1947     
1948           IF(.NOT.ALLOCATED(BLOCKED_CELL_AT)) THEN
1949              Allocate(  BLOCKED_CELL_AT  (DIMENSION_3) )
1950              BLOCKED_CELL_AT = .FALSE.
1951           ENDIF
1952     
1953           IF(.NOT.ALLOCATED(STANDARD_CELL_AT)) THEN
1954              Allocate(  STANDARD_CELL_AT  (DIMENSION_3) )
1955              STANDARD_CELL_AT = .TRUE.
1956           ENDIF
1957     
1958           IF(.NOT.ALLOCATED(CUT_CELL_AT)) THEN
1959              Allocate(  CUT_CELL_AT  (DIMENSION_3) )
1960              CUT_CELL_AT = .FALSE.
1961           ENDIF
1962     
1963           IF(.NOT.ALLOCATED(NUMBER_OF_NODES))  Allocate(  NUMBER_OF_NODES  (DIMENSION_3) )
1964           NUMBER_OF_NODES= 0
1965     
1966           IF(.NOT.ALLOCATED(CONNECTIVITY))     Allocate(  CONNECTIVITY  (DIMENSION_3,15) )
1967     
1968     
1969     ! This is a shoter version of Get_cut_cell_Flags
1970           DO IJK = IJKSTART3, IJKEND3
1971     
1972     !======================================================================
1973     !  Get coordinates of eight nodes
1974     !======================================================================
1975     
1976              CALL GET_CELL_NODE_COORDINATES(IJK,'SCALAR')
1977     
1978              I = I_OF(IJK)
1979              J = J_OF(IJK)
1980              K = K_OF(IJK)
1981     
1982              IF(NO_K) THEN   ! 2D case
1983     
1984                 INTERIOR_CELL_AT(IJK) = (     (I >= ISTART1 ).AND.(I <= IEND1 )  &
1985                                          .AND.(J >= JSTART1 ).AND.(J <= JEND1 ) )
1986     
1987              ELSE            ! 3D case
1988     
1989                 INTERIOR_CELL_AT(IJK) = (     (I >= ISTART1 ).AND.(I <= IEND1 )  &
1990                                          .AND.(J >= JSTART1 ).AND.(J <= JEND1 )  &
1991                                          .AND.(K >= KSTART1 ).AND.(K <= KEND1 ) )
1992     
1993              ENDIF
1994     
1995     
1996              IF(INTERIOR_CELL_AT(IJK)) THEN
1997     
1998     ! Set up the connectivity: 4 nodes in 2D, 8 nodes in 3D
1999                 IF(NO_K) THEN
2000                    NUMBER_OF_NODES(IJK) = 4
2001                    CONNECTIVITY(IJK,1) = IJK_OF_NODE(5)
2002                    CONNECTIVITY(IJK,2) = IJK_OF_NODE(6)
2003                    CONNECTIVITY(IJK,3) = IJK_OF_NODE(8)
2004                    CONNECTIVITY(IJK,4) = IJK_OF_NODE(7)
2005                 ELSE
2006                    NUMBER_OF_NODES(IJK) = 8
2007                    DO NODE = 1,8
2008                       CONNECTIVITY(IJK,NODE) = IJK_OF_NODE(NODE)
2009                    END DO
2010                 ENDIF
2011     
2012     ! If obstacles are defined, they will be flagged as blocked cells
2013     ! and will not be visible in the VTK files
2014                 IF(WALL_AT(IJK)) THEN
2015                    BLOCKED_CELL_AT(IJK)  = .TRUE.
2016                    CUT_CELL_AT(IJK)      = .FALSE.
2017                    STANDARD_CELL_AT(IJK) = .FALSE.
2018                 ENDIF
2019     
2020              ENDIF
2021     
2022           ENDDO
2023     
2024     
2025     !======================================================================
2026     !  parallel processing
2027     !======================================================================
2028           call SEND_RECEIVE_1D_LOGICAL(STANDARD_CELL_AT,2)
2029           call SEND_RECEIVE_1D_LOGICAL(BLOCKED_CELL_AT,2)
2030           call SEND_RECEIVE_1D_LOGICAL(CUT_CELL_AT,2)
2031     
2032           Allocate(  SHIFTED_CONNECTIVITY  (DIMENSION_3,15) )
2033     
2034           SHIFTED_CONNECTIVITY = CONNECTIVITY
2035     
2036     ! Replace local node index by global node index before gathering the array
2037           DO IJK = IJKSTART3,IJKEND3
2038              DO L=1,NUMBER_OF_NODES(IJK)
2039                 IF(CONNECTIVITY(IJK,L) <= IJKEND3) THEN
2040                    I = I_OF(CONNECTIVITY(IJK,L))
2041                    J = J_OF(CONNECTIVITY(IJK,L))
2042                    K = K_OF(CONNECTIVITY(IJK,L))
2043                    SHIFTED_CONNECTIVITY(IJK,L) = funijk_gl(I,J,K)
2044                 ENDIF
2045              ENDDO
2046           ENDDO
2047     
2048     ! Allocate, initialize and gather arrays
2049           IF(.NOT.GLOBAL_VAR_ALLOCATED) THEN
2050     
2051              IF (myPE == PE_IO) THEN
2052                 allocate (GLOBAL_I_OF(ijkmax3))
2053                 allocate (GLOBAL_J_OF(ijkmax3))
2054                 allocate (GLOBAL_K_OF(ijkmax3))
2055                 allocate (GLOBAL_CONNECTIVITY(ijkmax3,15))
2056                 allocate (GLOBAL_NUMBER_OF_NODES(ijkmax3))
2057                 allocate (GLOBAL_INTERIOR_CELL_AT(ijkmax3))
2058                 allocate (GLOBAL_BLOCKED_CELL_AT(ijkmax3))
2059                 allocate (GLOBAL_STANDARD_CELL_AT(ijkmax3))
2060                 allocate (GLOBAL_CUT_CELL_AT(ijkmax3))
2061     
2062              ELSE
2063                 allocate (GLOBAL_I_OF(1))
2064                 allocate (GLOBAL_J_OF(1))
2065                 allocate (GLOBAL_K_OF(1))
2066                 allocate (GLOBAL_CONNECTIVITY(1,15))
2067                 allocate (GLOBAL_NUMBER_OF_NODES(1))
2068                 allocate (GLOBAL_INTERIOR_CELL_AT(1))
2069                 allocate (GLOBAL_BLOCKED_CELL_AT(1))
2070                 allocate (GLOBAL_STANDARD_CELL_AT(1))
2071                 allocate (GLOBAL_CUT_CELL_AT(1))
2072              ENDIF
2073     
2074              GLOBAL_VAR_ALLOCATED = .TRUE.
2075     
2076           ENDIF
2077     
2078     
2079           GLOBAL_INTERIOR_CELL_AT = .FALSE.
2080           GLOBAL_BLOCKED_CELL_AT  = .FALSE.
2081           GLOBAL_CUT_CELL_AT      = .FALSE.
2082           GLOBAL_STANDARD_CELL_AT = .TRUE.
2083     
2084           call gather (I_OF,GLOBAL_I_OF,root)
2085           call gather (J_OF,GLOBAL_J_OF,root)
2086           call gather (K_OF,GLOBAL_K_OF,root)
2087           call gather (SHIFTED_CONNECTIVITY,GLOBAL_CONNECTIVITY,root)
2088           call gather (NUMBER_OF_NODES,GLOBAL_NUMBER_OF_NODES,root)
2089           call gather (INTERIOR_CELL_AT,GLOBAL_INTERIOR_CELL_AT,root)
2090           call gather (BLOCKED_CELL_AT,GLOBAL_BLOCKED_CELL_AT,root)
2091           call gather (STANDARD_CELL_AT,GLOBAL_STANDARD_CELL_AT,root)
2092           call gather (CUT_CELL_AT,GLOBAL_CUT_CELL_AT,root)
2093     
2094           deAllocate(  SHIFTED_CONNECTIVITY   )
2095     
2096     
2097     
2098     ! Count the number of cells
2099           GLOBAL_NUMBER_OF_NEW_POINTS = 0
2100     
2101           IF (myPE == PE_IO) THEN
2102     
2103              POLY_COUNTER = 0
2104     
2105              NUMBER_OF_CELLS = 0
2106     
2107              NUMBER_OF_CUT_CELLS = 0
2108     
2109              NUMBER_OF_BLOCKED_CELLS = 0
2110     
2111              NUMBER_OF_STANDARD_CELLS = 0
2112     
2113              DO IJK = 1, IJKMAX3
2114     
2115                 IF (GLOBAL_INTERIOR_CELL_AT(IJK)) THEN
2116     
2117                    NUMBER_OF_CELLS = NUMBER_OF_CELLS + 1
2118     
2119                    IF (GLOBAL_BLOCKED_CELL_AT(IJK))  NUMBER_OF_BLOCKED_CELLS  = NUMBER_OF_BLOCKED_CELLS + 1
2120                    IF (GLOBAL_STANDARD_CELL_AT(IJK)) NUMBER_OF_STANDARD_CELLS = NUMBER_OF_STANDARD_CELLS + 1
2121                    IF (GLOBAL_CUT_CELL_AT(IJK))      NUMBER_OF_CUT_CELLS = NUMBER_OF_CUT_CELLS + 1
2122     
2123                    IF (.NOT.GLOBAL_BLOCKED_CELL_AT(IJK))   POLY_COUNTER = POLY_COUNTER + GLOBAL_NUMBER_OF_NODES(IJK) + 1
2124     
2125                 ENDIF
2126     
2127              END DO
2128     
2129     ! There are no new points since there a no cut cells
2130              NUMBER_OF_POINTS = IJKMAX3
2131     
2132           ENDIF
2133     
2134           RETURN
2135     
2136     
2137           END SUBROUTINE SETUP_VTK_NO_CUTCELL
2138     
2139     
2140     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2141     !                                                                      C
2142     !  Module name: PRINT_GRID_STATISTICS                                  C
2143     !  Purpose: PRINT_GRID_STATISTICS ON SCREEN                            C
2144     !                                                                      C
2145     !                                                                      C
2146     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
2147     !  Reviewer:                                          Date:            C
2148     !                                                                      C
2149     !  Revision Number #                                  Date: ##-###-##  C
2150     !  Author: #                                                           C
2151     !  Purpose: #                                                          C
2152     !                                                                      C
2153     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
2154       SUBROUTINE PRINT_GRID_STATISTICS
2155     
2156           USE param
2157           USE param1
2158           USE parallel
2159           USE constant
2160           USE run
2161           USE toleranc
2162           USE geometry
2163           USE indices
2164           USE compar
2165           USE mpi_utility
2166           USE sendrecv
2167           USE quadric
2168           USE cutcell
2169           USE fldvar
2170           USE vtk
2171           USE functions
2172     
2173           IMPLICIT NONE
2174     
2175           INTEGER :: IJK
2176     
2177           INTEGER :: IERR
2178     
2179           DOUBLE PRECISION :: MIN_VOL, MAX_VOL, GLOBAL_MIN_VOL,GLOBAL_MAX_VOL
2180           DOUBLE PRECISION :: MIN_AYZ, MAX_AYZ, GLOBAL_MIN_AYZ,GLOBAL_MAX_AYZ
2181           DOUBLE PRECISION :: MIN_AXZ, MAX_AXZ, GLOBAL_MIN_AXZ,GLOBAL_MAX_AXZ
2182           DOUBLE PRECISION :: MIN_AXY, MAX_AXY, GLOBAL_MIN_AXY,GLOBAL_MAX_AXY
2183           DOUBLE PRECISION :: MIN_CUT, MAX_CUT, GLOBAL_MIN_CUT,GLOBAL_MAX_CUT
2184           DOUBLE PRECISION :: LOCAL_MIN_Q,LOCAL_MAX_Q, GLOBAL_MIN_Q,GLOBAL_MAX_Q
2185     
2186           IF (myPE == PE_IO) THEN
2187     
2188              IF(.NOT.GRID_INFO_PRINTED_ON_SCREEN) THEN
2189                 WRITE(*,5) 'GRID STATISTICS:'
2190                 WRITE(*,5) 'NUMBER OF CELLS          = ', NUMBER_OF_CELLS
2191                 WRITE(*,10)'NUMBER OF STANDARD CELLS = ', &
2192                             NUMBER_OF_STANDARD_CELLS,DBLE(NUMBER_OF_STANDARD_CELLS) / DBLE(NUMBER_OF_CELLS) * 100.0D0
2193                 WRITE(*,10)'NUMBER OF CUT CELLS      = ', &
2194                             NUMBER_OF_CUT_CELLS,DBLE(NUMBER_OF_CUT_CELLS) / DBLE(NUMBER_OF_CELLS) * 100.0D0
2195                 WRITE(*,10)'NUMBER OF BLOCKED CELLS  = ', &
2196                             NUMBER_OF_BLOCKED_CELLS,DBLE(NUMBER_OF_BLOCKED_CELLS) / DBLE(NUMBER_OF_CELLS) * 100.0D0
2197     
2198     5           FORMAT(1X,A,I8)
2199     10          FORMAT(1X,A,I8,' (',F6.2,' % of Total)')
2200     
2201              ENDIF
2202     
2203              GRID_INFO_PRINTED_ON_SCREEN = .TRUE.
2204     
2205           ENDIF
2206     
2207     
2208     !======================================================================
2209     !  Scalar Cell volumes and areas
2210     !======================================================================
2211     
2212           MIN_VOL =   LARGE_NUMBER
2213           MAX_VOL = - LARGE_NUMBER
2214           MIN_AYZ =   LARGE_NUMBER
2215           MAX_AYZ = - LARGE_NUMBER
2216           MIN_AXZ =   LARGE_NUMBER
2217           MAX_AXZ = - LARGE_NUMBER
2218           MIN_AXY =   LARGE_NUMBER
2219           MAX_AXY = - LARGE_NUMBER
2220     
2221           DO IJK = IJKSTART3, IJKEND3
2222              IF(STANDARD_CELL_AT(IJK)) THEN              ! STANDARD CELLS
2223     
2224                 MIN_VOL =   DMIN1(MIN_VOL,VOL(IJK))
2225                 MAX_VOL =   DMAX1(MAX_VOL,VOL(IJK))
2226                 MIN_AYZ =   DMIN1(MIN_AYZ,AYZ(IJK))
2227                 MAX_AYZ =   DMAX1(MAX_AYZ,AYZ(IJK))
2228                 MIN_AXZ =   DMIN1(MIN_AXZ,AXZ(IJK))
2229                 MAX_AXZ =   DMAX1(MAX_AXZ,AXZ(IJK))
2230                 MIN_AXY =   DMIN1(MIN_AXY,AXY(IJK))
2231                 MAX_AXY =   DMAX1(MAX_AXY,AXY(IJK))
2232     
2233              ENDIF
2234           END DO
2235     
2236           call global_min(MIN_VOL, GLOBAL_MIN_VOL,  PE_IO, ierr )
2237           call global_max(MAX_VOL, GLOBAL_MAX_VOL,  PE_IO, ierr )
2238           call global_min(MIN_AYZ, GLOBAL_MIN_AYZ,  PE_IO, ierr )
2239           call global_max(MAX_AYZ, GLOBAL_MAX_AYZ,  PE_IO, ierr )
2240           call global_min(MIN_AXZ, GLOBAL_MIN_AXZ,  PE_IO, ierr )
2241           call global_max(MAX_AXZ, GLOBAL_MAX_AXZ,  PE_IO, ierr )
2242           call global_min(MIN_AXY, GLOBAL_MIN_AXY,  PE_IO, ierr )
2243           call global_max(MAX_AXY, GLOBAL_MAX_AXY,  PE_IO, ierr )
2244     
2245           IF (myPE == PE_IO) THEN
2246              WRITE(UNIT_CUT_CELL_LOG,1000)  '################################################################'
2247              WRITE(UNIT_CUT_CELL_LOG,1000)  '                       CELLS STATISTICS                         '
2248              WRITE(UNIT_CUT_CELL_LOG,1000)  '################################################################'
2249              WRITE(UNIT_CUT_CELL_LOG,1000)  'SCALAR STANDARD CELLS:'
2250              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXY                   = ', GLOBAL_MIN_AXY,GLOBAL_MAX_AXY
2251              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXZ                   = ', GLOBAL_MIN_AXZ,GLOBAL_MAX_AXZ
2252              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AYZ                   = ', GLOBAL_MIN_AYZ,GLOBAL_MAX_AYZ
2253              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF VOLUME                = ', GLOBAL_MIN_VOL,GLOBAL_MAX_VOL
2254           ENDIF
2255     
2256     
2257           MIN_VOL =   LARGE_NUMBER
2258           MAX_VOL = - LARGE_NUMBER
2259           MIN_AYZ =   LARGE_NUMBER
2260           MAX_AYZ = - LARGE_NUMBER
2261           MIN_AXZ =   LARGE_NUMBER
2262           MAX_AXZ = - LARGE_NUMBER
2263           MIN_AXY =   LARGE_NUMBER
2264           MAX_AXY = - LARGE_NUMBER
2265     
2266           DO IJK = IJKSTART3, IJKEND3
2267              IF(CUT_CELL_AT(IJK)) THEN                   ! CUT CELLS
2268     
2269                 MIN_VOL =   DMIN1(MIN_VOL,VOL(IJK))
2270                 MAX_VOL =   DMAX1(MAX_VOL,VOL(IJK))
2271                 MIN_AYZ =   DMIN1(MIN_AYZ,AYZ(IJK))
2272                 MAX_AYZ =   DMAX1(MAX_AYZ,AYZ(IJK))
2273                 MIN_AXZ =   DMIN1(MIN_AXZ,AXZ(IJK))
2274                 MAX_AXZ =   DMAX1(MAX_AXZ,AXZ(IJK))
2275                 MIN_AXY =   DMIN1(MIN_AXY,AXY(IJK))
2276                 MAX_AXY =   DMAX1(MAX_AXY,AXY(IJK))
2277     
2278              ENDIF
2279           END DO
2280     
2281           call global_min(MIN_VOL, GLOBAL_MIN_VOL,  PE_IO, ierr )
2282           call global_max(MAX_VOL, GLOBAL_MAX_VOL,  PE_IO, ierr )
2283           call global_min(MIN_AYZ, GLOBAL_MIN_AYZ,  PE_IO, ierr )
2284           call global_max(MAX_AYZ, GLOBAL_MAX_AYZ,  PE_IO, ierr )
2285           call global_min(MIN_AXZ, GLOBAL_MIN_AXZ,  PE_IO, ierr )
2286           call global_max(MAX_AXZ, GLOBAL_MAX_AXZ,  PE_IO, ierr )
2287           call global_min(MIN_AXY, GLOBAL_MIN_AXY,  PE_IO, ierr )
2288           call global_max(MAX_AXY, GLOBAL_MAX_AXY,  PE_IO, ierr )
2289     
2290           IF (myPE == PE_IO) THEN
2291              WRITE(UNIT_CUT_CELL_LOG,1000)  'SCALAR CUT CELLS:'
2292              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXY                   = ', GLOBAL_MIN_AXY,GLOBAL_MAX_AXY
2293              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXZ                   = ', GLOBAL_MIN_AXZ,GLOBAL_MAX_AXZ
2294              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AYZ                   = ', GLOBAL_MIN_AYZ,GLOBAL_MAX_AYZ
2295              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF VOLUME                = ', GLOBAL_MIN_VOL,GLOBAL_MAX_VOL
2296              WRITE(UNIT_CUT_CELL_LOG,1010)  'NUMBER OF SMALL SCALAR CELLS   = ', NUMBER_OF_SMALL_CELLS
2297              WRITE(UNIT_CUT_CELL_LOG,1000)  '################################################################'
2298           ENDIF
2299     
2300     1000 FORMAT(A,E14.7,2X,E14.7)
2301     1010 FORMAT(A,I8)
2302     
2303     !======================================================================
2304     !  U-Momentum Cell volumes and areas
2305     !======================================================================
2306     
2307           MIN_VOL =   LARGE_NUMBER
2308           MAX_VOL = - LARGE_NUMBER
2309           MIN_AYZ =   LARGE_NUMBER
2310           MAX_AYZ = - LARGE_NUMBER
2311           MIN_AXZ =   LARGE_NUMBER
2312           MAX_AXZ = - LARGE_NUMBER
2313           MIN_AXY =   LARGE_NUMBER
2314           MAX_AXY = - LARGE_NUMBER
2315     
2316           DO IJK = IJKSTART3, IJKEND3
2317              IF(STANDARD_U_CELL_AT(IJK)) THEN              ! STANDARD CELLS
2318     
2319                 MIN_VOL =   DMIN1(MIN_VOL,VOL_U(IJK))
2320                 MAX_VOL =   DMAX1(MAX_VOL,VOL_U(IJK))
2321                 MIN_AYZ =   DMIN1(MIN_AYZ,AYZ_U(IJK))
2322                 MAX_AYZ =   DMAX1(MAX_AYZ,AYZ_U(IJK))
2323                 MIN_AXZ =   DMIN1(MIN_AXZ,AXZ_U(IJK))
2324                 MAX_AXZ =   DMAX1(MAX_AXZ,AXZ_U(IJK))
2325                 MIN_AXY =   DMIN1(MIN_AXY,AXY_U(IJK))
2326                 MAX_AXY =   DMAX1(MAX_AXY,AXY_U(IJK))
2327     
2328              ENDIF
2329           END DO
2330     
2331           call global_min(MIN_VOL, GLOBAL_MIN_VOL,  PE_IO, ierr )
2332           call global_max(MAX_VOL, GLOBAL_MAX_VOL,  PE_IO, ierr )
2333           call global_min(MIN_AYZ, GLOBAL_MIN_AYZ,  PE_IO, ierr )
2334           call global_max(MAX_AYZ, GLOBAL_MAX_AYZ,  PE_IO, ierr )
2335           call global_min(MIN_AXZ, GLOBAL_MIN_AXZ,  PE_IO, ierr )
2336           call global_max(MAX_AXZ, GLOBAL_MAX_AXZ,  PE_IO, ierr )
2337           call global_min(MIN_AXY, GLOBAL_MIN_AXY,  PE_IO, ierr )
2338           call global_max(MAX_AXY, GLOBAL_MAX_AXY,  PE_IO, ierr )
2339     
2340           IF (myPE == PE_IO) THEN
2341              WRITE(UNIT_CUT_CELL_LOG,1000)  'U-MOMENTUM STANDARD CELLS:'
2342              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXY                   = ', GLOBAL_MIN_AXY,GLOBAL_MAX_AXY
2343              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXZ                   = ', GLOBAL_MIN_AXZ,GLOBAL_MAX_AXZ
2344              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AYZ                   = ', GLOBAL_MIN_AYZ,GLOBAL_MAX_AYZ
2345              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF VOLUME                = ', GLOBAL_MIN_VOL,GLOBAL_MAX_VOL
2346           ENDIF
2347     
2348           MIN_VOL =   LARGE_NUMBER
2349           MAX_VOL = - LARGE_NUMBER
2350           MIN_AYZ =   LARGE_NUMBER
2351           MAX_AYZ = - LARGE_NUMBER
2352           MIN_AXZ =   LARGE_NUMBER
2353           MAX_AXZ = - LARGE_NUMBER
2354           MIN_AXY =   LARGE_NUMBER
2355           MAX_AXY = - LARGE_NUMBER
2356           MIN_CUT =   LARGE_NUMBER
2357           MAX_CUT = - LARGE_NUMBER
2358     
2359           DO IJK = IJKSTART3, IJKEND3
2360              IF(CUT_U_CELL_AT(IJK).AND.(.NOT.WALL_U_AT(IJK))) THEN      ! CUT CELLS
2361     
2362                 MIN_VOL =   DMIN1(MIN_VOL,VOL_U(IJK))
2363                 MAX_VOL =   DMAX1(MAX_VOL,VOL_U(IJK))
2364                 MIN_AYZ =   DMIN1(MIN_AYZ,AYZ_U(IJK))
2365                 MAX_AYZ =   DMAX1(MAX_AYZ,AYZ_U(IJK))
2366                 MIN_AXZ =   DMIN1(MIN_AXZ,AXZ_U(IJK))
2367                 MAX_AXZ =   DMAX1(MAX_AXZ,AXZ_U(IJK))
2368                 MIN_AXY =   DMIN1(MIN_AXY,AXY_U(IJK))
2369                 MAX_AXY =   DMAX1(MAX_AXY,AXY_U(IJK))
2370                 MIN_CUT =   DMIN1(MIN_CUT,AREA_U_CUT(IJK))
2371                 MAX_CUT =   DMAX1(MAX_CUT,AREA_U_CUT(IJK))
2372     
2373              ENDIF
2374           END DO
2375     
2376           call global_min(MIN_VOL, GLOBAL_MIN_VOL,  PE_IO, ierr )
2377           call global_max(MAX_VOL, GLOBAL_MAX_VOL,  PE_IO, ierr )
2378           call global_min(MIN_AYZ, GLOBAL_MIN_AYZ,  PE_IO, ierr )
2379           call global_max(MAX_AYZ, GLOBAL_MAX_AYZ,  PE_IO, ierr )
2380           call global_min(MIN_AXZ, GLOBAL_MIN_AXZ,  PE_IO, ierr )
2381           call global_max(MAX_AXZ, GLOBAL_MAX_AXZ,  PE_IO, ierr )
2382           call global_min(MIN_AXY, GLOBAL_MIN_AXY,  PE_IO, ierr )
2383           call global_max(MAX_AXY, GLOBAL_MAX_AXY,  PE_IO, ierr )
2384           call global_min(MIN_CUT, GLOBAL_MIN_CUT,  PE_IO, ierr )
2385           call global_max(MAX_CUT, GLOBAL_MAX_CUT,  PE_IO, ierr )
2386     
2387           IF (myPE == PE_IO) THEN
2388              WRITE(UNIT_CUT_CELL_LOG,1000)  'U-MOMENTUM CUT CELLS:'
2389              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXY                   = ', GLOBAL_MIN_AXY,GLOBAL_MAX_AXY
2390              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXZ                   = ', GLOBAL_MIN_AXZ,GLOBAL_MAX_AXZ
2391              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AYZ                   = ', GLOBAL_MIN_AYZ,GLOBAL_MAX_AYZ
2392              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF CUT AREA              = ', GLOBAL_MIN_CUT,GLOBAL_MAX_CUT
2393              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF VOLUME                = ', GLOBAL_MIN_VOL,GLOBAL_MAX_VOL
2394              WRITE(UNIT_CUT_CELL_LOG,1010)  'NUMBER OF U WALL CELLS         = ', NUMBER_OF_U_WALL_CELLS
2395              WRITE(UNIT_CUT_CELL_LOG,1000)  '################################################################'
2396           ENDIF
2397     !======================================================================
2398     !  V-Momentum Cell volumes and areas
2399     !======================================================================
2400     
2401     
2402           MIN_VOL =   LARGE_NUMBER
2403           MAX_VOL = - LARGE_NUMBER
2404           MIN_AYZ =   LARGE_NUMBER
2405           MAX_AYZ = - LARGE_NUMBER
2406           MIN_AXZ =   LARGE_NUMBER
2407           MAX_AXZ = - LARGE_NUMBER
2408           MIN_AXY =   LARGE_NUMBER
2409           MAX_AXY = - LARGE_NUMBER
2410     
2411           DO IJK = IJKSTART3, IJKEND3
2412              IF(STANDARD_V_CELL_AT(IJK)) THEN              ! STANDARD CELLS
2413     
2414                 MIN_VOL =   DMIN1(MIN_VOL,VOL_V(IJK))
2415                 MAX_VOL =   DMAX1(MAX_VOL,VOL_V(IJK))
2416                 MIN_AYZ =   DMIN1(MIN_AYZ,AYZ_V(IJK))
2417                 MAX_AYZ =   DMAX1(MAX_AYZ,AYZ_V(IJK))
2418                 MIN_AXZ =   DMIN1(MIN_AXZ,AXZ_V(IJK))
2419                 MAX_AXZ =   DMAX1(MAX_AXZ,AXZ_V(IJK))
2420                 MIN_AXY =   DMIN1(MIN_AXY,AXY_V(IJK))
2421                 MAX_AXY =   DMAX1(MAX_AXY,AXY_V(IJK))
2422     
2423              ENDIF
2424           END DO
2425     
2426           call global_min(MIN_VOL, GLOBAL_MIN_VOL,  PE_IO, ierr )
2427           call global_max(MAX_VOL, GLOBAL_MAX_VOL,  PE_IO, ierr )
2428           call global_min(MIN_AYZ, GLOBAL_MIN_AYZ,  PE_IO, ierr )
2429           call global_max(MAX_AYZ, GLOBAL_MAX_AYZ,  PE_IO, ierr )
2430           call global_min(MIN_AXZ, GLOBAL_MIN_AXZ,  PE_IO, ierr )
2431           call global_max(MAX_AXZ, GLOBAL_MAX_AXZ,  PE_IO, ierr )
2432           call global_min(MIN_AXY, GLOBAL_MIN_AXY,  PE_IO, ierr )
2433           call global_max(MAX_AXY, GLOBAL_MAX_AXY,  PE_IO, ierr )
2434     
2435           IF (myPE == PE_IO) THEN
2436              WRITE(UNIT_CUT_CELL_LOG,1000)  'V-MOMENTUM STANDARD CELLS:'
2437              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXY                   = ', GLOBAL_MIN_AXY,GLOBAL_MAX_AXY
2438              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXZ                   = ', GLOBAL_MIN_AXZ,GLOBAL_MAX_AXZ
2439              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AYZ                   = ', GLOBAL_MIN_AYZ,GLOBAL_MAX_AYZ
2440              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF VOLUME                = ', GLOBAL_MIN_VOL,GLOBAL_MAX_VOL
2441           ENDIF
2442     
2443           MIN_VOL =   LARGE_NUMBER
2444           MAX_VOL = - LARGE_NUMBER
2445           MIN_AYZ =   LARGE_NUMBER
2446           MAX_AYZ = - LARGE_NUMBER
2447           MIN_AXZ =   LARGE_NUMBER
2448           MAX_AXZ = - LARGE_NUMBER
2449           MIN_AXY =   LARGE_NUMBER
2450           MAX_AXY = - LARGE_NUMBER
2451           MIN_CUT =   LARGE_NUMBER
2452           MAX_CUT = - LARGE_NUMBER
2453     
2454           DO IJK = IJKSTART3, IJKEND3
2455              IF(CUT_V_CELL_AT(IJK).AND.(.NOT.WALL_V_AT(IJK))) THEN      ! CUT CELLS
2456     
2457                 MIN_VOL =   DMIN1(MIN_VOL,VOL_V(IJK))
2458                 MAX_VOL =   DMAX1(MAX_VOL,VOL_V(IJK))
2459                 MIN_AYZ =   DMIN1(MIN_AYZ,AYZ_V(IJK))
2460                 MAX_AYZ =   DMAX1(MAX_AYZ,AYZ_V(IJK))
2461                 MIN_AXZ =   DMIN1(MIN_AXZ,AXZ_V(IJK))
2462                 MAX_AXZ =   DMAX1(MAX_AXZ,AXZ_V(IJK))
2463                 MIN_AXY =   DMIN1(MIN_AXY,AXY_V(IJK))
2464                 MAX_AXY =   DMAX1(MAX_AXY,AXY_V(IJK))
2465                 MIN_CUT =   DMIN1(MIN_CUT,AREA_V_CUT(IJK))
2466                 MAX_CUT =   DMAX1(MAX_CUT,AREA_V_CUT(IJK))
2467     
2468              ENDIF
2469           END DO
2470     
2471           call global_min(MIN_VOL, GLOBAL_MIN_VOL,  PE_IO, ierr )
2472           call global_max(MAX_VOL, GLOBAL_MAX_VOL,  PE_IO, ierr )
2473           call global_min(MIN_AYZ, GLOBAL_MIN_AYZ,  PE_IO, ierr )
2474           call global_max(MAX_AYZ, GLOBAL_MAX_AYZ,  PE_IO, ierr )
2475           call global_min(MIN_AXZ, GLOBAL_MIN_AXZ,  PE_IO, ierr )
2476           call global_max(MAX_AXZ, GLOBAL_MAX_AXZ,  PE_IO, ierr )
2477           call global_min(MIN_AXY, GLOBAL_MIN_AXY,  PE_IO, ierr )
2478           call global_max(MAX_AXY, GLOBAL_MAX_AXY,  PE_IO, ierr )
2479           call global_min(MIN_CUT, GLOBAL_MIN_CUT,  PE_IO, ierr )
2480           call global_max(MAX_CUT, GLOBAL_MAX_CUT,  PE_IO, ierr )
2481     
2482           IF (myPE == PE_IO) THEN
2483              WRITE(UNIT_CUT_CELL_LOG,1000)  'V-MOMENTUM CUT CELLS:'
2484              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXY                   = ', GLOBAL_MIN_AXY,GLOBAL_MAX_AXY
2485              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXZ                   = ', GLOBAL_MIN_AXZ,GLOBAL_MAX_AXZ
2486              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AYZ                   = ', GLOBAL_MIN_AYZ,GLOBAL_MAX_AYZ
2487              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF CUT AREA              = ', GLOBAL_MIN_CUT,GLOBAL_MAX_CUT
2488              WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF VOLUME                = ', GLOBAL_MIN_VOL,GLOBAL_MAX_VOL
2489              WRITE(UNIT_CUT_CELL_LOG,1010)  'NUMBER OF V WALL CELLS         = ', NUMBER_OF_V_WALL_CELLS
2490              WRITE(UNIT_CUT_CELL_LOG,1000)  '################################################################'
2491           ENDIF
2492     
2493     !======================================================================
2494     !  W-Momentum Cell volumes and areas
2495     !======================================================================
2496     
2497     
2498           IF(DO_K) THEN
2499     
2500              MIN_VOL =   LARGE_NUMBER
2501              MAX_VOL = - LARGE_NUMBER
2502              MIN_AYZ =   LARGE_NUMBER
2503              MAX_AYZ = - LARGE_NUMBER
2504              MIN_AXZ =   LARGE_NUMBER
2505              MAX_AXZ = - LARGE_NUMBER
2506              MIN_AXY =   LARGE_NUMBER
2507              MAX_AXY = - LARGE_NUMBER
2508     
2509              DO IJK = IJKSTART3, IJKEND3
2510                 IF(STANDARD_W_CELL_AT(IJK)) THEN              ! STANDARD CELLS
2511     
2512                    MIN_VOL =   DMIN1(MIN_VOL,VOL_W(IJK))
2513                    MAX_VOL =   DMAX1(MAX_VOL,VOL_W(IJK))
2514                    MIN_AYZ =   DMIN1(MIN_AYZ,AYZ_W(IJK))
2515                    MAX_AYZ =   DMAX1(MAX_AYZ,AYZ_W(IJK))
2516                    MIN_AXZ =   DMIN1(MIN_AXZ,AXZ_W(IJK))
2517                    MAX_AXZ =   DMAX1(MAX_AXZ,AXZ_W(IJK))
2518                    MIN_AXY =   DMIN1(MIN_AXY,AXY_W(IJK))
2519                    MAX_AXY =   DMAX1(MAX_AXY,AXY_W(IJK))
2520     
2521                 ENDIF
2522              END DO
2523     
2524              call global_min(MIN_VOL, GLOBAL_MIN_VOL,  PE_IO, ierr )
2525              call global_max(MAX_VOL, GLOBAL_MAX_VOL,  PE_IO, ierr )
2526              call global_min(MIN_AYZ, GLOBAL_MIN_AYZ,  PE_IO, ierr )
2527              call global_max(MAX_AYZ, GLOBAL_MAX_AYZ,  PE_IO, ierr )
2528              call global_min(MIN_AXZ, GLOBAL_MIN_AXZ,  PE_IO, ierr )
2529              call global_max(MAX_AXZ, GLOBAL_MAX_AXZ,  PE_IO, ierr )
2530              call global_min(MIN_AXY, GLOBAL_MIN_AXY,  PE_IO, ierr )
2531              call global_max(MAX_AXY, GLOBAL_MAX_AXY,  PE_IO, ierr )
2532     
2533              IF (myPE == PE_IO) THEN
2534                 WRITE(UNIT_CUT_CELL_LOG,1000)  'W-MOMENTUM STANDARD CELLS:'
2535                 WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXY                   = ', GLOBAL_MIN_AXY,GLOBAL_MAX_AXY
2536                 WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXZ                   = ', GLOBAL_MIN_AXZ,GLOBAL_MAX_AXZ
2537                 WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AYZ                   = ', GLOBAL_MIN_AYZ,GLOBAL_MAX_AYZ
2538                 WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF VOLUME                = ', GLOBAL_MIN_VOL,GLOBAL_MAX_VOL
2539              ENDIF
2540     
2541              MIN_VOL =   LARGE_NUMBER
2542              MAX_VOL = - LARGE_NUMBER
2543              MIN_AYZ =   LARGE_NUMBER
2544              MAX_AYZ = - LARGE_NUMBER
2545              MIN_AXZ =   LARGE_NUMBER
2546              MAX_AXZ = - LARGE_NUMBER
2547              MIN_AXY =   LARGE_NUMBER
2548              MAX_AXY = - LARGE_NUMBER
2549              MIN_CUT =   LARGE_NUMBER
2550              MAX_CUT = - LARGE_NUMBER
2551     
2552              DO IJK = IJKSTART3, IJKEND3
2553                 IF(CUT_W_CELL_AT(IJK).AND.(.NOT.WALL_W_AT(IJK))) THEN      ! CUT CELLS
2554     
2555                    MIN_VOL =   DMIN1(MIN_VOL,VOL_W(IJK))
2556                    MAX_VOL =   DMAX1(MAX_VOL,VOL_W(IJK))
2557                    MIN_AYZ =   DMIN1(MIN_AYZ,AYZ_W(IJK))
2558                    MAX_AYZ =   DMAX1(MAX_AYZ,AYZ_W(IJK))
2559                    MIN_AXZ =   DMIN1(MIN_AXZ,AXZ_W(IJK))
2560                    MAX_AXZ =   DMAX1(MAX_AXZ,AXZ_W(IJK))
2561                    MIN_AXY =   DMIN1(MIN_AXY,AXY_W(IJK))
2562                    MAX_AXY =   DMAX1(MAX_AXY,AXY_W(IJK))
2563                    MIN_CUT =   DMIN1(MIN_CUT,AREA_W_CUT(IJK))
2564                    MAX_CUT =   DMAX1(MAX_CUT,AREA_W_CUT(IJK))
2565     
2566                 ENDIF
2567              END DO
2568     
2569              call global_min(MIN_VOL, GLOBAL_MIN_VOL,  PE_IO, ierr )
2570              call global_max(MAX_VOL, GLOBAL_MAX_VOL,  PE_IO, ierr )
2571              call global_min(MIN_AYZ, GLOBAL_MIN_AYZ,  PE_IO, ierr )
2572              call global_max(MAX_AYZ, GLOBAL_MAX_AYZ,  PE_IO, ierr )
2573              call global_min(MIN_AXZ, GLOBAL_MIN_AXZ,  PE_IO, ierr )
2574              call global_max(MAX_AXZ, GLOBAL_MAX_AXZ,  PE_IO, ierr )
2575              call global_min(MIN_AXY, GLOBAL_MIN_AXY,  PE_IO, ierr )
2576              call global_max(MAX_AXY, GLOBAL_MAX_AXY,  PE_IO, ierr )
2577              call global_min(MIN_CUT, GLOBAL_MIN_CUT,  PE_IO, ierr )
2578              call global_max(MAX_CUT, GLOBAL_MAX_CUT,  PE_IO, ierr )
2579     
2580              IF (myPE == PE_IO) THEN
2581                 WRITE(UNIT_CUT_CELL_LOG,1000)  'W-MOMENTUM CUT CELLS:'
2582                 WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXY                   = ', GLOBAL_MIN_AXY,GLOBAL_MAX_AXY
2583                 WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AXZ                   = ', GLOBAL_MIN_AXZ,GLOBAL_MAX_AXZ
2584                 WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF AYZ                   = ', GLOBAL_MIN_AYZ,GLOBAL_MAX_AYZ
2585                 WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF CUT AREA              = ', GLOBAL_MIN_CUT,GLOBAL_MAX_CUT
2586                 WRITE(UNIT_CUT_CELL_LOG,1000)  'RANGE OF VOLUME                = ', GLOBAL_MIN_VOL,GLOBAL_MAX_VOL
2587                 WRITE(UNIT_CUT_CELL_LOG,1010)  'NUMBER OF W WALL CELLS         = ', NUMBER_OF_W_WALL_CELLS
2588                 WRITE(UNIT_CUT_CELL_LOG,1000)  '################################################################'
2589              ENDIF
2590     
2591           ENDIF
2592     
2593     
2594     
2595           LOCAL_MIN_Q = MINVAL(Alpha_Ue_c)
2596           LOCAL_MAX_Q = MAXVAL(Alpha_Ue_c)
2597           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2598           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2599           IF (myPE == PE_IO)  WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF U-MOMENTUM Alpha_Ue_c = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2600     
2601           LOCAL_MIN_Q = MINVAL(Alpha_Un_c)
2602           LOCAL_MAX_Q = MAXVAL(Alpha_Un_c)
2603           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2604           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2605           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF U-MOMENTUM Alpha_Un_c = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2606     
2607           LOCAL_MIN_Q = MINVAL(Alpha_Ut_c)
2608           LOCAL_MAX_Q = MAXVAL(Alpha_Ut_c)
2609           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2610           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2611           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF U-MOMENTUM Alpha_Ut_c = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2612           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)
2613     
2614           LOCAL_MIN_Q = MINVAL(Theta_Ue)
2615           LOCAL_MAX_Q = MAXVAL(Theta_Ue)
2616           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2617           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2618           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF U-MOMENTUM Theta_Ue   = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2619     
2620           LOCAL_MIN_Q = MINVAL(Theta_Un)
2621           LOCAL_MAX_Q = MAXVAL(Theta_Un)
2622           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2623           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2624           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF U-MOMENTUM Theta_Un   = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2625     
2626           LOCAL_MIN_Q = MINVAL(Theta_Ut)
2627           LOCAL_MAX_Q = MAXVAL(Theta_Ut)
2628           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2629           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2630           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF U-MOMENTUM Theta_Ut   = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2631           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)
2632     
2633           LOCAL_MIN_Q = MINVAL(Theta_U_ne)
2634           LOCAL_MAX_Q = MAXVAL(Theta_U_ne)
2635           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2636           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2637           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF U-MOMENTUM Theta_U_ne = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2638     
2639           LOCAL_MIN_Q = MINVAL(Theta_U_te)
2640           LOCAL_MAX_Q = MAXVAL(Theta_U_te)
2641           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2642           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2643           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF U-MOMENTUM Theta_U_te = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2644           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)
2645     
2646           LOCAL_MIN_Q = MINVAL(NOC_U_E)
2647           LOCAL_MAX_Q = MAXVAL(NOC_U_E)
2648           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2649           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2650           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF U-MOMENTUM NOC_U_E    = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2651     
2652           LOCAL_MIN_Q = MINVAL(NOC_U_N)
2653           LOCAL_MAX_Q = MAXVAL(NOC_U_N)
2654           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2655           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2656           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF U-MOMENTUM NOC_U_N    = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2657     
2658           LOCAL_MIN_Q = MINVAL(NOC_U_T)
2659           LOCAL_MAX_Q = MAXVAL(NOC_U_T)
2660           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2661           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2662           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF U-MOMENTUM NOC_U_T    = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2663           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)
2664     
2665           LOCAL_MIN_Q = MINVAL(DELH_U)
2666           LOCAL_MAX_Q = MAXVAL(DELH_U)
2667           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2668           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2669           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF U-MOMENTUM DELH_U     = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2670           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  '################################################################'
2671     
2672     
2673     
2674           LOCAL_MIN_Q = MINVAL(Alpha_Ve_c)
2675           LOCAL_MAX_Q = MAXVAL(Alpha_Ve_c)
2676           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2677           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2678           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF V-MOMENTUM Alpha_Ve_c = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2679           LOCAL_MIN_Q = MINVAL(Alpha_Vn_c)
2680           LOCAL_MAX_Q = MAXVAL(Alpha_Vn_c)
2681           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2682           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2683           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF V-MOMENTUM Alpha_Vn_c = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2684           LOCAL_MIN_Q = MINVAL(Alpha_Vt_c)
2685           LOCAL_MAX_Q = MAXVAL(Alpha_Vt_c)
2686           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2687           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2688           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF V-MOMENTUM Alpha_Vt_c = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2689           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)
2690           LOCAL_MIN_Q = MINVAL(Theta_Ve)
2691           LOCAL_MAX_Q = MAXVAL(Theta_Ve)
2692           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2693           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2694           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF V-MOMENTUM Theta_Ve   = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2695           LOCAL_MIN_Q = MINVAL(Theta_Vn)
2696           LOCAL_MAX_Q = MAXVAL(Theta_Vn)
2697           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2698           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2699           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF V-MOMENTUM Theta_Vn   = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2700           LOCAL_MIN_Q = MINVAL(Theta_Vt)
2701           LOCAL_MAX_Q = MAXVAL(Theta_Vt)
2702           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2703           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2704           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF V-MOMENTUM Theta_Vt   = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2705           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)
2706           LOCAL_MIN_Q = MINVAL(Theta_V_ne)
2707           LOCAL_MAX_Q = MAXVAL(Theta_V_ne)
2708           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2709           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2710           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF V-MOMENTUM Theta_V_ne = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2711           LOCAL_MIN_Q = MINVAL(Theta_V_nt)
2712           LOCAL_MAX_Q = MAXVAL(Theta_V_nt)
2713           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2714           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2715           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF V-MOMENTUM Theta_V_nt = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2716           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)
2717           LOCAL_MIN_Q = MINVAL(NOC_V_E)
2718           LOCAL_MAX_Q = MAXVAL(NOC_V_E)
2719           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2720           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2721           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF V-MOMENTUM NOC_V_E    = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2722           LOCAL_MIN_Q = MINVAL(NOC_V_N)
2723           LOCAL_MAX_Q = MAXVAL(NOC_V_N)
2724           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2725           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2726           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF V-MOMENTUM NOC_V_N    = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2727           LOCAL_MIN_Q = MINVAL(NOC_V_T)
2728           LOCAL_MAX_Q = MAXVAL(NOC_V_T)
2729           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2730           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2731           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF V-MOMENTUM NOC_V_T    = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2732           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)
2733           LOCAL_MIN_Q = MINVAL(DELH_V)
2734           LOCAL_MAX_Q = MAXVAL(DELH_V)
2735           call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2736           call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2737           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF V-MOMENTUM DELH_V     = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2738           IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  '################################################################'
2739     
2740     
2741           IF(DO_K) THEN
2742     
2743              LOCAL_MIN_Q = MINVAL(Alpha_We_c)
2744              LOCAL_MAX_Q = MAXVAL(Alpha_We_c)
2745              call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2746              call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2747              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF W-MOMENTUM Alpha_We_c = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2748              LOCAL_MIN_Q = MINVAL(Alpha_Wn_c)
2749              LOCAL_MAX_Q = MAXVAL(Alpha_Wn_c)
2750              call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2751              call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2752              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF W-MOMENTUM Alpha_Wn_c = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2753              LOCAL_MIN_Q = MINVAL(Alpha_Wt_c)
2754              LOCAL_MAX_Q = MAXVAL(Alpha_Wt_c)
2755              call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2756              call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2757              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF W-MOMENTUM Alpha_Wt_c = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2758              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)
2759              LOCAL_MIN_Q = MINVAL(Theta_We)
2760              LOCAL_MAX_Q = MAXVAL(Theta_We)
2761              call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2762              call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2763              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF W-MOMENTUM Theta_We   = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2764              LOCAL_MIN_Q = MINVAL(Theta_Wn)
2765              LOCAL_MAX_Q = MAXVAL(Theta_Wn)
2766              call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2767              call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2768              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF W-MOMENTUM Theta_Wn   = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2769              LOCAL_MIN_Q = MINVAL(Theta_Wt)
2770              LOCAL_MAX_Q = MAXVAL(Theta_Wt)
2771              call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2772              call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2773              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF W-MOMENTUM Theta_Wt   = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2774              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)
2775              LOCAL_MIN_Q = MINVAL(Theta_W_te)
2776              LOCAL_MAX_Q = MAXVAL(Theta_W_te)
2777              call global_min(LOCAL_MAX_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2778              call global_max(LOCAL_MIN_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2779              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF W-MOMENTUM Theta_W_te = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2780              LOCAL_MIN_Q = MINVAL(Theta_W_tn)
2781              LOCAL_MAX_Q = MAXVAL(Theta_W_tn)
2782              call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2783              call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2784              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF W-MOMENTUM Theta_W_tn = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2785              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)
2786              LOCAL_MIN_Q = MINVAL(NOC_W_E)
2787              LOCAL_MAX_Q = MAXVAL(NOC_W_E)
2788              call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2789              call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2790              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF W-MOMENTUM NOC_W_E    = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2791              LOCAL_MIN_Q = MINVAL(NOC_W_N)
2792              LOCAL_MAX_Q = MAXVAL(NOC_W_N)
2793              call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2794              call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2795              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF W-MOMENTUM NOC_W_N    = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2796              LOCAL_MIN_Q = MINVAL(NOC_W_T)
2797              LOCAL_MAX_Q = MAXVAL(NOC_W_T)
2798              call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2799              call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2800              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF W-MOMENTUM NOC_W_T    = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2801              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)
2802              LOCAL_MIN_Q = MINVAL(DELH_W)
2803              LOCAL_MAX_Q = MAXVAL(DELH_W)
2804              call global_min(LOCAL_MIN_Q, GLOBAL_MIN_Q,  PE_IO, ierr )
2805              call global_max(LOCAL_MAX_Q, GLOBAL_MAX_Q,  PE_IO, ierr )
2806              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  ' RANGE OF W-MOMENTUM DELH_W     = ', GLOBAL_MIN_Q, GLOBAL_MAX_Q
2807              IF (myPE == PE_IO) WRITE(UNIT_CUT_CELL_LOG,1000)  '################################################################'
2808     
2809           ENDIF
2810     
2811           RETURN
2812     
2813           END SUBROUTINE PRINT_GRID_STATISTICS
2814     
2815     
2816     
2817     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2818     !                                                                      C
2819     !  Module name: CLEAN_GEOMETRY                                         C
2820     !  Purpose: Clean-up the list of point and only keep points            C
2821     !           that are used in the connectivity list.                    C
2822     !                                                                      C
2823     !  Author: Jeff Dietiker                              Date: 19-Dec-14  C
2824     !  Reviewer:                                          Date:            C
2825     !                                                                      C
2826     !  Revision Number #                                  Date: ##-###-##  C
2827     !  Author: #                                                           C
2828     !  Purpose: #                                                          C
2829     !                                                                      C
2830     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
2831     
2832           SUBROUTINE CLEAN_GEOMETRY
2833     
2834           USE param
2835           USE param1
2836           USE parallel
2837           USE constant
2838           USE run
2839           USE toleranc
2840           USE geometry
2841           USE indices
2842           USE compar
2843           USE mpi_utility
2844           USE sendrecv
2845           USE quadric
2846           USE cutcell
2847           USE fldvar
2848           USE vtk
2849           USE cdist
2850           USE functions
2851     
2852           IMPLICIT NONE
2853     
2854           INTEGER :: IJK,L
2855     
2856           INTEGER ::POINT_ID,IJKC
2857           INTEGER , ALLOCATABLE        ::  POINT_NEW_ID(:)
2858           INTEGER , ALLOCATABLE        ::  NEW_POINT_NEW_ID(:)
2859           LOGICAL , ALLOCATABLE        ::  KEEP_POINT(:)
2860           LOGICAL , ALLOCATABLE        ::  KEEP_NEW_POINT(:)
2861     
2862     
2863           IF (myPE == PE_IO.AND.(.NOT.BDIST_IO)) THEN
2864     
2865              IF(ALLOCATED(GLOBAL_CLEANED_CONNECTIVITY)) DEALLOCATE(GLOBAL_CLEANED_CONNECTIVITY)
2866              IF(ALLOCATED(KEEP_NEW_POINT))              DEALLOCATE (KEEP_NEW_POINT)
2867              IF(ALLOCATED(POINT_NEW_ID))                DEALLOCATE (POINT_NEW_ID)
2868              IF(ALLOCATED(NEW_POINT_NEW_ID))            DEALLOCATE (NEW_POINT_NEW_ID)
2869              IF(ALLOCATED(KEEP_POINT))                  DEALLOCATE (KEEP_POINT)
2870     
2871              ALLOCATE (GLOBAL_CLEANED_CONNECTIVITY(ijkmax3,15))
2872              ALLOCATE (KEEP_NEW_POINT(GLOBAL_NUMBER_OF_NEW_POINTS))
2873     
2874              ALLOCATE (POINT_NEW_ID(IJKMAX3))
2875              ALLOCATE (NEW_POINT_NEW_ID(IJKMAX3))
2876              ALLOCATE (KEEP_POINT(IJKMAX3))
2877     
2878     ! Step 1: Go through connectivity list and only keep points that are used.
2879     !         For background cell corners, assign KEEP_POINT = .TRUE.
2880     !         For cut cells, the new intersection points were called NEW_POINTS,
2881     !         so assign KEEP_NEW_POINT = .TRUE.
2882     !         A NEW_POINT had an IJK index larger than IJKMAX3
2883     
2884              KEEP_POINT = .FALSE.
2885              KEEP_NEW_POINT = .FALSE.
2886     
2887              DO IJK = 1,IJKMAX3
2888                 IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) THEN
2889                    DO L=1,GLOBAL_NUMBER_OF_NODES(IJK)
2890                       IJKC = GLOBAL_CONNECTIVITY(IJK,L)
2891                       IF(IJKC<=IJKMAX3) KEEP_POINT(IJKC) = .TRUE.
2892                       IF(IJKC>IJKMAX3) KEEP_NEW_POINT(IJKC-IJKMAX3) = .TRUE.
2893                    ENDDO
2894                 ENDIF
2895              END DO
2896     
2897     
2898     ! Step 2: Clean-up list of used points and store cleaned connectivity
2899              POINT_NEW_ID = -1
2900              NEW_POINT_NEW_ID = -1
2901              POINT_ID = 1
2902     ! This is for the background grid cell corners
2903              DO IJK = 1,IJKMAX3
2904                 IF(KEEP_POINT(IJK)) THEN
2905                    POINT_NEW_ID(IJK) = POINT_ID
2906                    POINT_ID = POINT_ID + 1
2907                 ENDIF
2908              END DO
2909     ! This is for the cut cell new corners
2910              DO IJK = 1,GLOBAL_NUMBER_OF_NEW_POINTS
2911                 IF(KEEP_NEW_POINT(IJK)) THEN
2912                    NEW_POINT_NEW_ID(IJK) = POINT_ID
2913                    POINT_ID = POINT_ID + 1
2914                 ENDIF
2915              END DO
2916     
2917     ! Update the true (clean) number of points
2918              NUMBER_OF_POINTS = POINT_ID - 1
2919     
2920     ! Now, store a list of coordinates for all used points
2921              IF(ALLOCATED(GLOBAL_COORDS_OF_POINTS)) DEALLOCATE(GLOBAL_COORDS_OF_POINTS)
2922     
2923              ALLOCATE(GLOBAL_COORDS_OF_POINTS(3,NUMBER_OF_POINTS))
2924     
2925              POINT_ID = 1
2926     ! This is for the background grid cell corners
2927              DO IJK = 1,IJKMAX3
2928                 IF(KEEP_POINT(IJK)) THEN
2929                    GLOBAL_COORDS_OF_POINTS(1:3,POINT_ID) = &
2930                         (/REAL(XG_E(GLOBAL_I_OF(IJK))),REAL(YG_N(GLOBAL_J_OF(IJK))),REAL(ZG_T(GLOBAL_K_OF(IJK)))/)
2931                    POINT_ID = POINT_ID + 1
2932                 ENDIF
2933              END DO
2934     ! This is for the cut cell new corners
2935              DO IJK = 1,GLOBAL_NUMBER_OF_NEW_POINTS
2936                 IF(KEEP_NEW_POINT(IJK)) THEN
2937                    NEW_POINT_NEW_ID(IJK) = POINT_ID
2938                    GLOBAL_COORDS_OF_POINTS(1:3,POINT_ID) = &
2939                         (/REAL(GLOBAL_X_NEW_POINT(IJK)),REAL(GLOBAL_Y_NEW_POINT(IJK)),REAL(GLOBAL_Z_NEW_POINT(IJK))/)
2940                    POINT_ID = POINT_ID + 1
2941                 ENDIF
2942              END DO
2943     
2944     
2945     ! Step 3: Shift connectivity with new point indices
2946              DO IJK = 1,IJKMAX3
2947                 IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) THEN
2948                    DO L=1,GLOBAL_NUMBER_OF_NODES(IJK)
2949                       IF(GLOBAL_CONNECTIVITY(IJK,L)<=IJKMAX3) THEN
2950                          GLOBAL_CLEANED_CONNECTIVITY(IJK,L) = POINT_NEW_ID(GLOBAL_CONNECTIVITY(IJK,L))
2951                       ELSE
2952                          GLOBAL_CLEANED_CONNECTIVITY(IJK,L) = NEW_POINT_NEW_ID(GLOBAL_CONNECTIVITY(IJK,L)-IJKMAX3)
2953                       ENDIF
2954                    ENDDO
2955                 ENDIF
2956              END DO
2957     
2958     
2959     
2960            ELSEIF(BDIST_IO) THEN
2961     
2962     
2963               IF(ALLOCATED(CLEANED_CONNECTIVITY))  DEALLOCATE (CLEANED_CONNECTIVITY)
2964               IF(ALLOCATED(KEEP_NEW_POINT))        DEALLOCATE (KEEP_NEW_POINT)
2965               IF(ALLOCATED(POINT_NEW_ID))          DEALLOCATE (POINT_NEW_ID)
2966               IF(ALLOCATED(NEW_POINT_NEW_ID))      DEALLOCATE (NEW_POINT_NEW_ID)
2967               IF(ALLOCATED(KEEP_POINT))            DEALLOCATE (KEEP_POINT)
2968     
2969               ALLOCATE (CLEANED_CONNECTIVITY(IJKEND3,15))
2970               ALLOCATE (KEEP_NEW_POINT(NUMBER_OF_NEW_POINTS))
2971     
2972               ALLOCATE (POINT_NEW_ID(IJKEND3))
2973               ALLOCATE (NEW_POINT_NEW_ID(IJKEND3))
2974               ALLOCATE (KEEP_POINT(IJKEND3))
2975     
2976     ! Step 1: Go through connectivity list and only keep points that are used.
2977     !         For background cell corners, assign KEEP_POINT = .TRUE.
2978     !         For cut cells, the new intersection points were called NEW_POINTS,
2979     !         so assign KEEP_NEW_POINT = .TRUE.
2980     !         A NEW_POINT had an IJK index larger than IJKMAX3
2981     
2982               KEEP_POINT = .FALSE.
2983               KEEP_NEW_POINT = .FALSE.
2984     
2985               DO IJK = 1,IJKEND3
2986                  IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) THEN
2987                     DO L=1,NUMBER_OF_NODES(IJK)
2988                        IJKC = CONNECTIVITY(IJK,L)
2989                        IF(IJKC<=IJKEND3) KEEP_POINT(IJKC) = .TRUE.
2990                        IF(IJKC>IJKEND3) KEEP_NEW_POINT(IJKC-IJKEND3) = .TRUE.
2991                     ENDDO
2992                  ENDIF
2993               END DO
2994     
2995     
2996     ! Step 2: Clean-up list of used points and store cleaned connectivity
2997               POINT_NEW_ID = -1
2998               NEW_POINT_NEW_ID = -1
2999               POINT_ID = 1
3000     ! This is for the background grid cell corners
3001               DO IJK = 1,IJKEND3
3002                  IF(KEEP_POINT(IJK)) THEN
3003                     POINT_NEW_ID(IJK) = POINT_ID
3004                     POINT_ID = POINT_ID + 1
3005                  ENDIF
3006               END DO
3007     ! This is for the cut cell new corners
3008               DO IJK = 1,NUMBER_OF_NEW_POINTS
3009                  IF(KEEP_NEW_POINT(IJK)) THEN
3010                     NEW_POINT_NEW_ID(IJK) = POINT_ID
3011                     POINT_ID = POINT_ID + 1
3012                  ENDIF
3013               END DO
3014     
3015     ! Update the true (clean) number of points
3016               NUMBER_OF_POINTS = POINT_ID - 1
3017     
3018     ! Now, store a list of coordinates for all used points
3019               IF(ALLOCATED(COORDS_OF_POINTS)) DEALLOCATE(COORDS_OF_POINTS)
3020     
3021               ALLOCATE(COORDS_OF_POINTS(NUMBER_OF_POINTS,3))
3022     
3023               POINT_ID = 1
3024     ! This is for the background grid cell corners
3025               DO IJK = 1,IJKEND3
3026                  IF(KEEP_POINT(IJK)) THEN
3027                     COORDS_OF_POINTS(POINT_ID,1:3) = &
3028                          (/REAL(XG_E(I_OF(IJK))),REAL(YG_N(J_OF(IJK))),REAL(ZG_T(K_OF(IJK)))/)
3029                     POINT_ID = POINT_ID + 1
3030                  ENDIF
3031               END DO
3032     ! This is for the cut cell new corners
3033               DO IJK = 1,NUMBER_OF_NEW_POINTS
3034                  IF(KEEP_NEW_POINT(IJK)) THEN
3035                     NEW_POINT_NEW_ID(IJK) = POINT_ID
3036                     COORDS_OF_POINTS(POINT_ID,1:3) = &
3037                          (/REAL(X_NEW_POINT(IJK)),REAL(Y_NEW_POINT(IJK)),REAL(Z_NEW_POINT(IJK))/)
3038                     POINT_ID = POINT_ID + 1
3039                  ENDIF
3040               END DO
3041     
3042     
3043     ! Step 3: Shift connectivity with new point indices
3044               DO IJK = 1,IJKEND3
3045                  IF (BELONGS_TO_VTK_SUBDOMAIN(IJK)) THEN
3046                     DO L=1,NUMBER_OF_NODES(IJK)
3047                        IF(CONNECTIVITY(IJK,L)<=IJKEND3) THEN
3048                           CLEANED_CONNECTIVITY(IJK,L) = POINT_NEW_ID(CONNECTIVITY(IJK,L))
3049                        ELSE
3050                           CLEANED_CONNECTIVITY(IJK,L) = NEW_POINT_NEW_ID(CONNECTIVITY(IJK,L)-IJKEND3)
3051                        ENDIF
3052                     ENDDO
3053                  ENDIF
3054               END DO
3055     
3056            ENDIF
3057     
3058           RETURN
3059     
3060           END SUBROUTINE CLEAN_GEOMETRY
3061     
3062     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
3063     !                                                                      C
3064     !  Module name: SETUP_VTK_REGION                                       C
3065     !                                                                      C
3066     !  Purpose: Filter the cells based on the VTK region bounds and        C
3067     !           set the flag BELONGS_TO_VTK_SUBDOMAIN(IJK) to .TRUE.       C
3068     !           to keep the cell.                                          C
3069     !                                                                      C
3070     !  Author: Jeff Dietiker                              Date: 19-Dec-14  C
3071     !  Reviewer:                                          Date:            C
3072     !                                                                      C
3073     !  Revision Number #                                  Date: ##-###-##  C
3074     !  Author: #                                                           C
3075     !  Purpose: #                                                          C
3076     !                                                                      C
3077     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
3078           SUBROUTINE SETUP_VTK_REGION
3079     
3080           USE cdist
3081           USE compar, only: mype, pe_io, ijkend3
3082           USE cutcell
3083           USE geometry
3084           USE indices, only: i_of, j_of, k_of
3085           USE vtk
3086     
3087           IMPLICIT NONE
3088     
3089           INTEGER :: IJK,I,J,K,I_E,I_W,J_N,J_S,K_T,K_B
3090           INTEGER :: NXS,NYS,NZS,NS,I_TMP,J_TMP,K_TMP
3091           INTEGER :: I_SLICE(DIM_I),J_SLICE(DIM_J),K_SLICE(DIM_K)
3092           DOUBLE PRECISION :: XE,XW,YS,YN,ZB,ZT
3093           DOUBLE PRECISION :: XSLICE,YSLICE,ZSLICE
3094           LOGICAL :: KEEP_XDIR,KEEP_YDIR,KEEP_ZDIR
3095     
3096     ! Get VTK region bounds
3097           XE = VTK_X_E(VTK_REGION)
3098           XW = VTK_X_W(VTK_REGION)
3099           YS = VTK_Y_S(VTK_REGION)
3100           YN = VTK_Y_N(VTK_REGION)
3101           ZB = VTK_Z_B(VTK_REGION)
3102           ZT = VTK_Z_T(VTK_REGION)
3103     
3104           NXS = VTK_NXS(VTK_REGION)
3105           NYS = VTK_NYS(VTK_REGION)
3106           NZS = VTK_NZS(VTK_REGION)
3107     
3108           CALL CALC_CELL (XMIN, VTK_X_W(VTK_REGION), DX, IMAX, I_W)
3109           I_W = I_W !+ 1
3110           CALL CALC_CELL (XMIN, VTK_X_E(VTK_REGION), DX, IMAX, I_E)
3111     
3112     
3113           CALL CALC_CELL (ZERO, VTK_Y_S(VTK_REGION), DY, JMAX, J_S)
3114           J_S = J_S !+ 1
3115           CALL CALC_CELL (ZERO, VTK_Y_N(VTK_REGION), DY, JMAX, J_N)
3116     
3117           IF (NO_K) THEN
3118              K_B = 1
3119              K_T = 1
3120           ELSE
3121              CALL CALC_CELL (ZERO, VTK_Z_B(VTK_REGION), DZ, KMAX, K_B)
3122              K_B = K_B !+ 1
3123              CALL CALC_CELL (ZERO, VTK_Z_T(VTK_REGION), DZ, KMAX, K_T)
3124           ENDIF
3125     
3126     ! get slice(s) location
3127           DO NS = 1,NXS
3128              XSLICE = XW + (XE-XW)/(NXS-1)*(NS-1)
3129              CALL CALC_CELL (XMIN, XSLICE, DX, IMAX, I_TMP)
3130              I_SLICE(NS) = MAX(MIN(I_TMP,IMAX1),IMIN1)
3131           ENDDO
3132     
3133           DO NS = 1,NYS
3134              YSLICE = YS + (YN-YS)/(NYS-1)*(NS-1)
3135              CALL CALC_CELL (ZERO, YSLICE, DY, JMAX, J_TMP)
3136              J_SLICE(NS) = MAX(MIN(J_TMP,JMAX1),JMIN1)
3137           ENDDO
3138     
3139           DO NS = 1,NZS
3140              ZSLICE = ZB + (ZT-ZB)/(NZS-1)*(NS-1)
3141              CALL CALC_CELL (ZERO, ZSLICE, DZ, KMAX, K_TMP)
3142              K_SLICE(NS) = MAX(MIN(K_TMP,KMAX1),KMIN1)
3143           ENDDO
3144     
3145           IF (myPE == PE_IO.AND.(.NOT.BDIST_IO)) THEN
3146     
3147              IF(ALLOCATED(BELONGS_TO_VTK_SUBDOMAIN)) DEALLOCATE(BELONGS_TO_VTK_SUBDOMAIN)
3148     
3149              ALLOCATE (BELONGS_TO_VTK_SUBDOMAIN(ijkmax3))
3150     
3151     ! Filter the cells based on the VTK region bounds and set the
3152     ! flag BELONGS_TO_VTK_SUBDOMAIN(IJK) to .TRUE. to keep the cell.
3153     
3154              BELONGS_TO_VTK_SUBDOMAIN = .FALSE.
3155              NUMBER_OF_VTK_CELLS      = 0
3156     
3157              DO IJK = 1,IJKMAX3
3158                 IF (GLOBAL_INTERIOR_CELL_AT(IJK))      THEN
3159                    IF (.NOT.GLOBAL_BLOCKED_CELL_AT(IJK)) THEN
3160                       I = GLOBAL_I_OF(IJK)
3161                       J = GLOBAL_J_OF(IJK)
3162                       K = GLOBAL_K_OF(IJK)
3163     
3164     ! X-direction
3165                       KEEP_XDIR=.FALSE.
3166                       IF(NXS==0) THEN
3167                          IF(I_W<=I.AND.I<=I_E) KEEP_XDIR=.TRUE.
3168                       ELSE
3169                          DO NS = 1,NXS
3170                             IF(I==I_SLICE(NS)) KEEP_XDIR=.TRUE.
3171                          ENDDO
3172                       ENDIF
3173     
3174     ! Y-direction
3175                       KEEP_YDIR=.FALSE.
3176                       IF(NYS==0) THEN
3177                          IF(J_S<=J.AND.J<=J_N) KEEP_YDIR=.TRUE.
3178                       ELSE
3179                          DO NS = 1,NYS
3180                             IF(J==J_SLICE(NS)) KEEP_YDIR=.TRUE.
3181                          ENDDO
3182                       ENDIF
3183     
3184     ! Z-direction
3185                       KEEP_ZDIR=.FALSE.
3186                       IF(NZS==0) THEN
3187                          IF(K_B<=K.AND.K<=K_T) KEEP_ZDIR=.TRUE.
3188                       ELSE
3189                          DO NS = 1,NZS
3190                             IF(K==K_SLICE(NS)) KEEP_ZDIR=.TRUE.
3191                          ENDDO
3192                       ENDIF
3193     
3194     ! Now combine
3195                       IF(KEEP_XDIR.AND.KEEP_YDIR.AND.KEEP_ZDIR) THEN
3196                          BELONGS_TO_VTK_SUBDOMAIN(IJK) = .TRUE.
3197                          NUMBER_OF_VTK_CELLS = NUMBER_OF_VTK_CELLS + 1
3198                       ENDIF
3199                    ENDIF
3200                 ENDIF
3201              END DO
3202     
3203           ELSE  ! BDIST_IO
3204     
3205              IF(ALLOCATED(BELONGS_TO_VTK_SUBDOMAIN)) DEALLOCATE(BELONGS_TO_VTK_SUBDOMAIN)
3206     
3207              ALLOCATE (BELONGS_TO_VTK_SUBDOMAIN(ijkend3))
3208     
3209     ! Filter the cells based on the VTK region bounds and set the
3210     ! flag BELONGS_TO_VTK_SUBDOMAIN(IJK) to .TRUE. to keep the cell.
3211     
3212              BELONGS_TO_VTK_SUBDOMAIN = .FALSE.
3213              NUMBER_OF_VTK_CELLS      = 0
3214     
3215              DO IJK = 1,IJKEND3
3216                 IF (INTERIOR_CELL_AT(IJK))      THEN
3217                    IF (.NOT.BLOCKED_CELL_AT(IJK)) THEN
3218                       I = I_OF(IJK)
3219                       J = J_OF(IJK)
3220                       K = K_OF(IJK)
3221     
3222     ! X-direction
3223                       KEEP_XDIR=.FALSE.
3224                       IF(NXS==0) THEN
3225                          IF(I_W<=I.AND.I<=I_E) KEEP_XDIR=.TRUE.
3226                       ELSE
3227                          DO NS = 1,NXS
3228                             IF(I==I_SLICE(NS)) KEEP_XDIR=.TRUE.
3229                          ENDDO
3230                       ENDIF
3231     
3232     ! Y-direction
3233                       KEEP_YDIR=.FALSE.
3234                       IF(NYS==0) THEN
3235                          IF(J_S<=J.AND.J<=J_N) KEEP_YDIR=.TRUE.
3236                       ELSE
3237                          DO NS = 1,NYS
3238                             IF(J==J_SLICE(NS)) KEEP_YDIR=.TRUE.
3239                          ENDDO
3240                       ENDIF
3241     
3242     ! Z-direction
3243                       KEEP_ZDIR=.FALSE.
3244                       IF(NZS==0) THEN
3245                          IF(K_B<=K.AND.K<=K_T) KEEP_ZDIR=.TRUE.
3246                       ELSE
3247                          DO NS = 1,NZS
3248                             IF(K==K_SLICE(NS)) KEEP_ZDIR=.TRUE.
3249                          ENDDO
3250                       ENDIF
3251     
3252     ! Now combine
3253                       IF(KEEP_XDIR.AND.KEEP_YDIR.AND.KEEP_ZDIR) THEN
3254                          BELONGS_TO_VTK_SUBDOMAIN(IJK) = .TRUE.
3255                          NUMBER_OF_VTK_CELLS = NUMBER_OF_VTK_CELLS + 1
3256                       ENDIF
3257                    ENDIF
3258                 ENDIF
3259              END DO
3260     
3261           ENDIF
3262     
3263           RETURN
3264     
3265           END SUBROUTINE SETUP_VTK_REGION
3266