
SUBROUTINE init_radx(params)

  INTERFACE

     SUBROUTINE input_radx(params)
       USE input_worldtube 
       USE input_exterior
       USE input_params
       real*8, dimension(100), intent(in) :: params
     END SUBROUTINE input_radx

     SUBROUTINE worldtube_init
       USE input_worldtube
       IMPLICIT NONE
     END SUBROUTINE worldtube_init

     SUBROUTINE exterior_init(params)
       USE sch_mpe
       IMPLICIT NONE
       real*8, dimension(100), intent(in) :: params
     END SUBROUTINE exterior_init

  END INTERFACE

  real*8, dimension(100) :: params

  CALL input_radx(params)
  CALL worldtube_init
  CALL exterior_init(params)

END SUBROUTINE init_radx

 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



!!$ RadX input subroutine
!!$ Version 1.6
!!$
!!$ Reads in values for input_worldtube variables from the file wt.dat

SUBROUTINE input_radx(params)

  USE input_worldtube 
  USE input_exterior
  USE input_params

  implicit none
  real*8, dimension(100),intent(in) :: params
  real*8, dimension(3)              :: dx

  OPEN (unit=1,file='radx.in')

  READ (1,501) nn
  READ (1,500) wt_r
  READ (1,500) router
  READ (1,500) cfext
  READ (1,501) df
  READ (1,501) lmin
  READ (1,501) lmax
  read (1,501) mmax
  read (1,501) interp_order
  read (1,501) evol_type

500 FORMAT(21x, f10.7)
501 FORMAT(21x, i7)

  CLOSE(1)

  dx(1) = (params(P_XMAX)-params(P_XMIN))/(params(P_NX)-1)
  dx(2) = (params(P_YMAX)-params(P_YMIN))/(params(P_NY)-1)
  dx(3) = (params(P_ZMAX)-params(P_ZMIN))/(params(P_NZ)-1)

  dr = params(P_COURANT)*minval(dx)/cfext

  rinner = wt_r
  nrr = int((router-rinner)/dr)

  t_lev = 0

  PRINT *,'nn =          ',nn
  PRINT *,'wt_r =        ',wt_r
  PRINT *,'cfext =       ',cfext
  PRINT *,'router =      ',router
  PRINT *,'nrr =         ',nrr
  PRINT *,'df =          ',df
  PRINT *,'dr =          ',dr

  RETURN

END SUBROUTINE input_radx



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!!$ Allocates memory for worldtube variables, and calls subroutines to
!!$    define the worldtube grid and evaluate the Cartesian -> spherical
!!$    Jacobian

SUBROUTINE worldtube_init

  USE input_worldtube

  IMPLICIT NONE

  INTERFACE

     SUBROUTINE wt_stereo_grid
       USE input_worldtube
       IMPLICIT NONE 
     END SUBROUTINE wt_stereo_grid

     SUBROUTINE wt_define
       USE input_worldtube
       IMPLICIT NONE
     END SUBROUTINE wt_define

     SUBROUTINE wt_jacob
       USE input_worldtube
       IMPLICIT NONE 
     END SUBROUTINE wt_jacob

  END INTERFACE

  REAL*8, DIMENSION(2,3,3,nn,nn)            :: wt_K, wt_L
  REAL*8                                    :: time

  ALLOCATE(q(nn,nn))
  ALLOCATE(p(nn,nn))

  ALLOCATE(wtan(nn,nn,4,3))
  ALLOCATE(wtas(nn,nn,4,3))

  ALLOCATE(jac(2,3,3,nn,nn))

  CALL wt_stereo_grid
  CALL wt_define 
  CALL wt_jacob

END SUBROUTINE worldtube_init



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



!!$ Defines the stereographic coordinate patches on the worldtube,
!!$    where z = q + i p is the complex expression for points in the
!!$    stereographic coordinate grid.
!!$
!!$ This subroutine was part of the Pitt/PSU extraction module, and
!!$    has been slightly modified


SUBROUTINE wt_stereo_grid

  USE input_worldtube

  IMPLICIT NONE 

  INTEGER                                 :: i

  dd = 2.0/(nn-5)

  DO i=1,nn 
     q(i,:) = -1.0 + (i-3) * dd 
     p(:,i) = -1.0 + (i-3) * dd 
  ENDDO

END SUBROUTINE wt_stereo_grid



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



!!$ Defines the cartesian coordinate values, relative to the interior grid
!!$    of stereographic coordinate points on the world tube.
!!$
!!$ This subroutine was originally part of the Pitt/PSU extraction module
!!$    and has been slightly modified.
!!$
!!$ Modified 12-19-96

SUBROUTINE wt_define

  USE input_worldtube

  IMPLICIT NONE

!!$ Northern Grid

  wtan(:,:,1,1) = 2.0*wt_r*q/(1.0+q**2+p**2)
  wtan(:,:,2,1) = 2.0*wt_r*p/(1.0+q**2+p**2)
  wtan(:,:,3,1) = wt_r*(1.0-q**2-p**2)/(1.0+q**2+p**2)

!!$ Southern Grid

  wtas(:,:,1,1) = 2.0*wt_r*q/(1.0+q**2+p**2)
  wtas(:,:,2,1) = -2.0*wt_r*p/(1.0+q**2+p**2)
  wtas(:,:,3,1) = -wt_r*(1.0-q**2-p**2)/(1.0+q**2+p**2)

END SUBROUTINE wt_define



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



!!$ Calculates the jacobian for the transformation of covariant tensor
!!$    components from a cartesian basis to a spherical basis on the 
!!$    worldtube.


SUBROUTINE wt_jacob

  USE input_worldtube

  IMPLICIT NONE 

  REAL*8,  DIMENSION(2,nn,nn)            :: x,y,z,cosp,sinp
  INTEGER                                :: i,j,axis

  axis = (nn+1)/2

  x(1,:,:) = wtan(:,:,1,1) 
  y(1,:,:) = wtan(:,:,2,1) 
  z(1,:,:) = wtan(:,:,3,1)

  x(2,:,:) = wtas(:,:,1,1) 
  y(2,:,:) = wtas(:,:,2,1) 
  z(2,:,:) = wtas(:,:,3,1)

  DO i=1,nn
     DO j=1,nn
        IF(i.EQ.axis.AND.j.EQ.axis) THEN
           cosp(:,i,j) = 1.
           sinp(:,i,j) = 0.
        ELSE
           cosp(:,i,j) = x(:,i,j)/SQRT( x(:,i,j)**2 + y(:,i,j)**2 )
           sinp(:,i,j) = y(:,i,j)/SQRT( x(:,i,j)**2 + y(:,i,j)**2 )
        ENDIF
     ENDDO
  ENDDO

  jac(:,1,1,:,:) = x(:,:,:)/wt_r
  jac(:,2,1,:,:) = y(:,:,:)/wt_r
  jac(:,3,1,:,:) = z(:,:,:)/wt_r
  jac(:,1,2,:,:) = z(:,:,:)*cosp(:,:,:)
  jac(:,2,2,:,:) = z(:,:,:)*sinp(:,:,:)
  jac(:,3,2,:,:) = -SQRT( x(:,:,:)**2 + y(:,:,:)**2 )
  jac(:,1,3,:,:) = -y(:,:,:)
  jac(:,2,3,:,:) = x(:,:,:)
  jac(:,3,3,:,:) = 0.d0

END SUBROUTINE wt_jacob



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



!!$ Calls subroutines which set up the exterior data structure and
!!$    initializes their values

SUBROUTINE exterior_init(params)

  USE sch_mpe

  IMPLICIT NONE

  INTERFACE

     SUBROUTINE setup_exterior(sch)
       USE sch_mpe
       USE alias
       USE input_exterior
       IMPLICIT NONE
       TYPE(GF3TL)                           :: sch
     END SUBROUTINE setup_exterior

     SUBROUTINE initial_exterior(sch,params)
       USE sch_mpe
       USE alias 
       USE input_exterior
       IMPLICIT NONE
       TYPE(GF3TL)                      :: sch
       real*8, dimension(100)           :: params
     END SUBROUTINE initial_exterior

  END INTERFACE

  real*8, dimension(100) :: params
  TYPE (GF3TL)  :: sch

  CALL setup_exterior(sch)

  CALL initial_exterior(sch,params)

END SUBROUTINE exterior_init



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



!!$ Exterior setup subroutine
!!$ Version 1.4
!!$
!!$ Allocates memory for exterior grid functions. These will be 
!!$    aliased by the module alias.

SUBROUTINE setup_exterior(sch)

  USE sch_mpe
  USE alias
  USE input_exterior

  IMPLICIT NONE

!!$ type(GF3TL) is defined in sch_mpe.f90

  TYPE(GF3TL)                           :: sch
  INTEGER                               :: i,nl,ll,numfcn

! Number of 3-time-level grid functions:

  numfcn = 9

  sch%n_gfcn  = numfcn
  ALLOCATE(sch%gfcn(sch%n_gfcn))

  DO i=1,numfcn
     sch%gfcn(i)%n_t = 3
     ALLOCATE(sch%gfcn(i)%tlev(sch%gfcn(i)%n_t))
  ENDDO

!!$  set up coordinate system

  ALLOCATE(sch%gfcn(1)%tlev(1)%c(1))

  DO i=2,numfcn
     sch%gfcn(i)%tlev(1)%c => sch%gfcn(1)%tlev(1)%c
  ENDDO

  sch%gfcn(1)%tlev(1)%c(1)%n_coord=nrr+1

  ALLOCATE(sch%gfcn(1)%tlev(1)%c(1)%coord(nrr+1))

  sch%gfcn(1)%tlev(1)%c(1)%coord(1) = rinner

  DO i=2,nrr+1
     sch%gfcn(1)%tlev(1)%c(1)%coord(i) = rinner + dr*(i-1)
  ENDDO

!!$ allocate memory for the odd-parity gridfunctions

  ALLOCATE(sch%gfcn(1)%tlev(1)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(1)%tlev(2)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(1)%tlev(3)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(2)%tlev(1)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(2)%tlev(2)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(2)%tlev(3)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(3)%tlev(1)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(3)%tlev(2)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(3)%tlev(3)%data(lmin:lmax,-mmax:mmax,nrr+1))

!!$ allocate memory for the even-parity gridfunctions

  ALLOCATE(sch%gfcn(4)%tlev(1)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(4)%tlev(2)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(4)%tlev(3)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(5)%tlev(1)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(5)%tlev(2)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(5)%tlev(3)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(6)%tlev(1)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(6)%tlev(2)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(6)%tlev(3)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(7)%tlev(1)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(7)%tlev(2)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(7)%tlev(3)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(8)%tlev(1)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(8)%tlev(2)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(8)%tlev(3)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(9)%tlev(1)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(9)%tlev(2)%data(lmin:lmax,-mmax:mmax,nrr+1))
  ALLOCATE(sch%gfcn(9)%tlev(3)%data(lmin:lmax,-mmax:mmax,nrr+1))

  RETURN
END SUBROUTINE setup_exterior



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



SUBROUTINE initial_exterior(sch,params)

!!$ this is initial data for time-symmetric data so all 
!!$ evolve variables are set to zero
!!$ This will, of course, be modified if different initial 
!!$ data is used


!!$ use the standard modules 

  USE sch_mpe
  USE alias 
  USE input_exterior

  IMPLICIT NONE

!!$ sch is the module containing the exterior variables 
!!$ GF3TL is the standard GCBBH data structure

  TYPE(GF3TL)                        :: sch
  real*8, dimension(100), intent(in) :: params
  REAL*8                             :: dt
  INTEGER                            :: i,n_par,l,m
  REAL*8                             :: oneh,zero,ae20
  EXTERNAL ae20

  oneh = 0.5d0
  dt = cfext*dr
  zero = 0.d0

!!$ setup the aliases for the code variables otherwise 
!!$ the code variables look like: 
!!$ sch%gfcn(5)%tlev(1)%data 

  CALL setup_alias(sch)

! Initialize all variables to zero

  ao_m1 = zero
  ao_   = zero
  ao_p1 = zero

  aot_m1 = zero
  aot_   = zero
  aot_p1 = zero

  aor_m1 = zero
  aor_   = zero
  aor_p1 = zero

  ae_m1 = zero
  ae_   = zero
  ae_p1 = zero

  aet_m1 = zero
  aet_   = zero
  aet_p1 = zero

  aer_m1 = zero
  aer_   = zero
  aer_p1 = zero

  he_m1 = zero
  he_   = zero
  he_p1 = zero

  het_m1 = zero
  het_   = zero
  het_p1 = zero

  her_m1 = zero
  her_   = zero
  her_p1 = zero


! Initial data is put into time level n+1 (initial time) and n
! these are shifted down before extraction/1D evolution

  DO i=1,nr
     ae_p1(2,0,i)  = ae20(zero,rr(i),params)
     aet_p1(2,0,i) = (ae20(dt,rr(i),params)-ae20(-dt,rr(i),params))/(2*dt)
     aer_p1(2,0,i) = (ae20(zero,rr(i)+dr,params)-ae20(zero,rr(i)-dr,params))/&
                  (2*dr)
  ENDDO

!!$! Initialize variables at t = -dt for leapfrog integration
!!$! aeh, etc contain ae at t_old
!!$
!!$  do i=2,nr-1
!!$     ae_(2,0,i)  = ae20(-dt,rr(i),params)
!!$     aet_(2,0,i) = (ae20(zero,rr(i),params)-ae20(-2*dt,rr(i),params))/(2*dt)
!!$     aer_(2,0,i) = (ae20(-dt,rr(i)+dr,params)-ae20(-dt,rr(i)-dr,params))/&
!!$                   (2*dr)
!!$  end do

  RETURN 

END SUBROUTINE initial_exterior
