
program vel_convert
use tau_z_conv_3d_mod
use sep3d_struct_mod
use sep
implicit none
type(sep3d)        :: in,out
integer            :: ndim
character(len=128) :: op
integer            :: nin(3),nout(3)
real               :: oin(3),din(3),oout(3),dout(3)
real,pointer       :: ina(:,:,:),outa(:,:,:),map(:,:,:),extra(:,:,:)
integer, external  :: sreed,srite
logical            :: outs


call init_sep3d("in",in,"INPUT")
call from_param("op",op)
ndim=sep3d_ndims(in)
if( ndim < 2 .or. ndim >3) call seperr("expecting 2 or 3-D data")
oin=0.;oin(1:ndim)=in%o(1:ndim)
din=1.;din(1:ndim)=in%d(1:ndim)
nin=1 ;nin(1:ndim)=in%n(1:ndim)
nout=nin
oout=oin
dout=din
call from_param("n1",nout(1))
call from_param("o1",oout(1))
call from_param("d1",dout(1))
allocate(ina (nin(1) ,nin(2) ,nin(3)))
allocate(map(nout(1),nout(2),nout(3)))
outs=.false.

if(4*product(nin)/=sreed("in",ina,size(ina)*4)) &
  call seperr("trouble reading in input  velocity")


if (op(1:5)=="vel_z" .or. op(1:7)=="field_t" ) then !convert velocity to time
  outs=.true.
  allocate(outa(nout(1),nout(2),nout(3)))
  if(.not. init_tau_z_conv_3d(nout,oout,dout,nin,oin,din)) &
    call seperr("trouble initializing conversion")
  if(.not. create_z_map_3d(ina,map)) call seperr("trouble creating z map")
  if ( op(1:7)=="field_t") then
    if(4*size(ina)/=sreed("field_in",ina,size(ina)*4)) &
      call seperr("trouble reading in field")
  end if
  if(.not. convert_field_tau_z_3d(ina,outa)) call seperr("trouble converting")
else 
  outs=.true.
  if(.not. init_tau_z_conv_3d(nin,oin,din,nout,oout,dout))&
    call seperr("trouble initializing conversion")
  if(.not. create_z_map_3d(ina,map)) call seperr("trouble creating tau map")
  if (op(1:5)=="vel_t") then !convert velocity to depth
    outs=.true.
    allocate(outa(nout(1),nout(2),nout(3)))
    if(.not. convert_field_z_tau_3d(ina,outa)) call seperr("trouble converting")
  else if (op(1:6)=="sigx_z") then !calculate sigmax in depth space
    allocate(outa(nin(1),nin(2),nin(3)))
    if(.not. calc_sigmax_3d(ina,outa)) call seperr("trouble calc sigmax")
  else if (op(1:6)=="sigy_z") then !calculate sigmay in depth space
    allocate(outa(nin(1),nin(2),nin(3)))
    if(.not. calc_sigmay_3d(ina,outa)) call seperr("trouble calc sigmay")
  else 
    outs=.true.
    allocate(outa(nout(1),nout(2),nout(3)))
    allocate(extra(nin(1),nin(2),nin(3)))
    if (op(1:5)=="sigx_t") then !calculate sigmax in tau space
      if(.not. calc_sigmax_3d(ina,extra)) call seperr("trouble calc sigmax")
      if(.not. convert_field_z_tau_3d(extra,outa))call seperr("trouble convert")
    else if (op(1:5)=="sigy_t") then !calculate sigmay in tau space
      if(.not. calc_sigmay_3d(ina,extra)) call seperr("trouble calc sigmay")
      if(.not. convert_field_z_tau_3d(extra,outa))call seperr("trouble convert")
    else 
      write(0,*) "op=",trim(op)
      call seperr("unknown operation")
    end if
  end if
end if
if(outs) then
  call init_sep3d(in,out,"OUTPUT")
  out%n(1)=nout(1)
  out%n(2)=nout(2)
  out%n(3)=nout(3)
  call sep3d_set_sep3d(out)
  call sep3d_write_description("out",out)
end if
  if(4*size(outa)/=srite("out",out,outa)) &
     call seperr("trouble writing out data")
end program







