program RWE2D_FM
  use sep
  use util_
  use RWE2D_WI_wemig
  use C2R2D
  use sep_timers_mod
  implicit none
  !! ax is the Cartesian length dimension associated with Velocity 
  !! az is the Cartesian depth  dimension associated with Velocity 
  !! aw is the frequency axis
  !! as is the shot axis
  !! ag is the offset axis in Ray space
  !! at is the extrapolation direction in Ray space
  !! ar is the number of references
  type(myaxis) :: ax,az,aw,ar,as,at,ag

  integer :: ix,iz,ii,iw2,ishot(1),iis,mxl,mxr,oldn,ig,it,numit,in,nw,izcut(1),ixcut(1),ib
  integer :: nref,nsx,nsz,norm,cnx,cnz,migapp,ptap,nzcut,nxcut,nblock
  real    :: forward,xmin,zmin,xmax,zmax,dxx,dzz,omin,omax,pmax,zcut,jscale,njs
  real    :: rms,rmsold,www,r0
  real,    allocatable,dimension(:,:)  :: Vel,gnorm,Velin,dVel,Vtmp,CDnew,CDold,dVelold
  real,    allocatable,dimension(:,:)  :: gradient,Steplength,gtmp,Dscale,Sfilter
  real,    pointer    ,dimension(:,:)  :: HessScale,Htmp
  real,    allocatable,dimension(:,:,:):: raysin
  real,    allocatable,dimension(:,:):: Dwfld
  complex, pointer    ,dimension(:,:)  :: data,swf,model,resid,inject
  complex, allocatable,dimension(:,:,:):: swftmp,datatmp,wfld,SGF,Hwfld

  logical                              :: verbose,kin,adj,wf,is1D,usecg,forwardmodel
  character(len=128) :: name

  !! Geometry
  integer, allocatable :: mask(:,:)
  real,    allocatable :: refs(:,:,:), fields(:,:,:)

  !! Timers
  integer :: t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15
  logical :: logic


  !call sep_init(SOURCE)
  call sep_init()
  call getch_add_string("head=/dev/null")
  call getch_add_string("noheader=y")
  call set_no_putch()
  call sep_begin_prog()

  write(0,*) 'STARTING PROGRAM'

  call from_param("forward",forward,0.) !! Forward scattering option
  call from_param("nref",nref,30)        !! # reference velocities
  call from_param("verbose",verbose,.false.)
  call from_param("kinematic",  kin,.true.)
  call from_param("norm",norm,1)
  call from_param("nsx",nsx,3)
  call from_param("nsz",nsz,3)
  call from_param("omin",omin,2000.)
  call from_param("omax",omax,2400.)
  call from_param("wf",wf,.false.)
  call from_param("pmax",pmax,0.01)
  call from_param("numit",numit,1)
  call from_param("ptap",ptap,10)
  call from_param("zcut",zcut,0.)
  call from_param("jscale",jscale,1.)
  call from_param("is1D",is1D,.true.)
  call from_param("usecg",usecg,.true.)
  call from_param("nblock",nblock,1)
  call from_param("forwardmodel",forwardmodel,.true.)
  !! Read in input dimensions
  name="swf"
  call ainaxis(ag,1,name) 
  call ainaxis(aw,2,name) 
  call ainaxis(as,3,name) 

  name="rays"
  call ainaxis(at,1,name)

  !!    Read in velocity file
  name="vel"
  call ainaxis(az,1,name)
  call ainaxis(ax,2,name)
  !! Do only 1 frequency
  nw = aw%n
  aw%n=1


  write(0,*) 'DEPTH',at%n,at%o,at%d
  write(0,*) 'H-POS',ag%n,ag%o,ag%d
  write(0,*) 'SHOTS',as%n,as%o,as%d
  write(0,*) 'FREQ ',aw%n,aw%o,aw%d
  write(0,*) 'VEL X',ax%n,ax%o,ax%d
  write(0,*) 'VEL Z',az%n,az%o,az%d


  !! . . Create offset axis
  call aouaxis(az,1,"image")
  call aouaxis(ax,2,"image")
  call auxputch("esize","i",4,"image")
  call auxputch("n3","i",nw,"image")
  call auxputch("d3","r",1.,"image")
  call auxputch("o3","r",0.,"image")

  !  call aouaxis(ag,1,"DH.H")
  !  call aouaxis(at,2,"DH.H")!

  !  call aouaxis(at,1,"gnorm.H")
  !  call aouaxis(ag,2,"gnorm.H")

  !  call auxputch("esize","i",8,"DH.H")
  !  call auxputch("n3","i",2,"DH.H")
  !  call auxputch("d3","r",1.,"DH.H")
  !  call auxputch("o3","r",0.,"DH.H")

  call aouaxis(az,1,"HessScale.H")
  call aouaxis(ax,2,"HessScale.H")
  call aouaxis(az,1,"Steplength.H")
  call aouaxis(ax,2,"Steplength.H")
  call auxputch("n3","i",numit,"Steplength.H")
  call auxputch("n4","i",nw,"Steplength.H")

  call aouaxis(az,1,"CDnew.H")
  call aouaxis(ax,2,"CDnew.H")
  call auxputch("n3","i",numit,"CDnew.H")
  call auxputch("n4","i",nw,"CDnew.H")

  call aouaxis(az,1,"Dvel.H")
  call aouaxis(ax,2,"Dvel.H")
  call auxputch("n3","i",numit,"Dvel.H")
  call auxputch("n4","i",nw,"Dvel.H")

  !! . . Output Imaged WF in RWE coordinates
  call aouaxis(at,1,"Rimage")
  call aouaxis(ag,2,"Rimage")
  call auxputch("n3","i",numit,"Rimage")
  call auxputch("n4","i",nw,"Rimage")
  call auxputch("esize","i",4,"Rimage")

  call aouaxis(at,1,"Dwfld.H")
  call aouaxis(ag,2,"Dwfld.H")
  call auxputch("n3","i",numit,"Dwfld.H")
  call auxputch("n4","i",nw   ,"Dwfld.H")

  call aouaxis(at,1,"Sfilter.H")
  call aouaxis(ag,2,"Sfilter.H")
  call auxputch("n3","i",numit,"Sfilter.H")
  call auxputch("n4","i",nw   ,"Sfilter.H")


  call aouaxis(at,1,"Perturb.H")
  call aouaxis(aw,2,"Perturb.H")
  call auxputch("n3","i",3 ,"Perturb.H")
  call auxputch("o3","r",0.,"Perturb.H")
  call auxputch("d3","r",1.,"Perturb.H")
  call auxputch("n4","i",as%n,"Perturb.H")
  call auxputch("n5","i",numit,"Perturb.H")
  call auxputch("n6","i",1+nw/nblock, "Perturb.H")

  call auxputch("esize","i",8,"Perturb.H")
  call auxputch("n3","i",nw,"wfld1.H")
  !! Input data files
  allocate( swftmp(ag%n,nw,as%n),datatmp(at%n,nw,as%n),Dscale(nw,as%n) )

  !! Along data line
  allocate(  model(at%n,nblock),data(at%n,nblock),inject(at%n,nblock),resid(at%n,nblock) )

  !! Extrapolation
  allocate( swf(ag%n,nblock) )

  !! Velocity and Step
  allocate( dVel(az%n,ax%n),Vtmp(az%n,ax%n),CDold(az%n,ax%n),CDnew(az%n,ax%n) )
  allocate( dVelold(az%n,ax%n),gtmp(az%n,ax%n),Sfilter(at%n,ag%n) )
  allocate( Dwfld(at%n,ag%n),HessScale(az%n,ax%n),Htmp(az%n,ax%n) )
  allocate( Vel(at%n,ag%n),Velin(az%n,ax%n),gradient(at%n,ag%n),Steplength(az%n,ax%n) )

  !! Wavefield
  allocate( Hwfld(ag%n,at%n,nblock),wfld(ag%n,at%n,nblock),SGF(ag%n,at%n,nblock))

  !! Geometry
  allocate( mask(ag%n ,at%n), refs(nref,6,at%n),fields(ag%n,6,at%n) )
  allocate( raysin(at%n,ag%n,2),gnorm(at%n,ag%n) )

  !! . . Read in Rays and Velocity
  call sep_read(raysin,"rays")
  call sep_read(Velin ,"vel ")
  write(0,*) 'VELOCITY MIN/MAX', minval(Velin),maxval(Velin)

  !! Where to cut
  izcut = minloc(abs(raysin(1,:,2) - zcut))
  nzcut= izcut(1)
  write(0,*) 'CUTTING AT ',nzcut,' of ', ag%n 

  !! . . Initialize Routines
  call geometry_init(at%n,ag%n,at%d,ag%d,kin,nref,verbose)
  call wemig_init(ag,at,aw,ar,forward,kin,nref,verbose,ptap,nzcut,nblock)
  call C2R_init(0,0,norm,ax%d,az%d,ax%n,az%n,ag%n,at%n,ax%o,az%o)
  call wftaper_init(omin,omax)

  !! read in Source and data wavefields
  swftmp=0.
  call sep_read(swftmp ,"swf")
  call sep_read(datatmp,"data")
!  write(0,*) swftmp(1:3,1:3,1)
!  do iw2=1,nw
!     do ii=1,as%n
!        write(0,*) maxval(abs(swftmp(:,iw2,ii))),maxval(abs(datatmp(:,iw2,ii)))
!        Dscale(iw2,ii) = jscale
!        datatmp(:,iw2,ii)=datatmp(:,iw2,ii)/&
!        maxval(abs(datatmp(:,iw2,ii)))/sqrt( dble(iw2-1)*aw%d+aw%d)
!        swftmp(:,iw2,ii)=  swftmp(:,iw2,ii)/&
!        maxval(abs( swftmp(:,iw2,ii)))/sqrt( dble(iw2-1)*aw%d+aw%d)
!        !  swftmp(:,iw2,ii)=jscale* swftmp(:,iw2,ii)/(abs(datatmp(1,iw2,ii)))
!        write(0,*) maxval(abs(swftmp(:,iw2,1))),maxval(abs(datatmp(:,iw2,1)))
!     end do
!  end do
!  datatmp = datatmp*jscale
!  swftmp =  swftmp*jscale
  CDold=0.;CDnew=0.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !! . . MAIN PROGRAM
  !!
  njs = jscale
  logic=init_sep_timers()
  logic=setup_next_timer("Frequencies",t1)
  logic=setup_next_timer("NL Interations",t2)
  logic=setup_next_timer("Shots",t3)
  logic=setup_next_timer("Vel update",t4)
  logic=setup_next_timer("StepCalc",t5)
  logic=setup_next_timer("geometry",t6)
  logic=setup_next_timer("Interp1",t7)
  logic=setup_next_timer("DataAssign",t8)
  logic=setup_next_timer("FP1",t9)
  logic=setup_next_timer("Reshape",t10)
  logic=setup_next_timer("Residual",t11)
  logic=setup_next_timer("BPResid",t12)
  logic=setup_next_timer("DVelcalc",t13)
  logic=setup_next_timer("Interp2",t14)
  logic=setup_next_timer("Output",t15)

  call C2R_run(Vel,raysin,Velin,gnorm,.false.)
  write(0,*) maxval(Vel)  ,minval(Vel)
  write(0,*) maxval(Velin),minval(Velin)

  call geometry_calc(raysin,Vel,fields,refs,mask,gnorm)

  do iw2=1,nw,nblock !! . . Loop over input frequencies

     www = 2.*3.1415926*( (iw2-1.)*aw%d+aw%o )

     !! . . Initalization and info
     swf = 0.;wfld=0.;inject=0.
     swf  =  swftmp(:,iw2:iw2+nblock,1)   

     call wemig(swf,inject,wfld,.true.,1,fields,refs,mask,gnorm,iw2)

  end do

  call sep_end_prog()
  call exit()
end program RWE2D_FM
