forked from reedmaxwell/EcoSLIM
-
Notifications
You must be signed in to change notification settings - Fork 0
/
vtk_write_points.f90
101 lines (95 loc) · 3.25 KB
/
vtk_write_points.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
! Copyright 2010 Reed M. Maxwell
!
! This file is part of EcoSLIM; Originally in SLIM-Fast.
!
! SLIM-Fast is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License.
!
! SLIM-Fast is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with SLIM-Fast in /src/gpl.txt. If not, see <http://www.gnu.org/licenses/>.
SUBROUTINE vtk_write_points(P,np_active, np,icycle,vtk_file, dx, dy,nx,ny,maxZ, dem)
REAL*8 :: P(:,:)
INTEGER :: icycle
INTEGER*4 :: np_active
INTEGER*4 :: np
REAL*8 :: dx
REAL*8 :: dy
INTEGER*4 :: nx
INTEGER*4 :: ny
REAL*8 :: maxZ
REAL*8 :: DEM(:,:)
CHARACTER (LEN=200) :: vtk_file
INTEGER*4 i,j,k, ijk, debug,l, nxyz,nxyzp1,Px, Py
CHARACTER*1 lf
CHARACTER*12 num1, num2, num3
CHARACTER*8 ctime
real*8 number,X,Clocx, Clocy, Dxl, Dxu, Dyl, Dyu
debug = 0
!
! Open File
!
Write(ctime,'(i8.8)') icycle
OPEN(15,FILE=trim(vtk_file)//'.'//ctime//'.vtk',FORM='unformatted', &
access='stream',convert='BIG_ENDIAN')
!
! Write header info
!
lf = char(10) ! line feed character
Write(15) "# vtk DataFile Version 2.0"//lf
Write(15) "EcoSLIM Points Output"//lf
Write(15) "BINARY"//lf
Write(15) "DATASET POLYDATA"//lf
write(num1, '(i12)') np_active
Write(15) "POINTS "//num1//" FLOAT"//lf
!write(15) ((real(P(j,i),kind=4), i=1,3), j=1,np_active) ! This forces the expected write order
do j =1, np_active
write(15) real(P(j,1:2),kind=4)
! find integer cell location
Px = floor(P(j,1) / dx)
Py = floor(P(j,2) / dy)
! Find each particle's factional cell location
Clocx = (P(j,1) - float(Px)*dx) / dx
Clocy = (P(j,2) - float(Py)*dy) / dy
! X is local adjustment for DEM
! linearly interpolate particle location in DEM cell
Dxl = DEM(Px+1,Py+1)
Dxu = DEM(Px+2,Py+1)
Dyl = DEM(Px+1,Py+1)
Dyu = DEM(Px+1,Py+2)
! correct for edges
if (Px == nx -1) Dxu = Dxl
if (Py == ny -1) Dyu = Dyl
X = ( ((1.0d0-Clocx)*Dxl &
+ Dxu*Clocx) + &
((1.0d0-Clocy)*Dyl &
+ Dyu*Clocy) ) / 2.0D0
! write new location
write(15) real(P(j,3)+X -maxZ,kind=4)
end do !!j
write(15) lf
write(15) "POINT_DATA "//num1//lf
write(15) "SCALARS Time float"//lf
Write(15) "LOOKUP_TABLE default"//lf
write(15) (real(P(j,4),kind=4), j=1,np_active)
write(15) lf
write(15) "SCALARS Mass float"//lf
Write(15) "LOOKUP_TABLE default"//lf
write(15) (real(P(j,6),kind=4), j=1,np_active)
write(15) lf
write(15) "SCALARS Source float"//lf
Write(15) "LOOKUP_TABLE default"//lf
write(15) (real(P(j,7),kind=4), j=1,np_active)
write(15) lf
write(15) "SCALARS ID float"//lf
Write(15) "LOOKUP_TABLE default"//lf
write(15) (real(P(j,11),kind=4), j=1,np_active)
write(15) lf
CLOSE(15)
RETURN
END SUBROUTINE vtk_write_points