




















































































































 


































































































































































 






















































































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                


















































        


      subroutine IDLinearWaves_TeukWaves(cctk_dim,cctk_gsh,cctk_lsh,cctk
     &_lbnd,cctk_ubnd,cctk_tile_min,cctk_tile_max,cctk_ash,cctk_alignmen
     &t,cctk_alignment_offset,cctk_from,cctk_to,cctk_bbox,cctk_delta_tim
     &e,cctk_time,cctk_delta_space,cctk_origin_space,cctk_levfac,cctk_le
     &voff,cctk_levoffdenom,cctk_timefac,cctk_convlevel,cctk_convfac,cct
     &k_nghostzones,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,betaz_p,betaz_p_p,coarse_dx,coarse_dy,coarse_dz,conformal_s
     &tate,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_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,p
     &siyy,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 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,m,ra,pi
      REAL*8 wave,wave2,wave3,wave4,wave5,wave6,wave7,wave8
      INTEGER ipacket,iparity
      INTEGER ingoing,outgoing
      REAL*8 kappa
c
      REAL*8 xp,yp,zp,rp
      
c
      REAL*8 teuk_grr,teuk_grt,teuk_grp,teuk_gtt,teuk_gtp,teuk_gpp
      REAL*8 teuk_hrr,teuk_hrt,teuk_hrp,teuk_htt,teuk_htp,teuk_hpp
      
c
c
c
      REAL*8 teuk_tp,teuk_tn
      REAL*8 tp2,tn2,fp,fn,fpa,fna,fpb,fnb,fpc,fnc,fpd,fnd,fpe,fne
      REAL*8 ca,cb,cc,ck,cl,frr,frt,frp,ftt1,ftt2,ftp,fpp1,fpp2
      REAL*8 cadot,cbdot,ccdot,ckdot,cldot
      REAL*8 drt,drp,dtt,dtp,dpp,sint,cost,sinp,cosp,tmp

c
      REAL*8 drx,dry,drz,dtx,dty,dtz,dpx,dpy,dpz
      
      CHARACTER*200 infoline
      INTEGER i,j,k
      INTEGER CCTK_Equals

      REAL*8 amplitude
      integer, parameter :: cctki_use_amplitude = kind(amplitude)
      REAL*8 wavecenter
      integer, parameter :: cctki_use_wavecenter = kind(wavecenter)
      REAL*8 wavelength
      integer, parameter :: cctki_use_wavelength = kind(wavelength)
      REAL*8 wavephi
      integer, parameter :: cctki_use_wavephi = kind(wavephi)
      REAL*8 wavepulse
      integer, parameter :: cctki_use_wavepulse = kind(wavepulse)
      REAL*8 wavetheta
      integer, parameter :: cctki_use_wavetheta = kind(wavetheta)
      integer*8 packet
      integer, parameter :: cctki_use_packet = kind(packet)
      integer*8 parity
      integer, parameter :: cctki_use_parity = kind(parity)
      integer*8 teuk_no_vee
      integer, parameter :: cctki_use_teuk_no_vee = kind(teuk_no_vee)
      integer*8 wavesgoing
      integer, parameter :: cctki_use_wavesgoing = kind(wavesgoing)
      INTEGER*4 mvalue
      integer, parameter :: cctki_use_mvalue = kind(mvalue)
      COMMON /IDLinearWavespriv/amplitude, wavecenter, wavelength, wavep
     &hi, wavepulse, wavetheta, packet, parity, teuk_no_vee, wavesgoing,
     & mvalue
      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
            
      pi    = 3.14159265358979
      kappa = sqrt(143.0d0/pi)/12288.0d0
      
      if (wavecenter.ne.0.) then
         call CCTK_Warn(0,70,"teukwaves.F77","IDLinearWaves","Teukwaves:
     & wavecenter must be zero for time symmetry")
      endif
 
      amp  = amplitude
      m    = mvalue
      wave = wavelength
      ra   = wavecenter
      
c
c
c
      if (CCTK_Equals(packet,'eppley') .eq. 1) then
         ipacket = 1
         call CCTK_Info("IDLinearWaves",'Teukolsky Packet = eppley') 
      elseif (CCTK_Equals(packet,'evans') .eq. 1) then
         ipacket = 2
         call CCTK_Info("IDLinearWaves",'Teukolsky Packet = evans')
c
c
c
         amp = amp*(4.0d0/3.0d0)*kappa
      elseif (CCTK_Equals(packet,'square') .eq. 1) then
         ipacket = 3
         call CCTK_Info("IDLinearWaves",'Teukolsky Packet = square')
      end if
      
      write(infoline,'(A13,G12.7)')
     &    ' amplitude = ',amp
      call CCTK_Info("IDLinearWaves",infoline)
      write(infoline,'(A11,G12.7)')
     &    ' m value = ',m
      call CCTK_Info("IDLinearWaves",infoline)
      write(infoline,'(A14,G12.7)')
     &    ' wavecenter = ',ra
      call CCTK_Info("IDLinearWaves",infoline)
      write(infoline,'(A14,G12.7)')
     &    ' wavelength = ',wave
      call CCTK_Info("IDLinearWaves",infoline)
      
      if (CCTK_Equals(parity,'even') .eq. 1) then    
         iparity = 0
      elseif (CCTK_Equals(parity,'odd') .eq. 1) then 
         iparity = 1
      endif
      
      if(CCTK_Equals(wavesgoing,'in') .eq. 1) then
         ingoing = 1
      endif
      
      if(CCTK_Equals(wavesgoing,'out') .eq. 1) then
         outgoing = 1
      endif
      
      if(CCTK_Equals(wavesgoing,'both') .eq. 1) then
         ingoing  = 1
         outgoing = 1
      endif
      
      if (ipacket.eq.1) then
         if (ingoing.eq.0.or.outgoing.eq.0) then
            call CCTK_Warn(4,130,"teukwaves.F77","IDLinearWaves",'Epply 
     &packet is only non-singular at the origin for ingoing-outgoing com
     &bination of waves')
         endif
      endif
      
      wave2 = wave*wave
      wave3 = wave2*wave
      wave4 = wave3*wave
      wave5 = wave4*wave
      wave6 = wave5*wave
      wave7 = wave6*wave
      wave8 = wave7*wave
         

c
c
c
      do k=1,cctk_lsh(3)
         do j=1,cctk_lsh(2)
            do i=1,cctk_lsh(1)
               
               xp = x(i,j,k)
               yp = y(i,j,k)
               zp = z(i,j,k)
               rp = r(i,j,k)
               
               teuk_tp = (cctk_time+rp-ra)
               teuk_tn = (cctk_time-rp+ra)
               
               tp2 = teuk_tp**2
               tn2 = teuk_tn**2

c
               fp = 0.0d0
               fpa = 0.0d0
               fpb = 0.0d0
               fpc = 0.0d0
               fpd = 0.0d0
               fpe = 0.0d0

               fn = 0.0d0
               fna = 0.0d0
               fnb = 0.0d0
               fnc = 0.0d0
               fnd = 0.0d0
               fne = 0.0d0
               
c
               if (ipacket.eq.1) then
c
                  if (ra.ne.0.and.ingoing.eq.1.and.outgoing.eq.1) then
                     teuk_tn = (cctk_time-rp-ra)
                     tn2 = teuk_tn**2
                  endif
                  
                  if (ingoing.eq.1) then
                     
                     tmp = exp(-tp2)
                     fp = amp * teuk_tp*tmp
                     fpa = amp * (1 - 2*tp2)*tmp
                     fpb = amp * teuk_tp*(4*tp2 - 6)*tmp
                     fpc = amp * (24*tp2 - 8*tp2*tp2 - 6)*tmp
                     fpd = amp * teuk_tp*(60 - 80*tp2 + 16*tp2*tp2)*tmp
                     fpe = amp * (
     &                    (60 - 360*tp2 + 240*tp2*tp2 - 32*tp2*tp2*tp2)*
     &tmp )
                     
                     if (outgoing.eq.0.and.(teuk_tp+wavepulse).le.0) 
     &                    then
                        
                        fp  = 0.0d0
                        fpa = 0.0d0
                        fpb = 0.0d0
                        fpc = 0.0d0
                        fpd = 0.0d0
                        fpe = 0.0d0
                        
                     endif
                     
                  endif
                  
                  if (outgoing.eq.1) then
                     
                     tmp = exp(-tn2)
                     fn = amp * teuk_tn*tmp
                     fna = amp * (1 - 2*tn2)*tmp
                     fnb = amp * teuk_tn*(4*tn2 - 6)*tmp
                     fnc = amp * (24*tn2 - 8*tn2*tn2 -6)*tmp
                     fnd = amp * teuk_tn*(60 - 80*tn2 + 16*tn2*tn2)*tmp
                     fne = amp * (
     &                    (60 - 360*tn2 + 240*tn2*tn2 - 32*tn2*tn2*tn2)*
     &tmp )
                     
                  endif
               endif
               
               
c
               if(ipacket.eq.2) then
                  
                  if (ingoing.eq.1) then
                     
                     tmp = 1.0d0 - tp2/wave2
                     
                     if (abs(teuk_tp).lt.wave) then
                        
                        fp  =  amp*wave5*tmp**6
                        fpa =  -12.0d0*amp*wave3*teuk_tp*tmp**5
                        fpb =  -12.0d0*amp*wave3*(
     &                       tmp**5 - 10.0d0*tp2/wave2*tmp**4 )
                        fpc =  120.0d0*amp*wave3*(
     &                       3.0d0*teuk_tp/wave2*tmp**4 
     &                       -8.0d0*teuk_tp*tp2/wave4*tmp**3 )
                        fpd = 360.0d0*amp*wave3*(
     &                       1.0d0/wave2*tmp**4 
     &                       -16.0d0*tp2/wave4*tmp**3
     &                       +16.0d0*tp2*tp2/wave6*tmp**2 )
                        fpe = 2880.0d0*amp*wave3*(
     &                       -5.0d0*teuk_tp/wave4*tmp**3
     &                       +20.0d0*teuk_tp*tp2/wave6*tmp**2
     &                       -8.0d0*teuk_tp*tp2*tp2/wave8*tmp )
                        
                     else
                        fp=0.
                        fpa=0.
                        fpb=0.
                        fpc=0.
                        fpd=0.
                        fpe=0.
                     endif
                     
                  endif
                  
                  if (outgoing.eq.1) then
                     
                     tmp = 1.0d0 - tn2/wave2
                     
                     if (abs(teuk_tn).lt.wave) then
                        
                        fn  =  amp*wave5*tmp**6
                        fna =  -12.0d0*amp*wave3*teuk_tn*tmp**5
                        fnb =  -12.0d0*amp*wave3*(
     &                       tmp**5 - 10.0d0*tn2/wave2*tmp**4 )
                        fnc =  120.0d0*amp*wave3*(
     &                       3.0d0*teuk_tn/wave2*tmp**4 
     &                       -8.0d0*teuk_tn*tn2/wave4*tmp**3 )
                        fnd = 360.0d0*amp*wave3*(
     &                       1.0d0/wave2*tmp**4 
     &                       -16.0d0*tn2/wave4*tmp**3
     &                       +16.0d0*tn2*tn2/wave6*tmp**2 )
                        fne = 2880.0d0*amp*wave3*(
     &                       -5.0d0*teuk_tn/wave4*tmp**3
     &                       +20.0d0*teuk_tn*tn2/wave6*tmp**2
     &                       -8.0d0*teuk_tn*tn2*tn2/wave8*tmp )
                        
                     else
                        fn=0.
                        fna=0.
                        fnb=0.
                        fnc=0.
                        fnd=0.
                        fne=0.
                     endif
                     
                  endif
                  
               endif
      
c
               if(ipacket.eq.3) then 
                  
                  call CCTK_Warn(4,298,"teukwaves.F77","IDLinearWaves",'
     &Need to calculate fpe,fne')
                  call CCTK_Warn(4,299,"teukwaves.F77","IDLinearWaves",'
     &Need conditionals for in out waves')
                  
                  if (abs(teuk_tp).lt.wave) then
                     fp = amp * (1-(tp2/wave2))**2
                     fpa = amp * (-4*teuk_tp)*(1-tp2/wave2)/wave2
                     fpb = amp * (8*tp2/wave4-4*(1-tp2/wave2)/wave2)
                     fpc = amp * 24*teuk_tp/wave4
                     fpd = amp * 24/wave4
                  else
                     fp=0.
                     fpa=0.
                     fpb=0.
                     fpc=0.
                     fpd=0.
                  endif
         
                  if (abs(teuk_tn).lt.wave) then
                     fn = amp * (1-(tn2/wave2))**2
                     fna = amp * (-4*teuk_tn)*(1-tn2/wave2)/wave2
                     fnb = amp * (8*tn2/wave4-4*(1-tn2/wave2)/wave2)
                     fnc = amp * 24*teuk_tn/wave4
                     fnd = amp * 24/wave4  
                  else
                     fn=0.
                     fna=0.
                     fnb=0.
                     fnc=0.
                     fnd=0.
                  endif
               endif
               
c
c
c
               ca = 3*( (fpb-fnb)/rp**3 -3*(fna+fpa)/rp**4 +3*(fp-fn)/rp
     &**5 )
               cb = -( -(fnc+fpc)/rp**2 +3*(fpb-fnb)/rp**3 -6*(fna+fpa)/
     &rp**4
     &              +6*(fp-fn)/rp**5 ) 
               cc = ( (fpd-fnd)/rp -2*(fnc+fpc)/rp**2 +9*(fpb-fnb)/rp**3
     &
     &              -21*(fna+fpa)/rp**4 +21*(fp-fn)/rp**5 )/4
               ck = (fpb-fnb)/rp**2-3*(fpa+fna)/rp**3+3*(fp-fn)/rp**4
               cl = -(fpc+fnc)/rp+2*(fpb-fnb)/rp**2-3*(fpa+fna)/rp**3+
     &              3*(fp-fn)/rp**4
               
c
c
c
c
c
               
               if(ingoing.eq.0.and.outgoing.eq.1) then
                  ca=-ca
                  cb=-cb
                  cc=-cc
                  ck=-ck
                  cl=-cl
               endif
               
               sint = (xp**2+yp**2)**0.5/rp
               cost = zp/rp
               sinp = yp/(xp**2+yp**2)**0.5
               cosp = xp/(xp**2+yp**2)**0.5
               
c
               if (m.eq.0) then
                  frr = 2-3*sint**2
                  frt = -3*sint*cost
                  frp = 0.
                  ftt1 = 3*sint**2
                  ftt2 = -1.
                  ftp = 0.
                  fpp1 = -ftt1
                  fpp2 = 3*sint**2-1
                  drt = 0.
                  drp = -4*cost*sint
                  dtt = 0.
                  dtp = -sint**2
                  dpp = 0.
               elseif (m.eq.-1) then
                  frr = 2*sint*cost*sinp
                  frt = (cost**2-sint**2)*sinp 
                  frp = cost*cosp 
                  ftt1 = -frr
                  ftt2 = 0.
                  ftp = sint*cosp 
                  fpp1 = -ftt1
                  fpp2 = ftt1
                  drt = -2*cost*sinp
                  drp = -2*(cost**2-sint**2)*cosp
                  dtt = -sint*sinp
                  dtp = -cost*sint*cosp
                  dpp = sint*sinp
               elseif (m.eq.-2) then
                  frr = sint**2*2*sinp*cosp
                  frt = sint*cost*2*sinp*cosp
                  frp = sint*(1-2*sinp**2)
                  ftt1 = (1+cost**2)*2*sinp*cosp
                  ftt2 = -2*sinp*cosp
                  ftp = cost*(1-2*cosp**2)
                  fpp1 = -ftt1
                  fpp2 = cost**2*2*sinp*cosp
                  drt = 8*sint*sinp*cosp
                  drp = 4*sint*cost*(1-2*sinp**2)
                  dtt = -4*cost*sinp*cosp
                  dtp = (2-sint**2)*(2*sinp**2-1)
                  dpp = 2*cost*2*sinp*cosp
               elseif (m.eq.1) then
                  frr = 2*sint*cost*cosp
                  frt = (cost**2-sint**2)*cosp
                  frp = -cost*sinp
                  ftt1 = -2*sint*cost*cosp
                  ftt2 = 0.
                  ftp = -sint*sinp
                  fpp1 = -ftt1
                  fpp2 = ftt1
                  drt = -2*cost*cosp
                  drp = 2*(cost**2-sint**2)*sinp
                  dtt = -sint*cosp
                  dtp = cost*sint*sinp
                  dpp = sint*cosp
               elseif (m.eq.2) then
                  frr = sint**2*(1-2*sinp**2)
                  frt = sint*cost*(1-2*sinp**2)
                  frp = -sint*2*sinp*cosp
                  ftt1 = (1+cost**2)*(1-2*sinp**2)
                  ftt2 = (2*sinp**2-1)
                  ftp = cost*2*sinp*cosp
                  fpp1 = -ftt1
                  fpp2 = cost**2*(1-2*sinp**2)
                  drt = 4*sint*(1-2*sinp**2)
                  drp = -4*sint*cost*2*sinp*cosp
                  dtt = -2*cost*(1-2*sinp**2)
                  dtp = (2-sint**2)*2*sinp*cosp
                  dpp = 2*cost*(1-2*sinp**2)
               else
                  call CCTK_Warn(0,433,"teukwaves.F77","IDLinearWaves",'
     &m should be one of -2,-1,0,1,2')
               endif
               
c
               if (iparity.eq.0) then    
                  teuk_grr = 1 + ca*frr
                  teuk_grt = cb*frt*rp
                  teuk_grp = cb*frp*(xp**2+yp**2)**0.5
                  teuk_gtt = rp**2*(1 + cc*ftt1 + ca*ftt2)
                  teuk_gtp = (ca - 2*cc)*ftp*rp*(xp**2+yp**2)**0.5
                  teuk_gpp = (xp**2+yp**2)*(1 + cc*fpp1 + ca*fpp2)
               else
                  teuk_grr = 1.
                  teuk_grt = ck*drt*rp
                  teuk_grp = ck*drp*rp*sint
                  teuk_gtt = (1+cl*dtt)*rp**2
                  teuk_gtp = cl*dtp*rp**2*sint
                  teuk_gpp = (1+cl*dpp)*rp**2*sint**2
               endif
               
c
               
c
               drx = xp/rp
               dry = yp/rp
               drz = zp/rp
      
               dtx = xp*zp/(rp**2*sqrt(xp**2+yp**2))
               dty = yp*zp/(rp**2*sqrt(xp**2+yp**2))
               dtz = (zp**2/rp**2-1)/sqrt(xp**2+yp**2)
               
               dpx = -yp/(xp**2+yp**2)
               dpy = xp/(xp**2+yp**2)
               dpz = 0.
               
c
               gxx(i,j,k) = drx**2*teuk_grr+2.*drx*dtx*teuk_grt
     &              +2.*drx*dpx*teuk_grp+dtx**2*teuk_gtt
     &              +2.*dtx*dpx*teuk_gtp+dpx**2*teuk_gpp
               
               gyy(i,j,k) = dry**2*teuk_grr+2.*dry*dty*teuk_grt
     &              +2.*dry*dpy*teuk_grp+dty**2*teuk_gtt
     &              +2.*dty*dpy*teuk_gtp+dpy**2*teuk_gpp
 
               gzz(i,j,k) = drz**2*teuk_grr+2.*drz*dtz*teuk_grt
     &              +2.*drz*dpz*teuk_grp+dtz**2*teuk_gtt
     &              +2.*dtz*dpz*teuk_gtp+dpz**2*teuk_gpp
               
               gxy(i,j,k) = drx*dry*teuk_grr+(drx*dty+dtx*dry)*teuk_grt
     &              +(drx*dpy+dpx*dry)*teuk_grp
     &              +dtx*dty*teuk_gtt+(dtx*dpy+dpx*dty)*teuk_gtp+
     &              dpx*dpy*teuk_gpp
               
               gxz(i,j,k) = drx*drz*teuk_grr+(drx*dtz+dtx*drz)*teuk_grt
     &              +(drx*dpz+dpx*drz)*teuk_grp
     &              +dtx*dtz*teuk_gtt+(dtx*dpz+dpx*dtz)*teuk_gtp+
     &              dpx*dpz*teuk_gpp
               
               gyz(i,j,k) = dry*drz*teuk_grr+(dry*dtz+dty*drz)*teuk_grt
     &              +(dry*dpz+dpy*drz)*teuk_grp
     &              +dty*dtz*teuk_gtt+(dty*dpz+dpy*dtz)*teuk_gtp+
     &              dpy*dpz*teuk_gpp
                


c
               cadot = 3*( (fpc-fnc)/rp**3 -3*(fnb+fpb)/rp**4 
     &              +3*(fpa-fna)/rp**5 )
               cbdot = -( -(fnd+fpd)/rp**2 +3*(fpc-fnc)/rp**3 
     &              -6*(fnb+fpb)/rp**4
     &              +6*(fpa-fna)/rp**5 ) 
               ccdot = ( (fpe-fne)/rp -2*(fnd+fpd)/rp**2 +9*(fpc-fnc)/rp
     &**3
     &              -21*(fnb+fpb)/rp**4 +21*(fpa-fna)/rp**5 )/4
               ckdot = (fpc-fnc)/rp**2-3*(fpb+fnb)/rp**3+3*(fpa-fna)/rp*
     &*4
               cldot = -(fpd+fnd)/rp+2*(fpc-fnc)/rp**2-3*(fpb+fnb)/rp**3
     &+
     &              3*(fpa-fna)/rp**4
               
c
               if (iparity.eq.0) then    
                  teuk_hrr = -0.5*(cadot*frr)
                  teuk_hrt = -0.5*(cbdot*frt*rp)
                  teuk_hrp = -0.5*(cbdot*frp*(xp**2+yp**2)**0.5)
                  teuk_htt = -0.5*(rp**2*(ccdot*ftt1 + cadot*ftt2))
                  teuk_htp = 
     &                 -0.5*((cadot - 2*ccdot)*ftp*rp*(xp**2+yp**2)**0.5
     &)
                  teuk_hpp = -0.5*((xp**2+yp**2)*(ccdot*fpp1 + cadot*fpp
     &2))
               else
                  teuk_hrr = 0.
                  teuk_hrt = -0.5*(ckdot*drt*rp)
                  teuk_hrp = -0.5*(ckdot*drp*rp*sint)
                  teuk_htt = -0.5*((cldot*dtt)*rp**2)
                  teuk_htp = -0.5*(cldot*dtp*rp**2*sint)
                  teuk_hpp = -0.5*((cldot*dpp)*rp**2*sint**2)
               endif
               
c
               kxx(i,j,k) = drx**2*teuk_hrr+2.*drx*dtx*teuk_hrt
     &              +2.*drx*dpx*teuk_hrp+dtx**2*teuk_htt
     &              +2.*dtx*dpx*teuk_htp+dpx**2*teuk_hpp
               
               kyy(i,j,k) = dry**2*teuk_hrr+2.*dry*dty*teuk_hrt
     &              +2.*dry*dpy*teuk_hrp+dty**2*teuk_htt
     &              +2.*dty*dpy*teuk_htp+dpy**2*teuk_hpp
               
               kzz(i,j,k) = drz**2*teuk_hrr+2.*drz*dtz*teuk_hrt
     &              +2.*drz*dpz*teuk_hrp+dtz**2*teuk_htt
     &              +2.*dtz*dpz*teuk_htp+dpz**2*teuk_hpp
               
               kxy(i,j,k) = drx*dry*teuk_hrr+(drx*dty+dtx*dry)*teuk_hrt
     &              +(drx*dpy+dpx*dry)*teuk_hrp
     &              +dtx*dty*teuk_htt+(dtx*dpy+dpx*dty)*teuk_htp+
     &              dpx*dpy*teuk_hpp
               
               kxz(i,j,k) = drx*drz*teuk_hrr+(drx*dtz+dtx*drz)*teuk_hrt
     &              +(drx*dpz+dpx*drz)*teuk_hrp
     &              +dtx*dtz*teuk_htt+(dtx*dpz+dpx*dtz)*teuk_htp+
     &              dpx*dpz*teuk_hpp
               
               kyz(i,j,k) = dry*drz*teuk_hrr+(dry*dtz+dty*drz)*teuk_hrt
     &              +(dry*dpz+dpy*drz)*teuk_hrp
     &              +dty*dtz*teuk_htt+(dty*dpz+dpy*dtz)*teuk_htp+
     &              dpy*dpz*teuk_hpp
               
               
               
            enddo
         enddo
      enddo

c
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

         else if((CCTK_Equals(conformal_storage,"factor+derivs+2nd deriv
     &s") .ne. 0)) then
            conformal_state = 3
         end if

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

                  psi(i,j,k)   = 1d0

                  if(conformal_state .gt. 1) then
                     psix(i,j,k)  = 0d0
                     psiy(i,j,k)  = 0d0
                     psiz(i,j,k)  = 0d0    
                  endif
                  if(conformal_state .gt. 2) then
                     psixy(i,j,k) = 0d0
                     psixz(i,j,k) = 0d0
                     psiyz(i,j,k) = 0d0  
                     psixx(i,j,k) = 0d0
                     psiyy(i,j,k) = 0d0
                     psizz(i,j,k) = 0d0
                  endif
               end do
            end do
         end do
      end if
      
      return
      end 
      

      
