module geodip1
!
!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 inv_tomo_mod
  use filters_mod
  use i_op_mod

  implicit none
!-----------------------------------------------------------
subroutine OW_geodip(g)

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

!MASK
  real,allocatable                      :: mask(:),mask_in(:,:)

  call owsteer_parms()
  slow0(frame+1:frame+model_param%nz_mod,:,:)=transpose(g)
  do i1=1,frame
    slow0(i1,:,:)=slow0(frame+1,:,:)
    slow0(frame+img%z%n+i1,:,:)=slow0(frame+img%z%n,:,:)
  end do
  newz=img%z%n+frame*2
!############
!ARRAY SIZES#
!############
  nelem_model=model_param%ndim_mod*img%z%n*img%xm%n*img%ym%n
!###############
!INIT OPERATORS#
!###############
  write(0,*) "initialize operators"
  call helix_chain_init(nelem_model)
  call polydiv_init(filt)
!#################
!RUN THE OPERATOR#
!#################
  if (exist_file("mask")) then
     call from_aux("mask","n2",temp2i)
     call from_aux("mask","n1",tempi)
     if (tempi+frame*2.ne.img%z%n .or. temp2i.ne.img%xm%n) then
        write(0,*) "mask",tempi,temp2i
        write(0,*) "model",img%z%n,img%xm%n
        call erexit("angles not the same size as model")
     end if
     allocate(mask((tempi+frame*2)*temp2i),mask_in(tempi,temp2i))
     call sreed('mask',mask_in,tempi*temp2i*4)
     do i2=1,temp2i
       do i1=1,tempi
         mask(i1+frame+(i2-1)*(2*frame+tempi))=mask_in(i1,i2)
       end do
       do i1=1,frame
         mask(i1+(i2-1)*(2*frame+tempi))=mask(frame+1+(i2-1)*(2*frame+tempi))
       end do
       do i1=1,frame
          mask(i1+frame+tempi+(i2-1)*(2*frame+tempi))= mask(frame+tempi+(i2-1)*(2*frame+tempi))
       end do
     end do
     deallocate(mask_in)
     call mask_init(mask)
  end if
  call inv_tomo (param,model_param,nelem_model,slow,slow0,i_op,epspoly)
  if(allocated(mask)) then
    deallocate(mask)
    call mask_close()
  end if
!#################
!WRITE RESULTS   #
!#################
  newz=newz-frame*2
!#################
! Change when 3D #
!#################
  a1=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=a1/b1*transpose(slow(frame+1:newz+frame,:,1))

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