module steepest_mod

use sep
use mpi_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
real, private, allocatable             :: vel(:),taper(:,:),pick1(:)
integer, private                       :: iter, rect1, rect2,rect3,realz,i_node,n_node,ierr

contains 

!---------------------------------------------------------------------
subroutine fgeval_init(vel_in,taper_in,iter_in,pick1_in)
 real                        :: vel_in(:),taper_in(:,:)
 real,optional               :: pick1_in(:)
 integer                     :: iter_in
 iter=iter_in
 allocate(vel(img%xm%n*img%ym%n*nz),taper(kxm%n,kym%n))
 vel=vel_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(vel))   deallocate(vel)
 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,i_r,ierr,i_node,n_node
real, allocatable                             :: vref(:,:,:)
double precision                              :: fun
logical                                       :: riteD
integer                                       :: ninx,nfnx,niny,nfny,mpierr,rc,shtblk,shtres,shtidx
integer                                       :: izb, initz,finz,status
complex                                       :: ws(kxm%n,kym%n,sou%w%n),wr(kxm%n,kym%n,sou%w%n)
real                                          :: bimg(kxm%n,nsecdimx,kym%n,nsecdimy,zblock%ns)
real                                          :: image(kxm%n,nsecdimx,kym%n,nsecdimy,zblock%ns)

       call sseek("rec",0,0)
       call sseek("sou",0,0)
       start=0.;finish=0.
       call sseek("bkgimg",0,0)
       call sseek("Dimage",0,0)

       ninx=tapertrx+padtracesx+1
       nfnx=rec%hx%n+padtracesx+tapertrx
       niny=tapertry+padtracesy+1
       nfny=rec%hy%n+padtracesy+tapertry

       call MPI_INIT(mpierr)
       if (mpierr .ne. MPI_SUCCESS) then
          write(0,*) 'Error starting MPI program. Terminating.'
          call MPI_ABORT(MPI_COMM_WORLD, rc, mpierr)
       end if

       call MPI_COMM_RANK(MPI_COMM_WORLD, i_node, ierr)
       call MPI_COMM_SIZE(MPI_COMM_WORLD, n_node, ierr)
       call owparm_bcast()

       if (i_node==0) then
         do i_r=1,n_node-1
           call MPI_Send(taper,kxm%n*kym%n,         MPI_REAL,i_r,100,MPI_COMM_WORLD,ierr)
           call MPI_Send(vel, img%xm%n*img%ym%n*nz,MPI_REAL,i_r,110,MPI_COMM_WORLD,ierr)
           if (allocated(pick1)) call MPI_Send(pick1, img%xm%n*img%ym%n,MPI_REAL,i_r,120,MPI_COMM_WORLD,ierr)
         end do
       else
         call MPI_Recv(taper,kxm%n*kym%n,         MPI_REAL,0,100,MPI_COMM_WORLD,status,ierr)
         call MPI_Recv(vel, img%xm%n*img%ym%n*nz,MPI_REAL,0,110,MPI_COMM_WORLD,status,ierr)
         if (allocated(pick1)) call MPI_Recv(pick1, img%xm%n*img%ym%n,MPI_REAL,i_r,120,MPI_COMM_WORLD,ierr)
       end if

       if (i_node==0) then
          if (lloyd) then 
             allocate(vref(nref,2,img%z%n))
             call lloyd_vel1(vel,vref)
             do i_r=1,n_node-1
                call MPI_Send(vref,nref*2*img%z%n, MPI_REAL,i_r,140,MPI_COMM_WORLD,ierr)
             end do
          end if
       else
          if (lloyd) then 
             call MPI_Recv(vref,nref*2*img%z%n, MPI_REAL,0,140,MPI_COMM_WORLD,status,ierr)
          end if
       end if
       start1=0.;finish1=0.;
       call owexec_time(start1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Computing bkg.image
       shtblk=int(rec%xd%n/(n_node-1))
       shtres=mod(rec%xd%n,(n_node-1))
       do idip=1,shtblk
          shtidx=n_node-1
          if (idip==shtblk) shtidx=shtres-1
          if (i_node==0) then
             do i_r=1,shtidx
                wr=0;   call gathers_io((idip-1)*(n_node-1)+i_r,wr(ninx:nfnx,niny:nfny,:),"rec","reed")
                ws=0;   call gathers_io((idip-1)*(n_node-1)+i_r,ws(ninx:nfnx,niny:nfny,:),"sou","reed")
                call MPI_Send(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,150,MPI_COMM_WORLD,ierr)
                call MPI_Send(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,160,MPI_COMM_WORLD,ierr)
             end do
          else
             call MPI_Recv(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,150,MPI_COMM_WORLD,status,ierr)
             call MPI_Recv(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,160,MPI_COMM_WORLD,status,ierr)
          end if
          call owwei_init()
          do izb=1,zblock%nb
             if (izb==1) then
                initz=2
                finz=0
             else
                initz=1
                finz=(izb-1)*zblock%blen(izb-1)
             end if
             if (lloyd) then
                call owwei(taper,vel,image,vref)
             else
                call owwei(taper,vel,image)
             end if
             call MPI_Reduce(image, bimg, kxm%n*nsecdimx*kym%n*nsecdimy*zblock%blen(izb), MPI_REAL,MPI_SUM,0,MPI_COMM_WORLD,ierr)
             if (i_node==0) then
                call image_io("bkgimg","reed",izb,finz,bimg)
                call image_io("bkgimg","rite",izb,finz,bimg)
             end if
          end do
          call owexec_time(finish)
          write(3,111) "Shot = ",idip," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do
       call owexec_time(finish1)
       !----------------------------------------------------------------  
       !! shutdown MPI
       call MPI_FINALIZE(mpierr) 

       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(vel,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,vel,vref)
          else
             call owwemva(deltaS,taper,vel)
          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,vel,pick1)
!
! real, optional             :: pick1(:),vel(:)
! 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,vel)
!       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,vel)
!       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
!---------------------------------------------------------------------
subroutine gathers_io(idip,GG,tag,op)
complex              :: GG(:,:,:)
character*3          :: tag
character*4          :: op
integer              :: idip,ireal

if(op.eq."reed") then 
  call sreed_window(tag,5,(/rec%hx%n,rec%hy%n,sou%w%n,rec%xd%n/),&
                          (/rec%hx%n,rec%hy%n,sou%w%n,       1/),&
                          (/       0,       0,      0,  idip-1/),&
                          (/       1,       1,      1,       1/),&
                          8,GG) 
else
  call srite_window(tag,5,(/rec%hx%n,rec%hy%n,sou%w%n,rec%xd%n/),&
                          (/rec%hx%n,rec%hy%n,sou%w%n,       1/),&
                          (/       0,       0,      0,  idip-1/),&
                          (/       1,       1,      1,       1/),&
                          8,GG) 
end if
end subroutine
!---------------------------------------------------------------------
end module
