576 lines
16 KiB
Fortran
576 lines
16 KiB
Fortran
MODULE Compact
|
|
use, intrinsic :: iso_fortran_env, only: real64
|
|
IMPLICIT NONE
|
|
|
|
REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
|
|
lyf,lys,wyf,wys, &
|
|
lzf,lzs,wzf,wzs
|
|
|
|
INTEGER :: nxc,nyc,nzc
|
|
|
|
REAL(KIND=8), PARAMETER :: ezero = 1.0e-14
|
|
|
|
CONTAINS
|
|
|
|
SUBROUTINE ludcmp(nx,ny,nz,xp,yp,zp)
|
|
INTEGER, INTENT(IN) :: nx,ny,nz
|
|
INTEGER, INTENT(IN) :: xp,yp,zp
|
|
INTEGER :: ierr
|
|
|
|
nxc=nx
|
|
nyc=ny
|
|
nzc=nz
|
|
|
|
CALL ludcmp_allocate(nx,ny,nz,xp,yp,zp)
|
|
|
|
CALL ludcmp_calculate(nx,ny,nz,xp,yp,zp)
|
|
|
|
END SUBROUTINE ludcmp
|
|
|
|
SUBROUTINE ludcmp_allocate(nx,ny,nz,xp,yp,zp)
|
|
INTEGER, INTENT(IN) :: nx,ny,nz
|
|
INTEGER, INTENT(IN) :: xp,yp,zp
|
|
INTEGER :: ierr
|
|
|
|
nxc=nx
|
|
nyc=ny
|
|
nzc=nz
|
|
|
|
! IF(nyc /= nzc) PRINT*,'ny should be equal nz'
|
|
|
|
! xp, yp, zp = 0 : periodic
|
|
ALLOCATE(lxf(nxc),STAT=ierr)
|
|
IF(ierr /= 0) THEN
|
|
PRINT*, 'work array for lud allocation failed'
|
|
STOP 1
|
|
ENDIF
|
|
ALLOCATE(lxs(nxc),STAT=ierr)
|
|
IF(ierr /= 0) THEN
|
|
PRINT*, 'work array for lud allocation failed'
|
|
STOP 1
|
|
ENDIF
|
|
IF(xp.eq.0) THEN
|
|
ALLOCATE(wxf(nxc),STAT=ierr)
|
|
IF(ierr /= 0) THEN
|
|
PRINT*, 'work array for lud allocation failed'
|
|
STOP 1
|
|
ENDIF
|
|
ALLOCATE(wxs(nxc),STAT=ierr)
|
|
IF(ierr /= 0) THEN
|
|
PRINT*, 'work array for lud allocation failed'
|
|
STOP 1
|
|
ENDIF
|
|
ENDIF
|
|
|
|
ALLOCATE(lyf(nyc),STAT=ierr)
|
|
IF(ierr /= 0) THEN
|
|
PRINT*, 'work array for lud allocation failed'
|
|
STOP 1
|
|
ENDIF
|
|
ALLOCATE(lys(nyc),STAT=ierr)
|
|
IF(ierr /= 0) THEN
|
|
PRINT*, 'work array for lud allocation failed'
|
|
STOP 1
|
|
ENDIF
|
|
IF(yp.eq.0) THEN
|
|
ALLOCATE(wyf(nyc),STAT=ierr)
|
|
IF(ierr /= 0) THEN
|
|
PRINT*, 'work array for lud allocation failed'
|
|
STOP 1
|
|
ENDIF
|
|
ALLOCATE(wys(nyc),STAT=ierr)
|
|
IF(ierr /= 0) THEN
|
|
PRINT*, 'work array for lud allocation failed'
|
|
STOP 1
|
|
ENDIF
|
|
ENDIF
|
|
|
|
ALLOCATE(lzf(nzc),STAT=ierr)
|
|
IF(ierr /= 0) THEN
|
|
PRINT*, 'work array for lud allocation failed'
|
|
STOP 1
|
|
ENDIF
|
|
ALLOCATE(lzs(nzc),STAT=ierr)
|
|
IF(ierr /= 0) THEN
|
|
PRINT*, 'work array for lud allocation failed'
|
|
STOP 1
|
|
ENDIF
|
|
IF(zp.eq.0) THEN
|
|
ALLOCATE(wzf(nzc),STAT=ierr)
|
|
IF(ierr /= 0) THEN
|
|
PRINT*, 'work array for lud allocation failed'
|
|
STOP 1
|
|
ENDIF
|
|
ALLOCATE(wzs(nzc),STAT=ierr)
|
|
IF(ierr /= 0) THEN
|
|
PRINT*, 'work array for lud allocation failed'
|
|
STOP 1
|
|
ENDIF
|
|
ENDIF
|
|
|
|
END SUBROUTINE ludcmp_allocate
|
|
|
|
SUBROUTINE ludcmp_deallocate(xp,yp,zp)
|
|
INTEGER, INTENT(IN) :: xp,yp,zp
|
|
|
|
! IF(nyc /= nzc) PRINT*,'ny should be equal nz'
|
|
|
|
! xp, yp, zp = 0 : periodic
|
|
DEALLOCATE(lxf)
|
|
DEALLOCATE(lxs)
|
|
IF(xp.eq.0) THEN
|
|
DEALLOCATE(wxf)
|
|
DEALLOCATE(wxs)
|
|
ENDIF
|
|
|
|
DEALLOCATE(lyf)
|
|
DEALLOCATE(lys)
|
|
IF(yp.eq.0) THEN
|
|
DEALLOCATE(wyf)
|
|
DEALLOCATE(wys)
|
|
ENDIF
|
|
|
|
DEALLOCATE(lzf)
|
|
DEALLOCATE(lzs)
|
|
IF(zp.eq.0) THEN
|
|
DEALLOCATE(wzf)
|
|
DEALLOCATE(wzs)
|
|
ENDIF
|
|
|
|
END SUBROUTINE ludcmp_deallocate
|
|
|
|
SUBROUTINE ludcmp_testalloc
|
|
|
|
IF (.not. ALLOCATED(lxf)) print *, "lxf not allocated"
|
|
IF (.not. ALLOCATED(lxs)) print *, "lxs not allocated"
|
|
IF (.not. ALLOCATED(wxf)) print *, "wxf not allocated"
|
|
IF (.not. ALLOCATED(wxs)) print *, "wxs not allocated"
|
|
|
|
IF (.not. ALLOCATED(lyf)) print *, "lyf not allocated"
|
|
IF (.not. ALLOCATED(lys)) print *, "lys not allocated"
|
|
IF (.not. ALLOCATED(wyf)) print *, "wyf not allocated"
|
|
IF (.not. ALLOCATED(wys)) print *, "wys not allocated"
|
|
|
|
IF (.not. ALLOCATED(lzf)) print *, "lzf not allocated"
|
|
IF (.not. ALLOCATED(lzs)) print *, "lzs not allocated"
|
|
IF (.not. ALLOCATED(wzf)) print *, "wzf not allocated"
|
|
IF (.not. ALLOCATED(wzs)) print *, "wzs not allocated"
|
|
|
|
END SUBROUTINE ludcmp_testalloc
|
|
|
|
SUBROUTINE ludcmp_calculate(nx,ny,nz,xp,yp,zp)
|
|
INTEGER, INTENT(IN) :: nx,ny,nz
|
|
INTEGER, INTENT(IN) :: xp,yp,zp
|
|
INTEGER :: ierr
|
|
|
|
nxc=nx
|
|
nyc=ny
|
|
nzc=nz
|
|
|
|
! CALL ludcmp_testalloc
|
|
|
|
! IF(nyc /= nzc) PRINT*,'ny should be equal nz'
|
|
|
|
! xp, yp, zp = 0 : periodic
|
|
IF(xp.eq.0) THEN
|
|
CALL p_lud(1,nxc)
|
|
ELSE
|
|
CALL nonp_lud(1,nxc)
|
|
ENDIF
|
|
|
|
IF(yp.eq.0) THEN
|
|
CALL p_lud(2,nyc)
|
|
ELSE
|
|
call nonp_lud(2,nyc)
|
|
ENDIF
|
|
|
|
IF(zp.eq.0) THEN
|
|
CALL p_lud(3,nzc)
|
|
ELSE
|
|
call nonp_lud(3,nzc)
|
|
ENDIF
|
|
|
|
END SUBROUTINE ludcmp_calculate
|
|
|
|
|
|
SUBROUTINE test_nonp_lud1(xx, coef)
|
|
INTEGER :: xx
|
|
REAL(KIND=8), DIMENSION(xx) :: aa
|
|
REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef
|
|
aa=3.
|
|
aa(1)=0.5 ; aa(2)=4.
|
|
aa(xx-1)=4. ; aa(xx)=0.5
|
|
|
|
CALL stdlu(aa,xx,coef)
|
|
END SUBROUTINE test_nonp_lud1
|
|
|
|
|
|
SUBROUTINE test_nonp_lud2(xx, coef)
|
|
INTEGER :: xx
|
|
REAL(KIND=8), DIMENSION(xx) :: aa
|
|
REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef
|
|
aa=5.5
|
|
aa(1)=2./11. ; aa(2)=10.
|
|
aa(xx-1)=10. ; aa(xx)=2./11.
|
|
|
|
CALL stdlu(aa,xx,coef)
|
|
END SUBROUTINE test_nonp_lud2
|
|
|
|
|
|
SUBROUTINE test_p_lud1(xx, coef1, coef2)
|
|
INTEGER :: xx
|
|
REAL(KIND=8) :: a
|
|
REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef1, coef2
|
|
a=3. ! first derivative
|
|
|
|
CALL ptdlu(a,xx,coef1,coef2) ! x-direction
|
|
END SUBROUTINE test_p_lud1
|
|
|
|
|
|
SUBROUTINE test_p_lud2(xx, coef1, coef2)
|
|
INTEGER :: xx
|
|
REAL(KIND=8) :: a
|
|
REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef1, coef2
|
|
a=11./2. ! second derivative
|
|
|
|
CALL ptdlu(a,xx,coef1,coef2) ! x-direction
|
|
END SUBROUTINE test_p_lud2
|
|
|
|
|
|
SUBROUTINE nonp_lud(xyz,xx)
|
|
INTEGER :: i,xyz,xx
|
|
REAL(KIND=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(KIND=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(KIND=8), INTENT(IN) :: a(n)
|
|
REAL(KIND=8), INTENT(OUT) :: l(n)
|
|
REAL(KIND=8) :: d
|
|
INTEGER :: i
|
|
l(1)=1.0d0/a(1)
|
|
DO i=2,n
|
|
d=a(i)-l(i-1)
|
|
l(i)=1.0d0/d
|
|
ENDDO
|
|
END SUBROUTINE stdlu
|
|
|
|
SUBROUTINE ptdlu(a,n,l,w)
|
|
INTEGER :: n
|
|
REAL(KIND=8), INTENT(IN) :: a
|
|
REAL(KIND=8), INTENT(OUT) :: l(n),w(n)
|
|
INTEGER :: i
|
|
REAL(KIND=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 rhs1np(n,h,x,dx,nd)
|
|
INTEGER,INTENT(IN) :: n,nd
|
|
REAL(KIND=8),INTENT(IN) :: h
|
|
REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
|
|
REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
|
|
INTEGER :: i,j
|
|
REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
|
|
|
|
h1=1.d0/h
|
|
|
|
r1=7.d0/3.d0
|
|
r2=1.d0/12.d0
|
|
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
|
|
|
|
END SUBROUTINE rhs1np
|
|
|
|
|
|
SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
|
|
INTEGER,INTENT(IN) :: n,nd,dir
|
|
REAL(KIND=8),INTENT(IN) :: h
|
|
REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
|
|
REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
|
|
INTEGER :: i,j
|
|
REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
|
|
|
|
CALL rhs1np (n,h,x,dx,nd)
|
|
|
|
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(KIND=8),INTENT(IN) :: h
|
|
REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
|
|
REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
|
|
INTEGER :: i,j
|
|
REAL(KIND=8) :: r1,r2,h1
|
|
|
|
|
|
! print *, "dfnonp received (nd,n)", nd, n
|
|
|
|
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(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r
|
|
REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l,w
|
|
INTEGER i,j
|
|
REAL(KIND=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(KIND=8),INTENT(IN) :: h
|
|
REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
|
|
REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
|
|
INTEGER :: i,j
|
|
REAL(KIND=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(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r
|
|
REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l
|
|
INTEGER i,j
|
|
REAL(KIND=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(KIND=8),INTENT(IN) :: h
|
|
REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
|
|
REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
|
|
INTEGER :: i,j
|
|
REAL(KIND=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 Compact
|