
!---------------------------------------------------------------------
      program btdv3

      include 'header3d.h'
      integer  i,niter,step,fstatus,n3
      double precision  navg,mflops
      external timer_read,verify
      double precision  tmax,timer_read
      logical  verified
      character  class

!---------------------------------------------------------------------
!      Root node reads input file (if it exists) else takes
!      defaults from parameters
!---------------------------------------------------------------------
      write (unit = *,fmt = 1000) 
      open (unit = 2,file = 'inputbt.data',status = 'old',iostat = fstat
     &us)
      if (fstatus .eq. 0) then
         write (unit = *,fmt = 233) 
233      format(' Reading from input file inputbt.data')
         read (unit = 2,fmt = *) niter
         read (unit = 2,fmt = *) dt
         read (unit = 2,fmt = *) grid_points(1),grid_points(2),grid_poin
     &ts(3)
         close (unit = 2)
      else  
         write (unit = *,fmt = 234) 
         niter = niter_default
         dt = dt_default
         grid_points(1) = problem_size
         grid_points(2) = problem_size
         grid_points(3) = problem_size
      endif  
234   format(' No input file inputbt.data. Using compiled defaults')
      write (unit = *,fmt = 1001) grid_points(1),grid_points(2),grid_poi
     &nts(3)
      write (unit = *,fmt = 1002) niter,dt
1000  format(//, ' NAS Parallel Benchmarks 3.3.1 - DVMH version',' - BT 
     &Benchmark ',/)
1001  format(' Size: ', i3, 'x', i3, 'x', i3)
1002  format(' Iterations: ', i3, '    dt: ', F10.6)
      if (grid_points(1) .gt. imax .or. grid_points(2) .gt. jmax .or. gr
     &id_points(3) .gt. kmax) then
         print *, (grid_points(i), i = 1,3)
         print *, ' Problem size too big for compiled array sizes'
         goto 999
      endif  
      open (unit = 2,file = 'inputStage',status = 'old',iostat = fstat
     &us)	
      if (fstatus .eq. 0) then
         read (unit = 2,fmt = *) stage_n
         close (unit = 2)
      else  
         stage_n = 0
      endif 
      write(*,*) 'stage = ', stage_n
	  
      call set_constants()     
      call initialize()	
      call exact_rhs()
	  
! ************* DO 2 iterations for touch all code 
       call adi_first
       call adi_first
       call initialize

      call timer_clear(1)
      call timer_start(1)
      do  step = 1,niter
         if (mod (step,20) .eq. 0 .or. step .eq. 1) then
            write (unit = *,fmt = 200) step
200         format(' Time step ', i8)
         endif  
         call adi()
      enddo  
      call timer_stop(1)
      tmax = timer_read (1)
      call verify(niter,class,verified)
      n3 = grid_points(1) * grid_points(2) * grid_points(3)
      navg = (grid_points(1) + grid_points(2) + grid_points(3)) / 3.0
      if (tmax .ne. 0.) then
         mflops = 1.0e-6 * float (niter) * (3478.8 * float (n3) - 17655.
     &7 * navg** 2 + 28023.7 * navg) / tmax
      else  
         mflops = 0.0
      endif  
      call print_results('BT',class,grid_points(1),grid_points(2),grid_p
     &oints(3),niter,tmax,mflops,'          floating point',verified,npb
     &version)

!       ,compiletime, cs1, cs2, cs3, cs4, cs5,cs6, '(none)')
999   continue
      end
	  
      subroutine adi_first()
      call compute_rhs()  
      call x_solve()  
      call y_solve()  
      call z_solve()
      return
      end
	  
      subroutine adi ()

!DVM$ interval 1
      call compute_rhs()
!DVM$ end interval	  
!DVM$ interval 11	  
      call x_solve()
!DVM$ end interval	  
!DVM$ interval 12	  
      call y_solve()
!DVM$ end interval	  
!DVM$ interval 13	  
      call z_solve()
!DVM$ end interval
      return
      end
