deleted m_calculate and added m_chemistry

This commit is contained in:
ignis 2019-08-30 16:26:24 +09:00
parent c50c487b54
commit cf620e1668
4 changed files with 50 additions and 319 deletions

View file

@ -1,279 +0,0 @@
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

43
code/m_chemistry.f90 Normal file
View file

@ -0,0 +1,43 @@
module m_chemistry
use m_parameters
implicit none
contains
real function rate_1step (yr, theta)
real, intent(in) :: yr
real, intent(in) :: theta
real :: y
real :: t_reduce
y=yr
if(yr.lt.0.) y=0.
if(yr.gt.1.) y=1.
t_reduce=theta
if(theta.lt.0.) t_reduce=0.
if(theta.gt.1.) t_reduce=1.
if (t_reduce.gt.c_ref) then
rate_1step = pre*y*exp(-ac/(1.+bc*t_reduce))
else if (t_reduce.le.c_cut) then
rate_1step = min_wr
else
rate_1step = &
((refwr-min_wr)*exp(prof_wr*(t_reduce-c_ref)) + min_wr - refwr*exp(prof_wr*(c_cut-c_ref))) &
/ (1.-exp(prof_wr*(c_cut-c_ref)))
endif
end function rate_1step
end module m_chemistry

View file

@ -12,24 +12,24 @@ endif
compiler = gfortran
ex : test.o ysolve.o m_parameters.o m_calculate.o m_compact.o
${compiler} -o ex test.o ysolve.o m_parameters.o m_calculate.o m_compact.o
ex : test.o ysolve.o m_chemistry.o m_parameters.o m_compact.o
${compiler} -o ex test.o ysolve.o m_chemistry.o m_parameters.o m_compact.o
test.o : test.f90 ysolve.mod m_compact.mod
${compiler} -c ${flags} test.f90
ysolve.o ysolve.mod : ysolve.f90 m_compact.mod m_parameters.mod
ysolve.o ysolve.mod : ysolve.f90 m_compact.mod m_chemistry.mod m_parameters.mod
${compiler} -c ${flags} ysolve.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
m_chemistry.o m_chemistry.mod : m_chemistry.f90 m_parameters.mod
${compiler} -c ${flags} m_chemistry.f90
test: ex
sh test.sh

View file

@ -1,6 +1,7 @@
MODULE ysolve
USE m_compact
USE m_parameters
USE m_chemistry
IMPLICIT NONE
@ -412,38 +413,4 @@
!------------------------------------------------------------------------
real function rate_1step (yr, theta)
real, intent(in) :: yr
real, intent(in) :: theta
real :: y
real :: t_reduce
y=yr
if(yr.lt.0.) y=0.
if(yr.gt.1.) y=1.
t_reduce=theta
if(theta.lt.0.) t_reduce=0.
if(theta.gt.1.) t_reduce=1.
if (t_reduce.gt.c_ref) then
rate_1step = pre*y*exp(-ac/(1.+bc*t_reduce))
else if (t_reduce.le.c_cut) then
rate_1step = min_wr
else
rate_1step = &
((refwr-min_wr)*exp(prof_wr*(t_reduce-c_ref)) + min_wr - refwr*exp(prof_wr*(c_cut-c_ref))) &
/ (1.-exp(prof_wr*(c_cut-c_ref)))
endif
end function rate_1step
END MODULE ysolve