1 files changed, 154 insertions, 0 deletions
diff --git a/src/qlm_boundary.F90 b/src/qlm_boundary.F90
new file mode 100644
index 0000000..33196bc
--- /dev/null
+++ b/src/qlm_boundary.F90
@@ -0,0 +1,154 @@
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+
+
+module qlm_boundary
+  use cctk
+  implicit none
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+  private
+  public set_boundary
+  
+  interface set_boundary
+     module procedure set_boundary_real
+     module procedure set_boundary_complex
+  end interface
+  
+contains
+  
+  subroutine set_boundary_real (CCTK_ARGUMENTS, hn, f, parity)
+    DECLARE_CCTK_ARGUMENTS
+    integer, intent(in) :: hn
+    CCTK_REAL           :: f(:,:)
+    integer, intent(in) :: parity
+    
+    CCTK_REAL, allocatable :: ff(:,:)
+    
+    integer :: ni, nj, gi, gj
+    
+    integer :: i, j, jj
+    
+    if (abs(parity) /= 1) then
+       call CCTK_WARN (0, "Parity must be +1 or -1")
+    end if
+    
+    
+    
+    ni = qlm_ntheta(hn)
+    nj = qlm_nphi(hn)
+    gi = qlm_nghoststheta(hn)
+    gj = qlm_nghostsphi(hn)
+    
+    
+    
+    f(       :gi, :nj) = -42
+    f(ni-gi+1:ni, :nj) = -42
+    f(:ni,        :gj) = -42
+    f(:ni, nj-gi+1:nj) = -42
+    
+    allocate (ff(gi, nj))
+    
+    ! lower theta
+    ff(:, :nj) = f(1+gi:2*gi, :nj)
+    
+    ! polar boundary condition at north pole
+    do j=1,nj/2
+       do i=1,gi
+          jj = j + (nj-2*gj) / 2
+          f(i,j) = parity * ff(gi-i+1,jj)
+       end do
+    end do
+    do j=nj/2+1,nj
+       do i=1,gi
+          jj = j - (nj-2*gj) / 2
+          f(i,j) = parity * ff(gi-i+1,jj)
+       end do
+    end do
+    
+    ! upper theta
+    ff(:, :nj) = f(ni-2*gi+1:ni-gi, :nj)
+    
+    ! polar boundary condition at south pole
+    do j=1,nj/2
+       do i=1,gi
+          jj = j + (nj-2*gj) / 2
+          f(ni-i+1,j) = parity * ff(i,jj)
+       end do
+    end do
+    do j=nj/2+1,nj
+       do i=1,gi
+          jj = j - (nj-2*gj) / 2
+          f(ni-i+1,j) = parity * ff(i,jj)
+       end do
+    end do
+    
+    deallocate (ff)
+    
+    
+    
+    allocate (ff(ni,gj))
+    
+    ! lower phi
+    ff(:ni, :) = f(:ni, nj-2*gj+1:nj-gj)
+    
+    ! periodic boundary at null meridian
+    do j=1,gj
+       do i=1,ni
+          f(i,j) = ff(i,j)
+       end do
+    end do
+    
+    ! upper phi
+    ff(:ni, :) = f(:ni, gj+1:2*gj)
+    
+    ! periodic boundary at null meridian
+    do j=1,gj
+       do i=1,ni
+          f(i,nj-j+1) = ff(i,gj-j+1)
+       end do
+    end do
+    
+    deallocate (ff)
+    
+    
+    
+    f(ni+1:, :nj) = 0
+    f(:, nj+1:) = 0
+    
+  end subroutine set_boundary_real
+  
+  
+  
+  subroutine set_boundary_complex (CCTK_ARGUMENTS, hn, f, parity)
+    DECLARE_CCTK_ARGUMENTS
+    integer, intent(in) :: hn
+    CCTK_COMPLEX        :: f(:,:)
+    integer, intent(in) :: parity
+    
+    CCTK_REAL :: fre(size(f,1),size(f,2))
+    CCTK_REAL :: fim(size(f,1),size(f,2))
+    
+    integer :: ni, nj, gi, gj
+    
+    ni = qlm_ntheta(hn)
+    nj = qlm_nphi(hn)
+    gi = qlm_nghoststheta(hn)
+    gj = qlm_nghostsphi(hn)
+    
+    f(       :gi, :nj) = 0
+    f(ni-gi+1:ni, :nj) = 0
+    f(:ni,        :gj) = 0
+    f(:ni, nj-gj+1:nj) = 0
+    
+    fre = real(f)
+    fim = aimag(f)
+    call set_boundary_real (CCTK_PASS_FTOF, hn, fre, parity)
+    call set_boundary_real (CCTK_PASS_FTOF, hn, fim, parity)
+    f = cmplx(fre, fim, kind(f))
+  end subroutine set_boundary_complex
+  
+end module qlm_boundary