BoundaryCondition Subroutine

  subroutine BoundaryCondition
    implicit none
    integer :: i,j,k

!$omp target teams distribute parallel do collapse(3)
    do k=1,kn-1
    do j=1,jn-1
    do i=1,mgn
      BsXend(i,j,k,1) =  d(ie-mgn+i,j,k)
      BsXend(i,j,k,2) = ei(ie-mgn+i,j,k)
      BsXend(i,j,k,nv1) = v1(ie-mgn+i,j,k)
      BsXend(i,j,k,nv2) = v2(ie-mgn+i,j,k)
      BsXend(i,j,k,nv3) = v3(ie-mgn+i,j,k)
      BsXend(i,j,k,6) = b1(ie-mgn+i,j,k)
      BsXend(i,j,k,7) = b2(ie-mgn+i,j,k)
      BsXend(i,j,k,8) = b3(ie-mgn+i,j,k)
      BsXend(i,j,k,9) = bp(ie-mgn+i,j,k)
      BsXend(i,j,k,10:nbc) = Xcomp(1:ncomp,ie-mgn+i,j,k)

      BsXstt(i,j,k,1) =  d(is+i-1,j,k)
      BsXstt(i,j,k,2) = ei(is+i-1,j,k)
      BsXstt(i,j,k,nv1) = v1(is+i-1,j,k)
      BsXstt(i,j,k,nv2) = v2(is+i-1,j,k)
      BsXstt(i,j,k,nv3) = v3(is+i-1,j,k)
      BsXstt(i,j,k,6) = b1(is+i-1,j,k)
      BsXstt(i,j,k,7) = b2(is+i-1,j,k)
      BsXstt(i,j,k,8) = b3(is+i-1,j,k)
      BsXstt(i,j,k,9) = bp(is+i-1,j,k)
      BsXstt(i,j,k,10:nbc) = Xcomp(1:ncomp,is+i-1,j,k)
    end do
    end do
    end do

!$omp target teams distribute parallel do collapse(3)
    do k=1,kn-1
    do i=1,in-1
    do j=1,mgn
      BsYend(i,j,k,1) =  d(i,je-mgn+j,k)
      BsYend(i,j,k,2) = ei(i,je-mgn+j,k)
      BsYend(i,j,k,nv1) = v1(i,je-mgn+j,k)
      BsYend(i,j,k,nv2) = v2(i,je-mgn+j,k)
      BsYend(i,j,k,nv3) = v3(i,je-mgn+j,k)
      BsYend(i,j,k,6) = b1(i,je-mgn+j,k)
      BsYend(i,j,k,7) = b2(i,je-mgn+j,k)
      BsYend(i,j,k,8) = b3(i,je-mgn+j,k)
      BsYend(i,j,k,9) = bp(i,je-mgn+j,k)
      BsYend(i,j,k,10:nbc) = Xcomp(1:ncomp,i,je-mgn+j,k)

      BsYstt(i,j,k,1) =  d(i,js+j-1,k)
      BsYstt(i,j,k,2) = ei(i,js+j-1,k)
      BsYstt(i,j,k,nv1) = v1(i,js+j-1,k)
      BsYstt(i,j,k,nv2) = v2(i,js+j-1,k)
      BsYstt(i,j,k,nv3) = v3(i,js+j-1,k)
      BsYstt(i,j,k,6) = b1(i,js+j-1,k)
      BsYstt(i,j,k,7) = b2(i,js+j-1,k)
      BsYstt(i,j,k,8) = b3(i,js+j-1,k)
      BsYstt(i,j,k,9) = bp(i,js+j-1,k)
      BsYstt(i,j,k,10:nbc) = Xcomp(1:ncomp,i,js+j-1,k)
    end do
    end do
    end do

!$omp target teams distribute parallel do collapse(3)
    do j=1,jn-1
    do i=1,in-1
    do k=1,mgn
      BsZend(i,j,k,1) =  d(i,j,ke-mgn+k)
      BsZend(i,j,k,2) = ei(i,j,ke-mgn+k)
      BsZend(i,j,k,nv1) = v1(i,j,ke-mgn+k)
      BsZend(i,j,k,nv2) = v2(i,j,ke-mgn+k)
      BsZend(i,j,k,nv3) = v3(i,j,ke-mgn+k)
      BsZend(i,j,k,6) = b1(i,j,ke-mgn+k)
      BsZend(i,j,k,7) = b2(i,j,ke-mgn+k)
      BsZend(i,j,k,8) = b3(i,j,ke-mgn+k)
      BsZend(i,j,k,9) = bp(i,j,ke-mgn+k)
      BsZend(i,j,k,10:nbc) = Xcomp(1:ncomp,i,j,ke-mgn+k)

      BsZstt(i,j,k,1) =  d(i,j,ks+k-1)
      BsZstt(i,j,k,2) = ei(i,j,ks+k-1)
      BsZstt(i,j,k,nv1) = v1(i,j,ks+k-1)
      BsZstt(i,j,k,nv2) = v2(i,j,ks+k-1)
      BsZstt(i,j,k,nv3) = v3(i,j,ks+k-1)
      BsZstt(i,j,k,6) = b1(i,j,ks+k-1)
      BsZstt(i,j,k,7) = b2(i,j,ks+k-1)
      BsZstt(i,j,k,8) = b3(i,j,ks+k-1)
      BsZstt(i,j,k,9) = bp(i,j,ks+k-1)
      BsZstt(i,j,k,10:nbc) = Xcomp(1:ncomp,i,j,ks+k-1)
    end do
    end do
    end do

    ! ----------------------
    ! Exchange / apply physical BCs (MPI on host, BC possibly on device)
    ! ----------------------
    call XbcSendRecv()
    call YbcSendRecv()
    call ZbcSendRecv()

    ! ----------------------
    ! Unpack from receive buffers on device
    ! ----------------------

!$omp target teams distribute parallel do collapse(3)
    do k=1,kn-1
    do j=1,jn-1
    do i=1,mgn
       d(i,j,k) = BrXstt(i,j,k,1)
      ei(i,j,k) = BrXstt(i,j,k,2)
      v1(i,j,k) = BrXstt(i,j,k,nv1)
      v2(i,j,k) = BrXstt(i,j,k,nv2)
      v3(i,j,k) = BrXstt(i,j,k,nv3)
      b1(i,j,k) = BrXstt(i,j,k,6)
      b2(i,j,k) = BrXstt(i,j,k,7)
      b3(i,j,k) = BrXstt(i,j,k,8)
      bp(i,j,k) = BrXstt(i,j,k,9)
      Xcomp(1:ncomp,i,j,k) = BrXstt(i,j,k,10:nbc)

       d(ie+i,j,k) = BrXend(i,j,k,1)
      ei(ie+i,j,k) = BrXend(i,j,k,2)
      v1(ie+i,j,k) = BrXend(i,j,k,nv1)
      v2(ie+i,j,k) = BrXend(i,j,k,nv2)
      v3(ie+i,j,k) = BrXend(i,j,k,nv3)
      b1(ie+i,j,k) = BrXend(i,j,k,6)
      b2(ie+i,j,k) = BrXend(i,j,k,7)
      b3(ie+i,j,k) = BrXend(i,j,k,8)
      bp(ie+i,j,k) = BrXend(i,j,k,9)
      Xcomp(1:ncomp,ie+i,j,k) = BrXend(i,j,k,10:nbc)
    end do
    end do
    end do

!$omp target teams distribute parallel do collapse(3)
    do k=1,kn-1
    do j=1,mgn
    do i=1,in-1
       d(i,j,k) = BrYstt(i,j,k,1)
      ei(i,j,k) = BrYstt(i,j,k,2)
      v1(i,j,k) = BrYstt(i,j,k,nv1)
      v2(i,j,k) = BrYstt(i,j,k,nv2)
      v3(i,j,k) = BrYstt(i,j,k,nv3)
      b1(i,j,k) = BrYstt(i,j,k,6)
      b2(i,j,k) = BrYstt(i,j,k,7)
      b3(i,j,k) = BrYstt(i,j,k,8)
      bp(i,j,k) = BrYstt(i,j,k,9)
      Xcomp(1:ncomp,i,j,k) = BrYstt(i,j,k,10:nbc)

       d(i,je+j,k) = BrYend(i,j,k,1)
      ei(i,je+j,k) = BrYend(i,j,k,2)
      v1(i,je+j,k) = BrYend(i,j,k,nv1)
      v2(i,je+j,k) = BrYend(i,j,k,nv2)
      v3(i,je+j,k) = BrYend(i,j,k,nv3)
      b1(i,je+j,k) = BrYend(i,j,k,6)
      b2(i,je+j,k) = BrYend(i,j,k,7)
      b3(i,je+j,k) = BrYend(i,j,k,8)
      bp(i,je+j,k) = BrYend(i,j,k,9)
      Xcomp(1:ncomp,i,je+j,k) = BrYend(i,j,k,10:nbc)
    end do
    end do
    end do

!$omp target teams distribute parallel do collapse(3)
    do k=1,mgn
    do j=1,jn-1
    do i=1,in-1
       d(i,j,k) = BrZstt(i,j,k,1)
      ei(i,j,k) = BrZstt(i,j,k,2)
      v1(i,j,k) = BrZstt(i,j,k,nv1)
      v2(i,j,k) = BrZstt(i,j,k,nv2)
      v3(i,j,k) = BrZstt(i,j,k,nv3)
      b1(i,j,k) = BrZstt(i,j,k,6)
      b2(i,j,k) = BrZstt(i,j,k,7)
      b3(i,j,k) = BrZstt(i,j,k,8)
      bp(i,j,k) = BrZstt(i,j,k,9)
      Xcomp(1:ncomp,i,j,k) = BrZstt(i,j,k,10:nbc)

       d(i,j,ke+k) = BrZend(i,j,k,1)
      ei(i,j,ke+k) = BrZend(i,j,k,2)
      v1(i,j,ke+k) = BrZend(i,j,k,nv1)
      v2(i,j,ke+k) = BrZend(i,j,k,nv2)
      v3(i,j,ke+k) = BrZend(i,j,k,nv3)
      b1(i,j,ke+k) = BrZend(i,j,k,6)
      b2(i,j,ke+k) = BrZend(i,j,k,7)
      b3(i,j,ke+k) = BrZend(i,j,k,8)
      bp(i,j,ke+k) = BrZend(i,j,k,9)
      Xcomp(1:ncomp,i,j,ke+k) = BrZend(i,j,k,10:nbc)
    end do
    end do
    end do

  end subroutine BoundaryCondition