program CombRadTran

  use sep
  !solver mods
  use huber_reg_mod_dp
  use huber_bfgs_b_dp
  !weights and regularizors
  use identity_mod_dp
  use weighthrt_dp
  !model operators
  use row_mod2_dp
  use velsimp6_dp
!  use linearsimp_mod_dp
  use planemod2
  implicit none
!  INTEGER, PARAMETER      :: dbl = SELECTED_REAL_KIND(15,307)
  INTEGER, PARAMETER      :: dbl = 8


!from HRT.f90
  integer                          :: nt, nx, ns, iadj, n2, niter, stat,its,ntr,ntrr
  integer                          :: mmax, l1, n3, verb, j, bound, em ,bm
  real(kind=dbl)                   :: gamd, gamx, eps, nu, mu, accuracy
  real                             :: spgamd, spgamx, speps, spnu, spmu, spaccuracy
  real                             :: o1,d1, o2,d2, o3, d3, x0,dx, s0,ds, llow, lup, hlow, hup
  real, dimension (:), allocatable :: spdata, spmodel
  real, dimension (:), pointer     :: spmask
  real(kind=dbl), dimension (:), allocatable :: data, model, recons, l, u
  real(kind=dbl), dimension (:), pointer     :: mask
  logical                          :: adj
!from LRT.f90
  integer                          :: weight, inversion,nq! , alias
  real                             :: pi!,fhi,flo!,gamd,gamx
  real                             :: q0,dq
  

  call sep_init ()                 ! start processing parameters
  call from_history ("n3", n3)
!  call from_history ("o3", o3)
!  call from_history ("d3", d3)

  call from_history ("o1", o1)
  call from_history ("d1", d1)
  call from_history ("n1", nt)
  call from_history ("n2", ntrr)
  call from_history ("n2",ntr,"geom")
  if(ntrr.ne.ntr) call seperr("wrong number of headers")

  call from_history ("o2", o2)
  call from_history ("d2", d2)
  nx=ntr
  x0=0.
  dx=1.
!  call from_param  ("nx", nx, n2)
!  call from_param  ("x0", x0, o2) 
!  call from_param  ("dx", dx, d2) 

  !!HRT stuff
  call from_param  ("ns", ns, 2*nx)
  call from_param  ("s0", s0) 
  call from_param  ("ds", ds)  
  !!direct arrival stuff
  call from_param  ("nq", nq, ns)
  call from_param  ("q0", q0, s0)
  call from_param  ("dq", dq, ds)
  !----------------------------------------------------------------
  call from_param   ("adj",iadj,1) ; adj = (iadj == 1)

  if (adj) then 
     call to_history("ns", ns)
     call to_history("s0", s0) 
     call to_history("ds", ds)  
     call to_history("label2","slowness(sec/km)")
     call to_history("nq", nq)
     call to_history("q0", q0) 
     call to_history("dq", dq)  
     n2 = ns+nq; o2 = s0; d2 = ds
  else
     call to_history ("label2","distance(km)")
     n2 = nx; o2 = x0; d2 = dx
  end if


  !!HRT model space
  call to_history ("o2", o2)
  call to_history ("d2", d2)
  call to_history ("n2", n2)
  call from_param ("verb",verb,2)
  call from_param ("niter", niter, 1)
  call from_param ("mmax",mmax,2)
  call from_param ("gamd",spgamd,1.)
  gamd=spgamd
  call from_param ("gamx",spgamx,1.)
  gamx=spgamx
  call from_param ("nu",spnu,0.9)
  nu=spnu
  call from_param ("mu",spmu,0.0001)
  mu=spmu
  call from_param ("eps",speps,0.)
  eps=speps
  call from_param ("l1",l1,0)
  call from_param ("bound_type",bound,0)
  call from_param ("accuracy",spaccuracy,1e-16)
  accuracy=spaccuracy
  call from_param ("lower_h",hlow,0.)
  call from_param ("upper_h",hup,1.)
  !!LRT model space
  call from_param   ("inversion",inversion,3)
  call from_param   ("adj",iadj,1) ; adj = (iadj == 1)
  call from_param   ("niter",niter,10)
  call from_param   ("lower_l",llow,0.)
  call from_param   ("upper_l",lup,1.)
  !----------------------------------------------------------------
  if (exist_file ("recons")) then
     call to_history ("o1", o1,"recons")
     call to_history ("d1", d1,"recons")
     call to_history ("n1", nt,"recons")
     call to_history ("o2", x0,"recons")
     call to_history ("d2", dx,"recons")
     call to_history ("n2", nx,"recons")
     call to_history ("o3", o3,"recons")
     call to_history ("d3", d3,"recons")
     call to_history ("n3", n3,"recons")
     call to_history ("label2","distance(km)","recons")
     allocate(recons(nt*nx))
  end if
  !----------------------------------------------------------------
  if (verb>0) write(0,*) ' '
  if (verb>0) write(0,*) ' '
  if (verb>0) write(0,*) '>> Valhall Direct Arrival & HRT in the Time Domain'
  if (verb>0) write(0,*) '>> '
  if (.not.adj) write(0,*) '>> BE SURE THAT HRT MODEL SPACE IS FIRST '
  if (verb>0) write(0,*) '>> '
  if (verb>0) write(0,*) '>> '


  if (verb>0) write(0,*) '>>  '
  if (verb>0) write(0,*) '>>  '
  if (verb>0) write(0,*) '>>      parameters'
  if (verb>0) write(0,*) '>> GEN - gamd  =:',gamd
  if (verb>0) write(0,*) '>> GEN - gamx  =:',gamx
  if (verb>0) write(0,*) '>> GEN - niter =:',niter
  if (verb>0) write(0,*) '>> GEN - eps   =:',eps
  if (verb>0) write(0,*) '>> GEN - l1    =:',l1
  if (verb>0) write(0,*) '>> '
  if (verb>0) write(0,*) '>> '
  if (verb>0) write(0,*) '>> DAT - x0  =:',x0
  if (verb>0) write(0,*) '>> DAT - dx  =:',dx
  if (verb>0) write(0,*) '>> DAT - nx  =:',nx
  if (verb>0) write(0,*) '>>     - '
  if (verb>0) write(0,*) '>>     - '
  if (verb>0) write(0,*) '>> HRT - s0  =:',s0
  if (verb>0) write(0,*) '>> HRT - ns  =:',ns
  if (verb>0) write(0,*) '>> HRT - ds  =:',ds
  if (verb>0) write(0,*) '>> HRT - verb=:',verb
  if (verb>0) write(0,*) '>>  '
  if (verb>0) write(0,*) '>>  ' 
  if (verb>0) write(0,*) '>> VDA - q0  =:',q0
  if (verb>0) write(0,*) '>> VDA - nq  =:',nq
  if (verb>0) write(0,*) '>> VDA - dq  =:',dq
  if (verb>0) write(0,*) '>> VDA - verb=:',verb
  if (verb>0) write(0,*) '>> VDA - inv =:',inversion
  if (verb>0) write(0,*) '>>  '
  if (verb>0) write(0,*) '>>  '



  allocate (  data(nt*nx),  model(nt*(ns+nq)), l(nt*(ns+nq)), u(nt*(ns+nq)))
  allocate (spdata(nt*nx),spmodel(nt*(ns+nq)))
  bm=size(u)
  em=nt*ns
  l(1:em)=hlow  ; l(em+1:bm)=llow
  u(1:em)=hup   ; u(em+1:bm)=lup
  if (exist_file("mask")) then
     allocate(spmask(nt*nx))
     allocate(mask(nt*nx))
  end if

  call row_init(o1,d1,x0,dx,q0,dq,nt,nx,nq,verb,s0,ds,ns)
  call huber_init_out(gamd,gamx,l1,eps,nt*(ns+nq),verb)
  do j=1,n3
     data=0.
     model=0.
     if (adj) then 
        data=0.
        model=0.

        if (verb>0) write(0,*) '>>>> Combo -- GATHER #',j,'/',n3

        if (exist_file("mask")) then
           mask=0.
           call sep_read(spmask,"mask")
           mask=spmask
           deallocate(spmask)

           call weight_init(mask)
           call sep_read (spdata)
           data=spdata

           data=data*mask
           if (bound.eq.3) then
              where (data .lt. 0.) data=0.
           elseif(bound.eq.1) then
              where (data .gt. 0.) data=0.
           end if
           if(gamd.eq.1) gamd=abs(maxval(data))
           call huber_init_in(data,gamd)
           call lbfgs (row_lop,identity_lop,huberfct,hubergdt,model,niter,accuracy,bound,l,u,verb,its=its)
           call huber_close()
           call sep_write (spmodel)

        else

           call sep_read (spdata)
           data=spdata
           if (bound.eq.1) then
              where (data .lt. 0.0000) data=0.0000
           elseif(bound.eq.3) then
              where (data .gt. 0.0000) data=0.0000
           end if
           call huber_init_in(data,gamd)
           call lbfgs (row_lop,identity_lop,huberfct,hubergdt,model,niter,accuracy,bound,l,u,verb,its=its)
           call huber_close()
           spmodel=model
           call sep_write (spmodel)

        end if

        if (exist_file("recons")) then
           stat = row_lop(.false.,.false.,model,data)           
           spdata=data
           call sep_write(spdata,"recons")
        end if

     else
        call sep_read (spmodel)
        model=spmodel
        stat = row_lop (.false.,.false.,model,data)
        spdata=data
        call sep_write(spdata)
     end if

  end do

  deallocate (model,data,spmodel,spdata,l,u)
  if (exist_file("mask")) deallocate(mask)
  if (exist_file ("recons")) deallocate(recons)
  if(adj) then
     write(0,*) "*************************"
     write(0,*) "Output has 2 model spaces packed into one."
     write(0,*) ' o2 & d2 are correct only for the HRT of n2=',ns
     write(0,*) '                              The output n2=',ns+nq
     write(0,*) 'All the real values (ns,nq,dq,oq) are in the history file.'
  end if
  write(0,*) 'Finished ========='
  if(adj) call to_history("iterations",its)
  call sep_close()                ! end processing parameters
  call exit (0)

end program CombRadTran
