!/*@@
!
! $Id: input.f90,v 1.18 1997/09/12 23:28:20 mijan Exp $
!
!@@*/

!/*@@
!   @routine read_parameters
!   @author Tomasz Haupt
!   @desc
!   reads parameters from a file
!   @enddescr
!   @var PIN
!   @vdesc logical number of the read unit (use 10 for file fort.10)
!   @vtype integer
!   @vio IN
!   @vcomment temporary quick solution; will be improved
!   @endvar 
!   @var params
!   @vdesc list of parameters to be read in
!   @vtype integer, dimension(100)
!   @vio OUT
!   @endvar 
!   @calledby main
!@@*/

subroutine  read_parameters(PIN, params)
use input_params
implicit none
integer PIN,i,inum
real*8, dimension(100) :: params

integer nparam
real*8,dimension(3) :: h

nparam=96

! ---------------------------
! Base Grid Parameters
! ---------------------------
! 1. nx
! 2. ny 
! 3. nz
! 4. dxmin
! 5. dxmax
! 6. dymin
! 7. dymax
! 8. dzmin
! 9. dzmax
! ---------------------------
! Evolution control
! ---------------------------
!10. CourantFactor
!11. BG_Steps
!12. StartTime
!13. Tolerance (main evolution loop)
!14. ExtBndryWidth
!15. MaxLev
!16. RefineFactor
!17. TruncErrTrsh  (threshold value of the truncatrion error)
!18. Parameter to be determined
!19. Outer Boundary Choice 
!         Dirichlet Condition     = 1.0 (Dependent on Initial data)
!         Sommerfeld Condition    = 2.0 (Dependent on params 66 and 67)
!         Periodic Condition      = 3.0 (Dependent on params 58)
!         Perturbative OB         = 4.0
!         Characteristic OB       = 5.0
!         Freezing Condition      = 6.0
!         Pert. Somm. OB          = 7.0 (Dependent on params 66 and 67)
!20. Inner Boundary Choice
!         No boundary: One-sided diff  = 0.0
!         Freezing in buffer zone      = 1.0 See params 48 and 68
!-------------------------------
! initial values
!-------------------------------
!21. type of initial conditions
!    1=iso
!    2=read in initial data from hdf files
!    3=linear waves
!    4=boosted ef
!    5=kerr
!    6=Harmonic Schwarzschild
!22. amplitude of linear wave
!23. width of linear wave
!24. initial radius of hole 1
!25. initial radius of hole 2
!26. initial x coord of hole 1 center
!27. initial y coord of hole 1 center
!28. initial z coord of hole 1 center
!29. initial x coord of hole 2 center
!30. initial y coord of hole 2 center
!31. initial z coord of hole 2 center
!32. G-code slicing method (Only if param34=1)
!33. G-code lapse exponent (Only if param34=1)
!34. lapse choice
!          1.0 Glapse,
!          2.0 Isotropic,
!          4.0 Boosted Eddington-Finkelstein,
!          5.0 Kerr,
!          6.0 Algebraic Kerr-Schild form 1
!          7.0 Algebraic Kerr-Schild form 2
!          8.0 Harmonic Schwarzschild Algebraic
!35. shift choice
!          0.0 zero shift,
!          1.0 periodic(space and time) shift see params(60,61,62,63,64,65,66)
!          2.0 Boosted Eddington-Finkelstein,
!          3.0 Kerr,
!          4.0 Areal
!          5.0 Harmonic Schwarzschild Algebraic
!          6.0 Constant in z-direction (see param 59)
!36.  initialization convergence criterion
!37.  apparent horizon  variable ilucgstp (unused)
!38.  apparent horizon stopping criteria  (unused)
!39.  apparent horizon variable epsilon (unused)
!40.  apparent horizon variable nst (unused)
!41.  apparent horizon variable nsp (unused)
!42.  turn on/off the apparent horizon locator (unused)
!43.  velocity of boost in z direction (Param21=4)
!44.  (unused)
!45.  Frequency of 3D hdf file dump (# of timesteps)
!46.  Evolution scheme, 0=leap, 1=crank 
!47.  (unused)
!48.  number of buffer zones inside horizon (approx.)
!49.  0=use inner-bndry pts for interps, 1=don't use them
!50.  width of inner-boundary smoothing domain
!51.  number of inner-boundary smoothing passes
!52.  Spin parameter in Kerr data [0-1) 
!53.  order of boundary derivative scheme (1 or 2)
!54.  order of interpolation (4th or 1st)
!55.  dissipation parameter epsilon [0-1]
!56.  Max number of CN iterations (only if p46=1)
!57.  Convergence criterion for CN scheme (only if p46=1)
!58.  Periodic boundary conditions stencil size (0=no periodic bc)
!59.  Value for shift in the z-direction (only if p35=6)
!60.  Value for beta0 in the x-direction (only if p35=1)
!61.  Value for beta0 in the y-direction (only if p35=1)
!62.  Value for beta0 in the z-direction (only if p35=1)
!63.  Frequency in time for periodic shift(only if p35=1)
!64.  Spatial frequency for periodic shift (L/2pi) (only if p35=1)
!65.  Phase of periodic shift in the z-direction (/pi) (only if p35=1)
!66.  Simple Sommerfeld toggle(1:=K_ij only, 2:=g_ij only, 3:=both)
!67.  Simple Sommerfeld falloff ( assumes nonzero p66, works only if p19=2,7)
!         for p66=1 or 2 power of r taken as is. ie: u = r^p g_ij or K_ij
!         for p66=3 gives r^p g_ij and r^(p+1) K_ij
!68.  Outermost radius for inner buffer freezing(in # of points)
!69.  Inner boundary Mask Geometry (0=sphere,1=cube)
!70.  Type of difference operator to use (0=reg, 1=back in y)
!71.  (unused)
!72.  (unused)
!73.  (unused)
!74.  (unused)
!75.  (unused)
!76.  (unused)
!77.  (unused)
!78.  (unused)
!79.  (unused)
!80.  (unused)
!81.  (unused)
!82.  (unused)
!83.  (unused)
!84.  (unused)
!85.  (unused)
!86.  (unused)
!87.  (unused)
!88.  (unused)
!89.  (unused)
!90.  (unused)
!91.  (unused)
!92.  (unused)
!93.  (unused)
!94.  (unused)
!95.  (unused)
!96.  (unused)
!-------------------------------
! Calculated at execution time
!-------------------------------
!97.  Time
!98. h (at current level) 
!99. dt (at current level)   

do i=1,nparam
   read(PIN,*) params(i), inum
   if(i /= inum) then
      write(*,*) "Input parameter ",i," is incorrect!"
      stop
   end if
enddo

!-------------------------------------------\
! For periodic b.c.'s, must add points       \
!  and redefine min/max values to ensure that \
!  periodic boundaries are where the user     /
!  thinks they are.                          /
!-------------------------------------------/

if(params(P_OUTER_BOUNDARY).eq. 3.0d0) then

   !---------------------\
   ! Save grid spacing    >
   !---------------------/
   
   h(1) = (params(P_XMAX)-params(P_XMIN))/(int(params(P_NX))-1)
   h(2) = (params(P_YMAX)-params(P_YMIN))/(int(params(P_NY))-1)
   h(3) = (params(P_ZMAX)-params(P_ZMIN))/(int(params(P_NZ))-1)
   
   !---------------------------------------------\
   ! Grow grid by 2s-1 zones in each direction    >
   !---------------------------------------------/
   
   params(P_NX) = params(P_NX) + 2*params(P_OB_PERIODIC)-1.d0
   params(P_NY) = params(P_NY) + 2*params(P_OB_PERIODIC)-1.d0
   params(P_NZ) = params(P_NZ) + 2*params(P_OB_PERIODIC)-1.d0
   
   !------------------------------------------------\
   ! Rescale min to lie on its former grid point - s >
   !------------------------------------------------/
   
   params(P_XMIN) = params(P_XMIN) - int(params(P_OB_PERIODIC))*h(1)
   params(P_YMIN) = params(P_YMIN) - int(params(P_OB_PERIODIC))*h(2)
   params(P_ZMIN) = params(P_ZMIN) - int(params(P_OB_PERIODIC))*h(3)
   
   !--------------------------------------------------\
   ! Rescale max to lie on its former grid point + s-1 >
   !--------------------------------------------------/
   
   params(P_XMAX) = params(P_XMAX) + (int(params(P_OB_PERIODIC))-1)*h(1)
   params(P_YMAX) = params(P_YMAX) + (int(params(P_OB_PERIODIC))-1)*h(2)
   params(P_ZMAX) = params(P_ZMAX) + (int(params(P_OB_PERIODIC))-1)*h(3)
   
endif

!-----------------------------------------------------------------------
! If we are using the perturbative outer boundaries then read in 
! necessary parameters.
!-----------------------------------------------------------------------
if ( (params(P_OUTER_BOUNDARY).eq.4).or.   &
     (params(P_OUTER_BOUNDARY).eq.7) ) then
   call init_radx(params)
end if


return
end
