module boundary_module
  use common_module
contains
  function bc(U0,dU,clock)
    implicit none
    real(DP),dimension(xlen,ylen,zlen,8) :: bc
    real(DP),intent(IN),dimension(xlen,ylen,zlen,8)::U0,dU
    real(DP),intent(IN) :: clock
    real(DP) :: V0(xlen,ylen,zlen),v8(8)
    integer :: cnt,s(3)
    bc(:,:,:,Sax)=U0(:,:,:,Sax)
    bc(:,:,:,Vax)=U0(:,:,:,Vax)
    bc(:,:,:,Bax)=U0(:,:,:,Bax)
    if (rkstep .lt. 0) then
       ! Energy
       v0=U0(:,:,:,Sax(2))+dU(:,:,:,Sax(2))
       call periodic_1(1,1,v0)    !xL-boundary
       call periodic_1(1,2,v0)    !xR-boundary
       call periodic_2(2,1,1,v0)  !yL-boundary
       call periodic_2(2,1,2,v0)  !yR-boundary
       call periodic_3(3,1,v0)    !zL-boundary
       call periodic_3(3,2,v0)    !zR-boundary
       bc(:,:,:,Sax(2))=v0
    end if
  contains
    subroutine loopcounter(count,sizes,org,subs)
      implicit none
      integer,intent(IN) :: count,sizes(:),org(:)
      integer,intent(OUT) :: subs(:)
      integer :: dim,s,quo(size(subs))
      dim=size(sizes)
      quo=count
      do s=2,dim; quo(s)=(quo(s-1)-1)/sizes(s-1)+1
      end do
      do s=1,dim; subs(s)=mod(quo(s)-1,sizes(s))+org(s)
      end do
    end subroutine loopcounter
    subroutine common_in_bc(axis,ends,ar,amar,abox,endc,kdlt,kdlt2,kdlt3)
      implicit none
      integer,intent(IN) :: axis,ends
      integer,intent(OUT) :: ar(2),amar(2),abox,endc(2)
      integer,intent(OUT) :: kdlt(3),kdlt2(3),kdlt3(3)
      kdlt=0**abs(axis-(/1,2,3/)); kdlt2=0**abs(mod(axis,3)+1-(/1,2,3/))
      kdlt3=1-kdlt-kdlt2
      ar=xr*kdlt(1)+yr*kdlt(2)+zr*kdlt(3)
      amar=xmar*kdlt(1)+ymar*kdlt(2)+zmar*kdlt(3)
      abox=xbox*kdlt(1)+ybox*kdlt(2)+zbox*kdlt(3)
      endc(2)=0**abs(ends-2) ; endc(1)=1-endc(2)
    end subroutine common_in_bc
    subroutine periodic_1(axis,ends,V0)
      implicit none
      real(DP),intent(INOUT),dimension(xlen,ylen,zlen)::V0
      integer,intent(IN) :: axis,ends
      integer :: endc(2),ir(2),jr(2),kr(2),cnt,szi,szj,szk
      integer :: kdlt(3),kdlt2(3),kdlt3(3),kji(3),kji1(3),ar(2),amar(2),abox
      call common_in_bc(axis,ends,ar,amar,abox,endc,kdlt,kdlt2,kdlt3)
      ! memo ----------------------------------------------------------!
      !endc=(/1,0/)...left boundary  i=1,...,xmar(1)                   !
      !   -> U(i,yr(1):yr(2),zr(1):zr(2),:)=U(xr(2)-xmar(1)+i,...,:)   !
      !endc=(/0,1/)...right boundary  i=1,...,xmar(2)                  !
      !   -> U(xr(2)+i,yr(1):yr(2),zr(1):zr(2),:)=U(xr(1)-1+i,...,:)   !
      !----------------------------------------------------------------!
      ir=((/1,xmar(ends)/)+xr(2)*endc(2))*kdlt(1)+xr*(1-kdlt(1))
      jr=((/1,ymar(ends)/)+yr(2)*endc(2))*kdlt(2)+yr*(1-kdlt(2))
      kr=((/1,zmar(ends)/)+zr(2)*endc(2))*kdlt(3)+zr*(1-kdlt(3))
      szi=ir(2)-ir(1)+1; szj=jr(2)-jr(1)+1; szk=kr(2)-kr(1)+1
      !$omp parallel do private(kji,kji1),default(shared)
      do cnt=1,szi*szj*szk
         call loopcounter(cnt,(/szk,szj,szi/),(/kr(1),jr(1),ir(1)/),kji)
         kji1=kji
         kji(4-axis)=(ar(2)-mod(amar(1)-kji1(4-axis),abox))*endc(1)&
              +(ar(1)+mod(kji1(4-axis)-ar(2)-1,abox))*endc(2)
         V0(kji1(3),kji1(2),kji1(1))=V0(kji(3),kji(2),kji(1))
      end do
      !$omp end parallel do
    end subroutine periodic_1

    subroutine periodic_2(axis1,axis2,ends,V0)
      implicit none
      real(DP),intent(INOUT),dimension(xlen,ylen,zlen)::V0
      integer,intent(IN) :: axis1,axis2,ends
      integer :: endc(2),ir(2),jr(2),kr(2),cnt,szi,szj,szk
      integer :: kdlt(3),kdlt2(3),kdlt3(3),kji(3),kji1(3),ar(2),amar(2),abox
      call common_in_bc(axis1,ends,ar,amar,abox,endc,kdlt,kdlt2,kdlt3)
      kdlt2=0**abs(axis2-(/1,2,3/)); kdlt3=(/1,1,1/)-kdlt-kdlt2
      ir=((/1,xmar(ends)/)+xr(2)*endc(2))*kdlt(1)&
           +(/1,xlen/)*kdlt2(1)+xr*kdlt3(1)
      jr=((/1,ymar(ends)/)+yr(2)*endc(2))*kdlt(2)&
           +(/1,ylen/)*kdlt2(2)+yr*kdlt3(2)
      kr=((/1,zmar(ends)/)+zr(2)*endc(2))*kdlt(3)&
           +(/1,zlen/)*kdlt2(3)+zr*kdlt3(3)
      szi=ir(2)-ir(1)+1; szj=jr(2)-jr(1)+1; szk=kr(2)-kr(1)+1
      ! memo ----------------------------------------------------------!
      !endc=(/1,0/)...left boundary  i=1,...,xmar(1)                   !
      !   -> U(i,yr(1):yr(2),zr(1):zr(2),:)=U(xr(2)-xmar(1)+i,...,:)   !
      !endc=(/0,1/)...right boundary  i=1,...,xmar(2)                  !
      !   -> U(xr(2)+i,yr(1):yr(2),zr(1):zr(2),:)=U(xr(1)-1+i,...,:)   !
      !----------------------------------------------------------------!
      !$omp parallel do private(kji,kji1),default(shared)
      do cnt=1,szi*szj*szk
         call loopcounter(cnt,(/szk,szj,szi/),(/kr(1),jr(1),ir(1)/),kji)
         kji1=kji
         kji(4-axis1)=(ar(2)-mod(amar(1)-kji1(4-axis1),abox))*endc(1)&
              +(ar(1)+mod(kji1(4-axis1)-ar(2)-1,abox))*endc(2)
         V0(kji1(3),kji1(2),kji1(1))=V0(kji(3),kji(2),kji(1))
      end do
      !$omp end parallel do
    end subroutine periodic_2

    subroutine periodic_3(axis,ends,V0)
      implicit none
      real(DP),intent(INOUT),dimension(xlen,ylen,zlen)::V0
      integer,intent(IN) :: axis,ends
      integer :: endc(2),ir(2),jr(2),kr(2),cnt,szi,szj,szk
      integer :: kdlt(3),kdlt2(3),kdlt3(3),kji(3),kji1(3),ar(2),amar(2),abox
      call common_in_bc(axis,ends,ar,amar,abox,endc,kdlt,kdlt2,kdlt3)
      ir=((/1,xmar(ends)/)+xr(2)*endc(2))*kdlt(1)+(/1,xlen/)*(1-kdlt(1))
      jr=((/1,ymar(ends)/)+yr(2)*endc(2))*kdlt(2)+(/1,ylen/)*(1-kdlt(2))
      kr=((/1,zmar(ends)/)+zr(2)*endc(2))*kdlt(3)+(/1,zlen/)*(1-kdlt(3))
      szi=ir(2)-ir(1)+1; szj=jr(2)-jr(1)+1; szk=kr(2)-kr(1)+1
      ! memo ----------------------------------------------------------!
      !endc=(/1,0/)...left boundary  i=1,...,xmar(1)                   !
      !   -> U(i,yr(1):yr(2),zr(1):zr(2),:)=U(xr(2)-xmar(1)+i,...,:)   !
      !endc=(/0,1/)...right boundary  i=1,...,xmar(2)                  !
      !   -> U(xr(2)+i,yr(1):yr(2),zr(1):zr(2),:)=U(xr(1)-1+i,...,:)   !
      !----------------------------------------------------------------!
      !$omp parallel do private(kji,kji1),default(shared)
      do cnt=1,szi*szj*szk
         call loopcounter(cnt,(/szk,szj,szi/),(/kr(1),jr(1),ir(1)/),kji)
         kji1=kji
         kji(4-axis)=(ar(2)-mod(amar(1)-kji1(4-axis),abox))*endc(1)&
              +(ar(1)+mod(kji1(4-axis)-ar(2)-1,abox))*endc(2)
         V0(kji1(3),kji1(2),kji1(1))=V0(kji(3),kji(2),kji(1))
      end do
      !$omp end parallel do
    end subroutine periodic_3
  end function bc
end module boundary_module