





















































































































 


































































































































































 
















































































































































































































 























































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                

















































































module NullNews_Bondi
  implicit none
contains

  subroutine NullNews_News2BondiNews(cctkGH, lsh, dt, qsize, Uo, Un, Uyo, Uyn,&
       zEvolo, zEvoln, zEvolh, News, NewsB, uBondio,&
       uBondi, uBondin, redshiftB, omega, Jh, eth_J,&
       beth_J, eth_U, beth_U, du_J, beta, patch_index,&
       deltao, deltan, nzeta, n_tmp1_cgfn, n_tmp1_cgfs,&
       n_tmp2_cgfn, n_tmp2_cgfs, n_tmp3_cgfn,&
       n_tmp3_cgfs, n_tmp1_rgfn, n_tmp1_rgfs,&
       n_tmp2_rgfn, n_tmp2_rgfs, starting, stereo_mask)              
    implicit none

    integer*8,                intent(in) :: cctkGH
    INTEGER*4,     dimension(2),  intent(in) :: lsh
    REAL*8,                   intent(in) :: dt, qsize

    INTEGER*4, dimension(lsh(1),lsh(2)), intent(in) :: stereo_mask

    COMPLEX*16, dimension (lsh(1), lsh(2),2), intent (in)    :: Uo, Un, Jh, eth_J, beth_J
    COMPLEX*16, dimension (lsh(1), lsh(2),2), intent (in)    :: eth_U, beth_u, du_J
    COMPLEX*16, dimension (lsh(1), lsh(2),2), intent (inout) :: Uyo, Uyn
    COMPLEX*16, dimension (lsh(1), lsh(2),2), intent (inout) :: zEvolo, zEvoln, zEvolh, nzeta
    COMPLEX*16, dimension (lsh(1), lsh(2),2), intent (in)    :: News
    COMPLEX*16, dimension (lsh(1), lsh(2),2), intent (out)   :: NewsB
    REAL*8,    dimension (lsh(1), lsh(2),2), intent (in)    :: omega, beta
    REAL*8,    dimension (lsh(1), lsh(2),2), intent (inout) :: uBondi, uBondin, uBondio
    REAL*8,    dimension (lsh(1), lsh(2),2), intent (inout) :: redshiftB
    REAL*8,    dimension (lsh(1), lsh(2),2), intent (inout) :: deltao, deltan
    INTEGER*4,     dimension (lsh(1), lsh(2),2), intent (inout) :: patch_index

    COMPLEX*16, dimension(lsh(1), lsh(2)), intent(inout) ::&
         n_tmp1_cgfn, n_tmp1_cgfs,&
         n_tmp2_cgfn, n_tmp2_cgfs,&
         n_tmp3_cgfn, n_tmp3_cgfs

    REAL*8,    dimension(lsh(1), lsh(2)), intent(inout) ::&
         n_tmp1_rgfn, n_tmp1_rgfs,&
         n_tmp2_rgfn, n_tmp2_rgfs 

    INTEGER*4, intent(inout) :: starting

    COMPLEX*16, dimension (:,:,:), allocatable, save :: Uh, Uyh, Po, Pm
    COMPLEX*16, dimension (:,:,:), allocatable, save :: delta_RHS, delta_RHS_B
    REAL*8,    dimension (:,:,:), allocatable, save :: delta, K, betaB, omegaB

    REAL*8, save :: pii = 3.1415926535897932385d0
    INTEGER*4 ::  i, j, ip, ipo, ipm
    REAL*8 qh, ph
    COMPLEX*16 :: ii  = (0.d0, 1.d0)
    INTEGER*4 :: p_this, p_other, l, ierror
    logical, save :: FirstTime = .true.
    
    integer :: N_in_arrays, N_out_arrays, N_interp_points, N_dims
    integer, save :: interp_handle, coord_system_handle, param_table_handle
    integer*8, dimension(2) :: interp_coords
    integer                     :: interp_coords_type_code
    INTEGER*4,     dimension(5) :: in_array_indices
    integer*8, dimension(5) :: out_arrays
    INTEGER*4,     dimension(5) :: out_array_type_codes
    integer,      dimension(10), save :: gf_indices


    COMPLEX*16, dimension(:), allocatable, save ::&
         Iarr1_this, Iarr1_other,&
         Iarr2_this, Iarr2_other,&
         Iarr4_this, Iarr4_other

    REAL*8, dimension(:), allocatable, save ::&
         Rarr3_this, Rarr3_other,&
         Rarr5_this, Rarr5_other,&
         qcoord_this, qcoord_other,&
         pcoord_this, pcoord_other

    INTEGER*4, dimension(:), allocatable, save ::&
         reindexi_this, reindexi_other,&
         reindexj_this, reindexj_other


    character(len=500) :: message

    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
  interface
  INTEGER*4 function IO_TruncateOutputFiles (GH)
  implicit none
  integer*8 GH
  end function IO_TruncateOutputFiles
  end interface
    INTEGER*4 linearized_news
  integer, parameter :: cctki_use_linearized_news = kind(linearized_news)
  COMMON /NullNewsrest/linearized_news
  REAL*8 eps_omega
  integer, parameter :: cctki_use_eps_omega = kind(eps_omega)
  INTEGER*4 DEBUG_skip_BondiNews
  integer, parameter :: cctki_use_DEBUG_skip_BondiNews = kind(DEBUG_skip_BondiNews)
  INTEGER*4 Jl_deriv_order
  integer, parameter :: cctki_use_Jl_deriv_order = kind(Jl_deriv_order)
  INTEGER*4 NewsSpinWeight
  integer, parameter :: cctki_use_NewsSpinWeight = kind(NewsSpinWeight)
  INTEGER*4 Ylm_at_Scri
  integer, parameter :: cctki_use_Ylm_at_Scri = kind(Ylm_at_Scri)
  INTEGER*4 compute_lin_strain
  integer, parameter :: cctki_use_compute_lin_strain = kind(compute_lin_strain)
  INTEGER*4 debug_output
  integer, parameter :: cctki_use_debug_output = kind(debug_output)
  INTEGER*4 filter_news
  integer, parameter :: cctki_use_filter_news = kind(filter_news)
  INTEGER*4 filter_omega
  integer, parameter :: cctki_use_filter_omega = kind(filter_omega)
  INTEGER*4 filter_scri_fields
  integer, parameter :: cctki_use_filter_scri_fields = kind(filter_scri_fields)
  INTEGER*4 first_order_scheme
  integer, parameter :: cctki_use_first_order_scheme = kind(first_order_scheme)
  INTEGER*4 interp_to_constant_uBondi
  integer, parameter :: cctki_use_interp_to_constant_uBondi = kind(interp_to_constant_uBondi)
  INTEGER*4 l_YlmScri
  integer, parameter :: cctki_use_l_YlmScri = kind(l_YlmScri)
  INTEGER*4 linearized_inertial_frame
  integer, parameter :: cctki_use_linearized_inertial_frame = kind(linearized_inertial_frame)
  INTEGER*4 m_YlmScri
  integer, parameter :: cctki_use_m_YlmScri = kind(m_YlmScri)
  INTEGER*4 mask_NewsB
  integer, parameter :: cctki_use_mask_NewsB = kind(mask_NewsB)
  INTEGER*4 mask_Psi4
  integer, parameter :: cctki_use_mask_Psi4 = kind(mask_Psi4)
  INTEGER*4 max_timelevels
  integer, parameter :: cctki_use_max_timelevels = kind(max_timelevels)
  INTEGER*4 news_interp_order
  integer, parameter :: cctki_use_news_interp_order = kind(news_interp_order)
  INTEGER*4 use_linearized_omega
  integer, parameter :: cctki_use_use_linearized_omega = kind(use_linearized_omega)
  INTEGER*4 write_spherical_harmonics
  integer, parameter :: cctki_use_write_spherical_harmonics = kind(write_spherical_harmonics)
  COMMON /NullNewspriv/eps_omega, DEBUG_skip_BondiNews, Jl_deriv_order, NewsSpinWeight, Ylm_at_Scri, compute_lin_strain, debug_outp&
  &ut, filter_news, filter_omega, filter_scri_fields, first_order_scheme, interp_to_constant_uBondi, l_YlmScri, linearized_inertial&
  &_frame, m_YlmScri, mask_NewsB, mask_Psi4, max_timelevels, news_interp_order, use_linearized_omega, write_spherical_harmonics
  REAL*8 CCTKH0
  integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
  REAL*8 CCTKH1
  integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
  REAL*8 CCTKH2
  integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
  integer*8 CCTKH3
  integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
  integer*8 CCTKH4
  integer, parameter :: cctki_use_CCTKH4 = kind(CCTKH4)
  INTEGER*4 cctk_itlast
  integer, parameter :: cctki_use_cctk_itlast = kind(cctk_itlast)
  INTEGER*4 CCTKH5
  integer, parameter :: cctki_use_CCTKH5 = kind(CCTKH5)
  COMMON /CACTUSrest/CCTKH0, CCTKH1, CCTKH2, CCTKH3, CCTKH4, cctk_itlast, CCTKH5
  REAL*8 null_rwt
  integer, parameter :: cctki_use_null_rwt = kind(null_rwt)
  REAL*8 CCTKH6
  integer, parameter :: cctki_use_CCTKH6 = kind(CCTKH6)
  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 N_radial_pts
  integer, parameter :: cctki_use_N_radial_pts = kind(N_radial_pts)
  INTEGER*4 CCTKH7
  integer, parameter :: cctki_use_CCTKH7 = kind(CCTKH7)
  COMMON /NULLGRIDrest/null_rwt, CCTKH6, N_ang_ev_outside_eq, N_ang_ghost_pts, N_ang_pts_inside_eq, N_ang_stencil_size, N_radial_pt&
  &s, CCTKH7

    

    
    

    if (FirstTime) then
       allocate (Uh(lsh(1), lsh(2),2), Uyh(lsh(1), lsh(2),2),&
            Po(lsh(1), lsh(2),2), Pm(lsh(1), lsh(2),2),&
            delta_RHS(lsh(1), lsh(2),2), delta_RHS_B(lsh(1), lsh(2),2),&
            delta(lsh(1), lsh(2),2), K(lsh(1), lsh(2),2),&
            betaB(lsh(1), lsh(2),2), omegaB(lsh(1), lsh(2),2) )

       allocate(Iarr1_this(lsh(1)*lsh(2)), Iarr1_other(lsh(1)*lsh(2)),&
            Iarr2_this(lsh(1)*lsh(2)), Iarr2_other(lsh(1)*lsh(2)),&
            Iarr4_this(lsh(1)*lsh(2)), Iarr4_other(lsh(1)*lsh(2)),&
            Rarr3_this(lsh(1)*lsh(2)), Rarr3_other(lsh(1)*lsh(2)),&
            Rarr5_this(lsh(1)*lsh(2)), Rarr5_other(lsh(1)*lsh(2)),&
            qcoord_this(lsh(1)*lsh(2)), qcoord_other(lsh(1)*lsh(2)),&
            pcoord_this(lsh(1)*lsh(2)), pcoord_other(lsh(1)*lsh(2)),&
            reindexi_this(lsh(1)*lsh(2)), reindexi_other(lsh(1)*lsh(2)),&
            reindexj_this(lsh(1)*lsh(2)), reindexj_other(lsh(1)*lsh(2))  )

       Uh=0; Uyh=0; Po=0; Pm=0; delta_RHS=0; delta_RHS_B=0
       delta=0; K=0; betaB=0; omegaB=0; Iarr1_this=0; Iarr1_other=0
       Iarr2_this=0; Iarr2_other=0; Iarr4_this=0; Iarr4_other=0
       Rarr3_this=0; Rarr3_other=0
       Rarr5_this=0; Rarr5_other=0; qcoord_this=0; qcoord_other=0
       pcoord_this=0; pcoord_other=0; reindexi_this=0; reindexi_other=0
       reindexj_this=0; reindexj_other=0

       FirstTime = .false.
       if (starting .eq. 1) then
          starting = 0
          Uyo = Uo
          
          do ip = 1, 2
             do j = 1, lsh(2)
                do i = 1, lsh(1)
                   if (dble(zEvolo(i,j,ip)*conjg(zEvolo(i,j,ip)))>1.0d0) then
                      patch_index(i,j,ip) = 1
                      deltao(i,j,ip) = pii -2.0d0 *atan2(dimag(zEvolo(i,j,ip)), dble(zEvolo(i,j,ip)))
                      zEvolo(i,j,ip) = 1.0d0 / zEvolo(i,j,ip)
                      Uyo(i,j,ip) = Uyo(i,j,ip) * (-zEvolo(i,j,ip) /conjg(zEvolo(i,j,ip)))
                   end if
                end do
             end do
          end do
       end if


       call CCTK_VarIndex(gf_indices(1),  "NullNews::n_tmp1_cgfn")
       call CCTK_VarIndex(gf_indices(2),  "NullNews::n_tmp1_cgfs")
       call CCTK_VarIndex(gf_indices(3),  "NullNews::n_tmp2_cgfn")
       call CCTK_VarIndex(gf_indices(4),  "NullNews::n_tmp2_cgfs")
       call CCTK_VarIndex(gf_indices(5),  "NullNews::n_tmp1_rgfn")
       call CCTK_VarIndex(gf_indices(6),  "NullNews::n_tmp1_rgfs")
       call CCTK_VarIndex(gf_indices(7),  "NullNews::n_tmp3_cgfn")
       call CCTK_VarIndex(gf_indices(8),  "NullNews::n_tmp3_cgfs")
       call CCTK_VarIndex(gf_indices(9),  "NullNews::n_tmp2_rgfn")
       call CCTK_VarIndex(gf_indices(10), "NullNews::n_tmp2_rgfs")

       
       
       
       
       
       

       interp_handle =-1
       param_table_handle = -1
       coord_system_handle = -1

       if(minval(gf_indices) < 0 ) then
          call CCTK_Warn(0,164,"NullNews_Bondi.F90","NullNews", "Error Obtaining VarIndex")
       end if
       call CCTK_InterpHandle (interp_handle, "generalized polynomial interpolation")
       if (interp_handle < 0) then
          call CCTK_Warn(0,168,"NullNews_Bondi.F90","NullNews","Interpolation operator not found")
       end if

       call CCTK_CoordSystemHandle (coord_system_handle, "stereo")
       if (coord_system_handle < 0) then
          call CCTK_Warn(0,173,"NullNews_Bondi.F90","NullNews","Coordinate system 'stereo' not registered")
       end if

       call  Util_TableCreateFromString(param_table_handle, &
             "order = " // char(ichar('0') + news_interp_order))
       if (param_table_handle < 0) then
          call  CCTK_Warn(-1,179,"NullNews_Bondi.F90","NullNews", "Can t create parameter table!")
       end if

    end if 

    


    Po = 1.0d0+ zEvolo*conjg(zEvolo)
    zEvoln = zEvolo + 0.25d0 * Po * Uyo * dt
    Uh = 0.5d0 * (Un + Uo)
    Pm = 1.0d0+ zEvoln*conjg(zEvoln)
    
    
    
    
    
    
    

    if (linearized_inertial_frame .eq. 0) then
       Uyh = 0

       do ip = 1, 2

          p_this = 0
          p_other = 0
          do j = 1, lsh(2)
             do i = 1, lsh(1)
                qh = dble(zEvoln(i,j,ip))
                ph = dimag(zEvoln(i,j,ip))
                ipm = mod( ip + patch_index(i,j,ip) - 1 ,2 ) +1
                if (abs(qh) > 1.2 .OR. abs(ph) > 1.2 ) then
                   CALL CCTK_Warn(0,212,"NullNews_Bondi.F90","NullNews","Abs qh or abs ph > 1.2")
                endif
              if (stereo_mask(i,j).eq.1) then
                if ( ipm .eq. ip ) then
                   p_this = p_this + 1
                   qcoord_this(p_this) = qh
                   pcoord_this(p_this) = ph
                   reindexi_this(p_this) = i
                   reindexj_this(p_this) = j
                else
                   p_other = p_other + 1
                   qcoord_other(p_other) = qh
                   pcoord_other(p_other) = ph
                   reindexi_other(p_other) = i
                   reindexj_other(p_other) = j
                end if
              end if
             end do
          end do
          ipo = mod( ip, 2 ) +1
          n_tmp1_cgfn = Uh(1:lsh(1), 1:lsh(2), ip)  
          n_tmp1_cgfs = Uh(1:lsh(1), 1:lsh(2), ipo)

          
          interp_coords(1)        = CCTK_PointerTo(qcoord_this)
          interp_coords(2)        = CCTK_PointerTo(pcoord_this)
          in_array_indices(1)     = gf_indices(1)
          out_arrays(1)           = CCTK_PointerTo(Iarr1_this)
          interp_coords_type_code = 130
          out_array_type_codes(1) = 140

          N_dims=2
          N_in_arrays=1
          N_out_arrays=1
          N_interp_points = p_this

          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(1:1),&
               N_out_arrays, out_array_type_codes(1:1), out_arrays(1:1))

          
          
          
          
          
          

          if (ierror .ne. 0 ) then
             call CCTK_Warn(0,263,"NullNews_Bondi.F90","NullNews", "News Interp Error")
          endif

          
          interp_coords(1)        = CCTK_PointerTo(qcoord_other)
          interp_coords(2)        = CCTK_PointerTo(pcoord_other)
          in_array_indices(2)     = gf_indices(2)
          out_arrays(2)           = CCTK_PointerTo(Iarr1_other)
          interp_coords_type_code = 130
          out_array_type_codes(2) = 140

          N_dims=2
          N_in_arrays=1
          N_out_arrays=1
          N_interp_points = p_other

          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(2:2),&
               N_out_arrays, out_array_type_codes(2:2), out_arrays(2:2))

          if (ierror .ne. 0 ) then
             call CCTK_Warn(0,287,"NullNews_Bondi.F90","NullNews", "News Interp Error")
          endif

          do l = 1, p_this
             Uyh(reindexi_this(l), reindexj_this(l), ip) = Iarr1_this(l)
          end do
          do l = 1, p_other
             Uyh(reindexi_other(l), reindexj_other(l), ip) = Iarr1_other(l)
          end do

       end do 

       
       

       zEvoln = zEvolo + .5*Pm * Uyh * dt
       zEvolh = 0.5d0 * (zEvoln + zEvolo)

       
       
       
       

       
       
       


       Uyn = 0

       do ip = 1, 2

          p_this = 0
          p_other = 0
          do j = 1, lsh(2)
             do i = 1, lsh(1)
                qh = dble(zEvoln(i,j,ip))
                ph = dimag(zEvoln(i,j,ip))
                if (abs(qh) > 1.2 .OR. abs(ph) > 1.2 ) then
                   write(message, *) "CFL_VIOALTION_NEWS: ", &
                        zEvoln(i,j,ip), zEvolo(i,j,ip), Uyh(i,j,ip) * dt
                   CALL CCTK_Warn(0,328,"NullNews_Bondi.F90","NullNews",trim(message))
                endif
                ipm = mod( ip + patch_index(i,j,ip) - 1 ,2 ) +1
              if (stereo_mask(i,j).eq.1) then
                if ( ipm .eq. ip ) then
                   p_this = p_this + 1
                   qcoord_this(p_this) = qh
                   pcoord_this(p_this) = ph
                   reindexi_this(p_this) = i
                   reindexj_this(p_this) = j
                else
                   p_other = p_other + 1
                   qcoord_other(p_other) = qh
                   pcoord_other(p_other) = ph
                   reindexi_other(p_other) = i
                   reindexj_other(p_other) = j
                end if
              end if
             end do
          end do
          ipo = mod( ip, 2 ) +1
          n_tmp1_cgfn = Un(1:lsh(1), 1:lsh(2), ip)  
          n_tmp1_cgfs = Un(1:lsh(1), 1:lsh(2), ipo)

          
          interp_coords(1)        = CCTK_PointerTo(qcoord_this)
          interp_coords(2)        = CCTK_PointerTo(pcoord_this)
          in_array_indices(1)     = gf_indices(1)
          out_arrays(1)           = CCTK_PointerTo(Iarr1_this)
          interp_coords_type_code = 130
          out_array_type_codes(1) = 140

          N_dims=2
          N_in_arrays=1
          N_out_arrays=1
          N_interp_points = p_this

          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(1:1),&
               N_out_arrays, out_array_type_codes(1:1), out_arrays(1:1))


          if (ierror .ne. 0 ) then
             call CCTK_Warn(0,374,"NullNews_Bondi.F90","NullNews", "News Interp Error")
          endif

          
          interp_coords(1)        = CCTK_PointerTo(qcoord_other)
          interp_coords(2)        = CCTK_PointerTo(pcoord_other)
          in_array_indices(2)     = gf_indices(2)
          out_arrays(2)           = CCTK_PointerTo(Iarr1_other)
          interp_coords_type_code = 130
          out_array_type_codes(2) = 140

          N_dims=2
          N_in_arrays=1
          N_out_arrays=1
          N_interp_points = p_other

          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(2:2),&
               N_out_arrays, out_array_type_codes(2:2), out_arrays(2:2))

          if (ierror .ne. 0 ) then
             call CCTK_Warn(0,398,"NullNews_Bondi.F90","NullNews", "News Interp Error")
          endif

          do l = 1, p_this
             Uyn(reindexi_this(l), reindexj_this(l), ip) = Iarr1_this(l)
          end do
          do l = 1, p_other
             Uyn(reindexi_other(l), reindexj_other(l), ip) = Iarr1_other(l)
          end do

       end do 




       K = dsqrt(1.0d0 + dble(Jh*conjg(Jh)))
       delta_RHS =  conjg(du_J)*Jh / (K + 1.0d0) + &
            .5d0 /(K + 1.0d0) * (Jh*(Uh*conjg(eth_J)+&
            conjg(Uh)*conjg(beth_J))) +&
            Jh*conjg(eth_U) + K*beth_U


       
       

       Uyh = 0
       delta_RHS_B=0
       betaB = 0
       NewsB = 0
       OmegaB = 1

       do ip = 1, 2
          p_this = 0
          p_other = 0
          do j = 1, lsh(2)
             do i = 1, lsh(1)
                qh = dble(zEvolh(i,j,ip))
                ph = dimag(zEvolh(i,j,ip))
                ipm = mod( ip + patch_index(i,j,ip) - 1 ,2 ) +1
              if (stereo_mask(i,j).eq.1) then
                if ( ipm .eq. ip ) then
                   p_this = p_this + 1
                   qcoord_this(p_this) = qh
                   pcoord_this(p_this) = ph
                   reindexi_this(p_this) = i
                   reindexj_this(p_this) = j
                else
                   p_other = p_other + 1
                   qcoord_other(p_other) = qh
                   pcoord_other(p_other) = ph
                   reindexi_other(p_other) = i
                   reindexj_other(p_other) = j
                end if
              end if
             end do
          end do
          ipo = mod( ip, 2 ) +1

          n_tmp1_cgfn = Uh(1:lsh(1), 1:lsh(2), ip)  
          n_tmp1_cgfs = Uh(1:lsh(1), 1:lsh(2), ipo)
          n_tmp2_cgfn = delta_RHS(1:lsh(1), 1:lsh(2), ip)
          n_tmp2_cgfs = delta_RHS(1:lsh(1), 1:lsh(2), ipo)
          n_tmp1_rgfn = beta(1:lsh(1), 1:lsh(2), ip)
          n_tmp1_rgfs = beta(1:lsh(1), 1:lsh(2), ipo)
          n_tmp3_cgfn = News(1:lsh(1), 1:lsh(2), ip)
          n_tmp3_cgfs = News(1:lsh(1), 1:lsh(2), ipo)
          n_tmp2_rgfn = omega(1:lsh(1), 1:lsh(2), ip)
          n_tmp2_rgfs = omega(1:lsh(1), 1:lsh(2), ipo)
          

          
          interp_coords(1)        = CCTK_PointerTo(qcoord_this)
          interp_coords(2)        = CCTK_PointerTo(pcoord_this)


          in_array_indices(1)     = gf_indices(1)
          in_array_indices(2)     = gf_indices(3)
          in_array_indices(3)     = gf_indices(7)
          out_arrays(1)           = CCTK_PointerTo(Iarr1_this)
          out_arrays(2)           = CCTK_PointerTo(Iarr2_this)
          out_arrays(3)           = CCTK_PointerTo(Iarr4_this)
          interp_coords_type_code = 130
          out_array_type_codes(1) = 140
          out_array_type_codes(2) = 140
          out_array_type_codes(3) = 140

          N_dims=2
          N_in_arrays=3
          N_out_arrays=3
          N_interp_points = p_this

          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(1:3),&
               N_out_arrays, out_array_type_codes(1:3), out_arrays(1:3))

          if (ierror .ne. 0 ) then
             call CCTK_Warn(0,497,"NullNews_Bondi.F90","NullNews", "News Interp Error")
          endif

          in_array_indices(1)     = gf_indices(5)
          in_array_indices(2)     = gf_indices(9)
          out_arrays(1)           = CCTK_PointerTo(Rarr3_this)
          out_arrays(2)           = CCTK_PointerTo(Rarr5_this)
          interp_coords_type_code = 130
          out_array_type_codes(1) = 130
          out_array_type_codes(2) = 130

          N_dims=2
          N_in_arrays=2
          N_out_arrays=2
          N_interp_points = p_this

          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(1:2),&
               N_out_arrays, out_array_type_codes(1:2), out_arrays(1:2))

          if (ierror .ne. 0 ) then
             call CCTK_Warn(0,521,"NullNews_Bondi.F90","NullNews", "News Interp Error")
          endif

          
          interp_coords(1)        = CCTK_PointerTo(qcoord_other)
          interp_coords(2)        = CCTK_PointerTo(pcoord_other)
          in_array_indices(1)     = gf_indices(2)
          in_array_indices(2)     = gf_indices(4)
          in_array_indices(3)     = gf_indices(8)
          out_arrays(1)           = CCTK_PointerTo(Iarr1_other)
          out_arrays(2)           = CCTK_PointerTo(Iarr2_other)
          out_arrays(3)           = CCTK_PointerTo(Iarr4_other)
          interp_coords_type_code = 130
          out_array_type_codes(1) = 140
          out_array_type_codes(2) = 140
          out_array_type_codes(3) = 140

          N_dims=2
          N_in_arrays=3
          N_out_arrays=3
          N_interp_points = p_other

          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(1:3),&
               N_out_arrays, out_array_type_codes(1:3), out_arrays(1:3))

          if (ierror .ne. 0 ) then
             call CCTK_Warn(0,551,"NullNews_Bondi.F90","NullNews", "News Interp Error")
          endif

          
          interp_coords(1)        = CCTK_PointerTo(qcoord_other)
          interp_coords(2)        = CCTK_PointerTo(pcoord_other)
          in_array_indices(1)     = gf_indices(6)
          in_array_indices(2)     = gf_indices(10)
          out_arrays(1)           = CCTK_PointerTo(Rarr3_other)
          out_arrays(2)           = CCTK_PointerTo(Rarr5_other)
          interp_coords_type_code = 130
          out_array_type_codes(1) = 130
          out_array_type_codes(2) = 130

          N_dims=2
          N_in_arrays=2
          N_out_arrays=2
          N_interp_points = p_other

          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(1:2),&
               N_out_arrays, out_array_type_codes(1:2), out_arrays(1:2))

          if (ierror .ne. 0 ) then
             call CCTK_Warn(0,578,"NullNews_Bondi.F90","NullNews", "News Interp Error")
          endif

          do l = 1, p_this
             i = reindexi_this(l)
             j = reindexj_this(l)
             Uyh(i, j, ip)         = Iarr1_this(l)
             delta_RHS_B(i, j, ip) = Iarr2_this(l)
             betaB(i, j, ip)       = Rarr3_this(l)
             NewsB(i, j, ip)       = Iarr4_this(l)
             OmegaB(i, j, ip)      = Rarr5_this(l)

          end do
          do l = 1, p_other
             i = reindexi_other(l)
             j = reindexj_other(l)
             Uyh(i, j, ip)         = Iarr1_other(l)
             delta_RHS_B(i, j, ip) = Iarr2_other(l)
             betaB(i, j, ip)       = Rarr3_other(l)
             NewsB(i, j, ip)       = Iarr4_other(l)
             OmegaB(i, j, ip)      = Rarr5_other(l)
          end do

       end do


       delta_RHS_B = .5d0 * delta_RHS_B  +  uyh*conjg(zEvolh) 
       deltan = deltao + dt * dimag(delta_RHS_B)

    else 

       zEvoln = nzeta
       Uyh = Uh
       Uyn = Un
       NewsB = News
       deltao  = 0
       betaB = beta
       omegaB = omega
       zEvoln = nzeta
       deltan = 0

    endif 

    redshiftB = omegaB * exp(2*betaB) 
    uBondin = uBondio + dt * omegaB * exp(2*betaB)
    uBondi = .5d0 * (uBondin + uBondio)


    delta = .5d0 * (deltan + deltao)
    NewsB = NewsB * ( cos(2.0d0*delta) - ii*sin(2.0d0*delta) )

    

    zEvolo = zEvoln
    Uyo = Uyn
    deltao = deltan
    uBondio = uBondin

    do ip = 1, 2
       do j = 1, lsh(2)
          do i = 1, lsh(1)
             if ( dble(zEvolo(i,j,ip) * conjg(zEvolo(i,j,ip))) > 1.0d0 ) then
                patch_index(i,j,ip) = mod(patch_index(i,j,ip) +1, 2)  
                
                deltao(i,j,ip) = deltao(i,j,ip) -2.0d0 *atan2(dimag(zEvolo(i,j,ip)), dble(zEvolo(i,j,ip)))&
                     + pii - 2 * pii * nint(deltao(i,j,ip)/(2.0d0*pii)) 
                zEvolo(i,j,ip) = 1.0d0 / zEvolo(i,j,ip)
                Uyo(i,j,ip) = Uyo(i,j,ip) * (-zEvolo(i,j,ip) /conjg(zEvolo(i,j,ip)))
              
             end if
          end do
       end do
    end do


  end subroutine NullNews_News2BondiNews


end module NullNews_Bondi
