program Veltran                       ! Velocity transform
  use sep
  use velsimpx
  use solver_tiny_mod
  use cgstep_mod

  implicit none

  integer                          :: nt, nx, velsimpx_lopns, iadj, iinv, ilint, n2, niter, stat, ns
  real                             :: o1,d1, o2,d2, x0,dx, s0,ds
  real, dimension(:), allocatable  :: data, model, error,resid
  logical                          :: adj, inv, lint

  call sep_init()
  call from_history( "n1", nt       )
  call from_history( "n2", n2       )
  call from_history( "o1", o1       )
  call from_history( "d1", d1       )
  call from_either(  "nx", nx, n2   )
  call from_either(  "ns", ns, 2*nx )
  call from_either(  "s0", s0       ) 
  call from_either(  "ds", ds       )  
  call from_either(  "x0", x0       )
  call from_either(  "dx", dx       ) 

  call from_param( "adj",  iadj, 1);   adj  = (iadj==1)
  call from_param( "inv",  iinv, 0);   inv  = (iinv==1)
  call from_param( "lint", ilint, 0);  lint = (ilint==1)

  if (adj) then 
     call to_history ("label2","slowness(sec/km)")
     n2 = ns; o2 = s0; d2 = ds
  else
     call to_history ("label2","distance(km)")
     n2 = nx; o2 = x0; d2 = dx
  end if
  call to_history ("o2", o2)
  call to_history ("d2", d2)
  call to_history ("n2", n2)
  
  call from_param( "niter", niter, 1)
  if (exist_file ("error")) call to_history ("n1",niter,"error")
  call sep_close()
  
  allocate( data(nt*nx), model(nt*ns), resid(nt*nx) )

  call velsimpx_init( d1, d1, x0, dx, s0, ds, nt, nx, ns, lint)

  if( adj) then  
     allocate (error (niter))
     call sep_read (data)
     if( inv) then		              ! Do LS inverse...	
        stat = velsimpx_lop( adj, .false., model, data)
!	To test different initial models --> null space
!	model=sin(model)
!	model=model+1000.
        call solver_tiny(m=model, d=data, &
        Fop=velsimpx_lop, stepper=cgstep, &
        niter=20,resd=resid)
     write (0,*) 'resid = ',sqrt(dot_product(resid,resid))
     else                                     ! ...or simple adjoint.
        stat = velsimpx_lop( adj, .false., model, data)
     end if
     call sep_write( model)
     if (exist_file ("error")) call sep_write (error/dot_product (data,data), "error")
     deallocate (error)
  else
     call sep_read (model)
     stat = velsimpx_lop (adj, .false., model, data)
     call sep_write (data)
  end if

  deallocate( model, data)

  call exit (0)
end program Veltran
