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

implicit none
integer, private                       :: iter,realz

contains 
!---------------------------------------------------------------------
subroutine steep_mpi_init()


 call init_files()
 allocate(taper(kxm%n,kxm%n));taper=0.
 if (wpar%i_mpi==0) then
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Set up the taper
    call owtaper(taper)
    do i_r=1,wpar%n_mpi-1
       call MPI_Send_(taper,   kxm%n*kym%n,         MPI_REAL,i_r,100,MPI_COMM_WORLD,ierr)
    end do
 else
    call MPI_Recv_(taper,   kxm%n*kym%n,         MPI_REAL,0,100,MPI_COMM_WORLD,status,ierr)
 end if

 shtblk=int(rec%xd%n/(wpar%n_mpi-1))
 shtres=mod(rec%xd%n,(wpar%n_mpi-1))
 shtblm = wpar%n_mpi-1



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()
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 fgeval(inode,f,gn,grad,riteD)


    if (wpar%i_mpi==0) then
       call sseek("rec",0,0)
       call sseek("sou",0,0)
       fn=0;gn=0;g_r=0;fn_r=0
       if
    else
       fn=0;g=0;
    end if
    do idip=1,shtblk ! Reading in/transfering data
       if (shtblk - idip +1 <= shtres) then
          shtidx = wpar%n_mpi
       else
          shtidx = wpar%n_mpi-1
       end if
       if (wpar%i_mpi==0) then
          do i_r=1,shtidx
             wr=0;   call gathers_io((idip-1)*shtblm+i_r,wr(ninx:nfnx,niny:nfny,:),"rec","reed")
             ws=0;   call gathers_io((idip-1)*shtblm+i_r,ws(ninx:nfnx,niny:nfny,:),"sou","reed")
             call MPI_Send_(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,120,MPI_COMM_WORLD,ierr)
             call MPI_Send_(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,130,MPI_COMM_WORLD,ierr)
          end do
          shtblm=shtidx
       else
          call MPI_Recv_(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,120,MPI_COMM_WORLD,status,ierr)
          call MPI_Recv_(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,130,MPI_COMM_WORLD,status,ierr)
       end if
       if (wpar%i_mpi > 0) then
          call fgeval_init(iters)
          call feval_exe(ws,wr,sngl(x),taper,fn,.true.)
          if (present(grad)) call geval(ws,wr,sngl(x),taper,g)
          
!          call MPI_Send_(fn,                   1, MPI_DOUBLE_PRECISION,0,140,MPI_COMM_WORLD,ierr)
!          call MPI_Send_(g, img%xm%n*img%ym%n*nz, MPI_DOUBLE_PRECISION,0,150,MPI_COMM_WORLD,ierr)
          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_r, 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_r
          f  = f + fn_r
          gn = gn + fn*fn
          call MPI_BCAST_(gn,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpierr) 
       end if
    end do

    if (present(grad)) then
     if (wpar%i_mpi==0) then
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       ! Smoothing the gradient
       write(0,*) "main => ",maxval(g),minval(g)
       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)
          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 to_history("n4",iter,"grad")
       call srite("grad",sngl(g),4*img%z%n*img%xm%n*img%ym%n)



end subroutine
!---------------------------------------------------------------------
subroutine feval_exe(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)

       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
