program OWcgwemva

use sep 
use OW_type
use OW_grabvel2_areal
use OW_exec_time
use OW_parms_areal
use OW_taper_areal
use OW_weicg_areal
use OW_dotst1
!!!!!!!!!
! solver and regularization
use solver_prc_mod
use solver_reg_mod
use solver_tiny_mod
use laplac2
use cgstep_mod
use igrad1
use igrad2
use OW_wghtmute
! vel
use lloyd_vel
! using dips to smooth gradient
!use killwave

implicit none

character*80                                  :: hor1,log,listagem
real                                          :: start, finish, start1, finish1, start2, finish2
real                                          :: a1,b1,vmin,vmax,stat1,v0,r
integer                                       :: ixs,ino,ni,i,st,ishot,t3,iz,j,k,iter,ix,ish,titer,i1,i2,i3,nl_iter
real, allocatable                             :: velint(:),taper(:),pick1(:)
real, allocatable                             :: deltaI(:,:,:),scatS(:),deltaS(:),illum(:,:),illumination(:,:)
real, allocatable                             :: R_image(:),B_image(:),tmp(:)
real, allocatable                             :: vel(:),vref(:,:,:),vel0(:),indip(:),velnew(:)
complex, allocatable                          :: recv(:,:),source(:,:),source1(:)
logical                                       :: cond,bounds
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 call sep_init()
 call sep_begin_prog()
 start=0.;finish=0.
 call owexec_time(start)
 call owparam()
 call from_param("log",log,"log")
 open(unit=2,file=log)
 call from_param("list",listagem,"listagem")
 open(unit=3,file=listagem)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Memory allocation
!
allocate(taper(kxm%n))
allocate(velint(v%z%n*v%xm%n),vel(img%z%n*kxm%n),vel0(img%z%n*img%xm%n),velnew(img%z%n*img%xm%n))
allocate(scatS(img%z%n*img%xm%n),deltaI(kxm%n,img%sh%n,img%z%n))

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Set up the taper
!
 call owtaper(tapertr,taper)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Read velocities
!
 call sreed("vel",velint,v%z%n*v%xm%n*4)
 call owgrab_velim(velint,vel0,ni)
 deallocate(velint)
 velnew=vel0
!----------------------------------------------
! Horizon with no velocity update
if (exist_file("horv0")) then
   call from_param('horv0',hor1)
   open(51,file=hor1,status='old')
   allocate(pick1(img%xm%n));pick1=9999.
   do i=1,img%xm%n
      read(51,*,end=1) pick1(i)
   end do
1 close(51)
end if

 iter=1
 call from_param('nl_iter',nl_iter)
 if (dot) then
    call owwemva_lop_init(taper,velnew)
    call owdotst_wemva1()
    call owwemva_close()
    if (extrap(1:4).eq."pspi") call owdotst_pspi(velnew)
    if (extrap(1:3).eq."ffd") call owdotst_ffd(velnew)
!    call owdot_screen(velnew)
    call owdotst_image()
 else
    do while(iter <= nl_iter)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Computing bkg.image
!
       start1=0.;finish1=0.;
       if (allocated(B_image)) deallocate(B_image) ;allocate(B_image(img%xm%n*nsecdim*img%z%n));B_image=0.

       call sseek("rec",0,0)
       call sseek("sou",0,0)
       call owexec_time(start1)
       if (not(areal)) then
          if (allocated(source1)) deallocate(source1); allocate(source1(sou%w%n))
          source1=0;
          if (allocated(source)) deallocate(source); allocate(source(rec%h%n,sou%w%n))
          source=0;
          call sreed("sou",source1,8*sou%w%n)
          source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
       end if
       do ishot=1,rec%xd%n
          start=0.;finish=0. 
          call owexec_time(start)
          if (areal) then
             if (allocated(source)) deallocate(source); allocate(source(sou%h%n,sou%w%n))
             source=0;
             call sreed("sou",source,8*sou%h%n*sou%w%n)
             vel=velnew
             ni=1
          else
             call owgrab_velsht(velnew,vel,ni,ishot)
          end if

          if (allocated(recv))   deallocate(recv);   allocate(recv(rec%h%n,rec%w%n))
          recv=0;deltaI=0

          call sreed("rec",recv,8*rec%h%n*rec%w%n)

          if (lloyd) then
             if (allocated(vref)) deallocate(vref);allocate(vref(nref,2,img%z%n))
             call lloyd_vel1(1/vel,vref)
             call owwei_init()
             call owwei(taper,recv,source,vel,deltaI,vref)
          else
             call owwei_init()
             call owwei(taper,recv,source,vel,deltaI)
          end if
          if (areal) then
             do k=1,img%z%n;do ish=1,nsecdim;do j=1,kxm%n;
                B_image((k-1)*img%xm%n*img%sh%n+(ish-1)*img%xm%n+j)= &
                              B_image((k-1)*img%xm%n*img%sh%n+(ish-1)*img%xm%n+j)+deltaI(j,ish,k)
             end do;end do;end do
          else
             do k=1,img%z%n;do ish=1,nsecdim;do j=ni,ni+kxm%n-1;
                B_image((k-1)*img%xm%n*img%sh%n+(ish-1)*img%xm%n+j)= &
                              B_image((k-1)*img%xm%n*img%sh%n+(ish-1)*img%xm%n+j)+deltaI(j-ni+1,ish,k)
             end do;end do;end do
          end if

          call owexec_time(finish)
          write(0,111) "EXTRAPOLATED WAVEFIELDS = ",ishot," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do

       call owexec_time(finish1)
       write(0,112) "COMPUTED BACKGROUND IMAGE IN ",(finish1-start1)/60," min OF ITERATION # ", iter
       call to_history("n4",iter,"bkgimg")
       call srite("bkgimg",B_image,4*img%z%n*nsecdim*img%xm%n)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Computing DSO
!
       call owdso(.false. , B_image)
       if (exist_file("horv0")) then
          allocate(R_image(img%xm%n*nsecdim*img%z%n));R_image=0.
          call wghtmute_lop_init(pick1)
          call wghtmute_lop3(.false.,.false.,B_image,R_image )
          call wghtmute_close()
          B_image=R_image
          deallocate(R_image)
       end if
       if (l1norm) then
          allocate(R_image(img%xm%n*nsecdim*img%z%n));R_image=0.
          call owsmooth(B_image,R_image)
          B_image=B_image*abs(B_image)/(R_image**R_image)
          if (allocated(R_image)) deallocate(R_image)
       end if
       call owdso(.true.  , B_image)
       call to_history("n4",iter,"Dimage")
       call srite("Dimage",B_image,4*img%z%n*nsecdim*img%xm%n)
       write(0,112) "COMPUTED PERTURBED  IMAGE IN ",(finish1-start1)/60," min OF ITERATION # ", iter

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Slowness update
!
       scatS=0.
!       call laplac2_init(img%xm%n,img%z%n)
!       call igrad2_init(img%xm%n,img%z%n)
       call owwemva_lop_init(taper,velnew)
       if (exist_file("horv0")) then
          call wghtmute_lop_init(pick1)
          if (cgsolv(1:4).eq."tiny") then
             call solver_tiny( m=scatS,d=B_image, Fop=owwemva_lop, stepper=cgstep, niter=niter)
          else if (cgsolv(1:3).eq."reg") then
             call solver_reg( m=scatS,d=B_image, Fop=owwemva_lop, Aop=identity_lop, stepper=cgstep, nAop=img%xm%n*img%z%n, &
                           Jop=wghtmute_lop,niter=niter,eps1=epscg,verb=.true.)!,rm0,err,resd,resm,mmov,rmov,verb)
          else if (cgsolv(1:4).eq."prec") then
              call solver_prc( m=scatS,d=B_image, Fop=owwemva_lop, Sop=laplac2_lop, stepper=cgstep, nSop=img%xm%n*img%z%n, &
                           Jop=wghtmute_lop,niter=niter,eps1=epscg,verb=.true.)!,rm0,err,resd,resm,mmov,rmov,verb)
          end if
          call wghtmute_close()
       else
          if (cgsolv(1:4).eq."tiny") then
             call solver_tiny( m=scatS,d=B_image, Fop=owwemva_lop, stepper=cgstep, niter=niter)
          else if (cgsolv(1:3).eq."reg") then
              call solver_reg( m=scatS,d=B_image, Fop=owwemva_lop, Aop=laplac2_lop, stepper=cgstep, nAop=img%xm%n*img%z%n, &
                           niter=niter,eps1=epscg,verb=.true.)!,rm0,err,resd,resm,mmov,rmov,verb)
          else if (cgsolv(1:4).eq."prec") then
!              call solver_prc( m=scatS,d=B_image, Fop=owwemva_lop, Sop=laplac2_lop, stepper=cgstep, nSop=img%xm%n*img%z%n, &
              call solver_prc( m=scatS,d=B_image, Fop=owwemva_lop, Sop=igrad1_lop, stepper=cgstep, nSop=img%xm%n*img%z%n, &
                           niter=niter,eps1=epscg,verb=.true.)!,rm0,err,resd,resm,mmov,rmov,verb)
          end if
       end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!       call owwemva_lop2(.true.,.false.,scatS,B_image)
!       if (exist_file("horv0")) then
!          allocate(tmp(img%z%n*img%xm%n))
!          call wghtmute_lop_init(pick1)
!          call wghtmute_lop2(.false.,.false.,scatS,tmp )
!          call wghtmute_close()
!          scatS=tmp
!          deallocate(tmp)
!       end if      
       if (exist_file("grad")) then
          call srite("grad",scatS,4*img%z%n*img%xm%n)
       end if
!       B_image=0.
!       call owwemva_lop2(.false.,.false.,scatS,B_image)
!       if (exist_file("Dimage")) then
!          call srite("Dimage",B_image,4*img%z%n*nsecdim*img%xm%n)
!       end if      
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       if (exist_file("horv0")) then
          call oweditvel_sngl(pick1,scatS)
       end if

       velnew=velnew+scatS

       call cgstep_close()
       call owwemva_close()
       if (exist_file("horv0")) then
          call oweditvel_sngl(pick1,velnew,vel0)
       end if
       call to_history("n4",iter,"Velout")
       call srite("Velout",velnew,4*img%z%n*img%xm%n)
       iter=iter+1

    end do
 end if
 111 format(a26,i4,a8,i4,a4,f6.2,a4)
 112 format(a29,f10.2,a20,i4)
 113 format(2(a8,i4),a8,a5,a16,f15.6)
 if (allocated(source)) deallocate(source)
 if (allocated(recv)) deallocate(recv)
 if (allocated(deltaS)) deallocate(deltaS)
 if (allocated(scatS)) deallocate(scatS)
 if (allocated(deltaI)) deallocate(deltaI)
 if (allocated(R_image)) deallocate(R_image)
 if (allocated(B_image)) deallocate(B_image)
 if (allocated(taper)) deallocate(taper)
 if (allocated(vel)) deallocate(vel)
 if (allocated(vref)) deallocate(vref)
 if (allocated(vel0)) deallocate(vel0)

 close(2)
 close(3)
11 format(a8,i4,a11,i5,a8,i4,2(a16,f20.10))

 call sep_end_prog()

end program OWcgwemva
