297 lines
5.5 KiB
Fortran
297 lines
5.5 KiB
Fortran
module m_calculate
|
|
|
|
use Compact
|
|
use m_parameters
|
|
|
|
implicit none
|
|
|
|
real*8, allocatable :: work(:,:,:,:)
|
|
|
|
real*8, allocatable, dimension(:,:) :: xsrc
|
|
real*8, allocatable, dimension(:,:) :: xdst
|
|
|
|
real*8, allocatable, dimension(:,:) :: rsrc
|
|
real*8, allocatable, dimension(:,:) :: rdst
|
|
|
|
real*8, allocatable, dimension(:,:,:) :: zsrc
|
|
real*8, allocatable, dimension(:,:,:) :: zdst
|
|
|
|
integer, parameter :: nb = 16
|
|
|
|
private :: work, nb
|
|
private :: xsrc, xdst, zsrc, zdst
|
|
|
|
contains
|
|
|
|
subroutine m_calculate_init
|
|
|
|
integer :: ierr
|
|
|
|
! allocate(work(nxp, nyp, nzp, 3), stat=ierr)
|
|
|
|
call ludcmp(nxp,nyp,nzp,1,0,0) ! 1,1,0
|
|
|
|
allocate(xsrc(nyp, nxp), stat=ierr)
|
|
allocate(xdst(nyp, nxp), stat=ierr)
|
|
|
|
allocate(rsrc(nxp, nzp), stat=ierr)
|
|
allocate(rdst(nxp, nzp), stat=ierr)
|
|
|
|
allocate(zsrc(nyp, nzp, nxp), stat=ierr)
|
|
allocate(zdst(nyp, nzp, nxp), stat=ierr)
|
|
|
|
end subroutine m_calculate_init
|
|
|
|
|
|
subroutine m_calculate_finalize
|
|
|
|
! deallocate(work)
|
|
|
|
deallocate(xsrc)
|
|
deallocate(xdst)
|
|
|
|
deallocate(rsrc)
|
|
deallocate(rdst)
|
|
|
|
deallocate(zsrc)
|
|
deallocate(zdst)
|
|
|
|
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
|
|
|
|
|
|
do k = 1,nzp
|
|
|
|
call tp2(xsrc, src(:,:,k), nyp, nxp)
|
|
|
|
call dfnonp(nxp, hxp, xsrc, xdst, nyp, 1)
|
|
|
|
call tp2(dst(:,:,k), xdst, nxp, nyp)
|
|
|
|
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
|
|
|
|
|
|
call tp2(zsrc, src, nyp*nzp, nxp)
|
|
|
|
do i = 1,nxp
|
|
|
|
call dfp(nzp, hzp, zsrc(:,:,i), zdst(:,:,i), nyp, 3)
|
|
|
|
end do
|
|
|
|
call tp2(dst, zdst, nxp, nyp*nzp)
|
|
|
|
end subroutine ddz
|
|
|
|
|
|
subroutine ddz_ref(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
|
|
do i = 1,nxp
|
|
rsrc(i,k) = src(i,j,k)
|
|
end do
|
|
end do
|
|
|
|
call dfp(nzp, hzp, rsrc, rdst, nxp, 3)
|
|
|
|
do k = 1,nzp
|
|
do i = 1,nxp
|
|
dst(i,j,k) = rdst(i,k)
|
|
end do
|
|
end do
|
|
|
|
end do
|
|
|
|
end subroutine ddz_ref
|
|
|
|
|
|
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
|
|
|
|
|
|
do k = 1,nzp
|
|
|
|
call tp2(xsrc, src(:,:,k), nyp, nxp)
|
|
|
|
call d2fnonp(nxp, hxp, xsrc, xdst, nyp, 1)
|
|
|
|
call tp2(dst(:,:,k), xdst, nxp, nyp)
|
|
|
|
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
|
|
|
|
|
|
call tp2(zsrc, src, nyp*nzp, nxp)
|
|
|
|
do i = 1,nxp
|
|
|
|
call d2fp(nzp, hzp, zsrc(:,:,i), zdst(:,:,i), nyp, 3)
|
|
|
|
end do
|
|
|
|
call tp2(dst, zdst, nxp, nyp*nzp)
|
|
|
|
end subroutine d2dz
|
|
|
|
|
|
subroutine d2dz_ref(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
|
|
do i = 1,nxp
|
|
rsrc(i,k) = src(i,j,k)
|
|
end do
|
|
end do
|
|
|
|
call d2fp(nzp, hzp, rsrc, rdst, nxp, 3)
|
|
|
|
do k = 1,nzp
|
|
do i = 1,nxp
|
|
dst(i,j,k) = rdst(i,k)
|
|
end do
|
|
end do
|
|
|
|
end do
|
|
|
|
end subroutine d2dz_ref
|
|
|
|
|
|
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,jj+nb-1
|
|
DO i=ii,ii+nb-1
|
|
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
|
|
|
|
end module m_calculate
|