program OWrenorm05

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                             :: vel(:),vref(:,:,:),vel0(:),velnew(:)
complex, allocatable                          :: recv(:,:),source(:,:),source1(:),scatS(:),deltaS(:)
complex, allocatable                          :: B_image(:)
real, allocatable                             :: P_image(:)
real,allocatable                              :: err(:),temp(:),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

 iter=1
 if (dot) then
    allocate(d1(img%xm%n*img%z%n*nsecdim),d2(img%xm%n*img%z%n*nsecdim))
    allocate(m1(img%xm%n*img%z%n),m2(img%xm%n*img%z%n))
    allocate(temp(img%xm%n*img%z%n))
    call random_number(temp)
    m1=cmplx(acos(temp),asin(temp))
    deallocate(temp)
    allocate(temp(img%xm%n*nsecdim*img%z%n))
    call random_number(temp);temp=2*temp-1
    d2=cmplx(acos(temp),0)
    m1=m1/1000;d2=d2/1000;m2=0;d1=0;
    call owwemva_lop_init(taper,velnew)
    call owwemva_init()
    stat=owwemva_ren(.false.,.false.,m1,d1)
    stat=owwemva_ren(.true.,.false.,m2,d2)
    dott(1)=sum(conjg(m1)*m2)
    dott(2)=sum(conjg(d1)*d2)
    write(0,*) dott
    call owwemva_close()
    deallocate(m1,m2,d1,d2,temp)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 else
    do while(iter <= nl_iter)
          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)
             if (allocated(dslow)) deallocate(dslow) ;allocate(dslow(img%xm%n*img%z%n));dslow=0.
             call sreed("dslow",dslow,4*img%xm%n*img%z%n)
             call owwemva_forward(dslow,P_image)
             call owwemva_close()
             if (exist_file("dkim")) call srite("dkim",P_image,4*img%z%n*img%xm%n*nsecdim)
             deallocate(dslow)
          end if
          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))
          scatS=0.
          call owwemva_init()
          call owwemva_lop_init(taper,velnew)
          allocate(mmov(img%xm%n*img%z%n,nitercg1),rmov(img%xm%n*img%z%n,nitercg1),err(nitercg1))
          call solver_regzl2_cplx1( m=scatS,d=B_image, Fop=owwemva_ren, Aop=identity_lop2, stepper=cgstep_cplx, &
                   nAop=img%xm%n*img%z%n,niter1=nitercg1,eps1=epscg1,verb=.true.,mmov=mmov,rmov=rmov,err=err)
          call cgstep_close()
          call owwemva_close()
          if (exist_file("Vrenorm")) call srite("Vrenorm",scatS,8*img%z%n*img%xm%n)
          deallocate(mmov,rmov,err)
          B_image=0.
          stat=owwemva_ren(.false.,.false.,scatS,B_image)
          if (exist_file("dimg")) call srite("dimg",real(B_image),4*img%z%n*img%xm%n*nsecdim)       
          deallocate(B_image)
          velnew=velnew+real(scatS(nz*img%xm%n*(i1-1)+1:nz*img%xm%n*i1))!/rec%w%n
          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(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)
11 format(a8,i4,a11,i5,a8,i4,2(a16,f20.10))

 call sep_end_prog()

end program
