program nlin_dir

use sep 
use OW_type
use OW_grabvel2_areal
use OW_exec_time
use OW_parms_areal
use OW_wei_areal
use steepest_mod
use base_spline_mod


implicit none

real                                          :: start, finish, start1, finish1, start2, finish2,alpha0
integer                                       :: ni,i,ishot,j,k,titer,stat1
integer                                       :: ir,is,iww
real, allocatable                             :: velint(:),vel0(:),pick1(:),velout(:)
double precision, allocatable                 :: x(:), g(:),fun(:),fn(:)
real, allocatable                             :: l(:), u(:)
integer                                       :: mpierr,i_r,ierr,status,shtida,ils,imin

real                                          :: lambda,norm_dx,minxx
double precision                              :: ft,f2,gt,gn,g2,ngrad0,ngrad1,teta,fmin
double precision, allocatable                 :: xn(:)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 call sep_init()
 call owparam()

 allocate(x(img%z%n*img%xm%n*img%ym%n));x=dble(0.)
 allocate(g(img%z%n*img%ym%n*img%xm%n));g=dble(0.)
 allocate(pick1(img%xm%n*img%ym%n))
 allocate(xn(img%z%n*img%xm%n*img%ym%n));xn=dble(0.)
 allocate(vel0(nzold*img%ym%n*img%xm%n));vel0=0
 allocate(velout(nzold*img%ym%n*img%xm%n));velout=0
 allocate(fun(max_iters));fun=0
 allocate(fn(nls));fn=0 
write(0,*) "Main: vectors allocated"

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Read velocities
    allocate(velint(v%z%n*v%xm%n*v%ym%n));velint=0
    call sreed("vel",velint,v%z%n*v%xm%n*v%ym%n*4)
    call owgrab_velim(velint,vel0,ni)
    deallocate(velint)
    !----------------------------------------------
    ! Horizon with no velocity update

    if (exist_file("horv0")) then
       open(51,file=hor1,status='old')
       do i=1,img%xm%n*img%ym%n
          read(51,*,end=111) pick1(i)
       end do
111    close(51)
    else
       pick1=0
    end if


write(0,*) "==================================="
write(0,*) "Beginning of the optimization loops"
write(0,*) "==================================="

!--------------------------------------------
 ! Initial velocity model = background velocity
 do k=nzmin,nzmax;do j=1,img%ym%n;do i=1,img%xm%n
    x((k-nzmin)*img%ym%n*img%xm%n+(j-1)*img%xm%n+i)=dble(vel0((k-1)*img%ym%n*img%xm%n+(j-1)*img%xm%n+i))
 end do;end do;end do

 call to_history("n4",1,"Velout")
 call srite("Velout",vel0,4*nzold*img%xm%n*img%ym%n)

 xn=x
 minxx=sngl(minval(x))
 if (bounds) then
    allocate(l(img%z%n*img%ym%n*img%xm%n));l=0.
    allocate(u(img%z%n*img%ym%n*img%xm%n));u=0.
    if (cbounds) then
       l=vmin
       u=vmax
    else
       l=p_lo*1/x
       u=p_up*1/x
    end if
 end if 

 call steep_mpi_init(pick1,minxx)
 if (scaling) call scale_init()
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Main code loop starts here  
!! 
 alpha0=alpha
 fun=0;gn=0;g=0.
 call fgeval_init(1,sngl(xn))
 call scale_init1(.true.)

 call fgeval(vel0=sngl(xn),ff=fun(1),gg=gn,grad=g,ngrad=ngrad0)
 write(3,*) 'Iteration = ',1, " fn = ",fun(1),0
 lambda=1.
 i=1
800 x=xn+lambda*g
    if (any(mask=(abs(1/x - 1/xn) > deltav/xn)).and.i < 40) then
        lambda=factor*lambda
        x=xn+lambda*g
        i=i+1
        if (i > 40) goto 1000
        goto 800
    end if
1000 g=g*lambda
 write(3,*) "lambda = ",lambda,i
 call to_history("n4",1,"grad_")
 call grad_io("grad_","rite",1,sngl(g))
 titer=2
 do iters =2, max_iters
    fn=0.;gn=0
    do ils=1,nls
       x = xn + alpha*ils*g
       write(3,*) ils,"Max/min vel.",int(maxval(1/x)),int(minval(1/x))
       call scale_init1(.false.)
       call fgeval(vel0=sngl(x),ff=fn(ils),gg=gn)
    end do
    if (wpar%i_mpi==0) write(0,*) iters," fn_vector = ", fn
    fmin = fn(1)
    imin = 1
    do ils=1,nls
       if (fn(ils).le.fmin) then
          fmin = fn(ils)
          imin = ils
       end if
    end do

    if (fmin <= fun(iters-1)) then
       x = xn + alpha*imin*g
       xn = x
       do k=nzmin,nzmax;do j=1,img%ym%n;do i=1,img%xm%n
          vel0((k-1)*img%ym%n*img%xm%n+(j-1)*img%xm%n+i)=sngl(xn((k-nzmin)*img%ym%n*img%xm%n+(j-1)*img%xm%n+i))
       end do;end do;end do
       call to_history("n4",iters,"Velout")
       call srite("Velout",vel0,4*nzold*img%xm%n*img%ym%n)
       g=0.
       call fgeval_init(iters)
       call scale_init1(.false.)
       call fgeval(vel0=sngl(xn),ff=fun(iters),gg=gn,grad=g,ngrad=ngrad1)
       g=g*lambda*fun(iters)/fun(1)
       call to_history("n4",titer,"grad_")
       call grad_io("grad_","rite",iters,sngl(g))
!       alpha=alpha * 0.5
       titer=titer+1
       ! ****************************************************************** 
       ! Line search complete 
       ! ****************************************************************** 

       write(3,*) 'Iteration = ',iters, " fn = ",fun(iters),alpha

    else
       alpha = alpha * 0.5
    end if

 end do

01  format(a9,i7,a6,i9)
02  format(a9,i7,2(a6,i9))
03  format(a9,i7)
04  format(a12,i4,a7,i11)
1212 close (3)
 if (scaling) call scale_close()
 if (allocated(x)) deallocate(x)
 if (allocated(g)) deallocate(g)
 if (allocated(xn)) deallocate(xn)
 if (allocated(vel0)) deallocate(vel0)
 close(3)
 call steep_mpi_close()
 call sep_close()
 stop
end program
