added m_calculate and m_parameters
This commit is contained in:
ignis
parent
cb68140f62
commit
b01be79870
4 changed files with 533 additions and 57 deletions
279
code/m_calculate.f90
Normal file
279
code/m_calculate.f90
Normal file
|
|
@ -0,0 +1,279 @@
|
|||
module m_calculate
|
||||
|
||||
use m_parameters
|
||||
use m_compact
|
||||
|
||||
implicit none
|
||||
|
||||
real*8, allocatable, dimension(:,:) :: xsrc
|
||||
real*8, allocatable, dimension(:,:) :: xdst
|
||||
|
||||
real*8, allocatable, dimension(:,:) :: rsrc
|
||||
real*8, allocatable, dimension(:,:) :: rdst
|
||||
|
||||
integer, parameter :: nb = BLOCKSIZE
|
||||
|
||||
private :: nb
|
||||
private :: xsrc, xdst, rsrc, rdst
|
||||
|
||||
contains
|
||||
|
||||
subroutine m_calculate_init
|
||||
|
||||
integer :: ierr
|
||||
|
||||
call ludcmp(nxp,nyp,nzp,1,0,0) ! 1,1,0
|
||||
|
||||
allocate(xsrc(nb, nxp), stat=ierr)
|
||||
allocate(xdst(nb, nxp), stat=ierr)
|
||||
|
||||
allocate(rsrc(nxp, nzp), stat=ierr)
|
||||
allocate(rdst(nxp, nzp), stat=ierr)
|
||||
|
||||
end subroutine m_calculate_init
|
||||
|
||||
|
||||
subroutine m_calculate_finalize
|
||||
|
||||
deallocate(xsrc)
|
||||
deallocate(xdst)
|
||||
|
||||
deallocate(rsrc)
|
||||
deallocate(rdst)
|
||||
|
||||
end subroutine m_calculate_finalize
|
||||
|
||||
|
||||
subroutine ddx(dst, src)
|
||||
|
||||
real*8, dimension(nxp,nyp,nzp), intent(in) :: src
|
||||
real*8, dimension(nxp,nyp,nzp), intent(out) :: dst
|
||||
|
||||
integer :: i, j ,k
|
||||
integer :: ju
|
||||
|
||||
|
||||
do k = 1,nzp
|
||||
do j = 1,nyp,nb
|
||||
|
||||
ju = min(j+nb-1,nyp)
|
||||
|
||||
call tp2(xsrc, src(:,j:ju,k), nb, nxp)
|
||||
|
||||
call dfnonp(nxp, hxp, xsrc, xdst, nb, 1)
|
||||
|
||||
call tp2(dst(:,j:ju,k), xdst, nxp, nb)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end subroutine ddx
|
||||
|
||||
|
||||
subroutine ddy(dst, src)
|
||||
|
||||
real*8, dimension(nxp,nyp,nzp), intent(in) :: src
|
||||
real*8, dimension(nxp,nyp,nzp), intent(out) :: dst
|
||||
|
||||
integer :: i, j ,k
|
||||
|
||||
|
||||
do k = 1,nzp
|
||||
|
||||
call dfp(nyp, hyp, src(:,:,k), dst(:,:,k), nxp, 2)
|
||||
|
||||
end do
|
||||
|
||||
end subroutine ddy
|
||||
|
||||
|
||||
subroutine ddz(dst, src)
|
||||
|
||||
real*8, dimension(nxp,nyp,nzp), intent(in) :: src
|
||||
real*8, dimension(nxp,nyp,nzp), intent(out) :: dst
|
||||
|
||||
integer :: i, j ,k
|
||||
|
||||
do j = 1,nyp
|
||||
|
||||
do k = 1,nzp
|
||||
rsrc(:,k) = src(:,j,k)
|
||||
end do
|
||||
|
||||
call dfp(nzp, hzp, rsrc, rdst, nxp, 3)
|
||||
|
||||
do k = 1,nzp
|
||||
dst(:,j,k) = rdst(:,k)
|
||||
end do
|
||||
|
||||
end do
|
||||
|
||||
end subroutine ddz
|
||||
|
||||
|
||||
subroutine d2dx(dst, src)
|
||||
|
||||
real*8, dimension(nxp,nyp,nzp), intent(in) :: src
|
||||
real*8, dimension(nxp,nyp,nzp), intent(out) :: dst
|
||||
|
||||
integer :: i, j ,k
|
||||
integer :: ju
|
||||
|
||||
|
||||
do k = 1,nzp
|
||||
do j = 1,nyp,nb
|
||||
|
||||
ju = min(j+nb-1,nyp)
|
||||
|
||||
call tp2(xsrc, src(:,j:ju,k), nb, nxp)
|
||||
|
||||
call d2fnonp(nxp, hxp, xsrc, xdst, nb, 1)
|
||||
|
||||
call tp2(dst(:,j:ju,k), xdst, nxp, nb)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end subroutine d2dx
|
||||
|
||||
|
||||
subroutine d2dy(dst, src)
|
||||
|
||||
real*8, dimension(nxp,nyp,nzp), intent(in) :: src
|
||||
real*8, dimension(nxp,nyp,nzp), intent(out) :: dst
|
||||
|
||||
integer :: i, j ,k
|
||||
|
||||
|
||||
do k = 1,nzp
|
||||
|
||||
call d2fp(nyp, hyp, src(:,:,k), dst(:,:,k), nxp, 2)
|
||||
|
||||
end do
|
||||
|
||||
end subroutine d2dy
|
||||
|
||||
|
||||
subroutine d2dz(dst, src)
|
||||
|
||||
real*8, dimension(nxp,nyp,nzp), intent(in) :: src
|
||||
real*8, dimension(nxp,nyp,nzp), intent(out) :: dst
|
||||
|
||||
integer :: i, j ,k
|
||||
|
||||
do j = 1,nyp
|
||||
|
||||
do k = 1,nzp
|
||||
rsrc(:,k) = src(:,j,k)
|
||||
end do
|
||||
|
||||
call d2fp(nzp, hzp, rsrc, rdst, nxp, 3)
|
||||
|
||||
do k = 1,nzp
|
||||
dst(:,j,k) = rdst(:,k)
|
||||
end do
|
||||
|
||||
end do
|
||||
|
||||
end subroutine d2dz
|
||||
|
||||
|
||||
subroutine tp(a, b, nx)
|
||||
! a(nb,nx) = transpose(b(nx,nb))
|
||||
|
||||
integer,intent(in) :: nx
|
||||
|
||||
real*8,intent(out) :: a(nb,nx)
|
||||
real*8,intent(in) :: b(nx,nb)
|
||||
|
||||
call tp2(a, b, nb, nx)
|
||||
|
||||
end subroutine tp
|
||||
|
||||
subroutine tp2 (a, b, n1, n2)
|
||||
! a = transpose(b)
|
||||
|
||||
implicit none
|
||||
|
||||
integer,intent(in) :: n1, n2
|
||||
real*8,intent(out) :: a(n1,n2)
|
||||
real*8,intent(in) :: b(n2,n1)
|
||||
integer :: i,j,ii,jj
|
||||
|
||||
DO jj=1,n2,nb
|
||||
DO ii=1,n1,nb
|
||||
|
||||
DO j=jj,min(jj+nb-1,n2)
|
||||
DO i=ii,min(ii+nb-1,n1)
|
||||
a(i,j) = b(j,i)
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
end subroutine tp2
|
||||
|
||||
real function rxn_rate (c)
|
||||
|
||||
real :: c
|
||||
|
||||
if(c.lt.0.) c=0.
|
||||
if(c.gt.1.) c=1.
|
||||
|
||||
if (c.le.c_cut) then
|
||||
|
||||
rxn_rate = min_wr
|
||||
|
||||
else if (c.gt.c_ref) then
|
||||
|
||||
rxn_rate = pre*(1.-c)*exp(-ac/(1.+bc*c))
|
||||
|
||||
else
|
||||
|
||||
rxn_rate = &
|
||||
((refwr-min_wr)*exp(prof_wr*(c-c_ref)) + min_wr - refwr*exp(prof_wr*(c_cut-c_ref))) &
|
||||
/ (1.-exp(prof_wr*(c_cut-c_ref)))
|
||||
|
||||
endif
|
||||
|
||||
end function rxn_rate
|
||||
|
||||
real function threshold_min_max (c, minc, maxc)
|
||||
|
||||
real :: c
|
||||
real :: minc, maxc
|
||||
|
||||
if ((c.lt.minc) .or. (c.gt.maxc)) then
|
||||
threshold_min_max = 0.
|
||||
else
|
||||
threshold_min_max = 1.0
|
||||
end if
|
||||
|
||||
end function threshold_min_max
|
||||
|
||||
real function positive (c)
|
||||
|
||||
real :: c
|
||||
|
||||
if (c > 0.0d0) then
|
||||
positive = c
|
||||
else
|
||||
positive = 0.
|
||||
end if
|
||||
|
||||
end function positive
|
||||
|
||||
real function negative (c)
|
||||
|
||||
real :: c
|
||||
|
||||
if (c < 0.0d0) then
|
||||
negative = c
|
||||
else
|
||||
negative = 0.
|
||||
end if
|
||||
|
||||
end function negative
|
||||
|
||||
end module m_calculate
|
||||
|
|
@ -1,14 +1,13 @@
|
|||
MODULE m_compact
|
||||
IMPLICIT NONE
|
||||
PRIVATE
|
||||
REAL, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
|
||||
lyf,lys,wyf,wys, &
|
||||
lzf,lzs,wzf,wzs
|
||||
! lyzf,lyzs,wyzf,wyzs
|
||||
|
||||
REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
|
||||
lyf,lys,wyf,wys, &
|
||||
lzf,lzs,wzf,wzs
|
||||
|
||||
INTEGER :: nxc,nyc,nzc
|
||||
REAL, PARAMETER :: ezero = 1.0e-14
|
||||
|
||||
PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp
|
||||
|
||||
REAL(KIND=8), PARAMETER :: ezero = 1.0e-14
|
||||
|
||||
CONTAINS
|
||||
|
||||
|
|
@ -20,6 +19,22 @@
|
|||
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
|
||||
|
|
@ -32,9 +47,6 @@
|
|||
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)
|
||||
|
|
@ -46,9 +58,6 @@
|
|||
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)
|
||||
|
|
@ -60,19 +69,140 @@
|
|||
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'
|
||||
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
|
||||
|
||||
! CALL x_lud
|
||||
! CALL yz_lud
|
||||
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
|
||||
|
||||
END SUBROUTINE ludcmp
|
||||
|
||||
SUBROUTINE nonp_lud(xyz,xx)
|
||||
INTEGER :: i,xyz,xx
|
||||
REAL, DIMENSION(xx) :: aa
|
||||
REAL(KIND=8), DIMENSION(xx) :: aa
|
||||
aa=3.
|
||||
aa(1)=0.5 ; aa(2)=4.
|
||||
aa(xx-1)=4. ; aa(xx)=0.5
|
||||
|
|
@ -92,7 +222,7 @@
|
|||
|
||||
SUBROUTINE p_lud(xyz,xx)
|
||||
INTEGER :: i,xyz,xx
|
||||
REAL :: a
|
||||
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
|
||||
|
|
@ -105,21 +235,23 @@
|
|||
|
||||
SUBROUTINE stdlu(a,n,l)
|
||||
INTEGER :: n
|
||||
REAL :: a(n),l(n)
|
||||
REAL :: d
|
||||
REAL(KIND=8), INTENT(IN) :: a(n)
|
||||
REAL(KIND=8), INTENT(OUT) :: l(n)
|
||||
REAL(KIND=8) :: d
|
||||
INTEGER :: i
|
||||
l(1)=1.0/a(1)
|
||||
l(1)=1.0d0/a(1)
|
||||
DO i=2,n
|
||||
d=a(i)-l(i-1)
|
||||
l(i)=1.0/d
|
||||
l(i)=1.0d0/d
|
||||
ENDDO
|
||||
END SUBROUTINE stdlu
|
||||
|
||||
SUBROUTINE ptdlu(a,n,l,w)
|
||||
INTEGER :: n
|
||||
REAL :: a,l(n),w(n)
|
||||
REAL(KIND=8), INTENT(OUT) :: a
|
||||
REAL(KIND=8), INTENT(OUT) :: l(n),w(n)
|
||||
INTEGER :: i
|
||||
REAL :: aa(n),d
|
||||
REAL(KIND=8) :: aa(n),d
|
||||
|
||||
DO i=1,n-1
|
||||
aa(i)=a
|
||||
|
|
@ -141,18 +273,19 @@
|
|||
l(n)=1./d
|
||||
END SUBROUTINE ptdlu
|
||||
|
||||
SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
|
||||
INTEGER,INTENT(IN) :: n,nd,dir
|
||||
REAL,INTENT(IN) :: h
|
||||
REAL,INTENT(IN),DIMENSION(nd,n) :: x
|
||||
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
|
||||
|
||||
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 :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
|
||||
REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
|
||||
|
||||
h1=1./h
|
||||
h1=1.d0/h
|
||||
|
||||
r1=7./3.
|
||||
r2=1./12.
|
||||
r1=7.d0/3.d0
|
||||
r2=1.d0/12.d0
|
||||
r3=3.
|
||||
a=-1.25
|
||||
b=1.
|
||||
|
|
@ -178,6 +311,20 @@
|
|||
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
|
||||
|
|
@ -185,11 +332,14 @@
|
|||
|
||||
SUBROUTINE dfp(n,h,x,dx,nd,dir)
|
||||
INTEGER,INTENT(IN) :: n,nd,dir
|
||||
REAL,INTENT(IN) :: h
|
||||
REAL,INTENT(IN),DIMENSION(nd,n) :: x
|
||||
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
|
||||
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 :: r1,r2,h1
|
||||
REAL(KIND=8) :: r1,r2,h1
|
||||
|
||||
|
||||
! print *, "dfnonp received (nd,n)", nd, n
|
||||
|
||||
h1=1./h
|
||||
r1=7./3.
|
||||
|
|
@ -221,10 +371,10 @@
|
|||
|
||||
SUBROUTINE ptdslv(r,n,l,w,nd)
|
||||
INTEGER,INTENT(IN) :: n,nd
|
||||
REAL,INTENT(INOUT),DIMENSION(nd,n) :: r
|
||||
REAL,INTENT(IN),DIMENSION(:) :: l,w
|
||||
REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r
|
||||
REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l,w
|
||||
INTEGER i,j
|
||||
REAL, DIMENSION(nd) :: sum
|
||||
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)
|
||||
|
|
@ -249,11 +399,13 @@
|
|||
|
||||
SUBROUTINE d2fp(n,h,x,dx,nd,dir)
|
||||
INTEGER,INTENT(IN) :: n,nd,dir
|
||||
REAL,INTENT(IN) :: h
|
||||
REAL,INTENT(IN),DIMENSION(nd,n) :: x
|
||||
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
|
||||
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 :: h2,r1,r2,t1,t2
|
||||
REAL(KIND=8) :: h2,r1,r2,t1,t2
|
||||
|
||||
|
||||
h2=1./(h*h)
|
||||
r1=6.
|
||||
r2=3./8.
|
||||
|
|
@ -312,10 +464,10 @@
|
|||
|
||||
SUBROUTINE tdslv(r,n,l,nd)
|
||||
INTEGER,INTENT(IN) :: n,nd
|
||||
REAL,INTENT(INOUT),DIMENSION(nd,n) :: r
|
||||
REAL,INTENT(IN),DIMENSION(:) :: l
|
||||
REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r
|
||||
REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l
|
||||
INTEGER i,j
|
||||
REAL t1
|
||||
REAL(KIND=8) t1
|
||||
DO j=1,nd
|
||||
r(j,1)=r(j,1)*l(1)
|
||||
ENDDO
|
||||
|
|
@ -334,11 +486,12 @@
|
|||
|
||||
SUBROUTINE d2fnonp(n,h,x,dx,nd,dir)
|
||||
INTEGER,INTENT(IN) :: n,nd,dir
|
||||
REAL,INTENT(IN) :: h
|
||||
REAL,INTENT(IN),DIMENSION(nd,n) :: x
|
||||
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
|
||||
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 :: h2,r1,r2,r3,a,b,c,e,t1,t2
|
||||
REAL(KIND=8) :: h2,r1,r2,r3,a,b,c,e,t1,t2
|
||||
|
||||
|
||||
h2=1./(h*h)
|
||||
r1=6.
|
||||
|
|
|
|||
36
code/m_parameters.f90
Normal file
36
code/m_parameters.f90
Normal file
|
|
@ -0,0 +1,36 @@
|
|||
module m_parameters
|
||||
|
||||
implicit none
|
||||
|
||||
integer :: nxp,nyp,nzp
|
||||
|
||||
real :: hxp,hyp,hzp
|
||||
|
||||
real :: l_0
|
||||
|
||||
integer :: startnum,endnum,skipnum,shiftnum
|
||||
|
||||
real :: scp,prp,lep,vis0p,rod
|
||||
|
||||
real :: prof_wr,min_wr,min_fsd,min_c,refwr
|
||||
|
||||
real :: pre,ac,bc,c_cut,c_ref
|
||||
|
||||
integer :: syp,eyp,twod
|
||||
|
||||
real :: SL_u
|
||||
|
||||
integer :: omitnum
|
||||
|
||||
integer, allocatable :: omit_t(:,:)
|
||||
|
||||
integer, allocatable :: file_dist(:)
|
||||
|
||||
integer, allocatable :: export_offset(:)
|
||||
|
||||
real, parameter :: pi=3.14159265358979323846
|
||||
real, parameter :: me=1.00e-20
|
||||
|
||||
contains
|
||||
|
||||
end module m_parameters
|
||||
|
|
@ -1,4 +1,6 @@
|
|||
flags = -Wall -O3 -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 -march=native
|
||||
BLOCKSIZE?=32
|
||||
|
||||
flags = -Wall -O3 -cpp -DBLOCKSIZE=$(BLOCKSIZE) -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 -march=native
|
||||
|
||||
ifdef CHECK
|
||||
flags += -fcheck=all
|
||||
|
|
@ -10,18 +12,24 @@ endif
|
|||
|
||||
compiler = gfortran
|
||||
|
||||
ex : test.o ysolve.o m_compact.o
|
||||
ex : test.o ysolve.o m_parameters.o m_calculate.o m_compact.o
|
||||
${compiler} -o ex test.o ysolve.o m_compact.o
|
||||
|
||||
test.o : test.f90 ysolve.mod compact.mod
|
||||
test.o : test.f90 ysolve.mod m_compact.mod
|
||||
${compiler} -c ${flags} test.f90
|
||||
|
||||
ysolve.o ysolve.mod : ysolve.f90 compact.mod
|
||||
ysolve.o ysolve.mod : ysolve.f90 m_compact.mod
|
||||
${compiler} -c ${flags} ysolve.f90
|
||||
|
||||
m_compact.o compact.mod : m_compact.f90
|
||||
m_compact.o m_compact.mod : m_compact.f90
|
||||
${compiler} -c ${flags} m_compact.f90
|
||||
|
||||
m_calculate.o m_calculate.mod : m_calculate.f90
|
||||
${compiler} -c ${flags} m_calculate.f90
|
||||
|
||||
m_parameters.o m_parameters.mod : m_parameters.f90
|
||||
${compiler} -c ${flags} m_parameters.f90
|
||||
|
||||
test: ex
|
||||
sh test.sh
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue