program conjgdir
 
use sep 
use OW_type
use OW_grabvel_areal
use OW_exec_time
use OW_parms_areal
use OW_taper_areal
use OW_wei_areal
use solver_cjg_mod
use steepest_mod
use lloyd_vel

implicit none

character*80                                  :: hor1
real                                          :: start, finish, start1, finish1, start2, finish2
real                                          :: vmin,vmax,p_up,p_lo,epsx
integer                                       :: ni,i,ishot,j,k,iters,titer,stat1
real, allocatable                             :: velint(:),taper(:),pick1(:)
real, allocatable                             :: vel0(:)
logical                                       :: cond,bounds,cbounds
double precision, allocatable                 :: x(:)
real, allocatable                             :: l(:), u(:)

integer                                       :: rect1, rect2, rect3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 call sep_init()
 call sep_begin_prog()
 start=0.;finish=0.
 call owexec_time(start)
 call owparam()

 call from_param("rect1",rect1,5)
 call from_param("rect2",rect2,5)
 call from_param("rect3",rect3,5)

 call owtri_init(rect1,rect2,rect3)
 allocate(taper(kxm%n));taper=0.
 allocate(velint(v%z%n*v%xm%n),vel0(img%z%n*img%xm%n))
 allocate(x(img%z%n*img%xm%n));x=dble(0.)

! Set up the taper
 call owtaper(tapertr,taper)

! Read velocities
 call sreed("vel",velint,v%z%n*v%xm%n*4)
 call owgrab_vel(velint,vel0,.true.,ni)
 deallocate(velint)
 x=dble(vel0)

 call from_param("epsx",epsx,1.)
 call from_param("bounds",bounds,.false.)
 if (bounds) then
    allocate(l(img%z%n*img%xm%n));l=0.
    allocate(u(img%z%n*img%xm%n));u=0.
    call from_param("p_up",p_up)
    call from_param("p_lo",p_lo)
    if (p_up >= 1.) p_up=1+p_up/100
    if (p_lo >= 1.) p_lo=1-p_lo/100
    l=p_lo*x
    u=p_up*x
 end if 
!----------------------------------------------
! Horizon with no velocity update
 if (exist_file("horv0")) then
    call from_param('horv0',hor1)
    open(51,file=hor1,status='old')
    allocate(pick1(img%xm%n));pick1=9999.
    do i=1,img%xm%n
       read(51,*,end=111) pick1(i)
    end do
 111 close(51)
 end if
 titer=30
 call slvcjg_init(img%z%n*img%xm%n*nsecdim,vel0,taper,pick1)
 if (bounds) then
    call solver_cjg(epsx,titer,l,u)
 else
    call solver_cjg(epsx,titer)
 end if
 if (allocated(taper)) deallocate(taper)
 if (allocated(pick1)) deallocate(pick1)
 if (allocated(vel0)) deallocate(vel0)
 if (allocated(x)) deallocate(x)
 call sep_end_prog()

end program
