386 lines
11 KiB
Fortran
386 lines
11 KiB
Fortran
MODULE m_compact
|
|
IMPLICIT NONE
|
|
PRIVATE
|
|
REAL*8, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
|
|
lyf,lys,wyf,wys, &
|
|
lzf,lzs,wzf,wzs
|
|
! lyzf,lyzs,wyzf,wyzs
|
|
INTEGER :: nxc,nyc,nzc
|
|
REAL*8, PARAMETER :: ezero = 1.0e-14
|
|
|
|
PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp
|
|
|
|
CONTAINS
|
|
|
|
SUBROUTINE ludcmp(nxx,nyy,nzz,xp,yp,zp)
|
|
INTEGER, INTENT(IN) :: nxx,nyy,nzz
|
|
INTEGER, INTENT(IN) :: xp,yp,zp
|
|
INTEGER :: ierr
|
|
|
|
nxc=nxx
|
|
nyc=nyy
|
|
nzc=nzz
|
|
! IF(nyc /= nzc) PRINT*,'ny should be equal nz'
|
|
|
|
! xp, yp, zp = 0 : periodic
|
|
ALLOCATE(lxf(nxc),STAT=ierr)
|
|
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
|
ALLOCATE(lxs(nxc),STAT=ierr)
|
|
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
|
IF(xp.eq.0) THEN
|
|
ALLOCATE(wxf(nxc),STAT=ierr)
|
|
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
|
ALLOCATE(wxs(nxc),STAT=ierr)
|
|
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
|
CALL p_lud(1,nxc)
|
|
ELSE
|
|
CALL nonp_lud(1,nxc)
|
|
ENDIF
|
|
|
|
ALLOCATE(lyf(nyc),STAT=ierr)
|
|
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
|
ALLOCATE(lys(nyc),STAT=ierr)
|
|
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
|
IF(yp.eq.0) THEN
|
|
ALLOCATE(wyf(nyc),STAT=ierr)
|
|
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
|
ALLOCATE(wys(nyc),STAT=ierr)
|
|
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
|
CALL p_lud(2,nyc)
|
|
ELSE
|
|
call nonp_lud(2,nyc)
|
|
ENDIF
|
|
|
|
ALLOCATE(lzf(nzc),STAT=ierr)
|
|
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
|
ALLOCATE(lzs(nzc),STAT=ierr)
|
|
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
|
IF(zp.eq.0) THEN
|
|
ALLOCATE(wzf(nzc),STAT=ierr)
|
|
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
|
ALLOCATE(wzs(nzc),STAT=ierr)
|
|
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
|
CALL p_lud(3,nzc)
|
|
ELSE
|
|
call nonp_lud(3,nzc)
|
|
ENDIF
|
|
|
|
! CALL x_lud
|
|
! CALL yz_lud
|
|
|
|
END SUBROUTINE ludcmp
|
|
|
|
SUBROUTINE nonp_lud(xyz,xx)
|
|
INTEGER :: i,xyz,xx
|
|
REAL*8, DIMENSION(xx) :: aa
|
|
aa=3.
|
|
aa(1)=0.5 ; aa(2)=4.
|
|
aa(xx-1)=4. ; aa(xx)=0.5
|
|
! first derivative
|
|
IF (xyz.eq.1) CALL stdlu(aa,xx,lxf) ! x-direction
|
|
IF (xyz.eq.2) CALL stdlu(aa,xx,lyf) ! y-direction
|
|
IF (xyz.eq.3) CALL stdlu(aa,xx,lzf) ! z-direction
|
|
aa=5.5
|
|
aa(1)=2./11. ; aa(2)=10.
|
|
aa(xx-1)=10. ; aa(xx)=2./11.
|
|
! second derivative
|
|
IF (xyz.eq.1) CALL stdlu(aa,xx,lxs) ! x-direction
|
|
IF (xyz.eq.2) CALL stdlu(aa,xx,lys) ! y-direction
|
|
IF (xyz.eq.3) CALL stdlu(aa,xx,lzs) ! z-direction
|
|
END SUBROUTINE nonp_lud
|
|
|
|
|
|
SUBROUTINE p_lud(xyz,xx)
|
|
INTEGER :: i,xyz,xx
|
|
REAL*8 :: a
|
|
a=3. ! first derivative
|
|
IF (xyz.eq.1) CALL ptdlu(a,xx,lxf,wxf) ! x-direction
|
|
IF (xyz.eq.2) CALL ptdlu(a,xx,lyf,wyf) ! y-direction
|
|
IF (xyz.eq.3) CALL ptdlu(a,xx,lzf,wzf) ! z-direction
|
|
a=11./2. ! second derivative
|
|
IF (xyz.eq.1) CALL ptdlu(a,xx,lxs,wxs) ! x-direction
|
|
IF (xyz.eq.2) CALL ptdlu(a,xx,lys,wys) ! y-direction
|
|
IF (xyz.eq.3) CALL ptdlu(a,xx,lzs,wzs) ! z-direction
|
|
END SUBROUTINE p_lud
|
|
|
|
SUBROUTINE stdlu(a,n,l)
|
|
INTEGER :: n
|
|
REAL*8 :: a(n),l(n)
|
|
REAL*8 :: d
|
|
INTEGER :: i
|
|
l(1)=1.0/a(1)
|
|
DO i=2,n
|
|
d=a(i)-l(i-1)
|
|
l(i)=1.0/d
|
|
ENDDO
|
|
END SUBROUTINE stdlu
|
|
|
|
SUBROUTINE ptdlu(a,n,l,w)
|
|
INTEGER :: n
|
|
REAL*8 :: a,l(n),w(n)
|
|
INTEGER :: i
|
|
REAL*8 :: aa(n),d
|
|
|
|
DO i=1,n-1
|
|
aa(i)=a
|
|
ENDDO
|
|
i=n-1
|
|
call stdlu(aa,i,l)
|
|
w(1)=1.0
|
|
DO i=2,n-2
|
|
w(i)=-l(i-1)*w(i-1)
|
|
ENDDO
|
|
w(n-1)=1.0-l(n-2)*w(n-2)
|
|
DO i=1,n-1
|
|
w(i)=w(i)*l(i)
|
|
ENDDO
|
|
d=a
|
|
DO i=1,n-1
|
|
d=d-w(i)*w(i)/l(i)
|
|
ENDDO
|
|
l(n)=1./d
|
|
END SUBROUTINE ptdlu
|
|
|
|
SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
|
|
INTEGER,INTENT(IN) :: n,nd,dir
|
|
REAL*8,INTENT(IN) :: h
|
|
REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
|
|
REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
|
|
INTEGER :: i,j
|
|
REAL*8 :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
|
|
|
|
h1=1./h
|
|
|
|
r1=7./3.
|
|
r2=1./12.
|
|
r3=3.
|
|
a=-1.25
|
|
b=1.
|
|
c=0.25
|
|
|
|
DO j=1,nd
|
|
dx(j,n-1)=x(j,n)-x(j,n-2)
|
|
dx(j,n)=-(a*x(j,n)+b*x(j,n-1)+c*x(j,n-2))
|
|
dx(j,1)=(a*x(j,1)+b*x(j,2)+c*x(j,3))
|
|
dx(j,2)=x(j,3)-x(j,1)
|
|
IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) dx(j,n)=0.
|
|
IF (x(j,1).eq.x(j,2).and.x(j,2).eq.x(j,3)) dx(j,1)=0.
|
|
dx(j,n-1)=dx(j,n-1)*h1*r3
|
|
dx(j,n)=dx(j,n)*h1
|
|
dx(j,1)=dx(j,1)*h1
|
|
dx(j,2)=dx(j,2)*h1*r3
|
|
ENDDO
|
|
|
|
DO i=3,n-2
|
|
DO j=1,nd
|
|
t1=x(j,i+1)-x(j,i-1)
|
|
t2=x(j,i+2)-x(j,i-2)
|
|
dx(j,i)=h1*(r1*t1+r2*t2)
|
|
ENDDO
|
|
ENDDO
|
|
IF (dir.eq.1) CALL tdslv(dx,n,lxf,nd) ! x-direction
|
|
IF (dir.eq.2) CALL tdslv(dx,n,lyf,nd) ! y-direction
|
|
IF (dir.eq.3) CALL tdslv(dx,n,lzf,nd) ! z-direction
|
|
END SUBROUTINE dfnonp
|
|
|
|
SUBROUTINE dfp(n,h,x,dx,nd,dir)
|
|
INTEGER,INTENT(IN) :: n,nd,dir
|
|
REAL*8,INTENT(IN) :: h
|
|
REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
|
|
REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
|
|
INTEGER :: i,j
|
|
REAL*8 :: r1,r2,h1
|
|
|
|
h1=1./h
|
|
r1=7./3.
|
|
r2=1./12.
|
|
|
|
DO j=1,nd
|
|
dx(j,n-1)=(r1*(x(j,n)-x(j,n-2))+r2*(x(j,1)-x(j,n-3)))
|
|
dx(j,n)=(r1*(x(j,1)-x(j,n-1))+r2*(x(j,2)-x(j,n-2)))
|
|
dx(j,1)=(r1*(x(j,2)-x(j,n))+r2*(x(j,3)-x(j,n-1)))
|
|
dx(j,2)=(r1*(x(j,3)-x(j,1))+r2*(x(j,4)-x(j,n)))
|
|
dx(j,n-1)=dx(j,n-1)*h1
|
|
dx(j,n)=dx(j,n)*h1
|
|
dx(j,1)=dx(j,1)*h1
|
|
dx(j,2)=dx(j,2)*h1
|
|
ENDDO
|
|
|
|
DO i=3,n-2
|
|
DO j=1,nd
|
|
dx(j,i)=(r1*(x(j,i+1)-x(j,i-1))+r2*(x(j,i+2)-x(j,i-2)))
|
|
dx(j,i)=dx(j,i)*h1
|
|
ENDDO
|
|
ENDDO
|
|
|
|
IF (dir.eq.1) CALL ptdslv(dx,n,lxf,wxf,nd) ! x-direction
|
|
IF (dir.eq.2) CALL ptdslv(dx,n,lyf,wyf,nd) ! y-direction
|
|
IF (dir.eq.3) CALL ptdslv(dx,n,lzf,wzf,nd) ! z-direction
|
|
|
|
END SUBROUTINE dfp
|
|
|
|
SUBROUTINE ptdslv(r,n,l,w,nd)
|
|
INTEGER,INTENT(IN) :: n,nd
|
|
REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
|
|
REAL*8,INTENT(IN),DIMENSION(:) :: l,w
|
|
INTEGER i,j
|
|
REAL*8, DIMENSION(nd) :: sum
|
|
DO j=1,nd
|
|
sum(j)=w(1)*r(j,1)
|
|
r(j,1)=r(j,1)*l(1)
|
|
ENDDO
|
|
DO i=2,n-1
|
|
DO j=1,nd
|
|
r(j,i)=r(j,i)-r(j,i-1)
|
|
sum(j)=sum(j)+w(i)*r(j,i)
|
|
r(j,i)=r(j,i)*l(i)
|
|
ENDDO
|
|
ENDDO
|
|
DO j=1,nd
|
|
r(j,n)=l(n)*(r(j,n)-sum(j))
|
|
r(j,n-1)=r(j,n-1)-w(n-1)*r(j,n)
|
|
ENDDO
|
|
DO i=n-2,1,-1
|
|
DO j=1,nd
|
|
r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r(j,n)
|
|
ENDDO
|
|
ENDDO
|
|
END SUBROUTINE ptdslv
|
|
|
|
SUBROUTINE d2fp(n,h,x,dx,nd,dir)
|
|
INTEGER,INTENT(IN) :: n,nd,dir
|
|
REAL*8,INTENT(IN) :: h
|
|
REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
|
|
REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
|
|
INTEGER :: i,j
|
|
REAL*8 :: h2,r1,r2,t1,t2
|
|
h2=1./(h*h)
|
|
r1=6.
|
|
r2=3./8.
|
|
DO j=1,nd
|
|
t1 = (x(j,n)-2.*x(j,n-1)+x(j,n-2))
|
|
t2 = (x(j,1)-2.*x(j,n-1)+x(j,n-3))
|
|
IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) t1=0.
|
|
IF (x(j,1).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-3)) t2=0.
|
|
dx(j,n-1)=(r1*t1+r2*t2)
|
|
|
|
t1 = (x(j,1)-2.*x(j,n)+x(j,n-1))
|
|
t2 = (x(j,2)-2.*x(j,n)+x(j,n-2))
|
|
IF (x(j,1).eq.x(j,n).and.x(j,n).eq.x(j,n-1)) t1=0.
|
|
IF (x(j,2).eq.x(j,n).and.x(j,n).eq.x(j,n-2)) t2=0.
|
|
! dx(j,n)=(r1*(x(j,1)-2.*x(j,n)+x(j,n-1)) &
|
|
! +r2*(x(j,2)-2.*x(j,n)+x(j,n-2)))
|
|
dx(j,n)=(r1*t1+r2*t2)
|
|
|
|
t1 = (x(j,2)-2.*x(j,1)+x(j,n))
|
|
t2 = (x(j,3)-2.*x(j,1)+x(j,n-1))
|
|
IF (x(j,2).eq.x(j,1).and.x(j,1).eq.x(j,n)) t1=0.
|
|
IF (x(j,3).eq.x(j,1).and.x(j,1).eq.x(j,n-1)) t2=0.
|
|
! dx(j,1)=(r1*(x(j,2)-2.*x(j,1)+x(j,n)) &
|
|
! +r2*(x(j,3)-2.*x(j,1)+x(j,n-1)))
|
|
dx(j,1)=(r1*t1+r2*t2)
|
|
|
|
t1 = (x(j,3)-2.*x(j,2)+x(j,1))
|
|
t2 = (x(j,4)-2.*x(j,2)+x(j,n))
|
|
IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) t1=0.
|
|
IF (x(j,4).eq.x(j,2).and.x(j,2).eq.x(j,n)) t2=0.
|
|
! dx(j,2)=(r1*(x(j,3)-2.*x(j,2)+x(j,1)) &
|
|
! +r2*(x(j,4)-2.*x(j,2)+x(j,n)))
|
|
dx(j,2)=(r1*t1+r2*t2)
|
|
|
|
dx(j,n-1)=dx(j,n-1)*h2
|
|
dx(j,n)=dx(j,n)*h2
|
|
dx(j,1)=dx(j,1)*h2
|
|
dx(j,2)=dx(j,2)*h2
|
|
ENDDO
|
|
DO i=3,n-2
|
|
DO j=1,nd
|
|
t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1))
|
|
t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2))
|
|
IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0.
|
|
IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0.
|
|
! dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) &
|
|
! +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2)))
|
|
dx(j,i)=(r1*t1+r2*t2)
|
|
dx(j,i)=dx(j,i)*h2
|
|
ENDDO
|
|
ENDDO
|
|
IF (dir.eq.1) CALL ptdslv(dx,n,lxs,wxs,nd) ! x-direction
|
|
IF (dir.eq.2) CALL ptdslv(dx,n,lys,wys,nd) ! y-direction
|
|
IF (dir.eq.3) CALL ptdslv(dx,n,lzs,wzs,nd) ! z-direction
|
|
END SUBROUTINE d2fp
|
|
|
|
SUBROUTINE tdslv(r,n,l,nd)
|
|
INTEGER,INTENT(IN) :: n,nd
|
|
REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
|
|
REAL*8,INTENT(IN),DIMENSION(:) :: l
|
|
INTEGER i,j
|
|
REAL*8 t1
|
|
DO j=1,nd
|
|
r(j,1)=r(j,1)*l(1)
|
|
ENDDO
|
|
DO i=2,n
|
|
DO j=1,nd
|
|
t1=r(j,i)-r(j,i-1)
|
|
r(j,i)=l(i)*t1
|
|
ENDDO
|
|
ENDDO
|
|
DO i=n-1,1,-1
|
|
DO j=1,nd
|
|
r(j,i)=r(j,i)-l(i)*r(j,i+1)
|
|
ENDDO
|
|
ENDDO
|
|
END SUBROUTINE tdslv
|
|
|
|
SUBROUTINE d2fnonp(n,h,x,dx,nd,dir)
|
|
INTEGER,INTENT(IN) :: n,nd,dir
|
|
REAL*8,INTENT(IN) :: h
|
|
REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
|
|
REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
|
|
INTEGER :: i,j
|
|
REAL*8 :: h2,r1,r2,r3,a,b,c,e,t1,t2
|
|
|
|
h2=1./(h*h)
|
|
r1=6.
|
|
r2=3./8.
|
|
r3=12.
|
|
a=13./11.
|
|
b=-27./11.
|
|
c=15./11.
|
|
e=-1./11.
|
|
|
|
DO j=1,nd
|
|
dx(j,1)=(a*x(j,1)+b*x(j,2)+c*x(j,3)+e*x(j,4))
|
|
dx(j,2)=(x(j,3)-2.*x(j,2)+x(j,1))
|
|
dx(j,n-1)=(x(j,n)-2.*x(j,n-1)+x(j,n-2))
|
|
dx(j,n)=(a*x(j,n)+b*x(j,n-1)+c*x(j,n-2)+e*x(j,n-3))
|
|
IF (x(j,1).eq.x(j,2).and.x(j,2).eq.x(j,3).and.x(j,3).eq.x(j,4)) dx(j,1)=0.
|
|
IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) dx(j,2)=0.
|
|
IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2).and.x(j,n-2).eq.x(j,n-3)) dx(j,n)=0.
|
|
IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) dx(j,n-1)=0.
|
|
dx(j,1)=dx(j,1)*h2
|
|
dx(j,2)=dx(j,2)*h2*r3
|
|
dx(j,n-1)=dx(j,n-1)*h2*r3
|
|
dx(j,n)=dx(j,n)*h2
|
|
ENDDO
|
|
DO i=3,n-2
|
|
DO j=1,nd
|
|
t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1))
|
|
t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2))
|
|
IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0.
|
|
IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0.
|
|
|
|
! dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) &
|
|
! +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2)))
|
|
dx(j,i)=(r1*t1+r2*t2)
|
|
dx(j,i)=dx(j,i)*h2
|
|
ENDDO
|
|
ENDDO
|
|
|
|
IF (dir.eq.1) CALL tdslv(dx,n,lxs,nd) ! x-direction
|
|
IF (dir.eq.2) CALL tdslv(dx,n,lys,nd) ! y-direction
|
|
IF (dir.eq.3) CALL tdslv(dx,n,lzs,nd) ! z-direction
|
|
|
|
END SUBROUTINE d2fnonp
|
|
|
|
END MODULE m_compact
|