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 OW_mpi

implicit none
integer, public, save                  :: iter,realz
real, private, allocatable             :: taper(:,:),pick1(:)

contains 
!---------------------------------------------------------------------
subroutine steep_mpi_init(pick_in)
integer                  :: i, i_r
real                     :: pick_in(:)

    allocate(pick1(img%xm%n*img%ym%n))
    pick1=pick_in
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Set up the taper
    allocate(taper(kxm%n,kym%n));taper=0.
    call owtaper(taper)

end subroutine
!---------------------------------------------------------------------
subroutine steep_mpi_close()
  if (allocated(taper)) deallocate(taper)
  if (allocated(pick1)) deallocate(pick1)
end subroutine
!---------------------------------------------------------------------
subroutine fgeval_init(iter_in)
 integer                     :: iter_in
 iter=iter_in

 if (nreal.ne.1) then
    realz=mod(iter,nreal)
 else
    realz=1
 end if

end subroutine

!---------------------------------------------------------------------
subroutine init_files()
integer            :: i,j,k
real,allocatable   :: a(:)

write(0,*) "init_files 00"
allocate(a(kxm%n*kym%n*nsecdimx*nsecdimy)); a=0
do i=1,img%z%n
   call srite("bkgimg",a,4*kxm%n*kym%n*nsecdimx*nsecdimy)
   call srite("grad_",a(1:kxm%n*kym%n),4*kxm%n*kym%n)
end do
deallocate(a)
write(0,*) "init_files 01"
end subroutine
!---------------------------------------------------------------------
subroutine fgeval(vel0,ff,gg,grad,riteD,wsave)

logical                            :: riteD,wsave
double precision,optional          :: grad(:)
double precision                   :: ff,gg
real                               :: vel0(:)
!---
double precision, allocatable      :: g(:),g_r(:),gspl(:)
integer                            :: ierr,idip,stat1,i_r,status
real                               :: start2,finish2

write(0,*) "fgeval 00"
    do idip=1,shtblk
       call feval_exe(idip,vel0,riteD,wsave)
    end do

write(0,*) "fgeval 01"
    ff=dble(0.);gg=ff
    call owdso(ff,riteD,pick1)
    gg = ff*ff
    if (wpar%i_mpi==0) write(0,*) "ff= ",ff," gg= ",gg

    if (present(grad)) then
       allocate(g(size(grad)));g=0
       allocate(g_r(size(grad)));g_r=0;
       do idip=1,shtblk
          g=0;g_r=0;
          call geval_exe(idip,vel0,g);
          call MPI_ALLREDUCE_(g, g_r, img%xm%n*img%ym%n*nz, MPI_DOUBLE_PRECISION, &
                          MPI_SUM,MPI_COMM_WORLD,ierr)
          grad = grad + g_r
       end do
       deallocate(g,g_r)
       if (wpar%i_mpi==0) write(0,*) "main => ",maxval(grad),minval(grad)
       if (axk%n.ne.99999 .and. ayk%n.ne.99999 .and. azk%n.ne.99999) then
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          ! B-spline representation of grad
          stat1 = base_spln3d_init(img%xm,img%ym,img%z,axk,ayk,azk)
          call owexec_time(start2)
          allocate(gspl(axk%n*ayk%n*azk%n))
          gspl=dble(0.)
          stat1 = base_spln3d_oper(.true. , .false., gspl, grad)
          grad=dble(0.)
          stat1 = base_spln3d_oper(.false. ,.false., gspl, grad)
          deallocate(gspl)
          call owexec_time(finish2)
         end if
         if (wpar%i_mpi==0) write(0,*) "main => ",maxval(grad),minval(grad)
         if (wpar%i_mpi==0) write(0,*) "SMOOTHED GRADIENT IN ",int((finish2-start2))," s"
         call oweditvel(pick1,grad)
         call to_history("n5",iter,"grad")
         call srite("grad",sngl(grad),4*img%z%n*img%xm%n*img%ym%n)
       end if

end subroutine
!---------------------------------------------------------------------
!subroutine feval_exe(idip,vel0,fun,riteD,wsave)
subroutine feval_exe(idip,vel0,riteD,wsave)

real                                          :: start, finish, start1, finish1, vel0(:)
integer                                       :: idip,ireal
real, allocatable                             :: vref(:,:,:),dcig(:,:,:,:,:)
double precision                              :: fun
logical                                       :: riteD,wsave
integer                                       :: izb,i_r,status,ierr

write(0,*) "feval_exe 00"
       if (wpar%i_mpi > 0) then
          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
       end if
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       ! Computing bkg.image
       call owexec_time(start)
       call owwei_init(idip,1)
write(0,*) "feval_exe 01"
       if (lloyd) then
          call owwei(taper=taper,vel=vel0,sref=vref,wsave=wsave)
       else
          call owwei(taper=taper,vel=vel0,wsave=wsave)
       end if
write(0,*) "feval_exe 02"

       call owexec_time(finish)
       write(3,112) "COMPUTED BACKGROUND IMAGE IN ",(finish-start)/60," min"
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       ! Computing DSO

!       fun=dble(0.)
!       call owdso(fun,riteD,pick1)
       if (wpar%i_mpi > 0) then
          if (allocated(vref)) deallocate(vref)
       end if

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

 real                                          :: start, finish, start1, finish1, start2, finish2,vel0(:)
 integer                                       :: ni,i,ishot,j,k,stat1,idip
 real,allocatable                              :: deltaS(:),scatS(:)
 real, allocatable                             :: vref(:,:,:)
 double precision                              :: g(:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Computing the gradient
!
       allocate(deltaS(img%z%n*kxm%n*kym%n),scatS(img%z%n*img%xm%n*img%ym%n))
       if (wpar%i_mpi > 0) then
          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
       end if
       start1=0.;finish1=0.;g=dble(0.)

       call owexec_time(start1)
       start1=0.;finish1=0.;deltaS=0.
       call owwei_init(idip,1)
       if (lloyd) then
          call owwemva(deltaS,taper,vel0,vref)
       else
          call owwemva(deltaS,taper,vel0)
       end if
       g=dble(deltaS)

       call owexec_time(finish1)
       if (wpar%i_mpi==0) write(0,112) "COMPUTED THE GRADIENT IN ",(finish1-start1)/60," min OF ITERATION # ", iter
       if (wpar%i_mpi > 0) then
          if (allocated(vref)) deallocate(vref)
       end if
 111 format(a26,i4,a8,i4,a4,f6.2,a4)
 112 format(a29,f10.2,a20,i4)
       deallocate(deltaS,scatS)
 end subroutine
!---------------------------------------------------------------------

real function norm(x)
  double precision                  :: x(:)
  integer                           :: i,j,k
  norm=0.

  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
              norm=norm+x((k-1)*img%xm%n+i)**2
           end if
        end do
     end do
  end do
  norm=sqrt(norm)
end function
!---------------------------------------------------------------------
end module
