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