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 = BLOCKSIZE

    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
        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_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
        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_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,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

end module m_calculate