






















































































































 


































































































































































 





































































































































































































 






















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                



















































































      subroutine setupbrilldata3D(cctk_dim,cctk_gsh,cctk_lsh,cctk_lbnd,c
     &ctk_ubnd,cctk_tile_min,cctk_tile_max,cctk_ash,cctk_alignment,cctk_
     &alignment_offset,cctk_from,cctk_to,cctk_bbox,cctk_delta_time,cctk_
     &time,cctk_delta_space,cctk_origin_space,cctk_levfac,cctk_levoff,cc
     &tk_levoffdenom,cctk_timefac,cctk_convlevel,cctk_convfac,cctk_nghos
     &tzones,cctk_iteration,cctkGH,cctk_ash1,cctk_ash2,cctk_ash3, brillM
     &linear,brillNsource,brillpsi, alp,alp_p,alp_p_p,betax,betax_p,beta
     &x_p_p,betay,betay_p,betay_p_p,betaz,betaz_p,betaz_p_p,coarse_dx,co
     &arse_dy,coarse_dz,conformal_state,dtalp,dtalp_p,dtalp_p_p,dtbetax,
     &dtbetax_p,dtbetax_p_p,dtbetay,dtbetay_p,dtbetay_p_p,dtbetaz,dtbeta
     &z_p,dtbetaz_p_p,dtlapse_state,dtshift_state,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,kxy,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)

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

c


      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 brillMlinear (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_brillMlinear = kind(brillMlinear)
      REAL*8 brillNsource (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_brillNsource = kind(brillNsource)
      REAL*8 brillpsi (cctk_ash1,cctk_ash2,cctk_ash3)
      integer, parameter :: cctki_use_brillpsi = kind(brillpsi)
      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 brill3d_d
      integer, parameter :: cctki_use_brill3d_d = kind(brill3d_d)
      REAL*8 brill3d_e
      integer, parameter :: cctki_use_brill3d_e = kind(brill3d_e)
      REAL*8 brill3d_m
      integer, parameter :: cctki_use_brill3d_m = kind(brill3d_m)
      REAL*8 brill3d_n
      integer, parameter :: cctki_use_brill3d_n = kind(brill3d_n)
      REAL*8 brill3d_phi0
      integer, parameter :: cctki_use_brill3d_phi0 = kind(brill3d_phi0)
      REAL*8 const_v0
      integer, parameter :: cctki_use_const_v0 = kind(const_v0)
      REAL*8 eppley_a
      integer, parameter :: cctki_use_eppley_a = kind(eppley_a)
      REAL*8 eppley_r0
      integer, parameter :: cctki_use_eppley_r0 = kind(eppley_r0)
      REAL*8 eppley_sigmar
      integer, parameter :: cctki_use_eppley_sigmar = kind(eppley_sigmar
     &)
      REAL*8 eppley_sigmarho
      integer, parameter :: cctki_use_eppley_sigmarho = kind(eppley_sigm
     &arho)
      REAL*8 exp_a
      integer, parameter :: cctki_use_exp_a = kind(exp_a)
      REAL*8 exp_rho0
      integer, parameter :: cctki_use_exp_rho0 = kind(exp_rho0)
      REAL*8 exp_sigmaz
      integer, parameter :: cctki_use_exp_sigmaz = kind(exp_sigmaz)
      REAL*8 gundlach_a
      integer, parameter :: cctki_use_gundlach_a = kind(gundlach_a)
      REAL*8 gundlach_r0
      integer, parameter :: cctki_use_gundlach_r0 = kind(gundlach_r0)
      REAL*8 gundlach_sigmar
      integer, parameter :: cctki_use_gundlach_sigmar = kind(gundlach_si
     &gmar)
      REAL*8 gundlach_sigmarho
      integer, parameter :: cctki_use_gundlach_sigmarho = kind(gundlach_
     &sigmarho)
      REAL*8 rhofudge
      integer, parameter :: cctki_use_rhofudge = kind(rhofudge)
      REAL*8 robin_inf
      integer, parameter :: cctki_use_robin_inf = kind(robin_inf)
      REAL*8 thresh
      integer, parameter :: cctki_use_thresh = kind(thresh)
      integer*8 bound
      integer, parameter :: cctki_use_bound = kind(bound)
      integer*8 q_function
      integer, parameter :: cctki_use_q_function = kind(q_function)
      integer*8 solver
      integer, parameter :: cctki_use_solver = kind(solver)
      INTEGER*4 eppley_b
      integer, parameter :: cctki_use_eppley_b = kind(eppley_b)
      INTEGER*4 eppley_c
      integer, parameter :: cctki_use_eppley_c = kind(eppley_c)
      INTEGER*4 exp_b
      integer, parameter :: cctki_use_exp_b = kind(exp_b)
      INTEGER*4 gundlach_b
      integer, parameter :: cctki_use_gundlach_b = kind(gundlach_b)
      INTEGER*4 gundlach_c
      integer, parameter :: cctki_use_gundlach_c = kind(gundlach_c)
      INTEGER*4 output_coeffs
      integer, parameter :: cctki_use_output_coeffs = kind(output_coeffs
     &)
      INTEGER*4 robin_falloff
      integer, parameter :: cctki_use_robin_falloff = kind(robin_falloff
     &)
      INTEGER*4 sor_maxit
      integer, parameter :: cctki_use_sor_maxit = kind(sor_maxit)
      COMMON /IDBrillDatapriv/brill3d_d, brill3d_e, brill3d_m, brill3d_n
     &, brill3d_phi0, const_v0, eppley_a, eppley_r0, eppley_sigmar, eppl
     &ey_sigmarho, exp_a, exp_rho0, exp_sigmaz, gundlach_a, gundlach_r0,
     & gundlach_sigmar, gundlach_sigmarho, rhofudge, robin_inf, thresh, 
     &bound, q_function, solver, eppley_b, eppley_c, exp_b, gundlach_b, 
     &gundlach_c, output_coeffs, robin_falloff, sor_maxit
      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 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 CCTKH9
      integer, parameter :: cctki_use_CCTKH9 = kind(CCTKH9)
      integer*8 CCTKH11
      integer, parameter :: cctki_use_CCTKH11 = kind(CCTKH11)
      integer*8 metric_type
      integer, parameter :: cctki_use_metric_type = kind(metric_type)
      integer*8 CCTKH13
      integer, parameter :: cctki_use_CCTKH13 = kind(CCTKH13)
      integer*8 CCTKH14
      integer, parameter :: cctki_use_CCTKH14 = kind(CCTKH14)
      INTEGER*4 CCTKH10
      integer, parameter :: cctki_use_CCTKH10 = kind(CCTKH10)
      INTEGER*4 CCTKH12
      integer, parameter :: cctki_use_CCTKH12 = kind(CCTKH12)
      INTEGER*4 CCTKH15
      integer, parameter :: cctki_use_CCTKH15 = kind(CCTKH15)
      COMMON /ADMBASErest/CCTKH0, CCTKH1, CCTKH2, CCTKH3, initial_data, 
     &CCTKH4, CCTKH5, CCTKH6, CCTKH7, CCTKH8, CCTKH9, CCTKH11, metric_ty
     &pe, CCTKH13, CCTKH14, CCTKH10, CCTKH12, CCTKH15
      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

      integer i,j,k

      REAL*8 x1,y1,z1,rho1,rho2
      REAL*8 phi,phif,e2q
      REAL*8 brillq,eps
      REAL*8 zero,one,two,three

c

      zero  = 0.0D0
      one   = 1.0D0
      two   = 2.0D0
      three = 3.0D0

c

      eps = cctk_delta_space(1)

      do k=1,cctk_lsh(3)
         do j=1,cctk_lsh(2)
            do i=1,cctk_lsh(1)

               x1 = x(i,j,k)
               y1 = y(i,j,k)
               z1 = z(i,j,k)

               rho2 = x1*x1 + y1*y1
               rho1 = sqrt(rho2)
               phi = phif(x1,y1)

               e2q  = exp(two*brillq(rho1,z1,phi))

c

               brillpsi(i,j,k) = one

c
c
c
c

c
c

c
c
c
c

               if (rho1.gt.rhofudge) then

                  gxx(i,j,k) = e2q + (one - e2q)*y1*y1/rho2
                  gyy(i,j,k) = e2q + (one - e2q)*x1*x1/rho2
                  gzz(i,j,k) = e2q
                  gxy(i,j,k) = - (one - e2q)*x1*y1/rho2

                  brillMlinear(i,j,k) = 0.25/e2q
     .                 *(brillq(rho1,z1+eps,phi) 
     .                 + brillq(rho1,z1-eps,phi) 
     .                 + brillq(rho1+eps,z1,phi) 
     .                 + brillq(rho1-eps,z1,phi) 
     .                 - 4.0*brillq(rho1,z1,phi))
     .                 / eps**2

                  brillMlinear(i,j,k) = brillMlinear(i,j,k) + 0.25/rho2
     .               *(three*0.25*(brillq(rho1,z1,phi+eps)
     .               - brillq(rho1,z1,phi-eps))**2
     .               + two*(brillq(rho1,z1,phi+eps)
     .               - two*brillq(rho1,z1,phi)
     .               + brillq(rho1,z1,phi-eps)))/eps**2

               else

                  gxx(i,j,k) = one 
                  gyy(i,j,k) = one 
                  gzz(i,j,k) = one
                  gxy(i,j,k) = zero

                  brillMlinear(i,j,k) = 0.25/e2q
     .                 *(brillq(rho1,z1+eps,phi) 
     .                 + brillq(rho1,z1-eps,phi) 
     .                 + two*brillq(eps,z1,phi) 
     .                 - two*brillq(rho1,z1,phi))
     .                 / eps**2

               end if

               gxz(i,j,k) = zero
               gyz(i,j,k) = zero

c
               brillNsource(i,j,k) = zero

            end do
         end do
      end do


c

c
c

      return
      end
