module steepest_mod

use sep
use OW_grabvel2_areal
use OW_exec_time
use OW_parms_areal
use OW_taper_areal
use OW_wei_areal
use lloyd_vel
use base_spline_mod
!use triangle2_mod
!use cgstep_mod0
use solver_tiny_mod

implicit none
real, private, allocatable             :: vel0(:),taper(:,:),pick1(:)
integer, private                       :: iter, rect1, rect2,rect3,realz,i_node,n_node,ierr

contains 

!---------------------------------------------------------------------
subroutine fgeval_init(vel0_in,taper_in,iter_in,pick1_in)
 real                        :: vel0_in(:),taper_in(:,:)
 real,optional               :: pick1_in(:)
 integer                     :: iter_in
 iter=iter_in
 allocate(vel0(img%xm%n*img%ym%n*nz),taper(kxm%n,kym%n))
 vel0=vel0_in
 taper=taper_in
 if (present(pick1_in)) then
    allocate(pick1(img%xm%n*img%ym%n))
    pick1=pick1_in
 end if
 if (nreal.ne.1) then
    realz=mod(iter,nreal)
 else
    realz=1
 end if

end subroutine
!---------------------------------------------------------------------
subroutine fgeval_close()
 if (allocated(vel0))   deallocate(vel0)
 if (allocated(taper)) deallocate(taper)
 if (allocated(pick1)) deallocate(pick1)
end subroutine
!---------------------------------------------------------------------
subroutine init_files()
call auxclose("Dimage"); call auxinout("Dimage"); 
call auxclose("bkgimg"); call auxinout("bkgimg"); 
call auxclose("bkws");   call auxinout("bkws"); 
call auxclose("bkwr");   call auxinout("bkwr"); 
end subroutine
!---------------------------------------------------------------------
subroutine feval(fun,riteD)

real                                          :: start, finish, start1, finish1
integer                                       :: idip,ireal
real, allocatable                             :: vref(:,:,:)
double precision                              :: fun
logical                                       :: riteD

       start=0.;finish=0.
       call sseek("bkgimg",0,0)
       call sseek("Dimage",0,0)
       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
       start1=0.;finish1=0.;
       call owexec_time(start1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Computing bkg.image
       call sseek("rec",0,0)
       call sseek("sou",0,0)
       do ireal=1,nreal
       do idip=1,rec%xd%n
          call owexec_time(start)
          if (lloyd) then
             call owwei_init(idip,ireal)
             call owwei(taper,vel0,vref)
          else
             call owwei_init(idip,ireal)
             call owwei(taper,vel0)
          end if
          call owexec_time(finish)
          write(3,111) "Shot = ",idip," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do
       end do
       call owexec_time(finish1)
       write(3,112) "COMPUTED BACKGROUND IMAGE IN ",(finish1-start1)/60," min"
       if (allocated(vref)) deallocate(vref)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Computing DSO
       fun=dble(0.)
       if (allocated(pick1)) then
          call owdso(fun,riteD,pick1)
       else
          call owdso(fun,riteD)
       end if

 111 format(a7,i4,a8,i4,a4,f6.2,a4)
 112 format(a29,f10.2,a4)
end subroutine
!---------------------------------------------------------------------
 subroutine geval(g)

 real                                          :: start, finish, start1, finish1, start2, finish2
 integer                                       :: ni,i,ishot,j,k,stat1
 real                                          :: deltaS(img%z%n*kxm%n*kym%n),scatS(img%z%n*img%xm%n*img%ym%n)
 real, allocatable                             :: vref(:,:,:)
 double precision                              :: g(:),gspl(img%z%n*img%xm%n*img%ym%n)

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

       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
       start1=0.;finish1=0.;g=dble(0.)

       call sseek("rec",0,0)
       call sseek("sou",0,0)
       call owexec_time(start1)
       start1=0.;finish1=0.;scatS=0.
       call owexec_time(start1)
       do ishot=1,rec%xd%n
          start=0.;finish=0.;deltaS=0.
          call owexec_time(start)
          call owwemva_init(ishot)
          if (lloyd) then
             call owwemva(deltaS,taper,vel0,vref)
          else
             call owwemva(deltaS,taper,vel0)
          end if
          scatS=scatS+deltaS
          call owexec_time(finish)
          write(3,111) "EXTRAPOLATED WAVEFIELDS = ",ishot," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do
       g=dble(-scatS)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Smoothing the gradient
       write(0,*) "main => ",maxval(g),minval(g)

!       if (exist_file("horv0")) call oweditvel(pick1,g)

       if (axkn.ne.99999 .or. axkn.ne.99999 .or. azkn.ne.99999) then
          stat1 = base_spln3d_init(img%xm,img%ym,img%z,axk,ayk,azk)
          call owexec_time(start2)
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          ! B-spline representation of grad
          gspl=dble(0.)
          stat1 = base_spln2d_oper(.true. , .false., gspl, g)
          g=dble(0.)
          stat1 = base_spln2d_oper(.false. ,.false., gspl, g)
          call owexec_time(finish2)
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
       end if
       write(0,*) "main => ",maxval(g),minval(g)
       write(0,112) "SMOOTHED GRADIENT IN ",(finish2-start2)/60
       if (exist_file("horv0")) call oweditvel(pick1,g)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       call owexec_time(finish1)
       write(0,112) "COMPUTED THE GRADIENT IN ",(finish1-start1)/60," min OF ITERATION # ", iter
       call to_history("n4",iter,"grad")
       call srite("grad",sngl(g),4*img%z%n*img%xm%n*img%ym%n)

       if (allocated(vref)) deallocate(vref)
 111 format(a26,i4,a8,i4,a4,f6.2,a4)
 112 format(a29,f10.2,a20,i4)
 end subroutine
!---------------------------------------------------------------------

!---------------------------------------------------------------------
! subroutine ow_smoovel(x,vel0,pick1)
!
! real, optional             :: pick1(:),vel0(:)
! double precision           :: x(:),ex,exs
!! double precision           :: gspl(axk%n*azk%n)
! real                       :: start2,finish2,slow(img%z%n*img%xm%n)
! real, allocatable          :: gspl(:)
! integer                    :: stat1,i,j,ix,iz
! 
!       call triangle2_init(rect1,rect3,img%xm%n,img%z%n)
!       if (present(pick1)) call oweditvel(pick1,x,vel0)
!       call owexec_time(start2)
!       stat1 = triangle2_lop(.true.,.false.,slow,sngl(x))
!       x=dble(slow)
!       call owexec_time(finish2)
!       if (present(pick1)) call oweditvel(pick1,x,vel0)
!       call triangle2_close()
!
!end subroutine
!---------------------------------------------------------------------
 subroutine owtri_init(rect1_in,rect2_in,rect3_in)
 integer           :: rect1_in,rect2_in,rect3_in
   rect1=rect1_in
   rect2=rect2_in
   rect3=rect3_in
 end subroutine
!---------------------------------------------------------------------
real function norm(x,pick1)
  double precision                  :: x(:)
  real,optional                     :: pick1(:)
  integer                           :: i,j,k
  norm=0.
  if (present(pick1)) then
     do i=1,img%xm%n
        do j=1,img%xm%n
           do k=1,img%z%n
              if ((k-1)*img%z%d+img%z%o > pick1(i)) then
!!!!!!!!!!!!!!!!
!
!resolver
                 norm=norm+x((k-1)*img%xm%n+i)**2
              end if
           end do
        end do
     end do
     norm=sqrt(norm)
  else
     norm = sqrt(sum(x*x))
  end if
end function
!---------------------------------------------------------------------
end module
