program DECO_
  !! depth conversion
  use sep
  use decomap_
  use util_

  implicit none

  type(myaxis) :: a1,a2,a3,a4,a5
  type(myaxis) :: r1,r2,r3
  type(myaxis) :: b1

  integer :: i3
 
  real, dimension(:,:,:), allocatable :: din,dou,rho,map
  character(len=128) :: name

  call sep_init()

  !! mapping file
  name="rho"
  call ainaxis(r1,1,name)
  call ainaxis(r2,2,name)
  call ainaxis(r3,3,name)

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

  !! output axis
  b1%l="out"
  call pinaxis(b1)
  call houaxis(b1,1)

  call sep_close()
  
  if(r1%n /= a1%n) call seperr("field_n1 /= vel_n1")
  if(r2%o              > a2%o              ) call seperr("velocity too small for tt")
  if(r2%o+r2%d*(r2%n-1)< a2%o+(a2%n-1)*a2%d) call seperr("velocity too big for field 2")
  
  call decomap_init(a1,a2,a3,a4,a5,b1)
  
  ! axis conversion mapping
  allocate(map(r1%n,r2%n,r3%n))
  allocate(rho(r1%n,r2%n,r3%n))

  call sreed("rho",rho,size(rho)*4)
  call createmap(rho,map)        

  ! data conversion
  allocate(din(a1%n,a2%n,a3%n))
  allocate(dou(b1%n,a2%n,a3%n))

  do i3=1,a3%n
     call sep_read(din)
     call convert(din,dou)
     call sep_write(dou)
  end do

  deallocate(din,dou)
     
end program DECO_
