program OWrenorm03
!
! Performs the experiment:
! 1 - computing Dsr using  the correct DU and Ds
!
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
!!!!!!!!!
! solver and regularization
use solver_incore_mod
use steppr_mod
use identity_mod
! vel
use lloyd_vel

implicit none

character*80                                  :: listagem
real                                          :: start, finish, start1, finish1, start2, finish2
real                                          :: a1,b1,vmin,vmax,stat1,v0,r,stat
integer                                       :: ixs,ino,ni,i,st,ishot,t3,iz,j,k,iter,ix,ish,titer,i1,i2,i3,mm
real, allocatable                             :: velint(:),taper(:)
real, allocatable                             :: deltaI(:,:,:)
real, allocatable                             :: vel(:),vref(:,:,:),vel0(:),velnew(:)
complex, allocatable                          :: recv(:,:),source(:,:),source1(:),scatS(:),deltaS(:)
complex, allocatable                          :: B_image(:)
real, allocatable                             :: P_image(:)
real,allocatable                              :: err(:),temp(:),temp0(:),dslow(:)
complex, allocatable                          :: mmov(:,:),rmov(:,:)
complex,allocatable                           :: m1(:),m2(:),d1(:),d2(:)
double precision                              :: dott(2)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 call sep_init()
 call sep_begin_prog()
 start=0.;finish=0.
 call owexec_time(start)
 call owparam()
 call from_param("list",listagem,"listagem")
 open(unit=3,file=listagem)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Set up the taper
!
 allocate(taper(kxm%n))
 call owtaper(tapertr,taper)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Read velocities
!
 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))
 call sreed("vel",velint,v%z%n*v%xm%n*4)
 call owgrab_velim(velint,vel0,ni)
 deallocate(velint)
 velnew=vel0

if (not(dot)) then
!------------------------------------------------
!
! Initializing IO of wavefields
!
 call ow_wvfld_io_init("bkws")
 call ow_wvfld_io_init("bkwr")
 call ow_wvfld_io_init("dkws")
 call ow_wvfld_io_init("dkwr")
end if

! 1 - computing Dsr using  the correct DU and Ds
    call owwemva_init()
    call owwemva_lop_init(taper,velnew) 
    if (allocated(temp))  deallocate(temp)    ;allocate(temp(v%z%n*v%xm%n));  temp=0.
    if (allocated(temp0)) deallocate(temp0)   ;allocate(temp0(img%xm%n*img%z%n)); temp0=0.
    call sreed("dslow",temp,4*v%z%n*v%xm%n)
    call owgrab_velim(temp,temp0,ni)
    deallocate(temp)
    call owwolfe(deltaS=temp0,n=1)
    if (allocated(deltaS)) deallocate(deltaS) ;allocate(deltaS(img%xm%n*img%z%n*rec%w%n));deltaS=0.
    do i3=1,rec%w%n
       do i2=1,img%z%n
          do i1=1,img%xm%n
             deltaS((i3-1)*img%z%n*img%xm%n+(i2-1)*img%xm%n+i1)=cmplx(temp0((i2-1)*img%xm%n+i1))
          end do
       end do
    end do
    deallocate(temp0)
    call owwemva_close()
    if (allocated(scatS)) deallocate(scatS) ;allocate(scatS(img%z%n*img%xm%n*rec%w%n))
    call auxclose("bkws")
    call auxclose("bkwr")
    call auxclose("dkws")
    call auxclose("dkwr")
    stat = owdeltaS_renorm_lop(.false.,.false.,deltaS,scatS)
    if (exist_file("Vrenorm")) call srite("Vrenorm",scatS,8*img%z%n*img%xm%n*rec%w%n)
    call solver_regzl2_cplx( m=deltaS,d=scatS, Fop=owdeltaS_renorm_lop, Aop=igrad1_lop, stepper=cgstep_cplx, &
             nAop=img%xm%n*img%z%n*rec%w%n,niter1=nitercg2,eps1=epscg2,verb=.true.)
    call auxclose("bkws")
    call auxclose("bkwr")
    call auxclose("dkws")
    call auxclose("dkwr")
    do i1=1,rec%w%n
          velnew=velnew+real(deltaS(nz*img%xm%n*(i1-1)+1:nz*img%xm%n*i1))/rec%w%n
    end do
    call cgstep_close()
    if (exist_file("grad")) call srite("grad",deltaS,8*img%z%n*img%xm%n*rec%w%n)
    call to_history("n4",iter,"Velout")
    call srite("Velout",velnew,4*img%z%n*img%xm%n)

    call ow_wvfld_io_init("bkws")
    call ow_wvfld_io_init("bkwr")
    call ow_wvfld_io_init("dkws")
    call ow_wvfld_io_init("dkwr")
    call owwemva_init()
    call owwemva_lop_init(taper,velnew) 
    call owwolfe(celtaS=scatS,n=rec%w%n)
    call owwemva_close()

 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(deltaS)) deallocate(deltaS)
 if (allocated(scatS)) deallocate(scatS)
 if (allocated(taper)) deallocate(taper)
 if (allocated(vel)) deallocate(vel)
 if (allocated(vref)) deallocate(vref)
 if (allocated(vel0)) deallocate(vel0)

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

 call sep_end_prog()

end program
