program OWrenorm01

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_reg1_mod
use solver_tiny_mod
use laplac2
use cgstep_mod
use identity_mod
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,stat
integer                                       :: ixs,ino,ni,i,st,ishot,t3,iz,j,k,iter,ix,ish,titer,i1,i2,i3,nl_iter,mm
real, allocatable                             :: velint(:),taper(:),pick1(:)
real, allocatable                             :: deltaI(:,:,:),scatS(:),deltaS(:)
real, allocatable                             :: R_image(:),B_image(:)
real, allocatable                             :: vel(:),vref(:,:,:),vel0(:),indip(:),velnew(:)
complex, allocatable                          :: recv(:,:),source(:,:),source1(:),cscatS(:)
logical                                       :: cond,bounds
real,allocatable                              :: m1(:),m2(:),d1(:),d2(:)
real, allocatable                             :: mmov(:,:),rmov(:,:),err(:)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 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))

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

if (not(dot).and.renorm) 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

 iter=1
 call from_param('nl_iter',nl_iter,1)
 if (dot) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    allocate(d1(img%xm%n*img%z%n*nsecdim),m1(img%xm%n*img%z%n*sou%w%n),d2(img%xm%n*img%z%n*nsecdim),m2(img%xm%n*img%z%n*sou%w%n))
    call random_number(m1);call random_number(d2)
    m1=m1/10;d2=(2*d2-1)/10;m2=0;d1=0;
    call owwemva_lop_init(taper,velnew)
    call owwemva_init()
    stat=owwemva_renorm_lop(.false.,.false.,m1,d1)
    stat=owwemva_renorm_lop(.true.,.false.,m2,d2)
    write(0,*) dot_product(m1,m2),dot_product(d1,d2)
    call owwemva_close()
    deallocate(m1,m2,d1,d2)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    allocate(d1(img%xm%n*img%z%n*nsecdim),m1(img%xm%n*img%z%n),d2(img%xm%n*img%z%n*nsecdim),m2(img%xm%n*img%z%n))
    call random_number(m1);call random_number(d2)
    m1=m1/10;d2=(2*d2-1)/10;m2=0;d1=0;
    call owwemva_lop_init(taper,velnew)
    call owwemva_init()
    stat=owwemva_lop(.false.,.false.,m1,d1)
    stat=owwemva_lop(.true.,.false.,m2,d2)
    write(0,*) dot_product(m1,m2),dot_product(d1,d2)
    call owwemva_close()
    deallocate(m1,m2,d1,d2)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    if (allocated(B_image)) deallocate(B_image) ;allocate(B_image(img%xm%n*nsecdim*img%z%n));B_image=0.
    call sreed("P_image",B_image,4*img%z%n*nsecdim*img%xm%n)
    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)
    stat=owwemva_renorm_lop(.true.,.false.,scatS,B_image)
    call to_history("n1",img%xm%n,"Vrenorm"); call to_history("d1",img%xm%d,"Vrenorm"); call to_history("o1",img%xm%o,"Vrenorm")  
    call to_history("n2",nz,"Vrenorm") ; call to_history("d2",img%z%d,"Vrenorm") ; call to_history("o2",zmin,"Vrenorm")
    call to_history("n3",1,"Vrenorm"); call to_history("d3",1,"Vrenorm"); call to_history("o3",0,"Vrenorm")  
    call to_history("n4",1,"Vrenorm"); call to_history("d4",1,"Vrenorm"); call to_history("o4",0,"Vrenorm")
    call to_history("n5",rec%w%n,"Vrenorm") ; call to_history("d5",rec%w%d,"Vrenorm") ; call to_history("o5",rec%w%o,"Vrenorm")
    call to_history("n6",iter,"Vrenorm"); call to_history("d6",1.,"Vrenorm"); call to_history("o6",1.,"Vrenorm")
    call to_history("esize",4,"Vrenorm")
    call srite("Vrenorm",scatS,4*img%z%n*img%xm%n*rec%w%n)
    deallocate(scatS,B_image)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 else
    do while(iter <= nl_iter)

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
    ! Computing bkg.image
    !
    if (exist_file("P_image")) then
       if (allocated(B_image)) deallocate(B_image) ;allocate(B_image(img%xm%n*nsecdim*img%z%n));B_image=0.
       call sreed("P_image",B_image,4*img%z%n*nsecdim*img%xm%n)
    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)
!       if (exist_file("horv0")) then
!          allocate(R_image(img%xm%n*nsecdim*img%z%n*sou%w%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
       call owdso(.true.  , B_image)
       write(0,112) "COMPUTED PERTURBED  IMAGE IN ",(finish1-start1)/60," min OF ITERATION # ", iter
    end if
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
    ! 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 (renorm) then
       if (not(exist_file("deltaS"))) then
          if (allocated(scatS)) deallocate(scatS) ;allocate(scatS(img%z%n*img%xm%n*sou%w%n))
!          if (allocated(cscatS)) deallocate(cscatS) ;allocate(cscatS(img%z%n*img%xm%n))
          scatS=0.
          call owwemva_init()
          call owwemva_lop_init(taper,velnew)
!          call laplac2_init(img%xm%n,img%z%n)
!          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_renorm_lop, stepper=cgstep, niter1=niter)
!             else if (cgsolv(1:3).eq."reg") then
!                call solver_reg( m=scatS,d=B_image, Fop=owwemva_renorm_lop, Aop=identity_lop, stepper=cgstep, nAop=img%xm%n*img%z%n*sou%w%n, &
!                           Jop=wghtmute_lop,niter1=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_renorm_lop, Sop=identity_lop, stepper=cgstep, nSop=img%xm%n*img%z%n*sou%w%n, &
!                           Jop=wghtmute_lop,niter1=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_renorm_lop, stepper=cgstep, niter1=niter)
!             else if (cgsolv(1:3).eq."reg") then
              allocate(mmov(img%xm%n*img%z%n*sou%w%n,niter),rmov(img%xm%n*img%z%n*sou%w%n,niter),err(niter))
              call solver_reg1( m=cscatS,d=B_image, Fop=owwemva_renorm_lop, Aop=identity_lop, stepper=cgstep, nAop=img%xm%n*img%z%n*sou%w%n, &
                           niter1=niter,eps1=epscg,verb=.true.,mmov=mmov,rmov=rmov,err=err)
!          call solver_reg1( m=scatS,d=B_image, Fop=owwemva_lop, Aop=identity_lop, stepper=cgstep, nAop=img%xm%n*img%z%n, &
!                           niter1=niter,eps1=epscg,verb=.true.,mmov=mmov,rmov=rmov,err=err)
!             else if (cgsolv(1:4).eq."prec") then
!                 call solver_prc( m=scatS,d=B_image, Fop=owwemva_renorm_lop, Sop=identity_lop, stepper=cgstep, nSop=img%xm%n*img%z%n*sou%w%n, &
!                           niter1=niter,eps1=epscg,verb=.true.)!,rm0,err,resd,resm,mmov,rmov,verb)
!             end if
!          end if
          call cgstep_close()
          call owwemva_close()
          if (exist_file("Vrenorm")) then
             call to_history("n1",img%xm%n,"Vrenorm"); call to_history("d1",img%xm%d,"Vrenorm"); call to_history("o1",img%xm%o,"Vrenorm")  
             call to_history("n2",nz,"Vrenorm") ; call to_history("d2",img%z%d,"Vrenorm") ; call to_history("o2",zmin,"Vrenorm")
             call to_history("n3",1,"Vrenorm"); call to_history("d3",1,"Vrenorm"); call to_history("o3",0,"Vrenorm")  
             call to_history("n4",1,"Vrenorm"); call to_history("d4",1,"Vrenorm"); call to_history("o4",0,"Vrenorm")
             call to_history("n5",rec%w%n,"Vrenorm") ; call to_history("d5",rec%w%d,"Vrenorm") ; call to_history("o5",rec%w%o,"Vrenorm")
             call to_history("n6",iter,"Vrenorm"); call to_history("d6",1.,"Vrenorm"); call to_history("o6",1.,"Vrenorm")
             call to_history("esize",4,"Vrenorm")
             call srite("Vrenorm",scatS,4*img%z%n*img%xm%n*rec%w%n)
          end if
          deallocate(mmov,rmov,err)
       else
          if (allocated(scatS)) deallocate(scatS) ;allocate(scatS(img%z%n*img%xm%n))
          call sreed("deltaS",scatS,4*img%z%n*img%xm%n)
       end if

       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       !
       ! Computing the perturbed wavefields from
       !  the renormalized deltaS
       !
       if (allocated(deltaS)) deallocate(deltaS) ;allocate(deltaS(img%z%n*img%xm%n))
       call owwemva_lop_init(taper,velnew)
       write(0,*) "renorm1 passei 00";          call owwfld_renorm(scatS)
!       write(0,*) "renorm1 passei 00";          call owwfld_renorm1(scatS)
       call auxclose("bkws")
       call auxclose("bkwr")
       call auxclose("dkws")
       call auxclose("dkwr")
       write(0,*) "renorm1 passei 02"
       call solver_reg( m=deltaS,d=scatS, Fop=owdeltaS_renorm_lop, Aop=identity_lop, stepper=cgstep, nAop=img%xm%n*img%z%n, &
                   niter1=niter,eps1=epscg,verb=.true.)!,rm0,err,resd,resm,mmov,rmov,verb)
       write(0,*) "renorm1 passei 03"
       call auxclose("bkws")
       call auxclose("bkwr")
       call auxclose("dkws")
       call auxclose("dkwr")
       write(0,*) "renorm1 passei 04";          velnew=velnew+deltaS
       write(0,*) "renorm1 passei 05";          call cgstep_close()
       write(0,*) "renorm1 passei 06"
       if (exist_file("grad")) call srite("grad",deltaS,4*img%z%n*img%xm%n)
    else
       if (allocated(scatS)) deallocate(scatS) ;allocate(scatS(img%z%n*img%xm%n))
       scatS=0.
       call owwemva_init()
       call owwemva_lop_init(taper,velnew)
       call laplac2_init(img%xm%n,img%z%n)
       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, niter1=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,niter1=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=identity_lop, stepper=cgstep, nSop=img%xm%n*img%z%n, &
                           Jop=wghtmute_lop,niter1=niter,eps1=epscg,verb=.true.)!,rm0,err,resd,resm,mmov,rmov,verb)
          end if
          call wghtmute_close()
       else
          allocate(mmov(img%xm%n*img%z%n,niter),rmov(img%xm%n*img%z%n,niter),err(niter))
          if (cgsolv(1:3).eq."reg") then
             if (lapl) then
write(0,*) "laplacian ", lapl
                call solver_reg1( m=scatS,d=B_image, Fop=owwemva_lop, Aop=laplac2_lop, stepper=cgstep, nAop=img%xm%n*img%z%n, &
                           niter1=niter,eps1=epscg,verb=.true.,mmov=mmov,rmov=rmov,err=err)
             else
                call solver_reg( m=scatS,d=B_image, Fop=owwemva_lop, Aop=identity_lop, stepper=cgstep, nAop=img%xm%n*img%z%n, &
                           niter1=niter,eps1=epscg,verb=.true.,mmov=mmov,rmov=rmov,err=err)
             end if
          else if (cgsolv(1:4).eq."prec") then
              call solver_prc( m=scatS,d=B_image, Fop=owwemva_lop, Sop=identity_lop, stepper=cgstep, nSop=img%xm%n*img%z%n, &
                           niter1=niter,eps1=epscg,verb=.true.,mmov=mmov,rmov=rmov,err=err)
          end if
          deallocate(mmov,rmov,err)
          write(0,*) "passei 000";          call cgstep_close()
          write(0,*) "passei 001";          call owwemva_close()
          write(0,*) "passei 002"
       end if
       write(0,*) "passei 00"
       velnew=velnew+scatS
       if (exist_file("grad")) call srite("grad",scatS,4*img%z%n*img%xm%n)
    end if

!       if (exist_file("horv0")) call oweditvel_sngl(pick1,scatS)
!       if (exist_file("horv0")) call oweditvel_sngl(pick1,velnew,vel0)

    call to_history("n4",iter,"Velout")
    call srite("Velout",velnew,4*img%z%n*img%xm%n)
    iter=iter+1
    write(0,*) "passei 04"
    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
