program OWrenorm04
!
! Performs the experiment:
! 1 - computing Ds using corrct Du and inverted Dsr

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(:,:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 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

!------------------------------------------------
! 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")

    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(temp0)

    if (allocated(B_image)) deallocate(B_image) ;allocate(B_image(img%xm%n*nsecdim*img%z%n));B_image=0.
    if (allocated(P_image)) deallocate(P_image) ;allocate(P_image(img%xm%n*nsecdim*img%z%n));P_image=0.
    if (exist_file("P_image")) then
       call sreed("P_image",P_image,4*img%z%n*nsecdim*img%xm%n)
    else
       call owwemva_init()
       call owwemva_lop_init(taper,velnew)
       call owwemva_forward(temp0,P_image)
       call owwemva_close()
       if (exist_file("dkim")) call srite("dkim",P_image,4*img%z%n*img%xm%n*nsecdim)
    end if
    deallocate(temp0)
    B_image=cmplx(P_image,0.)
    deallocate(P_image)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Computing renormalized deltaS
    if (exist_file("mmov")) call ow_solver_io_init("mmov")
    if (exist_file("rmov")) call ow_solver_io_init("rmov")
    if (exist_file("err")) call ow_err_io_init()
    if (allocated(scatS)) deallocate(scatS) ;allocate(scatS(img%z%n*img%xm%n*sou%w%n))
    scatS=0.
    call owwemva_init()
    call owwemva_lop_init(taper,velnew)
    allocate(mmov(img%xm%n*img%z%n*sou%w%n,nitercg1),rmov(img%xm%n*img%z%n*sou%w%n,nitercg1),err(nitercg1))
    call solver_regzl2_cplx( m=scatS,d=B_image, Fop=owwemva_lop, Aop=identity_lop2, stepper=cgstep_cplx, &
            nAop=img%xm%n*img%z%n*sou%w%n,niter1=nitercg1,eps1=epscg1,verb=.true.,mmov=mmov,rmov=rmov,err=err)
    call cgstep_close()
    if (exist_file("Vrenorm")) call srite("Vrenorm",scatS,8*img%z%n*img%xm%n*rec%w%n)
    deallocate(mmov,rmov,err)
    if (exist_file("dimg")) then
       B_image=0.
       stat=owwemva_lop(.false.,.false.,scatS,B_image)
       call srite("dimg",real(B_image),4*img%z%n*img%xm%n*nsecdim)       
       deallocate(B_image)
    end if
    call owwemva_close()
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Computing the perturbed wavefields from
    !  the renormalized deltaS
    if (allocated(deltaS)) deallocate(deltaS) ;allocate(deltaS(img%z%n*img%xm%n*rec%w%n))
    call owwemva_lop_init(taper,velnew)
    call auxclose("bkws")
    call auxclose("bkwr")
    call auxclose("dkws")
    call auxclose("dkwr")
    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 srite("Velout",velnew,4*img%z%n*img%xm%n)

 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(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(3)

 call sep_end_prog()

end program

!
!
!       else
!          start1=0.;finish1=0.;
!          if (allocated(B_image)) deallocate(B_image) ;allocate(B_image(img%xm%n*nsecdim*img%z%n));B_image=0.
!          if (allocated(deltaI)) deallocate(deltaI) ;allocate(deltaI(kxm%n,nsecdim,img%z%n));deltaI=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
!          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!          !
!          ! Computing DSO
!          !
!          call owdso(.false. , B_image)
!          call owdso(.true.  , B_image)
!          write(0,112) "COMPUTED PERTURBED  IMAGE IN ",(finish1-start1)/60," min OF ITERATION # ", iter
!       end if
