














































c














































































 


































































































































































 






























































































































































































 





















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                



















































































c

      subroutine IDAxiBrillBH(cctk_dim,cctk_gsh,cctk_lsh,cctk_lbnd,cctk_
     &ubnd,cctk_tile_min,cctk_tile_max,cctk_ash,cctk_alignment,cctk_alig
     &nment_offset,cctk_from,cctk_to,cctk_bbox,cctk_delta_time,cctk_time
     &,cctk_delta_space,cctk_origin_space,cctk_levfac,cctk_levoff,cctk_l
     &evoffdenom,cctk_timefac,cctk_convlevel,cctk_convfac,cctk_nghostzon
     &es,cctk_iteration,cctkGH,cctk_ash1,cctk_ash2,cctk_ash3, alp,alp_p,
     &alp_p_p,betax,betax_p,betax_p_p,betay,betay_p,betay_p_p,betaz,beta
     &z_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_st
     &ate,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,kx
     &y,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)

      use cctk

      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 :: kxx (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kxx = kind(kxx)
      REAL*8 :: kxy (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kxy = kind(kxy)
      REAL*8 :: kxz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kxz = kind(kxz)
      REAL*8 :: kyy (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kyy = kind(kyy)
      REAL*8 :: kyz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kyz = kind(kyz)
      REAL*8 :: kzz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_kzz = kind(kzz)
      REAL*8 :: alp (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_alp = kind(alp)
      REAL*8 :: gxx (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gxx = kind(gxx)
      REAL*8 :: gxy (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gxy = kind(gxy)
      REAL*8 :: gxz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gxz = kind(gxz)
      REAL*8 :: gyy (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gyy = kind(gyy)
      REAL*8 :: gyz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gyz = kind(gyz)
      REAL*8 :: gzz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_gzz = kind(gzz)
      REAL*8 :: betax (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_betax = kind(betax)
      REAL*8 :: betay (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_betay = kind(betay)
      REAL*8 :: betaz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_betaz = kind(betaz)
      REAL*8, iNteNt(iN) :: r (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_r = kind(r)
      REAL*8, iNteNt(iN) :: x (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_x = kind(x)
      REAL*8, iNteNt(iN) :: y (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_y = kind(y)
      REAL*8, iNteNt(iN) :: z (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_z = kind(z)
      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 :: psiy (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psiy = kind(psiy)
      REAL*8 :: psiz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psiz = kind(psiz)
      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 :: 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 :: psizz (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_psizz = kind(psizz)
      INTEGER*4 :: conformal_state
      integer, parameter :: cctki_use_conformal_state = kind(conformal_s
     &tate)
      characTer*8, intent(IN) :: alp_p
      integer, parameter :: cctki_use_alp_p = kind(alp_p)
      characTer*8, intent(IN) :: alp_p_p
      integer, parameter :: cctki_use_alp_p_p = kind(alp_p_p)
      characTer*8, intent(IN) :: betax_p
      integer, parameter :: cctki_use_betax_p = kind(betax_p)
      characTer*8, intent(IN) :: betax_p_p
      integer, parameter :: cctki_use_betax_p_p = kind(betax_p_p)
      characTer*8, intent(IN) :: betay_p
      integer, parameter :: cctki_use_betay_p = kind(betay_p)
      characTer*8, intent(IN) :: betay_p_p
      integer, parameter :: cctki_use_betay_p_p = kind(betay_p_p)
      characTer*8, intent(IN) :: betaz_p
      integer, parameter :: cctki_use_betaz_p = kind(betaz_p)
      characTer*8, intent(IN) :: betaz_p_p
      integer, parameter :: cctki_use_betaz_p_p = kind(betaz_p_p)
      characTer*8, intent(IN) :: coarse_dx
      integer, parameter :: cctki_use_coarse_dx = kind(coarse_dx)
      characTer*8, intent(IN) :: coarse_dy
      integer, parameter :: cctki_use_coarse_dy = kind(coarse_dy)
      characTer*8, intent(IN) :: coarse_dz
      integer, parameter :: cctki_use_coarse_dz = kind(coarse_dz)
      characTer*8, intent(IN) :: dtalp
      integer, parameter :: cctki_use_dtalp = kind(dtalp)
      characTer*8, intent(IN) :: dtalp_p
      integer, parameter :: cctki_use_dtalp_p = kind(dtalp_p)
      characTer*8, intent(IN) :: dtalp_p_p
      integer, parameter :: cctki_use_dtalp_p_p = kind(dtalp_p_p)
      characTer*8, intent(IN) :: dtbetax
      integer, parameter :: cctki_use_dtbetax = kind(dtbetax)
      characTer*8, intent(IN) :: dtbetax_p
      integer, parameter :: cctki_use_dtbetax_p = kind(dtbetax_p)
      characTer*8, intent(IN) :: dtbetax_p_p
      integer, parameter :: cctki_use_dtbetax_p_p = kind(dtbetax_p_p)
      characTer*8, intent(IN) :: dtbetay
      integer, parameter :: cctki_use_dtbetay = kind(dtbetay)
      characTer*8, intent(IN) :: dtbetay_p
      integer, parameter :: cctki_use_dtbetay_p = kind(dtbetay_p)
      characTer*8, intent(IN) :: dtbetay_p_p
      integer, parameter :: cctki_use_dtbetay_p_p = kind(dtbetay_p_p)
      characTer*8, intent(IN) :: dtbetaz
      integer, parameter :: cctki_use_dtbetaz = kind(dtbetaz)
      characTer*8, intent(IN) :: dtbetaz_p
      integer, parameter :: cctki_use_dtbetaz_p = kind(dtbetaz_p)
      characTer*8, intent(IN) :: dtbetaz_p_p
      integer, parameter :: cctki_use_dtbetaz_p_p = kind(dtbetaz_p_p)
      characTer*8, intent(IN) :: dtlapse_state
      integer, parameter :: cctki_use_dtlapse_state = kind(dtlapse_state
     &)
      characTer*8, intent(IN) :: dtshift_state
      integer, parameter :: cctki_use_dtshift_state = kind(dtshift_state
     &)
      characTer*8, intent(IN) :: gxx_p
      integer, parameter :: cctki_use_gxx_p = kind(gxx_p)
      characTer*8, intent(IN) :: gxx_p_p
      integer, parameter :: cctki_use_gxx_p_p = kind(gxx_p_p)
      characTer*8, intent(IN) :: gxy_p
      integer, parameter :: cctki_use_gxy_p = kind(gxy_p)
      characTer*8, intent(IN) :: gxy_p_p
      integer, parameter :: cctki_use_gxy_p_p = kind(gxy_p_p)
      characTer*8, intent(IN) :: gxz_p
      integer, parameter :: cctki_use_gxz_p = kind(gxz_p)
      characTer*8, intent(IN) :: gxz_p_p
      integer, parameter :: cctki_use_gxz_p_p = kind(gxz_p_p)
      characTer*8, intent(IN) :: gyy_p
      integer, parameter :: cctki_use_gyy_p = kind(gyy_p)
      characTer*8, intent(IN) :: gyy_p_p
      integer, parameter :: cctki_use_gyy_p_p = kind(gyy_p_p)
      characTer*8, intent(IN) :: gyz_p
      integer, parameter :: cctki_use_gyz_p = kind(gyz_p)
      characTer*8, intent(IN) :: gyz_p_p
      integer, parameter :: cctki_use_gyz_p_p = kind(gyz_p_p)
      characTer*8, intent(IN) :: gzz_p
      integer, parameter :: cctki_use_gzz_p = kind(gzz_p)
      characTer*8, intent(IN) :: gzz_p_p
      integer, parameter :: cctki_use_gzz_p_p = kind(gzz_p_p)
      characTer*8, intent(IN) :: kxx_p
      integer, parameter :: cctki_use_kxx_p = kind(kxx_p)
      characTer*8, intent(IN) :: kxx_p_p
      integer, parameter :: cctki_use_kxx_p_p = kind(kxx_p_p)
      characTer*8, intent(IN) :: kxy_p
      integer, parameter :: cctki_use_kxy_p = kind(kxy_p)
      characTer*8, intent(IN) :: kxy_p_p
      integer, parameter :: cctki_use_kxy_p_p = kind(kxy_p_p)
      characTer*8, intent(IN) :: kxz_p
      integer, parameter :: cctki_use_kxz_p = kind(kxz_p)
      characTer*8, intent(IN) :: kxz_p_p
      integer, parameter :: cctki_use_kxz_p_p = kind(kxz_p_p)
      characTer*8, intent(IN) :: kyy_p
      integer, parameter :: cctki_use_kyy_p = kind(kyy_p)
      characTer*8, intent(IN) :: kyy_p_p
      integer, parameter :: cctki_use_kyy_p_p = kind(kyy_p_p)
      characTer*8, intent(IN) :: kyz_p
      integer, parameter :: cctki_use_kyz_p = kind(kyz_p)
      characTer*8, intent(IN) :: kyz_p_p
      integer, parameter :: cctki_use_kyz_p_p = kind(kyz_p_p)
      characTer*8, intent(IN) :: kzz_p
      integer, parameter :: cctki_use_kzz_p = kind(kzz_p)
      characTer*8, intent(IN) :: kzz_p_p
      integer, parameter :: cctki_use_kzz_p_p = kind(kzz_p_p)
      characTer*8, intent(IN) :: shift_state
      integer, parameter :: cctki_use_shift_state = kind(shift_state)
      REAL*8 amp
      integer, parameter :: cctki_use_amp = kind(amp)
      REAL*8 error_tolerance
      integer, parameter :: cctki_use_error_tolerance = kind(error_toler
     &ance)
      REAL*8 eta0
      integer, parameter :: cctki_use_eta0 = kind(eta0)
      REAL*8 etamax
      integer, parameter :: cctki_use_etamax = kind(etamax)
      REAL*8 sigma
      integer, parameter :: cctki_use_sigma = kind(sigma)
      integer*8 interpolator_name
      integer, parameter :: cctki_use_interpolator_name = kind(interpola
     &tor_name)
      integer*8 interpolator_pars
      integer, parameter :: cctki_use_interpolator_pars = kind(interpola
     &tor_pars)
      integer*8 output_psi2D_file_name
      integer, parameter :: cctki_use_output_psi2D_file_name = kind(outp
     &ut_psi2D_file_name)
      INTEGER*4 debug
      integer, parameter :: cctki_use_debug = kind(debug)
      INTEGER*4 debug_i
      integer, parameter :: cctki_use_debug_i = kind(debug_i)
      INTEGER*4 debug_ii
      integer, parameter :: cctki_use_debug_ii = kind(debug_ii)
      INTEGER*4 debug_j
      integer, parameter :: cctki_use_debug_j = kind(debug_j)
      INTEGER*4 debug_jj
      integer, parameter :: cctki_use_debug_jj = kind(debug_jj)
      INTEGER*4 debug_k
      integer, parameter :: cctki_use_debug_k = kind(debug_k)
      INTEGER*4 generate_StaticConformal_metric
      integer, parameter :: cctki_use_generate_StaticConformal_metric = 
     &kind(generate_StaticConformal_metric)
      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 ne
      integer, parameter :: cctki_use_ne = kind(ne)
      INTEGER*4 nq
      integer, parameter :: cctki_use_nq = kind(nq)
      INTEGER*4 output_psi2D
      integer, parameter :: cctki_use_output_psi2D = kind(output_psi2D)
      COMMON /IDAxiBrillBHpriv/amp, error_tolerance, eta0, etamax, sigma
     &, interpolator_name, interpolator_pars, output_psi2D_file_name, de
     &bug, debug_i, debug_ii, debug_j, debug_jj, debug_k, generate_Stati
     &cConformal_metric, interpolation_order, n, ne, nq, output_psi2D
      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
      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
      interface
      subroutine ConfToPhysInPlace (nx, ny, nz, psi, gxx, gxy, gxz, gyy,
     & gyz, gzz)
      implicit none
      INTEGER*4 nx
      INTEGER*4 ny
      INTEGER*4 nz
      REAL*8 psi(*)
      REAL*8 gxx(*)
      REAL*8 gxy(*)
      REAL*8 gxz(*)
      REAL*8 gyy(*)
      REAL*8 gyz(*)
      REAL*8 gzz(*)
      end subroutine ConfToPhysInPlace
      end interface

      interface
      subroutine mgparm(m,ifd59,id5,id9,idi,idg,imx,jmx)
      implicit REAL*8 (a-h,o-z)
      end
      subroutine mg5 (idim,ilower,iupper,jdim,jlower,jupper,
     &     cc,cn,cs,cw,ce,u,rhs,
     &     id5,id9,idi,idg,ifmg,eps,rmax,ier)
      implicit REAL*8 (a-h,o-z)
      REAL*8 cc(idim,jdim),cn(idim,jdim),cs(idim,jdim),cw(idim,jdim),
     &     ce(idim,jdim),cnw(idim,jdim),cne(idim,jdim),csw(idim,jdim),
     &     cse(idim,jdim),u(idim,jdim),rhs(idim,jdim)
      end
      end interface

      REAL*8, parameter :: one = 1
      REAL*8, parameter :: eps = 1.0e-10
      REAL*8 axibheps, rmax, dq, deta
      integer levels,id5,id9,idi,idg,ier
      REAL*8, allocatable :: cc(:,:),ce(:,:),cw(:,:),cn(:,:),cs(:,:),
     $     rhs(:,:),psi2d(:,:),detapsi2d(:,:),dqpsi2d(:,:),
     $     detaetapsi2d(:,:),detaqpsi2d(:,:),dqqpsi2d(:,:)
      REAL*8, allocatable :: etagrd(:),qgrd(:)
      REAL*8, allocatable :: eta(:,:,:),abseta(:,:,:),sign_eta(:,:,:),
     $     q(:,:,:),phi(:,:,:)
      REAL*8, allocatable :: psi2dv(:,:,:),detapsi2dv(:,:,:),
     $     dqpsi2dv(:,:,:),detaetapsi2dv(:,:,:),detaqpsi2dv(:,:,:),
     $     dqqpsi2dv(:,:,:)
      REAL*8 error_at_this_grid_point,max_error_in_grid
      REAL*8 o1,o2,o3,o4,o5,o6,o7,o8,o9
      REAL*8 o10,o11,o12,o13,o14,o15,o16,o17,o18,o19
      REAL*8 o20,o21,o22,o23,o24,o25,o26,o27,o28,o29
      REAL*8 o30,o31,o32,o33,o34,o35,o36,o37,o38,o39
      REAL*8 o40,o41,o42,o43,o44,o45,o46,o47,o48,o49
      REAL*8 o50,o51,o52,o53,o54,o55,o56,o57,o58,o59
      REAL*8 o60,o61,o62,o63,o64,o65,o66,o67,o68,o69
      REAL*8 o70,o71,o72,o73,o74,o75,o76,o77,o78,o79
      REAL*8 o80,o81,o82,o83,o84,o85,o86,o87,o88,o89
      REAL*8 o90,o91,o92,o93,o94,o95,o96
      integer i22
      REAL*8 pi
      REAL*8 adm
      REAL*8 exp_mhalf_eta, psi3d
      INTEGER*4 nx,ny,nz
      integer i,j,k,nquads
      integer npoints,ierror
      integer neb, nqb
      integer posn
      integer io_status

      integer, parameter :: max_string_length = 500

      integer                         param_table_handle, interp_handle
      character(max_string_length) :: message_buffer

      integer                         fstring_length
      character(max_string_length) :: interpolator_name_fstring
      character(max_string_length) :: interpolator_pars_fstring
      character(max_string_length) :: output_psi2D_file_name_fstring

      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


      pi = 4.0d0*atan(1.0d0)

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

c
c
c

c
c
c
c
c
c
c
      call CCTK_FortranString(fstring_length,
     $                        interpolator_name,
     $                        interpolator_name_fstring)
      call CCTK_InterpHandle(interp_handle,
     $                       interpolator_name_fstring(1:fstring_length)
     &)
      if (interp_handle .lt. 0) then
         call CCTK_Error(129,"IDAxiBrillBH.F","IDAxiBrillBH","Cannot get
     & interpolator handle!  did you forgot to activate a suitable local
     &-interpolator thorn?")
c
      endif

c
c
c
c
c
c
      call CCTK_FortranString(fstring_length,
     $                        interpolator_pars,
     $                        interpolator_pars_fstring)
      if (fstring_length .eq. 0) then
         interpolator_pars_fstring
     $         = "order=" // char(ichar('0') + interpolation_order)
      endif
      call Util_TableCreateFromString
     $            (param_table_handle,
     $             interpolator_pars_fstring(1:fstring_length))
      if (param_table_handle .lt. 0) then
        write(message_buffer, '(A,I8)')
     $      'failed to create interpolator param table: error code ',
     $      param_table_handle
        call CCTK_Warn(0,153,"IDAxiBrillBH.F","IDAxiBrillBH", message_bu
     &ffer)
      endif

c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
      neb = ne+2
      nqb = nq+2

      allocate(          cc(neb,nqb))
      allocate(          ce(neb,nqb))
      allocate(          cw(neb,nqb))
      allocate(          cn(neb,nqb))
      allocate(          cs(neb,nqb))
      allocate(         rhs(neb,nqb))
      allocate(       psi2d(neb,nqb))
      allocate(   detapsi2d(neb,nqb))
      allocate(     dqpsi2d(neb,nqb))
      allocate(detaetapsi2d(neb,nqb))
      allocate(  detaqpsi2d(neb,nqb))
      allocate(    dqqpsi2d(neb,nqb))

      allocate(etagrd(neb))
      allocate(  qgrd(nqb))

      allocate(          eta(nx,ny,nz))
      allocate(       abseta(nx,ny,nz))
      allocate(     sign_eta(nx,ny,nz))
      allocate(            q(nx,ny,nz))
      allocate(          phi(nx,ny,nz))
      allocate(       psi2dv(nx,ny,nz))
      allocate(   detapsi2dv(nx,ny,nz))
      allocate(     dqpsi2dv(nx,ny,nz))
      allocate(detaetapsi2dv(nx,ny,nz))
      allocate(  detaqpsi2dv(nx,ny,nz))
      allocate(    dqqpsi2dv(nx,ny,nz))

c
      psi2d = 1.0d0
      detapsi2d = 0.0d0

      nquads = 2
      dq = nquads*0.5d0*pi/(nqb-2)
      deta = etamax/(neb-3)

      do j=1,nqb
        qgrd(j) = (j-1.5d0)*dq
        do i=1,neb
          etagrd(i) = (i-2)*deta

        o1 = dq**2
        o2 = 1/o1
        o3 = 1/dq
        o4 = tan(qgrd(j))
        o5 = 1/o4
        o6 = deta**2
        o7 = 1/o6
        o8 = sin(qgrd(j))
        o9 = o8**2
        o10 = 1/o9
        o11 = cos(qgrd(j))
        o12 = o11**2
        o13 = etagrd(i)**2
        o14 = sigma**2
        o15 = 1/o14
        o16 = -(o13*o15)
        o17 = -2.00000000000000d0*etagrd(i)*eta0*o15
        o18 = eta0**2
        o19 = -(o15*o18)
        o20 = o16 + o17 + o19
        o21 = exp(o20)
        i22 = -2 + n
c
c
c
c
c
c
        o23 = o8**i22
        o24 = n**2
        o25 = 2.00000000000000d0*etagrd(i)*eta0*o15
        o26 = o16 + o19 + o25
        o27 = exp(o26)
        o28 = o8**n
        o29 = o14**2
        o30 = 1/o29
        o31 = o4**2
        o32 = 1/o31
        o33 = 5.0000000000000d-1*etagrd(i)
        o34 = cosh(o33)
        Cn(i,j) = o2 + 5.0000000000000d-1*o3*o5
        Cs(i,j) = o2 - 5.0000000000000d-1*o3*o5
        Cw(i,j) = o7
        Cc(i,j) = -1.25000000000000d-1 - 1.25000000000000d-1*o10 - 2.000
     &
     &  00000000000d0*o2 - 2.50000000000000d-1*amp*n*o12*o21*o23 + 2.500
     &
     &  00000000000d-1*amp*o12*o21*o23*o24 - 2.50000000000000d-1*amp*n*o
     &
     &  12*o23*o27 + 2.50000000000000d-1*amp*o12*o23*o24*o27 - 2.5000000
     &
     &  0000000d-1*amp*n*o21*o28 - 5.0000000000000d-1*amp*o15*o21*o28 - 
     &
     &  2.50000000000000d-1*amp*n*o27*o28 - 5.0000000000000d-1*amp*o15*o
     &
     &  27*o28 + 2.00000000000000d0*amp*etagrd(i)*eta0*o21*o28*o30 + amp
     &*o13*o
     &  21*o28*o30 + amp*o18*o21*o28*o30 - 2.00000000000000d0*amp*etagrd
     &(i)*et
     &  a0*o27*o28*o30 + amp*o13*o27*o28*o30 + amp*o18*o27*o28*o30 + 1.2
     &
     &  5000000000000d-1*o32 - 2.00000000000000d0*o7
        Ce(i,j) = o7
        Rhs(i,j) = -2.50000000000000d-1*o34 + 2.50000000000000d-1*o10*o3
     &
     &  4 + 5.0000000000000d-1*amp*n*o12*o21*o23*o34 - 5.0000000000000d-
     &
     &  1*amp*o12*o21*o23*o24*o34 + 5.0000000000000d-1*amp*n*o12*o23*o27
     &
     &  *o34 - 5.0000000000000d-1*amp*o12*o23*o24*o27*o34 + 5.0000000000
     &
     &  000d-1*amp*n*o21*o28*o34 + amp*o15*o21*o28*o34 + 5.0000000000000
     &
     &  d-1*amp*n*o27*o28*o34 + amp*o15*o27*o28*o34 - 4.0000000000000d0*
     &
     &  amp*etagrd(i)*eta0*o21*o28*o30*o34 - 2.00000000000000d0*amp*o13*
     &o21*o2
     &  8*o30*o34 - 2.00000000000000d0*amp*o18*o21*o28*o30*o34 + 4.00000
     &
     &  00000000d0*amp*etagrd(i)*eta0*o27*o28*o30*o34 - 2.00000000000000
     &d0*amp
     &  *o13*o27*o28*o30*o34 - 2.00000000000000d0*amp*o18*o27*o28*o30*o3
     &
     &  4 - 2.50000000000000d-1*o32*o34

        enddo
      enddo
c
      do j=1,nqb
        ce(2,j)=ce(2,j)+cw(2,j)
        cw(2,j)=0.0d0

        cw(neb-1,j)=cw(neb-1,j)+ce(neb-1,j)
        cc(neb-1,j)=cc(neb-1,j)-deta*ce(neb-1,j)
        ce(neb-1,j)=0.0d0

      enddo
      do i=1,neb
        cc(i,2)=cc(i,2)+cs(i,2)
        cs(i,2)=0.0d0
        cc(i,nqb-1)=cc(i,nqb-1)+cn(i,nqb-1)
        cn(i,nqb-1)=0.0d0
      enddo

c
      axibheps = error_tolerance
      call CCTK_Info("IDAxiBrillBH","Calling axisymmetric solver")

      call mgparm (levels,5,id5,id9,idi,idg,neb,nqb)
      call mg5 (neb,2,neb-1,nqb,2,nqb-1,
     $          cc,cn,cs,cw,ce,psi2d,rhs,
     $          id5,id9,idi,idg,1,axibheps,rmax,ier)

      call CCTK_Info("IDAxiBrillBH","Solve complete")

c
c
c
      if(ier .ne. 0) then
         write(message_buffer, '(A,I8)')
     $      'failed to solve elliptic equation: ier=', ier
         call CCTK_Warn(0,254,"IDAxiBrillBH.F","IDAxiBrillBH", message_b
     &uffer)
      end if

      print *,'rmax = ',rmax
      print *,'axibheps = ',axibheps
      print *,'psi2d = ',maxval(psi2d),' ',minval(psi2d)

      max_error_in_grid = 0.0d0
      do j=2,nqb-1
        do i=2,neb-1
          error_at_this_grid_point = rhs(i,j)
     $                               - psi2d(i,j  )*cc(i,j)
     $                               - psi2d(i,j+1)*cn(i,j)
     $                               - psi2d(i,j-1)*cs(i,j)
     $                               - psi2d(i+1,j)*ce(i,j)
     $                               - psi2d(i-1,j)*cw(i,j)
          max_error_in_grid = max(max_error_in_grid,
     $                            abs(error_at_this_grid_point))
        enddo
      enddo
      print *,'Resulting eps =',max_error_in_grid

c
      do j=1,nqb
         psi2d(1  ,j) = psi2d(3,j)
         psi2d(neb,j) = - deta*psi2d(neb-1,j) + psi2d(neb-2,j)
      enddo

      do i=1,neb
        psi2d(i,1  ) = psi2d(i,2)
        psi2d(i,nqb) = psi2d(i,nqb-1)
      enddo

c
      do j=2,nqb-1
         do i=2,neb-1
            dqpsi2d  (i,j) = 0.5d0*(psi2d(i,j+1)-psi2d(i,j-1))/dq
            dqqpsi2d (i,j) = (psi2d(i,j+1)+psi2d(i,j-1)-2.0d0*psi2d(i,j)
     &)/dq**2
            detapsi2d(i,j) = sinh(0.5d0*etagrd(i))
     $                       + 0.5d0*(psi2d(i+1,j)-psi2d(i-1,j))/deta
            detaetapsi2d(i,j)
     $          = 0.5d0*cosh(0.5d0*etagrd(i))
     $            + (psi2d(i+1,j)+psi2d(i-1,j)-2.0d0*psi2d(i,j))/deta**2
     &
         enddo
      enddo

      do j=1,nqb
        detapsi2d(1,j)=-detapsi2d(3,j)
        detapsi2d(neb,j)=detapsi2d(neb-2,j) ! simplified

        detaetapsi2d(1,j)=detaetapsi2d(3,j)
        detaetapsi2d(neb,j)=detaetapsi2d(neb-2,j) ! simplified...

        dqqpsi2d(1,j)=dqqpsi2d(3,j)
        dqqpsi2d(neb,j)=dqqpsi2d(neb-2,j) ! simplified

        dqpsi2d(1,j)=dqpsi2d(3,j)
        dqpsi2d(neb,j)=-dq*dqpsi2d(neb-1,j)+dqpsi2d(neb-2,j)
      enddo

      do i=1,neb
        detapsi2d(i,1)=detapsi2d(i,2)
        detapsi2d(i,nqb)=detapsi2d(i,nqb-1)

        detaetapsi2d(i,1)=detaetapsi2d(i,2)
        detaetapsi2d(i,nqb)=detaetapsi2d(i,nqb-1)

        dqqpsi2d(i,1)=dqqpsi2d(i,2)
        dqqpsi2d(i,nqb)=dqqpsi2d(i,nqb-1)

        dqpsi2d(i,1)=-dqpsi2d(i,2)
        dqpsi2d(i,nqb)=-dqpsi2d(i,nqb-1)
      enddo

      do j=2,nqb-1
        do i=2,neb-1
          detaqpsi2d(i,j)=0.5d0*(detapsi2d(i,j+1)-detapsi2d(i,j-1))/dq
        enddo
      enddo

      do j=1,nqb
        detaqpsi2d(1,j)=-detaqpsi2d(3,j)
        detaqpsi2d(neb,j)=detaqpsi2d(neb-2,j) ! simplified
      enddo

      do i=1,ne
        detaqpsi2d(i,1)=-detaqpsi2d(i,2)
        detaqpsi2d(i,nqb)=-detaqpsi2d(i,nqb-1)
      enddo

      do j=1,nqb
        psi2d(:,j)=psi2d(:,j)+2.0d0*cosh(0.5d0*etagrd)
      enddo

      if (debug .ge. 6) then
         print *, '### 2-D grid results (= inputs to interpolation) ###'
     &
         print *, 'effective 2-D grid size: neb,nqb =', neb,nqb
         print *, 'debug_{ii,jj} =', debug_ii,debug_jj
         print *, 'at this 2-D grid point...'
         print *, '             eta =', etagrd(debug_ii)
         print *, '       theta [q] =',   qgrd(debug_jj)
         print *, '           psi2d =',        psi2d(debug_ii,debug_jj)
         print *, '       detapsi2d =',    detapsi2d(debug_ii,debug_jj)
         print *, '         dqpsi2d =',      dqpsi2d(debug_ii,debug_jj)
         print *, '    detaetapsi2d =', detaetapsi2d(debug_ii,debug_jj)
         print *, '      detaqpsi2d =',   detaqpsi2d(debug_ii,debug_jj)
         print *, '        dqqpsi2d =',     dqqpsi2d(debug_ii,debug_jj)
      endif

c
c
c
      if (output_psi2D .ne. 0) then
         write (message_buffer, '(A,A,A)')
     $         'writing 2D psi to "',
     $            output_psi2D_file_name_fstring(1:fstring_length),
     $         '"'
         call CCTK_Info("IDAxiBrillBH",message_buffer)

         call CCTK_FortranString(fstring_length,
     $                           output_psi2D_file_name,
     $                           output_psi2D_file_name_fstring)
         open (9, iostat=io_status, status='replace',
     $         file=output_psi2D_file_name_fstring)
         if (io_status .ne. 0) then 
            write (message_buffer, '(A,A,A,I8)')
     $            'error opening psi2D output file "',
     $            output_psi2D_file_name_fstring(1:fstring_length),
     $            '": io_status=', io_status
            call CCTK_Warn(0,383,"IDAxiBrillBH.F","IDAxiBrillBH", messag
     &e_buffer)
         endif

         write (9, '(a)')
     $         '# eta          theta (=q)     psi (2D)            psi (3
     &D)'
            do i = 1,neb
            exp_mhalf_eta = exp(-0.5d0 * etagrd(i))
               do j = 1,nqb
               psi3d = psi2d(i,j) * exp_mhalf_eta
               write (9, '(f12.8,a, f12.8,a, g20.10e3,a, g20.10e3)')
     $               etagrd(i), ' ', qgrd(j), ' ',
     $               psi2d(i,j), ' ', psi3d
               end do
            write (9, '(a)') ''
            end do

         close (9, iostat=io_status)
         if (io_status .ne. 0) then 
            write(message_buffer, '(A,I8)')
     $         'error closing psi2D output file: io_status=', io_status
            call CCTK_Warn(0,403,"IDAxiBrillBH.F","IDAxiBrillBH", messag
     &e_buffer)
         endif
      endif



c
c
c
c
c
c
c
c
c
c
c
c

      call CCTK_Info("IDAxiBrillBH","interpolating solution to xyz grid 
     &points")

      eta = log (r + eps)
      abseta = abs (eta)
      sign_eta = sign (one, eta)
      q = atan2 (sqrt ((x + eps)**2 + y**2), z)
      phi = atan2 (y, x + eps)

      if (debug .ge. 6) then
c
c
	 posn = (debug_i-1) + nx*(debug_j-1) + nx*ny*(debug_k-1)
	 print *, '### 3-D interpolation coordinates and inputs ###'
	 print *, '3-D grid size: n[xyz] =', nx,ny,nz
	 print *, 'debug_[ijk] =', debug_i,debug_j,debug_k, '==> posn =', posn
         print *, 'at this 3-D grid point, ...'
	 print *, '        x =',      x(debug_i,debug_j,debug_k)
	 print *, '        y =',      y(debug_i,debug_j,debug_k)
	 print *, '        z =',      z(debug_i,debug_j,debug_k)
	 print *, '      eta =',    eta(debug_i,debug_j,debug_k)
	 print *, '   abseta =', abseta(debug_i,debug_j,debug_k)
	 print *, 'theta [q] =',      q(debug_i,debug_j,debug_k)
	 print *, '      phi =',    phi(debug_i,debug_j,debug_k)
      endif

c
      npoints = nx*ny*nz

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

      if (debug .ge. 6) then
         print *, '### 2-D grid origin:   eta=', coord_origin(1)
         print *, '                     theta=', coord_origin(2)
         print *, '              delta:   eta=', coord_delta(1)
         print *, '                     theta=', coord_delta(2)
      end if

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

      in_array_dims(1) = neb
      in_array_dims(2) = nqb

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

      out_arrays(1) = CCTK_PointerTo(       psi2dv)
      out_arrays(2) = CCTK_PointerTo(   detapsi2dv)
      out_arrays(3) = CCTK_PointerTo(     dqpsi2dv)
      out_arrays(4) = CCTK_PointerTo(detaetapsi2dv)
      out_arrays(5) = CCTK_PointerTo(  detaqpsi2dv)
      out_arrays(6) = CCTK_PointerTo(    dqqpsi2dv)

             psi2dv = 1
         detapsi2dv = 0
           dqpsi2dv = 0
      detaetapsi2dv = 0
        detaqpsi2dv = 0
          dqqpsi2dv = 0

      call CCTK_InterpLocalUniform (ierror, 2,
     $                              interp_handle, param_table_handle,
     $                              coord_origin, coord_delta,
     $                              npoints, int(type_codes(1)), interp_
     &coords,
     $                              6, in_array_dims, type_codes, in_arr
     &ays,
     $                              6, type_codes, out_arrays)
      if (ierror < 0) then
        write(message_buffer, '(A,I8)')
     $      'error in interpolator: ierror=', ierror
        call CCTK_Warn(1,498,"IDAxiBrillBH.F","IDAxiBrillBH", message_bu
     &ffer)
      endif

      call Util_TableDestroy (ierror, param_table_handle)

      if (debug .ge. 6) then
         print *, '### interpolation results (at this 3-D grid point) ##
     &#'
	 print *, '   psi2dv =', psi2dv(debug_i,debug_j,debug_k)
      end if
c
c
c
c
      psi = psi2dv * exp (-0.5d0 * eta)
      detapsi2dv = sign_eta * detapsi2dv
      detaqpsi2dv = sign_eta * detaqpsi2dv

      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 = psi2dv(i,j,k)
        o4 = 1/o3
        o5 = cos(phi(i,j,k))
        o6 = cos(q(i,j,k))
        o7 = dqpsi2dv(i,j,k)
        o8 = -1.50000000000000d0*eta(i,j,k)
        o9 = exp(o8)
        o10 = detapsi2dv(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 = psi2dv(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 = detapsi2dv(i,j,k)
        o10 = -2.50000000000000d0*eta(i,j,k)
        o11 = exp(o10)
        o12 = dqqpsi2dv(i,j,k)
        o13 = sin(phi(i,j,k))
        o14 = o13**2
        o15 = detaqpsi2dv(i,j,k)
        o16 = sin(q(i,j,k))
        o17 = dqpsi2dv(i,j,k)
        o18 = detaetapsi2dv(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

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

        o1 = 2.00000000000000d0*phi(i,j,k)
        o2 = cos(o1)
        o3 = -eta0
        o4 = eta(i,j,k) + o3
        o5 = o4**2
        o6 = sigma**2
        o7 = 1/o6
        o8 = -(o5*o7)
        o9 = exp(o8)
        o10 = eta(i,j,k) + eta0
        o11 = o10**2
        o12 = -(o11*o7)
        o13 = exp(o12)
        o14 = o13 + o9
        o15 = sin(q(i,j,k))
        o16 = o15**n
        o17 = 2.00000000000000d0*amp*o14*o16
        o18 = exp(o17)
        o19 = -1.00000000000000d0 + o18
        o20 = eta(i,j,k)**2
        o21 = 2.00000000000000d0*o20
        o22 = eta0**2
        o23 = 2.00000000000000d0*o22
        o24 = eta(i,j,k)*o6
        o25 = o21 + o23 + o24
        o26 = -(o25*o7)
        o27 = exp(o26)
        o28 = 1/o15
        o29 = o20 + o22
        o30 = 2.00000000000000d0*o29*o7
        o31 = exp(o30)
        o32 = sin(phi(i,j,k))
        o33 = sin(o1)
        o34 = o11*o7
        o35 = exp(o34)
        o36 = cos(phi(i,j,k))
        o37 = o32**2
        o38 = cos(q(i,j,k))
        o39 = o38**2
        o40 = o5*o7
        o41 = exp(o40)
        o42 = o15**2
        o43 = o36**2
        o44 = -2.00000000000000d0*o19*o31*o6
        o45 = -2.00000000000000d0*eta(i,j,k)
        o46 = -2.00000000000000d0*eta0
        o47 = n*o6
        o48 = 4.0000000000000d0*eta(i,j,k)*eta0*o7
        o49 = exp(o48)
        o50 = -2.00000000000000d0*eta(i,j,k)*o49
        o51 = 2.00000000000000d0*eta0*o49
        o52 = n*o49*o6
        o53 = 2.00000000000000d0*q(i,j,k)
        o54 = cos(o53)
        o55 = 2.00000000000000d0*eta(i,j,k)
        o56 = 2.00000000000000d0*eta0
        o57 = 2.00000000000000d0*eta(i,j,k)*o49
        o58 = -2.00000000000000d0*eta0*o49
        o59 = o47 + o52 + o55 + o56 + o57 + o58
        o60 = o54*o59
        o61 = o45 + o46 + o47 + o50 + o51 + o52 + o60
        o62 = o17 + o40
        o63 = exp(o62)
        o64 = amp*o16*o61*o63
        o65 = o44 + o64
        o66 = -1.00000000000000d0 + o49
        o67 = -2.00000000000000d0*eta0*o66
        o68 = 1.00000000000000d0 + o49
        o69 = 2.00000000000000d0*eta(i,j,k)*o68
        o70 = n*o6*o68
        o71 = o67 + o69 + o70
        o72 = o56 + o6
        o73 = eta(i,j,k)*o72
        o74 = o20 + o22 + o73
        o75 = -(o7*o74)
        o76 = o17 + o75
        o77 = exp(o76)
        o78 = -eta(i,j,k)
        o79 = exp(o78)
        o80 = -2.00000000000000d0*o29*o7
        o81 = o17 + o80
        o82 = exp(o81)
        o83 = -1.00000000000000d0 + n
        o84 = o15**o83
        o85 = o35 + o41
        o86 = -(n*o39*o6*o85)
        o87 = eta(i,j,k)*o49
        o88 = -(eta0*o49)
        o89 = eta(i,j,k) + eta0 + o87 + o88
        o90 = 2.00000000000000d0*o41*o42*o89
        o91 = o86 + o90
        o92 = o17 + o26
        o93 = exp(o92)
        o94 = n*o39*o6*o85
        o95 = -2.00000000000000d0*o41*o42*o89
        o96 = o94 + o95
        gxx(i,j,k) = 5.0000000000000d-1*(1.00000000000000d0 + o18 + o19*
     &
     &  o2)
c
c
c
c
c
c
c
c
c
c
        gxy(i,j,k) = 5.0000000000000d-1*o19*o33
c
c
c
c
c
        gxz(i,j,k) = 0
c
c
c
        gyy(i,j,k) = 5.0000000000000d-1*(1.00000000000000d0 + o18 - o19*
     &
     &  o2)
c
c
c
c
c
        gyz(i,j,k) = 0
c
c
c
        gzz(i,j,k) = o18
c
c
c

               if (r(i,j,k) < 1.0d-10) then
                  
                  psi(i,j,k) = 1
                  
                  psix(i,j,k) = 0
                  psiy(i,j,k) = 0
                  psiz(i,j,k) = 0
                  
                  psixx(i,j,k) = 0
                  psixy(i,j,k) = 0
                  psixz(i,j,k) = 0
                  psiyy(i,j,k) = 0
                  psiyz(i,j,k) = 0
                  psizz(i,j,k) = 0
                  
                  gxx(i,j,k) = 1
                  gxy(i,j,k) = 0
                  gxz(i,j,k) = 0
                  gyy(i,j,k) = 1
                  gyz(i,j,k) = 0
                  gzz(i,j,k) = 1
                  
               end if

            enddo
         enddo
      enddo

c
      if (generate_StaticConformal_metric .eq. 0) then
         call CCTK_Info("IDAxiBrillBH","converting to physical metric")
         call ConfToPhysInPlace(nx, ny, nz,
     $				psi,
     $				gxx, gxy, gxz,
     $				     gyy, gyz,
     $				          gzz)
c
         conformal_state = 0
      else
c
         conformal_state = 3
      end if

c
      kxx = 0.0d0
      kxy = 0.0d0
      kxz = 0.0d0
      kyy = 0.0d0
      kyz = 0.0d0
      kzz = 0.0d0

      if (debug .ge. 6) then
         print *, '### final results (again at this 3-D grid point) ###'
     &
         if (conformal_state .gt. 0) then
            print *, '### ... conformal metric with'
	    print *, '   psi =', psi(debug_i,debug_j,debug_k)
         else
            print *, '### ... physical metric with'
         end if
	 print *, '   gxx =', gxx(debug_i,debug_j,debug_k)
	 print *, '   gxy =', gxy(debug_i,debug_j,debug_k)
	 print *, '   gxz =', gxz(debug_i,debug_j,debug_k)
	 print *, '   gyy =', gyy(debug_i,debug_j,debug_k)
	 print *, '   gyz =', gyz(debug_i,debug_j,debug_k)
	 print *, '   gzz =', gzz(debug_i,debug_j,debug_k)
      endif

c
c
c

c
      call CCTK_Info("IDAxiBrillBH","computing ADM mass")
      i = neb-15
      adm = 0.0d0
      do j=2,nqb-1
        adm = adm
     $        + (psi2d(i,j)-(psi2d(i+1,j)-psi2d(i-1,j))/deta)
     $          *exp(0.5d0*etagrd(i))
      enddo
      adm=adm/(nqb-2)
      print *,'ADM mass: ',adm

      if ((CCTK_Equals(initial_lapse,"schwarz") .ne. 0)) then
         write (*,*)"Initial with schwarzschild-like lapse"
         write (*,*)"using alp = (2r - adm)/(2r+adm)"
         alp = (2.0d0*r - adm)/(2.0d0*r+adm)
      endif

      deallocate(cc,ce,cw,cn,cs,rhs,psi2d,detapsi2d,dqpsi2d,
     $     detaetapsi2d,detaqpsi2d,dqqpsi2d,
     $     etagrd,qgrd,
     $     eta,abseta,sign_eta,q,phi,psi2dv,detapsi2dv,dqpsi2dv,
     $     detaetapsi2dv,detaqpsi2dv,dqqpsi2dv)

      return
      end
