module OW_parms_areal

use sep
use OW_type
use base_blkcal_mod

implicit none
character*12,public, save                                   :: sufx
character*6,public, save                                    :: extrap,cgsolv
character*7,public, save                                    :: imgcon
character*80,public, save                                   :: hor1, listagem
logical,public, save				:: adj,dot,icwind,impulse,dso,resmig,migration
logical,public, save				:: light,moddso,smooth,imgnorm,inversion,multsc
logical,public, save                                        :: lapl,geomspread,zweight,resmig_lin,areal
logical,public, save				:: timeshift,l1norm,envelope,born2,renorm,scaling,chop
logical,public, save                                        :: bounds,cbounds,onlygrad
type (axis) ,public, save				:: kxm, hx, hy, kym, w, xd
type (axis) ,public, save				:: xm, ym, z, shx, shy, kz, khx, khy
type (dat_grid)	 ,public, save			:: rec, sou
type (mod_grid)	 ,public, save			:: img
type (velo) ,public, save				:: v
type (block) ,public, save				:: wblock,zblock
integer	 ,public, save				:: max_iters, line_search,iters,secgr
integer	 ,public, save				:: windx,windh,windz,frame,tc1,tc2
integer	 ,public, save				:: freq,jhx,jhy,nx,ny,nref,iwb
integer	 ,public, save				:: tapertrx,padtracesx,tapertry,padtracesy,nreal
integer	 ,public, save				:: nz,nzmin,nzmax,nsecdimy,nsecdimx,nitercg,nzold,mmax
integer	 ,public, save				:: isp,nwbott,node,niter,n_wind,n_sampt,nxfilt,nh0,dipop
integer, public, save                                       :: ninx,nfnx,niny,nfny,mdso
real ,public, save					:: xmin,xmax,ymin,ymax,zmin,zmax,pi,dwbott,yline_
real ,public, save					:: osecdimx,dsecdimx,osecdimy,dsecdimy,thresh,maxvchg
real, public, save					:: vsal,rho,maxang,damp,beta,deltav,expgrad
real, public, save					:: xv_min,xv_max,yv_min,yv_max,trick,expdso,zwgt
real,public, save                                           :: vmin,vmax,p_up,p_lo,chophi,choplo
real,public, save                                           :: err,err1,epsteep,factor,lambda_min,tola,tol(2)
integer, public, save				:: nxw, nah,nax
real, public, save					:: dah,oah,dax,oax
type(axis), public, save               		          :: axk, ayk, azk
real, public, save                               		:: axkd, aykd, azkd, axko, ayko, azko
integer, public, save                            		:: axkn, aykn, azkn
! using lloyd
logical,public, save				:: lloyd,perc_start
real,public,save					:: min_r_pct,min_s_pct,slow_own
integer,public,save					:: niter_over
! DIPS
logical,public,save					:: dipflt
!logical,public,save				:: dips,verb,v3d,dave
!character*20,public,save				:: method
!integer,public,save				:: maxsize,nl_iter,niter_dip,nwind
!real, public,save					:: eps,eps2
real, public, save                                          :: epscg
!
integer, public, save                                       :: serg_na, serg_niter, dip_niter
real, public, save                                          :: serg_eps1, serg_eps2, serg_epsh, dip_eps
integer, public, save                                       :: shtblk,shtres,shtblm,nls
real, public, save                                          :: alpha

contains

subroutine owparam()

real						:: wbott,dd1,dd2,dd3,od1,od2,od3,dv1,dv2,ov1,ov2
real						:: deltat, d_sampt,dzref,ozref,zmin1
integer						:: iw,nt, i, m,nd1,nd2,nd3,nv1,nv2,nzref,jh,nreal1

if (wpar%i_mpi==0) then
 call from_param("dipflt",dipflt,.false.)
 if (dipflt) call from_param('frame',frame,50)
call from_param("epscg",epscg,1.)
call from_param("nitercg",nitercg,50)
call from_param("cgsolv",cgsolv,"      ")
if (cgsolv(1:3).eq."reg") call from_param("lapl",lapl,.true.)
     call from_param("oxk", axko, 0.); call from_param("dxk", axkd, 1.); call from_param("nxk", axkn, 99999 )
     call from_param("oyk", ayko, 0.); call from_param("dyk", aykd, 1.); call from_param("nyk", aykn, 99999 )
     call from_param("ozk", azko, 0.); call from_param("dzk", azkd, 1.); call from_param("nzk", azkn, 99999 )
     axk%n=axkn
     ayk%n=aykn
     azk%n=azkn
     axk%d=axkd
     ayk%d=aykd
     azk%d=azkd
     axk%o=axko
     ayk%o=ayko
     azk%o=azko

     write(0,*) axkn,aykn,azkn
     write(0,*) axkd,aykd,azkd
     write(0,*) axko,ayko,azko
     write(0,*) (axk%n.ne.99999 .and. ayk%n.ne.99999 .and. azk%n.ne.99999) 
!!!!!!!!!!!!!!!!!!
!
! for smoothing along dips
  call from_param("smooth",smooth,.true.)
  if (smooth) then
     call from_param("serg_na"    ,serg_na,8)
     call from_param("serg_epsh"  ,serg_epsh,0.01)
  ! Set smoothing factors for SERGEY
     call from_param("serg_niter",serg_niter,30) 
     call from_param("serg_eps1",serg_eps1,0.002)
     call from_param("serg_eps2",serg_eps2,0.002)
     call from_param("dip_niter",dip_niter,100)
     call from_param("dip_eps",dip_eps,10.)
!!!!!!!!!!!!!!!!!!
  end if
 call from_param("mmax",mmax,10)
 call from_param("sufx",sufx,"_")
 call from_param("dot",dot,.false.)
 call from_param("geomspread",geomspread,.true.)
 call from_param("adj",adj,.true.)
 call from_param("inversion",inversion,.false.)
 call from_param("dso",dso,.false.)
 call from_param("expgrad",expgrad,0.125)
 if (dso) call from_param("mdso",mdso,0)
 call from_param("l1norm",l1norm,.false.)
 call from_param("windx",windx,11)
 call from_param("windh",windh,3)
 call from_param("windz",windz,11)
 call from_param("extrapolator",extrap,"ffd")
 call from_param("zwgt",zwgt,99999.)
 if (extrap(1:3).eq."ffd") then
   call from_param("dipop",dipop,45)
   call from_param("trick",trick,0.14)
 end if
 if (dso) then
   call from_param("maxvchg",maxvchg,0.15)
   call from_param("moddso",moddso,.false.)
   call from_param("expdso",expdso,1.)
   call from_param("envelope",envelope,.false.)
   if (envelope) then
     call from_param("tc1",tc1,25)
     call from_param("tc2",tc2,5)
     call from_param("freq",freq,0)
   end if
 end if
 call from_param("born2",born2,.false.)
 call from_param("renorm",renorm,.false.)
 if (moddso) call from_param("beta",beta,0.1)
 call from_param("migration",migration,.false.)
 call from_param("areal",areal,.true.)
 call from_param("timeshift",timeshift,.false.)
 call from_param("imgcon",imgcon,"imgprim")
 call from_param("imgnorm",imgnorm,.false.)
if (inversion) then
   if (not(dso)) then
      call from_param("resmig",resmig,.false.)
      call from_param("rho",rho,1.)
      if (not(resmig)) then
         call from_param("resmig_lin",resmig_lin,.false.)
         if (not(resmig_lin)) then
            call seperr("How to compute the perturbed image?")
         end if
      end if
   end if
end if
 call from_param("impulse",impulse,.false.)
 call from_param("ymin",ymin,0.)
 call from_param("ymax",ymax,0.)
 call from_param("node",node,9999)
 call from_param("niter",niter,5)
 call from_param("nxfilt",nxfilt,99999)
 call from_param("wbott",wbott,0.)
 call from_param("traces_in_taper_x",tapertrx,0)
 call from_param("padtraces_x",padtracesx,tapertrx)
 call from_param("traces_in_taper_y",tapertry,0)
 call from_param("padtraces_y",padtracesy,tapertry)
 call from_param("maxang",maxang,89.)
 call from_param("damp", damp, 25000.)
 call from_param("deltav", deltav, 0.15)
 if (deltav > 1.) deltav=deltav/100
 call from_param("zweight",zweight,.false.)
 call from_param("scale",scaling,.true.)
 call from_param("onlygrad",onlygrad,.false.)
 call from_param("chop",chop,.true.)
 call from_param("multsc",multsc,.false.)
if (chop) then
   call from_param("chophi",chophi,.1)
   call from_param("choplo",choplo,.1)
   if (choplo > 1) choplo=choplo/100
   if (chophi > 1) chophi=chophi/100
end if
 call from_param("secgr",secgr,0)
 call from_param("nls",nls,3)
 call from_param("alpha",alpha,.1)

 call from_aux("rec","n1",rec%hx%n);	call from_aux("rec","d1",rec%hx%d);     call from_aux("rec","o1",rec%hx%o)
 call from_aux("rec","n2",rec%hy%n);	call from_aux("rec","d2",rec%hy%d);     call from_aux("rec","o2",rec%hy%o)
 call from_aux("rec","n3",rec%w%n);	call from_aux("rec","d3",rec%w%d);      call from_aux("rec","o3",rec%w%o)
 call from_aux("rec","n4",rec%xd%n);    call from_aux("rec","d4",rec%xd%d);     call from_aux("rec","o4",rec%xd%o)  
 call from_aux("rec","n5",nreal1)
 call from_aux("rec","o6",zmin1)

 call from_aux("sou","n1",sou%hx%n);	call from_aux("sou","d1",sou%hx%d);	call from_aux("sou","o1",sou%hx%o)
 call from_aux("sou","n2",sou%hy%n);	call from_aux("sou","d2",sou%hy%d);	call from_aux("sou","o2",sou%hy%o)
 call from_aux("sou","n3",sou%w%n);	call from_aux("sou","d3",sou%w%d);	call from_aux("sou","o3",sou%w%o)
 call from_aux("sou","n4",sou%xd%n);	call from_aux("sou","d4",sou%xd%d);	call from_aux("sou","o4",sou%xd%o)  
 call from_aux("sou","n5",nreal)
 call from_aux("sou","o6",zmin)

 call from_aux("vel","n1",v%xm%n);	call from_aux("vel","d1",v%xm%d);	call from_aux("vel","o1",v%xm%o)
 call from_aux("vel","n2",v%ym%n);	call from_aux("vel","d2",v%ym%d);	call from_aux("vel","o2",v%ym%o)
 call from_aux("vel","n3",v%z%n);	call from_aux("vel","d3",v%z%d);	call from_aux("vel","o3",v%z%o)
 
 if (zmin.ne.zmin1) call seperr("source and receiver recorded at different depth levels")
 call from_param("zmax",zmax,(v%z%n-1)*v%z%d+zmin)

if (extrap(1:4).eq."pspi") then
   call from_param("nrefvel",nref,1) 
   call from_param("lloyd",lloyd,.false.)
   if (lloyd) then
       call from_param("min_region_pct",min_r_pct,1.)
       min_r_pct=min_r_pct/100.
       !MINIMUM ALLOWED DEVIATION BETWEEN CENTROIDS
       call from_param("min_slow_dev",min_s_pct,9.)
       min_s_pct=min_s_pct/100.
       call from_param("perc_start",perc_start,.false.)
       call from_param("vel_own",slow_own,-1.)
       slow_own=1./slow_own
       !NUMBER OF ITERATIONS
       call from_param("niter_lloyd",niter_over,-1)
   else
       call from_param("vsal",vsal)
   end if
end if
 if (nreal1 .ne. nreal) then
    call seperr("# source realizations .ne. # receiver realizations")
 end if
 if (sou%w%n .ne. rec%w%n) then
    call seperr("# source freq. .ne. # receiver freq.")
 end if
 if (sou%xd%n .ne. rec%xd%n) then
    call seperr("# source planes. .ne. # receiver planes")
 end if
 if (sou%hx%n .ne. rec%hx%n) then
    call seperr("# source x-traces. .ne. # receiver x-traces")
 end if
 if (sou%hy%n .ne. rec%hy%n) then
    call seperr("# source y-traces. .ne. # receiver y-traces")
 end if
 if (sou%hx%d .ne. rec%hx%d) then
    call seperr("source x-trace interval .ne. areal shot x-trace interval")
 end if
 if (sou%hy%d .ne. rec%hy%d) then
    call seperr("source y-trace interval .ne. areal shot y-trace interval")
 end if
 if (sou%hx%o .ne. rec%hx%o) then
    call seperr("source x-trace origin .ne. areal shot x-trace origin")
 end if
 if (sou%hy%o .ne. rec%hy%o) then
    call seperr("source y-trace origin .ne. areal shot y-trace origin")
 end if
 if (sou%w%d .ne. rec%w%d) then
    call seperr("source freq.interval .ne. receiver freq.interval")
 end if
 if (sou%w%o .ne. rec%w%o) then
    call seperr("source freq.origin .ne. receiver freq.origin")
 end if

   img%xm%n=rec%hx%n+2*padtracesx+2*tapertrx
   img%xm%d=rec%hx%d
   img%xm%o=rec%hx%o-(padtracesx+tapertrx)*img%xm%d
   img%ym%n=rec%hy%n+2*padtracesy+2*tapertry
   img%ym%d=rec%hy%d
   img%ym%o=rec%hy%o-(padtracesy+tapertry)*img%ym%d

 xmin=img%xm%o
 xmax=((img%xm%n-1)+(padtracesx+tapertrx))*img%xm%d+img%xm%o
 if (xmin < v%xm%o ) call seperr("Min.x-image < Min.x-velocity")
 if (xmax > v%xm%o+(v%xm%n-1)*v%xm%d ) call seperr("Max.x-image > Max.x-velocity")

 ymin=img%ym%o
 ymax=((img%ym%n-1)+(padtracesy+tapertry))*img%ym%d+img%ym%o
 if (ymin < v%ym%o ) call seperr("Min.y-image < Min.y-velocity")
 if (ymax > v%ym%o+(v%ym%n-1)*v%ym%d ) call seperr("Max.y-image > Max.y-velocity")

 nx=(xmax-xmin)/img%xm%d+1.5
 ny=(ymax-ymin)/img%ym%d+1.5

if (inversion) then
   call from_param("xv_min",xv_min,img%xm%o+10*img%xm%d)
   call from_param("xv_max",xv_max,(img%xm%n-1)*img%xm%d+img%xm%o-10*img%xm%d)
   call from_param("yv_min",yv_min,img%ym%o+10*img%ym%d)
   call from_param("yv_max",yv_max,(img%ym%n-1)*img%ym%d+img%ym%o-10*img%ym%d)
end if

  if (img%xm%o < v%xm%o ) call seperr("Min.x-image space < Min.x-velocity")
  if ((img%xm%n-1)*img%xm%d+img%xm%o > v%xm%o+(v%xm%n-1)*v%xm%d ) call seperr("Max.x-image space > Max.x-velocity")
 
  if (img%ym%o < v%ym%o ) call seperr("Min.y-image space < Min.y-velocity")
  if ((img%ym%n-1)*img%ym%d+img%ym%o > v%ym%o+(v%ym%n-1)*v%ym%d ) call seperr("Max.y-image space > Max.y-velocity")
 
 if (zmin < v%z%o ) call seperr("Min.z-image < Min.z-velocity")
 if (zmax > v%z%o+(v%z%n-1)*v%z%d ) call seperr("Max.z-image > Max.z-velocity")

 nz=int((zmax-zmin)/v%z%d)+1.5
 nzmax=int((zmax-v%z%o)/v%z%d)+1.5
 nzmin=int((zmin-v%z%o)/v%z%d)+1.5
 nzold=nzmax
 img%z%o=zmin

 write(0,*) "*************************"
 write(0,*) "* ",nz,nzmax,nzmin," *"
 write(0,*) "*************************"


  call from_param("wbott",wbott,0.)

   call from_param("nhx",img%shx%n,1)
   call from_param("jhx",jhx,1)
   call from_param("dhx",img%shx%d,jhx*img%xm%d)
   call from_param("ohx",img%shx%o,-int(img%shx%n/2)*img%shx%d)
   call from_param("nhy",img%shy%n,1)
   call from_param("jhy",jhy,1)
   call from_param("dhy",img%shy%d,jhy*img%ym%d)
   call from_param("ohy",img%shy%o,-int(img%shy%n/2)*img%shy%d)

  img%z%n=nz
  img%z%o=zmin
  img%z%d=v%z%d

  wbott=int(wbott/v%z%d)*v%z%d
  nwbott=int(wbott/v%z%d)+1
  dwbott=wbott/2.

pi=4*atan(1.)

!--------------------------------------------------
! Depth vector
!
allocate(img%z%all(img%z%n))
do iwb = 1,nz
   img%z%all(iwb) = img%z%o + (iwb - 1) * v%z%d
end do

!--------------------------------------------------
! Wavenumber vectors
!
! BUG DANGER: img%z%n ==> nz !!!!!
!
!
!if (resmig.or.resmig_lin) then
!!!kz
   kz%n = img%z%n
   allocate(kz%all(kz%n))
   kz%all=0.
   if(kz%n/=1) then
     m = kz%n/2
     do i=1,m+1
        kz%all(i)=i-1
     end do
     do i=m+2,kz%n
        kz%all(i)=i-kz%n-1
     end do
   else
     kz%all=0. !! make sure k=0 if n=1
   end if
   kz%all = kz%all / (kz%n*img%z%d)      !! 1/length
   kz%all = 2 * PI * kz%all     !! radians
   kz%d = kz%all(2)-kz%all(1)

!!! khy
   khy%n = img%shy%n
   allocate(khy%all(khy%n))
   khy%all=0.
   if(khy%n/=1) then
     m = khy%n/2
     do i=1,m+1
        khy%all(i)=i-1
     end do
     do i=m+2,khy%n
        khy%all(i)=i-khy%n-1
     end do
   else
     khy%all=0. !! make sure k=0 if n=1
   end if
   khy%all = khy%all / (khy%n*img%shy%d)      !! 1/length
   khy%all = 2 * PI * khy%all     !! radians
   if (khy%n/=1) then
      khy%d = khy%all(2)-khy%all(1)
   else
      khy%d=1/img%shy%d
   end if

!!! khx
   khx%n = img%shx%n
   allocate(khx%all(khx%n))
   khx%all=0.
   if(khx%n/=1) then
     m = khx%n/2
     do i=1,m+1
        khx%all(i)=i-1
     end do
     do i=m+2,khx%n
        khx%all(i)=i-khx%n-1
     end do
   else
     khx%all=0. !! make sure k=0 if n=1
   end if
   khx%all = khx%all / (khx%n*img%shx%d)      !! 1/length
   khx%all = 2 * PI * khx%all     !! radians
   if (khx%n/=1) then
      khx%d = khx%all(2)-khx%all(1)
   else
      khx%d=1/img%shx%d
   end if

kxm%n = rec%hx%n+2*(tapertrx+padtracesx)
kxm%d = 1./(kxm%n*rec%hx%d)
kxm%o = 1./(2*rec%hx%d)-kxm%d
allocate(kxm%all(kxm%n))
do i=1,kxm%n
   kxm%all(i) = mod(kxm%o+(i-1)*kxm%d,1./rec%hx%d)-kxm%o
end do
kxm%all = 2*pi * kxm%all

kym%n = rec%hy%n+2*(tapertry+padtracesy)
kym%d = 1./(kym%n*rec%hy%d)
kym%o = 1./(2*rec%hy%d)-kym%d
allocate(kym%all(kym%n))
do i=1,kym%n
   kym%all(i) = mod(kym%o+(i-1)*kym%d,1./rec%hy%d)-kym%o
end do
kym%all = 2*pi * kym%all

   nsecdimx=img%shx%n
   dsecdimx=img%shx%d
   osecdimx=img%shx%o

   nsecdimy=img%shy%n
   dsecdimy=img%shy%d
   osecdimy=img%shy%o

!-------->
 call from_param("list",listagem,"listagem")
 open(unit=3,file=listagem)
 call from_param("line_search",line_search,3)
 call from_param("max_iters",max_iters,10)

!----------------------------------------------
!
! Migration output
!
if (migration) then
   call auxscr("bkgimg")
   call to_history("n1",img%xm%n,"bkgimg"); call to_history("d1",img%xm%d,"bkgimg"); call to_history("o1",img%xm%o,"bkgimg")  
   call to_history("n2",img%shx%n,"bkgimg") ; call to_history("d2",img%shx%d,"bkgimg") ; call to_history("o2",img%shx%o,"bkgimg")
   call to_history("n3",img%ym%n,"bkgimg"); call to_history("d3",img%ym%d,"bkgimg"); call to_history("o3",img%ym%o,"bkgimg")  
   call to_history("n4",img%shy%n,"bkgimg") ; call to_history("d4",img%shy%d,"bkgimg") ; call to_history("o4",img%shy%o,"bkgimg")
   call to_history("n5",nz,"bkgimg") ; call to_history("d5",img%z%d,"bkgimg") ; call to_history("o5",zmin,"bkgimg")
   call to_history("esize",4,"bkgimg")
end if

if (inversion) then
!write(0,*) "before grad"
!      call auxscr("grad")
!      call to_history("n1",img%xm%n,"grad"); call to_history("d1",img%xm%d,"grad"); call to_history("o1",img%xm%o,"grad")  
!      call to_history("n2",img%ym%n,"grad"); call to_history("d2",img%ym%d,"grad"); call to_history("o2",img%ym%o,"grad")  
!      call to_history("n3",nz,"grad") ; call to_history("d3",img%z%d,"grad") ; call to_history("o3",img%z%o,"grad")
!      call to_history("n4",max_iters,"grad"); call to_history("d4",1,"grad"); call to_history("o4",1,"grad")
!      call to_history("n5",1,"grad"); call to_history("d5",1,"grad"); call to_history("o5",0,"grad")  
!      call to_history("esize",4,"grad")

write(0,*) "before grad_"
      call auxscr("grad_")
      call to_history("n1",img%xm%n,"grad_"); call to_history("d1",img%xm%d,"grad_"); call to_history("o1",img%xm%o,"grad_")  
      call to_history("n2",img%ym%n,"grad_"); call to_history("d2",img%ym%d,"grad_"); call to_history("o2",img%ym%o,"grad_")  
      call to_history("n3",nz,"grad_") ; call to_history("d3",img%z%d,"grad_") ; call to_history("o3",img%z%o,"grad_")
      call to_history("n4",max_iters,"grad_"); call to_history("d4",1,"grad_"); call to_history("o4",0,"grad_")
      call to_history("n5",1,"grad_"); call to_history("d5",1,"grad_"); call to_history("o5",0,"grad_")  
      call to_history("esize",4,"grad_")

      if (exist_file("scale")) then
      call to_history("n1",img%xm%n,"scale"); call to_history("d1",img%xm%d,"scale"); call to_history("o1",img%xm%o,"scale")  
      call to_history("n2",img%ym%n,"scale"); call to_history("d2",img%ym%d,"scale"); call to_history("o2",img%ym%o,"scale")  
      call to_history("n3",nz,"scale") ; call to_history("d3",img%z%d,"scale") ; call to_history("o3",img%z%o,"scale")
      call to_history("n4",max_iters,"scale"); call to_history("d4",1,"scale"); call to_history("o4",0,"scale")
      call to_history("n5",1,"scale"); call to_history("d5",1,"scale"); call to_history("o5",0,"scale")  
      call to_history("esize",4,"scale")
      end if

write(0,*) "before velout"
      call auxscr("Velout")
      call to_history("n1",img%xm%n,"Velout"); call to_history("d1",img%xm%d,"Velout"); call to_history("o1",img%xm%o,"Velout")  
      call to_history("n2",img%ym%n,"Velout"); call to_history("d2",img%ym%d,"Velout"); call to_history("o2",img%ym%o,"Velout")  
      call to_history("n3",nzold,"Velout") ; call to_history("d3",img%z%d,"Velout") ; call to_history("o3",0,"Velout")
      call to_history("n4",max_iters,"Velout"); call to_history("d4",1,"Velout"); call to_history("o4",1,"Velout")
      call to_history("n5",1,"Velout"); call to_history("d5",1,"Velout"); call to_history("o5",0,"Velout")  
      call to_history("esize",4,"Velout")

!      call auxscr("vtemp")
!      call to_history("n1",img%xm%n,"vtemp"); call to_history("d1",img%xm%d,"vtemp"); call to_history("o1",img%xm%o,"vtemp")  
!      call to_history("n2",img%ym%n,"vtemp"); call to_history("d2",img%ym%d,"vtemp"); call to_history("o2",img%ym%o,"vtemp")  
!      call to_history("n3",nzold,"vtemp") ; call to_history("d3",img%z%d,"vtemp") ; call to_history("o3",0,"vtemp")
!      call to_history("n4",1,"vtemp"); call to_history("d4",1,"vtemp"); call to_history("o4",1,"vtemp")
!      call to_history("n5",1,"vtemp"); call to_history("d5",1,"vtemp"); call to_history("o5",0,"vtemp")  
!      call to_history("esize",4,"vtemp")

!      call auxscr("vreff")
!      call to_history("n1",nref,"vreff"); call to_history("d1",1.,"vreff"); call to_history("o1",1.,"vreff")  
!      call to_history("n2",2,"vreff"); call to_history("d2",img%ym%d,"vreff"); call to_history("o2",img%ym%o,"vreff")  
!      call to_history("n3",img%z%n,"vreff") ; call to_history("d3",img%z%d,"vreff") ; call to_history("o3",zmin,"vreff")
!      call to_history("n4",1,"vreff"); call to_history("d4",1,"vreff"); call to_history("o4",1,"vreff")
!      call to_history("n5",1,"vreff"); call to_history("d5",1,"vreff"); call to_history("o5",0,"vreff")  
!      call to_history("esize",4,"vreff")

!write(0,*) "before dimage"
!      call auxscr("Dimage")
!      call to_history("n1",img%xm%n,"Dimage"); call to_history("d1",img%xm%d,"Dimage"); call to_history("o1",img%xm%o,"Dimage")  
!      call to_history("n2",img%shx%n,"Dimage") ; call to_history("d2",img%shx%d,"Dimage") ; call to_history("o2",img%shx%o,"Dimage")
!      call to_history("n3",img%ym%n,"Dimage"); call to_history("d3",img%ym%d,"Dimage"); call to_history("o3",img%ym%o,"Dimage")  
!      call to_history("n4",img%shy%n,"Dimage") ; call to_history("d4",img%shy%d,"Dimage") ; call to_history("o4",img%shy%o,"Dimage")
!      call to_history("n5",nz,"Dimage") ; call to_history("d5",img%z%d,"Dimage") ; call to_history("o5",zmin,"Dimage")
!      call to_history("esize",4,"Dimage")

!write(0,*) "before bkgimg"
!      call auxscr("bkgimg")
!      call to_history("n1",img%xm%n,"bkgimg"); call to_history("d1",img%xm%d,"bkgimg"); call to_history("o1",img%xm%o,"bkgimg")  
!      call to_history("n2",img%shx%n,"bkgimg") ; call to_history("d2",img%shx%d,"bkgimg") ; call to_history("o2",img%shx%o,"bkgimg")
!      call to_history("n3",img%ym%n,"bkgimg"); call to_history("d3",img%ym%d,"bkgimg"); call to_history("o3",img%ym%o,"bkgimg")  
!      call to_history("n4",img%shy%n,"bkgimg") ; call to_history("d4",img%shy%d,"bkgimg") ; call to_history("o4",img%shy%o,"bkgimg")
!      call to_history("n5",nz,"bkgimg") ; call to_history("d5",img%z%d,"bkgimg") ; call to_history("o5",zmin,"bkgimg")
!      call to_history("esize",4,"bkgimg")

write(0,*) "before dimg"
   call auxscr("dimg")
   call to_history("n1",img%xm%n,"dimg"); call to_history("d1",img%xm%d,"dimg"); call to_history("o1",img%xm%o,"dimg")  
   call to_history("n2",img%shx%n,"dimg") ; call to_history("d2",img%shx%d,"dimg") ; call to_history("o2",img%shx%o,"dimg")
   call to_history("n3",img%ym%n,"dimg"); call to_history("d3",img%ym%d,"dimg"); call to_history("o3",img%ym%o,"dimg")  
   call to_history("n4",img%shy%n,"dimg") ; call to_history("d4",img%shy%d,"dimg") ; call to_history("o4",img%shy%o,"dimg")
   call to_history("n5",nz,"dimg") ; call to_history("d5",img%z%d,"dimg") ; call to_history("o5",zmin,"dimg")
   call to_history("n6",max_iters,"dimg") ; call to_history("d6",1,"dimg") ; call to_history("o6",1,"dimg")
   call to_history("esize",4,"dimg")

write(0,*) "before bimg"
   call auxscr("bimg")
   call to_history("n1",img%xm%n,"bimg"); call to_history("d1",img%xm%d,"bimg"); call to_history("o1",img%xm%o,"bimg")  
   call to_history("n2",img%shx%n,"bimg") ; call to_history("d2",img%shx%d,"bimg") ; call to_history("o2",img%shx%o,"bimg")
   call to_history("n3",img%ym%n,"bimg"); call to_history("d3",img%ym%d,"bimg"); call to_history("o3",img%ym%o,"bimg")  
   call to_history("n4",img%shy%n,"bimg") ; call to_history("d4",img%shy%d,"bimg") ; call to_history("o4",img%shy%o,"bimg")
   call to_history("n5",nz,"bimg") ; call to_history("d5",img%z%d,"bimg") ; call to_history("o5",zmin,"bimg")
   call to_history("n6",max_iters,"bimg") ; call to_history("d6",1,"bimg") ; call to_history("o6",1,"bimg")
   call to_history("esize",4,"bimg")
end if

!--------------------------------------------------
! Frequency vector
! no zero frequency component
!
rec%w%o   = 2*pi * rec%w%o
rec%w%d   = 2*pi * rec%w%d
allocate(rec%w%all(rec%w%n))
do iwb = 1,rec%w%n
   rec%w%all(iwb) = rec%w%o + (iwb - 1) * rec%w%d
end do

sou%w%o   = 2*pi * sou%w%o
sou%w%d   = 2*pi * sou%w%d
allocate(sou%w%all(sou%w%n))
do iwb = 1,sou%w%n
   sou%w%all(iwb) = sou%w%o + (iwb - 1) * sou%w%d
end do

!=====================
! Defining blocks
! w and z
call from_param("nwblk",wblock%ns,node)
call from_param("nzblk",zblock%ns,int(nz/3))
call base_blkcal(wblock, sou%w)
call base_blkcal(zblock, img%z)

!-------->
write(0,*) "Frequency blocking"
write(0,*) wblock%bres,wblock%ds,wblock%nb,wblock%nb_res
write(0,*) wblock%blen(:)
write(0,*) wblock%os(:)
write(0,*) wblock%indx(:)
!-------->
write(0,*) "Depth blocking"
write(0,*) zblock%bres,zblock%ds,zblock%nb,zblock%nb_res
write(0,*) zblock%blen(:)
write(0,*) zblock%os(:)
write(0,*) zblock%indx(:)
!----------------------------------------------
! Set bounds
 call from_param("bounds",bounds,.false.)
 if (bounds) then
    call from_param("cbounds",cbounds,.false.)
    if (cbounds) then
       call from_param("vmin",vmin)
       call from_param("vmax",vmax)
    else
       call from_param("p_up",p_up)
       call from_param("p_lo",p_lo)
       if (p_up >= 1.) p_up=1+p_up/100
       if (p_lo >= 1.) p_lo=1-p_lo/100
    end if
 end if 
!----------------------------------------------
!!!!!!!!!!!!!!!!!!
   ! Horizon with no velocity update
   if (exist_file("horv0")) then
      call from_param('horv0',hor1)
   end if

   ! ********************************************** 
   !  Initialization of code variables 
    iters=1
    err=1
    epsteep=1e-8 
    tol(1)=1e-8
    tol(2)=1e-5

   ! Line search variables 
    factor=0.5
    lambda_min=0.00001
   ! ********************************************** 
   !   Set tolerance for stopping criterium 
    err1=tiny(1.)
   ! ********************************************** 
    tola=0; 
   ! ********************************************** 

!<--------
   ninx=tapertrx+padtracesx+1
   nfnx=rec%hx%n+padtracesx+tapertrx
   niny=tapertry+padtracesy+1
   nfny=rec%hy%n+padtracesy+tapertry
!<--------
!----------------------------------------------------
write(0,*) "   PARM   "
write(0,*) "=========="
write(0,*) "PARAMETERS"
write(0,*) "=========="
write(0,*) //
write(0,*) "no.ref.velocities   = ",nref
write(0,*) "min.x-coord         = ",xmin
write(0,*) "max.x-coord         = ",xmax
write(0,*) "min.y-coord         = ",ymin
write(0,*) "max.y-coord         = ",ymax
write(0,*) "min.depth           = ",zmin
write(0,*) "max.depth           = ",zmax
write(0,*) "traces in taper x   = ",tapertrx
write(0,*) "padded traces   x   = ",padtracesx
write(0,*) "traces in taper y   = ",tapertry
write(0,*) "padded traces   y   = ",padtracesy
write(0,*) //
write(0,*) "============"
write(0,*) "HEADER WORDS"
write(0,*) "============"
write(0,*) //
write(0,*) "SHOTS"
write(0,*) "====="
write(0,*) " # xlines  = ", rec%hx%n
write(0,*) "xline intv = ", rec%hx%d
write(0,*) "xline orig = ", rec%hx%o
write(0,*) " # ilines  = ", rec%hy%n
write(0,*) "iline intv = ", rec%hy%d
write(0,*) "iline orig = ", rec%hy%o
write(0,*) "number of freq  = ", rec%w%n
write(0,*) "freq spacing    = ", rec%w%d/(2*pi)
write(0,*) "minimum freq    = ", rec%w%o/(2*pi)
write(0,*) "number of shots = ", rec%xd%n
write(0,*) "shot spacing    = ", rec%xd%d
write(0,*) "min.shot-x      = ", rec%xd%o
write(0,*) //
write(0,*) "SOURCE FUNCTION"
write(0,*) "==============="
write(0,*) " # xlines  = ", sou%hx%n
write(0,*) "xline intv = ", sou%hx%d
write(0,*) "xline orig = ", sou%hx%o
write(0,*) " # ilines  = ", sou%hy%n
write(0,*) "iline intv = ", sou%hy%d
write(0,*) "iline orig = ", sou%hy%o
write(0,*) "number of freq  = ", sou%w%n
write(0,*) "freq spacing    = ", sou%w%d/(2*pi)
write(0,*) "minimum freq    = ", sou%w%o/(2*pi)
write(0,*) "number of shots = ", sou%xd%n
write(0,*) "shot spacing    = ", sou%xd%d
write(0,*) "min.shot-x      = ", sou%xd%o
write(0,*) //
write(0,*) "VELOCITY MODEL"
write(0,*) "=============="
write(0,*) " # xlines   = ", v%xm%n
write(0,*) " xline intv = ", v%xm%d
write(0,*) " xline orig = ", v%xm%o
write(0,*) " # ilines   = ", v%ym%n
write(0,*) " iline intv = ", v%ym%d
write(0,*) " iline orig = ", v%ym%o
write(0,*) "no.depth samples= ", v%z%n
write(0,*) "depth interval  = ", v%z%d
write(0,*) "minimum depth   = ", v%z%o
write(0,*) //
write(0,*) "IMAGE SPACE"
write(0,*) "==========="
write(0,*) " # xlines   = ", img%xm%n
write(0,*) " xline intv = ", img%xm%d
write(0,*) " xline orig = ", img%xm%o
write(0,*) " # ilines   = ", img%ym%n
write(0,*) " iline intv = ", img%ym%d
write(0,*) " iline orig = ", img%ym%o
write(0,*) "no.depth samples= ", img%z%n
write(0,*) "depth interval  = ", img%z%d
write(0,*) "minimum depth   = ", img%z%o
write(0,*) "no.of x-offsets    = ", img%shx%n
write(0,*) "x-offset spacing   = ", img%shx%d
write(0,*) "x-offset origin     = ", img%shx%o
write(0,*) "no.of y-offsets    = ", img%shy%n
write(0,*) "y-offset spacing   = ", img%shy%d
write(0,*) "y-offset origin     = ", img%shy%o

write(0,*) //
write(0,*) "WAVENUMBER"
write(0,*) "=========="
write(0,*) "no.wavenumbers x = ", kxm%n
write(0,*) "waveno.spacing x = ", kxm%d
write(0,*) "waveno.origin  x = ", kxm%o
write(0,*) "no.wavenumbers y = ", kym%n
write(0,*) "waveno.spacing y = ", kym%d
write(0,*) "waveno.origin  y = ", kym%o
write(0,*) //

write(0,*) "*****************************************************"
end if
end subroutine

end module
