module geodip
!
!infile=grd.H outfile=grdflt.H angles=xhh.H steer_lags=15 niter=20 epspoly=20. tri_width=6350
!infile=grd.H outfile=grdflt.H angles=xhh.H steer_lags=40 15 niter=20 epspoly=40. tri_width=130000
!infile=grd.H outfile=grdflt.H angles=xhh.H steer_lags=30 niter=12 epspoly=80. tri_width=130000
!infile=grd.H outfile=grdflt.H angles=xhh.H mask=1.H steer_lags=40 20 niter=12 epspoly=160. tri_width=1000000
!
  use OW_type
  use OW_parms_areal
  use inv_tomo_mod
  use filters_mod
  use i_op_mod
  use process_types
  use OW_grabvel_areal
  use OW_steer_parms

  implicit none
  real,allocatable                             :: pick(:)
  contains
!-----------------------------------------------------------
 subroutine owgeodip_init(pick_in)
   real,dimension(:)                             :: pick_in
   allocate(pick(size(pick_in)))
   pick=pick_in
   pick=pick+float(frame)
   write(0,*) "max pick", int(maxval(pick)),"min pick", int(minval(pick))
 end subroutine
!-----------------------------------------------------------
subroutine OW_geodip(g)

  double precision                      :: g(img%xm%n,img%z%n)
  real                                  :: slow0(img%z%n+frame*2,img%xm%n,1)!img%ym%n)
  real                                  :: a1,b1,slow(img%z%n+frame*2,img%xm%n,1)!img%ym%n)
  integer                               :: ierr,nelem_ray,nelem_model,verb,newz
  integer                               :: n_h,n_mid,auxout,i1,i2,tempi,temp2i,temp3i,type,test_oper

!MASK
  real,allocatable                      :: mask(:)

  call owsteer_parms()
  newz=img%z%n+frame*2
  slow0(frame+1:frame+model_param%nz_mod,:,1)=transpose(sngl(g))

  do i1=1,frame
    slow0(i1,:,1)=slow0(frame+1,:,1)
    slow0(frame+img%z%n+i1,:,1)=slow0(frame+img%z%n,:,1)
  end do
!############
!ARRAY SIZES#
!############
  nelem_model=model_param%ndim_mod*img%z%n*img%xm%n*img%ym%n
!###############
!INIT OPERATORS#
!###############
  call owfilter_init()
  write(0,*) "initialize operators"
!  call helix_chain_init(nelem_model)
  call polydiv_init(filt)
!#################
!RUN THE OPERATOR#
!#################
  if (allocated(pick)) then
     allocate(mask(newz*img%xm%n))
     mask=1
     do i2=1,img%xm%n
        do i1=1,newz
           if (i1 < int(pick(i2))) then
              mask((i2-1)*newz+i1)=0.
           end if
        end do
     end do
     call mask_init(mask)
  end if
  call inv_tomo (slow,slow0,nelem_model)!,i_op)
  if(allocated(mask)) then
    deallocate(mask,pick)
    call mask_close()
  end if
   if (exist_file("maskout")) then
      call to_history("n1",newz,"maskout"); call to_history("d1",img%z%d,"maskout"); call to_history("o1",img%z%o-(frame*img%z%d),"maskout")  
      call to_history("n3",1,"maskout") ; call to_history("d3",1.,"maskout") ; call to_history("o3",0.,"maskout")
      call to_history("n2",img%xm%n,"maskout") ; call to_history("d2",img%xm%d,"maskout") ; call to_history("o2",img%xm%o,"maskout")
      call to_history("n4",1,"maskout"); call to_history("d4",1,"maskout"); call to_history("o4",0,"maskout")
      call to_history("n5",1,"maskout"); call to_history("d5",1,"maskout"); call to_history("o5",0,"maskout")  
      call to_history("esize",4,"maskout")
      call srite("maskout",slow,4*newz*img%xm%n)      
   end if
!#################
!WRITE RESULTS   #
!#################
  newz=newz-frame*2
!#################
! Change when 3D #
!#################
  a1=sngl(sqrt(sum(g*g)/size(g)))
  b1=0
  do i2=1,img%z%n
     do i1=1,img%z%n
        b1=b1+slow(frame+i2,i1,1)*slow(frame+i2,i1,1)
     end do
  end do
  b1=sqrt(b1/(img%z%n*img%xm%n))
  g=dble(a1/b1*transpose(slow(frame+1:newz+frame,:,1)))

  end subroutine

end module

!  allocate(tempb(nelem_model),tempc(nelem_model))
!  tempi=0
!  do i2=1,img%xm%n
!     do i1=1,img%z%n
!        tempi= tempi+1
!        tempb(tempi)=d_slow(i1,i2,1)
!     end do
!  end do
!  ierr=polydiv_lop(.false.,.false.,tempb,tempc)
!  tempi=0
!  do i2=1,img%xm%n
!     do i1=1,img%z%n
!        tempi= tempi+1
!        tempa(i1,i2,1)=tempc(tempi)
!     end do
!  end do
!  slow=slow0+tempa
!  slow=1./slow
!  deallocate(tempa,stat=ierr)
!  if (0.ne. putch('n1','d',img%z%n )) call erexit('Trouble writing  n1 to history  ')
!  if (0.ne. putch('n2','d',img%xm%n )) call erexit('Trouble writing  n2 to history  ')
!  if (0.ne. putch('n3','d',3 )) call erexit('Trouble writing  n3 to history  ')
!   allocate (tempa(2*frame+img%z%n,img%xm%n,img%ym%n))
