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
!use cgstep_mod
!use base_spline_real
!use solver_reg_mod

implicit none
integer, public, save                  :: iter,realz
real, private, allocatable             :: taper(:,:),pick1(:),vel(:),vel1(:)
real, private, allocatable             :: image(:,:,:,:,:),image0(:,:,:,:,:),g_ir(:),grad0(:)
real, private, allocatable             :: img1(:,:,:,:),img0(:,:,:,:),img_ir(:,:,:,:,:)
real, private                          :: minvel
double precision                       :: aa,bb
double precision, private, allocatable :: gspl(:)
!real, private, allocatable             :: gspl(:)
contains 
!---------------------------------------------------------------------
subroutine steep_mpi_init(pick_in,minvel_in)
integer                  :: i, i_r
real                     :: pick_in(:),minvel_in

    allocate(pick1(img%xm%n*img%ym%n))
    pick1=pick_in
    minvel=minvel_in
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Set up the taper
    allocate(taper(kxm%n,kym%n));taper=0.
    call owtaper(taper)
    allocate(image0(kxm%n,nsecdimx,kym%n,nsecdimy,img%z%n));image0=0
    allocate(image(kxm%n,nsecdimx,kym%n,nsecdimy,img%z%n));image=0
    allocate(g_ir(kxm%n*nz*kym%n));g_ir=0.
    allocate(grad0(kxm%n*nz*kym%n));grad0=0
    allocate(vel(kxm%n*nz*kym%n));vel=0
    allocate(vel1(kxm%n*nz*kym%n));vel1=0
    if (axk%n.ne.99999 .and. ayk%n.ne.99999 .and. azk%n.ne.99999) allocate(gspl(axk%n*ayk%n*azk%n))

end subroutine
!---------------------------------------------------------------------
subroutine steep_mpi_close()
  if (allocated(taper)) deallocate(taper)
  if (allocated(pick1)) deallocate(pick1)
  if (allocated(image)) deallocate(image)
  if (allocated(image0)) deallocate(image0)
  if (allocated(grad0)) deallocate(grad0)
  if (allocated(g_ir)) deallocate(g_ir)
  if (allocated(gspl)) deallocate(gspl)
  if (allocated(vel)) deallocate(vel)
  if (allocated(vel1)) deallocate(vel1)
end subroutine
!---------------------------------------------------------------------
subroutine fgeval_init(iter_in,vel1_in)
 integer                     :: iter_in
 real,optional               :: vel1_in(:)
 iter=iter_in

 if (present(vel1_in)) vel1=vel1_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(:)

allocate(a(kxm%n*kym%n*nsecdimx*nsecdimy)); a=0
call sseek("bkgimg",0,0)
call sseek("Dimage",0,0)
call sseek("grad",0,0)
do i=1,img%z%n
   call srite("bkgimg",a,4*kxm%n*kym%n*nsecdimx*nsecdimy)
   call srite("Dimage",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)
end subroutine
!---------------------------------------------------------------------
subroutine fgeval(vel0,ff,gg,grad,ngrad)

double precision,optional          :: grad(:),ngrad
double precision                   :: ff,gg,funcao
real                               :: vel0(:)
!---
integer                            :: ierr,idip,stat1,i_r,status,izb,iz, ihy,iy,ihx,ix
real                               :: start2,finish2
integer                            :: nx_max,nx_min,ny_max,ny_min

   nx_max=(xv_max-img%xm%o)/img%xm%d+1.5
   nx_min=(xv_min-img%xm%o)/img%xm%d+1.5
   ny_max=(yv_max-img%ym%o)/img%ym%d+1.5
   ny_min=(yv_min-img%ym%o)/img%ym%d+1.5
    image=0.
    do idip=1,shtblk
       call feval_exe(idip,vel0,image)
    end do
    if(wpar%i_mpi == 0) image=0.
    call MPI_BARRIER(MPI_COMM_WORLD,ierr)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    call owexec_time(start2)
    image0=0
    do iz=1,img%z%n
       call MPI_BARRIER(MPI_COMM_WORLD,ierr)
       call MPI_ALLREDUCE(image(:,:,:,:,iz),image0(:,:,:,:,iz),&
                kxm%n*nsecdimx*kym%n*nsecdimy,MPI_REAL,MPI_SUM,MPI_COMM_WORLD,ierr)
    end do

    call owexec_time(finish2)
    if (wpar%i_mpi==1) write(0,*) "FINISHED COPYING bimg ",int((finish2-start2))," s"

    if (scaling) then
       call ow_scale1(iter,image0)
    end if
 
    call srite("bimg",image0(:,:,int((yline_-img%ym%o)/img%ym%d)+1,:,:),4*kxm%n*nsecdimx*nsecdimy*nz)
    call MPI_BARRIER(MPI_COMM_WORLD,ierr)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    call owexec_time(start2)
    ff=0.;gg=0.
    if (mdso==1) then
       write(0,*) "MDSO"
       call owdso_nozb1(pick1,image0,ff)
    else if (mdso==2) then
       call owdso_nozb2(pick1,image0,ff)
    else
       write(0,*) "DSO"
       call owdso_nozb0(pick1,image0,ff)
    end if
    call MPI_BARRIER(MPI_COMM_WORLD,ierr)
    gg = ff*ff

    call srite("dimg",image0(:,:,int((yline_-img%ym%o)/img%ym%d)+1,:,:),4*kxm%n*nsecdimx*nsecdimy*nz)
    call MPI_BARRIER(MPI_COMM_WORLD,ierr)

    if (wpar%i_mpi==1) write(0,*) "@@4 "
    call owexec_time(finish2)
    if (wpar%i_mpi==1) write(0,*) "FINISHED DSO IN ",int((finish2-start2))," s"
    write(0,*) "ff= ",ff," gg= ",gg,wpar%i_mpi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    call MPI_BARRIER(MPI_COMM_WORLD,ierr)
    if (present(grad)) then
       grad0=0;grad=0.
       do idip=1,shtblk
          call geval_exe(idip,vel0,image0,grad0);
       end do
       call MPI_BARRIER(MPI_COMM_WORLD,ierr)
       if(wpar%i_mpi == 0) grad0=0.
       grad0=-grad0
       g_ir=0.
       call MPI_BARRIER(MPI_COMM_WORLD,ierr)
       call MPI_ALLREDUCE(grad0,g_ir,kxm%n*kym%n*nz,MPI_REAL,MPI_SUM,MPI_COMM_WORLD,ierr)
       grad=dble(g_ir)

       if (present(ngrad)) ngrad=norm(grad) 
       if (chop) then
          if (iter==1) then
             aa=(1-chophi)*maxval(grad)
             bb=(1-choplo)*minval(grad)
          end if
          where(grad>aa) grad=aa
          where(grad<bb) grad=bb
       end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Somente teste
       grad=sign(abs(grad)**expgrad,grad)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!       call MPI_BARRIER(MPI_COMM_WORLD,ierr)
       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)
              gspl=dble(0.)
              stat1 = base_spln3d_oper(.true. , .false., gspl, grad)
              grad=dble(0.)
              stat1 = base_spln3d_oper(.false. ,.false., gspl, grad)
              call owexec_time(finish2)
              write(3,*) "SMOOTHED GRADIENT IN ",int((finish2-start2))," s",wpar%i_mpi
       end if

       grad=(grad/maxval(abs(grad)))*minvel
       if (wpar%i_mpi==0) write(0,*) "03 ",minval(grad),maxval(grad)
       call oweditvel(pick1,grad)
       if (wpar%i_mpi==0) write(0,*) "05 ",minval(grad),maxval(grad)

    end if
    call MPI_BARRIER(MPI_COMM_WORLD,ierr)

end subroutine
!---------------------------------------------------------------------
subroutine feval_exe(idip,vel0,image)

real                                          :: start, finish, start1, finish1, vel0(:)
integer                                       :: idip,ireal
real, allocatable                             :: vref(:,:,:)
real                                          :: image(:,:,:,:,:)
integer                                       :: izb,i_r,status,ierr

       if (extrap(1:4).eq."pspi") then
          if (lloyd) then 
            allocate(vref(nref,2,img%z%n))
            call lloyd_vel1(vel0,vref)
!            if (wpar%i_mpi==0) call srite("vreff",vref,4*nref*2*img%z%n)
          end if
       end if
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       ! Computing bkg.image
       call owexec_time(start)
       call owwei_init(idip,1,iter)
       if (lloyd) then
          call owwei_nozb(taper=taper,vel=vel0,image=image,sref=vref)
       else
          call owwei_nozb(taper=taper,vel=vel0,image=image)
       end if
       call owexec_time(finish)
       write(3,112) "COMPUTED BACKGROUND IMAGE IN ",(finish-start)/60," min"
       if (allocated(vref)) deallocate(vref)

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

 real                                          :: start, finish, start1, finish1, start2, finish2,vel0(:)
 integer                                       :: ni,i,ishot,j,k,stat1,idip
 real, allocatable                             :: vref(:,:,:)
 real                                          :: image(:,:,:,:,:),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.;

       call owexec_time(start1)
       call owwei_init(idip,1,iter)
       if (lloyd) then
          call owwemva(taper,vel0,iter,image,g,vref)
       else
          call owwemva(taper,vel0,iter,image,g)
       end if

       call owexec_time(finish1)
       if (wpar%i_mpi==1) write(3,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
!---------------------------------------------------------------------

double precision function norm(x)
  double precision                  :: x(img%xm%n,img%ym%n,img%z%n)
  integer                           :: i,j,k
  norm=0.

  do i=1,img%z%n
     do j=1,img%ym%n
        do k=1,img%xm%n
           norm=norm+x(k,j,i)**2
        end do
     end do
  end do
  norm=sqrt(norm)
end function
!---------------------------------------------------------------------
end module
