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,optional            :: pick_in(:)

 if(wpar%i_mpi>0) then
    if (present(pick_in)) then
       allocate(pick1(img%xm%n*img%ym%n))
       pick1=pick_in
    end if
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Set up the taper
    allocate(taper(kxm%n,kym%n));taper=0.
    call owtaper(taper)
 end if
!write(0,*) wpar%i_mpi,"after pick,taper"

end subroutine
!---------------------------------------------------------------------
subroutine steep_mpi_close()
  deallocate(taper)
  if (allocated(pick1)) deallocate(taper,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

do i=1,kxm%n*kym%n*img%z%n
!   call srite("Velout",0.,4)
!   call srite("grad",0.,4)
   do j=1,nsecdimx*nsecdimy
      call srite("Dimage",0.,4)
      call srite("bkgimg",0.,4)
   end do
end do
end subroutine
!---------------------------------------------------------------------
subroutine fgeval(vel0,ff,gg,grad,riteD)

logical                            :: riteD
double precision,optional          :: grad(:)
double precision                   :: ff,gg
real                               :: vel0(:)
!---
double precision                   :: fn,gn,fn_r
double precision, allocatable      :: g(:),gspl(:)
complex                            :: ws(kxm%n,kym%n,sou%w%n),wr(kxm%n,kym%n,sou%w%n)
integer                            :: ierr,idip,stat1
real                               :: start2,finish2

    if (present(grad)) then
       allocate(g(size(grad)))
       g=0
    end if
    if (wpar%i_mpi==0) then
       call sseek("rec",0,0)
       call sseek("sou",0,0)
       fn=0;gn=0
    else
       fn=0
    end if
    do idip=1,shtblk
       call feval_exe(idip,vel0,fn,riteD)
       if (present(grad)) call geval_exe(idip,vel0,g)
       if (wpar%i_mpi > 0) then
          call MPI_REDUCE_(fn, fn_r, 1, MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)
          if (present(grad)) then
             call MPI_REDUCE_(g, g, img%xm%n*img%ym%n*nz, MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)
          end if
       else
          if (present(grad)) grad = grad + g
          ff = ff + fn_r
          gg = gg + ff*ff
       end if
    end do

    if (present(grad)) then
     if (wpar%i_mpi==0) then
       write(0,*) "main => ",maxval(grad),minval(grad)
       if (axkn.ne.99999 .or. axkn.ne.99999 .or. azkn.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(axkn*aykn*azkn))
          gspl=dble(0.)
          stat1 = base_spln2d_oper(.true. , .false., gspl, grad)
          grad=dble(0.)
          stat1 = base_spln2d_oper(.false. ,.false., gspl, grad)
          deallocate(gspl)
          call owexec_time(finish2)
       end if
       write(0,*) "main => ",maxval(grad),minval(grad)
       write(0,*) "SMOOTHED GRADIENT IN ",int((finish2-start2))," s"
       if (allocated(pick1)) call oweditvel(pick1,grad)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       call to_history("n4",iter,"grad")
       call srite("grad",sngl(grad),4*img%z%n*img%xm%n*img%ym%n)
       deallocate(g)
     end if
    end if

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

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

       call sseek("bkgimg",0,0)
       call sseek("Dimage",0,0)
       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,*) "before owwei ", wpar%i_mpi
       if (lloyd) then
          call owwei(taper,vel0,vref)
       else
          call owwei(taper,vel0)
       end if

       if (wpar%i_mpi==0) then
          call owexec_time(finish)
          call owexec_time(finish)
          write(3,112) "COMPUTED BACKGROUND IMAGE IN ",(finish-start)/60," min"
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          ! Computing DSO
          fun=dble(0.)
          if (allocated(pick1)) then
             call owdso(fun,riteD,pick1)
          else
             call owdso(fun,riteD)
          end if
       else
          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                                          :: deltaS(img%z%n*kxm%n*kym%n),scatS(img%z%n*img%xm%n*img%ym%n)
 real, allocatable                             :: vref(:,:,:)
 double precision                              :: g(:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! 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 owwei_init(ishot,1)
          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)

       write(0,*) "main => ",maxval(g),minval(g)
       call owexec_time(finish1)
       write(0,112) "COMPUTED THE GRADIENT IN ",(finish1-start1)/60," min OF ITERATION # ", iter
       if (allocated(vref)) deallocate(vref)
 111 format(a26,i4,a8,i4,a4,f6.2,a4)
 112 format(a29,f10.2,a20,i4)
 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
