#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "cctk_Functions.h"



subroutine SWTNS_outerboundary (CCTK_ARGUMENTS)
 implicit none
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_FUNCTIONS
  DECLARE_CCTK_PARAMETERS  

   integer :: i, j, k
!write(*,*) "Define Xeta"
!$OMP PARALLEL DO private(i,j,k)
  do k = 1, cctk_lsh(3)
     do j = 1, cctk_lsh(2)
        do i = 1, cctk_lsh(1)
		!grid_r=r(i,j,k)
		if(grid_r(i,j,k) == 0) then
		Xeta(i,j,k) = 0
		Xi(i,j,k) = 0 !inner boundary!		  
		else
		Xeta(i,j,k) = Xi(i,j,k) * exp(PHI(i,j,k))/(r(i,j,k)**2)
		end if
        end do
     end do
  end do
!$OMP END PARALLEL DO 

!updating TOV surface in PostStep

outerloop: do k=1,cctk_lsh(3)
 	do j=1,cctk_lsh(2)
  	      do i=1,cctk_lsh(1)
		!grid_r=r(i,j,k)
		 if (grid_r(i,j,k) > (TOV_surface -rprec) .AND. grid_r(i,j,k) < (TOV_surface + rprec) ) then
		 TOV_surface = TOV_surface + Xeta(i,j,k)
		 !write(*,*) "i = ", i
!		 write(*,*) "j ", j
!		 write(*,*) "k = ", k
!		 write(*,*) "r = ", r(i,j,K)
!		 write(*,*) "Xi(i,j,k) = ", Xi(i,j,k)	
!                write(*,*) "TOV_surface = ", TOV_surface
		exit outerloop
		end if
	      end do
	end do 
  end do outerloop

end subroutine SWTNS_outerboundary


subroutine SWTNS_RHS_outerboundary (CCTK_ARGUMENTS)
  implicit none
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_FUNCTIONS
  DECLARE_CCTK_PARAMETERS

  integer :: i, j, k



!$OMP PARALLEL DO private(i,j,k)
  do k = 1, cctk_lsh(3)
     do j = 1, cctk_lsh(2)
        do i = 1, cctk_lsh(1)
		  !grid_r=r(i,j,k)
		  !if(grid_r(i,j,k) < (TOV_surface-rprec)) then
		  !Xidot(i,j,k) = Pi(i,j,k)
		  !Pidot(i,j,k) = Wrinv(i,j,k) *(drP(i,j,k) * drXi(i,j,k)  + P(i,j,k) * drrXi(i,j,k) + Qr(i,j,k)* Xi(i,j,k))
		  if (grid_r(i,j,k) >= (TOV_surface-rprec) .AND. grid_r(i,j,k) <= (TOV_surface+ rprec)) then
		  drXi(i,j,k) = 0
		  Pidot(i,j,k) = Wrinv(i,j,k) * P(i,j,k) * drrXi(i,j,k) + Qr(i,j,k) * Wrinv(i,j,k) * Xi(i,j,k) !drXi=0
		  !write(*,*) "r = " , r(i,j,k)
                  else if(grid_r(i,j,k) > (TOV_surface + rprec)) then
		  drXi(i,j,k) = 0
		  drrXi(i,j,k) = 0
		  Xi(i,j,k) = 0
 		  Pi(i,j,k) = 0
		  Xidot(i,j,k) = 0
		  Pidot(i,j,k) = 0
		  Xeta(i,j,k) = 0
		  end if
        end do
     end do
  end do
!$OMP END PARALLEL DO


end subroutine SWTNS_RHS_outerboundary

subroutine SWTNS_boundaries (CCTK_ARGUMENTS)
  implicit none
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_FUNCTIONS
  DECLARE_CCTK_PARAMETERS  
  
  character :: fbound*1000
  integer   :: fboundlen
  integer   :: ierr
  !write(*,*) "SWTNS_boundaries"
   call CCTK_FortranString (fboundlen, bound, fbound)
   if (fboundlen > len(fbound)) call CCTK_WARN (0, "internal error")
   
    ierr = Boundary_SelectGroupForBC &
         (cctkGH, CCTK_ALL_FACES, +1, -1, "SWTNS::scalar", fbound)
    if (ierr/=0) call CCTK_WARN (0, "internal error")
 
    ierr = Boundary_SelectGroupForBC &
        (cctkGH, CCTK_ALL_FACES, +1, -1, "SWTNS::density", fbound)
   if (ierr/=0) call CCTK_WARN (0, "internal error")
   
end subroutine SWTNS_boundaries