





























































































































 


































































































































































 






























































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                



















































































c

      subroutine RotatingDBHIVP(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 byJ
      integer, parameter :: cctki_use_byJ = kind(byJ)
      REAL*8 eta0
      integer, parameter :: cctki_use_eta0 = kind(eta0)
      REAL*8 etamax
      integer, parameter :: cctki_use_etamax = kind(etamax)
      REAL*8 mm
      integer, parameter :: cctki_use_mm = kind(mm)
      REAL*8 sigma
      integer, parameter :: cctki_use_sigma = kind(sigma)
      INTEGER*4 interpolation_order
      integer, parameter :: cctki_use_interpolation_order = kind(interpo
     &lation_order)
      INTEGER*4 kerr_shift
      integer, parameter :: cctki_use_kerr_shift = kind(kerr_shift)
      INTEGER*4 kerr_slice
      integer, parameter :: cctki_use_kerr_slice = kind(kerr_slice)
      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)
      INTEGER*4 verbose
      integer, parameter :: cctki_use_verbose = kind(verbose)
      COMMON /RotatingDBHIVPpriv/amp, byJ, eta0, etamax, mm, sigma, inte
     &rpolation_order, kerr_shift, kerr_slice, neta, nphi, ntheta, verbo
     &se
      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(:,:,:),Ksq(:,:,:),
     $     psisph(:,:,:),psip(:,:,:),psipp(:,:,:),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,o65,o66,o67,o68,o69,o70,o71,o72,o73,o74,o75,o76,o77,
     $     o78,o79,o80,o81,o82,o83,o84,o85,o86,o87,o88,o89,o90,o91,o92,
     $     o93,o94,o95,o96,o97,o98,o99,o100,o101,o102,o103,o104,o105,
     $     o106,o107,o108,o109,o110,o111,o112,o113,o114,o115,o116,o117,
     $     o118,o119,o120,o121,o122,o123,o124
      real*8 t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,
     $     t11,t12,t13,t14,t15,t16
      real*8 gtil,dngtil,dnngtil,dnnngtil,dnnnngtil,dnnnnngtil,
     $     dnnnnnngtil,dnnnnnnngtil
      real*8 rhsmax,rmax,adm
      real*8,parameter :: rbh_tol = 1.0d-7,rbh_eps = 1.0d-10
      integer,parameter :: itmax = 30 
      real*8 pi
      integer :: ne,nq,np
      integer :: nx,ny,nz
      integer i,j,k,it,ier,nquads,nocts,ntries
      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


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
c
      print *,"Brill wave + Rotating Distorted BH solve"
      write(*,123)amp,eta0,sigma,byJ,mm
      print*,'etamax=',etamax
 123  format(1x, 'Pars: amp',f8.5,' eta0',f8.5,' sigma',f8.5,' byJ',f9.5
     &
     $     ,' mm',f8.5)
      
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),
     $     Ksq(ne,nq,np),psisph(ne,nq,np),psip(ne,nq,np),
     $     psipp(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
      psip = 0.
      psipp = 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
      do k = 1,np
         do j = 1,nq
            do i = 1,ne

        gtil = amp*(exp(-((-eta0 + etagrd(i))**2/sigma**2)) + exp(-((eta
     &0 +
     &   etagrd(i))**2/sigma**2)))
        dngtil = amp*((-2.00000000000000d0*(-eta0 + etagrd(i))*exp(-((-e
     &ta0
     &   + etagrd(i))**2/sigma**2)))/sigma**2 - (2.00000000000000d0*(eta
     &0 +
     &   etagrd(i))*exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma**2)
        dnngtil = amp*((4.0000000000000d0*(-eta0 + etagrd(i))**2*exp(-((
     &-et
     &  a0 + etagrd(i))**2/sigma**2)))/sigma**4 - (2.00000000000000d0*ex
     &p(-
     &  ((-eta0 + etagrd(i))**2/sigma**2)))/sigma**2 + (4.0000000000000d
     &0*(
     &  eta0 + etagrd(i))**2*exp(-((eta0 + etagrd(i))**2/sigma**2)))/sig
     &ma**4 
     &  - (2.00000000000000d0*exp(-((eta0 + etagrd(i))**2/sigma**2)))/si
     &gma
     &  **2)
        dnnngtil = amp*((-8.0000000000000d0*(-eta0 + etagrd(i))**3*exp(-
     &((-
     &  eta0 + etagrd(i))**2/sigma**2)))/sigma**6 + (1.20000000000000d1*
     &(-e
     &  ta0 + etagrd(i))*exp(-((-eta0 + etagrd(i))**2/sigma**2)))/sigma*
     &*4 - (
     &  8.0000000000000d0*(eta0 + etagrd(i))**3*exp(-((eta0 + etagrd(i))
     &**2/si
     &  gma**2)))/sigma**6 + (1.20000000000000d1*(eta0 + etagrd(i))*exp(
     &-((
     &  eta0 + etagrd(i))**2/sigma**2)))/sigma**4)
        dnnnngtil = amp*((1.60000000000000d1*(-eta0 + etagrd(i))**4*exp(
     &-((
     &  -eta0 + etagrd(i))**2/sigma**2)))/sigma**8 - (4.8000000000000d1*
     &(-e
     &  ta0 + etagrd(i))**2*exp(-((-eta0 + etagrd(i))**2/sigma**2)))/sig
     &ma**6 
     &  + (1.20000000000000d1*exp(-((-eta0 + etagrd(i))**2/sigma**2)))/s
     &igm
     &  a**4 + (1.60000000000000d1*(eta0 + etagrd(i))**4*exp(-((eta0 + e
     &tag
     &  rd(i))**2/sigma**2)))/sigma**8 - (4.8000000000000d1*(eta0 + etag
     &rd(i))
     &  **2*exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma**6 + (1.200000
     &000
     &  00000d1*exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma**4)
        dnnnnngtil = amp*((-3.2000000000000d1*(-eta0 + etagrd(i))**5*exp
     &(-(
     &  (-eta0 + etagrd(i))**2/sigma**2)))/sigma**10 + (1.60000000000000
     &d2*
     &  (-eta0 + etagrd(i))**3*exp(-((-eta0 + etagrd(i))**2/sigma**2)))/
     &sigma*
     &  *8 - (1.20000000000000d2*(-eta0 + etagrd(i))*exp(-((-eta0 + etag
     &rd(i))
     &  **2/sigma**2)))/sigma**6 - (3.2000000000000d1*(eta0 + etagrd(i))
     &**5
     &  *exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma**10 + (1.60000000
     &000
     &  000d2*(eta0 + etagrd(i))**3*exp(-((eta0 + etagrd(i))**2/sigma**2
     &)))/si
     &  gma**8 - (1.20000000000000d2*(eta0 + etagrd(i))*exp(-((eta0 + et
     &agr
     &  d(i))**2/sigma**2)))/sigma**6)
        dnnnnnngtil = amp*((6.4000000000000d1*(-eta0 + etagrd(i))**6*exp
     &(-(
     &  (-eta0 + etagrd(i))**2/sigma**2)))/sigma**12 - (4.8000000000000d
     &2*(
     &  -eta0 + etagrd(i))**4*exp(-((-eta0 + etagrd(i))**2/sigma**2)))/s
     &igma**
     &  10 + (7.2000000000000d2*(-eta0 + etagrd(i))**2*exp(-((-eta0 + et
     &agr
     &  d(i))**2/sigma**2)))/sigma**8 - (1.20000000000000d2*exp(-((-eta0
     & + 
     &  etagrd(i))**2/sigma**2)))/sigma**6 + (6.4000000000000d1*(eta0 + 
     &eta
     &  grd(i))**6*exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma**12 - (
     &4.8000
     &  000000000d2*(eta0 + etagrd(i))**4*exp(-((eta0 + etagrd(i))**2/si
     &gma**2
     &  )))/sigma**10 + (7.2000000000000d2*(eta0 + etagrd(i))**2*exp(-((
     &eta
     &  0 + etagrd(i))**2/sigma**2)))/sigma**8 - (1.20000000000000d2*exp
     &(-(
     &  (eta0 + etagrd(i))**2/sigma**2)))/sigma**6)
        dnnnnnnngtil = amp*((-1.28000000000000d2*(-eta0 + etagrd(i))**7*
     &exp
     &  (-((-eta0 + etagrd(i))**2/sigma**2)))/sigma**14 + (1.34400000000
     &000
     &  d3*(-eta0 + etagrd(i))**5*exp(-((-eta0 + etagrd(i))**2/sigma**2)
     &))/sig
     &  ma**12 - (3.3600000000000d3*(-eta0 + etagrd(i))**3*exp(-((-eta0 
     &+ e
     &  tagrd(i))**2/sigma**2)))/sigma**10 + (1.68000000000000d3*(-eta0 
     &+ e
     &  tagrd(i))*exp(-((-eta0 + etagrd(i))**2/sigma**2)))/sigma**8 - (1
     &.28000
     &  000000000d2*(eta0 + etagrd(i))**7*exp(-((eta0 + etagrd(i))**2/si
     &gma**2
     &  )))/sigma**14 + (1.34400000000000d3*(eta0 + etagrd(i))**5*exp(-(
     &(et
     &  a0 + etagrd(i))**2/sigma**2)))/sigma**12 - (3.3600000000000d3*(e
     &ta0
     &   + etagrd(i))**3*exp(-((eta0 + etagrd(i))**2/sigma**2)))/sigma**
     &10 + (
     &  1.68000000000000d3*(eta0 + etagrd(i))*exp(-((eta0 + etagrd(i))**
     &2/sigm
     &  a**2)))/sigma**8)

        o1 = mm**2
        o2 = 5.0000000000000d1*o1
        o3 = -8.0000000000000d1 + o2
        o4 = 1/o3
        o5 = 2.00000000000000d0*qgrd(j)
        o6 = cos(o5)
        o7 = 5.0000000000000d0*o6
        o8 = 3.00000000000000d0 + o7
        o9 = dnnnnngtil
        o10 = mm*phigrd(k)
        o11 = sin(o10)
        o12 = sin(qgrd(j))
        o13 = o12**2
        o14 = o13**2
        o15 = 2.50000000000000d1*o1
        o16 = -4.0000000000000d1 + o15
        o17 = 1/o16
        o18 = dnnngtil
        o19 = 2.00000000000000d1*o1
        o20 = -3.2000000000000d1 + o19
        o21 = 1/o20
        o22 = dngtil
        o23 = 5.0000000000000d0*o1
        o24 = -8.0000000000000d0 + o23
        o25 = 1/o24
        o26 = cos(o10)
        o27 = cos(qgrd(j))
        o28 = 1.03000000000000d2*o1
        o29 = -9.2800000000000d2 + o28
        o30 = o12*o13
        o31 = 3.5000000000000d1*o1
        o32 = -8.3000000000000d1 + o31
        o33 = -1.60000000000000d1 + o1
        o34 = dnnnnnnngtil
        o35 = 4.0000000000000d0*qgrd(j)
        o36 = cos(o35)
        o37 = dnnnngtil
        o38 = 1.00000000000000d1*o1
        o39 = -1.60000000000000d1 + o38
        o40 = 1/o39
        o41 = dnngtil
        o42 = dnnnnnngtil
        o43 = o27**2
        o44 = gtil
        o45 = o13*o14
        o46 = 1.07000000000000d2*o1
        o47 = 5.3000000000000d1*o1
        o48 = sin(o5)
        o49 = byJ**2
        o50 = o24**2
        o51 = 1/o50
        o52 = o26**2
        o53 = o37**2
        o54 = o1**2
        o55 = o6**2
        o56 = o36**2
        o57 = o39**2
        o58 = 1/o57
        o59 = o41**2
        o60 = o9**2
        o61 = o16**2
        o62 = 1/o61
        o63 = o18**2
        o64 = o20**2
        o65 = 1/o64
        o66 = o22**2
        o67 = o11**2
        o68 = o1*o54
        o69 = o3**2
        o70 = 1/o69
        o71 = o42**2
        o72 = o34**2
        o73 = o14**2
        o74 = o48**2
        o75 = o43**2
        o76 = o44**2
        t1 = -1.00800000000000d2*byJ*o13*o25*o26*o37 - 1.39500000000000d
     &1*
     &  byJ*o1*o13*o25*o26*o37 - 9.7440000000000d2*byJ*o13*o25*o26*o36*o
     &37 +
     &   1.08150000000000d2*byJ*o1*o13*o25*o26*o36*o37 + 1.3320000000000
     &0d
     &  3*byJ*o13*o26*o40*o41 - 4.6800000000000d2*byJ*o1*o13*o26*o40*o41
     & + 6
     &  .9720000000000d3*byJ*o13*o26*o36*o40*o41 - 2.94000000000000d3*by
     &J*o1
     &  *o13*o26*o36*o40*o41 - 9.6000000000000d2*byJ*o14*o26*o4*o42 + 6.
     &00
     &  00000000000d1*byJ*o1*o14*o26*o4*o42 - 1.80000000000000d2*byJ*o14
     &*o26
     &  *o44 + 7.2000000000000d2*byJ*o13*o26*o43*o44 - 2.70000000000000d
     &3*
     &  byJ*o14*o26*o43*o44 + 4.5000000000000d2*byJ*o26*o44*o45 + 1.8000
     &0000
     &  000000d1*o13*o49 - 3.7296000000000d3*o13*o25*o37*o40*o41*o52 + 7
     &
     &  .9425000000000d2*o1*o13*o25*o37*o40*o41*o52 - 5.5574400000000d4*
     &
     &  o13*o25*o36*o37*o40*o41*o52 + 2.21991000000000d4*o1*o13*o25*o36*
     &
     &  o37*o40*o41*o52 + 2.68800000000000d3*o14*o25*o37*o4*o42*o52 + 2.
     &
     &  04000000000000d2*o1*o14*o25*o37*o4*o42*o52 + 2.59840000000000d4*
     &
     &  o14*o25*o36*o37*o4*o42*o52 - 4.5080000000000d3*o1*o14*o25*o36*o3
     &
     &  7*o4*o42*o52 - 3.5520000000000d4*o14*o4*o40*o41*o42*o52 + 1.4700
     &
     &  0000000000d4*o1*o14*o4*o40*o41*o42*o52 - 1.85920000000000d5*o14*
     &
     &  o36*o4*o40*o41*o42*o52 + 9.0020000000000d4*o1*o14*o36*o4*o40*o41
     &
     &  *o42*o52
        t2 = -1.10400000000000d5*o14*o17*o18*o21*o22*o43*o52 + 2.9160000
     &
     &  0000000d4*o1*o14*o17*o18*o21*o22*o43*o52 + 5.0400000000000d2*o14
     &
     &  *o25*o37*o44*o52 + 6.9750000000000d1*o1*o14*o25*o37*o44*o52 + 4.
     &
     &  8720000000000d3*o14*o25*o36*o37*o44*o52 - 5.4075000000000d2*o1*o
     &
     &  14*o25*o36*o37*o44*o52 - 6.6600000000000d3*o14*o40*o41*o44*o52 +
     &
     &   2.34000000000000d3*o1*o14*o40*o41*o44*o52 - 3.4860000000000d4*o
     &
     &  14*o36*o40*o41*o44*o52 + 1.47000000000000d4*o1*o14*o36*o40*o41*o
     &
     &  44*o52 - 2.01600000000000d3*o13*o25*o37*o43*o44*o52 - 2.79000000
     &
     &  000000d2*o1*o13*o25*o37*o43*o44*o52 + 7.5600000000000d3*o14*o25*
     &
     &  o37*o43*o44*o52 + 1.04625000000000d3*o1*o14*o25*o37*o43*o44*o52 
     &
     &  - 1.94880000000000d4*o13*o25*o36*o37*o43*o44*o52 + 2.16300000000
     &
     &  000d3*o1*o13*o25*o36*o37*o43*o44*o52 + 7.3080000000000d4*o14*o25
     &
     &  *o36*o37*o43*o44*o52 - 8.1112500000000d3*o1*o14*o25*o36*o37*o43*
     &
     &  o44*o52 + 2.66400000000000d4*o13*o40*o41*o43*o44*o52 - 9.3600000
     &
     &  000000d3*o1*o13*o40*o41*o43*o44*o52 - 9.9900000000000d4*o14*o40*
     &
     &  o41*o43*o44*o52 + 3.5100000000000d4*o1*o14*o40*o41*o43*o44*o52 +
     &
     &   1.39440000000000d5*o13*o36*o40*o41*o43*o44*o52 - 5.880000000000
     &
     &  0d4*o1*o13*o36*o40*o41*o43*o44*o52 - 5.2290000000000d5*o14*o36*o
     &
     &  40*o41*o43*o44*o52 + 2.20500000000000d5*o1*o14*o36*o40*o41*o43*o
     &
     &  44*o52 - 1.92000000000000d4*o14*o4*o42*o43*o44*o52
        t3 = 1.20000000000000d3*o1*o14*o4*o42*o43*o44*o52 - 7.6800000000
     &
     &  000d3*o17*o21*o22*o34*o43*o45*o52 + 1.44000000000000d3*o1*o17*o2
     &
     &  1*o22*o34*o43*o45*o52 - 1.26000000000000d3*o25*o37*o44*o45*o52 -
     &
     &   1.74375000000000d2*o1*o25*o37*o44*o45*o52 - 1.21800000000000d4*
     &
     &  o25*o36*o37*o44*o45*o52 + 1.35187500000000d3*o1*o25*o36*o37*o44*
     &
     &  o45*o52 + 1.66500000000000d4*o40*o41*o44*o45*o52 - 5.85000000000
     &
     &  00d3*o1*o40*o41*o44*o45*o52 + 8.7150000000000d4*o36*o40*o41*o44*
     &
     &  o45*o52 - 3.6750000000000d4*o1*o36*o40*o41*o44*o45*o52 + 4.80000
     &
     &  00000000d3*o4*o42*o44*o45*o52 - 3.00000000000000d2*o1*o4*o42*o44
     &
     &  *o45*o52 + 7.2000000000000d4*o4*o42*o43*o44*o45*o52 - 4.50000000
     &
     &  00000d3*o1*o4*o42*o43*o44*o45*o52 + 1.41120000000000d2*o13*o51*o
     &
     &  52*o53 + 3.9060000000000d1*o1*o13*o51*o52*o53 + 2.72832000000000
     &
     &  d3*o13*o36*o51*o52*o53 + 7.4760000000000d1*o1*o13*o36*o51*o52*o5
     &
     &  3 + 1.81350000000000d2*o13*o25*o37*o40*o41*o52*o54 - 2.667000000
     &
     &  00000d2*o13*o25*o36*o37*o40*o41*o52*o54 - 2.32500000000000d1*o14
     &
     &  *o25*o37*o4*o42*o52*o54 + 1.80250000000000d2*o14*o25*o36*o37*o4*
     &
     &  o42*o52*o54 - 7.8000000000000d2*o14*o4*o40*o41*o42*o52*o54 - 4.9
     &
     &  000000000000d3*o14*o36*o4*o40*o41*o42*o52*o54 - 1.92000000000000
     &
     &  d3*o14*o17*o18*o21*o22*o43*o52*o54 - 6.0000000000000d1*o17*o21*o
     &
     &  22*o34*o43*o45*o52*o54
        t4 = 2.70281250000000d0*o13*o51*o52*o53*o54 - 4.1908125000000d1*
     &
     &  o13*o36*o51*o52*o53*o54 - 2.74432000000000d4*o13*o25*o37*o40*o41
     &
     &  *o52*o55 + 9.5432000000000d3*o1*o13*o25*o37*o40*o41*o52*o55 + 1.
     &
     &  14688000000000d4*o14*o25*o37*o4*o42*o52*o55 - 4.2560000000000d2*
     &
     &  o1*o14*o25*o37*o4*o42*o52*o55 - 1.20064000000000d5*o14*o4*o40*o4
     &
     &  1*o42*o52*o55 + 5.2304000000000d4*o1*o14*o4*o40*o41*o42*o52*o55 
     &
     &  - 1.99200000000000d6*o14*o17*o18*o21*o22*o43*o52*o55 + 2.0850000
     &
     &  0000000d6*o1*o14*o17*o18*o21*o22*o43*o52*o55 + 1.31072000000000d
     &
     &  3*o13*o51*o52*o53*o55 + 6.6560000000000d1*o1*o13*o51*o52*o53*o55
     &
     &   + 2.60000000000000d2*o13*o25*o37*o40*o41*o52*o54*o55 - 1.820000
     &
     &  00000000d1*o14*o25*o37*o4*o42*o52*o54*o55 - 2.80000000000000d3*o
     &
     &  14*o4*o40*o41*o42*o52*o54*o55 - 5.2500000000000d5*o14*o17*o18*o2
     &
     &  1*o22*o43*o52*o54*o55 + 8.4500000000000d-1*o13*o51*o52*o53*o54*o
     &
     &  55 - 1.88708800000000d5*o13*o25*o37*o40*o41*o52*o56 + 1.00521050
     &
     &  000000d5*o1*o13*o25*o37*o40*o41*o52*o56 + 1.31868800000000d4*o13
     &
     &  *o51*o52*o53*o56 - 2.92726000000000d3*o1*o13*o51*o52*o53*o56 - 8
     &
     &  .8322500000000d3*o13*o25*o37*o40*o41*o52*o54*o56 + 1.62450312500
     &
     &  000d2*o13*o51*o52*o53*o54*o56 + 2.46420000000000d4*o13*o52*o58*o
     &
     &  59 - 1.73160000000000d4*o1*o13*o52*o58*o59 + 2.57964000000000d5*
     &
     &  o13*o36*o52*o58*o59 - 1.99416000000000d5*o1*o13*o36*o52*o58*o59
        t5 = 3.04200000000000d3*o13*o52*o54*o58*o59 + 3.8220000000000d4*
     &
     &  o13*o36*o52*o54*o58*o59 + 1.43648000000000d5*o13*o52*o55*o58*o59
     &
     &   - 1.07200000000000d5*o1*o13*o52*o55*o58*o59 + 2.00000000000000d
     &
     &  4*o13*o52*o54*o55*o58*o59 + 6.7512200000000d5*o13*o52*o56*o58*o5
     &
     &  9 - 5.6938000000000d5*o1*o13*o52*o56*o58*o59 + 1.20050000000000d
     &
     &  5*o13*o52*o54*o56*o58*o59 - 3.07200000000000d2*byJ*o13*o25*o26*o
     &37
     &  *o6 - 7.8000000000000d0*byJ*o1*o13*o25*o26*o37*o6 + 3.2160000000
     &00
     &  0d3*byJ*o13*o26*o40*o41*o6 - 1.20000000000000d3*byJ*o1*o13*o26*o
     &40*o
     &  41*o6 - 1.34400000000000d3*byJ*o14*o26*o4*o42*o6 + 8.40000000000
     &00
     &  d1*byJ*o1*o14*o26*o4*o42*o6 - 2.03712000000000d4*o13*o25*o37*o40
     &*o
     &  41*o52*o6 + 5.8188000000000d3*o1*o13*o25*o37*o40*o41*o52*o6 - 1.
     &
     &  46540800000000d5*o13*o25*o36*o37*o40*o41*o52*o6 + 6.571880000000
     &
     &  0d4*o1*o13*o25*o36*o37*o40*o41*o52*o6 + 1.19552000000000d4*o14*o
     &
     &  25*o37*o4*o42*o52*o6 - 1.84000000000000d1*o1*o14*o25*o37*o4*o42*
     &
     &  o52*o6 + 3.6377600000000d4*o14*o25*o36*o37*o4*o42*o52*o6 - 6.311
     &
     &  2000000000d3*o1*o14*o25*o36*o37*o4*o42*o52*o6 - 1.35488000000000
     &
     &  d5*o14*o4*o40*o41*o42*o52*o6 + 5.7940000000000d4*o1*o14*o4*o40*o
     &
     &  41*o42*o52*o6 - 2.60288000000000d5*o14*o36*o4*o40*o41*o42*o52*o6
     &
     &   + 1.26028000000000d5*o1*o14*o36*o4*o40*o41*o42*o52*o6 + 9.50400
     &
     &  00000000d5*o14*o17*o18*o21*o22*o43*o52*o6
        t6 = -6.3960000000000d5*o1*o14*o17*o18*o21*o22*o43*o52*o6 + 1.53
     &
     &  600000000000d3*o14*o25*o37*o44*o52*o6 + 3.9000000000000d1*o1*o14
     &
     &  *o25*o37*o44*o52*o6 - 1.60800000000000d4*o14*o40*o41*o44*o52*o6 
     &
     &  + 6.0000000000000d3*o1*o14*o40*o41*o44*o52*o6 - 6.1440000000000d
     &
     &  3*o13*o25*o37*o43*o44*o52*o6 - 1.56000000000000d2*o1*o13*o25*o37
     &
     &  *o43*o44*o52*o6 + 2.30400000000000d4*o14*o25*o37*o43*o44*o52*o6 
     &
     &  + 5.8500000000000d2*o1*o14*o25*o37*o43*o44*o52*o6 + 6.4320000000
     &
     &  000d4*o13*o40*o41*o43*o44*o52*o6 - 2.40000000000000d4*o1*o13*o40
     &
     &  *o41*o43*o44*o52*o6 - 2.41200000000000d5*o14*o40*o41*o43*o44*o52
     &
     &  *o6 + 9.0000000000000d4*o1*o14*o40*o41*o43*o44*o52*o6 - 2.688000
     &
     &  00000000d4*o14*o4*o42*o43*o44*o52*o6 + 1.68000000000000d3*o1*o14
     &
     &  *o4*o42*o43*o44*o52*o6 + 3.8400000000000d4*o17*o21*o22*o34*o43*o
     &
     &  45*o52*o6 - 2.64000000000000d4*o1*o17*o21*o22*o34*o43*o45*o52*o6
     &
     &   - 3.8400000000000d3*o25*o37*o44*o45*o52*o6 - 9.7500000000000d1*
     &
     &  o1*o25*o37*o44*o45*o52*o6 + 4.0200000000000d4*o40*o41*o44*o45*o5
     &
     &  2*o6 - 1.50000000000000d4*o1*o40*o41*o44*o45*o52*o6 + 6.72000000
     &
     &  00000d3*o4*o42*o44*o45*o52*o6 - 4.2000000000000d2*o1*o4*o42*o44*
     &
     &  o45*o52*o6 + 1.00800000000000d5*o4*o42*o43*o44*o45*o52*o6 - 6.30
     &
     &  00000000000d3*o1*o4*o42*o43*o44*o45*o52*o6 + 8.6016000000000d2*o
     &
     &  13*o51*o52*o53*o6 + 1.40880000000000d2*o1*o13*o51*o52*o53*o6
        t7 = 8.3148800000000d3*o13*o36*o51*o52*o53*o6 - 7.1176000000000d
     &
     &  2*o1*o13*o36*o51*o52*o53*o6 + 5.6640000000000d2*o13*o25*o37*o40*
     &
     &  o41*o52*o54*o6 - 2.96800000000000d3*o13*o25*o36*o37*o40*o41*o52*
     &
     &  o54*o6 - 4.5550000000000d1*o14*o25*o37*o4*o42*o52*o54*o6 + 2.523
     &
     &  50000000000d2*o14*o25*o36*o37*o4*o42*o52*o54*o6 - 3.092000000000
     &
     &  00d3*o14*o4*o40*o41*o42*o52*o54*o6 - 6.8600000000000d3*o14*o36*o
     &
     &  4*o40*o41*o42*o52*o54*o6 + 6.9000000000000d4*o14*o17*o18*o21*o22
     &
     &  *o43*o52*o54*o6 + 1.50000000000000d3*o17*o21*o22*o34*o43*o45*o52
     &
     &  *o54*o6 + 3.02250000000000d0*o13*o51*o52*o53*o54*o6 - 2.34325000
     &
     &  000000d1*o13*o36*o51*o52*o53*o54*o6 + 1.18992000000000d5*o13*o52
     &
     &  *o58*o59*o6 - 8.6208000000000d4*o1*o13*o52*o58*o59*o6 + 6.228320
     &
     &  0000000d5*o13*o36*o52*o58*o59*o6 - 4.9504000000000d5*o1*o13*o36*
     &
     &  o52*o58*o59*o6 + 1.56000000000000d4*o13*o52*o54*o58*o59*o6 + 9.8
     &
     &  000000000000d4*o13*o36*o52*o54*o58*o59*o6 + 6.1952000000000d2*o1
     &
     &  4*o43*o51*o52*o60 - 1.93600000000000d2*o1*o14*o43*o51*o52*o60 + 
     &
     &  1.51250000000000d1*o14*o43*o51*o52*o54*o60 + 4.3059200000000d3*o
     &
     &  14*o43*o51*o52*o55*o60 - 9.5584000000000d2*o1*o14*o43*o51*o52*o5
     &
     &  5*o60 + 5.3045000000000d1*o14*o43*o51*o52*o54*o55*o60 - 3.266560
     &
     &  0000000d3*o14*o43*o51*o52*o6*o60 + 8.7296000000000d2*o1*o14*o43*
     &
     &  o51*o52*o6*o60 - 5.6650000000000d1*o14*o43*o51*o52*o54*o6*o60
        t8 = 1.47200000000000d4*o18*o34*o43*o45*o52*o62 - 2.968000000000
     &
     &  00d3*o1*o18*o34*o43*o45*o52*o62 + 1.28000000000000d2*o18*o34*o43
     &
     &  *o45*o52*o54*o62 - 5.3120000000000d4*o18*o34*o43*o45*o52*o6*o62 
     &
     &  + 2.57200000000000d4*o1*o18*o34*o43*o45*o52*o6*o62 - 1.400000000
     &
     &  00000d3*o18*o34*o43*o45*o52*o54*o6*o62 + 1.05800000000000d5*o14*
     &
     &  o43*o52*o62*o63 - 2.94400000000000d4*o1*o14*o43*o52*o62*o63 + 2.
     &
     &  04800000000000d3*o14*o43*o52*o54*o62*o63 + 1.37780000000000d6*o1
     &
     &  4*o43*o52*o55*o62*o63 - 1.16200000000000d6*o1*o14*o43*o52*o55*o6
     &
     &  2*o63 + 2.45000000000000d5*o14*o43*o52*o54*o55*o62*o63 - 7.63600
     &
     &  00000000d5*o14*o43*o52*o6*o62*o63 + 4.2824000000000d5*o1*o14*o43
     &
     &  *o52*o6*o62*o63 - 4.4800000000000d4*o14*o43*o52*o54*o6*o62*o63 +
     &
     &   2.88000000000000d4*o14*o43*o52*o65*o66 - 7.2000000000000d3*o1*o
     &
     &  14*o43*o52*o65*o66 + 4.5000000000000d2*o14*o43*o52*o54*o65*o66 +
     &
     &   7.2000000000000d5*o14*o43*o52*o55*o65*o66 - 9.0000000000000d5*o
     &
     &  1*o14*o43*o52*o55*o65*o66 + 2.81250000000000d5*o14*o43*o52*o54*o
     &
     &  55*o65*o66 - 2.88000000000000d5*o14*o43*o52*o6*o65*o66 + 2.16000
     &
     &  000000000d5*o1*o14*o43*o52*o6*o65*o66 - 2.25000000000000d4*o14*o
     &
     &  43*o52*o54*o6*o65*o66 - 2.10600000000000d4*o1*o14*o17*o18*o21*o2
     &
     &  2*o67 - 6.7562000000000d4*o1*o14*o17*o18*o22*o25*o67 - 7.0180000
     &
     &  000000d4*o1*o14*o17*o18*o22*o40*o67 + 1.92000000000000d3*o1*o14*
     &
     &  o17*o42*o43*o44*o67
        t9 = -3.3600000000000d3*o1*o17*o42*o43*o44*o45*o67 - 5.440000000
     &
     &  0000d2*o1*o17*o25*o27*o30*o37*o42*o48*o67 + 1.98400000000000d3*o
     &
     &  1*o17*o27*o30*o40*o41*o42*o48*o67 - 1.02000000000000d3*o1*o12*o2
     &
     &  5*o27*o37*o44*o48*o67 + 1.78500000000000d3*o1*o25*o27*o30*o37*o4
     &
     &  4*o48*o67 + 3.7200000000000d3*o1*o12*o27*o40*o41*o44*o48*o67 - 6
     &
     &  .5100000000000d3*o1*o27*o30*o40*o41*o44*o48*o67 + 4.625200000000
     &
     &  0d4*o14*o17*o18*o22*o25*o54*o67 + 6.6494000000000d4*o14*o17*o18*
     &
     &  o22*o40*o54*o67 - 1.20000000000000d2*o14*o17*o42*o43*o44*o54*o67
     &
     &   + 2.10000000000000d2*o17*o42*o43*o44*o45*o54*o67 + 1.2200000000
     &
     &  0000d2*o17*o25*o27*o30*o37*o42*o48*o54*o67 - 5.7200000000000d2*o
     &
     &  17*o27*o30*o40*o41*o42*o48*o54*o67 + 1.65000000000000d2*o12*o25*
     &
     &  o27*o37*o44*o48*o54*o67 - 2.88750000000000d2*o25*o27*o30*o37*o44
     &
     &  *o48*o54*o67 - 8.4000000000000d2*o12*o27*o40*o41*o44*o48*o54*o67
     &
     &   + 1.47000000000000d3*o27*o30*o40*o41*o44*o48*o54*o67 - 5.850000
     &
     &  0000000d4*o1*o14*o17*o18*o21*o22*o55*o67 - 1.42506000000000d5*o1
     &
     &  *o14*o17*o18*o22*o25*o55*o67 - 1.31652000000000d5*o1*o14*o17*o18
     &
     &  *o22*o40*o55*o67 + 9.2016000000000d4*o14*o17*o18*o22*o25*o54*o55
     &
     &  *o67 + 1.25082000000000d5*o14*o17*o18*o22*o40*o54*o55*o67 - 7.02
     &
     &  00000000000d4*o1*o14*o17*o18*o21*o22*o6*o67 - 1.96332000000000d5
     &
     &  *o1*o14*o17*o18*o22*o25*o6*o67 - 1.92408000000000d5*o1*o14*o17*o
     &
     &  18*o22*o40*o6*o67 - 1.29920000000000d3*o1*o17*o25*o27*o30*o37*o4
     &
     &  2*o48*o6*o67 + 1.33760000000000d4*o1*o17*o27*o30*o40*o41*o42*o48
     &
     &  *o6*o67
        t10 = -2.43600000000000d3*o1*o12*o25*o27*o37*o44*o48*o6*o67 + 4.
     &
     &  2630000000000d3*o1*o25*o27*o30*o37*o44*o48*o6*o67 + 2.5080000000
     &
     &  0000d4*o1*o12*o27*o40*o41*o44*o48*o6*o67 - 4.3890000000000d4*o1*
     &
     &  o27*o30*o40*o41*o44*o48*o6*o67 + 1.30692000000000d5*o14*o17*o18*
     &
     &  o22*o25*o54*o6*o67 + 1.82664000000000d5*o14*o17*o18*o22*o40*o54*
     &
     &  o6*o67 + 1.08400000000000d2*o17*o25*o27*o30*o37*o42*o48*o54*o6*o
     &
     &  67 - 6.5320000000000d3*o17*o27*o30*o40*o41*o42*o48*o54*o6*o67 + 
     &
     &  5.1000000000000d1*o12*o25*o27*o37*o44*o48*o54*o6*o67 - 8.9250000
     &
     &  000000d1*o25*o27*o30*o37*o44*o48*o54*o6*o67 - 1.06800000000000d4
     &
     &  *o12*o27*o40*o41*o44*o48*o54*o6*o67 + 1.86900000000000d4*o27*o30
     &
     &  *o40*o41*o44*o48*o54*o6*o67 + 4.0960000000000d3*o1*o14*o60*o62*o
     &
     &  67 - 5.1200000000000d2*o14*o54*o60*o62*o67 + 9.2160000000000d3*o
     &
     &  1*o14*o55*o60*o62*o67 - 1.15200000000000d3*o14*o54*o55*o60*o62*o
     &
     &  67 + 1.22880000000000d4*o1*o14*o6*o60*o62*o67 - 1.53600000000000
     &
     &  d3*o14*o54*o6*o60*o62*o67 + 2.63438000000000d5*o1*o14*o62*o63*o6
     &
     &  7 - 1.49158000000000d5*o14*o54*o62*o63*o67 + 5.8030200000000d5*o
     &
     &  1*o14*o55*o62*o63*o67 - 3.2531400000000d5*o14*o54*o55*o62*o63*o6
     &
     &  7 + 7.8147600000000d5*o1*o14*o6*o62*o63*o67 - 4.4056800000000d5*
     &
     &  o14*o54*o6*o62*o63*o67 + 6.8890000000000d3*o1*o14*o51*o66*o67 - 
     &
     &  5.8100000000000d3*o14*o51*o54*o66*o67 + 1.36890000000000d4*o1*o1
     &
     &  4*o51*o55*o66*o67
        t11 = -1.05300000000000d4*o14*o51*o54*o55*o66*o67 + 1.4641000000
     &
     &  0000d4*o1*o14*o58*o66*o67 - 1.69400000000000d4*o14*o54*o58*o66*o
     &
     &  67 + 2.52810000000000d4*o1*o14*o55*o58*o66*o67 - 2.8620000000000
     &
     &  0d4*o14*o54*o55*o58*o66*o67 + 1.94220000000000d4*o1*o14*o51*o6*o
     &
     &  66*o67 - 1.56600000000000d4*o14*o51*o54*o6*o66*o67 + 3.847800000
     &
     &  0000d4*o1*o14*o58*o6*o66*o67 - 4.4040000000000d4*o14*o54*o58*o6*
     &
     &  o66*o67 + 8.1000000000000d3*o1*o14*o65*o66*o67 + 2.2500000000000
     &
     &  0d4*o1*o14*o55*o65*o66*o67 + 2.70000000000000d4*o1*o14*o6*o65*o6
     &
     &  6*o67 - 7.4900000000000d3*o14*o17*o18*o22*o25*o67*o68 - 1.498000
     &
     &  00000000d4*o14*o17*o18*o22*o40*o67*o68 - 5.5000000000000d0*o17*o
     &
     &  25*o27*o30*o37*o42*o48*o67*o68 + 2.80000000000000d1*o17*o27*o30*
     &
     &  o40*o41*o42*o48*o67*o68 - 1.43100000000000d4*o14*o17*o18*o22*o25
     &
     &  *o55*o67*o68 - 2.86200000000000d4*o14*o17*o18*o22*o40*o55*o67*o6
     &
     &  8 - 2.07600000000000d4*o14*o17*o18*o22*o25*o6*o67*o68 - 4.152000
     &
     &  0000000d4*o14*o17*o18*o22*o40*o6*o67*o68 - 1.70000000000000d0*o1
     &
     &  7*o25*o27*o30*o37*o42*o48*o6*o67*o68 + 3.5600000000000d2*o17*o27
     &
     &  *o30*o40*o41*o42*o48*o6*o67*o68 + 1.60000000000000d1*o14*o60*o62
     &
     &  *o67*o68 + 3.6000000000000d1*o14*o55*o60*o62*o67*o68 + 4.8000000
     &
     &  000000d1*o14*o6*o60*o62*o67*o68 + 2.28980000000000d4*o14*o62*o63
     &
     &  *o67*o68 + 5.0562000000000d4*o14*o55*o62*o63*o67*o68
        t12 = 6.8052000000000d4*o14*o6*o62*o63*o67*o68 + 1.2250000000000
     &
     &  0d3*o14*o51*o66*o67*o68 + 2.02500000000000d3*o14*o51*o55*o66*o67
     &
     &  *o68 + 4.9000000000000d3*o14*o58*o66*o67*o68 + 8.1000000000000d3
     &
     &  *o14*o55*o58*o66*o67*o68 + 3.15000000000000d3*o14*o51*o6*o66*o67
     &
     &  *o68 + 1.26000000000000d4*o14*o58*o6*o66*o67*o68 + 1.42420000000
     &
     &  000d4*o1*o14*o60*o67*o70 - 1.90400000000000d3*o14*o54*o60*o67*o7
     &
     &  0 + 3.15540000000000d4*o1*o14*o55*o60*o67*o70 - 4.2480000000000d
     &
     &  3*o14*o54*o55*o60*o67*o70 + 4.2396000000000d4*o1*o14*o6*o60*o67*
     &
     &  o70 - 5.6880000000000d3*o14*o54*o6*o60*o67*o70 + 6.4000000000000
     &
     &  d1*o14*o60*o67*o68*o70 + 1.44000000000000d2*o14*o55*o60*o67*o68*
     &
     &  o70 + 1.92000000000000d2*o14*o6*o60*o67*o68*o70 + 5.120000000000
     &
     &  0d2*o1*o43*o45*o62*o67*o71 - 6.4000000000000d1*o43*o45*o54*o62*o
     &
     &  67*o71 + 2.00000000000000d0*o43*o45*o62*o67*o68*o71 + 1.28000000
     &
     &  000000d4*o45*o52*o70*o71 - 1.60000000000000d3*o1*o45*o52*o70*o71
     &
     &   + 5.0000000000000d1*o45*o52*o54*o70*o71 + 2.50880000000000d4*o4
     &
     &  5*o52*o55*o70*o71 - 3.13600000000000d3*o1*o45*o52*o55*o70*o71 + 
     &
     &  9.8000000000000d1*o45*o52*o54*o55*o70*o71 + 3.5840000000000d4*o4
     &
     &  5*o52*o6*o70*o71 - 4.4800000000000d3*o1*o45*o52*o6*o70*o71 + 1.4
     &
     &  0000000000000d2*o45*o52*o54*o6*o70*o71
        t13 = -1.20000000000000d4*o4*o42*o44*o52*o73 + 7.5000000000000d2
     &
     &  *o1*o4*o42*o44*o52*o73 - 1.68000000000000d4*o4*o42*o44*o52*o6*o7
     &
     &  3 + 1.05000000000000d3*o1*o4*o42*o44*o52*o6*o73 + 5.120000000000
     &
     &  0d2*o43*o52*o62*o72*o73 - 6.4000000000000d1*o1*o43*o52*o62*o72*o
     &
     &  73 + 2.00000000000000d0*o43*o52*o54*o62*o72*o73 - 1.054000000000
     &
     &  00d3*o1*o25*o37*o40*o41*o67*o74 + 1.44500000000000d2*o1*o51*o53*
     &
     &  o67*o74 + 4.0850000000000d2*o25*o37*o40*o41*o54*o67*o74 - 4.6750
     &
     &  000000000d1*o51*o53*o54*o67*o74 - 1.69708000000000d4*o1*o25*o37*
     &
     &  o40*o41*o55*o67*o74 + 8.2418000000000d2*o1*o51*o53*o55*o67*o74 +
     &
     &   7.5821000000000d3*o25*o37*o40*o41*o54*o55*o67*o74 - 3.451000000
     &
     &  0000d1*o51*o53*o54*o55*o67*o74 + 1.92200000000000d3*o1*o58*o59*o
     &
     &  67*o74 - 8.6800000000000d2*o54*o58*o59*o67*o74 + 8.7362000000000
     &
     &  d4*o1*o55*o58*o59*o67*o74 - 7.4404000000000d4*o54*o55*o58*o59*o6
     &
     &  7*o74 - 9.6232000000000d3*o1*o25*o37*o40*o41*o6*o67*o74 + 6.9020
     &
     &  000000000d2*o1*o51*o53*o6*o67*o74 + 4.7966000000000d3*o25*o37*o4
     &
     &  0*o41*o54*o6*o67*o74 - 1.26100000000000d2*o51*o53*o54*o6*o67*o74
     &
     &   + 2.59160000000000d4*o1*o58*o59*o6*o67*o74 - 1.68880000000000d4
     &
     &  *o54*o58*o59*o6*o67*o74 - 3.8500000000000d1*o25*o37*o40*o41*o67*
     &
     &  o68*o74 + 3.7812500000000d0*o51*o53*o67*o68*o74
        t14 = -1.51300000000000d2*o25*o37*o40*o41*o55*o67*o68*o74 + 3.61
     &
     &  25000000000d-1*o51*o53*o55*o67*o68*o74 + 9.8000000000000d1*o58*o
     &
     &  59*o67*o68*o74 + 1.58420000000000d4*o55*o58*o59*o67*o68*o74 - 5.
     &
     &  0140000000000d2*o25*o37*o40*o41*o6*o67*o68*o74 + 2.3375000000000
     &
     &  0d0*o51*o53*o6*o67*o68*o74 + 2.49200000000000d3*o58*o59*o6*o67*o
     &
     &  68*o74 - 3.6000000000000d3*o14*o43*o52*o76 + 4.5000000000000d2*o
     &
     &  45*o52*o76 + 2.25000000000000d4*o43*o45*o52*o76 + 1.800000000000
     &
     &  00d3*o1*o13*o43*o67*o76 - 6.3000000000000d3*o1*o14*o43*o67*o76 +
     &
     &   5.5125000000000d3*o1*o43*o45*o67*o76 - 2.25000000000000d3*o52*o
     &
     &  73*o76 + 2.81250000000000d3*o13*o52*o73*o76 - 3.3750000000000d4*
     &
     &  o43*o52*o73*o76 + 7.2000000000000d3*o13*o52*o75*o76 - 5.40000000
     &
     &  00000d4*o14*o52*o75*o76 + 1.01250000000000d5*o45*o52*o75*o76 - 1
     &
     &  .61920000000000d4*o14*o17*o18*o25*o43*o52*o9 + 4.7828000000000d3
     &
     &  *o1*o14*o17*o18*o25*o43*o52*o9 + 8.4480000000000d3*o14*o21*o22*o
     &
     &  25*o43*o52*o9 - 2.37600000000000d3*o1*o14*o21*o22*o25*o43*o52*o9
     &
     &   - 1.12640000000000d3*o17*o25*o34*o43*o45*o52*o9 + 2.46400000000
     &
     &  000d2*o1*o17*o25*o34*o43*o45*o52*o9 - 3.5200000000000d2*o14*o17*
     &
     &  o18*o25*o43*o52*o54*o9 + 1.65000000000000d2*o14*o21*o22*o25*o43*
     &
     &  o52*o54*o9
        t15 = -1.10000000000000d1*o17*o25*o34*o43*o45*o52*o54*o9 - 1.540
     &
     &  48000000000d5*o14*o17*o18*o25*o43*o52*o55*o9 + 8.2058000000000d4
     &
     &  *o1*o14*o17*o18*o25*o43*o52*o55*o9 + 1.11360000000000d5*o14*o21*
     &
     &  o22*o25*o43*o52*o55*o9 - 8.1960000000000d4*o1*o14*o21*o22*o25*o4
     &
     &  3*o52*o55*o9 - 7.2100000000000d3*o14*o17*o18*o25*o43*o52*o54*o55
     &
     &  *o9 + 7.7250000000000d3*o14*o21*o22*o25*o43*o52*o54*o55*o9 + 1.0
     &
     &  1120000000000d5*o14*o17*o18*o25*o43*o52*o6*o9 - 4.4447200000000d
     &
     &  4*o1*o14*o17*o18*o25*o43*o52*o6*o9 - 6.4512000000000d4*o14*o21*o
     &
     &  22*o25*o43*o52*o6*o9 + 3.8256000000000d4*o1*o14*o21*o22*o25*o43*
     &
     &  o52*o6*o9 + 2.96960000000000d3*o17*o25*o34*o43*o45*o52*o6*o9 - 5
     &
     &  .1520000000000d2*o1*o17*o25*o34*o43*o45*o52*o6*o9 + 4.5092000000
     &
     &  000d3*o14*o17*o18*o25*o43*o52*o54*o6*o9 - 4.4340000000000d3*o14*
     &
     &  o21*o22*o25*o43*o52*o54*o6*o9 + 2.06000000000000d1*o17*o25*o34*o
     &
     &  43*o45*o52*o54*o6*o9 + 1.06240000000000d4*o1*o14*o17*o22*o25*o67
     &
     &  *o9 - 7.1126000000000d4*o1*o14*o17*o18*o4*o67*o9 + 1.62000000000
     &
     &  000d3*o1*o14*o21*o22*o4*o67*o9 + 2.87980000000000d4*o1*o14*o22*o
     &
     &  4*o40*o67*o9 - 5.1440000000000d3*o14*o17*o22*o25*o54*o67*o9 + 3.
     &
     &  01060000000000d4*o14*o17*o18*o4*o54*o67*o9 - 1.85960000000000d4*
     &
     &  o14*o22*o4*o40*o54*o67*o9 + 2.24640000000000d4*o1*o14*o17*o22*o2
     &
     &  5*o55*o67*o9 - 1.52406000000000d5*o1*o14*o17*o18*o4*o55*o67*o9 +
     &
     &   4.5000000000000d3*o1*o14*o21*o22*o4*o55*o67*o9 + 5.628600000000
     &
     &  0d4*o1*o14*o22*o4*o40*o55*o67*o9
        t16 = -1.00440000000000d4*o14*o17*o22*o25*o54*o55*o67*o9 + 6.622
     &
     &  2000000000d4*o14*o17*o18*o4*o54*o55*o67*o9 - 3.5676000000000d4*o
     &
     &  14*o22*o4*o40*o54*o55*o67*o9 + 3.09120000000000d4*o1*o14*o17*o22
     &
     &  *o25*o6*o67*o9 - 2.08212000000000d5*o1*o14*o17*o18*o4*o6*o67*o9 
     &
     &  + 5.4000000000000d3*o1*o14*o21*o22*o4*o6*o67*o9 + 8.067600000000
     &
     &  0d4*o1*o14*o22*o4*o40*o6*o67*o9 - 1.44120000000000d4*o14*o17*o22
     &
     &  *o25*o54*o6*o67*o9 + 8.9304000000000d4*o14*o17*o18*o4*o54*o6*o67
     &
     &  *o9 - 5.1648000000000d4*o14*o22*o4*o40*o54*o6*o67*o9 - 5.2096000
     &
     &  000000d4*o1*o14*o18*o62*o67*o9 + 1.69520000000000d4*o14*o18*o54*
     &
     &  o62*o67*o9 - 1.16928000000000d5*o1*o14*o18*o55*o62*o67*o9 + 3.78
     &
     &  36000000000d4*o14*o18*o54*o55*o62*o67*o9 - 1.56096000000000d5*o1
     &
     &  *o14*o18*o6*o62*o67*o9 + 5.0652000000000d4*o14*o18*o54*o6*o62*o6
     &
     &  7*o9 + 2.80000000000000d2*o14*o17*o22*o25*o67*o68*o9 - 1.7120000
     &
     &  0000000d3*o14*o17*o18*o4*o67*o68*o9 + 1.12000000000000d3*o14*o22
     &
     &  *o4*o40*o67*o68*o9 + 5.4000000000000d2*o14*o17*o22*o25*o55*o67*o
     &
     &  68*o9 - 3.8160000000000d3*o14*o17*o18*o4*o55*o67*o68*o9 + 2.1600
     &
     &  0000000000d3*o14*o22*o4*o40*o55*o67*o68*o9 + 7.8000000000000d2*o
     &
     &  14*o17*o22*o25*o6*o67*o68*o9 - 5.1120000000000d3*o14*o17*o18*o4*
     &
     &  o6*o67*o68*o9 + 3.12000000000000d3*o14*o22*o4*o40*o6*o67*o68*o9 
     &
     &  - 8.5600000000000d2*o14*o18*o62*o67*o68*o9 - 1.90800000000000d3*
     &
     &  o14*o18*o55*o62*o67*o68*o9 - 2.55600000000000d3*o14*o18*o6*o62*o
     &
     &  67*o68*o9
        ksq(i,j,k) = t1 + t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 + 
     &
     &  t11 + t12 + t13 + t14 + t15 + t16
            enddo
         enddo
      enddo
c
c
c
       ntries = 0
       
      do it = 1,itmax+1
         
         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
                  rhs(i,j,k) =
     $                 (psip(i+1,j,k)-2.*psip(i,j,k)+psip(i-1,j,k))/
     $                 deta**2+(psip(i,j+1,k)-2.*psip(i,j,k)+
     $                 psip(i,j-1,k))/dq**2+0.5*(psip(i,j+1,k)-
     $                 psip(i,j-1,k))/(dq*tan(qgrd(j)))+
     $                 (psip(i,j,k+1)-2.*psip(i,j,k)+psip(i,j,k-1))/
     $                 (dphi**2*sin(qgrd(j))**2)-0.25*psip(i,j,k)+0.125*
     &
     $                 Ksq(i,j,k)/(psisph(i,j,k)+psip(i,j,k))**7
               enddo
            enddo
         enddo
c
         rhs = -rhs
         
c
c
c
         rhsmax = 0.
c
         do k = 2,np-1
            do j = 2,nq-1
               do i = 2,ne-1
                  rhsmax = max(rhsmax,abs(rhs(i,j,k)))
               enddo
            enddo
         enddo
         
         if (verbose == 1) then
            ntries = ntries+1
            print*,'-----------Retry', ntries
            print*,'residual =', rhsmax
            print*,'--------------------------------------------------'
         endif
         
         if (rhsmax.lt.rbh_tol) goto 110
c
c
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./(dphi**2*
     $                 sin(qgrd(j))**2)-0.25-0.875*Ksq(i,j,k)/(psisph(i,
     &j,k)+
     $                 psip(i,j,k))**8
c
                  ae(i,j,k) = 1./deta**2
c
                  aw(i,j,k) = 1./deta**2
c
                  an(i,j,k) = 1./dq**2 + 0.5/(tan(qgrd(j))*dq)
c
                  as(i,j,k) = 1./dq**2 - 0.5/(tan(qgrd(j))*dq)
c
                  aq(i,j,k) = 1./(sin(qgrd(j))**2*dphi**2)
c
                  ab(i,j,k) = 1./(sin(qgrd(j))**2*dphi**2)
               enddo
            enddo
         enddo
c
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
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
c
c
c
            do i = 3,ne-2
               ac(i,2,k) = ac(i,2,k) + as(i,2,k)
               as(i,2,k) = 0.
c
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
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
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
         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
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
c
            ac(ne-1,2,k) = ac(ne-1,2,k) + 4.*ae(ne-1,2,k)/(3.+deta) +
     &           as(ne-1,2,k)
            aw(ne-1,2,k) = aw(ne-1,2,k) - ae(ne-1,2,k)/(3.+deta)
            ae(ne-1,2,k) = 0.
            as(ne-1,2,k) = 0.
c
c
c
            ac(ne-1,nq-1,k) = ac(ne-1,nq-1,k) + 4.*
     &           ae(ne-1,nq-1,k)/(3.+deta) + an(ne-1,nq-1,k)
            aw(ne-1,nq-1,k) = aw(ne-1,nq-1,k) -
     &           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
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
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
c
            ac(ne-1,j,2) = ac(ne-1,j,2) + 4.*ae(ne-1,j,2)/(3.+deta) +
     &           ab(ne-1,j,2)
            aw(ne-1,j,2) = aw(ne-1,j,2) - ae(ne-1,j,2)/(3.+deta)
            ae(ne-1,j,2) = 0.
            ab(ne-1,j,2) = 0.
c
c
c
            ac(ne-1,j,np-1) = ac(ne-1,j,np-1) + 4.*
     &           ae(ne-1,j,np-1)/(3.+deta) + aq(ne-1,j,np-1)
            aw(ne-1,j,np-1) = aw(ne-1,j,np-1) -
     &           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
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
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
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
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
         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
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
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
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
c
         ac(ne-1,2,2) = ac(ne-1,2,2) + 4.*ae(ne-1,2,2)/(3.+deta) +
     &        as(ne-1,2,2) + ab(ne-1,2,2)
         aw(ne-1,2,2) = aw(ne-1,2,2) - 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
         ac(ne-1,2,np-1) = ac(ne-1,2,np-1) + 4.*ae(ne-1,2,np-1)/(3.+deta
     &) +
     &        as(ne-1,2,np-1) + aq(ne-1,2,np-1)
         aw(ne-1,2,np-1) = aw(ne-1,2,np-1) - 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
c
         ac(ne-1,nq-1,2) = ac(ne-1,nq-1,2) + 4.*ae(ne-1,nq-1,2)/(3.+deta
     &) +
     &        an(ne-1,nq-1,2) + ab(ne-1,nq-1,2)
         aw(ne-1,nq-1,2) = aw(ne-1,nq-1,2) - 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
c
         ac(ne-1,nq-1,np-1) = ac(ne-1,nq-1,np-1) +
     &        4.*ae(ne-1,nq-1,np-1)/(3.+deta) + an(ne-1,nq-1,np-1) +
     &        aq(ne-1,nq-1,np-1)
         aw(ne-1,nq-1,np-1) = aw(ne-1,nq-1,np-1) - 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 rbicgst3d(ac,ae,aw,an,as,aq,ab,psipp,rhs,rbh_eps,rmax,ier,
     &
     $        ne,nq,np)
         if (rmax.gt.1.0e-10) then
            write(*,*) '***WARNING: bicgst3d did not converge.'
         endif
         if (ier.eq.-1) then
            write(*,*) '***WARNING: ier=-1'
         endif
c
c
c
         do k = 1,np
            do j = 1,nq
               psipp(1,j,k) = psipp(3,j,k)
               psipp(ne,j,k) = (4.*psipp(ne-1,j,k)-psipp(ne-2,j,k))/(3.+
     &
     $              deta)
            enddo
            do i = 1,ne
               psipp(i,1,k) = psipp(i,2,k)
               psipp(i,nq,k) = psipp(i,nq-1,k)
            enddo
         enddo
         do j = 1,nq
            do i = 1,ne
               psipp(i,j,1) = psipp(i,j,2)
               psipp(i,j,np) = psipp(i,j,np-1)
            enddo
         enddo
c
c
c
         do k = 1,np
            do j = 1,nq
               do i = 1,ne
                  psip(i,j,k) = psip(i,j,k) + psipp(i,j,k)
               enddo
            enddo
         enddo
c
      enddo
c
      write(*,*) 'It did not converge.'
      stop
c
 110  continue
      
      write(*,*) '--------------------------------------------------'
      print*,'Converge at Residual = ',rhsmax

c
c
c
      do k = 1,np
         do j = 1,nq
            do i = 2,ne-1
               detapsisph(i,j,k)=0.5*(psip(i+1,j,k)-psip(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*(psip(i,j+1,k)-psip(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*(psip(i,j,k+1)-psip(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)=(psip(i+1,j,k)-2.*psip(i,j,k)+
     &              psip(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)=psip(:,j,k)+2.0*cosh(0.5*etagrd)
         enddo
      enddo
c
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,680,"RotatingDBHIVP.F","RotatingDBHIVP","Cannot
     & create parameter table for interpolator")
      endif

      call CCTK_InterpHandle (interp_handle, "uniform cartesian")
      if (interp_handle .lt. 0) then
        call CCTK_Warn(0,685,"RotatingDBHIVP.F","RotatingDBHIVP","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

        gxx(i,j,k) = 1.00000000000000d0
c
c
c
        gxy(i,j,k) = 0
c
c
c
        gxz(i,j,k) = 0
c
c
c
        gyy(i,j,k) = 1.00000000000000d0
c
c
c
        gyz(i,j,k) = 0
c
c
c
        gzz(i,j,k) = 1.00000000000000d0
c
c
c
            enddo
         enddo
      enddo
c
      do k= 1,nz
         do j= 1,ny
            do i= 1,nx 

        gtil = amp*(exp(-((-eta0 + eta(i,j,k))**2/sigma**2)) + exp(-((et
     &a0 +
     &   eta(i,j,k))**2/sigma**2)))
        dngtil = amp*((-2.00000000000000d0*(-eta0 + eta(i,j,k))*exp(-((-
     &eta0
     &   + eta(i,j,k))**2/sigma**2)))/sigma**2 - (2.00000000000000d0*(et
     &a0 +
     &   eta(i,j,k))*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma**2)
        dnngtil = amp*((4.0000000000000d0*(-eta0 + eta(i,j,k))**2*exp(-(
     &(-et
     &  a0 + eta(i,j,k))**2/sigma**2)))/sigma**4 - (2.00000000000000d0*e
     &xp(-
     &  ((-eta0 + eta(i,j,k))**2/sigma**2)))/sigma**2 + (4.0000000000000
     &d0*(
     &  eta0 + eta(i,j,k))**2*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/s
     &igma**4 
     &  - (2.00000000000000d0*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/s
     &igma
     &  **2)
        dnnngtil = amp*((-8.0000000000000d0*(-eta0 + eta(i,j,k))**3*exp(
     &-((-
     &  eta0 + eta(i,j,k))**2/sigma**2)))/sigma**6 + (1.20000000000000d1
     &*(-e
     &  ta0 + eta(i,j,k))*exp(-((-eta0 + eta(i,j,k))**2/sigma**2)))/sigm
     &a**4 - (
     &  8.0000000000000d0*(eta0 + eta(i,j,k))**3*exp(-((eta0 + eta(i,j,k
     &))**2/si
     &  gma**2)))/sigma**6 + (1.20000000000000d1*(eta0 + eta(i,j,k))*exp
     &(-((
     &  eta0 + eta(i,j,k))**2/sigma**2)))/sigma**4)
        dnnnngtil = amp*((1.60000000000000d1*(-eta0 + eta(i,j,k))**4*exp
     &(-((
     &  -eta0 + eta(i,j,k))**2/sigma**2)))/sigma**8 - (4.8000000000000d1
     &*(-e
     &  ta0 + eta(i,j,k))**2*exp(-((-eta0 + eta(i,j,k))**2/sigma**2)))/s
     &igma**6 
     &  + (1.20000000000000d1*exp(-((-eta0 + eta(i,j,k))**2/sigma**2)))/
     &sigm
     &  a**4 + (1.60000000000000d1*(eta0 + eta(i,j,k))**4*exp(-((eta0 + 
     &eta(i,j,k)
     &  )**2/sigma**2)))/sigma**8 - (4.8000000000000d1*(eta0 + eta(i,j,k
     &))
     &  **2*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma**6 + (1.20000
     &0000
     &  00000d1*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma**4)
        dnnnnngtil = amp*((-3.2000000000000d1*(-eta0 + eta(i,j,k))**5*ex
     &p(-(
     &  (-eta0 + eta(i,j,k))**2/sigma**2)))/sigma**10 + (1.6000000000000
     &0d2*
     &  (-eta0 + eta(i,j,k))**3*exp(-((-eta0 + eta(i,j,k))**2/sigma**2))
     &)/sigma*
     &  *8 - (1.20000000000000d2*(-eta0 + eta(i,j,k))*exp(-((-eta0 + eta
     &(i,j,k))
     &  **2/sigma**2)))/sigma**6 - (3.2000000000000d1*(eta0 + eta(i,j,k)
     &)**5
     &  *exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma**10 + (1.6000000
     &0000
     &  000d2*(eta0 + eta(i,j,k))**3*exp(-((eta0 + eta(i,j,k))**2/sigma*
     &*2)))/si
     &  gma**8 - (1.20000000000000d2*(eta0 + eta(i,j,k))*exp(-((eta0 + e
     &ta(i,j,k)
     &  )**2/sigma**2)))/sigma**6)
        dnnnnnngtil = amp*((6.4000000000000d1*(-eta0 + eta(i,j,k))**6*ex
     &p(-(
     &  (-eta0 + eta(i,j,k))**2/sigma**2)))/sigma**12 - (4.8000000000000
     &d2*(
     &  -eta0 + eta(i,j,k))**4*exp(-((-eta0 + eta(i,j,k))**2/sigma**2)))
     &/sigma**
     &  10 + (7.2000000000000d2*(-eta0 + eta(i,j,k))**2*exp(-((-eta0 + e
     &ta(i,j,k)
     &  )**2/sigma**2)))/sigma**8 - (1.20000000000000d2*exp(-((-eta0 + 
     &  eta(i,j,k))**2/sigma**2)))/sigma**6 + (6.4000000000000d1*(eta0 +
     & eta(i,j,k)
     &  )**6*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma**12 - (4.800
     &0
     &  000000000d2*(eta0 + eta(i,j,k))**4*exp(-((eta0 + eta(i,j,k))**2/
     &sigma**2
     &  )))/sigma**10 + (7.2000000000000d2*(eta0 + eta(i,j,k))**2*exp(-(
     &(eta
     &  0 + eta(i,j,k))**2/sigma**2)))/sigma**8 - (1.20000000000000d2*ex
     &p(-(
     &  (eta0 + eta(i,j,k))**2/sigma**2)))/sigma**6)
        dnnnnnnngtil = amp*((-1.28000000000000d2*(-eta0 + eta(i,j,k))**7
     &*exp
     &  (-((-eta0 + eta(i,j,k))**2/sigma**2)))/sigma**14 + (1.3440000000
     &0000
     &  d3*(-eta0 + eta(i,j,k))**5*exp(-((-eta0 + eta(i,j,k))**2/sigma**
     &2)))/sig
     &  ma**12 - (3.3600000000000d3*(-eta0 + eta(i,j,k))**3*exp(-((-eta0
     & + e
     &  ta(i,j,k))**2/sigma**2)))/sigma**10 + (1.68000000000000d3*(-eta0
     & + e
     &  ta(i,j,k))*exp(-((-eta0 + eta(i,j,k))**2/sigma**2)))/sigma**8 - 
     &(1.28000
     &  000000000d2*(eta0 + eta(i,j,k))**7*exp(-((eta0 + eta(i,j,k))**2/
     &sigma**2
     &  )))/sigma**14 + (1.34400000000000d3*(eta0 + eta(i,j,k))**5*exp(-
     &((et
     &  a0 + eta(i,j,k))**2/sigma**2)))/sigma**12 - (3.3600000000000d3*(
     &eta0
     &   + eta(i,j,k))**3*exp(-((eta0 + eta(i,j,k))**2/sigma**2)))/sigma
     &**10 + (
     &  1.68000000000000d3*(eta0 + eta(i,j,k))*exp(-((eta0 + eta(i,j,k))
     &**2/sigm
     &  a**2)))/sigma**8)

        o1 = mm**2
        o2 = 5.0000000000000d1*o1
        o3 = -8.0000000000000d1 + o2
        o4 = 1/o3
        o5 = 2.00000000000000d0*q(i,j,k)
        o6 = cos(o5)
        o7 = 5.0000000000000d0*o6
        o8 = 3.00000000000000d0 + o7
        o9 = dnnnnngtil
        o10 = dnnngtil
        o11 = -2.60000000000000d1*o10
        o12 = dngtil
        o13 = 2.50000000000000d1*o12
        o14 = o11 + o13 + o9
        o15 = -3.00000000000000d0*eta(i,j,k)
        o16 = exp(o15)
        o17 = sin(phi(i,j,k))
        o18 = o17**2
        o19 = mm*phi(i,j,k)
        o20 = sin(o19)
        o21 = sin(q(i,j,k))
        o22 = o21**2
        o23 = 2.50000000000000d1*o1
        o24 = -4.0000000000000d1 + o23
        o25 = 1/o24
        o26 = cos(phi(i,j,k))
        o27 = o26**2
        o28 = -1.60000000000000d1 + o1
        o29 = 3.00000000000000d0*o6
        o30 = 2.00000000000000d0 + o29
        o31 = 2.00000000000000d0*o28*o30*o9
        o32 = -1.07000000000000d2*o1
        o33 = 6.0900000000000d2*o6
        o34 = -1.59000000000000d2*o1*o6
        o35 = 4.0700000000000d2 + o32 + o33 + o34
        o36 = o10*o35
        o37 = 3.5000000000000d1*o1
        o38 = 5.0000000000000d0*o1
        o39 = -1.30000000000000d1 + o38
        o40 = 9.0000000000000d0*o39*o6
        o41 = -8.3000000000000d1 + o37 + o40
        o42 = 5.0000000000000d0*o12*o41
        o43 = o31 + o36 + o42
        o44 = o22**2
        o45 = cos(q(i,j,k))
        o46 = o45**2
        o47 = -3.00000000000000d0*mm*o14*o20*o22*o4*o8
        o48 = -(mm*o20*o22*o25*o43)
        o49 = o47 + o48
        o50 = -8.0000000000000d0 + o38
        o51 = 1/o50
        o52 = cos(o19)
        o53 = -9.0000000000000d0*o1
        o54 = 5.0000000000000d0*o50*o6
        o55 = 8.0000000000000d0 + o53 + o54
        o56 = -7.5000000000000d1*o12*o55
        o57 = 2.88000000000000d2*o9
        o58 = -1.20000000000000d1*o1*o9
        o59 = -9.4000000000000d1*o1
        o60 = -4.1500000000000d2*o6
        o61 = 1.75000000000000d2*o1*o6
        o62 = 1.15000000000000d2 + o59 + o60 + o61
        o63 = 4.0000000000000d0*o10*o62
        o64 = dnnnnnnngtil
        o65 = 3.2000000000000d1*o22*o64
        o66 = -2.00000000000000d0*o1*o22*o64
        o67 = -9.2800000000000d2*o22*o9
        o68 = 1.03000000000000d2*o1*o22*o9
        o69 = o57 + o58 + o63 + o65 + o66 + o67 + o68
        o70 = 2.00000000000000d0*o69
        o71 = o56 + o70
        o72 = sin(o5)
        o73 = dnnnnnngtil
        o74 = 3.2000000000000d1*o73
        o75 = -2.00000000000000d0*o1*o73
        o76 = -3.2000000000000d1*o6*o73
        o77 = 2.00000000000000d0*o1*o6*o73
        o78 = dnnnngtil
        o79 = -3.4000000000000d2*o78
        o80 = 5.5000000000000d1*o1*o78
        o81 = -8.1200000000000d2*o6*o78
        o82 = 1.70000000000000d1*o1*o6*o78
        o83 = 7.0000000000000d0*o1
        o84 = -2.09000000000000d2*o6
        o85 = 8.9000000000000d1*o1*o6
        o86 = -3.10000000000000d1 + o83 + o84 + o85
        o87 = dnngtil
        o88 = -2.00000000000000d1*o86*o87
        o89 = 7.0000000000000d0*o6
        o90 = 1.00000000000000d0 + o89
        o91 = gtil
        o92 = 7.5000000000000d1*o50*o90*o91
        o93 = o74 + o75 + o76 + o77 + o79 + o80 + o81 + o82 + o88 + o92
        o94 = -2.40000000000000d2*byJ
        o95 = 1.50000000000000d2*byJ*o1
        o96 = -1.15200000000000d3*o46*o52*o78
        o97 = 7.2000000000000d1*o1*o46*o52*o78
        o98 = 3.9000000000000d1*o1
        o99 = -2.68000000000000d2*o6
        o100 = 1.00000000000000d2*o1*o6
        o101 = 4.0000000000000d0*q(i,j,k)
        o102 = cos(o101)
        o103 = -5.8100000000000d2*o102
        o104 = 2.45000000000000d2*o1*o102
        o105 = -1.11000000000000d2 + o100 + o103 + o104 + o98 + o99
        o106 = -5.0000000000000d0*o105*o52*o87
        o107 = -1.92000000000000d2*o22*o46*o52*o73
        o108 = 1.20000000000000d1*o1*o22*o46*o52*o73
        o109 = 2.88000000000000d2*o22*o52*o78
        o110 = -1.80000000000000d1*o1*o22*o52*o78
        o111 = 5.5680000000000d3*o22*o46*o52*o78
        o112 = -6.1800000000000d2*o1*o22*o46*o52*o78
        o113 = 3.2000000000000d1*o44*o52*o73
        o114 = -2.00000000000000d0*o1*o44*o52*o73
        o115 = -9.2800000000000d2*o44*o52*o78
        o116 = 1.03000000000000d2*o1*o44*o52*o78
        o117 = o106 + o107 + o108 + o109 + o110 + o111 + o112 + o113 + o
     &
     &  114 + o115 + o116 + o94 + o95 + o96 + o97
        o118 = 2.00000000000000d0*o117
        o119 = 2.00000000000000d1*o6
        o120 = 3.5000000000000d1*o102
        o121 = 9.0000000000000d0 + o119 + o120
        o122 = 1.87500000000000d1*o121*o50*o52*o91
        o123 = o118 + o122
        o124 = o21*o22
        kxx(i,j,k) = mm*o16*o20*o25*o27*o43*o44 + o16*o27*o46*o49 - 1.00
     &
     &  000000000000d-1*o123*o16*o17*o22*o26*o51 - 1.00000000000000d-1*o
     &
     &  16*o17*o22*o26*o46*o51*o52*o71 + 3.00000000000000d0*mm*o14*o16*o
     &
     &  18*o20*o22*o4*o8 + 5.0000000000000d-2*mm*o16*o20*o21*o27*o45*o51
     &
     &  *o72*o93
        kxy(i,j,k) = mm*o16*o17*o20*o25*o26*o43*o44 + o16*o17*o26*o46*o4
     &
     &  9 - 5.0000000000000d-2*o123*o16*o18*o22*o51 + 5.0000000000000d-2
     &
     &  *o123*o16*o22*o27*o51 - 5.0000000000000d-2*o16*o18*o22*o46*o51*o
     &
     &  52*o71 + 5.0000000000000d-2*o16*o22*o27*o46*o51*o52*o71 - 3.0000
     &
     &  0000000000d0*mm*o14*o16*o17*o20*o22*o26*o4*o8 + 5.0000000000000d
     &
     &  -2*mm*o16*o17*o20*o21*o26*o45*o51*o72*o93
        kxz(i,j,k) = mm*o124*o16*o20*o25*o26*o43*o45 - o16*o21*o26*o45*o
     &
     &  49 - 5.0000000000000d-2*o123*o16*o17*o21*o45*o51 + 5.00000000000
     &
     &  00d-2*o124*o16*o17*o45*o51*o52*o71 - 2.50000000000000d-2*mm*o16*
     &
     &  o20*o22*o26*o51*o72*o93 + 2.50000000000000d-2*mm*o16*o20*o26*o46
     &
     &  *o51*o72*o93
        kyy(i,j,k) = mm*o16*o18*o20*o25*o43*o44 + o16*o18*o46*o49 + 1.00
     &
     &  000000000000d-1*o123*o16*o17*o22*o26*o51 + 1.00000000000000d-1*o
     &
     &  16*o17*o22*o26*o46*o51*o52*o71 + 3.00000000000000d0*mm*o14*o16*o
     &
     &  20*o22*o27*o4*o8 + 5.0000000000000d-2*mm*o16*o18*o20*o21*o45*o51
     &
     &  *o72*o93
        kyz(i,j,k) = mm*o124*o16*o17*o20*o25*o43*o45 - o16*o17*o21*o45*o
     &
     &  49 + 5.0000000000000d-2*o123*o16*o21*o26*o45*o51 - 5.00000000000
     &
     &  00d-2*o124*o16*o26*o45*o51*o52*o71 - 2.50000000000000d-2*mm*o16*
     &
     &  o17*o20*o22*o51*o72*o93 + 2.50000000000000d-2*mm*o16*o17*o20*o46
     &
     &  *o51*o72*o93
        kzz(i,j,k) = mm*o16*o20*o22*o25*o43*o46 + o16*o22*o49 - 5.000000
     &
     &  0000000d-2*mm*o16*o20*o21*o45*o51*o72*o93
            enddo
         enddo
      enddo
      
      kxx = kxx/psi**2
      kxy = kxy/psi**2
      kxz = kxz/psi**2
      kyy = kyy/psi**2
      kyz = kyz/psi**2
      kzz = kzz/psi**2

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
c
      print*,'Angular momentum parameter: a/m = ', byJ/adm**2
c
c
c
c
c

      deallocate(ac,ae,aw,an,as,aq,ab,rhs,Ksq,psisph,psip,psipp,
     $     detapsisph,dqpsisph,dphipsisph,detaetapsisph,detaqpsisph,
     $     detaphipsisph,dqqpsisph,dqphipsisph,dphiphipsisph)
      deallocate(etagrd,qgrd,phigrd)
      
      return
      end

      
      
