program OWlbfgsc_adj

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 OW_imageIO
use OW_tim2ang
use solver_lbfgsb_mod
use base_spline_mod
use fx2d
use lloyd_vel
! using dips to smooth gradient
use killwave
!use dip_prog

implicit none

character*80                                  :: hor1,log,listagem
real                                          :: start, finish, start1, finish1, start2, finish2
real                                          :: a1,b1,vmin,vmax,stat1,v0,r
integer                                       :: ixs,ino,ni,i,st,ishot,t3,iz,j,k,iter,ix,ish,titer,i1,i2,i3
real, allocatable                             :: velint(:),taper(:),pick1(:)
real, allocatable                             :: deltaI(:,:,:),scatS(:),deltaS(:),illum(:,:),illumination(:,:)
real, allocatable                             :: R_image(:,:,:),R_image1(:,:),B_image(:,:,:)
real, allocatable                             :: vel(:),vref(:,:,:),vel0(:),indip(:)
complex, allocatable                          :: recv(:,:),source(:,:),source1(:)
logical                                       :: cond,bounds
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Variables of L-BFGSC
!
 integer                       :: nmax, mmax ! nmax - dimension of the problem and mmax - max.# of ltd.memory corrections.
 
!     Declare the variables needed by the code.
!     A description of all these variables is given at the end of 
!     the driver.
 
 character*70                  :: task, csave
 logical                       :: lsave(4)
 integer                       :: iprint,subiter
 integer,allocatable           :: nbd(:), wwa(:)
 integer                       :: isave(44)
 double precision              :: fun, factr, pgtol
 double precision              :: dsave(29)
 double precision, allocatable :: x(:), l(:), u(:), g(:),wa(:),xspl(:),gspl(:)
 real, allocatable             :: spl(:),spl2(:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 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)

 if (axkn.ne.99999 .or. azkn.ne.99999) then
    stat1 = base_spln2d_init(img%xm,img%z,axk,azk)
    allocate(gspl(nmax));gspl=dble(0.)
!    stat1 = base_spln3d_init(img%xm,img%ym,img%z,axk,ayk,azk)
 end if

 if (dipflt) then
    allocate(indip(img%xm%n*1*img%z%n))
!    allocate(indip(img%xm%n*img%ym%n*img%z%n))
    call sreed("indip",indip,4*img%xm%n*1*img%z%n)
 end if

 call from_param("log",log,"log")
 open(unit=2,file=log)
 call from_param("list",listagem,"listagem")
 open(unit=3,file=listagem)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Memory allocation
!

allocate(l(nmax));l=dble(0.)
allocate(u(nmax));u=dble(0.)
allocate(nbd(nmax));nbd=0
allocate(wa(2*mmax*nmax+4*nmax+12*mmax*mmax+12*mmax));wa=dble(0.)
allocate(wwa(3*nmax));wwa=0
allocate(deltaS(img%z%n*kxm%n));deltaS=0.
allocate(scatS(img%z%n*img%xm%n));scatS=0.
allocate(deltaI(kxm%n,nsecdim,img%z%n));deltaI=0.
allocate(B_image(img%xm%n,nsecdim,img%z%n));B_image=0.
allocate(taper(kxm%n));allocate(g(img%z%n*img%xm%n));g=dble(0.)
allocate(velint(v%z%n*v%xm%n),vel(img%z%n*kxm%n),vel0(img%z%n*img%xm%n))
allocate(spl2(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)
 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)

!     We wish to have output at every iteration.
 iprint = 99

! Tolerances and constraints
 factr  = 1.0d0+3
 pgtol  = 1.0d-16
 call from_param("bounds",bounds,.false.)
 if (bounds) then
    call from_param("vmin",vmin)
    call from_param("vmax",vmax)
    l=dble(1/vmax)
    u=dble(1/vmin)
    nbd=2
 else
    nbd=0
 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=1) pick1(i)
   end do
1 close(51)
end if
!----------------------------------------------
! Initial velocity model = background velocity
 x=dble(vel0)

 iter=0
 subiter=0
 titer = 1
 task = 'START'
 g=dble(0.)
 cond=.true.

 do while(cond .and. iter <= niter)
     cond=.false.
     call setulb(nmax,mmax,x,l,u,nbd,fun,g,factr,pgtol,wa,wwa,task,iprint,csave,lsave,isave,dsave)
     if (task(1:2) .eq. 'FG') then
!        the minimization routine has returned to request the
!        function fun and gradient g values at the current x.
!        Compute function value fun for the sample problem.
!
!     x is the updated slowness
!
       cond=.true.

       if (exist_file("horv0")) then
          call oweditvel(pick1,x,vel0)
       end if

       subiter=subiter+1

       call to_history("n4",titer,"Velout")
       call srite("Velout",sngl(1/x),4*img%z%n*img%xm%n)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Computing the gradient
!

       start1=0.;finish1=0.;g=dble(0.)

       call sseek("rec",0,0)
       call sseek("sou",0,0)
       call owexec_time(start1)
       if (not(areal)) then
          if (allocated(source1)) deallocate(source1); allocate(source1(sou%w%n))
          source1=0;
          if (allocated(source)) deallocate(source); allocate(source(rec%h%n,sou%w%n))
          source=0;
          call sreed("sou",source1,8*sou%w%n)
          source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
       end if

       start1=0.;finish1=0.;scatS=0.
       call owexec_time(start1)
       do ishot=1,rec%xd%n
          start=0.;finish=0. 
          call owexec_time(start)
          if (areal) then
             if (allocated(source)) deallocate(source); allocate(source(sou%h%n,sou%w%n))
             source=0;
             call sreed("sou",source,8*sou%h%n*sou%w%n)
             vel=x
             ni=1
          else
             call owgrab_vel(sngl(x),vel,.false.,ni,ishot)
          end if
             
          if (allocated(recv))   deallocate(recv);   allocate(recv(rec%h%n,rec%w%n))
          recv=0;deltaS=0
          call sreed("rec",recv,8*rec%h%n*rec%w%n)

          call owwemva_init(ishot)
          if (lloyd) then
             call owwemva_adj(taper,recv,source,vel,fun,deltaS,deltaI,vref)
          else
             call owwemva_adj(taper,recv,source,vel,fun,deltaI,deltaS)
          end if
          if (areal) then
             scatS=scatS+deltaS
             B_image=B_image+deltaI
          else
             do j=ni,ni+kxm%n-1;do k=1,img%z%n
                scatS((k-1)*img%xm%n+j)=scatS((k-1)*img%xm%n+j)+deltaS((k-1)*kxm%n+j-ni+1)
             end do;end do
             B_image(ni:ni+kxm%n-1,:,:)=B_image(ni:ni+kxm%n-1,:,:)+deltaI
          end if

          call owexec_time(finish)
          write(3,111) "EXTRAPOLATED WAVEFIELDS = ",ishot," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do

       call to_history("n4",titer,"Dimage")
       call srite("Dimage",B_image,4*img%z%n*nsecdim*img%xm%n)

       g=-dble(scatS)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Scaling/editing the gradient
!          spl2=sngl(g)
!          stat1 = fx2d_lop (img%z%n,1,img%xm%n,5,50,5,int(1.5*img%z%n/(5+1)),&
!                           int(1.5*img%xm%n/(50-5+1)),40.,spl2,spl)
!          g=dble(spl)
!       g = g/maxval(abs(g))!*minval(abs(x))
       write(0,*) "main => ",maxval(g),minval(g)

       if (axkn.ne.99999 .or. azkn.ne.99999) then
          call owexec_time(start2)
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          ! B-spline representation of grad
          gspl=dble(0.)
          stat1 = base_spln2d_oper(.true. , .false., gspl, g)
          stat1 = base_spln2d_oper(.false. ,.false., gspl, g)
          call owexec_time(finish2)
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
       end if

       if (dipflt) then
!          call ow_dips(sngl(g),indip)
          a1=sqrt(sum(g*g)/size(g))
          call ow_kill(g,indip)
          b1=sqrt(sum(g*g)/size(g))
          g=g*a1/b1
       end if

       write(3,112) "SMOOTHED GRADIENT IN ",(finish2-start2)/60

!       if (titer < 30 ) then
!          g=g*(1-exp(-float(titer)/damp))
!       else
!          g=g*(1-exp(-float(titer)))
!       end if
!
!
       if (exist_file("horv0")) then
          call oweditvel(pick1,g)
       end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       call owexec_time(finish1)
       write(3,112) "COMPUTED THE GRADIENT IN ",(finish1-start1)/60," min OF ITERATION # ", iter
       call to_history("n4",titer,"grad")
       call srite("grad",sngl(g),4*img%z%n*img%xm%n)

       write(2,113) " iter = ",iter," titer= ",titer," task = ",task(1:5)," Obj.function = ", sngl(fun)
       titer=titer+1
       cond=.true.
    end if

    if (task(1:5) .eq. 'NEW_X')  then
       iter=iter+1
       subiter=0
       cond=.true.
    end if

 end do
 111 format(a26,i4,a8,i4,a4,f6.2,a4)
 112 format(a29,f10.2,a20,i4)
 113 format(2(a8,i4),a8,a5,a16,f20.6)
 if (allocated(l)) deallocate(l)
 if (allocated(u)) deallocate(u)
 if (allocated(nbd)) deallocate(nbd)
 if (axkn.ne.99999 .or. azkn.ne.99999) then
    if (allocated(xspl)) deallocate(xspl)
    if (allocated(gspl)) deallocate(gspl)
 end if
 if (allocated(wa)) deallocate(wa)
 if (allocated(wwa)) deallocate(wwa)
 if (allocated(source)) deallocate(source)
 if (allocated(recv)) deallocate(recv)
 if (allocated(deltaS)) deallocate(deltaS)
 if (allocated(scatS)) deallocate(scatS)
 if (allocated(deltaI)) deallocate(deltaI)
 if (allocated(R_image)) deallocate(R_image)
 if (allocated(R_image1)) deallocate(R_image1)
 if (allocated(B_image)) deallocate(B_image)
 if (allocated(taper)) deallocate(taper)
 if (allocated(vel)) deallocate(vel)
 if (allocated(vref)) deallocate(vref)
 if (allocated(vel0)) deallocate(vel0)
 if (allocated(illum)) deallocate(illum)
 if (allocated(illumination)) deallocate(illumination)

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

 call sep_end_prog()

end program OWlbfgsc_adj
