! 
!
module res_o2a_mod

  use adj_mod

  implicit none
  real, private		 			:: o1,d1,o2,d2,o3,d3,oo2,do2,drs
  integer, private				:: n1,n2,n3,no2

  contains

  subroutine res_o2a_init(drs_in,n1_in,d1_in,o1_in,n2_in,d2_in,o2_in,n3_in,d3_in,o3_in,& 
	  no2_in,do2_in,oo2_in)
    real 			:: o1_in,o2_in,o3_in,d1_in,d2_in,d3_in,do2_in,oo2_in, drs_in
    integer			:: n1_in,n2_in,n3_in,no2_in

    n1 = n1_in;    n2 = n2_in;    n3 = n3_in
    o1 = o1_in;    o2 = o2_in;    o3 = o3_in
    d1 = d1_in;    d2 = d2_in;    d3 = d3_in
    no2=no2_in;    do2=do2_in;    oo2=oo2_in
    drs=drs_in
  end subroutine

  function res_o2a_lop ( adj, add, xx, yy) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    real,dimension(:,:)  :: xx,yy 
!    call adjnull (adj,add,xx,yy )
    call res_o2a_lop2(adj,add,xx,yy )
    stat=0
  end function 

  subroutine res_o2a_lop2(adj, add, xx, yy)
   integer 		:: ix,iz,it,ih, xstart, xend
   real			:: amp,t,z,z0,h, amax, xx(:,:), yy(:,:)
   real			:: dzdh,cosin,a,cosin15,tp
   real 		:: pi
   logical		:: add,adj

  pi=4*atan(1.) 
!  call adjnull( adj, add, xx,n1*n2, yy,n1*n2)


  do ig = 1,n2
    gama(ig) = o2 + (ig-1)*d2
  do ih= 1, no2
    h = do2 * (ih-1) + oo2				! h = subsurf.offset
    xstart = 1-ih;   xstart = max0( 1, xstart)
    xend = n2-ih;    xend = min0( n2, xend)
!    do ig= xstart, xend
!    write(2,*) ih,xstart, xend

    do iz= 1, n1
      zm = o1 + d1 * (iz-1) 			! z = depth
      z=sqrt(zm**2-h**2+h*drs)			! impulse response in SODCIG
      dzdh = (drs - 2*h)/(2*z)			! dz/dh of impulse response
      z0 = z - h*dzdh				! stacking curve in angle

      it = 1.5 + (t - o1) / d1			

      if ( it > n1 ) then
        exit
      end if

      a = 2*z / t
      cosin = cos(amax*pi/180.)
      cosin15 = cos((amax-15.)*pi/180.)
      amp = (2*z / t) * sqrt( 2*n1*d1 / t )
!
! Incude a linear taper
!
      if ( cosin15 <= a) then
        amp = amp
      else if ((a < cosin15) .and. (a > cosin)) then
        tp = (amax - acos(a)*180./pi)/15.
        amp = amp * tp
      else
        amp = 0.
      end if

      if (amp .ne. 0.) then
        if ( adj ) then
            xx(iz,ix )=xx(iz,ix )+yy(it,ix+ih)*amp
        else
            yy(it,ix+ih)=yy(it,ix+ih)+xx(iz,ix )*amp
        end if
      end if
    end do
  end do
  end do
  close(2)
  end subroutine
end module
