program Cod_
  !! convert to original depth
  !! after Stolt residual migration
  use sep
  use codmap_
  use util_

  !! complex input/output
  !! input:  z-hx-hy-mx-my
  !! output: z-hx-hy-mx-my
  !! need parameter rho

  implicit none

  type(myaxis) :: a1,a2,a3,a4,a5
  integer      :: i3,i4,i5
  real         :: rho
  logical      :: verb
  integer      :: nn,ii

  complex, dimension(:,:,:), allocatable :: din,dou
  real,    dimension(:,:,:), allocatable :: rin,rou
  real,    dimension(:,:,:), allocatable :: iin,iou

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

  !! ratio
  call from_param("rho",rho,1.)
  rho = 1./rho

  !! 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

  call sep_close()

  if(verb) write(0,*) "init COD map"
  call codmap_init(a1,a2,a3,rho)

  allocate(din(a1%n,a2%n,a3%n),rin(a1%n,a2%n,a3%n),iin(a1%n,a2%n,a3%n))
  allocate(dou(a1%n,a2%n,a3%n),rou(a1%n,a2%n,a3%n),iou(a1%n,a2%n,a3%n))

  nn=a5%n*a4%n*a3%n

  if(verb) write(0,*) "execute COD for rho=",rho
  do ii=1,nn
     if(verb) write(0,*) ii,nn
     call sep_read(din)
     rin =  real(din)
     iin = aimag(din)
     call codmap_execute(rin,rou)
     call codmap_execute(iin,iou)
     dou=cmplx(rou,iou)
     call sep_write(dou)
  end do

  deallocate(din,dou,rin,iin,rou,iou)
  call exit(0)
end program COD_
