module m_calculate

  use, intrinsic :: iso_fortran_env, only: real64
  use Compact
  use m_parameters

  implicit none

    real(real64), allocatable, dimension(:,:) :: xsrc
    real(real64), allocatable, dimension(:,:) :: xdst

    real(real64), allocatable, dimension(:,:) :: rsrc
    real(real64), 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)
     if (ierr /= 0) then
         write(0,*) "Error: allocation of xsrc failed on process", myid
         call MPI_ABORT(MPI_COMM_TASK, 1, mpi_err)
     end if

     allocate(xdst(nb, nxp), stat=ierr)
     if (ierr /= 0) then
         write(0,*) "Error: allocation of xdst failed on process", myid
         call MPI_ABORT(MPI_COMM_TASK, 1, mpi_err)
     end if

     allocate(rsrc(nxp, nzp), stat=ierr)
     if (ierr /= 0) then
         write(0,*) "Error: allocation of rsrc failed on process", myid
         call MPI_ABORT(MPI_COMM_TASK, 1, mpi_err)
     end if

     allocate(rdst(nxp, nzp), stat=ierr)
     if (ierr /= 0) then
         write(0,*) "Error: allocation of rdst failed on process", myid
         call MPI_ABORT(MPI_COMM_TASK, 1, mpi_err)
     end if

     end subroutine m_calculate_init


     subroutine m_calculate_finalize

     deallocate(xsrc)
     deallocate(xdst)

     deallocate(rsrc)
     deallocate(rdst)

     end subroutine m_calculate_finalize


     subroutine ddx1d(dst, src)

        real*8, dimension(1,nxp), intent(in) :: src
        real*8, dimension(1,nxp), intent(out) :: dst

        integer :: i, j ,k
        integer :: ju

        call dfnonp(nxp, hxp, src, dst, 1, 1)

     end subroutine ddx1d


     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 d2dx1d(dst, src)

        real*8, dimension(nxp), intent(in) :: src
        real*8, dimension(nxp), intent(out) :: dst

        integer :: i, j ,k
        integer :: ju

        call d2fnonp(nxp, hxp, src, dst, 1, 1)

     end subroutine d2dx1d


     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(real64) function rxn_rate (c)

      real(real64) :: c

        if(c.lt.0._real64) c=0._real64
        if(c.gt.1._real64) c=1._real64

        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(real64) function threshold_min_max (c, minc, maxc)

      real(real64) :: c
      real(real64) :: minc, maxc

        if ((c.lt.minc) .or. (c.gt.maxc)) then
            threshold_min_max = 0._real64
        else
            threshold_min_max = 1.0_real64
        end if

     end function threshold_min_max

     real(real64) function positive (c)

      real(real64) :: c

        if (c > 0.0_real64) then
            positive = c
        else
            positive = 0._real64
        end if

     end function positive

     real(real64) function negative (c)

      real(real64) :: c

        if (c < 0.0_real64) then
            negative = c
        else
            negative = 0._real64
        end if

     end function negative

end module m_calculate