# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!x_t,z_t,map_t<
! Prec_tomo
!
! Usage
! Prec_tomo
!
! INPUT PARAMETERS
! verb - integer [0] verbose
! dot_test - integer [0] whether or not to run the dot test
! niter - integer [1] number of itteratopms to run
! Vel0 - sepfile Starting velocity
! R_traj - sepfile Ray trajectory
! R_rpar - sepfile Ray parameters
! R_mod - sepfile Ray Model
! R_time_ref-sepfile Times to the reflection
! R_time_tot-sepfile Total time of the ray
! R_time_0 -sepfile Modeled matched times
! R_time -sepfile Correct travel times
! lap_lags -integer [0] number of lags to use for the laplacian
! steer_lags-integer [3] number of steering lags
! tri_width -integer [2.5] width of triangle to use
! lap_factor-real [0.] weight for the laplacian portion
! steer_factor-real [.999] weight for the steering factor portion
! angles -sepfile [otional if lap=1] angles for steering filters
! eps -real [1.] epsilon parameter for operator
! weight -sepfile weight for bad times
! frame -integer [50] number of elements in frame (top and bottom to
! avoid helix wrap-around)
! taumap -sepfile tau map file
! dtau -real sampling in tau
! ntau -integer number of samples in tau
!
! OUTPUT PARAMETERS
! Vel_tomo - sepfile Model, etc
!
! DESCRIPTION
!
!
!
!>
program Depth_tomo
  use model_io_mod
  use ray_io_mod
  use inv_tomo_mod
  use tomo_oper_mod
  use filters_mod
  implicit none
  type (tomo_param_type) param
  type (model_param_type) model_param
  type (ray_obj_type) ray
  type(cube) :: map,x_t,x_z
  real, allocatable, dimension (:,:) :: time0_ray,d_time_ray,time_ray
  real, allocatable, dimension (:,:,:) :: d_slow,slow,tempa
  real, allocatable, dimension (:) :: tempb,tempc
  complex, pointer, dimension (:,:,:) :: hold_vel
  real, pointer, dimension (:,:,:) :: slow0,sigma0,d_sigma,sigma
  real, pointer :: weight_in(:)
  character(len=128) :: model_tag
  integer ierr,nelem_ray,nelem_model,verb,frame
  integer :: n_h,n_mid,auxout,i1,i2,tempi,temp2i,temp3i,type,test_oper
  real :: mid_0,h_0,d_h,d_mid,rone,decay,background,epsilon,scale_lap&
    &,limit
  real :: dtau
  real, pointer :: taumap(:,:)
  real :: o(3),d(3)
!FILTER PARAMETERS
  type(filter) :: filt
  integer :: lap_lags,steer_lags,ntau
  real :: tri_width,lap_factor,steer_factor,eps
  real,allocatable :: angles(:),tangle(:,:)
!READ IN PARAMETERS
  integer getch,putch,auxpar,auxputch
  call initpar()
  call doc('/homes/sep/claudio/research.ln/Prst_ExpRefl/Tese/src/areal&
    &/wemva/sep97_bob/Prec_tomo.rs90')
    if (0>= getch('new','d',param%new )) then
      param%new = 0
    end if
    if (0.ne.putch('Read  from par:  new ','d',param%new)) then
      call erexit('Trouble writing param%new to history')
    end if
    write(0,*) 'new: from par is obsolete, please replace with from&
      & param'
    if (0>= getch('verb','d',verb )) then
      verb = 0
    end if
    if (0.ne.putch('Read  from par:  verb ','d',verb)) then
      call erexit('Trouble writing verb to history')
    end if
    write(0,*) 'verb: from par is obsolete, please replace with from&
      & param'
    if (0>= getch('frame','d',frame )) then
      frame = 50
    end if
    if (0.ne.putch('Read  from par:  frame ','d',frame)) then
      call erexit('Trouble writing frame to history')
    end if
    write(0,*) 'frame: from par is obsolete, please replace with from&
      & param'
    if (0>= getch('verb_solver','d',param%verb_solver )) then
      param%verb_solver = verb
    end if
    if (0.ne.putch('Read  from par:  verb_solver ','d',param%verb_solver&
      &)) then
      call erexit('Trouble writing param%verb_solver to history')
    end if
    write(0,*) 'verb_solver: from par is obsolete, please replace with&
      & from param'
    if (0>= getch('type','d',type )) then
      type = 1
    end if
    if (0.ne.putch('Read  from par:  type ','d',type)) then
      call erexit('Trouble writing type to history')
    end if
    write(0,*) 'type: from par is obsolete, please replace with from&
      & param'
  param%tau_domain=0
  param%approx=0
  param%test_oper=0
  param%prec_tau=1
    if (0>= getch('dot_test','d',param%dot_test )) then
      param%dot_test = 0
    end if
    if (0.ne.putch('Read  from par:  dot_test ','d',param%dot_test))&
      & then
      call erexit('Trouble writing param%dot_test to history')
    end if
    write(0,*) 'dot_test: from par is obsolete, please replace with from&
      & param'
    if (0>= getch('test_oper','d',param%test_oper )) then
      param%test_oper = 0
    end if
    if (0.ne.putch('Read  from par:  test_oper ','d',param%test_oper))&
      & then
      call erexit('Trouble writing param%test_oper to history')
    end if
    write(0,*) 'test_oper: from par is obsolete, please replace with&
      & from param'
    if (0>= getch('niter','d',param%niter )) then
      param%niter = 1
    end if
    if (0.ne.putch('Read  from par:  niter ','d',param%niter)) then
      call erexit('Trouble writing param%niter to history')
    end if
    write(0,*) 'niter: from par is obsolete, please replace with from&
      & param'
    if (0>= getch('ndim_mod','d',model_param%ndim_mod )) then
      model_param%ndim_mod = 1
    end if
    if (0.ne.putch('Read  from par:  ndim_mod ','d',model_param%ndim_mod&
      &)) then
      call erexit('Trouble writing model_param%ndim_mod to history')
    end if
    write(0,*) 'ndim_mod: from par is obsolete, please replace with from&
      & param'
    if (0>= getch('ndim_ray','d',ray%ndim_ray )) then
      ray%ndim_ray = 2
    end if
    if (0.ne.putch('Read  from par:  ndim_ray ','d',ray%ndim_ray))&
      & then
      call erexit('Trouble writing ray%ndim_ray to history')
    end if
    write(0,*) 'ndim_ray: from par is obsolete, please replace with from&
      & param'
  if (auxin("Vel0") .eq. -1) then
    call erexit("Must supply Vel0")
  end if
    if (0>= auxpar('n1','d',model_param%nv_mod ,'Vel0')) then
      call erexit('Trouble reading  n1 from aux Vel0 ')
    end if
    if (0.ne.putch('Read  from aux: Vel0_n1 ','d',model_param%nv_mod))&
      & then
      call erexit('Trouble writing model_param%nv_mod to history')
    end if
    if (0>= auxpar('n2','d',model_param%nx_mod ,'Vel0')) then
      call erexit('Trouble reading  n2 from aux Vel0 ')
    end if
    if (0.ne.putch('Read  from aux: Vel0_n2 ','d',model_param%nx_mod))&
      & then
      call erexit('Trouble writing model_param%nx_mod to history')
    end if
    if (0>= auxpar('n3','d',model_param%ny_mod ,'Vel0')) then
      model_param%ny_mod = 1
    end if
    if (0.ne.putch('Read  from aux: Vel0_n3 ','d',model_param%ny_mod))&
      & then
      call erexit('Trouble writing model_param%ny_mod to history')
    end if
    if (0>= auxpar('o1','f',model_param%v0_mod ,'Vel0')) then
      call erexit('Trouble reading  o1 from aux Vel0 ')
    end if
    if (0.ne.putch('Read  from aux: Vel0_o1 ','f',model_param%v0_mod))&
      & then
      call erexit('Trouble writing model_param%v0_mod to history')
    end if
    if (0>= auxpar('o2','f',model_param%x0_mod ,'Vel0')) then
      call erexit('Trouble reading  o2 from aux Vel0 ')
    end if
    if (0.ne.putch('Read  from aux: Vel0_o2 ','f',model_param%x0_mod))&
      & then
      call erexit('Trouble writing model_param%x0_mod to history')
    end if
    if (0>= auxpar('o3','f',model_param%y0_mod ,'Vel0')) then
      model_param%y0_mod = 0.
    end if
    if (0.ne.putch('Read  from aux: Vel0_o3 ','f',model_param%y0_mod))&
      & then
      call erexit('Trouble writing model_param%y0_mod to history')
    end if
    if (0>= auxpar('d1','f',model_param%dv_mod ,'Vel0')) then
      call erexit('Trouble reading  d1 from aux Vel0 ')
    end if
    if (0.ne.putch('Read  from aux: Vel0_d1 ','f',model_param%dv_mod))&
      & then
      call erexit('Trouble writing model_param%dv_mod to history')
    end if
    if (0>= auxpar('d2','f',model_param%dx_mod ,'Vel0')) then
      call erexit('Trouble reading  d2 from aux Vel0 ')
    end if
    if (0.ne.putch('Read  from aux: Vel0_d2 ','f',model_param%dx_mod))&
      & then
      call erexit('Trouble writing model_param%dx_mod to history')
    end if
    if (0>= auxpar('d3','f',model_param%dy_mod ,'Vel0')) then
      model_param%dy_mod = 1.
    end if
    if (0.ne.putch('Read  from aux: Vel0_d3 ','f',model_param%dy_mod))&
      & then
      call erexit('Trouble writing model_param%dy_mod to history')
    end if
  if (auxin("R_traj") .eq. -1) then
    call erexit("Must supply R_traj")
  end if
    if (0>= auxpar('n1','d',ray%nt_ray ,'R_traj')) then
      call erexit('Trouble reading  n1 from aux R_traj ')
    end if
    if (0.ne.putch('Read  from aux: R_traj_n1 ','d',ray%nt_ray)) then
      call erexit('Trouble writing ray%nt_ray to history')
    end if
    if (0>= auxpar('n2','d',n_h ,'R_traj')) then
      call erexit('Trouble reading  n2 from aux R_traj ')
    end if
    if (0.ne.putch('Read  from aux: R_traj_n2 ','d',n_h)) then
      call erexit('Trouble writing n_h to history')
    end if
    if (0>= auxpar('n3','d',n_mid ,'R_traj')) then
      call erexit('Trouble reading  n3 from aux R_traj ')
    end if
    if (0.ne.putch('Read  from aux: R_traj_n3 ','d',n_mid)) then
      call erexit('Trouble writing n_mid to history')
    end if
  ray%n_ray=n_h*n_mid
    if (0>= auxpar('n4','d',ray%nref_ray ,'R_traj')) then
      call erexit('Trouble reading  n4 from aux R_traj ')
    end if
    if (0.ne.putch('Read  from aux: R_traj_n4 ','d',ray%nref_ray))&
      & then
      call erexit('Trouble writing ray%nref_ray to history')
    end if
    if (0>= auxpar('o1','f',ray%t0_ray ,'R_traj')) then
      ray%t0_ray = 0.
    end if
    if (0.ne.putch('Read  from aux: R_traj_o1 ','f',ray%t0_ray)) then
      call erexit('Trouble writing ray%t0_ray to history')
    end if
    if (0>= auxpar('o2','f',h_0 ,'R_traj')) then
      call erexit('Trouble reading  o2 from aux R_traj ')
    end if
    if (0.ne.putch('Read  from aux: R_traj_o2 ','f',h_0)) then
      call erexit('Trouble writing h_0 to history')
    end if
    if (0>= auxpar('o3','f',mid_0 ,'R_traj')) then
      call erexit('Trouble reading  o3 from aux R_traj ')
    end if
    if (0.ne.putch('Read  from aux: R_traj_o3 ','f',mid_0)) then
      call erexit('Trouble writing mid_0 to history')
    end if
    if (0>= auxpar('d1','f',ray%dt_ray ,'R_traj')) then
      ray%dt_ray = 0.
    end if
    if (0.ne.putch('Read  from aux: R_traj_d1 ','f',ray%dt_ray)) then
      call erexit('Trouble writing ray%dt_ray to history')
    end if
    if (0>= auxpar('d2','f',d_h ,'R_traj')) then
      call erexit('Trouble reading  d2 from aux R_traj ')
    end if
    if (0.ne.putch('Read  from aux: R_traj_d2 ','f',d_h)) then
      call erexit('Trouble writing d_h to history')
    end if
    if (0>= auxpar('d3','f',d_mid ,'R_traj')) then
      call erexit('Trouble reading  d3 from aux R_traj ')
    end if
    if (0.ne.putch('Read  from aux: R_traj_d3 ','f',d_mid)) then
      call erexit('Trouble writing d_mid to history')
    end if
  param%tau_domain=1
  if (auxout("Vel_tomo") .eq. -1) then
    call erexit("Must supply R_traj")
  end if
  if (((param%niter >= 0).and.(param%test_oper .eq. 0)).or.(param%test_oper&
    & .eq. -1)) then
      if (0.ne. auxputch('esize','d',4 ,'Vel_tomo')) then
        call erexit('Trouble writing  esize to aux Vel_tomo ')
      end if
      if (0.ne. auxputch('n1','d',model_param%nv_mod ,'Vel_tomo'))&
        & then
        call erexit('Trouble writing  n1 to aux Vel_tomo ')
      end if
      if (0.ne. auxputch('n2','d',model_param%nx_mod ,'Vel_tomo'))&
        & then
        call erexit('Trouble writing  n2 to aux Vel_tomo ')
      end if
      if (0.ne. auxputch('n3','d',model_param%ny_mod ,'Vel_tomo'))&
        & then
        call erexit('Trouble writing  n3 to aux Vel_tomo ')
      end if
      if (0.ne. auxputch('o1','f',model_param%v0_mod ,'Vel_tomo'))&
        & then
        call erexit('Trouble writing  o1 to aux Vel_tomo ')
      end if
      if (0.ne. auxputch('o2','f',model_param%x0_mod ,'Vel_tomo'))&
        & then
        call erexit('Trouble writing  o2 to aux Vel_tomo ')
      end if
      if (0.ne. auxputch('o3','f',model_param%y0_mod ,'Vel_tomo'))&
        & then
        call erexit('Trouble writing  o3 to aux Vel_tomo ')
      end if
      if (0.ne. auxputch('d1','f',model_param%dv_mod ,'Vel_tomo'))&
        & then
        call erexit('Trouble writing  d1 to aux Vel_tomo ')
      end if
      if (0.ne. auxputch('d2','f',model_param%dx_mod ,'Vel_tomo'))&
        & then
        call erexit('Trouble writing  d2 to aux Vel_tomo ')
      end if
      if (0.ne. auxputch('d3','f',model_param%dy_mod ,'Vel_tomo'))&
        & then
        call erexit('Trouble writing  d3 to aux Vel_tomo ')
      end if
  end if
  allocate (sigma0(model_param%nv_mod,model_param%nx_mod,model_param%ny_mod&
    &),stat=ierr)
  if (ierr .ne. 0) then
    call seperr("Problems with allocation of sigma0")
  end if
  allocate (d_sigma(model_param%nv_mod,model_param%nx_mod,model_param%ny_mod&
    &),stat=ierr)
  if (ierr .ne. 0) then
    call seperr("Problems with allocation of d_sigma")
  end if
  allocate (sigma(model_param%nv_mod,model_param%nx_mod,model_param%ny_mod&
    &),stat=ierr)
  if (ierr .ne. 0) then
    call seperr("Problems with allocation of sigma")
  end if
  if (verb.eq.1) then
    write(0,*) "Done with initial parameter reading"
  end if
!ALLOCATE MODEL
  allocate (hold_vel(model_param%nv_mod,model_param%nx_mod,model_param%ny_mod&
    &),stat=ierr)
  allocate (slow0(model_param%nv_mod+frame*2,model_param%nx_mod&
    &,model_param%ny_mod),stat=ierr)
  if (ierr .ne. 0) then
    call seperr("Problems with allocation of slow0")
  end if
  allocate (tempa(2*frame+model_param%nv_mod,model_param%nx_mod&
    &,model_param%ny_mod),stat=ierr)
  if (ierr .ne. 0) then
    call seperr("Problems with allocation of tempa")
  end if
  allocate (d_slow(model_param%nv_mod+frame*2,model_param%nx_mod&
    &,model_param%ny_mod),stat=ierr)
  if (ierr .ne. 0) then
    call seperr("Problems with allocation of _slow")
  end if
  allocate (slow(model_param%nv_mod+frame*2,model_param%nx_mod&
    &,model_param%ny_mod),stat=ierr)
  if (ierr .ne. 0) then
    call seperr("Problems with allocation of slow")
  end if
  if (verb.eq.1) then
    write(0,*) "Done with initial array allocation"
  end if
!ALLOCATE RAYS
  call alloc_ray (param, ray)
  if (verb.eq.1) then
    write(0,*) "Done with array array allocation"
  end if
  allocate (time0_ray(ray%n_ray,ray%nref_ray),stat=ierr)
  if (ierr .ne. 0) then
    call seperr("Problems with allocation of time0_ray")
  end if
  allocate (d_time_ray(ray%n_ray,ray%nref_ray),stat=ierr)
  if (ierr .ne. 0) then
    call seperr("Problems with allocation of d_time_ray")
  end if
  allocate (time_ray(ray%n_ray,ray%nref_ray),stat=ierr)
  if (ierr .ne. 0) then
    call seperr("Problems with allocation of time_ray")
  end if
  if (verb.eq.1) then
    write(0,*) "Done with array allocation 2"
  end if
!########
  if (size(hold_vel)*8 .ne. sreed("Vel0",hold_vel,size(hold_vel)*8))&
    & then
    call erexit("trouble reading in vel0")
  end if
  slow0(frame+1:frame+model_param%nv_mod,:,:)=real(hold_vel)
  sigma0(frame+1:frame+model_param%nv_mod,:,:)=aimag(hold_vel)
  do i1=1,frame
    slow0(i1,:,:)=slow0(frame+1,:,:)
    slow0(frame+model_param%nv_mod+i1,:,:)=slow0(frame+model_param%nv_mod&
      &,:,:)
  end do
  if (verb.eq.1) then
    write(0,*)'Read Vel0'
  end if
  model_param%v0_mod=model_param%v0_mod-model_param%dv_mod*frame
  model_param%nv_mod=model_param%nv_mod+frame*2
  x_z%o=(/model_param%v0_mod,model_param%x0_mod,model_param%y0_mod/)
  x_z%d=(/model_param%dv_mod,model_param%dx_mod,model_param%dy_mod/)
  x_z%n=(/model_param%nv_mod,model_param%nx_mod,model_param%ny_mod/)
  where(slow0.ne.0.)
    slow0=1./slow0
end where
  if (verb.eq.1) then
    write(0,*) "Done reading velocity model"
  end if
  call read_ray ( param , ray )
  if (verb.eq.1) then
    write(0,*)'Done reading rays'
  end if
  call read_time_0_ray ( ray, time0_ray )
  if (verb.eq.1) then
    write(0,*)'Don reading 0 ray'
  end if
!##################
!CONSTRUCT FILTERS#
!##################
    if (0>= getch('lap_lags','d',lap_lags )) then
      lap_lags = 0
    end if
    if (0.ne.putch('Read  from par:  lap_lags ','d',lap_lags)) then
      call erexit('Trouble writing lap_lags to history')
    end if
    write(0,*) 'lap_lags: from par is obsolete, please replace with from&
      & param'
    if (0>= getch('steer_lags','d',steer_lags )) then
      steer_lags = 0
    end if
    if (0.ne.putch('Read  from par:  steer_lags ','d',steer_lags))&
      & then
      call erexit('Trouble writing steer_lags to history')
    end if
    write(0,*) 'steer_lags: from par is obsolete, please replace with&
      & from param'
    if (0>= getch('tri_width','f',tri_width )) then
      tri_width = 2.5
    end if
    if (0.ne.putch('Read  from par:  tri_width ','f',tri_width)) then
      call erexit('Trouble writing tri_width to history')
    end if
    write(0,*) 'tri_width: from par is obsolete, please replace with&
      & from param'
    if (0>= getch('lap_factor','f',lap_factor )) then
      lap_factor = 0.
    end if
    if (0.ne.putch('Read  from par:  lap_factor ','f',lap_factor))&
      & then
      call erexit('Trouble writing lap_factor to history')
    end if
    write(0,*) 'lap_factor: from par is obsolete, please replace with&
      & from param'
    if (0>= getch('steer_factor','f',steer_factor )) then
      steer_factor = .999
    end if
    if (0.ne.putch('Read  from par:  steer_factor ','f',steer_factor))&
      & then
      call erexit('Trouble writing steer_factor to history')
    end if
    write(0,*) 'steer_factor: from par is obsolete, please replace with&
      & from param'
    if (0>= getch('eps','f',eps )) then
      eps = 1.
    end if
    if (0.ne.putch('Read  from par:  eps ','f',eps)) then
      call erexit('Trouble writing eps to history')
    end if
    write(0,*) 'eps: from par is obsolete, please replace with from&
      & param'
  write(0,*) "width lap_factor lap_lags",tri_width,lap_factor,lap_lags
  write(0,*) "steerlags steer_factor eps",steer_lags,steer_factor,eps
  call init_filters(lap_scale=lap_factor,steer_scale=steer_factor&
    &,lap_lags=lap_lags,steer_lags=steer_lags,triangle_width=tri_width&
    &,run_triangle=.true.)
  if (-1.eq.auxin('angles')) then
    if (steer_lags.eq.0) then
      filt=new_steering((/model_param%nv_mod,model_param%nx_mod/))
    else
      call erexit("if steer_lag .ne. 0 angles must be provided")
    end if
  else
      if (0>= auxpar('n1','d',tempi ,'angles')) then
        call erexit('Trouble reading  n1 from aux angles ')
      end if
      if (0.ne.putch('Read  from aux: angles_n1 ','d',tempi)) then
        call erexit('Trouble writing tempi to history')
      end if
      if (0>= auxpar('n2','d',temp2i ,'angles')) then
        call erexit('Trouble reading  n2 from aux angles ')
      end if
      if (0.ne.putch('Read  from aux: angles_n2 ','d',temp2i)) then
        call erexit('Trouble writing temp2i to history')
      end if
    if (tempi+frame*2.ne.model_param%nv_mod .or. temp2i.ne.model_param%nx_mod&
      &) then
      write(0,*) "Angles",tempi,temp2i
      write(0,*) "model",model_param%nv_mod,model_param%nx_mod
      call erexit("angles not the same size as model")
    end if
    allocate(angles((tempi+frame*2)*temp2i),tangle(tempi,temp2i))
    if (sreed('angles',tangle,tempi*temp2i*4).ne.tempi*temp2i*4) then
      call erexit("trouble reading in angles ")
    end if
    ierr=0
    do i2=1,temp2i
      do i1=1,tempi
        angles(i1+frame+(i2-1)*(2*frame+tempi))=tangle(i1,i2)
      end do
      do i1=1,frame
        angles(i1+(i2-1)*(2*frame+tempi))=angles(frame+1+(i2-1)*(2*frame&
          &+tempi))
      end do
      do i1=1,frame
        angles(i1+frame+tempi+(i2-1)*(2*frame+tempi))= angles(frame&
          &+tempi+(i2-1)*(2*frame+tempi))
      end do
    end do
    filt=new_steering((/model_param%nv_mod,model_param%nx_mod/),angles&
      &,1)
    deallocate(angles,tangle)
  end if
  if (verb.eq.1) then
    write(0,*) "Done constructing angles"
  end if
!###########
!WEIGHT IN##
!###########
  if (-1.eq.auxin("weight")) then
    call erexit("must supply weight file")
  end if
    if (0>= auxpar('n1','d',tempi ,'weight')) then
      call erexit('Trouble reading  n1 from aux weight ')
    end if
    if (0.ne.putch('Read  from aux: weight_n1 ','d',tempi)) then
      call erexit('Trouble writing tempi to history')
    end if
    if (0>= auxpar('n2','d',temp2i ,'weight')) then
      call erexit('Trouble reading  n2 from aux weight ')
    end if
    if (0.ne.putch('Read  from aux: weight_n2 ','d',temp2i)) then
      call erexit('Trouble writing temp2i to history')
    end if
    if (0>= auxpar('n3','d',temp3i ,'weight')) then
      call erexit('Trouble reading  n3 from aux weight ')
    end if
    if (0.ne.putch('Read  from aux: weight_n3 ','d',temp3i)) then
      call erexit('Trouble writing temp3i to history')
    end if
  allocate(weight_in(tempi*temp2i*temp3i))
  ierr= sreed("weight",weight_in,size(weight_in)*4)
!###################
!TAU PRECONDITIONER#
!###################
    if (0>= getch('ntau','d',ntau )) then
      call erexit('Trouble reading  ntau from par  ')
    end if
    if (0.ne.putch('Read  from par:  ntau ','d',ntau)) then
      call erexit('Trouble writing ntau to history')
    end if
    write(0,*) 'ntau: from par is obsolete, please replace with from&
      & param'
    if (0>= getch('dtau','f',dtau )) then
      call erexit('Trouble reading  dtau from par  ')
    end if
    if (0.ne.putch('Read  from par:  dtau ','f',dtau)) then
      call erexit('Trouble writing dtau to history')
    end if
    write(0,*) 'dtau: from par is obsolete, please replace with from&
      & param'
  if (-1.eq.auxin("taumap")) then
    call erexit("must supply tau map")
  end if
    if (0>= auxpar('n1','d',tempi ,'taumap')) then
      call erexit('Trouble reading  n1 from aux taumap ')
    end if
    if (0.ne.putch('Read  from aux: taumap_n1 ','d',tempi)) then
      call erexit('Trouble writing tempi to history')
    end if
    if (0>= auxpar('n2','d',temp2i ,'taumap')) then
      call erexit('Trouble reading  n2 from aux taumap ')
    end if
    if (0.ne.putch('Read  from aux: taumap_n2 ','d',temp2i)) then
      call erexit('Trouble writing temp2i to history')
    end if
  o1=0.
  d1=1.
    if (0>= auxpar('o1','f',o(1) ,'taumap')) then
      call erexit('Trouble reading  o1 from aux taumap ')
    end if
    if (0.ne.putch('Read  from aux: taumap_o1 ','f',o(1))) then
      call erexit('Trouble writing o(1) to history')
    end if
    if (0>= auxpar('o2','f',o(2) ,'taumap')) then
      call erexit('Trouble reading  o2 from aux taumap ')
    end if
    if (0.ne.putch('Read  from aux: taumap_o2 ','f',o(2))) then
      call erexit('Trouble writing o(2) to history')
    end if
    if (0>= auxpar('d1','f',d(1) ,'taumap')) then
      call erexit('Trouble reading  d1 from aux taumap ')
    end if
    if (0.ne.putch('Read  from aux: taumap_d1 ','f',d(1))) then
      call erexit('Trouble writing d(1) to history')
    end if
    if (0>= auxpar('d2','f',d(2) ,'taumap')) then
      call erexit('Trouble reading  d2 from aux taumap ')
    end if
    if (0.ne.putch('Read  from aux: taumap_d2 ','f',d(2))) then
      call erexit('Trouble writing d(2) to history')
    end if
  allocate(taumap(tempi,temp2i))
  if (tempi*temp2i*4.ne.sreed("taumap",taumap,size(taumap)*4)) then
    call erexit("trouble reading taumap")
  end if
  map%n=(/tempi,temp2i,1/)
  map%o=o
  map%d=d
  x_t%n=(/ntau+2*frame,x_z%n(2:3)/)
  x_t%o=(/-dtau*frame,x_z%o(2:3)/)
  x_t%d=(/dtau,x_z%d(2:3)/)
  call init_map(x_t,z_t,map_t,taumap)
!############
!ARRAY SIZES#
!############
  nelem_ray=ray%n_ray*ray%nref_ray
  nelem_model=model_param%ndim_mod*model_param%nv_mod*model_param%nx_mod&
    &*model_param%ny_mod
  if (verb.eq.1) then
    write(0,*)'nelem_ray=',nelem_ray
    write(0,*)'nelem_model=',nelem_model
  end if
!###############
!INIT OPERATORS#
!###############
  if (verb.eq.1) then
    write(0,*) "initialize operators",size(sigma0)
  end if
  call tomo_oper_init (param,model_param,ray,sigma0,slow0)
  if (verb.eq.1) then
    write(0,*) "before chain"
  end if
  call helix_chain_init(nelem_model)
  if (verb.eq.1) then
    write(0,*) "before polydiv",size(filt%coeffs)
  end if
  call polydiv_init(filt)
  if (verb.eq.1) then
    write(0,*) "after polydif"
  end if
!#################
!RUN THE OPERATOR#
!#################
  if (param%dot_test .eq. 1) then
    if (verb.eq.1) then
      write(0,*) "before dottest"
    end if
    call inv_tomo (param,model_param,ray, nelem_ray,nelem_model, &
      & d_time_ray,d_slow,sigma0,slow0,tau_chain_oper,eps,weight_in)
  else
    if (verb.eq.1) then
      write(0,*)'Invert or Backproject'
    end if
    call read_time_ray ( ray, time_ray )
    d_time_ray=time_ray-time0_ray
    if (verb.eq.1) then
      write(0,*) "size of time_ray",size(time_ray)
      write(0,*)'minval(d_time_ray)=',minval(d_time_ray)
      write(0,*)'maxval(d_time_ray)=',maxval(d_time_ray)
      write(0,*)'minval(d_slow)=',minval(d_slow)
      write(0,*)'maxval(d_slow)=',maxval(d_slow)
    end if
    write(0,*) "calling inv tomo"
    call inv_tomo(param,model_param,ray,nelem_ray,nelem_model,d_time_ray&
      &,d_slow,sigma0,slow0,helix_chain2_oper,eps,weight_in)
    if (verb.eq.1) then
      write(0,*) "Done with inversion"
    end if
    allocate(tempb(nelem_model),tempc(nelem_model))
    tempi=0
    do i2=1,model_param%nx_mod
      do i1=1,model_param%nv_mod
        tempi= tempi+1
        tempb(tempi)=d_slow(i1,i2,1)
      end do
    end do
    write(0,*) "tempb",minval(tempb),maxval(tempb)
    if (verb.eq.1) then
      write(0,*) "before polydiv",tempi,size(d_slow)
    end if
    ierr=polydiv_lop(.false.,.false.,tempb,tempc)
    write(0,*) "tempc",minval(tempc),maxval(tempc)
    if (verb.eq.1) then
      write(0,*) "before conversion2",size(tempa),size(tempc)
    end if
    tempi=0
    do i2=1,model_param%nx_mod
      do i1=1,model_param%nv_mod
        tempi= tempi+1
        tempa(i1,i2,1)=tempc(tempi)
      end do
    end do
    write(0,*) "tempa",minval(tempa),maxval(tempa)
    write(0,*) "slow0",minval(slow0),maxval(slow0)
    if (verb.eq.1) then
      write(0,*) "before addition",tempi,size(slow),size(tempa),size&
        &(tempc)
    end if
    slow=slow0+tempa
    write(0,*) "slo0",minval(slow),maxval(slow)
    slow=1./slow
    deallocate(tempa,stat=ierr)
    if (ierr.ne.0) then
      call erexit("trouble deallocating helix arrays")
    end if
      if (0.ne. putch('n1','d',model_param%nv_mod )) then
        call erexit('Trouble writing  n1 to history  ')
      end if
      if (0.ne. putch('n2','d',model_param%nx_mod )) then
        call erexit('Trouble writing  n2 to history  ')
      end if
      if (0.ne. putch('n3','d',3 )) then
        call erexit('Trouble writing  n3 to history  ')
      end if
    model_param%nv_mod=model_param%nv_mod-frame*2
    write(0,*) model_param%v0_mod,"after2"
    write(0,*) size(slow,1),size(slow,2),minval(slow),maxval(slow)
    call write_model ('Vel_tomo', param , model_param, slow(frame&
      &+1:model_param%nv_mod+frame,:,:))
! call write_model ('Vel_tomo', param , model_param, slow)
  end if
  if (verb.eq.1) then
    write(0,*) "Done with tomo"
  end if
end program
