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(:,:,:) :: 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)

     allocate(xsrc(nyp, nxp), stat=ierr)
     allocate(xdst(nyp, nxp), 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(zsrc)
     deallocate(zdst)

     end subroutine m_calculate_finalize


     subroutine ddx(dst, src)

        real*8, dimension(nxp,nyp,nzp) :: src, 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) :: src, 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) :: src, dst

        integer :: i, j ,k


        call tp2(zsrc, src, nyp, nzp*nxp)

        do k = 1,nzp

        call dfp(nzp, hzp, zsrc, zdst, nyp, 3)

        end do

        call tp2(dst, zdst, nxp, nyp*nzp)

     end subroutine ddz

     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.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