program DirectRadon
  use sep
  use solver_tiny_brad
  use cgstep_mod
  use planemod2
  use dottest

  implicit none

  INTEGER, PARAMETER      :: dbl = 4
  logical    :: adj,add,dotme
  integer                          :: nt, nx, n2, niter, stat,its,ntr,ntrr,jniter,previts,st
  integer                          :: n3, j,i,nq
  real                             :: o1,d1, o2,d2, o3, d3, q0,dq, blck
  real(kind=dbl), dimension (:), allocatable :: data, model, recons
  real,dimension(2)::dot1,dot2

  call sep_init ()                 ! start processing parameters
  call from_param  ("adj",adj,.true.)
  call from_history ("n3", n3)
  !  call from_history ("o3", o3)
  !  call from_history ("d3", d3)
  call from_param("dot",dotme,.false.)
  call from_history ("o1", o1)
  call from_history ("d1", d1)
  call from_history ("n1", nt)
  call from_aux('geom',"n2",ntr)
  if(adj) then
     call from_history ("n2", ntrr)
     if(ntrr.ne.ntr) call seperr("wrong number of headers")
  end if
  call from_history ("o2", o2)
  call from_history ("d2", d2)

  nx=ntr
  write(0,*)' ntr,nx',ntr,nx

  write(0,*)' ------- Adj:',adj
  !!direct arrival stuff
  call from_param  ("nq", nq, 150)
  call from_param  ("q0", q0, 0.000645161)
  call from_param  ("dq", dq, 0.000000574)
  !----------------------------------------------------------------

  if (adj) then 
     call to_history("label2","slowness")
     call to_history("nq", nq)
     call to_history("q0", q0) 
     call to_history("dq", dq)  
     n2=nq; o2=q0 ; d2=dq
  else
     !     write(0,*) '>> BE SURE THAT HRT MODEL SPACE IS FIRST '
     call to_history ("label2","trace")
     call from_aux("geom","o2", o2)
     call from_aux("geom","n2", n2)
     call from_aux("geom","d2", d2)
  end if


  !!HRT model space
  call to_history ("o2", o2)
  call to_history ("d2", d2)
  call to_history ("n2", n2)
  call from_param ("niter", niter, 10) 
  call from_param ('jniter',jniter,niter)
  if (exist_file ("recons")) then
     call auxout("recons")
     call to_history ("o1", o1,"recons")
     call to_history ("d1", d1,"recons")
     call to_history ("n1", nt,"recons")
     call to_history ("label1","time(s)","recons")

     call to_history ("o2", 0.,"recons")
     call to_history ("d2", 1.,"recons")
     call to_history ("n2", ntr,"recons")

     call to_history ("o3", 0.,"recons")
     call to_history ("d3", 1.,"recons")
     call to_history ("n3", n3,"recons")
     call to_history ("label2","traces","recons")
     call to_history("n4",1,"recons")
  end if
  write(0,*)' nt,nx,nq',nt,nx,nq
  allocate (  data(nt*nx),  model(nt*(nq)))
  call planemod_init(o1,d1,q0,dq,nt,nx,nq)

  if (dotme) then
     j=size(model)
     i=size(data)
     write(0,*) 'Dot test'
     write(0,*) ' -------'
     !     write(0,*)' row_lop'
     !     call dot_test(row_lop,j,i,dot1,dot2)
     !     write(0,*)' vel_lop'
     !     call dot_test(velsimp6_lop,j/2,i,dot1,dot2)
     write(0,*) ' plane_lop'
     call dot_test(planemod_op,j,i,dot1,dot2)
     call seperr('Work???')
  end if
  blck=real(niter)/real(jniter)
  call solver_init(blck)
  do j=1,n3
     data=0.
     model=0.
     invert:if (adj) then 
        write(0,*) 'inverting data to mixed model domain'
        write(0,*) '>>>> Combo -- GATHER #',j,'/',n3
        call sep_read (data)
        write(0,'(a,f4.1)') ' Iteration blocks: ', blck
        !        call solver(model,data,row_lop,cgstep,add,niter,jniter)
        call solver(model,data,planemod_op,cgstep,add,niter,jniter)

        call to_history("iterations",niter+previts)
        if (exist_file("recons")) then
           st = planemod_op(.false.,.false.,model,data)
           call sep_write(data,"recons")
        end if
        call sep_write(model)
     else
        write(0,*)' reading model'
        call sep_read(model)
        st = planemod_op(.false.,.false.,model,data)
        call sep_write(data)
     end if invert
  end do
  call cgstep_close ( )
  call sep_close()
  write(0,*) 'Finished ========='
  call exit (0)

end program DirectRadon
