module oc_sos_mod
  !!=======================
  !! Paul Sava    (C) 2000
  !! paul@sep.stanford.edu
  !! oclib v2.0
  !!=======================
  use oc_global_mod
  use oc_file_mod
  use oc_adjnull_mod
  use util_

  implicit none

  complex, dimension(:,:,:,:), pointer, private :: xx,yy
  type(myaxis),                         private :: amx,amy,ahx,ahy,az
  integer :: ihx0,ihy0
  integer :: bites

contains
  !----------------------------------------------------------------
  subroutine oc_sos_init(t_)
    character(len=128), intent(in) :: t_
    integer                        :: er
    type(fileinfo)                 :: file

    call oc_allocateheader(file,t_)
!!$    call oc_infofile(file)

    amx%n=file%n(1); amx%o=file%o(1); amx%d=file%d(1)
    amy%n=file%n(2); amy%o=file%o(2); amy%d=file%d(2)
    ahx%n=file%n(3); ahx%o=file%o(3); ahx%d=file%d(3)
    ahy%n=file%n(4); ahy%o=file%o(4); ahy%d=file%d(4)
    az%n =file%n(5); az%o =file%o(5); az%d =file%d(5)

!!$    call printaxis(amx)
!!$    call printaxis(amy)
!!$    call printaxis(ahx)
!!$    call printaxis(ahy)
!!$    call printaxis(az )

    allocate(xx(amx%n,amy%n,ahx%n,ahy%n),stat=er)
    if(er/=0) call seperr('cant allocate xx')
    allocate(yy(amx%n,amy%n,ahx%n,ahy%n),stat=er)
    if(er/=0) call seperr('cant allocate yy')

    ihx0=1+nint((0-ahx%o)/ahx%d)
    ihy0=1+nint((0-ahy%o)/ahy%d)

    bites = amx%n*amy%n*ahx%n*ahy%n *8
!write(0,*) bites
  end subroutine oc_sos_init
  !----------------------------------------------------------------
  integer function oc_sos(adj,add, x_,yy_, op1,op2,op3,op4,op5,op6,op7,op8,op9) result(stat)
    logical,            intent(in) :: adj,add 
    character(len=128), intent(in) :: x_,yy_
    integer, external,  optional   :: op1,op2,op3,op4,op5,op6,op7,op8,op9
    integer :: iz, imx,imy,ihx,ihy
    integer :: stat

!write(0,*) "prepare to seek"
    call sseek( x_,0,0)
    call sseek(yy_,0,0)
!write(0,*) "done seek"
    call oc_adjnull(adj,add,x_,yy_)
!write(0,*) "done adjnull"

    if(adj) then ! adjoint
!       write(0,*) "ADJOINT"
       do iz=1,az%n
          call sep_read(yy,yy_)
          call sep_read(xx, x_)
          stat=sseek_block(x_,iz-1,bites,0)

          xx=xx+yy
          do imy=1,amy%n
             do imx=1,amx%n
                xx(imx,imy,:,:) = &
                xx(imx,imy,:,:) - yy(imx,imy,ihx0,ihy0)
             end do
          end do
          call sep_write(xx,x_)
       end do
    else         ! forward
!       write(0,*) "FORWARD"
       do iz=1,az%n
          call sep_read(xx, x_)
          call sep_read(yy,yy_)
          stat=sseek_block(yy_,iz-1,bites,0)

          yy=yy+xx
          do imy=1,amy%n
             do imx=1,amx%n
                yy(imx,imy,ihx0,ihy0) = &
                yy(imx,imy,ihx0,ihy0) - sum(xx(imx,imy,:,:))
             end do
          end do
          call sep_write(yy,yy_)
       end do
    end if

    call sseek( x_,0,0)
    call sseek(yy_,0,0)
    stat=0
  end function oc_sos
  !----------------------------------------------------------------
end module oc_sos_mod





