diff --git a/code/m_calculate.f90 b/code/m_calculate.f90 deleted file mode 100644 index d5a2a38..0000000 --- a/code/m_calculate.f90 +++ /dev/null @@ -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 diff --git a/code/m_chemistry.f90 b/code/m_chemistry.f90 new file mode 100644 index 0000000..28d11a2 --- /dev/null +++ b/code/m_chemistry.f90 @@ -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 diff --git a/code/makefile b/code/makefile index 4ad6d7b..11d7845 100644 --- a/code/makefile +++ b/code/makefile @@ -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 diff --git a/code/ysolve.f90 b/code/ysolve.f90 index 6c6743e..82be8ef 100644 --- a/code/ysolve.f90 +++ b/code/ysolve.f90 @@ -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