program OW_3dip
! Determine dip by linear Radon transform

  use sep
  use OW_type
  use OW_sltpar
  use OW_3diprad
  use OW_io
  use solver_tiny_mod
  use solver_prc_mod 
  use igrad1
  use cgstep_mod

  implicit none

  real, allocatable, dimension (:)		:: in,out,dip,temp,mask,curv
  integer					:: stat,ix,ih,iz,ind(1)

  call sep_init()
  call sep_begin_prog()
  call owsltpar()
  allocate(in(sog%z%n*sog%xm%n))

open(unit=2,file="lixo")

  call sreed("in",in,4*sog%z%n*sog%xm%n)

  if (exist_file("mask")) then
     allocate(mask(sog%z%n*sog%xm%n))
     call sreed("mask",mask,4*sog%z%n*sog%xm%n)
     where(mask >= maskval)
        mask=1
     elsewhere
        mask=0.
     end where
  end if

  if (plunge) then
     if (proc) then
        allocate(out(sog%z%n*nax*sog%xm%n))
        stat = ow3diprad(.true.,.false.,out,in)
        if (test) then
           stat = ow3diprad(.false.,.false.,out,in)
           call owio(in,0)
        else if (initdip) then
           allocate(dip(sog%z%n*sog%xm%n))
           allocate(temp(nax))
           do ix=1,sog%xm%n
              do iz=1,sog%z%n
                 do ih=1,nax
                    if (exist_file("mask")) then
                       temp(ih)=out((ix-1)*sog%z%n*nax+(iz-1)+(ih-1)*sog%z%n+1)*mask((ix-1)*sog%z%n+iz)
                    else
                       temp(ih)=out((ix-1)*sog%z%n*nax+(iz-1)+(ih-1)*sog%z%n+1)
                    end if
                 end do
                 ind=maxloc(abs(temp))
                 if(ind(1)==1) ind(1)=nax/2+1
                 dip((ix-1)*sog%z%n+iz)=(ind(1)-1)*dax+oax
              end do
           end do
           call owio(dip,1)
           deallocate(dip,temp)
        else
           call owio(out,2)
        end if
     else
        allocate(out(sog%z%n*sog%xm%n))
        if (tin)  call solver_tiny(m=out, d=in, Fop=ow3diprad, stepper=cgstep, niter=niter)
        if (prc)   call solver_prc( m=out,d=in, Fop=ow3diprad, Sop=igrad1_lop, stepper=cgstep, &
                                    nSop=size(out), niter=niter,eps=eps, verb=.true.)
         call owio(out,1)
     end if
  else if (fold) then
     allocate(dip(sog%z%n*sog%xm%n))
     call sreed("dip",dip,4*sog%z%n*sog%xm%n)
     if (proc) then
        allocate(out(sog%z%n*nrad*sog%xm%n))
        stat = ow3curvrad_init(dip)
        stat = ow3curvrad(.true.,.false.,out,in)
        stat = ow3curvrad_close()
        if (test) then
           stat = ow3curvrad_init(dip)
           stat = ow3curvrad(.false.,.false.,out,in)
           stat = ow3curvrad_close()
           call owio(in,3)
        else if (initcurv) then
           allocate(curv(sog%z%n*sog%xm%n))
           allocate(temp(nrad))
           do ix=1,sog%xm%n
              do iz=1,sog%z%n
                 do ih=1,nrad
                    if (exist_file("mask")) then
                       temp(ih)=out((ix-1)*sog%z%n*nrad+(iz-1)+(ih-1)*sog%z%n+1)*mask((ix-1)*sog%z%n+iz)
                    else
                       temp(ih)=out((ix-1)*sog%z%n*nrad+(iz-1)+(ih-1)*sog%z%n+1)
                    end if
                 end do
                 ind=maxloc(abs(temp))
                 if(ind(1)==1) ind(1)=nax/2
                 curv((ix-1)*sog%z%n+iz)=(ind(1)-1)*dax+oax
              end do
           end do
           call owio(curv,4)
           deallocate(curv,temp)
        else
           call owio(out,5)
        end if
     else
!        allocate(out(sog%z%n*sog%xm%n))
!        if (tin)  call solver_tiny(m=out, d=in, Fop=ow3diprad, stepper=cgstep, niter=niter)
!        if (prc)  call solver_prc( m=out,d=in, Fop=ow3diprad, Sop=igrad1_lop, stepper=cgstep, &
!                                    nSop=size(out), niter=niter,eps=eps, verb=.true.)
!         call owio(out,1)
     end if
  end if
  deallocate(in,out)
  if (exist_file("mask")) then
     deallocate(mask)
  end if

close(2)
  call sep_end_prog()
end program 
