program OWlnsrch

use sep 
use OW_type
use OW_grabvel_areal
use OW_exec_time
use OW_parms_areal
!use fgevaluation1
use fgevaluation
use lloyd_vel
! using dips to smooth velocity
!use killwave
!use geodip 
!use OW_steer_parms
!use dip_prog

implicit none

character*80                                  :: hor1,log,listagem,listwolfe
real                                          :: start, finish, start1, finish1, start2, finish2
real                                          :: r,p_up,p_lo
integer                                       :: iter,titer,ni,i
real, allocatable                             :: velint(:),pick1(:)
real, allocatable                             :: B_image(:,:,:),vref(:,:,:),vel0(:),R_image(:,:,:),z0(:,:)
logical                                       :: cond,bounds,cbounds
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Variables of Line search
!
 integer                                   :: nmax, mmax ! nmax - dimension of the problem and mmax - max.# of ltd.memory corrections.
 double precision                          :: ALPHA,FVAL,DFVAL
 integer                                   :: NFE,NGE
 
 double precision, allocatable             :: x(:), g(:), g0(:), l(:), u(:), xplus(:)
 integer, allocatable                      :: nbd(:)
 double precision                          :: fun,alpha0
 double precision                          :: factr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 call sep_init()
 call sep_begin_prog()
 start=0.;finish=0.
 call owexec_time(start)
 call owparam()

 nmax=img%z%n*img%xm%n!*img%ym%n
 call from_param('mmax',mmax,10)
 call from_param("log",log,"log")
 open(unit=2,file=log)
 call from_param("list",listagem,"listagem")
 open(unit=3,file=listagem)
 call from_param("listwolfe",listwolfe,"listwolfe")
 open(unit=4,file=listwolfe)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Memory allocation
!

!allocate(l(nmax));l=dble(0.)
!allocate(u(nmax));u=dble(0.)
allocate(g(img%z%n*img%xm%n))
allocate(g0(img%z%n*img%xm%n))
allocate(velint(v%z%n*v%xm%n),vel0(img%z%n*img%xm%n))
allocate(x(img%z%n*img%xm%n))
allocate(xplus(img%z%n*img%xm%n))
allocate(B_image(img%xm%n,nsecdim,img%z%n))
allocate(R_image(img%xm%n,nsecdim,img%z%n))

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Read velocities
!
 call sreed("vel",velint,v%z%n*v%xm%n*4)
 call owgrab_vel(velint,vel0,.true.,ni)
 if (extrap(1:4).eq."pspi") then
    if(lloyd) then
      allocate(vref(nref,2,img%z%n))
      call lloyd_vel1(vel0,vref)
    end if
 end if
 deallocate(velint)

! Tolerances and constraints
 factr  = tiny(1.)
!----------------------------------------------
! Initial velocity model = background velocity
 x=vel0

!----------------------------------------------
! 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=1) pick1(i)
   end do
1 close(51)
end if

 iter=1
 titer = 1
 g=0.

 do while(iter <= niter)
     write(4,*) '********************'
     write(4,*) "ITERATION #",iter
     write(4,*) '********************'
     if (exist_file("horv0")) then
        call oweditvel(pick1,x,vel0)
     end if
     call to_history("n4",iter,"Velout")
     call srite("Velout",sngl(1/x),4*img%z%n*img%xm%n)
     B_image=0;fun=0d+0;g0=0d+0;R_image=0;xplus=0.
     if (exist_file("horv0")) then
        call wolfe_init(pick1)
     end if
     call owfeval(x,fun,R_image,B_image)!,0,iter)
     call to_history("n4",iter,"bkgimg")
     call srite("bkgimg",B_image,4*img%z%n*nsecdim*img%xm%n)
     if (envelope) then
        allocate(z0(img%xm%n,img%z%n))
        do i=1,img%z%n
           z0(:,i)=B_image(:,nsecdim/2,i)
!write(0,*) maxval(abs(z0(:,i)))
        end do       
        call owgeval(x,g0,R_image,z0)!,lsiter,iter)
        deallocate(z0)
     else
        call owgeval(x,g0,R_image)!,0,iter)
     end if
     call to_history("n4",iter,"grad")
     call srite("grad",sngl(g0),4*img%z%n*img%xm%n)
     alpha0=sqrt(dot_product(x,x)/dot_product(g0,g0))/10
     call wolfe_search2(x,g0,alpha0,dble(.001),dble(.9),fun,-g0,iter,xplus,alpha,fval,g,nfe,nge)
     if (exist_file("horv0")) then
        call wolfe_close()
     end if
     x = xplus
     iter=iter+1
 end do

 if (allocated(l)) deallocate(l)
 if (allocated(u)) deallocate(u)
 if (allocated(x)) deallocate(x)
 if (allocated(g)) deallocate(g)
 if (allocated(g0)) deallocate(g0)
 if (allocated(nbd)) deallocate(nbd)
 if (allocated(B_image)) deallocate(B_image)
 if (allocated(R_image)) deallocate(R_image)
 if (allocated(vref)) deallocate(vref)
 if (allocated(vel0)) deallocate(vel0)
 if (allocated(pick1)) deallocate(pick1)

 close(2)
 close(3)
 close(4)
11 format(a8,i4,a11,i5,a8,i4,2(a16,f20.10))

 call sep_end_prog()

end program OWlnsrch
