Non-local boundary conditions

slab.f90

! Apply various non-local boundary conditions to a 3d grid
module slab
  implicit none

  integer, parameter :: nghosts = 2

contains

  ! Given an input array that is defined in the interior (on the dots
  ! and stars), fill the boundaries by rotating the array by 90
  ! degrees in either direction. The pivot point is the top left dot.
  ! The array values marked with a star end up being used to fill in
  ! the boundaries, the array values marked with a dot end up unused.
  !
  ! Formula: f(-x,y,z) = f(y,x,z)
  !          f(x,-y,z) = f(y,x,z)
  !
  ! Array shape:
  !   1 2 2 2 2
  !   3 . . . .
  !   3 . * * *
  !   3 . * * *
  !   3 . * * *
  subroutine rotating90(dst, src)
    double precision, intent(inout) :: dst(:,:,:)
    double precision, intent(in) :: src(:,:,:)
    integer :: i, j, k
    
    if (size(dst) /= size(src)) stop
    if (size(dst,1) /= size(dst,2)) stop

    i0 = nghosts
    j0 = nghosts

    ! Fill region marked "1" by rotating array by 180 degrees:
    forall (k = 1 : size(a,3))
       forall (j = 1 : nghosts)
          forall (i = 1 : nghosts)
             dst(i0-i, j0-j, k) = src(i0+i, j0+j, k)
          end forall
       end forall
    end forall

    ! Fill region marked "2" by rotating array 90 degrees counterclockwise:
    forall (k = 1 : size(a,3))
       forall (j = 1 : nghosts)
          forall (i = 1 : size(a,1)-nghosts-1)
             dst(i0+i, j0-j, k) = src(i0+j, j0+i, k)
          end forall
       end forall
    end forall

    ! Fill region marked "3" by rotating array 90 degrees clockwise:
    forall (k = 1 : size(a,3))
       forall (j = 1 : size(a,2)-nghosts-1)
          forall (i = 1 : nghosts)
             dst(i0-i, j0+j, k) = src(i0+j, j0+i, k)
          end forall
       end forall
    end forall
  end subroutine rotating90

  ! Given an input array that is defined in the interior (on the dots
  ! and stars), fill the boundaries by rotating the array by 180
  ! degrees in either direction. The pivot point is the centre dot.
  ! The array values marked with a star end up being used to fill in
  ! the boundaries, the array values marked with a dot end up unused.
  !
  ! Formula: f(-x,y,z) = f(x,-y,z)
  !
  ! Array shape:
  !   1 . * * *
  !   1 . * * *
  !   1 . * * *
  !   1 . * * *
  !   1 . * * *
  subroutine rotating180(dst, src)
    double precision, intent(inout) :: dst(:,:,:)
    double precision, intent(in) :: src(:,:,:)
    integer :: i, j, k
    
    if (size(dst) /= size(src)) stop
    if (mod(size(dst,2), 2) /= 1) stop

    i0 = nghosts
    j0 = nghosts

    ! Fill region marked "1" by rotating array by 180 degrees:
    forall (k = 1 : size(a,3))
       forall (j = 1 : size(a,2))
          forall (i = 1 : nghosts)
             dst(i0-i, j, k) = src(i0+i, size(a,2)+1-j, k)
          end forall
       end forall
    end forall
  end subroutine rotating180

  ! Given an input array that is defined in the interior (on the
  ! stars), fill the boundaries by cycling the array by 1/2 its length
  ! in the x direction and wrapping around in the y direction.
  !
  ! Formula: f(-theta,phi,z) = f(theta,phi+pi,z)
  !
  ! Array shape:
  !   1 1 1 2 2 2
  !   . . . . . .
  !   * * * * * *
  !   * * * * * *
  !   * * * * * *
  !   . . . . . .
  !   3 3 3 4 4 4
  subroutine polar
    double precision, intent(inout) :: dst(:,:,:)
    double precision, intent(in) :: src(:,:,:)
    integer :: i, j, k
    
    if (size(dst) /= size(src)) stop
    if (mod(size(dst,1), 2) /= 0) stop

    i0 = nghosts
    j0 = nghosts

    ! Fill regions marked "3" and "4":
    forall (k = 1 : size(a,3))
       forall (j = 1 : nghosts)
          forall (i = 1 : size(a,1))
             ii = mod(i + size(a,1) - 1, size(a,1)) + 1
             jj = j + size(a,2) - 2*nghosts
             dst(i, j0, k) = src(1 + mod(ii, jj, k)
          end forall
       end forall
    end forall

    ! ditto for regions 3 and 4

  end subroutine polar

end module slab

Should Line 47 forall (i = 1 : size(a,1)-nghosts-1) read forall (i = 1 : size(a,1)-nghosts)?

Yes, probably. Please treat this as pseudo-code to illustrate what the boundary conditions do. I didn't test the code thoroughly, and there are likely errors.

For rotating90, I assume you do want the ghost cells at the higher end of the domain being filled with rotated data too. Is that correct?

The boundary conditions should be applied at the lower i and j boundaries only, i.e. on 2 of the 4 boundaries.

The 180 degree rotating boundaries should be applied at the lower i boundary only.