














































c







      






































































 


































































































































































 






























































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                




















































































      
c

      subroutine DistortedBHIVP(cctk_dim,cctk_gsh,cctk_lsh,cctk_lbnd,cct
     &k_ubnd,cctk_tile_min,cctk_tile_max,cctk_ash,cctk_alignment,cctk_al
     &ignment_offset,cctk_from,cctk_to,cctk_bbox,cctk_delta_time,cctk_ti
     &me,cctk_delta_space,cctk_origin_space,cctk_levfac,cctk_levoff,cctk
     &_levoffdenom,cctk_timefac,cctk_convlevel,cctk_convfac,cctk_nghostz
     &ones,cctk_iteration,cctkGH,cctk_ash1,cctk_ash2,cctk_ash3, abseta,d
     &etaetapsi3d,detaphipsi3d,detapsi3d,detaqpsi3d,dphiphipsi3d,dphipsi
     &3d,dqphipsi3d,dqpsi3d,dqqpsi3d,eta,phi,psi3d,q,sign_eta, alp,alp_p
     &,alp_p_p,betax,betax_p,betax_p_p,betay,betay_p,betay_p_p,betaz,bet
     &az_p,betaz_p_p,coarse_dx,coarse_dy,coarse_dz,conformal_state,dtalp
     &,dtalp_p,dtalp_p_p,dtbetax,dtbetax_p,dtbetax_p_p,dtbetay,dtbetay_p
     &,dtbetay_p_p,dtbetaz,dtbetaz_p,dtbetaz_p_p,dtlapse_state,dtshift_s
     &tate,gxx,gxx_p,gxx_p_p,gxy,gxy_p,gxy_p_p,gxz,gxz_p,gxz_p_p,gyy,gyy
     &_p,gyy_p_p,gyz,gyz_p,gyz_p_p,gzz,gzz_p,gzz_p_p,kxx,kxx_p,kxx_p_p,k
     &xy,kxy_p,kxy_p_p,kxz,kxz_p,kxz_p_p,kyy,kyy_p,kyy_p_p,kyz,kyz_p,kyz
     &_p_p,kzz,kzz_p,kzz_p_p,psi,psix,psixx,psixy,psixz,psiy,psiyy,psiyz
     &,psiz,psizz,r,shift_state,x,y,z)

      implicit none


      INTEGER cctk_dim
      integer, parameter :: cctki_use_cctk_dim = kind(cctk_dim)
      INTEGER cctk_gsh (cctk_dim)
      integer, parameter :: cctki_use_cctk_gsh = kind(cctk_gsh)
      INTEGER cctk_lsh (cctk_dim)
      integer, parameter :: cctki_use_cctk_lsh = kind(cctk_lsh)
      INTEGER cctk_lbnd (cctk_dim)
      integer, parameter :: cctki_use_cctk_lbnd = kind(cctk_lbnd)
      INTEGER cctk_ubnd (cctk_dim)
      integer, parameter :: cctki_use_cctk_ubnd = kind(cctk_ubnd)
      INTEGER cctk_tile_min (cctk_dim)
      integer, parameter :: cctki_use_cctk_tile_min = kind(cctk_tile_min
     &)
      INTEGER cctk_tile_max (cctk_dim)
      integer, parameter :: cctki_use_cctk_tile_max = kind(cctk_tile_max
     &)
      INTEGER cctk_ash (cctk_dim)
      integer, parameter :: cctki_use_cctk_ash = kind(cctk_ash)
      INTEGER cctk_alignment
      integer, parameter :: cctki_use_cctk_alignment = kind(cctk_alignme
     &nt)
      INTEGER cctk_alignment_offset
      integer, parameter :: cctki_use_cctk_alignment_offset = kind(cctk_
     &alignment_offset)
      INTEGER cctk_from (cctk_dim)
      integer, parameter :: cctki_use_cctk_from = kind(cctk_from)
      INTEGER cctk_to (cctk_dim)
      integer, parameter :: cctki_use_cctk_to = kind(cctk_to)
      INTEGER cctk_bbox (2*cctk_dim)
      integer, parameter :: cctki_use_cctk_bbox = kind(cctk_bbox)
      REAL*8 cctk_delta_time
      integer, parameter :: cctki_use_cctk_delta_time = kind(cctk_delta_
     &time)
      REAL*8 cctk_time
      integer, parameter :: cctki_use_cctk_time = kind(cctk_time)
      REAL*8 cctk_delta_space (cctk_dim)
      integer, parameter :: cctki_use_cctk_delta_space = kind(cctk_delta
     &_space)
      REAL*8 cctk_origin_space (cctk_dim)
      integer, parameter :: cctki_use_cctk_origin_space = kind(cctk_orig
     &in_space)
      INTEGER cctk_levfac (cctk_dim)
      integer, parameter :: cctki_use_cctk_levfac = kind(cctk_levfac)
      INTEGER cctk_levoff (cctk_dim)
      integer, parameter :: cctki_use_cctk_levoff = kind(cctk_levoff)
      INTEGER cctk_levoffdenom (cctk_dim)
      integer, parameter :: cctki_use_cctk_levoffdenom = kind(cctk_levof
     &fdenom)
      INTEGER cctk_timefac
      integer, parameter :: cctki_use_cctk_timefac = kind(cctk_timefac)
      INTEGER cctk_convlevel
      integer, parameter :: cctki_use_cctk_convlevel = kind(cctk_convlev
     &el)
      INTEGER cctk_convfac
      integer, parameter :: cctki_use_cctk_convfac = kind(cctk_convfac)
      INTEGER cctk_nghostzones (cctk_dim)
      integer, parameter :: cctki_use_cctk_nghostzones = kind(cctk_nghos
     &tzones)
      INTEGER cctk_iteration
      integer, parameter :: cctki_use_cctk_iteration = kind(cctk_iterati
     &on)
      integer*8 cctkGH
      integer, parameter :: cctki_use_cctkGH = kind(cctkGH)
      INTEGER cctk_ash1
      integer, parameter :: cctki_use_cctk_ash1 = kind(cctk_ash1)
      INTEGER cctk_ash2
      integer, parameter :: cctki_use_cctk_ash2 = kind(cctk_ash2)
      INTEGER cctk_ash3
      integer, parameter :: cctki_use_cctk_ash3 = kind(cctk_ash3)
      REAL*8 abseta (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_abseta = kind(abseta)
      REAL*8 detaetapsi3d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_detaetapsi3d = kind(detaetapsi3d)
      REAL*8 detaphipsi3d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_detaphipsi3d = kind(detaphipsi3d)
      REAL*8 detapsi3d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_detapsi3d = kind(detapsi3d)
      REAL*8 detaqpsi3d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_detaqpsi3d = kind(detaqpsi3d)
      REAL*8 dphiphipsi3d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dphiphipsi3d = kind(dphiphipsi3d)
      REAL*8 dphipsi3d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dphipsi3d = kind(dphipsi3d)
      REAL*8 dqphipsi3d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dqphipsi3d = kind(dqphipsi3d)
      REAL*8 dqpsi3d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dqpsi3d = kind(dqpsi3d)
      REAL*8 dqqpsi3d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dqqpsi3d = kind(dqqpsi3d)
      REAL*8 eta (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_eta = kind(eta)
      REAL*8 phi (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_phi = kind(phi)
      REAL*8 psi3d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psi3d = kind(psi3d)
      REAL*8 q (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_q = kind(q)
      REAL*8 sign_eta (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_sign_eta = kind(sign_eta)
      REAL*8 alp (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_alp = kind(alp)
      REAL*8 alp_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_alp_p = kind(alp_p)
      REAL*8 alp_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_alp_p_p = kind(alp_p_p)
      REAL*8 betax (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_betax = kind(betax)
      REAL*8 betax_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_betax_p = kind(betax_p)
      REAL*8 betax_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_betax_p_p = kind(betax_p_p)
      REAL*8 betay (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_betay = kind(betay)
      REAL*8 betay_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_betay_p = kind(betay_p)
      REAL*8 betay_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_betay_p_p = kind(betay_p_p)
      REAL*8 betaz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_betaz = kind(betaz)
      REAL*8 betaz_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_betaz_p = kind(betaz_p)
      REAL*8 betaz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_betaz_p_p = kind(betaz_p_p)
      REAL*8 coarse_dx
      integer, parameter :: cctki_use_coarse_dx = kind(coarse_dx)
      REAL*8 coarse_dy
      integer, parameter :: cctki_use_coarse_dy = kind(coarse_dy)
      REAL*8 coarse_dz
      integer, parameter :: cctki_use_coarse_dz = kind(coarse_dz)
      INTEGER*4 conformal_state
      integer, parameter :: cctki_use_conformal_state = kind(conformal_s
     &tate)
      REAL*8 dtalp (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dtalp = kind(dtalp)
      REAL*8 dtalp_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dtalp_p = kind(dtalp_p)
      REAL*8 dtalp_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dtalp_p_p = kind(dtalp_p_p)
      REAL*8 dtbetax (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dtbetax = kind(dtbetax)
      REAL*8 dtbetax_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dtbetax_p = kind(dtbetax_p)
      REAL*8 dtbetax_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dtbetax_p_p = kind(dtbetax_p_p)
      REAL*8 dtbetay (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dtbetay = kind(dtbetay)
      REAL*8 dtbetay_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dtbetay_p = kind(dtbetay_p)
      REAL*8 dtbetay_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dtbetay_p_p = kind(dtbetay_p_p)
      REAL*8 dtbetaz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dtbetaz = kind(dtbetaz)
      REAL*8 dtbetaz_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dtbetaz_p = kind(dtbetaz_p)
      REAL*8 dtbetaz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dtbetaz_p_p = kind(dtbetaz_p_p)
      INTEGER*4 dtlapse_state
      integer, parameter :: cctki_use_dtlapse_state = kind(dtlapse_state
     &)
      INTEGER*4 dtshift_state
      integer, parameter :: cctki_use_dtshift_state = kind(dtshift_state
     &)
      REAL*8 gxx (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gxx = kind(gxx)
      REAL*8 gxx_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gxx_p = kind(gxx_p)
      REAL*8 gxx_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gxx_p_p = kind(gxx_p_p)
      REAL*8 gxy (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gxy = kind(gxy)
      REAL*8 gxy_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gxy_p = kind(gxy_p)
      REAL*8 gxy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gxy_p_p = kind(gxy_p_p)
      REAL*8 gxz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gxz = kind(gxz)
      REAL*8 gxz_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gxz_p = kind(gxz_p)
      REAL*8 gxz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gxz_p_p = kind(gxz_p_p)
      REAL*8 gyy (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gyy = kind(gyy)
      REAL*8 gyy_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gyy_p = kind(gyy_p)
      REAL*8 gyy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gyy_p_p = kind(gyy_p_p)
      REAL*8 gyz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gyz = kind(gyz)
      REAL*8 gyz_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gyz_p = kind(gyz_p)
      REAL*8 gyz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gyz_p_p = kind(gyz_p_p)
      REAL*8 gzz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gzz = kind(gzz)
      REAL*8 gzz_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gzz_p = kind(gzz_p)
      REAL*8 gzz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gzz_p_p = kind(gzz_p_p)
      REAL*8 kxx (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kxx = kind(kxx)
      REAL*8 kxx_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kxx_p = kind(kxx_p)
      REAL*8 kxx_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kxx_p_p = kind(kxx_p_p)
      REAL*8 kxy (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kxy = kind(kxy)
      REAL*8 kxy_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kxy_p = kind(kxy_p)
      REAL*8 kxy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kxy_p_p = kind(kxy_p_p)
      REAL*8 kxz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kxz = kind(kxz)
      REAL*8 kxz_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kxz_p = kind(kxz_p)
      REAL*8 kxz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kxz_p_p = kind(kxz_p_p)
      REAL*8 kyy (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kyy = kind(kyy)
      REAL*8 kyy_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kyy_p = kind(kyy_p)
      REAL*8 kyy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kyy_p_p = kind(kyy_p_p)
      REAL*8 kyz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kyz = kind(kyz)
      REAL*8 kyz_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kyz_p = kind(kyz_p)
      REAL*8 kyz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kyz_p_p = kind(kyz_p_p)
      REAL*8 kzz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kzz = kind(kzz)
      REAL*8 kzz_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kzz_p = kind(kzz_p)
      REAL*8 kzz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kzz_p_p = kind(kzz_p_p)
      REAL*8 psi (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psi = kind(psi)
      REAL*8 psix (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psix = kind(psix)
      REAL*8 psixx (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psixx = kind(psixx)
      REAL*8 psixy (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psixy = kind(psixy)
      REAL*8 psixz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psixz = kind(psixz)
      REAL*8 psiy (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psiy = kind(psiy)
      REAL*8 psiyy (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psiyy = kind(psiyy)
      REAL*8 psiyz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psiyz = kind(psiyz)
      REAL*8 psiz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psiz = kind(psiz)
      REAL*8 psizz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psizz = kind(psizz)
      REAL*8 r (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_r = kind(r)
      INTEGER*4 shift_state
      integer, parameter :: cctki_use_shift_state = kind(shift_state)
      REAL*8 x (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_x = kind(x)
      REAL*8 y (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_y = kind(y)
      REAL*8 z (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_z = kind(z)
      REAL*8 amp
      integer, parameter :: cctki_use_amp = kind(amp)
      REAL*8 c
      integer, parameter :: cctki_use_c = kind(c)
      REAL*8 eta0
      integer, parameter :: cctki_use_eta0 = kind(eta0)
      REAL*8 etamax
      integer, parameter :: cctki_use_etamax = kind(etamax)
      REAL*8 sigma
      integer, parameter :: cctki_use_sigma = kind(sigma)
      REAL*8 tolerance
      integer, parameter :: cctki_use_tolerance = kind(tolerance)
      INTEGER*4 interpolation_order
      integer, parameter :: cctki_use_interpolation_order = kind(interpo
     &lation_order)
      INTEGER*4 n
      integer, parameter :: cctki_use_n = kind(n)
      INTEGER*4 neta
      integer, parameter :: cctki_use_neta = kind(neta)
      INTEGER*4 nphi
      integer, parameter :: cctki_use_nphi = kind(nphi)
      INTEGER*4 ntheta
      integer, parameter :: cctki_use_ntheta = kind(ntheta)
      COMMON /DistortedBHIVPpriv/amp, c, eta0, etamax, sigma, tolerance,
     & interpolation_order, n, neta, nphi, ntheta
      integer*8 CCTKH0
      integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
      integer*8 CCTKH1
      integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
      integer*8 CCTKH2
      integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
      integer*8 CCTKH3
      integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
      integer*8 initial_data
      integer, parameter :: cctki_use_initial_data = kind(initial_data)
      integer*8 CCTKH4
      integer, parameter :: cctki_use_CCTKH4 = kind(CCTKH4)
      integer*8 CCTKH5
      integer, parameter :: cctki_use_CCTKH5 = kind(CCTKH5)
      integer*8 initial_lapse
      integer, parameter :: cctki_use_initial_lapse = kind(initial_lapse
     &)
      integer*8 CCTKH6
      integer, parameter :: cctki_use_CCTKH6 = kind(CCTKH6)
      integer*8 CCTKH7
      integer, parameter :: cctki_use_CCTKH7 = kind(CCTKH7)
      integer*8 CCTKH8
      integer, parameter :: cctki_use_CCTKH8 = kind(CCTKH8)
      integer*8 CCTKH10
      integer, parameter :: cctki_use_CCTKH10 = kind(CCTKH10)
      integer*8 metric_type
      integer, parameter :: cctki_use_metric_type = kind(metric_type)
      integer*8 CCTKH12
      integer, parameter :: cctki_use_CCTKH12 = kind(CCTKH12)
      integer*8 CCTKH13
      integer, parameter :: cctki_use_CCTKH13 = kind(CCTKH13)
      INTEGER*4 CCTKH9
      integer, parameter :: cctki_use_CCTKH9 = kind(CCTKH9)
      INTEGER*4 CCTKH11
      integer, parameter :: cctki_use_CCTKH11 = kind(CCTKH11)
      INTEGER*4 CCTKH14
      integer, parameter :: cctki_use_CCTKH14 = kind(CCTKH14)
      COMMON /ADMBASErest/CCTKH0, CCTKH1, CCTKH2, CCTKH3, initial_data, 
     &CCTKH4, CCTKH5, initial_lapse, CCTKH6, CCTKH7, CCTKH8, CCTKH10, me
     &tric_type, CCTKH12, CCTKH13, CCTKH9, CCTKH11, CCTKH14
      integer*8 conformal_storage
      integer, parameter :: cctki_use_conformal_storage = kind(conformal
     &_storage)
      COMMON /STATICCONFORMALrest/conformal_storage
      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

      real*8 :: deta,dq,dphi
      real*8, allocatable :: ac(:,:,:),ae(:,:,:),aw(:,:,:),an(:,:,:),
     $     as(:,:,:),aq(:,:,:),ab(:,:,:),rhs(:,:,:),
     $     qf(:,:,:),qfetaeta(:,:,:),qfqq(:,:,:),qfphi(:,:,:),
     $     qfphiphi(:,:,:),
     $     psisph(:,:,:),psiprim(:,:,:),detapsisph(:,:,:),
     $     dqpsisph(:,:,:),dphipsisph(:,:,:),detaetapsisph(:,:,:),
     $     detaqpsisph(:,:,:),detaphipsisph(:,:,:),dqqpsisph(:,:,:),
     $     dqphipsisph(:,:,:),dphiphipsisph(:,:,:)
      real*8, allocatable :: etagrd(:),qgrd(:),phigrd(:)
      real*8 o1,o2,o3,o4,o5,o6,o7,o8,o9,o10,o11,o12,o13,o14,o15,o16,o17,
     &
     $     o18,o19,o20,o21,o22,o23,o24,o25,o26,o27,o28,o29,o30,o31,o32,
     $     o33,o34,o35,o36,o37,o38,o39,o40,o41,o42,o43,o44,o45,o46,o47,
     $     o48,o49,o50,o51,o52,o53,o54,o55,o56,o57,o58,o59,o60,o61,o62,
     $     o63,o64
      real*8 rmax,adm
      real*8 dbh_eps
      real*8 pi
      integer :: ne, nq, np
      integer :: nx,ny,nz
      integer i,j,k,ier,nquads,nocts
      integer npoints,ierror
      integer make_conformal_derivs

      integer       param_table_handle, interp_handle
      character(30) options_string
      REAL*8,    dimension(3)  :: coord_origin, coord_delta
      integer*8, dimension(3)  :: interp_coords
      INTEGER*4,     dimension(3 ) :: in_array_dims
      integer*8, dimension(10) :: in_arrays, out_arrays
      INTEGER*4,     dimension(10), parameter :: type_codes = 130

      dbh_eps = tolerance

c

      if ((CCTK_Equals(metric_type, "static conformal") .ne. 0)) then

         conformal_state = 1
      
         if((CCTK_Equals(conformal_storage,"factor+derivs") .ne. 0)) the
     &n
         
            conformal_state = 2
            make_conformal_derivs = 1

          else if ((CCTK_Equals(conformal_storage,"factor+derivs+2nd der
     &ivs") .ne. 0)) then
          
            conformal_state = 3
            make_conformal_derivs = 1

          endif
          
       else
       
         make_conformal_derivs = 0
       
       endif
      
      pi = 4.0d0*atan(1.0d0)

c
      ne = neta
      nq = ntheta
      np = nphi

c
      nx = cctk_lsh(1)
      ny = cctk_lsh(2)
      nz = cctk_lsh(3)      
      
c
      
      print *,"Brill wave + Distorted BH solve"
      write(*,123)amp,eta0,c,sigma,n
      print*,'etamax=',etamax
 123  format(1x, 'Pars: amp',f8.5,' eta0',f8.5,' c',f8.5,' sigma',f8.5,'
     & n',i3)
      
c
c
c
c
      ne = ne + 2
      nq = nq + 2
      np = np + 2
c
      allocate(ac(ne,nq,np),ae(ne,nq,np),aw(ne,nq,np),an(ne,nq,np),
     $     as(ne,nq,np),aq(ne,nq,np),ab(ne,nq,np),rhs(ne,nq,np),
     $     qf(ne,nq,np),qfetaeta(ne,nq,np),qfqq(ne,nq,np),
     $     qfphi(ne,nq,np),qfphiphi(ne,nq,np),
     $     psisph(ne,nq,np),psiprim(ne,nq,np),detapsisph(ne,nq,np),
     $     dqpsisph(ne,nq,np),dphipsisph(ne,nq,np),
     $     detaetapsisph(ne,nq,np),detaqpsisph(ne,nq,np),
     $     detaphipsisph(ne,nq,np),dqqpsisph(ne,nq,np),
     $     dqphipsisph(ne,nq,np),dphiphipsisph(ne,nq,np))
      allocate(etagrd(ne),qgrd(nq),phigrd(np))
c
c
c
      nocts = 4
      nquads = 2      
      dphi = nocts*0.5*pi/(np-2)
      dq = nquads*0.5*pi/(nq-2)
      deta = etamax/(ne-3)
      
      do k = 1,np
         phigrd(k) = (k-1.5)*dphi
      enddo
      do j = 1,nq
         qgrd(j) = (j-1.5)*dq
      enddo
      do i=1,ne
         etagrd(i) = (i-2)*deta
      enddo
c
c
c
      do k = 1,np
         do j = 1,nq
            do i = 1,ne

	qf(i,j,k) = amp*(1.00000000000000d0+c*cos(phigrd(k))**2)*
     $	(exp(-((etagrd(i)-eta0)**2/sigma**2))+exp(-((etagrd(i)+eta0)**2
     $	/sigma**2)))*sin(qgrd(j))**n

	qfetaeta(i,j,k) = amp*(1.00000000000000d0+c*cos(phigrd(k))
     $	**2)*((4.0000000000000d0*(etagrd(i)-eta0)**2*exp(-((etagrd(i)-eta
     &0)
     $	**2/sigma**2)))/sigma**4-(2.00000000000000d0*exp(-((etagrd(i)
     $	-eta0)**2/sigma**2)))/sigma**2+(4.0000000000000d0*(etagrd(i)+
     $	eta0)**2*exp(-((etagrd(i)+eta0)**2/sigma**2)))/sigma**4-(2.00
     $	000000000000d0*exp(-((etagrd(i)+eta0)**2/sigma**2)))/sigma**2
     $	)*sin(qgrd(j))**n

	qfqq(i,j,k) = amp*(-1.00000000000000d0+n)*n*(1.00000000
     $	000000d0+c*cos(phigrd(k))**2)*cos(qgrd(j))**2*(exp(-((etagrd(i)-e
     &ta0)**2/
     $	sigma**2))+exp(-((etagrd(i)+eta0)**2/sigma**2)))*sin(qgrd(j))**(-
     &2
     $	+n)-amp*n*(1.00000000000000d0+c*cos(phigrd(k))
     $	**2)*(exp(-((etagrd(i)-eta0)**2/sigma**2))+exp(-((etagrd(i)+eta0)
     &**2
     $	/sigma**2)))*sin(qgrd(j))**n

	qfphi(i,j,k) = -2.00000000000000d0*amp*c*cos(phigrd(k))*(exp(-
     $	((etagrd(i)-eta0)**2/sigma**2))+exp(-((etagrd(i)+eta0)**2/sigma**
     &2))
     $	)*sin(phigrd(k))*sin(qgrd(j))**n

	qfphiphi(i,j,k) = -2.00000000000000d0*amp*c*cos(phigrd(k))**2*
     $	(exp(-((etagrd(i)-eta0)**2/sigma**2))+exp(-((etagrd(i)+eta0)**2/
     $	sigma**2)))*sin(qgrd(j))**n+2.00000000000000d0*amp*c*(exp(-((
     $	etagrd(i)-eta0)**2/sigma**2))+exp(-((etagrd(i)+eta0)**2/sigma**2)
     &))*
     $	sin(phigrd(k))**2*sin(qgrd(j))**n

            enddo
         enddo
      enddo
      
c
c
c
      psiprim = 0.
      
      do k = 1,np
         do j = 1,nq
            do i = 1,ne
               psisph(i,j,k) = 2.*cosh(0.5*etagrd(i))
            enddo
         enddo
      enddo
c
c
c
      ac = 0.
      ae = 0.
      aw = 0.
      an = 0.
      as = 0.
      aq = 0.
      ab = 0.
      rhs = 0.
c

      do k = 2,np-1
         do j = 2,nq-1
            do i = 2,ne-1
               ac(i,j,k) = -2./deta**2-2./dq**2-2.*exp(2.*
     &              qf(i,j,k))/(dphi**2*
     &              sin(qgrd(j))**2)+0.25*
     &              (qfetaeta(i,j,k)+qfqq(i,j,k)+2.*
     &              exp(2.*qf(i,j,k))*qfphiphi(i,j,k)/
     &              sin(qgrd(j))**2+3.*exp(2.*
     &              qf(i,j,k))*qfphi(i,j,k)**2/
     &              sin(qgrd(j))**2-1.)
               ae(i,j,k) = 1./deta**2
               aw(i,j,k) = 1./deta**2
               an(i,j,k) = 1./dq**2+0.5/(dq*tan(qgrd(j)))
               
               as(i,j,k) = 1./dq**2-0.5/(dq*tan(qgrd(j)))
               
               aq(i,j,k) = exp(2.*qf(i,j,k))/(dphi**2*
     &              sin(qgrd(j))**2)+exp(2.*
     &              qf(i,j,k))*qfphi(i,j,k)/(dphi*
     &              sin(qgrd(j))**2)
               
               ab(i,j,k) = exp(2.*qf(i,j,k))/(dphi**2*
     &              sin(qgrd(j))**2)-exp(2.*
     &              qf(i,j,k))*qfphi(i,j,k)/(dphi*
     &              sin(qgrd(j))**2)
               
               rhs(i,j,k) = -0.25*(qfetaeta(i,j,k)+
     &              qfqq(i,j,k)+2.*exp(2.*qf(i,j,k))*
     &              qfphiphi(i,j,k)/sin(qgrd(j))**2+3.*
     &              exp(2.*qf(i,j,k))*qfphi(i,j,k)**2/
     &              sin(qgrd(j))**2)*psisph(i,j,k)
            enddo
         enddo
      enddo
c
c
c
c
      do k = 3,np-2
         do j = 3,nq-2
            ae(2,j,k) = ae(2,j,k) + aw(2,j,k)
            aw(2,j,k) = 0.
c
c
            ac(ne-1,j,k) = ac(ne-1,j,k)+4.*ae(ne-1,j,k)/(3.+deta)
            aw(ne-1,j,k) = aw(ne-1,j,k) - ae(ne-1,j,k)/(3.+deta)
            ae(ne-1,j,k) = 0.
         enddo
      enddo
c
c
      do k = 3,np-2
         do i = 3,ne-2
            ac(i,2,k) = ac(i,2,k) + as(i,2,k)
            as(i,2,k) = 0.
c
c
            ac(i,nq-1,k) = ac(i,nq-1,k) + an(i,nq-1,k)
            an(i,nq-1,k) = 0.
         enddo
      enddo
c
c
      do j = 3,nq-2
         do i = 3,ne-2
            ac(i,j,2) = ac(i,j,2) + ab(i,j,2)
            ab(i,j,2) = 0.
c
c
            ac(i,j,np-1) = ac(i,j,np-1) + aq(i,j,np-1)
            aq(i,j,np-1) = 0.
         enddo
      enddo
c
c
c
c
      do k = 3,np-2
         ae(2,2,k) = ae(2,2,k) + aw(2,2,k)
         ac(2,2,k) = ac(2,2,k) + as(2,2,k)
         aw(2,2,k) = 0.
         as(2,2,k) = 0.
c
c
         aw(ne-1,2,k) = aw(ne-1,2,k) - ae(ne-1,2,k)/(3.+deta)
         ac(ne-1,2,k) = ac(ne-1,2,k) + as(ne-1,2,k) +
     &        4.*ae(ne-1,2,k)/(3.+deta)
         ae(ne-1,2,k) = 0.
         as(ne-1,2,k) = 0.
c
c
         ae(2,nq-1,k) = ae(2,nq-1,k) + aw(2,nq-1,k)
         ac(2,nq-1,k) = ac(2,nq-1,k) + an(2,nq-1,k)
         aw(2,nq-1,k) = 0.
         an(2,nq-1,k) = 0.
c
c
         aw(ne-1,nq-1,k) = aw(ne-1,nq-1,k) - ae(ne-1,nq-1,k)/
     &        (3.+deta)
         ac(ne-1,nq-1,k) = ac(ne-1,nq-1,k) + an(ne-1,nq-1,k) +
     &        4.*ae(ne-1,nq-1,k)/(3.+deta)
         ae(ne-1,nq-1,k) = 0.
         an(ne-1,nq-1,k) = 0.
      enddo
c
c
      do j = 3,nq-2
         ae(2,j,2) = ae(2,j,2) + aw(2,j,2)
         ac(2,j,2) = ac(2,j,2) + ab(2,j,2)
         aw(2,j,2) = 0.
         ab(2,j,2) = 0.
c
c
         aw(ne-1,j,2) = aw(ne-1,j,2) - ae(ne-1,j,2)/(3.+deta)
         ac(ne-1,j,2) = ac(ne-1,j,2) + ab(ne-1,j,2) +
     &        4.*ae(ne-1,j,2)/(3.+deta)
         ae(ne-1,j,2) = 0.
         ab(ne-1,j,2) = 0.
c
c
         ae(2,j,np-1) = ae(2,j,np-1) + aw(2,j,np-1)
         ac(2,j,np-1) = ac(2,j,np-1) + aq(2,j,np-1)
         aw(2,j,np-1) = 0.
         aq(2,j,np-1) = 0.
c
c
         aw(ne-1,j,np-1) = aw(ne-1,j,np-1) - ae(ne-1,j,np-1)/
     &        (3.+deta)
         ac(ne-1,j,np-1) = ac(ne-1,j,np-1) + aq(ne-1,j,np-1) +
     &        4.*ae(ne-1,j,np-1)/(3.+deta)
         ae(ne-1,j,np-1) = 0.
         aq(ne-1,j,np-1) = 0.
      enddo
c
c
      do i = 3,ne-2
         ac(i,2,2) = ac(i,2,2) + as(i,2,2) + ab(i,2,2)
         as(i,2,2) = 0.
         ab(i,2,2) = 0.
c
c
         ac(i,nq-1,2) = ac(i,nq-1,2) + an(i,nq-1,2) +
     &        ab(i,nq-1,2)
         an(i,nq-1,2) = 0.
         ab(i,nq-1,2) = 0.
c
c
         ac(i,2,np-1) = ac(i,2,np-1) + as(i,2,np-1) +
     &        aq(i,2,np-1)
         as(i,2,np-1) = 0.
         aq(i,2,np-1) = 0.
c
c
         ac(i,nq-1,np-1) = ac(i,nq-1,np-1) + an(i,nq-1,np-1) +
     &        aq(i,nq-1,np-1)
         an(i,nq-1,np-1) = 0.
         aq(i,nq-1,np-1) = 0.
      enddo
c
c
c
c
      ae(2,2,2) = ae(2,2,2) + aw(2,2,2)
      ac(2,2,2) = ac(2,2,2) + as(2,2,2) + ab(2,2,2)
      aw(2,2,2) = 0.
      as(2,2,2) = 0.
      ab(2,2,2) = 0.
c
c
      aw(ne-1,2,2) = aw(ne-1,2,2) - ae(ne-1,2,2)/(3.+deta)
      ac(ne-1,2,2) = ac(ne-1,2,2) + as(ne-1,2,2) + ab(ne-1,2,2) +
     &     4.*ae(ne-1,2,2)/(3.+deta)
      ae(ne-1,2,2) = 0.
      as(ne-1,2,2) = 0.
      ab(ne-1,2,2) = 0.
c
c
      ae(2,nq-1,2) = ae(2,nq-1,2) + aw(2,nq-1,2)
      ac(2,nq-1,2) = ac(2,nq-1,2) + an(2,nq-1,2) + ab(2,nq-1,2)
      aw(2,nq-1,2) = 0.
      an(2,nq-1,2) = 0.
      ab(2,nq-1,2) = 0.
c
c
      ae(2,2,np-1) = ae(2,2,np-1) + aw(2,2,np-1)
      ac(2,2,np-1) = ac(2,2,np-1) + as(2,2,np-1) + aq(2,2,np-1)
      aw(2,2,np-1) = 0.
      as(2,2,np-1) = 0.
      aq(2,2,np-1) = 0.
c
c
      aw(ne-1,nq-1,2) = aw(ne-1,nq-1,2) - ae(ne-1,nq-1,2)/(3.+deta)
      ac(ne-1,nq-1,2) = ac(ne-1,nq-1,2) + an(ne-1,nq-1,2) + ab(ne-1,nq-1
     &,2) +
     &     4.*ae(ne-1,nq-1,2)/(3.+deta)
      ae(ne-1,nq-1,2) = 0.
      an(ne-1,nq-1,2) = 0.
      ab(ne-1,nq-1,2) = 0.
c
c
      aw(ne-1,2,np-1) = aw(ne-1,2,np-1) - ae(ne-1,2,np-1)/(3.+deta)
      ac(ne-1,2,np-1) = ac(ne-1,2,np-1) + as(ne-1,2,np-1) + aq(ne-1,2,np
     &-1) +
     &     4.*ae(ne-1,2,np-1)/(3.+deta)
      ae(ne-1,2,np-1) = 0.
      as(ne-1,2,np-1) = 0.
      aq(ne-1,2,np-1) = 0.
c
c
      ae(2,nq-1,np-1) = ae(2,nq-1,np-1) + aw(2,nq-1,np-1)
      ac(2,nq-1,np-1) = ac(2,nq-1,np-1) + an(2,nq-1,np-1) + aq(2,nq-1,np
     &-1)
      aw(2,nq-1,np-1) = 0.
      an(2,nq-1,np-1) = 0.
      aq(2,nq-1,np-1) = 0.
c
c
      aw(ne-1,nq-1,np-1) = aw(ne-1,nq-1,np-1) - ae(ne-1,nq-1,np-1)/(3.+d
     &eta)
      ac(ne-1,nq-1,np-1) = ac(ne-1,nq-1,np-1) + an(ne-1,nq-1,np-1) +
     &     aq(ne-1,nq-1,np-1) + 4.*ae(ne-1,nq-1,np-1)/(3.+deta)
      ae(ne-1,nq-1,np-1) = 0.
      an(ne-1,nq-1,np-1) = 0.
      aq(ne-1,nq-1,np-1) = 0.
c
c
c
      call bicgst3d(ac,ae,aw,an,as,aq,ab,psiprim,rhs,dbh_eps,rmax,ier,
     $     ne,nq,np)
c
      if (rmax.gt.dbh_eps) then
         write(*,*) '***WARNING: bicgst3d did not converge.'
      endif
      if (ier.eq.-1) then
         write(*,*) '***WARNING: ier=-1'
      endif
      print *,'psiprim = ',maxval(psiprim),' ',minval(psiprim)
c
c
c
      do k = 1,np
         do j = 1,nq
            psiprim(1,j,k) = psiprim(3,j,k)
            psiprim(ne,j,k) = (4.*psiprim(ne-1,j,k)-psiprim(ne-2,j,k))/
     $           (3.+deta)
         enddo
      enddo
      do k = 1,np
         do i = 1,ne
            psiprim(i,1,k) = psiprim(i,2,k)
            psiprim(i,nq,k) = psiprim(i,nq-1,k)
         enddo
      enddo
      do j = 1,nq
         do i = 1,ne
            psiprim(i,j,1) = psiprim(i,j,2)
            psiprim(i,j,np) = psiprim(i,j,np-1)
         enddo
      enddo
c
c
c
c

      do k = 1,np
         do j = 1,nq
            do i = 2,ne-1
               detapsisph(i,j,k)=0.5*(psiprim(i+1,j,k)-psiprim(i-1,j,k))
     &
     $              /deta + sinh(0.5*etagrd(i))
            enddo
            detapsisph(1,j,k) = -detapsisph(3,j,k)
         enddo
      enddo
c
      do k = 1,np
         do j = 2,nq-1
            do i = 1,ne
               dqpsisph(i,j,k)=0.5*(psiprim(i,j+1,k)-psiprim(i,j-1,k))/
     $              dq
            enddo
         enddo
         do i = 1,ne
            dqpsisph(i,1,k) = -dqpsisph(i,2,k)
            dqpsisph(i,nq,k) = -dqpsisph(i,nq-1,k)
         enddo
      enddo
c
      do k = 2,np-1
         do j = 1,nq
            do i = 1,ne
               dphipsisph(i,j,k)=0.5*(psiprim(i,j,k+1)-psiprim(i,j,k-1))
     &
     $              /dphi
            enddo
         enddo
      enddo
      do j = 1,nq
         do i = 1,ne
            dphipsisph(i,j,1) = -dphipsisph(i,j,2)
            dphipsisph(i,j,np) = -dphipsisph(i,j,np-1)
         enddo
      enddo
c
      do k = 1,np
         do j = 1,nq
            do i = 2,ne-1
               detaetapsisph(i,j,k)=(psiprim(i+1,j,k)-2.*psiprim(i,j,k)+
     &
     &              psiprim(i-1,j,k))/deta**2+sqrt(0.25)*
     &              cosh(0.5*etagrd(i))
            enddo
            detaetapsisph(1,j,k) = detaetapsisph(3,j,k)
         enddo
      enddo
c
      do k = 1,np
         do j = 2,nq-1
            do i = 1,ne
               detaqpsisph(i,j,k)=0.5*(detapsisph(i,j+1,k)-detapsisph(i,
     &
     $              j-1,k))/dq
            enddo
         enddo
         do i = 1,ne
            detaqpsisph(i,1,k) = -detaqpsisph(i,2,k)
            detaqpsisph(i,nq,k) = -detaqpsisph(i,nq-1,k)
         enddo
      enddo
c
      do k = 2,np-1
         do j = 1,nq
            do i = 1,ne
               detaphipsisph(i,j,k)=0.5*(detapsisph(i,j,k+1)-detapsisph(
     &
     $              i,j,k-1))/dphi
            enddo
         enddo
      enddo
      do j = 1,nq
         do i = 1,ne
            detaphipsisph(i,j,1) = -detaphipsisph(i,j,2)
            detaphipsisph(i,j,np) = -detaphipsisph(i,j,np-1)
         enddo
      enddo
c
      do k = 1,np
         do j = 2,nq-1
            do i = 1,ne
               dqqpsisph(i,j,k)=0.5*(dqpsisph(i,j+1,k)-dqpsisph(i,j-1,k)
     &)/
     $              dq
            enddo
         enddo
         do i = 1,ne
            dqqpsisph(i,1,k) = dqqpsisph(i,2,k)
            dqqpsisph(i,nq,k) = dqqpsisph(i,nq-1,k)
         enddo
      enddo
c
      do k = 2,np-1
         do j = 1,nq
            do i = 1,ne
               dqphipsisph(i,j,k)=0.5*(dqpsisph(i,j,k+1)-dqpsisph(i,j,k-
     &1))/
     $              dphi
            enddo
         enddo
      enddo
      do j = 1,nq
         do i = 1,ne
            dqphipsisph(i,j,1) = -dqphipsisph(i,j,2)
            dqphipsisph(i,j,np) = -dqphipsisph(i,j,np-1)
         enddo
      enddo
c
      do k = 2,np-1
         do j = 1,nq
            do i = 1,ne
               dphiphipsisph(i,j,k)=0.5*(dphipsisph(i,j,k+1)-
     $              dphipsisph(i,j,k-1))/dphi
            enddo
         enddo
      enddo
      do j = 1,nq
         do i = 1,ne
            dphiphipsisph(i,j,1) = dphiphipsisph(i,j,2)
            dphiphipsisph(i,j,np) = dphiphipsisph(i,j,np-1)
         enddo
      enddo
c
      do k = 1,np
         do j = 1,nq
            psisph(:,j,k)=psiprim(:,j,k)+2.0*cosh(0.5*etagrd)
         enddo
      enddo
      
c
c
c
c

      eta = 0.5d0*log(x**2+y**2+z**2)
      abseta = abs(eta)
      q = atan2(sqrt(x**2+y**2),z)
      phi = atan2(y,x)

      do k = 1,nz
         do j = 1,ny
            do i = 1,nx

               if(eta(i,j,k) .lt. 0)then
                  sign_eta(i,j,k) = -1.0d0
               else
                  sign_eta(i,j,k) = 1.0d0
               endif

            enddo
         enddo
      enddo
      
      
      npoints = nx*ny*nz

c
      options_string = "order = " // char(ichar('0') + interpolation_ord
     &er)
      call Util_TableCreateFromString (param_table_handle, options_strin
     &g)
      if (param_table_handle .lt. 0) then
        call CCTK_Warn(0,606,"DistortedBHIVP.F","DistortedBHIVP","Cannot
     & create parameter table for interpolator")
      endif

      call CCTK_InterpHandle (interp_handle, "uniform cartesian")
      if (interp_handle .lt. 0) then
        call CCTK_Warn(0,611,"DistortedBHIVP.F","DistortedBHIVP","Cannot
     & get handle for interpolation ! Forgot to activate an implementati
     &on providing interpolation operators ??")
      endif

c
      coord_origin(1) = etagrd(1)
      coord_origin(2) = qgrd(1)
      coord_origin(3) = phigrd(1)-pi
      coord_delta(1)  = etagrd(2) - etagrd(1)
      coord_delta(2)  = qgrd(2) - qgrd(1)
      coord_delta(3)  = phigrd(2) - phigrd(1)

      interp_coords(1) = CCTK_PointerTo(abseta)
      interp_coords(2) = CCTK_PointerTo(q)
      interp_coords(3) = CCTK_PointerTo(phi)

      in_array_dims(1) = ne
      in_array_dims(2) = nq
      in_array_dims(3) = np

      in_arrays(1) = CCTK_PointerTo(psisph)
      in_arrays(2) = CCTK_PointerTo(detapsisph)
      in_arrays(3) = CCTK_PointerTo(dqpsisph)
      in_arrays(4) = CCTK_PointerTo(dphipsisph)
      in_arrays(5) = CCTK_PointerTo(detaetapsisph)
      in_arrays(6) = CCTK_PointerTo(detaqpsisph)
      in_arrays(7) = CCTK_PointerTo(detaphipsisph)
      in_arrays(8) = CCTK_PointerTo(dqqpsisph)
      in_arrays(9) = CCTK_PointerTo(dqphipsisph)
      in_arrays(10) = CCTK_PointerTo(dphiphipsisph)

      out_arrays(1) = CCTK_PointerTo(psi3d)
      out_arrays(2) = CCTK_PointerTo(detapsi3d)
      out_arrays(3) = CCTK_PointerTo(dqpsi3d)
      out_arrays(4) = CCTK_PointerTo(dphipsi3d)
      out_arrays(5) = CCTK_PointerTo(detaetapsi3d)
      out_arrays(6) = CCTK_PointerTo(detaqpsi3d)
      out_arrays(7) = CCTK_PointerTo(detaphipsi3d)
      out_arrays(8) = CCTK_PointerTo(dqqpsi3d)
      out_arrays(9) = CCTK_PointerTo(dqphipsi3d)
      out_arrays(10) = CCTK_PointerTo(dphiphipsi3d)

      call CCTK_InterpLocalUniform (ierror, 3,
     $                              interp_handle, param_table_handle,
     $                              coord_origin, coord_delta,
     $                              npoints, type_codes(1), interp_coord
     &s,
     $                              10, in_array_dims, type_codes, in_ar
     &rays,
     $                              10, type_codes, out_arrays)

      psi = psi3d*exp(-0.5*eta)
      detapsi3d = sign_eta*detapsi3d
      detaqpsi3d = sign_eta*detaqpsi3d
      detaphipsi3d = sign_eta*detaphipsi3d

      do k=1,nz
         do j=1,ny
            do i=1,nx
c

        o1 = 5.0000000000000d-1*eta(i,j,k)
        o2 = exp(o1)
        o3 = psi3d(i,j,k)
        o4 = 1/o3
        o5 = cos(phi(i,j,k))
        o6 = cos(q(i,j,k))
        o7 = dqpsi3d(i,j,k)
        o8 = -1.50000000000000d0*eta(i,j,k)
        o9 = exp(o8)
        o10 = dphipsi3d(i,j,k)
        o11 = sin(phi(i,j,k))
        o12 = sin(q(i,j,k))
        o13 = 1/o12
        o14 = detapsi3d(i,j,k)
        psix(i,j,k) = o2*o4*(-(o10*o11*o13*o9) + o12*o14*o5*o9 - 5.00000
     &
     &  00000000d-1*o12*o3*o5*o9 + o5*o6*o7*o9)
        psiy(i,j,k) = o2*o4*(o11*o12*o14*o9 - 5.0000000000000d-1*o11*o12
     &
     &  *o3*o9 + 1.00000000000000d0*o10*o13*o5*o9 + o11*o6*o7*o9)
        psiz(i,j,k) = o2*o4*(o14*o6*o9 - 5.0000000000000d-1*o3*o6*o9 - o
     &
     &  12*o7*o9)
c

        o1 = 5.0000000000000d-1*eta(i,j,k)
        o2 = exp(o1)
        o3 = psi3d(i,j,k)
        o4 = 1/o3
        o5 = cos(phi(i,j,k))
        o6 = o5**2
        o7 = cos(q(i,j,k))
        o8 = o7**2
        o9 = detapsi3d(i,j,k)
        o10 = -2.50000000000000d0*eta(i,j,k)
        o11 = exp(o10)
        o12 = dqqpsi3d(i,j,k)
        o13 = detaphipsi3d(i,j,k)
        o14 = sin(phi(i,j,k))
        o15 = dphipsi3d(i,j,k)
        o16 = o14**2
        o17 = sin(q(i,j,k))
        o18 = o17**2
        o19 = 1/o18
        o20 = dphiphipsi3d(i,j,k)
        o21 = detaqpsi3d(i,j,k)
        o22 = dqpsi3d(i,j,k)
        o23 = detaetapsi3d(i,j,k)
        o24 = tan(q(i,j,k))		
        o25 = o24**2	
        o26 = 1/o25
        o27 = dqphipsi3d(i,j,k)	
        o28 = 1/o24
        psixx(i,j,k) = o2*o4*(1.00000000000000d0*o11*o16*o19*o20 + 1.000
     &
     &  00000000000d0*o11*o16*o22*o28 - 5.0000000000000d-1*o11*o16*o3 - 
     &
     &  2.00000000000000d0*o11*o13*o14*o5 + 2.00000000000000d0*o11*o14*o
     &
     &  15*o5 + 1.00000000000000d0*o11*o14*o15*o19*o5 + 1.00000000000000
     &
     &  d0*o11*o14*o15*o26*o5 - 2.00000000000000d0*o11*o14*o27*o28*o5 + 
     &
     &  o11*o18*o23*o6 + 7.5000000000000d-1*o11*o18*o3*o6 + 2.0000000000
     &
     &  0000d0*o11*o17*o21*o6*o7 - 3.00000000000000d0*o11*o17*o22*o6*o7 
     &
     &  + o11*o12*o6*o8 - 5.0000000000000d-1*o11*o3*o6*o8 + o11*o16*o9 -
     &
     &   2.00000000000000d0*o11*o18*o6*o9 + o11*o6*o8*o9)
        psixy(i,j,k) = o2*o4*(-(o11*o13*o16) + 1.50000000000000d0*o11*o1
     &
     &  5*o16 + 1.00000000000000d0*o11*o15*o16*o26 - o11*o16*o27*o28 - o
     &
     &  11*o14*o19*o20*o5 + o11*o14*o18*o23*o5 - o11*o14*o22*o28*o5 + 5.
     &
     &  0000000000000d-1*o11*o14*o3*o5 + 7.5000000000000d-1*o11*o14*o18*
     &
     &  o3*o5 + o11*o13*o6 - 5.0000000000000d-1*o11*o15*o6 - o11*o15*o19
     &
     &  *o6 + 1.00000000000000d0*o11*o27*o28*o6 + 2.00000000000000d0*o11
     &
     &  *o14*o17*o21*o5*o7 - 3.00000000000000d0*o11*o14*o17*o22*o5*o7 + 
     &
     &  o11*o12*o14*o5*o8 - 5.0000000000000d-1*o11*o14*o3*o5*o8 - o11*o1
     &
     &  4*o5*o9 - 2.00000000000000d0*o11*o14*o18*o5*o9 + o11*o14*o5*o8*o
     &
     &  9)
        psixz(i,j,k) = o2*o4*(o11*o14*o27 - o11*o13*o14*o28 + 5.00000000
     &
     &  00000d-1*o11*o14*o15*o28 - o11*o18*o21*o5 + 1.50000000000000d0*o
     &
     &  11*o18*o22*o5 - o11*o12*o17*o5*o7 + o11*o17*o23*o5*o7 + 1.250000
     &
     &  00000000d0*o11*o17*o3*o5*o7 + o11*o21*o5*o8 - 1.50000000000000d0
     &
     &  *o11*o22*o5*o8 - 3.00000000000000d0*o11*o17*o5*o7*o9)
        psiyy(i,j,k) = o2*o4*(o11*o16*o18*o23 + 7.5000000000000d-1*o11*o
     &
     &  16*o18*o3 + 2.00000000000000d0*o11*o13*o14*o5 - 2.00000000000000
     &
     &  d0*o11*o14*o15*o5 - o11*o14*o15*o19*o5 - o11*o14*o15*o26*o5 + 2.
     &
     &  00000000000000d0*o11*o14*o27*o28*o5 + 1.00000000000000d0*o11*o19
     &
     &  *o20*o6 + 1.00000000000000d0*o11*o22*o28*o6 - 5.0000000000000d-1
     &
     &  *o11*o3*o6 + 2.00000000000000d0*o11*o16*o17*o21*o7 - 3.000000000
     &
     &  00000d0*o11*o16*o17*o22*o7 + o11*o12*o16*o8 - 5.0000000000000d-1
     &
     &  *o11*o16*o3*o8 - 2.00000000000000d0*o11*o16*o18*o9 + o11*o6*o9 +
     &
     &   o11*o16*o8*o9)
        psiyz(i,j,k) = o2*o4*(-(o11*o14*o18*o21) + 1.50000000000000d0*o1
     &
     &  1*o14*o18*o22 - o11*o27*o5 + 1.00000000000000d0*o11*o13*o28*o5 -
     &
     &   5.0000000000000d-1*o11*o15*o28*o5 - o11*o12*o14*o17*o7 + o11*o1
     &
     &  4*o17*o23*o7 + 1.25000000000000d0*o11*o14*o17*o3*o7 + o11*o14*o2
     &
     &  1*o8 - 1.50000000000000d0*o11*o14*o22*o8 - 3.00000000000000d0*o1
     &
     &  1*o14*o17*o7*o9)
        psizz(i,j,k) = o2*o4*(o11*o12*o18 - 5.0000000000000d-1*o11*o18*o
     &
     &  3 - 2.00000000000000d0*o11*o17*o21*o7 + 3.00000000000000d0*o11*o
     &
     &  17*o22*o7 + o11*o23*o8 + 7.5000000000000d-1*o11*o3*o8 + o11*o18*
     &
     &  o9 - 2.00000000000000d0*o11*o8*o9)
            enddo
         enddo
      enddo
      
c
c
      
c
c
c
      do k=1,nz
         do j=1,ny
            do i=1,nx

        o1 = 2.00000000000000d0*phi(i,j,k)
        o2 = cos(o1)
        o3 = cos(phi(i,j,k))
        o4 = o3**2
        o5 = c*o4
        o6 = 1.00000000000000d0 + o5
        o7 = -eta0
        o8 = eta(i,j,k) + o7
        o9 = o8**2
        o10 = sigma**2
        o11 = 1/o10
        o12 = -(o11*o9)
        o13 = exp(o12)
        o14 = eta(i,j,k) + eta0
        o15 = o14**2
        o16 = -(o11*o15)
        o17 = exp(o16)
        o18 = o13 + o17
        o19 = sin(q(i,j,k))
        o20 = o19**n
        o21 = 2.00000000000000d0*amp*o18*o20*o6
        o22 = exp(o21)
        o23 = -1.00000000000000d0 + o22
        o24 = -eta(i,j,k)
        o25 = exp(o24)
        o26 = sin(phi(i,j,k))
        o27 = o26**2
        o28 = 1/o19
        o29 = sin(o1)
        o30 = cos(q(i,j,k))
        o31 = o30**2
        o32 = -1.00000000000000d0 + n
        o33 = o19**o32
        o34 = amp*n*o18*o22*o3*o31*o33*o6
        o35 = 2.00000000000000d0*amp*c*o18*o22*o27*o3*o33
        o36 = 4.0000000000000d0*eta(i,j,k)*eta0*o11
        o37 = exp(o36)
        o38 = eta(i,j,k)*o37
        o39 = -(eta0*o37)
        o40 = eta(i,j,k) + eta0 + o38 + o39
        o41 = o16 + o21
        o42 = exp(o41)
        o43 = 1.00000000000000d0 + n
        o44 = o19**o43
        o45 = -2.00000000000000d0*amp*o11*o3*o40*o42*o44*o6
        o46 = -2.00000000000000d0*amp*c*o18*o22*o26*o33*o4
        o47 = amp*n*o18*o22*o26*o31*o33*o6
        o48 = -2.00000000000000d0*amp*o11*o26*o40*o42*o44*o6
        o49 = c*o2
        o50 = 2.00000000000000d0 + c + o49
        o51 = -1.00000000000000d0 + o37
        o52 = -2.00000000000000d0*eta0*o51
        o53 = 1.00000000000000d0 + o37
        o54 = 2.00000000000000d0*eta(i,j,k)*o53
        o55 = n*o10*o53
        o56 = o52 + o54 + o55
        o57 = o16 + o21 + o24
        o58 = exp(o57)
        o59 = o11*o9
        o60 = exp(o59)
        o61 = o11*o15
        o62 = exp(o61)
        o63 = o60 + o62
        o64 = o19**2
        gxx(i,j,k) = 5.0000000000000d-1*(1.00000000000000d0 + o22 + o2*o
     &
     &  23)
c
c
c
c
c
c
c
c
c
c
c
        gxy(i,j,k) = 5.0000000000000d-1*o23*o29
c
c
c
c
c
c
c
c
c
c
        gxz(i,j,k) = 0
c
c
c
        gyy(i,j,k) = 5.0000000000000d-1*(1.00000000000000d0 + o22 - o2*o
     &
     &  23)
c
c
c
c
c
c
c
c
c
c
c
        gyz(i,j,k) = 0
c
c
c
        gzz(i,j,k) = o22
c
c
c
c
c
c

            enddo
         enddo
      enddo
      
c
      kxx = 0.0d0
      kxy = 0.0d0
      kxz = 0.0d0
      kyy = 0.0d0
      kyz = 0.0d0
      kzz = 0.0d0
      
c
c
      i = ne-15
      adm = 0.0
      do k=2,np-1
         do j=2,nq-1
            adm=adm+(psisph(i,j,k)-(psisph(i+1,j,k)-psisph(i-1,j,k))/
     $           deta)*exp(0.5*etagrd(i))
         enddo
      enddo
      adm=adm/(nq-2)/(np-2)
      print *,'ADM mass: ',adm
      if (CCTK_Equals(initial_lapse,"schwarz")==1) then 
         write (*,*)"Initial with schwarzschild-like lapse"
         write (*,*)"using alp = (2.*r - adm)/(2.*r+adm)."
         alp = (2.*r - adm)/(2.*r+adm)
      endif
      
      deallocate(ac,ae,aw,an,as,aq,ab,rhs,
     $     qf,qfetaeta,qfqq,qfphi,qfphiphi,
     $     psisph,psiprim,detapsisph,dqpsisph,dphipsisph,detaetapsisph,
     $     detaqpsisph,detaphipsisph,dqqpsisph,dqphipsisph,
     $     dphiphipsisph)
      deallocate(etagrd,qgrd,phigrd)
      
      return
      end
