# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!
!
!
!
!
!
module dip_constraint_mod
  use adj_mod
  use param_check_mod
  use tomo_space_conv_mod
  use ray_database_mod
  use traceback_mod
  use basic_util_mod
  implicit none
  logical, private, save :: verb
  type(sep3d), private, save :: dip_pts_s
  real, private, save :: oang,dang
  real, pointer, private, save :: mask(:,:,:),d0(:),slow0(:,:,:)
  integer, private,save :: nt1,nt2,nt3,izero,n_dat
  real,private, save :: dt(3)
  contains
  logical function init_dip_constraint(ndat)
    integer :: npts,ndat
    real, pointer :: v(:,:,:),dip(:,:,:)
    integer :: n(3),np(3),na(3)
    real :: o(3),d(3),op(3),dp(3)
    type(sep3d) :: dip_s,vel_s
    init_dip_constraint=.false.
    call set_n_o_d_3d(dip_pts_s,n,o,d)
    na=n
    if (n(3).ne.1) then
      write(0,*) "dip constraint only written for 2-D for now"
      return
    end if
!READ IN THE FIELDS
    allocate(mask(n(1),n(2),n(3)),v(n(1),n(2),n(3)))
    call init_sep3d("dip_xz",dip_s,"INPUT")
    call set_n_o_d_3d(dip_s,np,op,dp)
    allocate(dip(np(1),np(2),np(3)))
    call sreed("dip_con_mask",mask,size(mask)*4)
    call sreed("dip_con_vals",v,size(mask)*4)
    call sseek("dip_xz",0,0)
    call sreed("dip_xz",dip,size(dip)*4)
    call sep3d_clean(dip_s)
    write(0,*) "size",size(v,1),size(v,2)
    write(0,*) "size",size(mask,1),size(mask,2)
    write(0,*) "size",size(dip,1),size(dip,2)
    call srite("da.H",v,size(v)*4)
    call srite("db.H",mask,size(mask)*4)
    call srite("dc.H",dip,size(dip)*4)
!CONVERT TO TAU
    if (.not. convert_field_tomo_space(mask)) then
      write(0,*) "    ERROR: converting mask from z to tau"
      return
    end if
    if (.not. convert_field_tomo_space(v)) then
      write(0,*) "    ERROR: converting dip constraint values from z to&
        & tau"
      return
    end if
    if (.not. field_z_to_tau(dip)) then
      write(0,*) "    ERROR: converting dip values from z to tau"
      return
    end if
    where(mask<1.)
      mask=0.
end where
    where(mask>1.)
      mask=1.
end where
    allocate(d0(size(v)))
    d0=reshape(mask*(v-dip),(/size(d0)/))
    write(0,*) "msa",size(mask,1),size(mask,2)
    call srite("d1.H",v,size(v)*4)
    call srite("d1.H",mask,size(mask)*4)
    call srite("d1.H",d0,size(d0)*4)
    call add_mem(size(d0)+2*size(mask))
    if (.not. vel_tomo_structure(vel_s)) then
      write(0,*) "trouble getting velo tomo struct"
      return
    end if
    call set_n_o_d_3d(vel_s,n,o,dt)
    nt1=n(1)
    nt2=n(2)
    nt3=n(3)
    allocate(slow0(nt1,nt2,nt3))
    if (.not. vel_tomo_3d(slow0)) then
      write(0,*) 'trouble grabbing velocity'
      return
    end if
    slow0=1./slow0
    izero=floor(-o(1)/dt(1)+1.5)
    n_dat=size(mask)
    ndat=size(mask)
    deallocate(v)
    allocate(v(na(1),na(2),na(3)))
    call sreed("abcde.H",v,size(v)*4)
    write(0,*) "CHECK1THIS",na
    if (.not. convert_vel_tomo_space(v)) then
      write(0,*) "    ERROR: converting mask from z to tau"
      return
    end if
    write(0,*) "CHECK2THIS",size(v,1),size(v,2)
    call srite("mtest3.H",v,size(dip)*4)
    v=1./v -slow0
    dip=0
    call dip_constraint_op3d(.false.,v,dip)
    call srite("mtest3.H",v,size(dip)*4)
    call srite("mtest3.H",dip,size(dip)*4)
    v=0
    call dip_constraint_op3d(.true.,v,dip)
    call srite("mtest3.H",v,size(dip)*4)
    write(0,*) "CHECK2THIS",size(v,1),size(v,2)
    deallocate(v,dip)
    init_dip_constraint=.true.
  end function
  integer function dip_constraint_op(adj,add,model,data)
    logical, intent (in) :: adj,add
    real :: model(:),data(:)
!else call srite("ttt.H",model,size(model)*4)
    call adjnull(adj,add,model,data)
    if (adj) then
      call srite("ttt.H",model,size(model)*4)
    end if
    call dip_constraint_op3d(adj,model,data)
    dip_constraint_op=0
    if (adj) then
      call srite("ttt.H",model,size(model)*4)
    end if
!else call srite("ttt.H",data,size(model)*4)
  end function
  subroutine dip_constraint_op3d(adj,slow,data)
    logical, intent(in) :: adj
    real :: slow(nt1,nt2,nt3),data(nt1,nt2,nt3)
    real(kind=wp) :: dz(nt1,nt2,nt3)
    integer :: ia,i1,i2,i3
    dz=0
    if (.not. adj) then
      do i1=izero+1,nt1
        dz(i1,:,:)=dz(i1-1,:,:)+dt(1)*.5*slow(i1-1,:,:)/slow0(i1-1,:,:&
          &)**2
      end do
      data(:,2:nt2-1,:)= data(:,2:nt2-1,:)+((dz(:,3:nt2,:)-dz(:,1:nt2&
        &-2,:))/(2*dt(2)))*mask(:,2:nt2-1,:)
! data(:,1,:)+=((dz(:,2,:)-dz(:,1,:))/(dt(2)))*mask(:,1,:)
! data(:,nt2,:)+=((dz(:,nt2,:)-dz(:,nt2-1,:))/(dt(2)))*mask(:,nt2,:)
    else
      dz(:,3:nt2,:) = dz(:,3:nt2,:) + data(:,2:nt2-1,:)/(2*dt(2))&
        &*mask(:,2:nt2-1,:)
      dz(:,1:nt2-2,:)= dz(:,1:nt2-2,:)- data(:,2:nt2-1,:)/(2*dt(2))&
        &*mask(:,2:nt2-1,:)
! dz(:,2,:)+= data(:,1,:)/(dt(2))*mask(:,1,:)
! dz(:,1,:)-= data(:,1,:)/(dt(2))*mask(:,1,:)
! dz(:,nt2,:)+= data(:,nt2,:)/(dt(2))*mask(:,nt2,:)
! dz(:,nt2-1,:)-= data(:,nt2,:)/(dt(2))*mask(:,nt2,:)
      do i3=1,nt3
        do i2=1,nt2
          do i1=izero+1,nt1
            do ia=i1-1,izero,-1
              slow(ia,i2,i3)= slow(ia,i2,i3)+.5*dt(1)*dz(i1,i2&
                &,i3)/slow0(ia,i2,i3)**2
            end do
          end do
        end do
      end do
    end if
    call srite("uuu.H",real(dz),size(dz)*4)
  end subroutine
  logical function get_extra_d_dip_con(d)
    real, pointer :: d(:)
    get_extra_d_dip_con=.false.
    allocate(d(size(d0)))
    d=d0
    get_extra_d_dip_con=.true.
  end function
  logical function check_dip_constraint_params()
    integer, external :: auxin
    logical :: check(10)
    type(sep3d) :: junk_s,vel_s
    integer :: ierr
    check_dip_constraint_params=.false.
    call from_param("verb",verb,.false.)
    check=.true.
    if (.not. vel_structure(vel_s)) then
      write(0,*) "trouble grabbing velocity tomo structure"
      return
    end if
    if (model_space(1:2).ne."ds") then
      write(0,*) "dip constraint only works when model space is ds"
      return
    end if
    if (-1.eq.auxin("dip_con_mask")) then
      write(0,*) "When doing dip constraing dip_con_mask must be&
        & specified"
      return
    end if
    call auxclose("dip_con_mask")
    call init_sep3d("dip_con_mask",dip_pts_s,"INPUT")
    if (-1.eq.auxin("dip_con_vals")) then
      write(0,*) "When doing dip constraing dip_con_vals must be&
        & specified"
      return
    end if
    call auxclose("dip_con_vals")
    call init_sep3d("dip_con_vals",junk_s,"INPUT")
    if (.not. sep3d_conform(junk_s,vel_s) .or. .not. sep3d_conform&
      &(junk_s,dip_pts_s)) then
      write(0,*) "dip_con_vals and dip_con_mask and  velocity must be&
        & the same shape"
      return
    end if
    call sep3d_clean(junk_s)
    check_dip_constraint_params=.true.
  end function
end module
