





















































































































 


































































































































































 

























































































































































































 





















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                

















































































module NullInterp_InterpUtil
  use cctk
  implicit none

contains

  subroutine NullInterp_Util_cinterp(cctkGH, lsh_, Jn, Js, gf_indices, spin, stereo_mask)

    use NullGrid_Vars

    implicit none

    interface
  integer*8 function CCTK_PointerTo (var)
  implicit none
  type(*),dimension(..),TARGET :: var
  end function CCTK_PointerTo
  end interface
  interface
  integer function CCTK_Equals (arg1, arg2)
  implicit none
  integer*8 arg1
  character(*) arg2
  end function CCTK_Equals
  integer function CCTK_MyProc (cctkGH)
  implicit none
  integer*8 cctkGH
  end function CCTK_MyProc
  integer function CCTK_nProcs (cctkGH)
  implicit none
  integer*8 cctkGH
  end function CCTK_nProcs
  integer function CCTK_IsThornActive (name)
  implicit none
  character(*) name
  end function CCTK_IsThornActive
  integer*8 function CCTK_NullPointer ()
  implicit none
  end function CCTK_NullPointer
  end interface
    integer*8 stereo_patch_type
  integer, parameter :: cctki_use_stereo_patch_type = kind(stereo_patch_type)
  INTEGER*4 deriv_accuracy
  integer, parameter :: cctki_use_deriv_accuracy = kind(deriv_accuracy)
  INTEGER*4 poison_test
  integer, parameter :: cctki_use_poison_test = kind(poison_test)
  COMMON /NullInterprest/stereo_patch_type, deriv_accuracy, poison_test
  REAL*8 mask_testinterp_tolerance
  integer, parameter :: cctki_use_mask_testinterp_tolerance = kind(mask_testinterp_tolerance)
  REAL*8 poison_value
  integer, parameter :: cctki_use_poison_value = kind(poison_value)
  integer*8 interpolation_name
  integer, parameter :: cctki_use_interpolation_name = kind(interpolation_name)
  INTEGER*4 extended_guard_zone
  integer, parameter :: cctki_use_extended_guard_zone = kind(extended_guard_zone)
  INTEGER*4 interpolation_order
  integer, parameter :: cctki_use_interpolation_order = kind(interpolation_order)
  INTEGER*4 skip_interpolation
  integer, parameter :: cctki_use_skip_interpolation = kind(skip_interpolation)
  INTEGER*4 use_edge_check
  integer, parameter :: cctki_use_use_edge_check = kind(use_edge_check)
  COMMON /NullInterppriv/mask_testinterp_tolerance, poison_value, interpolation_name, extended_guard_zone, interpolation_order, ski&
  &p_interpolation, use_edge_check
  REAL*8 CCTKH1
  integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
  REAL*8 CCTKH2
  integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
  INTEGER*4 N_ang_ev_outside_eq
  integer, parameter :: cctki_use_N_ang_ev_outside_eq = kind(N_ang_ev_outside_eq)
  INTEGER*4 N_ang_ghost_pts
  integer, parameter :: cctki_use_N_ang_ghost_pts = kind(N_ang_ghost_pts)
  INTEGER*4 N_ang_pts_inside_eq
  integer, parameter :: cctki_use_N_ang_pts_inside_eq = kind(N_ang_pts_inside_eq)
  INTEGER*4 N_ang_stencil_size
  integer, parameter :: cctki_use_N_ang_stencil_size = kind(N_ang_stencil_size)
  INTEGER*4 CCTKH0
  integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
  INTEGER*4 CCTKH3
  integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
  COMMON /NULLGRIDrest/CCTKH1, CCTKH2, N_ang_ev_outside_eq, N_ang_ghost_pts, N_ang_pts_inside_eq, N_ang_stencil_size, CCTKH0, CCTKH3

    
    integer*8, intent(in) :: cctkGH
    INTEGER*4, dimension(2), intent(in) :: lsh_
    COMPLEX*16, dimension(lsh_(1),lsh_(2)), intent(inout) :: Jn, Js
    INTEGER*4, dimension(lsh_(1),lsh_(2)), intent(in) :: stereo_mask
    INTEGER*4, dimension(2), intent(in) :: gf_indices
    INTEGER*4, intent(in) :: spin
    INTEGER*4 :: i, j, k

    logical, save :: FirstTime = .true.
    REAL*8, dimension(:), allocatable, save :: qcoord, pcoord
    COMPLEX*16, dimension(:), allocatable, save :: Iarr1, Iarr2
    INTEGER*4, dimension(:), allocatable, save :: reindexi, reindexj
    INTEGER*4, save :: tot_nbp, nbp, mask_count

    
    integer :: ierror, N_in_arrays, N_out_arrays, N_interp_points, N_dims
    integer, save :: interp_handle, coord_system_handle
    integer, save :: param_table_handle = -1
    integer*8, dimension(2) :: interp_coords
    integer :: interp_coords_type_code
    INTEGER*4, dimension(2) :: in_array_indices
    INTEGER*4, dimension(2) :: out_array_type_codes
    integer*8, dimension(2) ::  out_arrays


    if (skip_interpolation.ne.0) then
       if (FirstTime) call CCTK_Info("NullInterp","will skip interpolation");
       FirstTime = .false.
       do i = 1, lsh_(1)
          do j = 1, lsh_(2)
             if (stereo_mask(i,j).eq.1) then
                if (poison_test.ne.0) then
                   Jn(i,j) = poison_value
                   Js(i,j) = poison_value
                else
                   Jn(i,j) = 0
                   Js(i,j) = 0
                end if
             endif
          end do
       end do
       return;
    end if

    
    if ( FirstTime ) then
       FirstTime = .false.

       mask_count = count(stereo_mask.eq.1)
       tot_nbp = mask_count
       allocate(qcoord(tot_nbp), pcoord(tot_nbp), reindexi(tot_nbp),&
               reindexj(tot_nbp), Iarr1(tot_nbp), Iarr2(tot_nbp))

       nbp = 0
       do i = 1, lsh_(1)
          do j = 1, lsh_(2)
             if (stereo_mask(i,j).eq.1) then
                nbp = nbp + 1
                qcoord(nbp) = dble(qs(i,j)/(qs(i,j)**2 + ps(i,j)**2))
                pcoord(nbp) = dble(-ps(i,j)/(qs(i,j)**2 + ps(i,j)**2))
                reindexi(nbp) = i
                reindexj(nbp) = j
             endif
          end do
       end do

       if ( nbp > tot_nbp ) then
          call CCTK_Warn(0,89,"NullInterp_InterpUtil.F90","NullInterp","Boundary point determination error")
       endif

    end if 

    if (param_table_handle < 0) then
       call  Util_TableCreateFromString(param_table_handle,&
            "order = " // char(ichar('0') + interpolation_order))

       if (param_table_handle < 0) then
          call  CCTK_Warn(-1,99,"NullInterp_InterpUtil.F90","NullInterp", "can t create parameter table!")
       end if
    end if

    

    interp_handle =-1
    coord_system_handle = -1


    if ((CCTK_Equals(interpolation_name, "Lagrange") .ne. 0)) then
       call CCTK_InterpHandle (interp_handle, "generalized polynomial interpolation")
    else if ((CCTK_Equals(interpolation_name, "Hermite") .ne. 0)) then
       call CCTK_InterpHandle (interp_handle, "Hermite polynomial interpolation")
    end if

    if (interp_handle < 0) then
       call CCTK_Warn(0,116,"NullInterp_InterpUtil.F90","NullInterp","Interpolation operator not found")
    end if

    call CCTK_CoordSystemHandle (coord_system_handle, "stereo")
    if (coord_system_handle < 0) then
       call CCTK_Warn(0,121,"NullInterp_InterpUtil.F90","NullInterp","Coordinate system 'stereo' not registered")
    end if

    interp_coords(1)        = CCTK_PointerTo(qcoord)
    interp_coords(2)        = CCTK_PointerTo(pcoord)
    interp_coords_type_code = 130
    in_array_indices        = gf_indices
    out_arrays(1)           = CCTK_PointerTo(Iarr1)
    out_arrays(2)           = CCTK_PointerTo(Iarr2)
    out_array_type_codes    = 140

    N_dims=2
    N_in_arrays=2
    N_out_arrays=2
    N_interp_points = nbp

    call CCTK_InterpGridArrays (ierror, cctkGH, N_dims,&
         interp_handle, param_table_handle, coord_system_handle,&
         N_interp_points, interp_coords_type_code,&
         interp_coords,&
         N_in_arrays, in_array_indices,&
         N_out_arrays, out_array_type_codes, out_arrays)

    if (ierror < 0) then
       call CCTK_Warn(1,145,"NullInterp_InterpUtil.F90","NullInterp","Interpolation error")
    endif

    do k = 1, nbp
       i = reindexi(k) 
       j = reindexj(k) 

       if (spin .ne. 0) then
          Jn(i,j) = Iarr1(k) * (-zz(i,j)/conjg(zz(i,j)))**spin
          Js(i,j) = Iarr2(k) * (-zz(i,j)/conjg(zz(i,j)))**spin
       else
          Jn(i,j) = Iarr1(k)
          Js(i,j) = Iarr2(k)
       endif

    end do

  end subroutine NullInterp_Util_cinterp

  subroutine NullInterp_Util_rinterp(cctkGH, lsh_, Jn, Js, gf_indices, stereo_mask)
 
    use NullGrid_Vars

    implicit none

    interface
  integer*8 function CCTK_PointerTo (var)
  implicit none
  type(*),dimension(..),TARGET :: var
  end function CCTK_PointerTo
  end interface
  interface
  integer function CCTK_Equals (arg1, arg2)
  implicit none
  integer*8 arg1
  character(*) arg2
  end function CCTK_Equals
  integer function CCTK_MyProc (cctkGH)
  implicit none
  integer*8 cctkGH
  end function CCTK_MyProc
  integer function CCTK_nProcs (cctkGH)
  implicit none
  integer*8 cctkGH
  end function CCTK_nProcs
  integer function CCTK_IsThornActive (name)
  implicit none
  character(*) name
  end function CCTK_IsThornActive
  integer*8 function CCTK_NullPointer ()
  implicit none
  end function CCTK_NullPointer
  end interface
    integer*8 stereo_patch_type
  integer, parameter :: cctki_use_stereo_patch_type = kind(stereo_patch_type)
  INTEGER*4 deriv_accuracy
  integer, parameter :: cctki_use_deriv_accuracy = kind(deriv_accuracy)
  INTEGER*4 poison_test
  integer, parameter :: cctki_use_poison_test = kind(poison_test)
  COMMON /NullInterprest/stereo_patch_type, deriv_accuracy, poison_test
  REAL*8 mask_testinterp_tolerance
  integer, parameter :: cctki_use_mask_testinterp_tolerance = kind(mask_testinterp_tolerance)
  REAL*8 poison_value
  integer, parameter :: cctki_use_poison_value = kind(poison_value)
  integer*8 interpolation_name
  integer, parameter :: cctki_use_interpolation_name = kind(interpolation_name)
  INTEGER*4 extended_guard_zone
  integer, parameter :: cctki_use_extended_guard_zone = kind(extended_guard_zone)
  INTEGER*4 interpolation_order
  integer, parameter :: cctki_use_interpolation_order = kind(interpolation_order)
  INTEGER*4 skip_interpolation
  integer, parameter :: cctki_use_skip_interpolation = kind(skip_interpolation)
  INTEGER*4 use_edge_check
  integer, parameter :: cctki_use_use_edge_check = kind(use_edge_check)
  COMMON /NullInterppriv/mask_testinterp_tolerance, poison_value, interpolation_name, extended_guard_zone, interpolation_order, ski&
  &p_interpolation, use_edge_check
  REAL*8 CCTKH1
  integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
  REAL*8 CCTKH2
  integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
  INTEGER*4 N_ang_ev_outside_eq
  integer, parameter :: cctki_use_N_ang_ev_outside_eq = kind(N_ang_ev_outside_eq)
  INTEGER*4 N_ang_ghost_pts
  integer, parameter :: cctki_use_N_ang_ghost_pts = kind(N_ang_ghost_pts)
  INTEGER*4 N_ang_pts_inside_eq
  integer, parameter :: cctki_use_N_ang_pts_inside_eq = kind(N_ang_pts_inside_eq)
  INTEGER*4 N_ang_stencil_size
  integer, parameter :: cctki_use_N_ang_stencil_size = kind(N_ang_stencil_size)
  INTEGER*4 CCTKH0
  integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
  INTEGER*4 CCTKH3
  integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
  COMMON /NULLGRIDrest/CCTKH1, CCTKH2, N_ang_ev_outside_eq, N_ang_ghost_pts, N_ang_pts_inside_eq, N_ang_stencil_size, CCTKH0, CCTKH3

    integer*8, intent(in) :: cctkGH
    INTEGER*4, dimension(2), intent(in) :: lsh_
    REAL*8, dimension(lsh_(1),lsh_(2)), intent(inout) :: Jn, Js
    INTEGER*4, dimension(lsh_(1),lsh_(2)), intent(in) :: stereo_mask
    INTEGER*4, dimension(2), intent(in) :: gf_indices

    INTEGER*4 :: i, j, k
    logical, save :: FirstTime = .true.
    REAL*8, dimension(:), allocatable, save :: qcoord, pcoord
    REAL*8, dimension(:), allocatable, save :: Iarr1, Iarr2
    INTEGER*4, dimension(:), allocatable, save :: reindexi, reindexj
    INTEGER*4, save :: tot_nbp, nbp, mask_count

    
    integer :: ierror, N_in_arrays, N_out_arrays, N_interp_points, N_dims
    integer, save :: interp_handle, coord_system_handle
    integer, save :: param_table_handle = -1
    integer*8, dimension(2) :: interp_coords
    integer :: interp_coords_type_code
    INTEGER*4, dimension(2) :: in_array_indices
    INTEGER*4, dimension(2) :: out_array_type_codes
    integer*8, dimension(2) ::  out_arrays

    if (skip_interpolation.ne.0) then
       if (FirstTime) call CCTK_Info("NullInterp","will skip interpolation");
       FirstTime = .false.
       do i = 1, lsh_(1)
          do j = 1, lsh_(2)
             if (stereo_mask(i,j).eq.1) then
                if (poison_test.ne.0) then
                   Jn(i,j) = poison_value
                   Js(i,j) = poison_value
                else
                   Jn(i,j) = 0
                   Js(i,j) = 0
                end if
             endif
          end do
       end do
       return;
    end if

    if ( FirstTime ) then
       FirstTime = .false.

       mask_count = count(stereo_mask.eq.1)
       tot_nbp = mask_count 

       allocate(qcoord(tot_nbp), pcoord(tot_nbp), reindexi(tot_nbp), &
               reindexj(tot_nbp), Iarr1(tot_nbp), Iarr2(tot_nbp))

       nbp = 0
       do i = 1, lsh_(1) 
          do j = 1, lsh_(2) 
             if (stereo_mask(i,j).eq.1) then
                nbp = nbp + 1
                qcoord(nbp) = dble(qs(i,j)/(qs(i,j)**2 + ps(i,j)**2))
                pcoord(nbp) = dble(-ps(i,j)/(qs(i,j)**2 + ps(i,j)**2))
                reindexi(nbp) = i
                reindexj(nbp) = j
             endif
          end do
       end do

       if ( nbp > tot_nbp ) then
          call CCTK_Warn(0,238,"NullInterp_InterpUtil.F90","NullInterp","Boundary point determination error")
       end if

    end if 

    if (param_table_handle < 0) then
       call  Util_TableCreateFromString(param_table_handle,&
            "order = " // char(ichar('0') + interpolation_order))

       if (param_table_handle < 0) then
          call  CCTK_Warn(-1,248,"NullInterp_InterpUtil.F90","NullInterp", "can t create parameter table!")
       end if
    end if

    

    interp_handle =-1
    coord_system_handle = -1


    if ((CCTK_Equals(interpolation_name, "Lagrange") .ne. 0)) then
       call CCTK_InterpHandle (interp_handle, "generalized polynomial interpolation")
    else if ((CCTK_Equals(interpolation_name, "Hermite") .ne. 0)) then
       call CCTK_InterpHandle (interp_handle, "Hermite polynomial interpolation")
    end if

    if (interp_handle < 0) then
       call CCTK_Warn(0,265,"NullInterp_InterpUtil.F90","NullInterp","Interpolation operator not found")
    end if

    call CCTK_CoordSystemHandle (coord_system_handle, "stereo")
    if (coord_system_handle < 0) then
       call CCTK_Warn(0,270,"NullInterp_InterpUtil.F90","NullInterp","Coordinate system 'stereo' not registered")
    end if

    interp_coords(1)        = CCTK_PointerTo(qcoord)
    interp_coords(2)        = CCTK_PointerTo(pcoord)
    interp_coords_type_code = 130
    in_array_indices        = gf_indices
    out_arrays(1)           = CCTK_PointerTo(Iarr1)
    out_arrays(2)           = CCTK_PointerTo(Iarr2)
    out_array_type_codes    = 130

    N_dims=2
    N_in_arrays=2
    N_out_arrays=2
    N_interp_points = nbp

    call CCTK_InterpGridArrays (ierror, cctkGH, N_dims,&
         interp_handle, param_table_handle, coord_system_handle,&
         N_interp_points, interp_coords_type_code,&
         interp_coords,&
         N_in_arrays, in_array_indices,&
         N_out_arrays, out_array_type_codes, out_arrays)

    if (ierror < 0) then
       call CCTK_Warn(1,294,"NullInterp_InterpUtil.F90","NullInterp","Interpolation error")
    endif

    do k = 1, nbp
       i = reindexi(k) 
       j = reindexj(k) 
       Jn(i,j) = Iarr1(k)
       Js(i,j) = Iarr2(k)
    end do

  end subroutine NullInterp_Util_rinterp

  subroutine NullInterp_Util_3cinterp(cctkGH, lsh_, J1n, J1s, J2n, J2s, J3n, J3s, gf_indices, spin1, spin2, spin3, stereo_mask)

    use NullGrid_Vars

    implicit none

    interface
  integer*8 function CCTK_PointerTo (var)
  implicit none
  type(*),dimension(..),TARGET :: var
  end function CCTK_PointerTo
  end interface
  interface
  integer function CCTK_Equals (arg1, arg2)
  implicit none
  integer*8 arg1
  character(*) arg2
  end function CCTK_Equals
  integer function CCTK_MyProc (cctkGH)
  implicit none
  integer*8 cctkGH
  end function CCTK_MyProc
  integer function CCTK_nProcs (cctkGH)
  implicit none
  integer*8 cctkGH
  end function CCTK_nProcs
  integer function CCTK_IsThornActive (name)
  implicit none
  character(*) name
  end function CCTK_IsThornActive
  integer*8 function CCTK_NullPointer ()
  implicit none
  end function CCTK_NullPointer
  end interface
    integer*8 stereo_patch_type
  integer, parameter :: cctki_use_stereo_patch_type = kind(stereo_patch_type)
  INTEGER*4 deriv_accuracy
  integer, parameter :: cctki_use_deriv_accuracy = kind(deriv_accuracy)
  INTEGER*4 poison_test
  integer, parameter :: cctki_use_poison_test = kind(poison_test)
  COMMON /NullInterprest/stereo_patch_type, deriv_accuracy, poison_test
  REAL*8 mask_testinterp_tolerance
  integer, parameter :: cctki_use_mask_testinterp_tolerance = kind(mask_testinterp_tolerance)
  REAL*8 poison_value
  integer, parameter :: cctki_use_poison_value = kind(poison_value)
  integer*8 interpolation_name
  integer, parameter :: cctki_use_interpolation_name = kind(interpolation_name)
  INTEGER*4 extended_guard_zone
  integer, parameter :: cctki_use_extended_guard_zone = kind(extended_guard_zone)
  INTEGER*4 interpolation_order
  integer, parameter :: cctki_use_interpolation_order = kind(interpolation_order)
  INTEGER*4 skip_interpolation
  integer, parameter :: cctki_use_skip_interpolation = kind(skip_interpolation)
  INTEGER*4 use_edge_check
  integer, parameter :: cctki_use_use_edge_check = kind(use_edge_check)
  COMMON /NullInterppriv/mask_testinterp_tolerance, poison_value, interpolation_name, extended_guard_zone, interpolation_order, ski&
  &p_interpolation, use_edge_check
  REAL*8 CCTKH1
  integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
  REAL*8 CCTKH2
  integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
  INTEGER*4 N_ang_ev_outside_eq
  integer, parameter :: cctki_use_N_ang_ev_outside_eq = kind(N_ang_ev_outside_eq)
  INTEGER*4 N_ang_ghost_pts
  integer, parameter :: cctki_use_N_ang_ghost_pts = kind(N_ang_ghost_pts)
  INTEGER*4 N_ang_pts_inside_eq
  integer, parameter :: cctki_use_N_ang_pts_inside_eq = kind(N_ang_pts_inside_eq)
  INTEGER*4 N_ang_stencil_size
  integer, parameter :: cctki_use_N_ang_stencil_size = kind(N_ang_stencil_size)
  INTEGER*4 CCTKH0
  integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
  INTEGER*4 CCTKH3
  integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
  COMMON /NULLGRIDrest/CCTKH1, CCTKH2, N_ang_ev_outside_eq, N_ang_ghost_pts, N_ang_pts_inside_eq, N_ang_stencil_size, CCTKH0, CCTKH3

    
    integer*8,               intent(in) :: cctkGH
    INTEGER*4,     dimension(2), intent(in) :: lsh_
    INTEGER*4,                   intent(in) :: spin1, spin2, spin3

    COMPLEX*16, dimension(lsh_(1),lsh_(2)), intent(inout) :: J1n, J1s,&
         J2n, J2s, J3n, J3s
    INTEGER*4, dimension(lsh_(1),lsh_(2)), intent(in) :: stereo_mask
    INTEGER*4, dimension(6), intent(in) :: gf_indices
    INTEGER*4 :: i, j, k

    logical, save :: FirstTime = .true.
    REAL*8, dimension(:), allocatable, save :: qcoord, pcoord
    COMPLEX*16, dimension(:), allocatable, save :: Iarr1, Iarr2,&
         Iarr3, Iarr4, Iarr5, Iarr6
    INTEGER*4, dimension(:), allocatable, save :: reindexi, reindexj
    INTEGER*4, save :: tot_nbp, nbp, mask_count

    

    
    integer :: ierror, N_in_arrays, N_out_arrays, N_interp_points, N_dims
    integer, save :: interp_handle, param_table_handle, coord_system_handle
    integer*8, dimension(2) :: interp_coords
    integer :: interp_coords_type_code
    INTEGER*4, dimension(6) :: in_array_indices
    INTEGER*4, dimension(6) :: out_array_type_codes
    integer*8, dimension(6) ::  out_arrays

    if (skip_interpolation.ne.0) then
       if (FirstTime) call CCTK_Info("NullInterp","will skip interpolation");
       FirstTime = .false.
       do i = 1, lsh_(1)
          do j = 1, lsh_(2)
             if (stereo_mask(i,j).eq.1) then
                if (poison_test.ne.0) then
                   J1n(i,j) = poison_value
                   J1s(i,j) = poison_value
                else
                   J1n(i,j) = 0
                   J1s(i,j) = 0
                end if
             endif
          end do
       end do
       return;
    end if

    if ( FirstTime ) then
       FirstTime = .false.

       mask_count = count(stereo_mask.eq.1)
       tot_nbp = mask_count
       allocate(qcoord(tot_nbp), pcoord(tot_nbp),&
               reindexi(tot_nbp), reindexj(tot_nbp),&
               Iarr1(tot_nbp), Iarr2(tot_nbp),&
               Iarr3(tot_nbp), Iarr4(tot_nbp),&
               Iarr5(tot_nbp), Iarr6(tot_nbp))

       nbp = 0
       do i = 1, lsh_(1)
          do j = 1, lsh_(2)
             if (stereo_mask(i,j).eq.1) then
                nbp = nbp + 1
                qcoord(nbp) = qs(i,j)/(qs(i,j)**2 + ps(i,j)**2)
                pcoord(nbp) = -ps(i,j)/(qs(i,j)**2 + ps(i,j)**2)
                reindexi(nbp) = i
                reindexj(nbp) = j
             endif
          end do
       end do

       if ( nbp > tot_nbp ) then
          call CCTK_Warn(0,388,"NullInterp_InterpUtil.F90","NullInterp","Boundary point determination error")
       endif

       param_table_handle = -1
       call  Util_TableCreateFromString(param_table_handle,&
            "order = " // char(ichar('0') + interpolation_order))

       if (param_table_handle < 0) then
          call  CCTK_Warn(-1,396,"NullInterp_InterpUtil.F90","NullInterp", "can t create parameter table!")
       end if

    end if 

    

















    

    interp_handle =-1
    coord_system_handle = -1


    if ((CCTK_Equals(interpolation_name, "Lagrange") .ne. 0)) then
       call CCTK_InterpHandle (interp_handle, "generalized polynomial interpolation")
    else if ((CCTK_Equals(interpolation_name, "Hermite") .ne. 0)) then
       call CCTK_InterpHandle (interp_handle, "Hermite polynomial interpolation")
    end if

    if (interp_handle < 0) then
       call CCTK_Warn(0,432,"NullInterp_InterpUtil.F90","NullInterp","Interpolation operator not found")
    end if

    call CCTK_CoordSystemHandle (coord_system_handle, "stereo")
    if (coord_system_handle < 0) then
       call CCTK_Warn(0,437,"NullInterp_InterpUtil.F90","NullInterp","Coordinate system 'stereo' not registered")
    end if

    interp_coords(1)        = CCTK_PointerTo(qcoord)
    interp_coords(2)        = CCTK_PointerTo(pcoord)
    interp_coords_type_code = 130
    in_array_indices        = gf_indices
    out_arrays(1)           = CCTK_PointerTo(Iarr1)
    out_arrays(2)           = CCTK_PointerTo(Iarr2)
    out_arrays(3)           = CCTK_PointerTo(Iarr3)
    out_arrays(4)           = CCTK_PointerTo(Iarr4)
    out_arrays(5)           = CCTK_PointerTo(Iarr5)
    out_arrays(6)           = CCTK_PointerTo(Iarr6)
    out_array_type_codes    = 140

    N_dims=2
    N_in_arrays=6
    N_out_arrays=6
    N_interp_points = nbp

    call CCTK_InterpGridArrays (ierror, cctkGH, N_dims,&
         interp_handle, param_table_handle, coord_system_handle,&
         N_interp_points, interp_coords_type_code,&
         interp_coords,&
         N_in_arrays, in_array_indices,&
         N_out_arrays, out_array_type_codes, out_arrays)

    if (ierror < 0) then
       call CCTK_Warn(1,465,"NullInterp_InterpUtil.F90","NullInterp","Interpolation error")
    endif

    do k = 1, nbp
       i = reindexi(k) 
       j = reindexj(k) 

       if (spin1 .ne. 0) then
          J1n(i,j) = Iarr1(k) * (-zz(i,j)/conjg(zz(i,j)))**spin1
          J1s(i,j) = Iarr2(k) * (-zz(i,j)/conjg(zz(i,j)))**spin1
       else
          J1n(i,j) = Iarr1(k)
          J1s(i,j) = Iarr2(k)
       endif

       if (spin2 .ne. 0) then
          J2n(i,j) = Iarr3(k) * (-zz(i,j)/conjg(zz(i,j)))**spin2
          J2s(i,j) = Iarr4(k) * (-zz(i,j)/conjg(zz(i,j)))**spin2
       else
          J2n(i,j) = Iarr3(k)
          J2s(i,j) = Iarr4(k)
       endif

       if (spin3 .ne. 0) then
          J3n(i,j) = Iarr5(k) * (-zz(i,j)/conjg(zz(i,j)))**spin3
          J3s(i,j) = Iarr6(k) * (-zz(i,j)/conjg(zz(i,j)))**spin3
       else
          J3n(i,j) = Iarr5(k)
          J3s(i,j) = Iarr6(k)
       endif


       
       
       
       
       
       

    end do

  end subroutine NullInterp_Util_3cinterp

end module NullInterp_InterpUtil
