!
!Mon Jan 13 21:20:40 PST 2003
!
module calc_sigma_mod
  use adj_mod
  use double_mod
  use sep3d_struct_mod
  use basic_util_mod
  implicit none
  real, private, save, pointer :: sigma0(:,:,:),sigmax(:,:,:),sigmay(:,:,:),slow0(:,:,:)
  real, private, save :: d(3),o(3)
  logical, private, save :: do_x
  integer, private, save :: nh1,nh2,nh3,n(3)
 
  contains


  logical function init_calc_sigma(slow_s)
    type(sep3d) :: slow_s,sigmax_s,sigmay_s
    integer,external ::sreed

    init_calc_sigma=.false.
    call set_n_o_d_3d(slow_s,n,o,d)
    nh1=n(1);nh2=n(2);nh3=n(3)

    allocate(slow0(nh1,nh2,nh3))
    call sseek("in",0,0)
    if(4*product(n)/=sreed("in",slow0,size(slow0)*4)) then
      write(0,*) "trouble reading in slowness0"
      return
    end if

    call init_sep3d("sigma_x",sigmax_s,"INPUT")
    allocate(sigmax(nh1,nh2,nh3))
    if(4*product(n)/=sreed("sigma_x",sigmax,size(sigmax)*4)) then
      write(0,*) "trouble reading in sigmax"
      return
    end if
  
    if(nh3/=1) then
      call init_sep3d("sigma_y",sigmay_s,"INPUT")
      allocate(sigmay(nh1,nh2,nh3))
      if(4*product(n)/=sreed("sigma_y",sigmay,size(sigmay)*4)) then
        write(0,*) "trouble reading in sigmax"
        return
      end if
    end if
    init_calc_sigma=.true.
  end function


  integer function calc_delta_sigmay(adj,add,model,data)
    logical,intent(in) :: adj,add
    real :: model(:),data(:)
    do_x=.false.;sigma0=>sigmay
    calc_delta_sigmay=0
    call adjnull(adj,add,model,data)
    call calc_delta_sigma_3d(adj,add,model,data)
  end function
   
  integer function calc_delta_sigmax(adj,add,model,data)
    logical,intent(in) :: adj,add
    real :: model(:),data(:)
    do_x=.true.;sigma0=>sigmax
    calc_delta_sigmax=0
    call adjnull(adj,add,model,data)
    call calc_delta_sigma_3d(adj,add,model,data)
  end function
   






  subroutine calc_delta_sigma_3d(adj,add,slow,sigma)
    logical :: adj,add
    real    :: slow(nh1,nh2,nh3),sigma(nh1,nh2,nh3)
    real(kind=wp)   :: dz(nh1,nh2,nh3)
    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



    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
    if(adj) then
      do  i3=1,nh3
        do  i2=1,nh2
          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
