program OW_rhopanel
! rHO PANEL FOR VELOCITY PICKING

  use sep
  use OW_type
  use OW_parms_rho
  use OW_halfder
  use OW_rho
  use solver_reg_mod
  use solver_reg_mod1
  use solver_prc_mod
  use igrad1
  use causint
  use cgstep_mod

  implicit none

  real, allocatable, dimension (:)		:: in,out,in1,out1
  real, allocatable, dimension (:,:)              :: dips
  integer					:: stat,iy

  call sep_init()
  call sep_begin_prog()
  call owrhopar()
  allocate(in(img%z%n*img%ah%n*img%xm%n));in=0
  allocate(out(img%z%n*img%rho%n*img%xm%n));out=0
  if (op(1:3)=="dot") then
     allocate(in1(img%z%n*img%ah%n*img%xm%n))
     allocate(out1(img%z%n*img%rho%n*img%xm%n))
     call random_number(in1)
     in1=(2*in1-1)/1000
     out1=cos(in1)
     stat = ow_rho_lop(.true.,.false.,out,in1)
     stat = ow_rho_lop(.false.,.false.,out1,in)
     write(0,*) dot_product(in,in1),dot_product(out,out1)
     deallocate(in1,out1)
  else if (op(1:3)=="mod") then
     call sreed("in",out,4*img%z%n*img%rho%n*img%xm%n)
!     call halfdifa(out)
     stat = ow_rho_lop(.false.,.false.,out,in)
     call srite("out",in,4*img%z%n*img%ah%n*img%xm%n)
  else if (op(1:3)=="opt") then
     do iy=1,img%ym%n
        call sreed("in",in,4*img%z%n*img%ah%n*img%xm%n)
        if (geodip) then
           allocate(dips(img%z%n,img%xm%n))
           call sreed("dips",dips,4*img%z%n*img%xm%n)
        end if
        call halfdifa(in)
        if (method(1:3)=="reg") then
           call solver_reg(m=out,d=in, Fop=ow_rho_lop, Aop=igrad1_lop, stepper=cgstep, &
                nAop=img%z%n*img%rho%n*img%xm%n, niter=niter,eps1=eps, verb=.true.)
        else if (method(1:3)=="gcg") then
           call solver_reg1(m=out,d=in, Fop=ow_rho_lop, Aop=igrad1_lop, stepper=cgstep, &
                nAop=img%z%n*img%rho%n*img%xm%n, niter=niter,eps1=eps, verb=.true.)
        else if (method(1:3)=="prc") then
           call solver_prc(m=out,d=in, Fop=ow_rho_lop, Sop=causint_lop, stepper=cgstep, &
                nSop=img%z%n*img%rho%n*img%xm%n, niter=niter,eps1=eps, verb=.true.)
        end if
        call srite("out",out,4*img%z%n*img%rho%n*img%xm%n)
        if (exist_file("dips")) deallocate(dips)
     end do
  end if
  deallocate(in,out)
  call sep_end_prog()
end program 
