














































c














































































 


































































































































































 






























































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                



















































































c

      subroutine IDAxiOddBrillBH(cctk_dim,cctk_gsh,cctk_lsh,cctk_lbnd,cc
     &tk_ubnd,cctk_tile_min,cctk_tile_max,cctk_ash,cctk_alignment,cctk_a
     &lignment_offset,cctk_from,cctk_to,cctk_bbox,cctk_delta_time,cctk_t
     &ime,cctk_delta_space,cctk_origin_space,cctk_levfac,cctk_levoff,cct
     &k_levoffdenom,cctk_timefac,cctk_convlevel,cctk_convfac,cctk_nghost
     &zones,cctk_iteration,cctkGH,cctk_ash1,cctk_ash2,cctk_ash3, abseta,
     &detaetapsi2d,detapsi2d,detaqpsi2d,dqpsi2d,dqqpsi2d,eta,phi,psi2d,q
     &,sign_eta, alp,alp_p,alp_p_p,betax,betax_p,betax_p_p,betay,betay_p
     &,betay_p_p,betaz,betaz_p,betaz_p_p,coarse_dx,coarse_dy,coarse_dz,c
     &onformal_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,dtl
     &apse_state,dtshift_state,gxx,gxx_p,gxx_p_p,gxy,gxy_p,gxy_p_p,gxz,g
     &xz_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,kxy,kxy_p,kxy_p_p,kxz,kxz_p,kxz_p_p,kyy,kyy_p,k
     &yy_p_p,kyz,kyz_p,kyz_p_p,kzz,kzz_p,kzz_p_p,psi,psix,psixx,psixy,ps
     &ixz,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 detaetapsi2d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_detaetapsi2d = kind(detaetapsi2d)
      REAL*8 detapsi2d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_detapsi2d = kind(detapsi2d)
      REAL*8 detaqpsi2d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_detaqpsi2d = kind(detaqpsi2d)
      REAL*8 dqpsi2d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dqpsi2d = kind(dqpsi2d)
      REAL*8 dqqpsi2d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_dqqpsi2d = kind(dqqpsi2d)
      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 psi2d (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psi2d = kind(psi2d)
      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 amp_me
      integer, parameter :: cctki_use_amp_me = kind(amp_me)
      REAL*8 byJ
      integer, parameter :: cctki_use_byJ = kind(byJ)
      REAL*8 eta0
      integer, parameter :: cctki_use_eta0 = kind(eta0)
      REAL*8 eta0_me
      integer, parameter :: cctki_use_eta0_me = kind(eta0_me)
      REAL*8 etamax
      integer, parameter :: cctki_use_etamax = kind(etamax)
      REAL*8 sigma
      integer, parameter :: cctki_use_sigma = kind(sigma)
      REAL*8 sigma_me
      integer, parameter :: cctki_use_sigma_me = kind(sigma_me)
      REAL*8 xi
      integer, parameter :: cctki_use_xi = kind(xi)
      INTEGER*4 brandt_seidel
      integer, parameter :: cctki_use_brandt_seidel = kind(brandt_seidel
     &)
      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 n_me
      integer, parameter :: cctki_use_n_me = kind(n_me)
      INTEGER*4 neta
      integer, parameter :: cctki_use_neta = kind(neta)
      INTEGER*4 ntheta
      integer, parameter :: cctki_use_ntheta = kind(ntheta)
      INTEGER*4 sergio
      integer, parameter :: cctki_use_sergio = kind(sergio)
      INTEGER*4 verbose
      integer, parameter :: cctki_use_verbose = kind(verbose)
      COMMON /IDAxiOddBrillBHpriv/amp, amp_me, byJ, eta0, eta0_me, etama
     &x, sigma, sigma_me, xi, brandt_seidel, interpolation_order, n, n_m
     &e, neta, ntheta, sergio, verbose
      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
      REAL*8, allocatable :: ac(:,:),ae(:,:),aw(:,:),an(:,:),as(:,:),
     $     rhs(:,:),qfetaeta(:,:),qfqq(:,:),ksq(:,:), 
     $     psi2dv(:,:),dpsi2dv(:,:),ddpsi2dv(:,:),
     $     detapsisph(:,:),dqpsisph(:,:),detaetapsisph(:,:),
     $     detaqpsisph(:,:),dqqpsisph(:,:),
     $     etagrd(:),qgrd(:)
      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
      REAL*8 t1,t2,t3,t4
      REAL*8 gtil,dngtil,dnngtil,dnnngtil,dnnnngtil,dnnnnngtil,
     $     dnnnnnngtil,dnnnnnnngtil
      REAL*8 rhsmax,rmax,odd_get2d,adm,Jmom,mass,a,r_iso,rBL,rbh_eps
      REAL*8,parameter :: rbh_tol = 1.0d-7
      integer,parameter :: itmax = 100 
      REAL*8 pi,zero,one,two
      REAL*8 :: xmin,xmax,ymin,ymax,zmin,zmax
      integer :: ne, nq
      integer :: nx,ny,nz
      integer :: i,j,k,it,ier,nquads,ntries
      integer :: npoints,handle,ierror
      integer make_conformal_derivs

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


      rbh_eps = 1.0d-10
      pi = 4.0d0*atan(1.0d0)
c
      ne = neta
      nq = ntheta

c
      nx = cctk_lsh(1)
      ny = cctk_lsh(2)
      nz = cctk_lsh(3)      

c
c
      if (brandt_seidel == 1) then
         write(*,*) '"Solve Bradt-Seidel Initial Data Sets"'
      endif

      write(*,*) 'Brill Wave: amplitude for K :',amp
      write(*,*) 'Brill Wave: center for K:',eta0
      write(*,*) 'Brill Wave: sigma for K:',sigma
      write(*,*) 'Bowen & York Momenta:  :',byJ
      write(*,*) 'Brill Wave: sin^n theta for K:',n

      write(*,*) 'Brill Wave: amplitude for g :',amp_me
      write(*,*) 'Brill Wave: center for g:',eta0_me
      write(*,*) 'Brill Wave: sigma for g:',sigma_me
      write(*,*) 'Brill Wave: sin^n theta for g:',n_me


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

c

c
c
c
      ne = ne+2
      nq = nq+2
c
      allocate(ac(ne,nq),ae(ne,nq),aw(ne,nq),an(ne,nq),as(ne,nq),
     $     rhs(ne,nq),qfetaeta(ne,nq),qfqq(ne,nq),ksq(ne,nq),
     $     psi2dv(ne,nq),dpsi2dv(ne,nq),
     $     ddpsi2dv(ne,nq),detapsisph(ne,nq),dqpsisph(ne,nq),
     $     detaetapsisph(ne,nq),detaqpsisph(ne,nq),dqqpsisph(ne,nq),
     $     etagrd(ne),qgrd(nq))
c
c
c
      nquads = 2      
      dq = nquads*0.5*pi/(nq-2)
      deta = etamax/(ne-3)
      
      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 j = 1,nq
          do i = 1,ne

        o1 = -xi
        o2 = cos(qgrd(j))
        o3 = xi*o2
        o4 = 1.0 + o1 + o3
        o5 = -eta0_me
        o6 = etagrd(i) + o5
        o7 = o6**2
        o8 = sigma_me**2
        o9 = o8**2
        o10 = 1/o9
        o11 = 1/o8
        o12 = -(o11*o7)
        o13 = exp(o12)
        o14 = eta0_me + etagrd(i)
        o15 = o14**2
        o16 = -(o11*o15)
        o17 = exp(o16)
        o18 = sin(qgrd(j))
        o19 = o18**n_me
        o20 = o2**2
        o21 = o13 + o17
        qfetaeta(i,j) = amp_me*o19*o4*(-2.0*o11*o13 - 2.0*o11*o17 + 4.0*
     &o10*o15*o
     &  17 + 4.0*o10*o13*o7)
        qfqq(i,j) = -(amp_me*xi*n_me*o19*o2*o21) - amp_me*xi*(1.0 + n_me
     &)*o19*o2*o21 -
     &   amp_me*n_me*o19*o21*o4 + amp_me*(-1.0 + n_me)*n_me*o18**(-2 + n
     &)*o20*o21*o
     &  4
          enddo
       enddo


c
c
c
      dpsi2dv = 0.
      ddpsi2dv = 0.
      do j = 1,nq
         do i = 1,ne
            psi2dv(i,j) = 2.*cosh(0.5*etagrd(i))
         enddo
      enddo
c
c
c
      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)
            if(brandt_seidel == 1) then

        o1 = cos(qgrd(j))
        o2 = 1.0
        o3 = o2**2
        o4 = 1/o3
        o5 = sin(qgrd(j))
        o6 = 2.0 + n
        o7 = o5**2
        o8 = -(o6*o7)
        o9 = 1.0 + n + o8
        o10 = o3**2
        o11 = o10**2
        o12 = dngtil**2
        o13 = o1**2
        o14 = gtil**2
        o15 = o9**2
c
c
c
c
c
c
        ksq(i,j) = (2.0*o5**(2*(
     $	-2 + n))*(o14*o15 + o12*o13*o7))/(o10*o11)
	    else if (sergio == 1) then

        mass = sqrt(2.0 + sqrt(4.0 + byJ**2
     &  ))
        a = byJ/mass
        rBL = 5.0e-1*exp(etagrd(i))*(1.0 + (1.0*(a
     &   + mass)*exp(-etagrd(i)))/sqrt(-a**2 + mass**2))*sqrt(-a**2 + ma
     &ss**2)
     &  *(1.0 + (1.0*(a + mass)*exp(-etagrd(i)))/sqrt(a**2 + m
     &  ass**2))

        o1 = a**2
        o2 = -2.00000000000000d0*mass*rBL
        o3 = rBL**2
        o4 = o1 + o2 + o3
        o5 = sqrt(o4)
        o6 = cos(qgrd(j))
        o7 = o6**2
        o8 = o1*o7
        o9 = o3 + o8
        o10 = o9**2
        o11 = 1/o10
        o12 = 1.0!psi2dv(eta,q)
        o13 = o12**2
        o14 = 1/o13
        o15 = sin(qgrd(j))
        o16 = o15**2
        o17 = o1 + o3
        o18 = 2.00000000000000d0*o17*o3
        o19 = -o1
        o20 = o19 + o3
        o21 = o20*o9
        o22 = o18 + o21
        o23 = mass**2
        o24 = o10**2
        o25 = 1/o24
        o26 = o22**2
        o27 = o13**2
        o28 = o27**2
        o29 = o27*o28
        o30 = 1/o29
        o31 = o1**2
        o32 = o16**2
c
c
c
c
c
c
c
        ksq(i,j) = 2.00000000000000d0*o1*o16*o23*o25*o26*o30 + 8.0000000
     &
     &  000000d0*o1*o23*o25*o3*o30*o31*o32*o4*o7
            else

        o1 = cos(qgrd(j))
        o2 = 2.00000000000000d0*qgrd(j)
        o3 = cos(o2)
        o4 = dnnnnngtil
        o5 = sin(qgrd(j))
        o6 = o5**2
        o7 = o5*o6
        o8 = dnnngtil
        o9 = dngtil
        o10 = dnnnnnnngtil
        o11 = o6**2
        o12 = 4.0000000000000d0*qgrd(j)
        o13 = cos(o12)
        o14 = dnnnngtil
        o15 = dnngtil
        o16 = o1**2
        o17 = dnnnnnngtil
        o18 = o11*o6
        o19 = gtil
        o20 = byJ**2
        o21 = o14**2
        o22 = o3**2
        o23 = o13**2
        o24 = o15**2
        o25 = o4**2
        o26 = o8**2
        o27 = o9**2
        o28 = o16**2
        o29 = o17**2
        o30 = o10**2
        o31 = o11**2
        o32 = o31*o6
        o33 = o19**2
        t1 = 2.88000000000000d1*byJ*o11*o16*o17 + 1.00800000000000d1*o11
     &*o
     &  14*o16*o17 + 9.7440000000000d1*o11*o13*o14*o16*o17 - 6.660000000
     &
     &  0000d1*o11*o15*o16*o17 - 3.4860000000000d2*o11*o13*o15*o16*o17 -
     &
     &   4.8000000000000d0*byJ*o17*o18 - 1.68000000000000d0*o14*o17*o18 
     &- 
     &  1.62400000000000d1*o13*o14*o17*o18 + 1.11000000000000d1*o15*o17*
     &
     &  o18 + 5.8100000000000d1*o13*o15*o17*o18 - 1.80000000000000d2*byJ
     &*o
     &  11*o19 - 6.3000000000000d1*o11*o14*o19 - 6.0900000000000d2*o11*o
     &
     &  13*o14*o19 + 4.1625000000000d2*o11*o15*o19 + 2.17875000000000d3*
     &
     &  o11*o13*o15*o19 - 2.70000000000000d3*byJ*o11*o16*o19 - 9.4500000
     &00
     &  0000d2*o11*o14*o16*o19 - 9.1350000000000d3*o11*o13*o14*o16*o19 +
     &
     &   6.2437500000000d3*o11*o15*o16*o19 + 3.2681250000000d4*o11*o13*o
     &
     &  15*o16*o19 + 4.5000000000000d2*byJ*o18*o19 + 1.57500000000000d2*
     &o1
     &  4*o18*o19 + 1.52250000000000d3*o13*o14*o18*o19 - 1.0406250000000
     &
     &  0d3*o15*o18*o19 - 5.4468750000000d3*o13*o15*o18*o19 - 2.40000000
     &
     &  000000d2*o16*o17*o18*o19 + 9.6800000000000d0*o11*o16*o25
        t2 = 6.7280000000000d1*o11*o16*o22*o25 + 6.6125000000000d1*o11*o
     &
     &  16*o26 + 8.6112500000000d2*o11*o16*o22*o26 + 2.81250000000000d1*
     &
     &  o11*o16*o27 + 7.0312500000000d2*o11*o16*o22*o27 + 5.760000000000
     &
     &  0d2*o11*o17*o19*o28 - 2.16000000000000d3*o17*o18*o19*o28 + 1.152
     &
     &  00000000000d1*o18*o28*o29 + 3.07200000000000d1*o11*o14*o16*o17*o
     &
     &  3 - 1.60800000000000d2*o11*o15*o16*o17*o3 - 5.1200000000000d0*o1
     &
     &  4*o17*o18*o3 + 2.68000000000000d1*o15*o17*o18*o3 - 1.92000000000
     &
     &  000d2*o11*o14*o19*o3 + 1.00500000000000d3*o11*o15*o19*o3 - 2.880
     &
     &  00000000000d3*o11*o14*o16*o19*o3 + 1.50750000000000d4*o11*o15*o1
     &
     &  6*o19*o3 + 4.8000000000000d2*o14*o18*o19*o3 - 2.51250000000000d3
     &
     &  *o15*o18*o19*o3 - 5.1040000000000d1*o11*o16*o25*o3 - 4.772500000
     &
     &  0000d2*o11*o16*o26*o3 - 2.81250000000000d2*o11*o16*o27*o3 + 2.40
     &
     &  000000000000d1*o17*o19*o31 + 7.2000000000000d2*o16*o17*o19*o31 -
     &
     &   3.8400000000000d0*o16*o29*o31 + 3.2000000000000d-1*o16*o30*o31 
     &
     &  - 6.0000000000000d1*o17*o19*o32 + 3.2000000000000d-1*o29*o32 - 3
     &
     &  .6000000000000d3*o11*o16*o33
        t3 = 4.5000000000000d2*o18*o33 + 2.25000000000000d4*o16*o18*o33 
     &
     &  - 5.4000000000000d4*o11*o28*o33 + 1.01250000000000d5*o18*o28*o33
     &
     &   - 2.25000000000000d3*o31*o33 - 3.3750000000000d4*o16*o31*o33 + 
     &
     &  2.81250000000000d3*o32*o33 - 3.5200000000000d0*o10*o16*o18*o4 + 
     &
     &  9.2800000000000d0*o10*o16*o18*o3*o4 + 1.26000000000000d1*byJ*o14
     &*o
     &  6 + 1.21800000000000d2*byJ*o13*o14*o6 - 8.3250000000000d1*byJ*o1
     &5*o6
     &   - 4.3575000000000d2*byJ*o13*o15*o6 - 2.91375000000000d1*o14*o15
     &*o
     &  6 - 4.3417500000000d2*o13*o14*o15*o6 + 7.2000000000000d2*byJ*o16
     &*o
     &  19*o6 + 2.52000000000000d2*o14*o16*o19*o6 + 2.43600000000000d3*o
     &
     &  13*o14*o16*o19*o6 - 1.66500000000000d3*o15*o16*o19*o6 - 8.715000
     &
     &  0000000d3*o13*o15*o16*o19*o6 + 1.80000000000000d1*o20*o6 + 2.205
     &
     &  00000000000d0*o21*o6 + 4.2630000000000d1*o13*o21*o6 - 2.14400000
     &
     &  000000d2*o14*o15*o22*o6 + 2.04800000000000d1*o21*o22*o6 - 1.4742
     &
     &  8750000000d3*o14*o15*o23*o6 + 2.06045000000000d2*o21*o23*o6
        t4 = 9.6257812500000d1*o24*o6 + 1.00767187500000d3*o13*o24*o6 + 
     &
     &  5.6112500000000d2*o22*o24*o6 + 2.63719531250000d3*o23*o24*o6 + 3
     &
     &  .8400000000000d1*byJ*o14*o3*o6 - 2.01000000000000d2*byJ*o15*o3*o
     &6 - 
     &  1.59150000000000d2*o14*o15*o3*o6 - 1.14485000000000d3*o13*o14*o1
     &
     &  5*o3*o6 + 7.6800000000000d2*o14*o16*o19*o3*o6 - 4.0200000000000d
     &
     &  3*o15*o16*o19*o3*o6 + 1.34400000000000d1*o21*o3*o6 + 1.299200000
     &
     &  00000d2*o13*o21*o3*o6 + 4.6481250000000d2*o24*o3*o6 + 2.43293750
     &
     &  000000d3*o13*o24*o3*o6 + 7.2000000000000d3*o28*o33*o6 + 9.200000
     &
     &  0000000d0*o10*o16*o18*o8 - 3.3200000000000d1*o10*o16*o18*o3*o8 -
     &
     &   5.0600000000000d1*o11*o16*o4*o8 - 4.8140000000000d2*o11*o16*o22
     &
     &  *o4*o8 + 3.16000000000000d2*o11*o16*o3*o4*o8 - 6.0000000000000d0
     &
     &  *o10*o16*o18*o9 + 3.00000000000000d1*o10*o16*o18*o3*o9 + 3.30000
     &
     &  00000000d1*o11*o16*o4*o9 + 4.3500000000000d2*o11*o16*o22*o4*o9 -
     &
     &   2.52000000000000d2*o11*o16*o3*o4*o9 - 8.6250000000000d1*o11*o16
     &
     &  *o8*o9 - 1.55625000000000d3*o11*o16*o22*o8*o9 + 7.4250000000000d
     &
     &  2*o11*o16*o3*o8*o9
        ksq(i,j) = t1 + t2 + t3 + t4
            endif
         enddo
      enddo
c
c
c

      ntries = 0

      do it = 1,itmax+1
         
         ac = 0.
         ae = 0.
         aw = 0.
         an = 0.
         as = 0.
         rhs = 0.

c
         do j = 2,nq-1
            do i = 2,ne-1
               rhs(i,j) =
     $              (dpsi2dv(i+1,j)-2.*dpsi2dv(i,j)+dpsi2dv(i-1,j))/
     $              deta**2+(dpsi2dv(i,j+1)-2.*dpsi2dv(i,j)+
     $              dpsi2dv(i,j-1))/dq**2+0.5*(dpsi2dv(i,j+1)-
     $              dpsi2dv(i,j-1))/(dq*tan(qgrd(j)))
     $              +0.25*dpsi2dv(i,j)*
     $              (qfetaeta(i,j)+qfqq(i,j)-1.0d0)+
     $              0.125*ksq(i,j)/(psi2dv(i,j)+dpsi2dv(i,j))**7
            enddo
         enddo
c
         rhs = -rhs
c
c
c
         rhsmax = 0.
c
         do j = 2,nq-1
            do i = 2,ne-1
               rhsmax = max(rhsmax,abs(rhs(i,j)))
            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 j = 2,nq-1
            do i = 2,ne-1
               ac(i,j) = -2./deta**2-2./dq**2-sin(qgrd(j))**2+0.25
     $              *(qfetaeta(i,j)+qfqq(i,j) - 1.0d0)
     $              -0.875*ksq(i,j)/(psi2dv(i,j)+ 
     $              dpsi2dv(i,j))**8
c
               ae(i,j) = 1./deta**2
c
               aw(i,j) = 1./deta**2
c
               an(i,j) = 1./dq**2 + 0.5/(tan(qgrd(j))*dq)
c
               as(i,j) = 1./dq**2 - 0.5/(tan(qgrd(j))*dq)
c
            enddo
         enddo
         
c
c
c
c
c
         do j = 3,nq-2
            ae(2,j) = ae(2,j) + aw(2,j)
            aw(2,j) = 0.
c
c
c
            ac(ne-1,j) = ac(ne-1,j) + 4.*ae(ne-1,j)/(3.+deta)
            aw(ne-1,j) = aw(ne-1,j) - ae(ne-1,j)/(3.+deta)
            ae(ne-1,j) = 0.
         enddo
c
c
c
         do i = 3,ne-2
            ac(i,2) = ac(i,2) + as(i,2)
            as(i,2) = 0.
c
c
c
            ac(i,nq-1) = ac(i,nq-1) + an(i,nq-1)
            an(i,nq-1) = 0.
         enddo
         
c
c
c
         ae(2,2) = ae(2,2) + aw(2,2)
         ac(2,2) = ac(2,2) + as(2,2)
         aw(2,2) = 0.
         as(2,2) = 0.
c
c
c
         ae(2,nq-1) = ae(2,nq-1) + aw(2,nq-1)
         ac(2,nq-1) = ac(2,nq-1) + an(2,nq-1)
         aw(2,nq-1) = 0.
         an(2,nq-1) = 0.
c
c
c
         ac(ne-1,2) = ac(ne-1,2) + 4.*ae(ne-1,2)/(3.+deta) + 
     $        as(ne-1,2)
         aw(ne-1,2) = aw(ne-1,2) - ae(ne-1,2)/(3.+deta)
         ae(ne-1,2) = 0.
         as(ne-1,2) = 0.
c
c
c
         ac(ne-1,nq-1) = ac(ne-1,nq-1) + 4.*ae(ne-1,nq-1)/(3.+deta) +
     $        an(ne-1,nq-1)
         aw(ne-1,nq-1) = aw(ne-1,nq-1) - ae(ne-1,nq-1)/(3.+deta)
         ae(ne-1,nq-1) = 0.
         an(ne-1,nq-1) = 0.
c
         call bicgst2d(ac,ae,aw,an,as,ddpsi2dv,rhs,rbh_eps
     $        ,rmax,ier,ne,nq)
         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 j = 1,nq
            ddpsi2dv(1,j) = ddpsi2dv(3,j)
            ddpsi2dv(ne,j) = (4.*ddpsi2dv(ne-1,j)-ddpsi2dv(ne-2,j))/(3.+
     &deta)
         enddo
         do i = 1,ne
            ddpsi2dv(i,1) = ddpsi2dv(i,2)
            ddpsi2dv(i,nq) = ddpsi2dv(i,nq-1)
         enddo
         do j = 1,nq
            do i = 1,ne
               ddpsi2dv(i,j) = ddpsi2dv(i,j)
               ddpsi2dv(i,j) = ddpsi2dv(i,j)
            enddo
         enddo
c
c
c
         do j = 1,nq
            do i = 1,ne
               dpsi2dv(i,j) = dpsi2dv(i,j) + ddpsi2dv(i,j)
            enddo
         enddo
         
c
      enddo
c
      write(*,*) 'It did not converge.'
      stop
c
 110  continue
      write(*,*) '--------------------------------------------------'
      print*,'Converge at Residual = ',rhsmax
c
c
c
      do j = 1,nq
         do i = 2,ne-1
            detapsisph(i,j)=0.5*(dpsi2dv(i+1,j)-dpsi2dv(i-1,j))/
     $           deta + sinh(0.5*etagrd(i))
         enddo
         detapsisph(1,j) = -detapsisph(3,j)
      enddo
c
      do j = 2,nq-1
         do i = 1,ne
            dqpsisph(i,j)=0.5*(dpsi2dv(i,j+1)-dpsi2dv(i,j-1))/dq
         enddo
      enddo
      do i = 1,ne
         dqpsisph(i,1) = -dqpsisph(i,2)
         dqpsisph(i,nq) = -dqpsisph(i,nq-1)
      enddo
      
c
      do j = 1,nq
         do i = 2,ne-1
            detaetapsisph(i,j)=(dpsi2dv(i+1,j)-2.*dpsi2dv(i,j)+dpsi2dv(i
     &-1,j))/
     $           deta**2 + sqrt(0.25)*cosh(0.5*etagrd(i))
         enddo
         detaetapsisph(1,j) = detaetapsisph(3,j)
      enddo
      
c
      do j = 2,nq-1
         do i = 1,ne
            detaqpsisph(i,j)=0.5*(detapsisph(i,j+1)-
     $           detapsisph(i,j-1))/dq
         enddo
      enddo
      do i = 1,ne
         detaqpsisph(i,1) = -detaqpsisph(i,2)
         detaqpsisph(i,nq) = -detaqpsisph(i,nq-1)
      enddo
c
      do j = 2,nq-1
         do i = 1,ne
            dqqpsisph(i,j)=0.5*(dqpsisph(i,j+1)-dqpsisph(i,j-1))/dq
         enddo
      enddo
      do i = 1,ne
         dqqpsisph(i,1) = dqqpsisph(i,2)
         dqqpsisph(i,nq) = dqqpsisph(i,nq-1)
      enddo
c
      do j = 1,nq
         psi2dv(:,j)=dpsi2dv(:,j)+2.0*cosh(0.5*etagrd)
      enddo
c
c
c
c
         
      eta = 0.5 * 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
               else
                  sign_eta(i,j,k) = 1
               endif
            enddo
         enddo
      enddo

     
c
c

      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,450,"IDAxiOddBrillBH.F","IDAxiOddBrillBH","Cann
     &ot create parameter table for interpolator")
      endif

      call CCTK_InterpHandle (interp_handle, "uniform cartesian")
      if (interp_handle .lt. 0) then
        call CCTK_Warn(0,455,"IDAxiOddBrillBH.F","IDAxiOddBrillBH","Cann
     &ot get handle for interpolation ! Forgot to activate an implementa
     &tion providing interpolation operators ??")
      endif

c
      coord_origin(1) = etagrd(1)
      coord_origin(2) = qgrd(1)
      coord_delta(1)  = etagrd(2) - etagrd(1)
      coord_delta(2)  = qgrd(2) - qgrd(1)

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

      in_array_dims(1) = ne; in_array_dims(2) = nq

      in_arrays(1) = CCTK_PointerTo(psi2dv)
      in_arrays(2) = CCTK_PointerTo(detapsisph)
      in_arrays(3) = CCTK_PointerTo(dqpsisph)
      in_arrays(4) = CCTK_PointerTo(detaetapsisph)
      in_arrays(5) = CCTK_PointerTo(detaqpsisph)
      in_arrays(6) = CCTK_PointerTo(dqqpsisph)

      out_arrays(1) = CCTK_PointerTo(psi2d)
      out_arrays(2) = CCTK_PointerTo(detapsi2d)
      out_arrays(3) = CCTK_PointerTo(dqpsi2d)
      out_arrays(4) = CCTK_PointerTo(detaetapsi2d)
      out_arrays(5) = CCTK_PointerTo(detaqpsi2d)
      out_arrays(6) = CCTK_PointerTo(dqqpsi2d)

      call CCTK_InterpLocalUniform (ierror, 2,
     $                              interp_handle, param_table_handle,
     $                              coord_origin, coord_delta,
     $                              npoints, type_codes(1), interp_coord
     &s,
     $                              6, in_array_dims, type_codes, in_arr
     &ays,
     $                              6, type_codes, out_arrays)

      psi = psi2d * exp(-0.5 * eta)
      detapsi2d = sign_eta * detapsi2d
      detaqpsi2d = sign_eta * detaqpsi2d
  
      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 = psi2d(i,j,k)
        o4 = 1/o3
        o5 = cos(phi(i,j,k))
        o6 = cos(q(i,j,k))
        o7 = dqpsi2d(i,j,k)
        o8 = -1.50000000000000d0*eta(i,j,k)
        o9 = exp(o8)
        o10 = detapsi2d(i,j,k)
        o11 = sin(q(i,j,k))
        o12 = sin(phi(i,j,k))
        psix(i,j,k) = o2*o4*(o10*o11*o5*o9 - 5.0000000000000d-1*o11*o3*o
     &
     &  5*o9 + o5*o6*o7*o9)
        psiy(i,j,k) = o2*o4*(o10*o11*o12*o9 - 5.0000000000000d-1*o11*o12
     &
     &  *o3*o9 + o12*o6*o7*o9)
        psiz(i,j,k) = o2*o4*(o10*o6*o9 - 5.0000000000000d-1*o3*o6*o9 - o
     &
     &  11*o7*o9)
               
c

        o1 = 5.0000000000000d-1*eta(i,j,k)
        o2 = exp(o1)
        o3 = psi2d(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 = detapsi2d(i,j,k)
        o10 = -2.50000000000000d0*eta(i,j,k)
        o11 = exp(o10)
        o12 = dqqpsi2d(i,j,k)
        o13 = sin(phi(i,j,k))
        o14 = o13**2
        o15 = detaqpsi2d(i,j,k)
        o16 = sin(q(i,j,k))
        o17 = dqpsi2d(i,j,k)
        o18 = detaetapsi2d(i,j,k)
        o19 = o16**2
        o20 = tan(q(i,j,k))
        o21 = 1/o20
        psixx(i,j,k) = o2*o4*(1.00000000000000d0*o11*o14*o17*o21 - 5.000
     &
     &  0000000000d-1*o11*o14*o3 + o11*o18*o19*o6 + 7.5000000000000d-1*o
     &
     &  11*o19*o3*o6 + 2.00000000000000d0*o11*o15*o16*o6*o7 - 3.00000000
     &
     &  000000d0*o11*o16*o17*o6*o7 + o11*o12*o6*o8 - 5.0000000000000d-1*
     &
     &  o11*o3*o6*o8 + o11*o14*o9 - 2.00000000000000d0*o11*o19*o6*o9 + o
     &
     &  11*o6*o8*o9)
        psixy(i,j,k) = o2*o4*(o11*o13*o18*o19*o5 - o11*o13*o17*o21*o5 + 
     &
     &  5.0000000000000d-1*o11*o13*o3*o5 + 7.5000000000000d-1*o11*o13*o1
     &
     &  9*o3*o5 + 2.00000000000000d0*o11*o13*o15*o16*o5*o7 - 3.000000000
     &
     &  00000d0*o11*o13*o16*o17*o5*o7 + o11*o12*o13*o5*o8 - 5.0000000000
     &
     &  000d-1*o11*o13*o3*o5*o8 - o11*o13*o5*o9 - 2.00000000000000d0*o11
     &
     &  *o13*o19*o5*o9 + o11*o13*o5*o8*o9)
        psixz(i,j,k) = o2*o4*(-(o11*o15*o19*o5) + 1.50000000000000d0*o11
     &
     &  *o17*o19*o5 - o11*o12*o16*o5*o7 + o11*o16*o18*o5*o7 + 1.25000000
     &
     &  000000d0*o11*o16*o3*o5*o7 + o11*o15*o5*o8 - 1.50000000000000d0*o
     &
     &  11*o17*o5*o8 - 3.00000000000000d0*o11*o16*o5*o7*o9)
        psiyy(i,j,k) = o2*o4*(o11*o14*o18*o19 + 7.5000000000000d-1*o11*o
     &
     &  14*o19*o3 + 1.00000000000000d0*o11*o17*o21*o6 - 5.0000000000000d
     &
     &  -1*o11*o3*o6 + 2.00000000000000d0*o11*o14*o15*o16*o7 - 3.0000000
     &
     &  0000000d0*o11*o14*o16*o17*o7 + o11*o12*o14*o8 - 5.0000000000000d
     &
     &  -1*o11*o14*o3*o8 - 2.00000000000000d0*o11*o14*o19*o9 + o11*o6*o9
     &
     &   + o11*o14*o8*o9)
        psiyz(i,j,k) = o2*o4*(-(o11*o13*o15*o19) + 1.50000000000000d0*o1
     &
     &  1*o13*o17*o19 - o11*o12*o13*o16*o7 + o11*o13*o16*o18*o7 + 1.2500
     &
     &  0000000000d0*o11*o13*o16*o3*o7 + o11*o13*o15*o8 - 1.500000000000
     &
     &  00d0*o11*o13*o17*o8 - 3.00000000000000d0*o11*o13*o16*o7*o9)
        psizz(i,j,k) = o2*o4*(o11*o12*o19 - 5.0000000000000d-1*o11*o19*o
     &
     &  3 - 2.00000000000000d0*o11*o15*o16*o7 + 3.00000000000000d0*o11*o
     &
     &  16*o17*o7 + o11*o18*o8 + 7.5000000000000d-1*o11*o3*o8 + o11*o19*
     &
     &  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.0*phi(i,j,k)
        o2 = cos(o1)
        o3 = -xi
        o4 = cos(q(i,j,k))
        o5 = o4*xi
        o6 = 1.0 + o3 + o5
        o7 = -eta0_me
        o8 = eta(i,j,k) + o7
        o9 = o8**2
        o10 = sigma_me**2
        o11 = 1/o10
        o12 = -(o11*o9)
        o13 = exp(o12)
        o14 = eta(i,j,k) + eta0_me
        o15 = o14**2
        o16 = -(o11*o15)
        o17 = exp(o16)
        o18 = o13 + o17
        o19 = sin(q(i,j,k))
        o20 = o19**n_me
        o21 = 2.0*amp_me*o18*o20*o6
        o22 = exp(o21)
        o23 = -1.0 + o22
        gxx(i,j,k) = 5.d-1*(1.0 + o22 + o2*o23)
        gxy(i,j,k) = o23*cos(phi(i,j,k))*sin(phi(i,j,k))
        gxz(i,j,k) = 0
        gyy(i,j,k) = 5.d-1*(1.0 + o22 - o2*o23)
        gyz(i,j,k) = 0
        gzz(i,j,k) = o22
            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)
               if (brandt_seidel == 1) then

        o1 = cos(phi(i,j,k))
        o2 = cos(q(i,j,k))
        o3 = o2**2
        o4 = -3.00000000000000d0*eta(i,j,k)
        o5 = exp(o4)
        o6 = sin(phi(i,j,k))
        o7 = sin(q(i,j,k))
        o8 = -1.00000000000000d0 + n
        o9 = o7**o8
        o10 = 2.00000000000000d0 + n
        o11 = o7**2
        o12 = -(o10*o11)
        o13 = 1.00000000000000d0 + n + o12
        o14 = o1**2
        o15 = o6**2
        o16 = o7**n
        o17 = -2.00000000000000d0 + n
        o18 = o7**o17
        kxx(i,j,k) = -2.00000000000000d0*gtil*o1*o13*o5*o6*o9 + 2.000000
     &
     &  00000000d0*dngtil*o1*o3*o5*o6*o9
        kxy(i,j,k) = gtil*o13*o14*o5*o9 - gtil*o13*o15*o5*o9 - dngtil*o1
     &
     &  4*o3*o5*o9 + dngtil*o15*o3*o5*o9
        kxz(i,j,k) = -(dngtil*o16*o2*o5*o6) - gtil*o13*o18*o2*o5*o6
        kyy(i,j,k) = 2.00000000000000d0*gtil*o1*o13*o5*o6*o9 - 2.0000000
     &
     &  0000000d0*dngtil*o1*o3*o5*o6*o9
        kyz(i,j,k) = dngtil*o1*o16*o2*o5 + gtil*o1*o13*o18*o2*o5
        kzz(i,j,k) = 0
	       else if (sergio == 1) then

        mass = sqrt(2.0 + sqrt(4.0 + byJ**2
     &  ))
        a = byJ/mass
        rBL = 5.0e-1*exp(eta(i,j,k))*(1.0 + (1.0*(a
     &   + mass)*exp(-eta(i,j,k)))/sqrt(-a**2 + mass**2))*sqrt(-a**2 + m
     &ass**2)
     &  *(1.0 + (1.*(a + mass)*exp(-eta(i,j,k)))/sqrt(a**2 + m
     &  ass**2))

        o1 = a**2
        o2 = a*o1
        o3 = -2.00000000000000d0*mass*rBL
        o4 = rBL**2
        o5 = o1 + o3 + o4
        o6 = sqrt(o5)
        o7 = cos(phi(i,j,k))
        o8 = cos(q(i,j,k))
        o9 = o8**2
        o10 = o1*o9
        o11 = o10 + o4
        o12 = o11**2
        o13 = 1/o12
        o14 = -3.00000000000000d0*eta(i,j,k)
        o15 = exp(o14)
        o16 = sin(phi(i,j,k))
        o17 = sin(q(i,j,k))
        o18 = o17**2
        o19 = o1 + o4
        o20 = 2.00000000000000d0*o19*o4
        o21 = -o1
        o22 = o21 + o4
        o23 = o11*o22
        o24 = o20 + o23
        o25 = o7**2
        o26 = o16**2
        o27 = o17*o18
        kxx(i,j,k) = -2.00000000000000d0*a*mass*o13*o15*o16*o18*o24*o7 +
     &
     &   4.0000000000000d0*mass*o13*o15*o16*o18*o2*o6*o7*o9*rBL
        kxy(i,j,k) = a*mass*o13*o15*o18*o24*o25 - a*mass*o13*o15*o18*o24
     &
     &  *o26 - 2.00000000000000d0*mass*o13*o15*o18*o2*o25*o6*o9*rBL + 2.
     &
     &  00000000000000d0*mass*o13*o15*o18*o2*o26*o6*o9*rBL
        kxz(i,j,k) = -(a*mass*o13*o15*o16*o17*o24*o8) - 2.00000000000000
     &
     &  d0*mass*o13*o15*o16*o2*o27*o6*o8*rBL
        kyy(i,j,k) = 2.00000000000000d0*a*mass*o13*o15*o16*o18*o24*o7 - 
     &
     &  4.0000000000000d0*mass*o13*o15*o16*o18*o2*o6*o7*o9*rBL
        kyz(i,j,k) = a*mass*o13*o15*o17*o24*o7*o8 + 2.00000000000000d0*m
     &
     &  ass*o13*o15*o2*o27*o6*o7*o8*rBL
        kzz(i,j,k) = 0
               else

        o1 = cos(phi(i,j,k))
        o2 = cos(q(i,j,k))
        o3 = o2**2
        o4 = -3.00000000000000d0*eta(i,j,k)
        o5 = exp(o4)
        o6 = sin(phi(i,j,k))
        o7 = sin(q(i,j,k))
        o8 = o7**2
        o9 = dnnnnngtil
        o10 = -3.6000000000000d1*o9
        o11 = dnnnnnnngtil
        o12 = -4.0000000000000d0*o11*o8
        o13 = 1.16000000000000d2*o8*o9
        o14 = dnnngtil
        o15 = -4.1500000000000d2*o8
        o16 = 1.50000000000000d2 + o15
        o17 = o14*o16
        o18 = dngtil
        o19 = 5.0000000000000d0*o8
        o20 = -2.00000000000000d0 + o19
        o21 = 7.5000000000000d1*o18*o20
        o22 = o10 + o12 + o13 + o17 + o21
        o23 = 6.0000000000000d1*byJ
        o24 = dnnnnnngtil
        o25 = 3.00000000000000d0*o24
        o26 = 2.00000000000000d0*q(i,j,k)
        o27 = cos(o26)
        o28 = 4.0000000000000d0*o24*o27
        o29 = 4.0000000000000d0*q(i,j,k)
        o30 = cos(o29)
        o31 = -7.0000000000000d0*o24*o30
        o32 = dnnnngtil
        o33 = 2.10000000000000d1*o32
        o34 = 6.4000000000000d1*o27*o32
        o35 = 2.03000000000000d2*o30*o32
        o36 = o23 + o25 + o28 + o31 + o33 + o34 + o35
        o37 = 4.0000000000000d0*o36
        o38 = 2.68000000000000d2*o27
        o39 = 5.8100000000000d2*o30
        o40 = 1.11000000000000d2 + o38 + o39
        o41 = dnngtil
        o42 = -5.0000000000000d0*o40*o41
        o43 = 2.00000000000000d1*o27
        o44 = 3.5000000000000d1*o30
        o45 = 9.0000000000000d0 + o43 + o44
        o46 = gtil
        o47 = 7.5000000000000d1*o45*o46
        o48 = o37 + o42 + o47
        o49 = o1**2
        o50 = o6**2
        o51 = o7*o8
        kxx(i,j,k) = -2.00000000000000d-1*o1*o22*o3*o5*o6*o8 - 2.5000000
     &
     &  0000000d-2*o1*o48*o5*o6*o8
        kxy(i,j,k) = 1.00000000000000d-1*o22*o3*o49*o5*o8 + 1.2500000000
     &
     &  0000d-2*o48*o49*o5*o8 - 1.00000000000000d-1*o22*o3*o5*o50*o8 - 1
     &
     &  .25000000000000d-2*o48*o5*o50*o8
        kxz(i,j,k) = 1.00000000000000d-1*o2*o22*o5*o51*o6 - 1.2500000000
     &
     &  0000d-2*o2*o48*o5*o6*o7
        kyy(i,j,k) = 2.00000000000000d-1*o1*o22*o3*o5*o6*o8 + 2.50000000
     &
     &  000000d-2*o1*o48*o5*o6*o8
        kyz(i,j,k) = -1.00000000000000d-1*o1*o2*o22*o5*o51 + 1.250000000
     &
     &  00000d-2*o1*o2*o48*o5*o7
        kzz(i,j,k) = 0
               endif
            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 j=2,nq-1
         adm=adm+(psi2dv(i,j)-(psi2dv(i+1,j)-psi2dv(i-1,j))/
     $        deta)*exp(0.5*etagrd(i))
      enddo
      
      adm=adm/(nq-2)
      print *,'ADM mass: ',adm

c
      print*,'Angular momentum parameter: a/m = ', byJ/(adm**2)
c
c
c
c
c
      
c
c
      
      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,rhs,qfetaeta,qfqq,
     $     ksq,psi2dv,dpsi2dv,ddpsi2dv,
     $     detapsisph,dqpsisph,detaetapsisph,detaqpsisph,dqqpsisph,
     $     etagrd,qgrd)
      
      return
      end
      
