






















































































































 


































































































































































 
























































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                





































module NullConstr_Util

contains



subroutine x2r_derivs(dr_dx,d2r_dx2,dr_dxh,d2r_dx2h,null_nx,xb,xbh)
  use NullGrid_Vars, only: rwt
  implicit none
  INTEGER*4, intent(in) :: null_nx
  REAL*8, dimension(null_nx), intent(in) :: xb, xbh
  REAL*8, dimension(null_nx), intent(out) :: dr_dx,d2r_dx2,dr_dxh,d2r_dx2h
  
  

  dr_dx=rwt/(1-xb)**2
  d2r_dx2=2*rwt/(1-xb)**3

  dr_dxh=rwt/(1-xbh)**2
  d2r_dx2h=2*rwt/(1-xbh)**3
 
  
 
  dr_dx(null_nx)=rwt/(1-xb(null_nx-1))**2
  d2r_dx2(null_nx)=2*rwt/(1-xb(null_nx-1))**3

  dr_dxh(null_nx)=rwt/(1-xbh(null_nx-1))**2
  d2r_dx2h(null_nx)=2*rwt/(1-xbh(null_nx-1))**3
  
end subroutine x2r_derivs

subroutine real_derivs (aa_n,aa_p,dt,dx,null_nx,null_lsh,i,dr_dx,d2r_dx2, &
  aa_00,aa_01,aa_02,aa_03,aa_04,aa_11,aa_12,aa_13,aa_14, &
  aa_22,aa_23,aa_24,aa_33,aa_34)
  use NullInterp
  implicit none
  
  INTEGER*4, parameter :: izero = 0
  integer, parameter :: ik = kind(izero)

  INTEGER*4, intent(in)  :: null_nx,null_lsh(2),i
  REAL*8, intent(in) :: dt,dx,dr_dx,d2r_dx2
  REAL*8, dimension(null_lsh(1),null_lsh(2),null_nx), intent(in) :: aa_n,aa_p
  REAL*8, dimension(null_lsh(1),null_lsh(2)), intent(out) :: aa_00,aa_01,aa_04,aa_11,aa_14
  COMPLEX*16, dimension(null_lsh(1),null_lsh(2)), intent(out)  :: aa_02,aa_03,aa_12,aa_13,aa_22, &  
  aa_23,aa_24,aa_33,aa_34
  REAL*8, dimension(null_lsh(1),null_lsh(2)) :: aa_0x,aa_xx,aa_x4
 
  if (i.le.1.or.i.gt.null_nx) call CCTK_Warn(0,53,"NullConstr_Util.F90","NullConstr", "invalid index")

  aa_00=(aa_n(:,:,i)+aa_p(:,:,i))/2
  aa_04=(aa_n(:,:,i)-aa_p(:,:,i))/dt

  if (i.eq.null_nx) then 
     aa_0x=(aa_n(:,:,i)-aa_n(:,:,i-1)+aa_p(:,:,i)-aa_p(:,:,i-1))/(2*dx)
     aa_01=aa_0x/dr_dx
     aa_x4=(aa_n(:,:,i)-aa_n(:,:,i-1)-aa_p(:,:,i)+aa_p(:,:,i-1))/(dt*dx)
     aa_14=aa_x4/dr_dx
     aa_xx=(aa_n(:,:,(i-1)+1)-2*aa_n(:,:,(i-1))+aa_n(:,:,(i-1)-1) &
     +aa_p(:,:,(i-1)+1)-2*aa_p(:,:,(i-1))+aa_p(:,:,(i-1)-1))/(2*dx**2)
     aa_11=aa_xx/dr_dx**2-aa_0x*d2r_dx2/dr_dx**3
  else
     aa_0x=(aa_n(:,:,i+1)-aa_n(:,:,i-1)+aa_p(:,:,i+1)-aa_p(:,:,i-1))/(4*dx)
     aa_01=aa_0x/dr_dx
     aa_x4=(aa_n(:,:,i+1)-aa_n(:,:,i-1)-aa_p(:,:,i+1)+aa_p(:,:,i-1))/(2*dt*dx)
     aa_14=aa_x4/dr_dx
     aa_xx=(aa_n(:,:,i+1)-2*aa_n(:,:,i)+aa_n(:,:,i-1) &
     +aa_p(:,:,i+1)-2*aa_p(:,:,i)+aa_p(:,:,i-1))/(2*dx**2)
     aa_11=aa_xx/dr_dx**2-aa_0x*d2r_dx2/dr_dx**3
  end if

  call NullInterp_d1(aa_02, dcmplx(aa_00,0), 0_ik, 1_ik) 
  call NullInterp_d1(aa_03, dcmplx(aa_00,0), 0_ik, -1_ik) 
  call NullInterp_d1(aa_12, dcmplx(aa_01,0), 0_ik, 1_ik) 
  call NullInterp_d1(aa_13, dcmplx(aa_01,0), 0_ik, -1_ik) 
  call NullInterp_d1(aa_24, dcmplx(aa_04,0), 0_ik, 1_ik) 
  call NullInterp_d1(aa_34, dcmplx(aa_04,0), 0_ik, -1_ik) 
  call NullInterp_d1(aa_22, aa_02, 1_ik, 1_ik) 
  call NullInterp_d1(aa_33, aa_03, -1_ik, -1_ik) 
  call NullInterp_d1(aa_23, aa_03, -1_ik, 1_ik) 


end subroutine real_derivs

subroutine cmplx_derivs &
  (aa_n,aa_p,dt,dx,null_nx,null_lsh,i,spin_wt,dr_dx,d2r_dx2,half_grid, &
  aa_00,aa_01,aa_02,aa_03,aa_04,aa_11,aa_12,aa_13,aa_14, &
  aa_22,aa_23,aa_24,aa_33,aa_34, &
  aab_00,aab_01,aab_02,aab_03,aab_04,aab_11,aab_12,aab_13,aab_14, &
  aab_22,aab_23,aab_24,aab_33,aab_34)
  use NullInterp
  implicit none
  
  INTEGER*4, parameter :: izero = 0
  integer, parameter :: ik = kind(izero)

  INTEGER*4, intent(in)  :: null_nx,null_lsh(2),i,spin_wt
  LOGICAL, intent(in) :: half_grid
  REAL*8, intent(in) :: dt,dx,dr_dx,d2r_dx2

  COMPLEX*16, dimension(null_lsh(1),null_lsh(2),null_nx), intent(in) :: aa_n,aa_p
  COMPLEX*16, dimension(null_lsh(1),null_lsh(2)), intent(out) :: aa_00,aa_01,aa_04,aa_11,aa_14
  COMPLEX*16, dimension(null_lsh(1),null_lsh(2)), intent(out) :: aa_02,aa_03,aa_12,aa_13,aa_22, &
  aa_23,aa_24,aa_33,aa_34
  COMPLEX*16, dimension(null_lsh(1),null_lsh(2)), intent(out) :: aab_00,aab_01,aab_04,aab_11,aab_14
  COMPLEX*16, dimension(null_lsh(1),null_lsh(2)), intent(out) :: aab_02,aab_03,aab_12,aab_13,aab_22, &
  aab_23,aab_24,aab_33,aab_34
  COMPLEX*16, dimension(null_lsh(1),null_lsh(2)):: aa_0x,aa_x4,aa_xx

  if (i.le.1.or.i.gt.null_nx) call CCTK_Warn(0,122,"NullConstr_Util.F90","NullConstr", "invalid index")

  if(half_grid) then
     aa_00=(aa_n(:,:,i)+aa_p(:,:,i)+aa_n(:,:,i-1)+aa_p(:,:,i-1))/4
     aa_04=(aa_n(:,:,i)-aa_p(:,:,i)+aa_n(:,:,i-1)-aa_p(:,:,i-1))/(2*dt)
     aa_0x=(aa_n(:,:,i)-aa_n(:,:,i-1)+aa_p(:,:,i)-aa_p(:,:,i-1))/(2*dx)
     aa_x4=(aa_n(:,:,i)-aa_n(:,:,i-1)- &
     aa_p(:,:,i)+aa_p(:,:,i-1))/(dt*dx)
     if (i.eq.null_nx) then 
        aa_xx=(aa_n(:,:,i-1+1)-aa_n(:,:,i-1)-aa_n(:,:,i-1-1)+aa_n(:,:,i-1-2) &
        +aa_p(:,:,i-1+1)-aa_p(:,:,i-1)-aa_p(:,:,i-1-1)+aa_p(:,:,i-1-2) )/(4*dx**2)
     else
        aa_xx=(aa_n(:,:,i+1)-aa_n(:,:,i)-aa_n(:,:,i-1)+aa_n(:,:,i-2) &
        +aa_p(:,:,i+1)-aa_p(:,:,i)-aa_p(:,:,i-1)+aa_p(:,:,i-2) )/(4*dx**2)
     end if
  else
     aa_00=(aa_n(:,:,i)+aa_p(:,:,i))/2
     aa_04=(aa_n(:,:,i)-aa_p(:,:,i))/dt
     if (i.eq.null_nx) then  
        aa_0x=(aa_n(:,:,i)-aa_n(:,:,i-1)+aa_p(:,:,i)-aa_p(:,:,i-1))/(2*dx)
        aa_x4=(aa_n(:,:,i)-aa_n(:,:,i-1)- &
        aa_p(:,:,i)+aa_p(:,:,i-1))/(dt*dx)
        aa_xx=(aa_n(:,:,i-1+1)-2*aa_n(:,:,i-1)+aa_n(:,:,i-1-1) &
        +aa_p(:,:,i-1+1)-2*aa_p(:,:,i-1)+aa_p(:,:,i-1-1))/(2*dx**2)
     else
        aa_0x=(aa_n(:,:,i+1)-aa_n(:,:,i-1)+aa_p(:,:,i+1)-aa_p(:,:,i-1))/(4*dx)
        aa_x4=(aa_n(:,:,i+1)-aa_n(:,:,i-1)- &
        aa_p(:,:,i+1)+aa_p(:,:,i-1))/(2*dt*dx)
        aa_xx=(aa_n(:,:,i+1)-2*aa_n(:,:,i)+aa_n(:,:,i-1) &
        +aa_p(:,:,i+1)-2*aa_p(:,:,i)+aa_p(:,:,i-1))/(2*dx**2)
     end if
  end if
  aa_01=aa_0x/dr_dx
  aa_14=aa_x4/dr_dx
  aa_11=aa_xx/dr_dx**2-aa_0x*d2r_dx2/dr_dx**3
  call NullInterp_d1(aa_02, aa_00, spin_wt, 1_ik) 
  call NullInterp_d1(aa_03, aa_00, spin_wt, -1_ik) 
  call NullInterp_d1(aa_12, aa_01, spin_wt, 1_ik) 
  call NullInterp_d1(aa_13, aa_01, spin_wt, -1_ik) 
  call NullInterp_d1(aa_24, aa_04, spin_wt, 1_ik) 
  call NullInterp_d1(aa_34, aa_04, spin_wt, -1_ik) 
  call NullInterp_d1(aa_22, aa_02, spin_wt+1, 1_ik) 
  call NullInterp_d1(aa_33, aa_03, spin_wt-1, -1_ik) 
  call NullInterp_d1(aa_23, aa_03, spin_wt-1, 1_ik)
  
  aab_00=conjg(aa_00)
  aab_01=conjg(aa_01)
  aab_04=conjg(aa_04)
  aab_14=conjg(aa_14)
  aab_11=conjg(aa_11)
  aab_02=conjg(aa_03)
  aab_03=conjg(aa_02)
  aab_12=conjg(aa_13)
  aab_13=conjg(aa_12)
  aab_24=conjg(aa_34)
  aab_34=conjg(aa_24)
  aab_22=conjg(aa_33)
  aab_33=conjg(aa_22)
  call NullInterp_d1(aab_23, aab_03, -spin_wt-1, 1_ik)


end subroutine cmplx_derivs












end module NullConstr_Util
