File: RELATIVE:/../../../mfix.git/model/des/vtp_mod.f

1           MODULE vtp
2     
3           use mpi_utility
4           use cdist
5     
6           use desmpi
7           use mpi_comm_des
8           use error_manager
9     
10           IMPLICIT NONE
11     
12           INTEGER, PRIVATE :: GLOBAL_CNT
13           INTEGER, PRIVATE :: LOCAL_CNT
14     
15           INTEGER :: DES_UNIT = 2000
16     
17     ! file unit for ParaView *.pvd data
18           INTEGER, PARAMETER :: PVD_UNIT = 2050
19     
20     ! formatted file name
21           CHARACTER(LEN=64) :: FNAME_VTP
22     
23           INTERFACE VTP_WRITE_DATA
24              MODULE PROCEDURE VTP_WRITE_DP1
25              MODULE PROCEDURE VTP_WRITE_DP2
26              MODULE PROCEDURE VTP_WRITE_I1
27           END INTERFACE
28     
29           CONTAINS
30     
31     !``````````````````````````````````````````````````````````````````````!
32     ! Subroutine: VTP_WRITE_DP1                                            !
33     !                                                                      !
34     ! Purpose: Collect and write 1D double percision arrays to the VTP     !
35     ! file. This routine is designed to collect the data for parallel and  !
36     ! serial runs. This routine also manages the distribted IO case.       !
37     !``````````````````````````````````````````````````````````````````````!
38           SUBROUTINE VTP_WRITE_DP1(NAME, DATA)
39     
40           CHARACTER(len=*), INTENT(in) :: NAME
41           DOUBLE PRECISION, INTENT(in) :: DATA(:)
42     
43           INTEGER :: LC, PC
44     
45           IF(bDist_IO) THEN
46     
47              WRITE(DES_UNIT,1000) NAME
48     
49              PC = 1
50              DO LC = 1, MAX_PIP
51                 IF(PC > PIP) EXIT
52                 IF(IS_NONEXISTENT(LC)) CYCLE
53                 PC = PC+1
54                 IF(IS_GHOST(LC) .OR. IS_ENTERING_GHOST(LC) .OR. IS_EXITING_GHOST(LC)) CYCLE
55                 WRITE(DES_UNIT, 1001,ADVANCE="NO") real(DATA(LC))
56              ENDDO
57              WRITE(DES_UNIT,1002)
58     
59           ELSE
60     
61              allocate (dProcBuf(LOCAL_CNT) )
62              allocate (dRootBuf(GLOBAL_CNT))
63     
64              CALL DES_GATHER(DATA)
65     
66              IF(myPE == PE_IO) THEN
67                 WRITE(DES_UNIT,1000) NAME
68                 DO LC=1, GLOBAL_CNT
69                    WRITE(DES_UNIT,1001,ADVANCE="NO") real(drootbuf(LC))
70                 ENDDO
71                 WRITE(DES_UNIT,1002)
72              ENDIF
73     
74              deallocate(dProcBuf, dRootBuf)
75     
76           ENDIF
77     
78      1000 FORMAT('<DataArray type="Float32" Name="',A,'" format="ascii">')
79      1001 FORMAT(ES14.6,1X)
80      1002 FORMAT('</DataArray>')
81     
82           END SUBROUTINE VTP_WRITE_DP1
83     
84     !``````````````````````````````````````````````````````````````````````!
85     ! Subroutine: VTP_WRITE_DP2                                            !
86     !                                                                      !
87     ! Purpose: Collect and write 2D double percision arrays to the VTP     !
88     ! file. This routine is designed to collect the data for parallel and  !
89     ! serial runs. This routine also manages the distribted IO case.       !
90     !``````````````````````````````````````````````````````````````````````!
91           SUBROUTINE VTP_WRITE_DP2(NAME, DATA)
92     
93           CHARACTER(len=*), INTENT(in) :: NAME
94           DOUBLE PRECISION, INTENT(in) :: DATA(:,:)
95     
96           DOUBLE PRECISION, ALLOCATABLE :: ltemp_array(:,:)
97     
98           CHARACTER(len=16) :: NOC
99           INTEGER :: LB, UB
100           INTEGER :: PC, LC1, LC2
101     
102           LB = LBOUND(DATA,1)
103           UB = UBOUND(DATA,1)
104           NOC=''; WRITE(NOC,*) (UB-LB)+1
105     
106           IF(bDist_IO) THEN
107     
108              WRITE(DES_UNIT,1000) NAME, trim(adjustl(NOC))
109     
110              PC = 1
111              DO LC1 = 1, MAX_PIP
112                 IF(PC > PIP) EXIT
113                 IF(IS_NONEXISTENT(LC1)) CYCLE
114                 PC = PC+1
115                 IF(IS_GHOST(LC1) .OR. IS_ENTERING_GHOST(LC1) .OR. IS_EXITING_GHOST(LC1)) CYCLE
116                 DO LC2=LB, UB
117                    WRITE(DES_UNIT,1001,ADVANCE="NO") real(DATA(LC1,LC2))
118                 ENDDO
119              ENDDO
120              WRITE(DES_UNIT,1002)
121     
122           ELSE
123     
124              allocate (dProcBuf(LOCAL_CNT) )
125              allocate (dRootBuf(GLOBAL_CNT))
126              allocate (ltemp_array((UB-LB)+1,GLOBAL_CNT))
127     
128              DO LC1 = LB, UB
129                 CALL DES_GATHER(DATA(LC1,:))
130                 ltemp_array(LC1,:) = drootbuf(:)
131              ENDDO
132     
133              IF(myPE == PE_IO) THEN
134                 WRITE(DES_UNIT,1000) NAME, trim(adjustl(NOC))
135                 DO LC1=1, GLOBAL_CNT
136                    DO LC2=LB, UB
137                       WRITE(DES_UNIT,1001,ADVANCE="NO") &
138                          real(ltemp_array(LC2,LC1))
139                    ENDDO
140                 ENDDO
141                 WRITE(DES_UNIT,1002)
142              ENDIF
143     
144              deallocate (dProcBuf, dRootBuf, ltemp_array)
145     
146           ENDIF
147     
148     
149      1000 FORMAT('<DataArray type="Float32" Name="',A,'" NumberOf',        &
150              'Components="',A,'" format="ascii">')
151      1001 FORMAT(ES14.6,1X)
152      1002 FORMAT('</DataArray>')
153     
154           END SUBROUTINE VTP_WRITE_DP2
155     
156     
157     
158     !``````````````````````````````````````````````````````````````````````!
159     ! Subroutine: VTP_WRITE_I1                                             !
160     !                                                                      !
161     ! Purpose: Collect and write 1D integer arrays to the VTP file. This   !
162     ! routine is designed to collect the data for parallel and serial      !
163     ! runs. This routine also manages the distribted IO case.              !
164     !``````````````````````````````````````````````````````````````````````!
165           SUBROUTINE VTP_WRITE_I1(NAME, DATA)
166     
167           CHARACTER(len=*), INTENT(in) :: NAME
168           INTEGER, INTENT(in) :: DATA(:)
169     
170           INTEGER :: LC, PC
171     
172           IF(bDist_IO) THEN
173     
174              WRITE(DES_UNIT,1000) NAME
175     
176              PC = 1
177              DO LC = 1, MAX_PIP
178                 IF(PC > PIP) EXIT
179                 IF(IS_NONEXISTENT(LC)) CYCLE
180                 PC = PC+1
181                 IF(IS_GHOST(LC) .OR. IS_ENTERING_GHOST(LC) .OR. IS_EXITING_GHOST(LC)) CYCLE
182                 WRITE(DES_UNIT, 1001,ADVANCE="NO") DATA(LC)
183              ENDDO
184              WRITE(DES_UNIT,1002)
185     
186           ELSE
187     
188              allocate (iProcBuf(LOCAL_CNT) )
189              allocate (iRootBuf(GLOBAL_CNT))
190     
191              CALL DES_GATHER(DATA)
192     
193              IF(myPE == PE_IO) THEN
194                 WRITE(DES_UNIT,1000) NAME
195                 DO LC=1, GLOBAL_CNT
196                    WRITE(DES_UNIT,1001,ADVANCE="NO") irootbuf(LC)
197                 ENDDO
198                 WRITE(DES_UNIT,1002)
199              ENDIF
200     
201              deallocate(iProcBuf, iRootBuf)
202     
203           ENDIF
204     
205      1000 FORMAT('<DataArray type="Float32" Name="',A,'" format="ascii">')
206      1001 FORMAT(I10,1X)
207      1002 FORMAT('</DataArray>')
208     
209           END SUBROUTINE VTP_WRITE_I1
210     
211     
212     !``````````````````````````````````````````````````````````````````````!
213     ! Subroutine: VTP_WRITE_ELEMENT                                        !
214     !                                                                      !
215     ! Purpose: Write a string to the VTP file. It masks the need to check  !
216     ! the logical before flushing.                                         !
217     !``````````````````````````````````````````````````````````````````````!
218           SUBROUTINE VTP_WRITE_ELEMENT(ELEMENT)
219     
220           CHARACTER(len=*), INTENT(in) :: ELEMENT
221     
222           IF(bDist_IO .OR. myPE == PE_IO) &
223              WRITE(DES_UNIT,"(A)") ELEMENT
224     
225           RETURN
226           END SUBROUTINE VTP_WRITE_ELEMENT
227     
228     
229     
230     !``````````````````````````````````````````````````````````````````````!
231     ! Subroutine: VTP_OPEN_FILE                                            !
232     !                                                                      !
233     ! Purpose: This routine opens the VTP file and calcualtes the offsets  !
234     ! for dmp data collection.                                             !
235     !``````````````````````````````````````````````````````````````````````!
236           SUBROUTINE VTP_OPEN_FILE(NoPc)
237     
238     ! Modules
239     !-----------------------------------------------
240           USE run, only: RUN_TYPE, RUN_NAME
241     
242           IMPLICIT NONE
243     
244           CHARACTER(len=*) :: NoPc
245     
246           INTEGER :: NumberOfPoints
247     
248     ! Variables related to gather
249           integer lgathercnts(0:numpes-1), lproc
250     
251     ! check whether an error occurs in opening a file
252           INTEGER :: IOS
253     ! Integer error flag.
254           INTEGER :: IER
255     
256     ! logical used for testing is the data file already exists
257           LOGICAL :: EXISTS_VTP
258     ! status of the vtp file to be written
259           CHARACTER(LEN=8) :: STATUS_VTP
260     
261     
262     ! Initial the global count.
263           GLOBAL_CNT = 10
264     ! Calculate the number of 'real' particles on the local process.
265           LOCAL_CNT = PIP - iGHOST_CNT
266     
267     ! Distributed IO
268           IF(bDIST_IO) THEN
269              NumberOfPoints = LOCAL_CNT
270              WRITE(NoPc,"(I10.10)") NumberOfPoints
271     
272              WRITE(fname_vtp,'(A,"_DES",I4.4,"_",I5.5,".vtp")') &
273                 trim(run_name), vtp_findex, mype
274     
275     ! Serial IO
276           ELSE
277     
278     ! Calculate the total number of particles system-wide.
279              call global_sum(LOCAL_CNT, GLOBAL_CNT)
280              NumberOfPoints = GLOBAL_CNT
281              WRITE(NoPc,"(I10.10)") NumberOfPoints
282     
283     ! Set the send count from the local process.
284              igath_sendcnt = LOCAL_CNT
285     
286     ! Collect the number of particles on each rank.all ranks.
287              lgathercnts = 0
288              lgathercnts(myPE) = LOCAL_CNT
289              call global_sum(lgathercnts,igathercnts)
290     
291     ! Calculate the rank displacements.
292              idispls(0) = 0
293              DO lPROC = 1,NUMPEs-1
294                 idispls(lproc) = idispls(lproc-1) + igathercnts(lproc-1)
295              ENDDO
296     
297     ! set the file name and unit number and open file
298              WRITE(fname_vtp,'(A,"_DES_",I5.5,".vtp")') &
299                 trim(run_name), vtp_findex
300           ENDIF
301     
302           IER = 0
303           IF(bDIST_IO .OR. myPE == PE_IO) THEN
304     
305     ! The file should be new but could exist due to restarting.
306              STATUS_VTP = 'NEW'
307     ! Check to see if the file already exists.
308              INQUIRE(FILE=FNAME_VTP,EXIST=EXISTS_VTP)
309     ! The given file should not exist if the run type is NEW.
310              IF(EXISTS_VTP)THEN
311     ! The VTP should never exist for a NEW run.
312                 IF(RUN_TYPE == 'NEW')THEN
313                    IER = 1
314     ! The file may exist during a RESTART.
315                 ELSE
316                    STATUS_VTP = 'REPLACE'
317                 ENDIF
318              ENDIF
319     
320     ! Open the file and record any erros.
321              IF(IER == 0) THEN
322                 OPEN(CONVERT='BIG_ENDIAN',UNIT=DES_UNIT, FILE=FNAME_VTP,                        &
323                    STATUS=STATUS_VTP, IOSTAT=IOS)
324                 IF(IOS /= 0) IER = 2
325              ENDIF
326           ENDIF
327     
328           CALL GLOBAL_ALL_MAX(IER)
329     
330           IF(IER /= 0) THEN
331              CALL INIT_ERR_MSG("VTP_MOD --> OPEN_VTP")
332              WRITE(ERR_MSG, 1100) IER
333              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
334           ENDIF
335     
336      1100 FORMAT('Error 1100: Unable to open VTP file. This could be ',    &
337              'caused by a VTP',/'file with the same file name already ',   &
338              'existing. or an error code',/' returned by the OPEN ',       &
339              'function.'/'Error code: ',I2,4x,'Aborting.')
340     
341     
342           END SUBROUTINE VTP_OPEN_FILE
343     
344     
345     
346     !......................................................................!
347     ! SUBROUTINE: VTP_CLOSE_FILE                                           !
348     !                                                                      !
349     ! Purpose: This routine closes the vtp file.                           !
350     !``````````````````````````````````````````````````````````````````````!
351           SUBROUTINE VTP_CLOSE_FILE
352     
353     
354           VTP_FINDEX=VTP_FINDEX+1
355     
356           IF(bDist_io .OR. (myPE .eq.pe_IO)) CLOSE(des_unit)
357     
358     
359           END SUBROUTINE VTP_CLOSE_FILE
360     
361     
362     !......................................................................!
363     ! SUBROUTINE: ADD_VTP_TO_PVD                                           !
364     !                                                                      !
365     ! Purpose: This routine opens the pvd file.                            !
366     !``````````````````````````````````````````````````````````````````````!
367           SUBROUTINE ADD_VTP_TO_PVD
368     
369           use run, only: RUN_TYPE, RUN_NAME
370     
371     !-----------------------------------------------
372     ! Local Variables
373     !-----------------------------------------------
374     ! Index position of desired character
375           INTEGER IDX_f, IDX_b
376     ! logical used for testing is the data file already exists
377           LOGICAL :: EXISTS_PVD
378     ! Generic input limited to 256 characters
379           CHARACTER(LEN=256) INPUT
380     
381     ! formatted file name
382           CHARACTER(LEN=64) :: FNAME_PVD = ''
383     
384           LOGICAL, SAVE :: FIRST_PASS = .TRUE.
385     
386     ! IO Status flag
387           INTEGER :: IOS
388     
389     ! Variables related to gather
390           integer :: IER
391     
392     !-----------------------------------------------
393     
394           CALL INIT_ERR_MSG('VTP_MOD --> ADD_VTP_TO_PVD')
395     
396     ! Initialize the error flag.
397           IER = 0
398     
399     ! Obtain the file name and open the pvd file
400           FNAME_PVD = TRIM(RUN_NAME)//'_DES.pvd'
401     
402     ! The PVD file is only written by PE_IO with serial IO.
403           IF(myPE == PE_IO .AND. .NOT.bDist_IO) THEN
404     
405     ! Check to see if the file already exists.
406              INQUIRE(FILE=FNAME_PVD,EXIST=EXISTS_PVD)
407     
408              IF(FIRST_PASS) THEN
409     
410     ! Open the "NEW" file and write the necessary header information.
411                 IF(RUN_TYPE /= 'RESTART_1')THEN
412     
413     ! The file exists but first_pass is also true so most likely an existing
414     ! file from an earlier/other run is present in the directory. Exit to
415     ! prevent accidently overwriting the existing file.
416                    IF(EXISTS_PVD) THEN
417                       IER = 1
418                    ELSE
419                       OPEN(CONVERT='BIG_ENDIAN',UNIT=PVD_UNIT,FILE=FNAME_PVD,STATUS='NEW')
420                       WRITE(PVD_UNIT,"(A)")'<?xml version="1.0"?>'
421                       WRITE(PVD_UNIT,"(A)")'<VTKFile type="Collection" &
422                          &version="0.1" byte_order="LittleEndian">'
423                       WRITE(PVD_UNIT,"(3X,'<Collection>')")
424     ! write two generic lines that will be removed later
425                       WRITE(PVD_UNIT,"('SCRAP LINE 1')")
426                       WRITE(PVD_UNIT,"('SCRAP LINE 2')")
427                    ENDIF
428     
429     ! This is the first pass of a restart run. Extra care is needed to make
430     ! sure that the pvd file is ready to accept new data.
431                 ELSE ! a restart run
432                    IF(EXISTS_PVD) THEN
433     ! Open the file at the beginning.
434                       OPEN(CONVERT='BIG_ENDIAN',UNIT=PVD_UNIT,FILE=FNAME_PVD,&
435                          POSITION="REWIND",STATUS='OLD',IOSTAT=IOS)
436                       IF(IOS /= 0) IER = 2
437                    ELSE ! a pvd file does not exist
438                       IER = 3
439                    ENDIF
440     
441                    IF(IER == 0) THEN
442     ! Loop over the entries in the PVD file, looking for a match to the
443     ! file that is being written. If no match is found, the data will be
444     ! appended to the end of the pvd file, otherwise, the old data will
445     ! be over-written.
446                       DO
447     ! Read in the entires of the PVD file.
448                          READ(PVD_UNIT,"(A)",IOSTAT=IOS)INPUT
449                          IF(IOS > 0) THEN
450                             IER = 4
451                             EXIT
452                          ELSEIF(IOS<0)THEN
453     ! The end of the pvd file has been reached without finding an entry
454     ! matching the current record. Exit the loop.
455                             BACKSPACE(PVD_UNIT)
456                             EXIT
457                          ENDIF
458     ! Find the first instances of file=" and "/> in the read data.
459                          IDX_f = INDEX(INPUT,'file="')
460                          IDX_b = INDEX(INPUT,'"/>')
461     ! Skip rows that do not contain file data
462                          IF(IDX_f == 0 .AND. IDX_b == 0) CYCLE
463     ! Truncate the file name from the read data
464                          WRITE (INPUT,"(A)") INPUT(IDX_f+6:IDX_b-1)
465     ! If the file name matches the current VTP record, break the loop to
466     ! over-write this record.
467                          IF(TRIM(FNAME_VTP) == TRIM(INPUT)) THEN
468                             BACKSPACE(PVD_UNIT)
469                             EXIT
470                          ENDIF
471                       ENDDO
472                    ENDIF ! No errors
473                 ENDIF ! run_type new or restart
474     ! Identify that the files has been created and opened for next pass
475                 FIRST_PASS = .FALSE.
476     
477              ELSE ! not FIRST_PASS
478                 OPEN(CONVERT='BIG_ENDIAN',UNIT=PVD_UNIT,FILE=FNAME_PVD,&
479                    POSITION="APPEND",STATUS='OLD',IOSTAT=IOS)
480                 IF (IOS /= 0) IER = 2
481              ENDIF
482     
483           ENDIF ! if myPE == PE_IO and not distributed IO
484     
485     
486           CAlL GLOBAL_ALL_SUM(IER)
487           IF(IER /= 0) THEN
488              SELECT CASE(IER)
489              CASE(1); WRITE(ERR_MSG,1101) trim(FNAME_PVD)
490              CASE(2); WRITE(ERR_MSG,1102) trim(FNAME_PVD)
491              CASE(3); WRITE(ERR_MSG,1103) trim(FNAME_PVD)
492              CASE(4); WRITE(ERR_MSG,1104) trim(FNAME_PVD)
493              CASE DEFAULT; WRITE(ERR_MSG,1105) trim(FNAME_PVD)
494              END SELECT
495              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
496           ENDIF
497     
498      1101 FORMAT('Error 1101: A PVD file was detected in the run ',        &
499              'directory which should',/'not exist for a NEW run.',/        &
500              'File: ',A)
501     
502      1102 FORMAT('Error 1102: Fatal error status returned while OPENING ', &
503              'PVD file.',/'File: ', A)
504     
505      1103 FORMAT('Error 1103: PVD file MISSING from run directory.',/      &
506              'File: ',A)
507     
508      1104 FORMAT('Error 1104: Fatal error status returned while READING ', &
509              'PVD file.',/'File: ', A)
510     
511      1105 FORMAT('Error 1105:: Fatal unclassified error when processing ', &
512              'PVD file.',/'File: ', A)
513     
514     
515     ! If there were no errors, updated the file.
516           IF(myPE == PE_IO .AND. .NOT.bDist_IO) THEN
517     
518     ! Remove the last two lines written so that additional data can be added
519              BACKSPACE(PVD_UNIT)
520              BACKSPACE(PVD_UNIT)
521     
522     ! Write the data to the file
523              WRITE(PVD_UNIT,"(6X,A,A,A,A,A,A,A)")&
524              '<DataSet timestep="',TRIM(iVal(S_TIME)),'" ',&
525              'group="" part="0" ',& ! necessary file data
526              'file="',TRIM(FNAME_VTP),'"/>' ! file name of vtp
527     
528     ! Write the closing tags
529              WRITE(PVD_UNIT,"(3X,A)")'</Collection>'
530              WRITE(PVD_UNIT,"(A)")'</VTKFile>'
531     
532              CLOSE(PVD_UNIT)
533           ENDIF
534     
535           CALL FINL_ERR_MSG
536     
537     ! Return to the calling routine
538           RETURN
539     
540           END SUBROUTINE ADD_VTP_TO_PVD
541     
542     
543     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
544     !                                                                      C
545     !  Module name: WRITE_VTP_FILE                                         C
546     !  Purpose: Writes particles data in VTK format (Polydata VTP)         C
547     !           Binary format                                              C
548     !                                                                      C
549     !  Author: Jeff Dietiker                              Date: 11-Feb-15  C
550     !  Reviewer:                                          Date:            C
551     !                                                                      C
552     !  Revision Number #                                  Date: ##-###-##  C
553     !  Author: #                                                           C
554     !  Purpose: #                                                          C
555     !                                                                      C
556     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
557           SUBROUTINE WRITE_VTP_FILE(LCV,MODE)
558     
559           USE vtk, only: DIMENSION_VTK, VTK_DEFINED, FRAME
560           USE vtk, only: VTK_REGION,VTK_DEFINED,VTK_DATA
561           USE vtk, only: VTK_PART_DIAMETER,VTK_PART_VEL,VTK_PART_USR_VAR,VTK_PART_TEMP
562           USE vtk, only: VTK_PART_ANGULAR_VEL,VTK_PART_ORIENTATION
563           USE vtk, only: VTK_PART_X_S, VTK_PART_COHESION
564           USE vtk, only: TIME_DEPENDENT_FILENAME,VTU_FRAME_UNIT,VTU_FRAME_FILENAME
565           USE vtk, only: VTK_DBG_FILE
566           USE output, only: FULL_LOG
567     
568           IMPLICIT NONE
569           INTEGER :: L,N,LCV
570     
571           INTEGER :: PASS
572           INTEGER :: WRITE_HEADER = 1
573           INTEGER :: WRITE_DATA   = 2
574           INTEGER :: MODE   ! MODE = 0 : Write regular VTK region file
575                             ! MODE = 1 : Write debug   VTK region file (VTK_DBG_FILE = .TRUE.)
576     
577           VTK_REGION = LCV
578     ! There is nothing to write if we are not in a defined vtk region
579           IF(.NOT.VTK_DEFINED(VTK_REGION)) RETURN
580     
581           IF(VTK_DATA(LCV)/='P') RETURN
582           IF(MODE==0.AND.(VTK_DBG_FILE(LCV))) RETURN
583           IF(MODE==1.AND.(.NOT.VTK_DBG_FILE(LCV))) RETURN
584     
585           CALL SETUP_VTK_REGION_PARTICLES
586     
587           CALL OPEN_VTP_FILE_BIN(MODE)
588     
589     ! Only open pvd file when there are particles in vtk region
590           IF(GLOBAL_CNT>0.AND.MODE==0) CALL OPEN_PVD_FILE
591     
592     ! First pass write the data header.
593     ! Second pass writes the data (appended binary format).
594     
595           DO PASS=WRITE_HEADER,WRITE_DATA
596     
597     
598              CALL WRITE_GEOMETRY_IN_VTP_BIN(PASS)
599     
600              IF(VTK_PART_DIAMETER(VTK_REGION)) &
601                 CALL WRITE_SCALAR_IN_VTP_BIN('Diameter',2.0D0*DES_RADIUS,PASS)
602     
603              IF(VTK_PART_VEL(VTK_REGION)) &
604                 CALL WRITE_VECTOR_IN_VTP_BIN('Velocity',DES_VEL_NEW,PASS)
605     
606              IF(VTK_PART_ANGULAR_VEL(VTK_REGION)) &
607                 CALL WRITE_VECTOR_IN_VTP_BIN('Angular_velocity', OMEGA_NEW,PASS)
608     
609              IF(PARTICLE_ORIENTATION) THEN
610                 IF(VTK_PART_ORIENTATION(VTK_REGION)) &
611                    CALL WRITE_VECTOR_IN_VTP_BIN('Orientation', ORIENTATION,PASS)
612              ENDIF
613     
614              DO N=1, DES_USR_VAR_SIZE
615                 IF(VTK_PART_USR_VAR(VTK_REGION,N)) &
616                   CALL WRITE_SCALAR_IN_VTP_BIN('User Defined Var',DES_USR_VAR(N,:),PASS)
617              ENDDO
618     
619              IF(ENERGY_EQ.AND.VTK_PART_TEMP(VTK_REGION)) &
620                CALL WRITE_SCALAR_IN_VTP_BIN('Temperature', DES_T_s_NEW,PASS)
621     
622              IF(ANY_SPECIES_EQ) THEN
623                 DO N=1, DIMENSION_N_S
624                    IF(VTK_PART_X_s(VTK_REGION,N)) &
625                      CALL WRITE_SCALAR_IN_VTP_BIN(trim(iVar('X_s',N)), DES_X_s(:,N),PASS)
626                 ENDDO
627              ENDIF
628     
629           IF(USE_COHESION.AND.VTK_PART_COHESION(VTK_REGION)) &
630              CALL WRITE_SCALAR_IN_VTP_BIN('CohesiveForce', PostCohesive,PASS)
631     
632           ENDDO ! PASS LOOP, EITHER HEADER OR DATA
633     
634     
635           CALL CLOSE_VTP_FILE_BIN(MODE)
636     
637     ! Only update pvd file when there are particles in vtk region
638           IF(GLOBAL_CNT>0.AND.MODE==0) CALL UPDATE_AND_CLOSE_PVD_FILE
639     
640           call MPI_barrier(MPI_COMM_WORLD,mpierr)
641     
642     ! Update Frames
643           IF (myPE == PE_IO.AND.TIME_DEPENDENT_FILENAME) THEN
644              OPEN(CONVERT='BIG_ENDIAN',UNIT = VTU_FRAME_UNIT, FILE = TRIM(VTU_FRAME_FILENAME))
645              DO L = 1, DIMENSION_VTK
646                 IF(VTK_DEFINED(L)) WRITE(VTU_FRAME_UNIT,*) L,FRAME(L)
647              ENDDO
648              CLOSE(VTU_FRAME_UNIT)
649           ENDIF
650     
651          IF (FULL_LOG.AND.myPE == PE_IO) WRITE(*,20)' DONE.'
652     
653     20    FORMAT(A,1X/)
654           RETURN
655     
656           END SUBROUTINE WRITE_VTP_FILE
657     
658     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
659     !                                                                      C
660     !  Module name: OPEN_VTP_FILE                                          C
661     !  Purpose: Open a vtp file and writes the header                      C
662     !           Binary format                                              C
663     !                                                                      C
664     !  Author: Jeff Dietiker                              Date: 11-Feb-15  C
665     !  Reviewer:                                          Date:            C
666     !                                                                      C
667     !  Revision Number #                                  Date: ##-###-##  C
668     !  Author: #                                                           C
669     !  Purpose: #                                                          C
670     !                                                                      C
671     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
672       SUBROUTINE OPEN_VTP_FILE_BIN(MODE)
673     
674           USE run, only: TIME
675           USE output, only: FULL_LOG
676           USE vtk, only: TIME_DEPENDENT_FILENAME, VTU_FRAME_FILENAME, VTU_FRAME_UNIT
677           USE vtk, only: RESET_FRAME_AT_TIME_ZERO,PVTU_FILENAME,PVTU_UNIT,BUFFER,END_REC
678           USE vtk, only: DIMENSION_VTK, VTK_DEFINED,  FRAME,  VTK_REGION
679           USE vtk, only: VTU_FILENAME, VTK_FILEBASE, VTU_DIR, VTU_UNIT
680     
681           IMPLICIT NONE
682           LOGICAL :: VTU_FRAME_FILE_EXISTS, NEED_TO_WRITE_VTP
683           INTEGER :: ISTAT,BUFF1,BUFF2,L
684           INTEGER :: MODE   ! MODE = 0 : Write regular VTK region file
685                             ! MODE = 1 : Write debug   VTK region file (VTK_DBG_FILE = .TRUE.)
686     
687     
688           IF(BDIST_IO) THEN
689              NEED_TO_WRITE_VTP = (LOCAL_CNT>0)
690           ELSE
691              NEED_TO_WRITE_VTP = (MyPE==0.AND.GLOBAL_CNT>0)
692           ENDIF
693     
694     ! Only open the file from head node when not using distributed I/O
695           IF (myPE /= PE_IO.AND.(.NOT.BDIST_IO)) RETURN
696     
697           IF(TIME_DEPENDENT_FILENAME) THEN
698              INQUIRE(FILE=VTU_FRAME_FILENAME,EXIST=VTU_FRAME_FILE_EXISTS)
699              IF(VTU_FRAME_FILE_EXISTS) THEN
700                 OPEN(CONVERT='BIG_ENDIAN',UNIT = VTU_FRAME_UNIT, FILE = TRIM(VTU_FRAME_FILENAME))
701                 DO L = 1, DIMENSION_VTK
702                    IF(VTK_DEFINED(L)) THEN
703                       READ(VTU_FRAME_UNIT,*)BUFF1,BUFF2
704                       FRAME(L)=BUFF2
705                    ENDIF
706                 ENDDO
707                 CLOSE(VTU_FRAME_UNIT)
708              ENDIF
709              IF(RESET_FRAME_AT_TIME_ZERO.AND.TIME==ZERO) THEN
710                 DO L = 1, DIMENSION_VTK
711                    IF(L==VTK_REGION) FRAME(L)=-1
712                 ENDDO
713              ENDIF
714              DO L = 1, DIMENSION_VTK
715                 IF(L==VTK_REGION) FRAME(L) = FRAME(L) + 1
716              ENDDO
717           ENDIF
718     
719     ! For distributed I/O, define the file name for each processor that owns particles
720           IF (BDIST_IO) THEN
721              IF (LOCAL_CNT>0) THEN
722                 IF(TIME_DEPENDENT_FILENAME.AND.MODE==0) THEN
723                    WRITE(VTU_FILENAME,20) TRIM(VTK_FILEBASE(VTK_REGION)),FRAME(VTK_REGION),MYPE
724                 ELSE
725                    WRITE(VTU_FILENAME,25) TRIM(VTK_FILEBASE(VTK_REGION)),MYPE
726                 ENDIF
727              ENDIF
728           ELSE
729              IF(MYPE.EQ.PE_IO) THEN
730                 IF(TIME_DEPENDENT_FILENAME.AND.MODE==0) THEN
731                    WRITE(VTU_FILENAME,30) TRIM(VTK_FILEBASE(VTK_REGION)),FRAME(VTK_REGION)
732                 ELSE
733                    WRITE(VTU_FILENAME,35) TRIM(VTK_FILEBASE(VTK_REGION))
734                 ENDIF
735              END IF
736           END IF
737     
738     ! Add the VTU directory path if necessary
739     
740           IF (NEED_TO_WRITE_VTP) THEN
741              IF(TRIM(VTU_DIR)/='.') VTU_FILENAME='./'//TRIM(VTU_DIR)//'/'//VTU_FILENAME
742           ENDIF
743     
744     ! Echo
745           IF (FULL_LOG) THEN
746              IF (.NOT.BDIST_IO) THEN
747                 WRITE(*,10,ADVANCE='NO')' WRITING VTP FILE : ', TRIM(VTU_FILENAME),' .'
748              ELSE
749                 IF(myPE==PE_IO) WRITE(*,15,ADVANCE='NO')' EACH PROCESOR IS WRITING ITS OWN VTP FILE.'
750              ENDIF
751           ENDIF
752     
753     ! Open File
754     
755           IF (NEED_TO_WRITE_VTP) THEN
756     
757              VTU_UNIT = 678
758              OPEN(CONVERT='BIG_ENDIAN',UNIT     = VTU_UNIT,           &
759                   FILE     = TRIM(VTU_FILENAME), &
760                   FORM     = 'UNFORMATTED',      &  ! works with gfortran 4.3.4 and ifort 10.1 but may not be supported by all compilers
761                                                     ! use 'BINARY' if 'UNFORMATTED' is not supported
762                   ACCESS   = 'STREAM',           &  ! works with gfortran 4.3.4 and ifort 10.1 but may not be supported by all compilers
763                                                     ! use 'SEQUENTIAL' if 'STREAM' is not supported
764                   ACTION   = 'WRITE',            &
765                   IOSTAT=ISTAT)
766     
767     
768              IF (ISTAT /= 0) THEN
769                 IF(DMP_LOG) WRITE(UNIT_LOG, 1001) VTU_FILENAME, VTU_UNIT,VTU_DIR
770                 IF(FULL_LOG.AND.myPE == PE_IO) WRITE(*, 1001) VTU_FILENAME, VTU_UNIT,VTU_DIR
771                 CALL MFIX_EXIT(myPE)
772              ENDIF
773     
774     
775     1001 FORMAT(/1X,70('*')//, ' From: OPEN_VTP_FILE',/,' Message: ',          &
776                 'Error opening vtp file. Terminating run.',/10X,          &
777                 'File name:  ',A,/10X,                                         &
778                 'VTP_UNIT :  ',i4, /10X,                                       &
779                 'PLEASE VERIFY THAT VTU_DIR EXISTS: ', A, &
780                 /1X,70('*')/)
781     
782     
783     ! Write file Header
784              BUFFER='<?xml version="1.0"?>'
785              WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
786     
787              WRITE(BUFFER,*)'<!-- Time =',TIME,' sec. -->'
788              WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
789     
790              BUFFER='<VTKFile type="PolyData" version="0.1" byte_order="BigEndian">'
791              WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
792     
793              BUFFER='  <PolyData>'
794              WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
795     
796     
797           ENDIF
798     ! For distributed I/O, open .p))vtp file that combines all *.vtp files for a given FRAME
799     ! this is a simple ASCII file
800     
801           IF (myPE == PE_IO.AND.BDIST_IO.AND.GLOBAL_CNT>0) THEN
802     
803              IF(TIME_DEPENDENT_FILENAME.AND.MODE==0) THEN
804                 WRITE(PVTU_FILENAME,40) TRIM(VTK_FILEBASE(VTK_REGION)),FRAME(VTK_REGION)
805              ELSE
806                 WRITE(PVTU_FILENAME,45) TRIM(VTK_FILEBASE(VTK_REGION))
807              ENDIF
808     
809              IF(TRIM(VTU_DIR)/='.') PVTU_FILENAME='./'//TRIM(VTU_DIR)//'/'//PVTU_FILENAME
810     
811              OPEN(CONVERT='BIG_ENDIAN',UNIT = PVTU_UNIT, FILE = TRIM(PVTU_FILENAME))
812     
813              WRITE(PVTU_UNIT,100) '<?xml version="1.0"?>'
814              WRITE(PVTU_UNIT,110) '<!-- Time =',TIME,' sec. -->'
815              WRITE(PVTU_UNIT,120) '<VTKFile type="PPolyData"',&
816                       ' version="0.1" byte_order="BigEndian">'
817     
818              WRITE(PVTU_UNIT,100) '  <PPolyData GhostLevel="0">'
819              WRITE(PVTU_UNIT,100) '      <PPoints>'
820              WRITE(PVTU_UNIT,100) '        <PDataArray type="Float32" Name="coordinates" NumberOfComponents="3" &
821                   &format="appended" offset=" 0" />'
822              WRITE(PVTU_UNIT,100) '      </PPoints>'
823              WRITE(PVTU_UNIT,100) ''
824              WRITE(PVTU_UNIT,100) '      <PPointData Scalars="Diameter" Vectors="Velocity">'
825     
826           ENDIF
827     
828     100   FORMAT(A)
829     110   FORMAT(A,E14.7,A)
830     120   FORMAT(A,A)
831     10    FORMAT(/1X,3A)
832     15    FORMAT(/1X,A)
833     20    FORMAT(A,"_",I4.4,"_",I5.5,".vtp")
834     25    FORMAT(A,"_",I5.5,".vtp")
835     30    FORMAT(A,"_",I4.4,".vtp")
836     35    FORMAT(A,".vtp")
837     40    FORMAT(A,"_",I4.4,".pvtp")
838     45    FORMAT(A,".pvtp")
839     
840           RETURN
841     
842           END SUBROUTINE OPEN_VTP_FILE_BIN
843     
844     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
845     !                                                                      C
846     !  Module name: WRITE_GEOMETRY_IN_VTP_BIN                              C
847     !  Purpose: Write Geometry and connectivity in a vtu file              C
848     !           Binary format                                              C
849     !                                                                      C
850     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
851     !  Reviewer:                                          Date:            C
852     !                                                                      C
853     !  Revision Number #                                  Date: ##-###-##  C
854     !  Author: #                                                           C
855     !  Purpose: #                                                          C
856     !                                                                      C
857     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
858       SUBROUTINE WRITE_GEOMETRY_IN_VTP_BIN(PASS)
859     
860           USE vtk, only: NUMBER_OF_POINTS,BUFFER, VTU_UNIT,END_REC,VTU_OFFSET,BELONGS_TO_VTK_SUBDOMAIN
861     
862           IMPLICIT NONE
863     
864           REAL(c_float) :: float
865           INTEGER(c_int) :: int
866     
867           INTEGER ::     nbytes_vector
868           INTEGER ::     offset_xyz
869     
870           INTEGER :: PASS
871           INTEGER :: WRITE_HEADER = 1
872           INTEGER :: WRITE_DATA   = 2
873     
874           DOUBLE PRECISION, ALLOCATABLE :: ltemp_array(:,:)  ! local
875           DOUBLE PRECISION, ALLOCATABLE :: gtemp_array(:,:)  ! global
876     
877           INTEGER :: LB, UB
878           INTEGER :: PC, LC1, LC2
879     
880     ! Loop through all particles and kee a list of particles belonging to a VTK region
881     
882     ! Since the data is appended (i.e., written after all tags), the
883     ! offset, in number of bytes must be specified.  The offset includes
884     ! the size of the data for each field, plus the size of the integer
885     ! that stores the number of bytes.  this is why the offset of a field
886     ! equals the offset of the previous field plus sizeof(int) plus the
887     ! number of bytes of the field.
888     
889     ! Next, the actual data is written for the geometry (PASS=WRITE_DATA)
890     ! The DATA is converted to single precision to save memory.
891     
892           IF (.NOT.BDIST_IO) THEN
893     ! The number of points in the pvd file is the global number of particles
894     ! computed from SETUP_VTK_REGION_PARTICLES
895     
896              NUMBER_OF_POINTS = GLOBAL_CNT
897     
898     ! Number of bytes of position field (vector,3 components)
899              nbytes_vector       = NUMBER_OF_POINTS * 3 * sizeof(float)
900     
901     ! Offset of each field
902              offset_xyz = 0
903     
904              IF(PASS==WRITE_HEADER) THEN
905                 IF(myPE == PE_IO) THEN
906     
907                    WRITE(BUFFER,*)'    <Piece NumberOfPoints="',NUMBER_OF_POINTS, &
908                          '"  NumberOfVerts="0" NumberOfLines ="0" NumberOfStrips="0" NumberOfPolys="0" >'
909                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
910     
911                    WRITE(BUFFER,*)'      <Points>'
912                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
913     
914                    WRITE(BUFFER,*)'        <DataArray type="Float32" Name="coordinates" NumberOfComponents="3" &
915                                            &format="appended" offset="',offset_xyz,'" />'
916                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
917     
918                    WRITE(BUFFER,*)'      </Points>'
919                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
920     
921                    WRITE(BUFFER,*)'<PointData Scalars="Diameter" Vectors="Velocity"> '!preparing pointData
922                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
923     
924     ! calculate offset for next field
925                    VTU_offset = offset_xyz + sizeof(int) + nbytes_vector
926     
927                 ENDIF
928     
929              ELSEIF(PASS==WRITE_DATA) THEN
930     
931                 IF(myPE == PE_IO) THEN
932     
933                    WRITE(BUFFER,*)'      </PointData>'
934                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
935     
936                    WRITE(BUFFER,*)'      <Verts> </Verts>'
937                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
938     
939                    WRITE(BUFFER,*)'      <Lines> </Lines>'
940                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
941     
942                    WRITE(BUFFER,*)'      <Strips> </Strips>'
943                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
944     
945                    WRITE(BUFFER,*)'      <Polys> </Polys>'
946                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
947     
948                    WRITE(BUFFER,*)'    </Piece>'
949                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
950     
951                    WRITE(BUFFER,*)'  </PolyData>'
952                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
953     
954                    WRITE(BUFFER,*)'  <AppendedData encoding="raw">'
955                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
956     
957     
958     ! Starting raw binary data with an underscore
959     
960                    WRITE(BUFFER,*)'_'
961                    WRITE(VTU_UNIT)TRIM(BUFFER)
962     
963     ! Number of bytes for X,Y,Z coordinates
964                 WRITE(VTU_UNIT) nbytes_vector
965     
966     
967              ENDIF
968     
969              LB = LBOUND(DES_POS_NEW,1) ! This should always be 1
970              UB = UBOUND(DES_POS_NEW,1) ! This should always be 2
971     
972              ALLOCATE (dProcBuf(LOCAL_CNT) )
973              ALLOCATE (dRootBuf(GLOBAL_CNT))
974              ALLOCATE (ltemp_array((UB-LB)+1,LOCAL_CNT))
975              ALLOCATE (gtemp_array((UB-LB)+1,GLOBAL_CNT))
976     
977     ! Pack particle coordinates in a temporary local array
978              PC = 0
979              DO LC1 = 1, MAX_PIP
980                 IF(BELONGS_TO_VTK_SUBDOMAIN(LC1)) THEN
981                    PC =PC + 1
982                    DO LC2=LB, UB
983                       ltemp_array(LC2,PC) = DES_POS_NEW(LC2,LC1)
984                    ENDDO
985                 ENDIF
986                 IF(PC==LOCAL_CNT) EXIT
987              ENDDO
988     
989     ! For each coordinate (x,y, and z), gather the local list to global temporary array
990              DO LC1 = LB, UB
991                 dprocbuf(1:LOCAL_CNT)=ltemp_array(LC1,1:LOCAL_CNT)
992                 CALL desmpi_gatherv(ptype=2)
993                 gtemp_array(LC1,:) = drootbuf(:)
994              ENDDO
995     
996     ! Write the list of coordinates
997              IF(myPE == PE_IO) THEN
998                 DO LC1=1, GLOBAL_CNT
999                    DO LC2=LB, UB
1000                       WRITE(VTU_UNIT)  real(gtemp_array(LC2,LC1))
1001                    ENDDO
1002                 ENDDO
1003              ENDIF
1004     
1005              deallocate (dProcBuf, dRootBuf, ltemp_array,gtemp_array)
1006     
1007     
1008              ENDIF
1009     
1010     
1011           ELSEIF(BDIST_IO.AND.LOCAL_CNT>0) THEN
1012     
1013              IF(LOCAL_CNT==0) RETURN
1014     ! The number of points in the pvd file is the local number of particles
1015     ! computed from SETUP_VTK_REGION_PARTICLES
1016     
1017              NUMBER_OF_POINTS = LOCAL_CNT
1018     
1019     ! Number of bytes of position field (vector,3 components)
1020              nbytes_vector       = NUMBER_OF_POINTS * 3 * sizeof(float)
1021     
1022     ! Offset of each field
1023              offset_xyz = 0
1024     
1025              IF(PASS==WRITE_HEADER) THEN
1026     
1027                 WRITE(BUFFER,*)'    <Piece NumberOfPoints="',NUMBER_OF_POINTS, &
1028                       '"  NumberOfVerts="0" NumberOfLines ="0" NumberOfStrips="0" NumberOfPolys="0" >'
1029                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1030     
1031                 WRITE(BUFFER,*)'      <Points>'
1032                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1033     
1034                 WRITE(BUFFER,*)'        <DataArray type="Float32" Name="coordinates" NumberOfComponents="3" &
1035                                         &format="appended" offset="',offset_xyz,'" />'
1036                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1037     
1038                 WRITE(BUFFER,*)'      </Points>'
1039                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1040     
1041                 WRITE(BUFFER,*)'<PointData Scalars="Diameter" Vectors="Velocity"> '!preparing pointData
1042                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1043     
1044     ! calculate offset for next field
1045                 VTU_offset = offset_xyz + sizeof(int) + nbytes_vector
1046     
1047     
1048              ELSEIF(PASS==WRITE_DATA) THEN
1049     
1050                 WRITE(BUFFER,*)'      </PointData>'
1051                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1052     
1053                 WRITE(BUFFER,*)'      <Verts> </Verts>'
1054                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1055     
1056                 WRITE(BUFFER,*)'      <Lines> </Lines>'
1057                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1058     
1059                 WRITE(BUFFER,*)'      <Strips> </Strips>'
1060                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1061     
1062                 WRITE(BUFFER,*)'      <Polys> </Polys>'
1063                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1064     
1065                 WRITE(BUFFER,*)'    </Piece>'
1066                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1067     
1068                 WRITE(BUFFER,*)'  </PolyData>'
1069                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1070     
1071                 WRITE(BUFFER,*)'  <AppendedData encoding="raw">'
1072                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1073     
1074     ! Starting raw binary data with an underscore
1075     
1076                 WRITE(BUFFER,*)'_'
1077                 WRITE(VTU_UNIT)TRIM(BUFFER)
1078     
1079     ! Number of bytes for X,Y,Z coordinates
1080                 WRITE(VTU_UNIT) nbytes_vector
1081     
1082                 LB = LBOUND(DES_POS_NEW,1) ! This should always be 1
1083                 UB = UBOUND(DES_POS_NEW,1) ! This should always be 2
1084     
1085                 ALLOCATE (ltemp_array((UB-LB)+1,LOCAL_CNT))
1086     
1087     ! Pack particle coordinates in a temporary local array
1088                 PC = 0
1089                 DO LC1 = 1, MAX_PIP
1090                    IF(BELONGS_TO_VTK_SUBDOMAIN(LC1)) THEN
1091                       PC =PC + 1
1092                       DO LC2=LB, UB
1093                          ltemp_array(LC2,PC) = DES_POS_NEW(LC2,LC1)
1094                       ENDDO
1095                    ENDIF
1096                    IF(PC==LOCAL_CNT) EXIT
1097                 ENDDO
1098     
1099     
1100     ! Write the list of coordinates
1101                 DO LC1=1, LOCAL_CNT
1102                    DO LC2=LB, UB
1103                       WRITE(VTU_UNIT)  real(ltemp_array(LC2,LC1))
1104                    ENDDO
1105                 ENDDO
1106     
1107                 deallocate (ltemp_array)
1108     
1109              ENDIF
1110     
1111           ENDIF
1112     
1113           RETURN
1114     
1115           END SUBROUTINE WRITE_GEOMETRY_IN_VTP_BIN
1116     
1117     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1118     !                                                                      C
1119     !  Module name: WRITE_SCALAR_IN_VTP_BIN                                C
1120     !  Purpose: Write Scalar variable in a vtp file                        C
1121     !           Binary format                                              C
1122     !                                                                      C
1123     !  Author: Jeff Dietiker                              Date: 11-Feb-15  C
1124     !  Reviewer:                                          Date:            C
1125     !                                                                      C
1126     !  Revision Number #                                  Date: ##-###-##  C
1127     !  Author: #                                                           C
1128     !  Purpose: #                                                          C
1129     !                                                                      C
1130     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1131       SUBROUTINE WRITE_SCALAR_IN_VTP_BIN(VAR_NAME,VAR,PASS)
1132     
1133           USE vtk, only: BUFFER,VTU_OFFSET,VTU_UNIT,PVTU_UNIT
1134           USE vtk, only: END_REC,BELONGS_TO_VTK_SUBDOMAIN
1135           USE output, only: FULL_LOG
1136     
1137           IMPLICIT NONE
1138           INTEGER :: I,LC1,PC
1139     
1140           CHARACTER (*) :: VAR_NAME
1141           DOUBLE PRECISION, INTENT(in) :: VAR(:)
1142     
1143           REAL(c_float) :: float
1144     
1145           INTEGER :: nbytes_scalar
1146     
1147           INTEGER :: PASS
1148           INTEGER :: WRITE_HEADER = 1
1149           INTEGER :: WRITE_DATA   = 2
1150     
1151           IF (.NOT.BDIST_IO) THEN
1152     
1153     ! Number of bytes for each scalar field
1154              nbytes_scalar = GLOBAL_CNT * sizeof(float)
1155     
1156              IF(PASS==WRITE_HEADER) THEN
1157     
1158     ! Remove possible white space with underscore
1159                 DO I = 1,LEN_TRIM(VAR_NAME)
1160                    IF(VAR_NAME(I:I) == ' ') VAR_NAME(I:I) = '_'
1161                 ENDDO
1162     
1163     ! For each scalar, write a tag, with corresponding offset
1164                 WRITE(BUFFER,90)'        <DataArray type="Float32" Name="', &
1165                      TRIM(VAR_NAME),'" format="appended" offset="',VTU_offset,'" />'
1166                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1167     
1168     ! Prepare the offset for the next field
1169                 VTU_offset = VTU_offset + sizeof(float) + nbytes_scalar
1170     
1171     
1172              ELSEIF(PASS==WRITE_DATA) THEN
1173     
1174                allocate (dProcBuf(LOCAL_CNT) )
1175                allocate (dRootBuf(GLOBAL_CNT))
1176     
1177     ! Pack scalar list in a local buffer before gathering to root
1178                 PC = 0
1179                 DO LC1 = 1, MAX_PIP
1180                    IF(BELONGS_TO_VTK_SUBDOMAIN(LC1)) THEN
1181                       PC =PC + 1
1182                       dProcBuf(PC) = VAR(LC1)
1183                    ENDIF
1184                    IF(PC==LOCAL_CNT) EXIT
1185                 ENDDO
1186     
1187     ! Gather local buffer to root
1188              CALL desmpi_gatherv(ptype=2)
1189     
1190     ! Write the data, always preceded by its size in number of bytes
1191     ! Write root buffer to file
1192              WRITE(VTU_UNIT) nbytes_scalar
1193     
1194              IF(myPE == PE_IO) THEN
1195                 DO LC1=1, GLOBAL_CNT
1196                    WRITE(VTU_UNIT)  real(drootBuf(LC1))
1197                 ENDDO
1198              ENDIF
1199     
1200              deallocate (dProcBuf, dRootBuf)
1201     
1202     
1203              ENDIF
1204     
1205     
1206           ELSEIF(BDIST_IO.AND.LOCAL_CNT>0) THEN
1207     
1208     ! Number of bytes for each scalar field
1209              nbytes_scalar = LOCAL_CNT * sizeof(float)
1210     
1211     ! Remove possible white space with underscore
1212              DO I = 1,LEN_TRIM(VAR_NAME)
1213                 IF(VAR_NAME(I:I) == ' ') VAR_NAME(I:I) = '_'
1214              ENDDO
1215     
1216              IF(PASS==WRITE_HEADER) THEN
1217     
1218     ! For each scalar, write a tag, with corresponding offset
1219                 WRITE(BUFFER,90)'        <DataArray type="Float32" Name="', &
1220                      TRIM(VAR_NAME),'" format="appended" offset="',VTU_offset,'" />'
1221                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1222     
1223     ! Prepare the offset for the next field
1224                 VTU_offset = VTU_offset + sizeof(float) + nbytes_scalar
1225     
1226     
1227              ELSEIF(PASS==WRITE_DATA) THEN
1228     
1229                 allocate (dProcBuf(LOCAL_CNT) )
1230     
1231     ! Pack scalar list in a local buffer before writing in file
1232                 PC = 0
1233                 DO LC1 = 1, MAX_PIP
1234                    IF(BELONGS_TO_VTK_SUBDOMAIN(LC1)) THEN
1235                       PC =PC + 1
1236                       dProcBuf(PC) = VAR(LC1)
1237                    ENDIF
1238                    IF(PC==LOCAL_CNT) EXIT
1239                 ENDDO
1240     
1241     ! Write the data, always preceded by its size in number of bytes
1242     ! Write root buffer to file
1243                 WRITE(VTU_UNIT) nbytes_scalar
1244     
1245                 DO LC1=1, LOCAL_CNT
1246                    WRITE(VTU_UNIT)  real(dProcBuf(LC1))
1247                 ENDDO
1248     
1249                 deallocate (dProcBuf)
1250     
1251                 IF (myPE == PE_IO) THEN       ! Update pvtu file with variable name
1252                    WRITE(PVTU_UNIT,90) '        <PointArray type="Float32" Name="', &
1253                         TRIM(VAR_NAME),'" format="appended" offset="',VTU_offset,'" />'
1254                 ENDIF
1255     
1256     
1257              ENDIF
1258     
1259     
1260           ENDIF
1261     
1262     
1263           IF (PASS==WRITE_DATA.AND.FULL_LOG.AND.myPE == PE_IO) WRITE(*,10,ADVANCE='NO')'.'
1264     
1265     10    FORMAT(A)
1266     90    FORMAT(A,A,A,I12,A)
1267     
1268           RETURN
1269     
1270           END SUBROUTINE WRITE_SCALAR_IN_VTP_BIN
1271     
1272     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1273     !                                                                      C
1274     !  Module name: WRITE_VECTOR_IN_VTP                                    C
1275     !  Purpose: Write Vector variable in a vtp file                        C
1276     !           Binary format                                              C
1277     !                                                                      C
1278     !  Author: Jeff Dietiker                              Date: 11-Feb-15  C
1279     !  Reviewer:                                          Date:            C
1280     !                                                                      C
1281     !  Revision Number #                                  Date: ##-###-##  C
1282     !  Author: #                                                           C
1283     !  Purpose: #                                                          C
1284     !                                                                      C
1285     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1286       SUBROUTINE WRITE_VECTOR_IN_VTP_BIN(VAR_NAME,VAR,PASS)
1287     
1288           USE vtk, only: BUFFER,VTU_OFFSET,VTU_UNIT,PVTU_UNIT
1289           USE vtk, only: END_REC,BELONGS_TO_VTK_SUBDOMAIN
1290           USE output, only: FULL_LOG
1291     
1292           IMPLICIT NONE
1293     
1294           CHARACTER (*) :: VAR_NAME
1295           DOUBLE PRECISION, INTENT(in) :: VAR(:,:)
1296     
1297           REAL(c_float) :: float
1298     
1299           INTEGER :: nbytes_vector
1300     
1301           INTEGER :: PASS
1302           INTEGER :: WRITE_HEADER = 1
1303           INTEGER :: WRITE_DATA   = 2
1304     
1305           DOUBLE PRECISION, ALLOCATABLE :: ltemp_array(:,:)  ! local
1306           DOUBLE PRECISION, ALLOCATABLE :: gtemp_array(:,:)  ! global
1307     
1308           INTEGER :: LB, UB
1309           INTEGER :: PC, LC1, LC2
1310     
1311           IF (.NOT.BDIST_IO) THEN
1312     
1313     ! Number of bytes for each vector field
1314              nbytes_vector = GLOBAL_CNT * 3 * sizeof(float)
1315     
1316              IF(PASS==WRITE_HEADER) THEN
1317     ! For each vector, write a tag, with corresponding offset
1318     
1319                 WRITE(BUFFER,90)'        <DataArray type="Float32" Name="', &
1320                      TRIM(VAR_NAME),'"  NumberOfComponents="3" format="appended" offset="',VTU_offset,'" />'
1321                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1322     
1323     ! Prepare the offset for the next field
1324                 VTU_offset = VTU_offset + sizeof(float) + nbytes_vector
1325     
1326     
1327              ELSEIF(PASS==WRITE_DATA) THEN
1328     
1329                 LB = LBOUND(VAR,1) ! This should always be 1
1330                 UB = UBOUND(VAR,1) ! This should always be 2
1331     
1332                 ALLOCATE (dProcBuf(LOCAL_CNT) )
1333                 ALLOCATE (dRootBuf(GLOBAL_CNT))
1334                 ALLOCATE (ltemp_array((UB-LB)+1,LOCAL_CNT))
1335                 ALLOCATE (gtemp_array((UB-LB)+1,GLOBAL_CNT))
1336     
1337     ! For each vector component, pack component list in a local array
1338                 PC = 0
1339                 DO LC1 = 1, MAX_PIP
1340                    IF(BELONGS_TO_VTK_SUBDOMAIN(LC1)) THEN
1341                       PC =PC + 1
1342                       DO LC2=LB, UB
1343                          ltemp_array(LC2,PC) = VAR(LC2,LC1)
1344                       ENDDO
1345                    ENDIF
1346                    IF(PC==LOCAL_CNT) EXIT
1347                 ENDDO
1348     
1349     
1350     ! For each component, gather the local list to global temporary array
1351              DO LC1 = LB, UB
1352                 dprocbuf(1:LOCAL_CNT)=ltemp_array(LC1,1:LOCAL_CNT)
1353                 CALL desmpi_gatherv(ptype=2)
1354                 gtemp_array(LC1,:) = drootbuf(:)
1355              ENDDO
1356     
1357     ! Write the data, always preceded by its size in number of bytes
1358              IF(myPE == PE_IO) THEN
1359                 WRITE(VTU_UNIT) nbytes_vector
1360                 DO LC1=1, GLOBAL_CNT
1361                    DO LC2=LB, UB
1362                       WRITE(VTU_UNIT)  real(gtemp_array(LC2,LC1))
1363                    ENDDO
1364                 ENDDO
1365              ENDIF
1366     
1367              deallocate (dProcBuf, dRootBuf, ltemp_array,gtemp_array)
1368     
1369     
1370     
1371              ENDIF
1372     
1373     
1374           ELSEIF(BDIST_IO.AND.LOCAL_CNT>0) THEN
1375     
1376     ! Number of bytes for each vector field
1377              nbytes_vector = LOCAL_CNT * 3 * sizeof(float)
1378     
1379              IF(PASS==WRITE_HEADER) THEN
1380     ! For each vector, write a tag, with corresponding offset
1381     
1382                 WRITE(BUFFER,90)'        <DataArray type="Float32" Name="', &
1383                      TRIM(VAR_NAME),'"  NumberOfComponents="3" format="appended" offset="',VTU_offset,'" />'
1384                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1385     
1386     ! Prepare the offset for the next field
1387                 VTU_offset = VTU_offset + sizeof(float) + nbytes_vector
1388     
1389     
1390              ELSEIF(PASS==WRITE_DATA) THEN
1391     
1392                 LB = LBOUND(VAR,1) ! This should always be 1
1393                 UB = UBOUND(VAR,1) ! This should always be 2
1394     
1395                 ALLOCATE (ltemp_array((UB-LB)+1,LOCAL_CNT))
1396     
1397     ! For each vector component, pack component list in a local array
1398                 PC = 0
1399                 DO LC1 = 1, MAX_PIP
1400                    IF(BELONGS_TO_VTK_SUBDOMAIN(LC1)) THEN
1401                       PC =PC + 1
1402                       DO LC2=LB, UB
1403                          ltemp_array(LC2,PC) = VAR(LC2,LC1)
1404                       ENDDO
1405                    ENDIF
1406                    IF(PC==LOCAL_CNT) EXIT
1407                 ENDDO
1408     
1409     
1410     ! Write the data, always preceded by its size in number of bytes
1411                 WRITE(VTU_UNIT) nbytes_vector
1412                 DO LC1=1, LOCAL_CNT
1413                    DO LC2=LB, UB
1414                       WRITE(VTU_UNIT)  real(ltemp_array(LC2,LC1))
1415                    ENDDO
1416                 ENDDO
1417     
1418                 deallocate (ltemp_array)
1419     
1420     
1421                 IF (myPE == PE_IO) THEN       ! Update pvtu file with variable name
1422                    WRITE(PVTU_UNIT,90)'        <PointArray type="Float32" Name="', &
1423                      TRIM(VAR_NAME),'"  NumberOfComponents="3" format="appended" offset="',VTU_offset,'" />'
1424                 ENDIF
1425     
1426     
1427              ENDIF
1428     
1429           ENDIF
1430     
1431     
1432           IF (PASS==WRITE_DATA.AND.FULL_LOG.AND.myPE == PE_IO) WRITE(*,10,ADVANCE='NO')'.'
1433     
1434     10    FORMAT(A)
1435     90    FORMAT(A,A,A,I12,A)
1436     
1437           RETURN
1438     
1439           END SUBROUTINE WRITE_VECTOR_IN_VTP_BIN
1440     
1441     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1442     !                                                                      C
1443     !  Module name: CLOSE_VTP_FILE_BIN                                     C
1444     !  Purpose: Close a vtp file                                           C
1445     !           Binary format                                              C
1446     !                                                                      C
1447     !  Author: Jeff Dietiker                              Date: 11-Feb-15  C
1448     !  Reviewer:                                          Date:            C
1449     !                                                                      C
1450     !  Revision Number #                                  Date: ##-###-##  C
1451     !  Author: #                                                           C
1452     !  Purpose: #                                                          C
1453     !                                                                      C
1454     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1455       SUBROUTINE CLOSE_VTP_FILE_BIN(MODE)
1456     
1457           USE vtk, only: BUFFER,VTU_UNIT,END_REC,PVTU_UNIT,TIME_DEPENDENT_FILENAME
1458           USE vtk, only: VTK_REGION,VTK_FILEBASE,FRAME
1459     
1460           IMPLICIT NONE
1461     
1462           INTEGER:: N
1463           CHARACTER (LEN=32)  :: VTU_NAME
1464           INTEGER, DIMENSION(0:numPEs-1) :: ALL_PART_CNT
1465           INTEGER :: IERR
1466           INTEGER :: MODE   ! MODE = 0 : Write regular VTK region file
1467                             ! MODE = 1 : Write debug   VTK region file (VTK_DBG_FILE = .TRUE.)
1468     
1469     
1470           IF((myPE == PE_IO.AND.(.NOT.BDIST_IO)).OR.(BDIST_IO.AND.LOCAL_CNT>0)) THEN
1471     
1472     ! Write last tags and close the vtp file
1473           WRITE(BUFFER,110)'  </AppendedData>'
1474           WRITE(VTU_UNIT)END_REC//TRIM(BUFFER)//END_REC
1475     
1476           WRITE(BUFFER,110)'</VTKFile>'
1477           WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1478     
1479           CLOSE(VTU_UNIT)
1480     
1481           ENDIF
1482     
1483     ! Update pvtu file and close
1484     
1485           IF(BDIST_IO)  THEN
1486              CALL allgather_1i (LOCAL_CNT,ALL_PART_CNT,IERR)
1487     
1488              IF (myPE == PE_IO.AND.GLOBAL_CNT>0) THEN
1489                 WRITE(PVTU_UNIT,100) '      </PPointData>'
1490     
1491                 DO N = 0,NumPEs-1
1492                    IF(ALL_PART_CNT(N)>0) THEN
1493                       IF(TIME_DEPENDENT_FILENAME.AND.MODE==0) THEN
1494                          WRITE(VTU_NAME,20) TRIM(VTK_FILEBASE(VTK_REGION)),FRAME(VTK_REGION),N
1495                       ELSE
1496                          WRITE(VTU_NAME,25) TRIM(VTK_FILEBASE(VTK_REGION)),N
1497                       ENDIF
1498     
1499                       WRITE(PVTU_UNIT,110) '      <Piece Source="',TRIM(VTU_NAME),'"/>'
1500                    ENDIF
1501                 ENDDO
1502     
1503     
1504                 WRITE(PVTU_UNIT,100) '  </PPolyData>'
1505                 WRITE(PVTU_UNIT,100) '</VTKFile>'
1506                 CLOSE(PVTU_UNIT)
1507              ENDIF
1508           ENDIF
1509     
1510     20    FORMAT(A,"_",I4.4,"_",I5.5,".vtp")
1511     25    FORMAT(A,"_",I5.5,".vtp")
1512     
1513     100   FORMAT(A)
1514     110   FORMAT(A,A,A)
1515     
1516           RETURN
1517     
1518           END SUBROUTINE CLOSE_VTP_FILE_BIN
1519     
1520     
1521     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1522     !                                                                      C
1523     !  Module name: SETUP_VTK_REGION_PARTICLES                             C
1524     !                                                                      C
1525     !  Purpose: Filter the particles  based on the VTK region bounds and   C
1526     !           set the flag BELONGS_TO_VTK_SUBDOMAIN to .TRUE.            C
1527     !           to keep the particle.                                      C
1528     !                                                                      C
1529     !  Author: Jeff Dietiker                              Date: 11-Feb-15  C
1530     !  Reviewer:                                          Date:            C
1531     !                                                                      C
1532     !  Revision Number #                                  Date: ##-###-##  C
1533     !  Author: #                                                           C
1534     !  Purpose: #                                                          C
1535     !                                                                      C
1536     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1537           SUBROUTINE SETUP_VTK_REGION_PARTICLES
1538     
1539           USE vtk, only: VTK_REGION
1540           USE vtk, only: VTK_X_E, VTK_X_W, VTK_Y_S, VTK_Y_N, VTK_Z_B, VTK_Z_T
1541           USE vtk, only: VTK_NXS, VTK_NYS, VTK_NZS
1542           USE vtk, only: VTK_SLICE_TOL, VTK_SELECT_MODE
1543           USE vtk, only: BELONGS_TO_VTK_SUBDOMAIN
1544           USE discretelement, only: MAX_PIP,PIP,DES_POS_NEW
1545     
1546           IMPLICIT NONE
1547     
1548           INTEGER :: PC,LC1
1549           INTEGER :: NXS,NYS,NZS,NS
1550           INTEGER :: X_SLICE(DIM_I),Y_SLICE(DIM_J),Z_SLICE(DIM_K)
1551           DOUBLE PRECISION :: XE,XW,YS,YN,ZB,ZT
1552           DOUBLE PRECISION :: XP,YP,ZP,XP1,YP1,ZP1,XP2,YP2,ZP2,R
1553     
1554           DOUBLE PRECISION :: SLICE_TOL
1555           LOGICAL :: KEEP_XDIR,KEEP_YDIR,KEEP_ZDIR
1556     
1557     ! Variables related to gather
1558           integer lgathercnts(0:numpes-1), lproc
1559     
1560           CHARACTER(LEN=1) :: SELECT_PARTICLE_BY
1561     
1562     ! Get VTK region bounds
1563           XE = VTK_X_E(VTK_REGION)
1564           XW = VTK_X_W(VTK_REGION)
1565           YS = VTK_Y_S(VTK_REGION)
1566           YN = VTK_Y_N(VTK_REGION)
1567           ZB = VTK_Z_B(VTK_REGION)
1568           ZT = VTK_Z_T(VTK_REGION)
1569     
1570           NXS = VTK_NXS(VTK_REGION)
1571           NYS = VTK_NYS(VTK_REGION)
1572           NZS = VTK_NZS(VTK_REGION)
1573     
1574           SLICE_TOL = VTK_SLICE_TOL(VTK_REGION)
1575     
1576           SELECT_PARTICLE_BY = VTK_SELECT_MODE(VTK_REGION)
1577     
1578     ! get slice(s) location
1579           DO NS = 1,NXS
1580              X_SLICE(NS) = XW + (XE-XW)/(NXS-1)*(NS-1)
1581           ENDDO
1582     
1583           DO NS = 1,NYS
1584              Y_SLICE(NS) = YS + (YN-YS)/(NYS-1)*(NS-1)
1585           ENDDO
1586     
1587           DO NS = 1,NZS
1588              Z_SLICE(NS) = ZB + (ZT-ZB)/(NZS-1)*(NS-1)
1589           ENDDO
1590     
1591     
1592     
1593     ! Loop through all particles on local rank and keep a list of particles
1594     ! belonging to VTK region
1595     
1596           IF(ALLOCATED(BELONGS_TO_VTK_SUBDOMAIN)) DEALLOCATE(BELONGS_TO_VTK_SUBDOMAIN)
1597           ALLOCATE(BELONGS_TO_VTK_SUBDOMAIN(MAX_PIP))
1598     
1599           BELONGS_TO_VTK_SUBDOMAIN = .FALSE.
1600     
1601           LOCAL_CNT = 0
1602           PC = 1
1603           DO LC1 = 1, MAX_PIP
1604              IF(PC > PIP) EXIT
1605              IF(IS_NONEXISTENT(LC1)) CYCLE
1606              PC = PC+1
1607              IF(IS_GHOST(LC1) .OR. IS_ENTERING_GHOST(LC1) .OR. IS_EXITING_GHOST(LC1)) CYCLE
1608     
1609              SELECT CASE(SELECT_PARTICLE_BY)
1610                 CASE('C')  ! Particle center must be inside vtk region
1611     
1612                    XP = DES_POS_NEW(1,LC1)
1613                    YP = DES_POS_NEW(2,LC1)
1614                    ZP = DES_POS_NEW(3,LC1)
1615     
1616     ! X-direction
1617                    KEEP_XDIR=.FALSE.
1618                    IF(NXS==0) THEN
1619                       IF(XW<=XP.AND.XP<=XE) KEEP_XDIR=.TRUE.
1620                    ELSE
1621                       DO NS = 1,NXS
1622                          IF((X_SLICE(NS)-SLICE_TOL)<=XP.AND.XP<=(X_SLICE(NS)+SLICE_TOL)) KEEP_XDIR=.TRUE.
1623                       ENDDO
1624                    ENDIF
1625     
1626     ! Y-direction
1627                    KEEP_YDIR=.FALSE.
1628                    IF(NYS==0) THEN
1629                       IF(YS<=YP.AND.YP<=YN) KEEP_YDIR=.TRUE.
1630                    ELSE
1631                       DO NS = 1,NYS
1632                          IF((Y_SLICE(NS)-SLICE_TOL)<=YP.AND.YP<=(Y_SLICE(NS)+SLICE_TOL)) KEEP_YDIR=.TRUE.
1633                       ENDDO
1634                    ENDIF
1635     
1636     ! Z-direction
1637                    KEEP_ZDIR=.FALSE.
1638                    IF(NZS==0) THEN
1639                       IF(ZB<=ZP.AND.ZP<=ZT) KEEP_ZDIR=.TRUE.
1640                    ELSE
1641                       DO NS = 1,NZS
1642                          IF((Z_SLICE(NS)-SLICE_TOL)<=ZP.AND.ZP<=(Z_SLICE(NS)+SLICE_TOL)) KEEP_ZDIR=.TRUE.
1643                       ENDDO
1644                    ENDIF
1645     
1646     
1647                 CASE('P')  ! Entire particle must be inside vtk region
1648     
1649                    R = DES_RADIUS(LC1)
1650     
1651                    XP1 = DES_POS_NEW(1,LC1) - R
1652                    YP1 = DES_POS_NEW(2,LC1) - R
1653                    ZP1 = DES_POS_NEW(3,LC1) - R
1654     
1655                    XP2 = DES_POS_NEW(1,LC1) + R
1656                    YP2 = DES_POS_NEW(2,LC1) + R
1657                    ZP2 = DES_POS_NEW(3,LC1) + R
1658     
1659     ! X-direction
1660                    KEEP_XDIR=.FALSE.
1661                    IF(NXS==0) THEN
1662                       IF(XW<=XP1.AND.XP2<=XE) KEEP_XDIR=.TRUE.
1663                    ELSE
1664                       DO NS = 1,NXS
1665                          IF((X_SLICE(NS)-SLICE_TOL)<=XP1.AND.XP2<=(X_SLICE(NS)+SLICE_TOL)) KEEP_XDIR=.TRUE.
1666                       ENDDO
1667                    ENDIF
1668     
1669     ! Y-direction
1670                    KEEP_YDIR=.FALSE.
1671                    IF(NYS==0) THEN
1672                       IF(YS<=YP1.AND.YP2<=YN) KEEP_YDIR=.TRUE.
1673                    ELSE
1674                       DO NS = 1,NYS
1675                          IF((Y_SLICE(NS)-SLICE_TOL)<=YP1.AND.YP2<=(Y_SLICE(NS)+SLICE_TOL)) KEEP_YDIR=.TRUE.
1676                       ENDDO
1677                    ENDIF
1678     
1679     ! Z-direction
1680                    KEEP_ZDIR=.FALSE.
1681                    IF(NZS==0) THEN
1682                       IF(ZB<=ZP1.AND.ZP2<=ZT) KEEP_ZDIR=.TRUE.
1683                    ELSE
1684                       DO NS = 1,NZS
1685                          IF((Z_SLICE(NS)-SLICE_TOL)<=ZP1.AND.ZP2<=(Z_SLICE(NS)+SLICE_TOL)) KEEP_ZDIR=.TRUE.
1686                       ENDDO
1687                    ENDIF
1688     
1689     
1690                 CASE('I')  ! Particle must be inside or intersect the edge of the vtk region
1691     
1692                    R = DES_RADIUS(LC1)
1693     
1694                    XP1 = DES_POS_NEW(1,LC1) - R
1695                    YP1 = DES_POS_NEW(2,LC1) - R
1696                    ZP1 = DES_POS_NEW(3,LC1) - R
1697     
1698                    XP2 = DES_POS_NEW(1,LC1) + R
1699                    YP2 = DES_POS_NEW(2,LC1) + R
1700                    ZP2 = DES_POS_NEW(3,LC1) + R
1701     
1702     ! X-direction
1703                    KEEP_XDIR=.FALSE.
1704                    IF(NXS==0) THEN
1705                       IF(.NOT.(XE<=XP1.OR.XP2<=XW)) KEEP_XDIR=.TRUE.
1706                    ELSE
1707                       DO NS = 1,NXS
1708                          IF(.NOT.((X_SLICE(NS)+SLICE_TOL)<=XP1.OR.XP2<=(X_SLICE(NS)-SLICE_TOL))) KEEP_XDIR=.TRUE.
1709                       ENDDO
1710                    ENDIF
1711     
1712     ! Y-direction
1713                    KEEP_YDIR=.FALSE.
1714                    IF(NYS==0) THEN
1715                       IF(.NOT.(YN<=YP1.OR.YP2<=YS)) KEEP_YDIR=.TRUE.
1716                    ELSE
1717                       DO NS = 1,NYS
1718                          IF(.NOT.((Y_SLICE(NS)+SLICE_TOL)<=YP1.OR.YP2<=(Y_SLICE(NS)-SLICE_TOL))) KEEP_YDIR=.TRUE.
1719                       ENDDO
1720                    ENDIF
1721     
1722     ! Z-direction
1723                    KEEP_ZDIR=.FALSE.
1724                    IF(NZS==0) THEN
1725                       IF(.NOT.(ZT<=ZP1.OR.ZP2<=ZB)) KEEP_ZDIR=.TRUE.
1726                    ELSE
1727                       DO NS = 1,NZS
1728                          IF(.NOT.((Z_SLICE(NS)+SLICE_TOL)<=ZP1.OR.ZP2<=(Z_SLICE(NS)-SLICE_TOL))) KEEP_ZDIR=.TRUE.
1729                       ENDDO
1730                    ENDIF
1731     
1732     
1733                 CASE DEFAULT
1734                       print*,'should not be here'
1735              END SELECT
1736     
1737     ! Now combine
1738              IF(KEEP_XDIR.AND.KEEP_YDIR.AND.KEEP_ZDIR) THEN
1739                 BELONGS_TO_VTK_SUBDOMAIN(LC1) = .TRUE.
1740                 LOCAL_CNT = LOCAL_CNT + 1
1741              ENDIF
1742           ENDDO ! particle loop
1743     
1744     
1745     ! Calculate the total number of particles system-wide.
1746           call global_sum(LOCAL_CNT, GLOBAL_CNT)
1747     
1748     ! No need to set the send/reccv when using distributed IO
1749           IF (BDIST_IO) RETURN
1750     ! Set the send count from the local process.
1751           igath_sendcnt = LOCAL_CNT
1752     
1753     ! Collect the number of particles on each rank.all ranks.
1754           lgathercnts = 0
1755           lgathercnts(myPE) = LOCAL_CNT
1756           call global_sum(lgathercnts,igathercnts)
1757     
1758     ! Calculate the rank displacements.
1759           idispls(0) = 0
1760           DO lPROC = 1,NUMPEs-1
1761              idispls(lproc) = idispls(lproc-1) + igathercnts(lproc-1)
1762           ENDDO
1763     
1764           RETURN
1765     
1766           END SUBROUTINE SETUP_VTK_REGION_PARTICLES
1767     
1768           END MODULE VTP
1769