program Reverse_
  !! Paul Sava (paul@sep.stanford.edu)
  !! 3/11/04
  !! reverse axes 1 or 2
  !! real or complex input/output

  use sep
  use util_

  implicit none

  type(myaxis) :: a1,a2,a3,a4,a5
  integer      :: i1,i2,i3,i4,i5
  logical      :: verb
  integer      :: nn,ii
  integer      :: esize
  integer      :: axis

  complex, dimension(:,:,:), allocatable :: cin,cou
  real,    dimension(:,:,:), allocatable :: rin,rou

  integer :: ntot      !! total n2...n5
  integer :: ninc,iinc !! in core
  integer :: nblk,iblk !! no of blocks

  call sep_init()
  call from_param("verb",verb,.false.)
  call from_param("axis",axis,1)

  if(axis>2 .or. axis<1) call seperr("can only reverse axis 1 and 2")
  call from_history("esize",esize,4)

  !! input axes
  call hinaxis(a1,1)
  call hinaxis(a2,2)
  call hinaxis(a3,3)
  call hinaxis(a4,4)
  call hinaxis(a5,5)

  if(verb) then
     call printaxis(a1)
     call printaxis(a2)
     call printaxis(a3)
     call printaxis(a4)
     call printaxis(a5)
  end if

  !! how many slices to read in core
  call from_param("ninc",ninc,1)
  ntot = a3%n*a4%n*a5%n
  if(verb) then
     write(0,*) "ntot=",ntot
     write(0,*) "ninc=",ninc
  end if

  if(ntot/ninc * ninc == ntot) then
     nblk = ntot/ninc
  else
     if(verb) write(0,*) "WARNING: ninc does not divide ntot; reseting..."
     do while(ninc>=1 .and. ntot/ninc * ninc /= ntot)
        ninc=ninc-1
     end do
     nblk = ntot/ninc
     if(verb) write(0,*) "ninc=",ninc
  end if
  if(verb) write(0,*) "nblk=",nblk

  call sep_close()


  if(verb) write(0,*) "Reverse"
  if(esize==4) then
     allocate(rin(a1%n,a2%n,ninc),rou(a1%n,a2%n,ninc))

     if(axis==1) then
        do iblk=1,nblk
           if(verb) write(0,*) iblk,'/',nblk
           call sep_read(rin)
           rou(1:a1%n,:,:) = rin(a1%n:1:-1,:,:) 
           call sep_write(rou)
        end do
     end if

     if(axis==2) then
        do iblk=1,nblk
           if(verb) write(0,*) iblk,'/',nblk
           call sep_read(rin)
           rou(:,1:a2%n,:) = rin(:,a2%n:1:-1,:)
           call sep_write(rou)
        end do
     end if

     deallocate(rin,rou)
  else
     allocate(cin(a1%n,a2%n,ninc),cou(a1%n,a2%n,ninc))

     if(axis==1) then
        do iblk=1,nblk
           if(verb) write(0,*) iblk,'/',nblk
           call sep_read(cin)
           cou(1:a1%n,:,:) = cin(a1%n:1:-1,:,:) 
           call sep_write(cou)
        end do
     end if

     if(axis==2) then
        do iblk=1,nblk
           if(verb) write(0,*) iblk,'/',nblk
           call sep_read(cin)
           cou(:,1:a2%n,:) = cin(:,a2%n:1:-1,:)
           call sep_write(cou)
        end do
     end if

     deallocate(cin,cou)
  end if

  call exit(0)
end program REVERSE_
