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

      subroutine ivp_conformal_trans(CCTK_FARGUMENTS)
      implicit none
      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS
      INTEGER CCTK_Equals

      INTEGER, DIMENSION(3) :: sym
c     local vars
      integer i,j,k,bbox_save(6)
      CCTK_REAL divw,one,temp

      CCTK_REAL finf, npow, cw, velocity, dx, dy, dz
      integer interp_code, ierr
      integer octant, bitant, quadrant

      one = 1.0d0
      octant = (CCTK_Equals(domain,"octant"))
      bitant = (CCTK_Equals(domain,"bitant"))
      quadrant = (CCTK_Equals(domain,"quadrant"))

c     construct the A^TWIDLE (using CONFORMAL inverse metric!)
      kxx = 0.0d0
      kxy = 0.0d0
      kxz = 0.0d0
      kyy = 0.0d0
      kyz = 0.0d0
      kzz = 0.0d0

c     dx settings
      dx = cctk_delta_space(1)
      dy = cctk_delta_space(2)
      dz = cctk_delta_space(3)

c changing stencil_size to cctk_dim.
c Malcolm: cctk_dim is the dimensionality of the problem
c Is there anything wrong with just setting this to 1 for
c the reconstruction?

      do k = 2, cctk_lsh(3)-1
      do j = 2, cctk_lsh(2)-1
      do i = 2, cctk_lsh(1)-1
         divw = ( (sqrt(ivp_det(i+1,j,k))*ivp_wx(i+1,j,k) -
     .             sqrt(ivp_det(i-1,j,k))*ivp_wx(i-1,j,k))/(2.0d0*dx) +
     .            (sqrt(ivp_det(i,j+1,k))*ivp_wy(i,j+1,k) -
     .             sqrt(ivp_det(i,j-1,k))*ivp_wy(i,j-1,k))/(2.0d0*dy) +
     .            (sqrt(ivp_det(i,j,k+1))*ivp_wz(i,j,k+1) -
     .             sqrt(ivp_det(i,j,k-1))*ivp_wz(i,j,k-1))/(2.0d0*dz))/
     .          sqrt(ivp_det(i,j,k))
         kxx(i,j,k) = -
     .      ivp_wx(i,j,k) * (ivp_uxx(i+1,j,k)-ivp_uxx(i-1,j,k))/(2.0d0*dx) -
     .      ivp_wy(i,j,k) * (ivp_uxx(i,j+1,k)-ivp_uxx(i,j-1,k))/(2.0d0*dy) -
     .      ivp_wz(i,j,k) * (ivp_uxx(i,j,k+1)-ivp_uxx(i,j,k-1))/(2.0d0*dz) +
     .      ivp_uxx(i,j,k)* (ivp_wx(i+1,j,k) -ivp_wx(i-1,j,k ))/(2.0d0*dx) +
     .      ivp_uxy(i,j,k)* (ivp_wx(i,j+1,k) -ivp_wx(i,j-1,k ))/(2.0d0*dy) +
     .      ivp_uxz(i,j,k)* (ivp_wx(i,j,k+1) -ivp_wx(i,j,k-1 ))/(2.0d0*dz) +
     .      ivp_uxx(i,j,k)* (ivp_wx(i+1,j,k) -ivp_wx(i-1,j,k ))/(2.0d0*dx) +
     .      ivp_uxy(i,j,k)* (ivp_wx(i,j+1,k) -ivp_wx(i,j-1,k ))/(2.0d0*dy) +
     .      ivp_uxz(i,j,k)* (ivp_wx(i,j,k+1) -ivp_wx(i,j,k-1 ))/(2.0d0*dz) -
     .      (2.0d0/3.0d0) * divw * ivp_uxx(i,j,k)
         kyy(i,j,k) = -
     .      ivp_wx(i,j,k) * (ivp_uyy(i+1,j,k)-ivp_uyy(i-1,j,k))/(2.0d0*dx) -
     .      ivp_wy(i,j,k) * (ivp_uyy(i,j+1,k)-ivp_uyy(i,j-1,k))/(2.0d0*dy) -
     .      ivp_wz(i,j,k) * (ivp_uyy(i,j,k+1)-ivp_uyy(i,j,k-1))/(2.0d0*dz) +
     .      ivp_uxy(i,j,k)* (ivp_wy(i+1,j,k) -ivp_wy(i-1,j,k ))/(2.0d0*dx) +
     .      ivp_uyy(i,j,k)* (ivp_wy(i,j+1,k) -ivp_wy(i,j-1,k ))/(2.0d0*dy) +
     .      ivp_uyz(i,j,k)* (ivp_wy(i,j,k+1) -ivp_wy(i,j,k-1 ))/(2.0d0*dz) +
     .      ivp_uxy(i,j,k)* (ivp_wy(i+1,j,k) -ivp_wy(i-1,j,k ))/(2.0d0*dx) +
     .      ivp_uyy(i,j,k)* (ivp_wy(i,j+1,k) -ivp_wy(i,j-1,k ))/(2.0d0*dy) +
     .      ivp_uyz(i,j,k)* (ivp_wy(i,j,k+1) -ivp_wy(i,j,k-1 ))/(2.0d0*dz) -
     .      (2.0d0/3.0d0) * divw * ivp_uyy(i,j,k)
         kzz(i,j,k) = -
     .      ivp_wx(i,j,k) * (ivp_uzz(i+1,j,k)-ivp_uzz(i-1,j,k))/(2.0d0*dx) -
     .      ivp_wy(i,j,k) * (ivp_uzz(i,j+1,k)-ivp_uzz(i,j-1,k))/(2.0d0*dy) -
     .      ivp_wz(i,j,k) * (ivp_uzz(i,j,k+1)-ivp_uzz(i,j,k-1))/(2.0d0*dz) +
     .      ivp_uxz(i,j,k)* (ivp_wz(i+1,j,k) -ivp_wz(i-1,j,k ))/(2.0d0*dx) +
     .      ivp_uyz(i,j,k)* (ivp_wz(i,j+1,k) -ivp_wz(i,j-1,k ))/(2.0d0*dy) +
     .      ivp_uzz(i,j,k)* (ivp_wz(i,j,k+1) -ivp_wz(i,j,k-1 ))/(2.0d0*dz) +
     .      ivp_uxz(i,j,k)* (ivp_wz(i+1,j,k) -ivp_wz(i-1,j,k ))/(2.0d0*dx) +
     .      ivp_uyz(i,j,k)* (ivp_wz(i,j+1,k) -ivp_wz(i,j-1,k ))/(2.0d0*dy) +
     .      ivp_uzz(i,j,k)* (ivp_wz(i,j,k+1) -ivp_wz(i,j,k-1 ))/(2.0d0*dz) -
     .      (2.0d0/3.0d0) * divw * ivp_uzz(i,j,k)
         kxy(i,j,k) = -
     .      ivp_wx(i,j,k) * (ivp_uxy(i+1,j,k)-ivp_uxy(i-1,j,k))/(2.0d0*dx) -
     .      ivp_wy(i,j,k) * (ivp_uxy(i,j+1,k)-ivp_uxy(i,j-1,k))/(2.0d0*dy) -
     .      ivp_wz(i,j,k) * (ivp_uxy(i,j,k+1)-ivp_uxy(i,j,k-1))/(2.0d0*dz) +
     .      ivp_uxx(i,j,k)* (ivp_wy(i+1,j,k) -ivp_wy(i-1,j,k ))/(2.0d0*dx) +
     .      ivp_uxy(i,j,k)* (ivp_wy(i,j+1,k) -ivp_wy(i,j-1,k ))/(2.0d0*dy) +
     .      ivp_uxz(i,j,k)* (ivp_wy(i,j,k+1) -ivp_wy(i,j,k-1 ))/(2.0d0*dz) +
     .      ivp_uxy(i,j,k)* (ivp_wx(i+1,j,k) -ivp_wx(i-1,j,k ))/(2.0d0*dx) +
     .      ivp_uyy(i,j,k)* (ivp_wx(i,j+1,k) -ivp_wx(i,j-1,k ))/(2.0d0*dy) +
     .      ivp_uyz(i,j,k)* (ivp_wx(i,j,k+1) -ivp_wx(i,j,k-1 ))/(2.0d0*dz) -
     .      (2.0d0/3.0d0) * divw * ivp_uxy(i,j,k)
         kxz(i,j,k) = -
     .      ivp_wx(i,j,k) * (ivp_uxz(i+1,j,k)-ivp_uxz(i-1,j,k))/(2.0d0*dx) -
     .      ivp_wy(i,j,k) * (ivp_uxz(i,j+1,k)-ivp_uxz(i,j-1,k))/(2.0d0*dy) -
     .      ivp_wz(i,j,k) * (ivp_uxz(i,j,k+1)-ivp_uxz(i,j,k-1))/(2.0d0*dz) +
     .      ivp_uxx(i,j,k)* (ivp_wz(i+1,j,k) -ivp_wz(i-1,j,k ))/(2.0d0*dx) +
     .      ivp_uxy(i,j,k)* (ivp_wz(i,j+1,k) -ivp_wz(i,j-1,k ))/(2.0d0*dy) +
     .      ivp_uxz(i,j,k)* (ivp_wz(i,j,k+1) -ivp_wz(i,j,k-1 ))/(2.0d0*dz) +
     .      ivp_uxz(i,j,k)* (ivp_wx(i+1,j,k) -ivp_wx(i-1,j,k ))/(2.0d0*dx) +
     .      ivp_uyz(i,j,k)* (ivp_wx(i,j+1,k) -ivp_wx(i,j-1,k ))/(2.0d0*dy) +
     .      ivp_uzz(i,j,k)* (ivp_wx(i,j,k+1) -ivp_wx(i,j,k-1 ))/(2.0d0*dz) -
     .      (2.0d0/3.0d0) * divw * ivp_uxz(i,j,k)
         kyz(i,j,k) = -
     .      ivp_wx(i,j,k) * (ivp_uyz(i+1,j,k)-ivp_uyz(i-1,j,k))/(2.0d0*dx) -
     .      ivp_wy(i,j,k) * (ivp_uyz(i,j+1,k)-ivp_uyz(i,j-1,k))/(2.0d0*dy) -
     .      ivp_wz(i,j,k) * (ivp_uyz(i,j,k+1)-ivp_uyz(i,j,k-1))/(2.0d0*dz) +
     .      ivp_uxy(i,j,k)* (ivp_wz(i+1,j,k) -ivp_wz(i-1,j,k ))/(2.0d0*dx) +
     .      ivp_uyy(i,j,k)* (ivp_wz(i,j+1,k) -ivp_wz(i,j-1,k ))/(2.0d0*dy) +
     .      ivp_uyz(i,j,k)* (ivp_wz(i,j,k+1) -ivp_wz(i,j,k-1 ))/(2.0d0*dz) +
     .      ivp_uxz(i,j,k)* (ivp_wy(i+1,j,k) -ivp_wy(i-1,j,k ))/(2.0d0*dx) +
     .      ivp_uyz(i,j,k)* (ivp_wy(i,j+1,k) -ivp_wy(i,j-1,k ))/(2.0d0*dy) +
     .      ivp_uzz(i,j,k)* (ivp_wy(i,j,k+1) -ivp_wy(i,j,k-1 ))/(2.0d0*dz) -
     .      (2.0d0/3.0d0) * divw * ivp_uyz(i,j,k)
c        now, add in astar to get A^twidle
         kxx(i,j,k) = kxx(i,j,k) + ivp_astarxx(i,j,k)
         kxy(i,j,k) = kxy(i,j,k) + ivp_astarxy(i,j,k)
         kxz(i,j,k) = kxz(i,j,k) + ivp_astarxz(i,j,k)
         kyy(i,j,k) = kyy(i,j,k) + ivp_astaryy(i,j,k)
         kyz(i,j,k) = kyz(i,j,k) + ivp_astaryz(i,j,k)
         kzz(i,j,k) = kzz(i,j,k) + ivp_astarzz(i,j,k)
      enddo
      enddo
      enddo

c     Need to apply some boundary condition on the re-constructed Kij

      if(CCTK_Equals(ivp_conf_trans_bound,"static").ne.0) then
         continue
      elseif(CCTK_Equals(ivp_conf_trans_bound,"robin").ne.0) then
#ifdef MAHC
         write(*,*)'hey, thorn MAHC finally converted! Great, now fix the'
         write(*,*)'segments in IVP in order to use MAHC.'
c$$$         finf = 0.0d0
c$$$         velocity = 1.0d0
c$$$         npow = 1.0d0
c$$$         cw = 1.0d0
c$$$         interp_code = 2
c$$$          call mahc_rad_onebound(ARENA,kxx,kxx,finf,velocity,npow,cw,dt,
c$$$     .      one,one,one,r,x,y,z,interp_code)        
c$$$          call mahc_rad_onebound(ARENA,kxy,kxy,finf,velocity,npow,cw,dt,
c$$$     .      one,one,one,r,x,y,z,interp_code)      
c$$$          call mahc_rad_onebound(ARENA,kxz,kxz,finf,velocity,npow,cw,dt,
c$$$     .      one,one,one,r,x,y,z,interp_code)        
c$$$          call mahc_rad_onebound(ARENA,kyy,kyy,finf,velocity,npow,cw,dt,
c$$$     .      one,one,one,r,x,y,z,interp_code)        
c$$$          call mahc_rad_onebound(ARENA,kyz,kyz,finf,velocity,npow,cw,dt,
c$$$     .      one,one,one,r,x,y,z,interp_code)        
c$$$          call mahc_rad_onebound(ARENA,kzz,kzz,finf,velocity,npow,cw,dt,
c$$$     .      one,one,one,r,x,y,z,interp_code)        
#else
          write(*,*)'Sorry: Robin b.c. after the conf. transformation'
          write(*,*)'requires thorn MAHC.  Either recompile or add'
          write(*,*)'these conditions by hand.'
          write(*,*)'or set ivp_conf_trans_bound = static'
          STOP
#endif
      endif

      call CCTK_SyncGroup(cctkGH, "einstein::curv")

c     should apply inner boundary conditions here
      bbox_save = cctk_bbox
      sym(1) = 1
      sym(2) = 1
      sym(3) = 1
      if(octant.eq.1) then
         cctk_bbox(2) = 0
         cctk_bbox(4) = 0
         cctk_bbox(6) = 0
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxx')
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kyy')
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kzz')
         sym(1) = -1
         sym(2) = -1
         sym(3) = 1
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxy')
         sym(1) = -1
         sym(2) = 1
         sym(3) = -1
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxz')
         sym(1) = 1
         sym(2) = -1
         sym(3) = -1
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kyz')
      endif
      if(quadrant.eq.1)then
         cctk_bbox(2) = 0
         cctk_bbox(4) = 0
         cctk_bbox(5) = 0
         cctk_bbox(6) = 0
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxx')
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kyy')
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kzz')
         sym(1) = -1
         sym(2) = -1
         sym(3) = 1
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxy')
         sym(1) = -1
         sym(2) = 1
         sym(3) = -1
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxz')
         sym(1) = 1
         sym(2) = -1
         sym(3) = -1
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kyz')
      endif
      if(bitant.eq.1)then
         cctk_bbox(1) = 0
         cctk_bbox(2) = 0
         cctk_bbox(3) = 0
         cctk_bbox(4) = 0
         cctk_bbox(6) = 0
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxx')
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kyy')
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kzz')
         sym(1) = -1
         sym(2) = -1
         sym(3) = 1
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxy')
         sym(1) = -1
         sym(2) = 1
         sym(3) = -1
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxz')
         sym(1) = 1
         sym(2) = -1
         sym(3) = -1
         call SetCartSymVN(ierr,cctkGH,sym,'einstein::kyz')
      endif
      cctk_bbox = bbox_save

c     get rid of any numerically induced tr(A).  Watch out!  We are
c     using conformal metric to subtract out trace!
      if(enforce_traceless_Atwiddle.eq.1) then
         do k=1,cctk_lsh(3)
         do j=1,cctk_lsh(1)
         do i=1,cctk_lsh(1)

c           calculate trace(Atwiddle), store in temp
            temp = (kxx(i,j,k)*gxx(i,j,k) + kyy(i,j,k)*gyy(i,j,k) +
     .              kzz(i,j,k)*gzz(i,j,k) + 2.0d0*kxy(i,j,k)*gxy(i,j,k) +
     .              2.0d0*kxz(i,j,k)*gxz(i,j,k) + 
     .              2.0d0*kyz(i,j,k)*gyz(i,j,k))
            kxx(i,j,k) = kxx(i,j,k) - ivp_uxx(i,j,k)*temp/3.0d0
            kxy(i,j,k) = kxy(i,j,k) - ivp_uxy(i,j,k)*temp/3.0d0
            kxz(i,j,k) = kxz(i,j,k) - ivp_uxz(i,j,k)*temp/3.0d0
            kyy(i,j,k) = kyy(i,j,k) - ivp_uyy(i,j,k)*temp/3.0d0
            kyz(i,j,k) = kyz(i,j,k) - ivp_uyz(i,j,k)*temp/3.0d0
            kzz(i,j,k) = kzz(i,j,k) - ivp_uzz(i,j,k)*temp/3.0d0

         enddo
         enddo
         enddo
      endif

c     apply conformal factor to conformal metric to get physical metric:
      gxx = (ivp_psi**4) * gxx
      gxy = (ivp_psi**4) * gxy
      gxz = (ivp_psi**4) * gxz
      gyy = (ivp_psi**4) * gyy
      gyz = (ivp_psi**4) * gyz
      gzz = (ivp_psi**4) * gzz

c     construct upper physical extrinsic curvature:
      kxx = (ivp_psi**(-10))*kxx + 
     .   (ivp_psi**(-4)) * ivp_uxx * trK / 3.0d0
      kxy = (ivp_psi**(-10))*kxy + 
     .   (ivp_psi**(-4)) * ivp_uxy * trK / 3.0d0
      kxz = (ivp_psi**(-10))*kxz + 
     .   (ivp_psi**(-4)) * ivp_uxz * trK / 3.0d0
      kyy = (ivp_psi**(-10))*kyy + 
     .   (ivp_psi**(-4)) * ivp_uyy * trK / 3.0d0
      kyz = (ivp_psi**(-10))*kyz + 
     .   (ivp_psi**(-4)) * ivp_uyz * trK / 3.0d0
      kzz = (ivp_psi**(-10))*kzz + 
     .   (ivp_psi**(-4)) * ivp_uzz * trK / 3.0d0

c     evolution expects the lowered physical extrinsic curvature:
c     (use astar as temps)
      ivp_astarxx = 
     &  gxx**2*kxx+2.d0*gxx*gxy*kxy+2.d0*gxx*gxz*kxz+gxy**2*kyy+2.d0*gxy
     &  *gxz*kyz+gxz**2*kzz
      ivp_astarxy =
     &  gxx*gxy*kxx+gxy**2*kxy+gxx*gyy*kxy+gxy*gxz*kxz+gxx*gyz*kxz+gxy*g
     &  yy*kyy+gxz*gyy*kyz+gxy*gyz*kyz+gxz*gyz*kzz
      ivp_astarxz =
     &  gxx*gxz*kxx+gxy*gxz*kxy+gxx*gyz*kxy+gxz**2*kxz+gxx*gzz*kxz+gxy*g
     &  yz*kyy+gxz*gyz*kyz+gxy*gzz*kyz+gxz*gzz*kzz
      ivp_astaryy =
     &  gxy**2*kxx+2.d0*gxy*gyy*kxy+2.d0*gxy*gyz*kxz+gyy**2*kyy+2.d0*gyy
     &  *gyz*kyz+gyz**2*kzz
      ivp_astaryz =
     &  gxy*gxz*kxx+gxz*gyy*kxy+gxy*gyz*kxy+gxz*gyz*kxz+gxy*gzz*kxz+gyy*
     &  gyz*kyy+gyz**2*kyz+gyy*gzz*kyz+gyz*gzz*kzz
      ivp_astarzz =
     &  gxz**2*kxx+2.d0*gxz*gyz*kxy+2.d0*gxz*gzz*kxz+gyz**2*kyy+2.d0*gyz
     &  *gzz*kyz+gzz**2*kzz
      kxx = ivp_astarxx
      kxy = ivp_astarxy
      kxz = ivp_astarxz
      kyy = ivp_astaryy
      kyz = ivp_astaryz
      kzz = ivp_astarzz



      return
      end