!/*@@
!
! $Id: stub.f90,v 1.23 1997/05/15 16:54:49 mijan Exp haupt $
!
!@@*/

!/*@@
!   @routine main
!   @author T.Haupt
!   @calls BO_Integrate
!   @calls initialize
!   @calls read_parameters
!@@*/

program main
  use input_params
  use parameters
  
  ! This module provides I/O that is easily customizable through
  ! an input file
  
  implicit none

  real*8, dimension(100) :: params
  integer PIN
  integer nx,ny,nz,nst,nsp,stopnow

  integer, dimension(3) :: sh
  real*8, dimension(6)  :: b_box
  real*8 h,h1,h2,h3,dt
  real*8,dimension(:,:,:),allocatable :: s,sc


  integer Level
  integer t 
  integer iter     

  integer gft_close_all

  INTERFACE 

! ******************************************************* write interfaces.

     
     subroutine initialize(params, anm, an, anp, &
	bnm, bn, bnp, cnm, cn, cnp, dnm, dn, dnp, phi, temp)
       use input_params
       real*8, dimension(100)     :: params
       real*8, dimension(:,:,:)   :: anm, an, anp, bnm, bn, bnp, &
          cnm, cn, cnp, dnm, dn, dnp, phi, temp
       
     end subroutine initialize

     subroutine tstep(params, anm, an, anp, bnm, bn, bnp, &
          cnm, cn, cnp, dnm, dn, dnp, phim, phi, phip, temp)
       use input_params
       real*8, dimension(100)     :: params
         real*8, dimension(:,:,:)   :: anm, an, anp, bnm, bn, bnp, &
            cnm, cn, cnp, dnm, dn, dnp, phim, phi, phip, temp
       
     end subroutine tstep


  END INTERFACE



! -----------------------------------------------------------------------
! Sect. 1.0:  Initialize MPI 
! -----------------------------------------------------------------------


! Sect. 2.0 Set parameters
! -----------------------------------------------------------------------

  PIN = 10
  open(PIN,file='adm.in')
  open(33,file='adm.out')
  call read_parameters(PIN, params)

! Sect. 3.0  Set DAGH structures
! -----------------------------------------------------------------------

!  3.1 Base grid geometry
 

  sh(1) = int(params(P_NX))
  sh(2) = int(params(P_NY))
  sh(3) = int(params(P_NZ))
  b_box(1) = params(P_XMIN)
  b_box(2) = params(P_XMAX)
  b_box(3) = params(P_YMIN)
  b_box(4) = params(P_YMAX)
  b_box(5) = params(P_ZMIN)
  b_box(6) = params(P_ZMAX)
  write(6,32) b_box(1),b_box(2),b_box(3),b_box(4),b_box(5),b_box(6)
  write(33,32) b_box(1),b_box(2),b_box(3),b_box(4),b_box(5),b_box(6)
32 format ('Bb',1P,6e10.2)

  params(P_LOCAL_LOWER_X)=0
  params(P_LOCAL_LOWER_Y)=0
  params(P_LOCAL_LOWER_Z)=0
  params(P_LOCAL_UPPER_X)=params(P_NX)-1
  params(P_LOCAL_UPPER_Y)=params(P_NY)-1
  params(P_LOCAL_UPPER_Z)=params(P_NZ)-1


!  3.2 Spatial and temporal resolution of the base grid

  h1 = (b_box(2) - b_box(1))/(sh(1)-1)
  h2 = (b_box(4) - b_box(3))/(sh(2)-1)
  h3 = (b_box(6) - b_box(5))/(sh(3)-1)

  write(6,*) 'shape: ',sh(1),sh(2),sh(3)
  write(33,*)'shape: ',sh(1),sh(2),sh(3)

  if(h1 < h2) then
    h = h1
  else
    h = h2
  endif
  if(h3 < h) h = h3

  params(P_CURRENTTIME) = 0.d0
  params(P_GRIDSPACING) = h
  params(P_TIMESTEP) = params(P_COURANT)*h   !CourantFactor*h

  write(6,*) 'spatial  resolution: ',params(P_GRIDSPACING)
  write(33,*)'spatial  resolution: ',params(P_GRIDSPACING)
  write(6,*) 'temporal resolution: ',params(P_TIMESTEP)
  write(33,*)'temporal resolution: ',params(P_TIMESTEP)
  
!  3.3 Grid hierarchy

!  3.4 Grid functions

  nx=sh(1)
  ny=sh(2)
  nz=sh(3)

  allocate(anm(nx,ny,nz), an(nx,ny,nz), anp(nx,ny,nz), &
       bnm(nx,ny,nz), bn(nx,ny,nz), bnp(nx,ny,nz), &
       cnm(nx,ny,nz), cn(nx,ny,nz), cnp(nx,ny,nz), &
       dnm(nx,ny,nz), dn(nx,ny,nz), dnp(nx,ny,nz), &
       temp(nx,ny,nz), phi(nx,ny,nz), temp(nx,ny,nz))

  phip => temp
  phim => temp

! ******************************************************* allocate arrays.

  ! -----------------------------------------------------------------------
  ! Sect. 4.0  User's stuff
  ! -----------------------------------------------------------------------

  !---------------------------------------------\
  ! Read the input file for the outlist module.  >
  !---------------------------------------------/
!  call ol_parser("outlist.in",params)
  
! BEGIN_COMPUTE

! Sect. 5.0 Initialize data

! ******************************************************* initialize data

  !--------------------------------------\
  ! Initial data for dynamical variables. >
  !--------------------------------------/
  call initialize(params, anm, an, anp, bnm, bn, bnp, cnm, cn, cnp, &
       dnm, dn, dnp, phi, temp)

  !-------------------------\
  ! Output the initial data  >
  !-------------------------/
  
!  call ol_output(0,m,d,nl,params)


! -----------------------------------------------------------------------
! Sect. 6.0 Integrate
! -----------------------------------------------------------------------

  t=0
  iter=0
  write(6,*) "Beginning time index: ",t
  write(33,*) "Beginning time index: ",t


  do while (iter < int(params(P_NSTEPS))) 
     
     call tstep(params, anm, an, anp, bnm, bn, bnp, cnm, cn, cnp, &
	dnm, dn, dnp, phim, phi, phip, temp)

     params(P_CURRENTTIME) = t * params(P_TIMESTEP)

     t = t+1
     iter= iter +1 

    !---------------------------------------------\
    ! Output the most recently computed timestep   >
    !---------------------------------------------/
!    call ol_output(t,m,d,nl,params)

     
  enddo

! Sect. 7.0 Finalize

! END_COMPUTE 

!  t = gft_close_all()
  close(33)

end program main
