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) :: eth(xlen,ylen,zlen),V0(xlen,ylen,zlen),v8(8)
    integer :: cnt,s(3),cyc


    ! Mass density
    v0=dU(:,:,:,Sax(1))
    call dirichlet_free_1(1,1,v0)   !xL-boundary
    call dirichlet_free_1(1,2,v0)   !xR-boundary
    if (ylen.ge.2) then
       call dirichlet_free_2(2,1,1,v0)   !yL-boundary
       call dirichlet_free_2(2,1,2,v0)   !yR-boundary
    end if
    if (zlen.ge.2) then
       call dirichlet_free_3(3,1,v0)   !zL-boundary
       call dirichlet_free_3(3,2,v0)   !zR-boundary
    end if
    bc(:,:,:,Sax(1))=U0(:,:,:,Sax(1))+v0

    ! Momentum-x
    v0=dU(:,:,:,Vax(1))
    call dirichlet_free_1(1,1,v0)   !xL-boundary
    call dirichlet_free_1(1,2,v0)   !xR-boundary
    if (ylen.ge.2) then
       call dirichlet_free_2(2,1,1,v0)   !yL-boundary
       call dirichlet_free_2(2,1,2,v0)   !yR-boundary
    end if
    if (zlen.ge.2) then
       call dirichlet_free_3(3,1,v0)   !zL-boundary
       call dirichlet_free_3(3,2,v0)   !zR-boundary
    end if
    bc(:,:,:,Vax(1))=U0(:,:,:,Vax(1))+v0

    ! Momentum-y
    v0=dU(:,:,:,Vax(2))
    call dirichlet_free_1(1,1,v0)   !xL-boundary
    call dirichlet_free_1(1,2,v0)   !xR-boundary
    if (ylen.ge.2) then
       call dirichlet_free_2(2,1,1,v0)   !yL-boundary
       call dirichlet_free_2(2,1,2,v0)   !yR-boundary
    end if
    if (zlen.ge.2) then
       call dirichlet_free_3(3,1,v0)   !zL-boundary
       call dirichlet_free_3(3,2,v0)   !zR-boundary
    end if
    bc(:,:,:,Vax(2))=U0(:,:,:,Vax(2))+v0


    ! Momentum-z
    v0=dU(:,:,:,Vax(3))
    call dirichlet_free_1(1,1,v0)   !xL-boundary
    call dirichlet_free_1(1,2,v0)   !xR-boundary
    if (ylen.ge.2) then
       call dirichlet_free_2(2,1,1,v0)   !yL-boundary
       call dirichlet_free_2(2,1,2,v0)   !yR-boundary
    end if
    if (zlen.ge.2) then
       call dirichlet_free_3(3,1,v0)   !zL-boundary
       call dirichlet_free_3(3,2,v0)   !zR-boundary
    end if
    bc(:,:,:,Vax(3))=U0(:,:,:,Vax(3))+v0


    ! Magnetic Field-x
    v0=dU(:,:,:,Bax(1))
    call dirichlet_free_1(1,1,v0)   !xL-boundary
    call dirichlet_free_1(1,2,v0)   !xR-boundary
    if (ylen.ge.2) then
       call dirichlet_free_2(2,1,1,v0)   !yL-boundary
       call dirichlet_free_2(2,1,2,v0)   !yR-boundary
    end if
    if (zlen.ge.2) then
       call dirichlet_free_3(3,1,v0)   !zL-boundary
       call dirichlet_free_3(3,2,v0)   !zR-boundary
    end if
    bc(:,:,:,Bax(1))=U0(:,:,:,Bax(1))+v0

    ! Magnetic Field-y
    v0=dU(:,:,:,Bax(2))
    call dirichlet_free_1(1,1,v0)   !xL-boundary
    call dirichlet_free_1(1,2,v0)   !xR-boundary
    if (ylen.ge.2) then
       call dirichlet_free_2(2,1,1,v0)   !yL-boundary
       call dirichlet_free_2(2,1,2,v0)   !yR-boundary
    end if
    if (zlen.ge.2) then
       call dirichlet_free_3(3,1,v0)   !zL-boundary
       call dirichlet_free_3(3,2,v0)   !zR-boundary
    end if
    bc(:,:,:,Bax(2))=U0(:,:,:,Bax(2))+v0


    ! Magnetic Field-z
    v0=dU(:,:,:,Bax(3))
    call dirichlet_free_1(1,1,v0)   !xL-boundary
    call dirichlet_free_1(1,2,v0)   !xR-boundary
    if (ylen.ge.2) then
       call dirichlet_free_2(2,1,1,v0)   !yL-boundary
       call dirichlet_free_2(2,1,2,v0)   !yR-boundary
    end if
    if (zlen.ge.2) then
       call dirichlet_free_3(3,1,v0)   !zL-boundary
       call dirichlet_free_3(3,2,v0)   !zR-boundary
    end if
    bc(:,:,:,Bax(3))=U0(:,:,:,Bax(3))+v0

    ! Energy
    v0=dU(:,:,:,Sax(2))
    call dirichlet_free_1(1,1,v0)   !xL-boundary
    call dirichlet_free_1(1,2,v0)   !xR-boundary
    if (ylen.ge.2) then
       call dirichlet_free_2(2,1,1,v0)   !yL-boundary
       call dirichlet_free_2(2,1,2,v0)   !yR-boundary
    end if
    if (zlen.ge.2) then
       call dirichlet_free_3(3,1,v0)   !zL-boundary
       call dirichlet_free_3(3,2,v0)   !zR-boundary
    end if
    bc(:,:,:,Sax(2))=U0(:,:,:,Sax(2))+v0
  contains
    subroutine dirichlet_free_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                                  !
      !   -> U(i,yr(1):yr(2),zr(1):zr(2),:)=U(xr(1),...,:)           !
      !            where i=1,...,xmar(1)                             !
      !endc=(/0,1/)...right boundary  i=1,...,xmar(2)                !
      !   -> U(xr(2)+i,yr(1):yr(2),zr(1):zr(2),:)=U(xr(2),...,:)     !
      !--------------------------------------------------------------!    
      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(ends)
         V0(kji1(3),kji1(2),kji1(1))=V0(kji(3),kji(2),kji(1))
      end do
      !$omp end parallel do
    end subroutine dirichlet_free_1
    subroutine dirichlet_free_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
      !$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(ends)
         V0(kji1(3),kji1(2),kji1(1))=V0(kji(3),kji(2),kji(1))
      end do
      !$omp end parallel do
    end subroutine dirichlet_free_2
    subroutine dirichlet_free_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
      !$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(ends)
         V0(kji1(3),kji1(2),kji1(1))=V0(kji(3),kji(2),kji(1))
      end do
      !$omp end parallel do
    end subroutine dirichlet_free_3
    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
  end function bc
end module boundary_module