# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!
!Mon Jan 13 21:20:40 PST 2003
!
module calc_sigma_mod
  use adj_mod
  implicit none
  real, private, save, pointer :: sigma(:,:,:),sigmax(:,:,:),sigmay(:,:&
    &,:),slow0(:,:,:)
  real, private, save :: d(3)
  logical, private, save :: do_x
  integer, private, save :: n1h,n2h,n3h,n(3)
  contains
  subroutine calc_delta_sigmay(adj,add,model,data)
    logical :: adj,add
    real :: model(n1h,n2h,n3h),data(n1h,n2,n3h)
    if (.not. add) then
      if (adj) then
        model=0
      else
        data=0
      end if
    end if
    do_x=.false.
    sigma=>sigmay
    call calc_delta_sigma_3d(adj,add,model,data)
  end subroutine
  subroutine calc_delta_sigmax(adj,add,model,data)
    logical :: adj,add
    real :: model(n1h,n2h,n3h),data(n1h,n2,n3h)
    if (.not. add) then
      if (adj) then
        model=0
      else
        data=0
      end if
    end if
    do_x=.true.
    sigma=>sigmax
    call calc_delta_sigma_3d(adj,add,model,data)
  end subroutine
  logical function init_calc_sigma(s0,sx,sy,ds)
    real, target, dimension(:,:,:) :: s0,sx,sy
    real :: ds(3)
    d=ds
    sigmax=>sx
    sigmay=>sy
    slow0=>s0
    n1h=size(slow0,1)
    n3h=size(slow0,2)
    n2h=size(slow0,3)
    n=(/n1h,n2h,n3h/)
    init_calc_sigma=.true.
  end function
  integer function wrap_delta_sig(adj,add,model,data)
    logical :: adj,add
    real :: model(:),data(:)
    call adjnull(adj,add,model,data)
    do_x=.true.
    sigma=>sigmax
    wrap_delta_isg=calc_delta_sigma(adj,add,model,data)
  end function
  subroutine calc_delta_sigma_3d(adj,add,slow,sigma)
    logical :: adj,add
    real :: slow(n1h,n2h,n3h),sigma(n1h,n2h,n3h)
    real(kind=wp) :: dz(n1h,n2h,n3h)
    integer :: i1,i2,one,ia,izero,i3
    izero=floor(-o(1)/d(1)+1.5)
!First take care of the first term which is just a scaling of the
!slowness pertubation
    if (.not. adj) then
      sigma(izero,:,:)= sigma(izero,:,:)+2*slow(izero:,:,:)*sigma0&
        &(izero:,:,:)/slow0(izero:,:,:)
    else
      slow(izero:,:,:)= slow(izero:,:,:)+2*sigma(izero:,:,:)*sigma0&
        &(izero:,:,:)/slow0(izero:,:,:)
    end if
!
!#now we calcute a dz array which is the derivative of the integrand
!OLD
    dz=0.
    if (.not.adj) then
      dz(izero,:,:)=0
      do i1=izero+1,n(1)
        dz(i1,:,:)=dz(i1-1,:,:)-d(1)*.5*slow(i1-1,:,:)/slow0(i1-1,:,:)&
          &**2
      end do
    end if
!now we apply our changes to sigma
    if (do_x) then
      do i3=1,n(3)
        do i1=izero,n(1)
          do i2=2,n(2)-1
            if (.not. adj) then
              sigma(i1,i2,i3)=sigma(i1,i2,i3)-slow0(i1,i2,i3)*(dz(i1,i2&
                &+1,i3)-dz(i1,i2-1,i3 ))/d(2)
            else
              dz(i1,i2+1,i3)= dz(i1,i2+1,i3)-slow0(i1,i2,i3)*sigma&
                &(i1,i2,i3)/d(2)
              dz(i1,i2-1,i3)= dz(i1,i2-1,i3)+slow0(i1,i2,i3)*sigma&
                &(i1,i2,i3)/d(2)
            end if
          end do
          if (.not. adj) then
            sigma(i1,1,i3)= sigma(i1,1,i3)-2.*slow0(i1,1,i3)/d(2)*&
              &(dz(i1,2,i3)-dz(i1,1,i3))
            sigma(i1,n(2),i3)= sigma(i1,n(2),i3)-2.*slow0(i1,n(2)&
              &,i3)/d(2)*(dz(i1,n(2),i3)-dz(i1,n(2)-1,i3))
          else
            dz(i1,2,i3)= dz(i1,2,i3)-2*slow0(i1,1,i3)/d(2) * sigma&
              &(i1,1,i3)
            dz(i1,1,i3)= dz(i1,1,i3)+2*slow0(i1,1,i3)/d(2) * sigma&
              &(i1,1,i3)
            dz(i1,n(2),i3)= dz(i1,n(2),i3)-2*slow0(i1,n(2),i3)/d(2&
              &) * sigma(i1,n(2),i3)
            dz(i1,n(2)-1,i3)= dz(i1,n(2)-1,i3)+2*slow0(i1,n(2),i3)&
              &/d(2) * sigma(i1,n(2),i3)
          end if
        end do
      end do
    else
      do i2=1,n(2)
        do i1=izero,n(1)
          do i3=2,n(3)-1
            if (.not. adj) then
              sigma(i1,i2,i3)=sigma(i1,i2,i3)-slow0(i1,i2,i3)*(dz(i1,i2&
                &,i3+1)-dz(i1,i2,i3-1))/d(3)
            else
              dz(i1,i2,i3+1)= dz(i1,i2,i3+1)-slow0(i1,i2,i3)*sigma&
                &(i1,i2,i3)/d(3)
              dz(i1,i2,i3-1)= dz(i1,i2,i3-1)+slow0(i1,i2,i3)*sigma&
                &(i1,i2,i3)/d(3)
            end if
          end do
          if (.not. adj) then
            sigma(i1,i2,1)= sigma(i1,i2,1)-2.*slow0(i1,i2,1)/d(3)*&
              &(dz(i1,i2,2)-dz(i1,i2,1))
            sigma(i1,i2,n(3))= sigma(i1,i2,n(3))-2.*slow0(i1,i2,n(3&
              &))/d(3)*(dz(i1,i2,n(3))-dz(i1,i2,n(3) -1))
          else
            dz(i1,i2,2)= dz(i1,i2,2)-2*slow0(i1,i2,1)/d(3) * sigma&
              &(i1,i2,1)
            dz(i1,i2,1)= dz(i1,i2,1)+2*slow0(i1,i2,1)/d(3) * sigma&
              &(i1,i2,1)
            dz(i1,i2,n(3))= dz(i1,i2,n(3))-2*slow0(i1,i2,n(3))/d(3&
              &) * sigma(i1,i2,n(3))
            dz(i1,i2,n(3)-1)= dz(i1,i2,n(3)-1)+2*slow0(i1,i2,n(3))&
              &/d(3) * sigma(i1,i2,n(3))
          end if
        end do
      end do
    end if
!slow=dz
    if (adj) then
      do i3=1,n3h
        do i2=1,n2h
          do i1=izero+1,n(1)
            do ia=i1-1,izero,-1
              slow(ia,i2,i3)= slow(ia,i2,i3)-.5*d(1)*dz(i1,i2,i3&
                &)/slow0(ia,i2,i3)**2
            end do
          end do
        end do
      end do
    end if
  end subroutine
end module
