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