module biopermute

  implicit none
  private

  public  :: bio_permute, bio_depermute, npr_task_perm
  private :: bio_mpi_task_perm, bio_mpi_task_deperm

  real(8), allocatable, save, public :: eco_o(:,:), ben_o(:,:)

contains

!===============================================================================

  subroutine bio_permute ( mi, ni, ki, n2d, n3d, msrf, mcol, kh, ptc, ben, nl, &
                           nu, nben )
  !  wrapper for permuting eco and benthos arrays as needed

    implicit none

    integer(4), intent(in)    :: mi, ni, ki, n2d, n3d, nl, nu,nben
    integer(4), intent(in)    :: msrf(0:,0:), mcol(0:), kh(0:)
    real(8),    intent(inout) :: ptc(0:,1:)
    real(8),    intent(inout) :: ben(0:,1:)

    integer(4) :: neco

    neco = (nu-nl+1)
    call bio_mpi_task_perm (mi,ni,ki,n2d,n3d,neco,nben,msrf,mcol,kh,ptc,ben,nl)

  end subroutine bio_permute

!===============================================================================

  subroutine bio_depermute( mi, ni, ki, n2d, n3d, neco, nben, nl, msrf, mcol,  &
                            kh, ptc, ben)
  !  wrapper for depermuting eco and benthos arrays as needed

    implicit none

    integer(4), intent(in) :: mi, ni, ki, n2d, n3d, neco, nben, nl
    integer(4), intent(in) :: msrf(0:,0:), mcol(0:), kh(0:)
    real(8),    intent(in) :: ptc(0:,1:)
    real(8),    intent(in) :: ben(0:,1:)

    logical :: do_alloc

    do_alloc = .true.
    if (allocated(eco_o)) then
      do_alloc = .false.
      if (n3d > size(eco_o(:,1)) - 1) then
        deallocate(eco_o)
        do_alloc = .true.
      endif
    endif
    if (do_alloc) allocate( eco_o(0:n3d,neco) )
    do_alloc = .true.
    if (allocated(ben_o)) then
      do_alloc = .false.
      if (n2d > size(ben_o(:,1)) - 1) then
        deallocate(ben_o)
        do_alloc = .true.
      endif
    endif
    if (do_alloc) allocate( ben_o(0:n2d,nben) )

    call bio_mpi_task_deperm( mi, ni, ki, neco, nben, nl, msrf, mcol, kh, ptc, &
                              ben)

  end subroutine bio_depermute

!===============================================================================

  subroutine bio_mpi_task_perm(mi,ni,ki,n2d,n3d,neco,nben,msrf,mcol,kh,ptc,ben,&
                               nl)
  !  mpi task permute eco and benthos arrays to full domain

    implicit none

    integer(4), intent(in)    :: mi,ni,ki,n2d,n3d,neco,nben,nl
    integer(4), intent(in)    :: msrf(0:,0:), mcol(0:), kh(0:)
    real(8),    intent(inout) :: ptc(0:,1:)
    real(8),    intent(inout) :: ben(0:,1:)

    integer(4) :: i, j, k, mo, mc, msf
    real(8)    :: eco_n(1:n3d,1:neco), ben_n(1:n2d,1:nben)

    !  mo: index in old permutation
    !  mc: index in new permutation
    k  = 1
    mo = 0
    do j=1,ni
      do i=1,mi
        mc = msrf(i,j)
        if (mc <= 0) cycle
        mo = mo + 1
        eco_n(mc,1:neco) = ptc(mo,nl:nl-1+neco)
        ben_n(mc,1:nben) = ben(mo,1:nben)
      enddo
    enddo
    do k=2,ki
      do j=1,ni
        do i=1,mi
          msf = msrf(i,j)
          if (k > kh(msf)) cycle
          mc = mcol(msf)+k-2
          mo = mo + 1
          eco_n(mc,1:neco) = ptc(mo,nl:nl-1+neco)
        enddo
      enddo
    enddo

    ben(1:n2d,1:nben)       = ben_n(1:n2d,1:nben)
    ptc(1:n3d,nl:nl-1+neco) = eco_n(1:n3d,1:neco)

  end subroutine bio_mpi_task_perm

!===============================================================================

  subroutine npr_task_perm(mi,ni,n2d,msrf,npr)
  !  mpi task permute npr array to full domain

    implicit none

    integer(4), intent(in)    :: mi,ni,n2d
    integer(4), intent(in)    :: msrf(0:,0:)
    real(8),    intent(inout) :: npr(0:)

    integer(4) :: i, j, k, mo, mc, msf
    real(8)    :: npr_n(1:n2d)

    !  mo: index in old permutation
    !  mc: index in new permutation
    mo = 0
    do j=1,ni
      do i=1,mi
        mc = msrf(i,j)
        if (mc <= 0) cycle
        mo = mo + 1
        npr_n(mc) = npr(mo)
      enddo
    enddo

    npr(1:n2d)       = npr_n(1:n2d)

  end subroutine npr_task_perm

!===============================================================================

  subroutine bio_mpi_task_deperm(mi, ni, ki, neco, nben, nl, msrf, mcol, kh,   &
                                 ptc, ben)
  !  mpi task de-permute eco and benthos arrays to full domain

    use exits, only : exitme

    implicit none

    integer(4), intent(in) :: mi, ni, ki, neco, nben, nl
    integer(4), intent(in) :: msrf(0:,0:), mcol(0:), kh(0:)
    real(8),    intent(in) :: ptc(0:,1:)
    real(8),    intent(in) :: ben(0:,1:)

    integer(4) :: i, j, k, mc, mn, msf

    if (.not.allocated(eco_o)) then
      call exitme('Error in bio_mpi_task_deperm: eco_o not allocated')
    endif
    if (.not.allocated(ben_o)) then
      call exitme('Error in bio_mpi_task_deperm: ben_o not allocated')
    endif

    ! now, let's permute data arrays:
    eco_o(0,1:neco) = ptc(0,nl:nl-1+neco)
    mc = 0
    do k=1,ki
      do j=1,ni
        do i=1,mi
          msf = msrf(i,j)
          if (k > kh(msf)) cycle
          if (k == 1) then
            mn = msf
          else
            mn = mcol(msf)+k-2
          endif
          mc = mc + 1
          eco_o(mc,1:neco) = ptc(mn,nl:nl-1+neco)
        enddo
      enddo
    enddo
    ben_o(0,1:nben) = ben(0,1:nben)
    mc = 0
    k  = 1
    do j=1,ni
      do i=1,mi
        mn = msrf(i,j)
        if (mn == 0) cycle
        mc = mc + 1
        ben_o(mc,1:nben) = ben(mn,1:nben)
      enddo
    enddo

  end subroutine bio_mpi_task_deperm

!===============================================================================

end module biopermute
