program OWlbfgsc

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_dotst
use lbfgsb_mod

implicit none

logical						:: icwind1
real						:: start, finish, start1, finish1, o3,a,b,drho
integer						:: ixs,ino,ni,i,st,idip,n3,t3,iz,j,k,iter,ix,ish
real, allocatable					:: velint(:,:),taper(:),deltaSint(:,:)
real, allocatable					:: deltaI(:,:,:),R_image(:,:,:),B_image(:,:,:)
real, allocatable					:: scatS(:,:)
real, allocatable					:: vel(:,:),vref(:,:,:)
complex, allocatable				:: recv(:,:),source(:,:),source1(:,:),tfld3(:,:,:)
real,allocatable					:: deltaS(:,:),temp(:,:,:)
logical						:: cond

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! 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*60                  :: task, csave
 logical                       :: lsave(4)
 integer                       :: n, m, iprint
 integer,allocatable           :: nbd(:), iwa(:)
 integer                       :: isave(44)
 double precision              :: f, factr, pgtol
 double precision              :: dsave(29)
 double precision, allocatable :: x(:), l(:), u(:), g(:),wa(:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 call sep_init()
 call sep_begin_prog()
 open(unit=2,file='lixo')
 open(unit=3,file='lixo_vel')

 start=0.;finish=0.
 call owexec_time(start)
 call owparam()

 nmax=v%z%n*v%xm%n
 mmax=5

 allocate(x(nmax),l(nmax),u(nmax),g(nmax),nbd(nmax))
 allocate(wa(2*mmax*nmax+4*nmax+12*mmax*mmax+12*mmax))
 allocate(iwa(3*nmax))

 allocate(source1(rec%h%n,sou%w%n))
 allocate(source(sou%h%n,sou%w%n))
 allocate(recv(rec%h%n,rec%w%n))
 allocate(deltaS(img%z%n,img%xm%n));deltaS=0.
 allocate(scatS(img%z%n,img%xm%n));scatS=0.
 allocate(deltaI(img%z%n,img%sh%n,img%xm%n));deltaI=0.
 allocate(R_image(img%z%n,img%sh%n,img%xm%n));R_image=0.
 allocate(B_image(img%z%n,img%sh%n,img%xm%n));B_image=0.
 allocate(taper(img%xm%n))
 call owtaper(tapertr,taper)

 allocate(velint(v%z%n,v%xm%n),vel(img%z%n,img%xm%n))
 allocate(vref(nref,2,v%z%n))
 call sreed("vel",velint,v%z%n*v%xm%n*4)
 if(exist_file("vref")) then
    call sreed("vref",vref,v%z%n*2*nref*4)
 end if
 call owgrab_vel(velint,vel,ni)
 deallocate(velint)

 if (resmig.or.resmig_lin) then
    allocate(tfld3(img%z%n,img%sh%n,img%xm%n))
    call owfftw3(tfld3)
 end if

 call owexec_time(finish)
 write(0,'(8x,a,f10.2,a)')' Init.phase - Total elapsed time:', (finish - start),' (sec)'

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

! Tolerances and constraints
 factr  = 1.0d+3
 pgtol  = 1.0d-16
 l=dble(1/2.05)
 u=dble(1/1.95)
 nbd=0

 drho=1-rho
 write(0,*) " ***** DRHO ******     " , drho

! Initial velocity model = background velocity
 do j=1,img%xm%n
    do k=1,img%z%n
       x((j-1)*img%z%n+k)=dble(1/vel(k,j))
    end do
 end do
 iter=0
 task = 'START'

 cond=.true.
 icwind1=.false.

 do while(cond .or. iter <= niter)

    cond=.false.

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

       cond=.true.

       do j=1,img%xm%n
          do k=1,img%z%n
             vel(k,j)=1/x((j-1)*img%z%n+k)
          end do
       end do
       if (exist_file("Velout")) then
          call srite("Velout",vel(:,tapertr+padtraces+1:img%xm%n-tapertr-padtraces), &
                    4*img%z%n*(img%xm%n-2*(tapertr+padtraces)))
       end if

! Computing bkg.image
       start1=0.;finish1=0.;
       B_image=0.;R_image=0.;
       call sseek("in",0,0)
       call sseek("sou",0,0)
       call owexec_time(start1)
       if (sou%xd%n.eq.1) then
          call sreed("sou",source,8*sou%h%n*sou%w%n)
       end if
       do idip=1,rec%xd%n
          start=0.;finish=0. 
          call owexec_time(start)
          recv=0;deltaI=0
          if (sou%xd%n.eq.1) then
             source1=0;
             ino=int(((idip-1)*rec%xd%d+rec%xd%o-img%xm%o- &
                 (padtraces+tapertr)*img%xm%d)/img%xm%d)+1.5
             source1(ino,:)=source(1,:)
          else
             source=0;
             call sreed("sou",source,8*sou%h%n*sou%w%n)
             source1=source
          end if
          call sreed("in",recv,8*rec%h%n*rec%w%n)
          if (exist_file("vref")) then
             call owwei(taper,recv,source1,vel,deltaI,vref)
          else
             call owwei(taper,recv,source1,vel,deltaI)
          end if
          B_image=B_image+deltaI
          call owexec_time(finish)
          write(0,111) "EXTRAPOLATED WAVEFIELDS = ",idip," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do

       call owexec_time(finish1)
       write(0,112) "COMPUTED BACKGROUND IMAGE IN ",(finish1-start1)/60," min OF ITERATION # ", iter

       if (exist_file("bkgimg")) then
          call srite("bkgimg",real(B_image(:,:,tapertr+padtraces+1:img%xm%n-tapertr-padtraces)), &
                    4*img%z%n*img%sh%n*(img%xm%n-2*(tapertr+padtraces)))
       end if

! Computing the perturbed image

       if (dso) then
          call owdso(.false., B_image, R_image)
          B_image=0.
          call owdso(.true. , B_image, R_image)
          if (exist_file("Dimage")) then
             call srite("Dimage",real(B_image(:,:,tapertr+padtraces+1:img%xm%n-tapertr-padtraces)), &
                       4*img%z%n*img%sh%n*(img%xm%n-2*(tapertr+padtraces)))
          end if
       else if (resmig_lin) then
          tfld3=B_image
           do ix=1,kxm%n
             do ish=1,img%sh%n
                tfld3(:,ish,ix)=-cmplx(0.,1.)*tfld3(:,ish,ix)*img%z%all
             end do
          end do
          call owfft3Fexe(tfld3)
          tfld3=tfld3/sqrt(1.*kz%n*ksh%n*kxm%n)
          call ow_resmig_lin(tfld3, rho)
          call owfft3Iexe(tfld3)
          tfld3=tfld3/sqrt(1.*kz%n*ksh%n*kxm%n)
          B_image=real(tfld3)*drho
          if (exist_file("Dimage")) then
             call srite("Dimage",real(B_image(:,:,tapertr+padtraces+1:img%xm%n-tapertr-padtraces)), &
                       4*img%z%n*img%sh%n*(img%xm%n-2*(tapertr+padtraces)))
          end if
       else if (resmig) then
          tfld3=B_image
          call owfft3Fexe(tfld3)
          tfld3=tfld3/sqrt(1.*kz%n*ksh%n*kxm%n)
          call ow_resmig(tfld3, rho)
          call owfft3Iexe(tfld3)
          tfld3=tfld3/sqrt(1.*kz%n*ksh%n*kxm%n)
          B_image=B_image-tfld3
          if (exist_file("Dimage")) then
             call srite("Dimage",real(B_image(:,:,tapertr+padtraces+1:img%xm%n-tapertr-padtraces)), &
                       4*img%z%n*img%sh%n*(img%xm%n-2*(tapertr+padtraces)))
          end if
       end if  
 
       if (impulse) then
          allocate(temp(img%z%n,img%sh%n,2))
          temp(:,:,1)=B_image(:,:,int(img%xm%n/2))
          B_image=0
          B_image(:,:,int(img%xm%n/2))=temp(:,:,1)
          deallocate(temp)
       end if

       call sseek("in",0,0)
       call sseek("sou",0,0)
       if (sou%xd%n.eq.1) then
          call sreed("sou",source,8*sou%h%n*sou%w%n)
       end if

! Computing the gradient
       start1=0.;finish1=0.;scatS=0.
       call owexec_time(start1)
       if (icwind) then
          icwind1=icwind
          icwind=.false.
       end if
       do idip=1,rec%xd%n
          start=0.;finish=0. 
          call owexec_time(start)
          recv=0;deltaS=0
          if (sou%xd%n.eq.1) then
             source1=0;
             ino=int(((idip-1)*rec%xd%d+rec%xd%o-img%xm%o- &
                 (padtraces+tapertr)*img%xm%d)/img%xm%d)+1.5
             source1(ino,:)=source(1,:)
          else
             source=0;source1=0.
             call sreed("sou",source,8*sou%h%n*sou%w%n)
             source1=source
          end if
          call sreed("in",recv,8*rec%h%n*rec%w%n)
          if (exist_file("vref")) then
             call owwemva(.true.,taper,recv,source1,vel,deltaS,B_image,vref)
          else
             call owwemva(.true.,taper,recv,source1,vel,deltaS,B_image)
          end if
          scatS=scatS+deltaS
          call owexec_time(finish)
          write(0,111) "EXTRAPOLATED WAVEFIELDS = ",idip," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do
       if (icwind1) then
          icwind=.true.
       end if

       call owexec_time(finish1)
       write(0,112) "COMPUTED THE GRADIENT IN ",(finish1-start1)/60," min OF ITERATION # ", iter

       if (exist_file("grad")) then
          call srite("grad",scatS(:,tapertr+padtraces+1:img%xm%n-tapertr-padtraces), &
                    4*img%z%n*(img%xm%n-2*(tapertr+padtraces)))
       end if


! Formating the gradient and
! computing the objective function
       f=dble(0.)
       do i=1,img%xm%n
          do k=1,img%z%n
             g((i-1)*img%z%n+k)=(scatS(k,i))
             do j=1,img%sh%n
                f = f + (B_image(k,j,i)*B_image(k,j,i))
             end do
          end do
       end do

       write(2,*) "Sum gradient = ",sum(abs(g)) , "   Obj.function = ", f
 
    end if

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

 end do
 111 format(a26,i4,a8,i4,a4,f6.2,a4)
 112 format(a29,f10.2,a20,i4)
 close(2)
 close(3)
 call srite("out",scatS(:,tapertr+padtraces+1:img%xm%n-tapertr+padtraces), &
           4*img%z%n*(img%xm%n-2*(tapertr+padtraces)))

 deallocate(x,l,u,g,nbd)
 deallocate(wa)
 deallocate(iwa)
 deallocate(source)
 deallocate(recv)
 deallocate(deltaS)
 deallocate(scatS)
 deallocate(deltaI)
 deallocate(R_image)
 deallocate(B_image)
 deallocate(taper)
 deallocate(vel)
 deallocate(vref)
 if (resmig) then
    deallocate(tfld3)
    call owfftdestroy3()
 end if

 call sep_end_prog()

end program OWlbfgsc
