Page 430 - Handbook Of Multiphase Flow Assurance

Computer code (Makogon, 1994, 1997) 429

x3(i)=x3(i)+0.5D0
C shift the block by sqrt(2)/2 unit cell to locate edge at y=0
y1(i)=y1(i)+0.5D0*DSQRT(2.D0)
y2(i)=y2(i)+0.5D0*DSQRT(2.D0)
y3(i)=y3(i)+0.5D0*DSQRT(2.D0)
end do
C end of slicing of the basic [110] block
endif

C correct the upper surface position
C top boundary is strict, i.e. z.GE.0. is shifted.
do i=1,npts
z1(i)=z1(i)+(diag-B)/unitside
z2(i)=z2(i)+(diag-B)/unitside
z3(i)=z3(i)+(diag-B)/unitside
if(z1(i).ge.0.D0) then
z1(i)=z1(i)-diag/unitside
z2(i)=z2(i)-diag/unitside
z3(i)=z3(i)-diag/unitside
endif
end do

C add more unit cells in height if needed
nadd=0
CC=C
do while (CC.gt.diag)
nadd=nadd+1
do i=1,npts
x1(i+nadd*npts)=x1(i)
y1(i+nadd*npts)=y1(i)
z1(i+nadd*npts)=z1(i)-diag/unitside*nadd
x2(i+nadd*npts)=x2(i)
y2(i+nadd*npts)=y2(i)
z2(i+nadd*npts)=z2(i)-diag/unitside*nadd
x3(i+nadd*npts)=x3(i)
y3(i+nadd*npts)=y3(i)
z3(i+nadd*npts)=z3(i)-diag/unitside*nadd
end do
CC=CC-diag
end do
npts=npts+npts*nadd

C Correct the slice thickness to requested value. Adjust lower surface.
C Lower boundary is relaxed i.e. z.LT.-C is removed
do i=1,npts