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

subroutine SWTNS_calc_rhs (CCTK_ARGUMENTS)
  implicit none
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_FUNCTIONS
  DECLARE_CCTK_PARAMETERS
  
   CCTK_INT, parameter :: izero = 0
   CCTK_INT :: i, j, k
   CCTK_REAL :: pio
    
   integer, parameter :: ik = kind (izero)
   
   integer :: na, nb, nc
   
   CCTK_REAL, dimension(:,:,:), allocatable &
        :: TOVold, Babo, drBabo, dxBabo, dyBabo, dzBabo, dxXi, dyXi, dzXi, &
            dxxXi, dyyXi, dzzXi, dxyXi, dxzXi, dyzXi

   na = cctk_lsh(1); nb = cctk_lsh(2); nc = cctk_lsh(3)
 !   write(*,*)"allocate RHS..."
  
    pio=3.14159265359
    !write(*,*) "Surface = " , TOV_surface
 allocate( Babo(na, nb, nc), &
	   drBabo(na, nb, nc), &
	   dxBabo(na, nb, nc), &
	   dyBabo(na, nb, nc), &
	   dzBabo(na, nb, nc), &
	   dxXi(na, nb, nc), &
           dyXi(na, nb, nc), &
           dzXi(na, nb, nc), &
           dxxXi(na, nb, nc), &
           dxyXi(na, nb, nc), &
           dxzXi(na, nb, nc), &
           dyyXi(na, nb, nc), &
           dyzXi(na, nb, nc), &
           dzzXi(na, nb, nc), &
	   TOVold(na,nb,nc)  )
 

!Ableitung drXi
           call globalDiff_gv (cctkGH, 0_ik, Xi, dxXi, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33,-1_ik)
           call globalDiff_gv (cctkGH, 1_ik, Xi, dyXi, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
           call globalDiff_gv (cctkGH, 2_ik, Xi, dzXi, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
             
           !drXi = dxXi*dxdr + dyXi*dydr + dzXi*dzdr
	   drXi = dxXi*dxdr + dyXi*dydr + dzXi*dzdr


           !Babo = P *drXi  
!2te Ableitung drrXi
!aus erster....
	   call globalDiff_gv (cctkGH, 0_ik, drXi, dxxXi, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33,-1_ik)
           call globalDiff_gv (cctkGH, 1_ik, drXi, dyyXi, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
           call globalDiff_gv (cctkGH, 2_ik, drXi, dzzXi, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
             
           !drXi = dxXi*dxdr + dyXi*dydr + dzXi*dzdr
	   drrXi = dxxXi*dxdr + dyyXi*dydr + dzzXi*dzdr


          
!!! Wave Eq.
       	Xidot = Pi
        Pidot = Wrinv *(drP * drXi  + P * drrXi + Qr* Xi)	
	!drBabo (drP * drXi  + P * drrXi)	

!        Pidot = 1/Wr * (drPr * drXi + Pr * drrXi) + Qr/Wr *Xi | original



  deallocate ( Babo, drBabo, dxBabo, dyBabo, dzBabo, dxXi, dyXi, dzXi, &
              dxxXi, dxyXi, dxzXi, dyyXi, dyzXi, dzzXi, TOVold )
  
end subroutine SWTNS_calc_rhs