# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module avomodel  
! computes avo model
  use adj_mod

  implicit none
  real, private		 			:: o2, d2
  real,private,dimension(:),allocatable		:: mask
  integer, private				:: n1,n2,n3

  contains

  subroutine avomodel_init(n1_in,n2_in,n3_in,o2_in,d2_in,mask_in)
!  subroutine avomodel_init(n1_in,n2_in,o2_in,d2_in)
    real 			:: o2_in, d2_in
    integer			:: n1_in, n2_in,n3_in
    real,dimension(n1_in*n2_in)	:: mask_in

    n1 = n1_in
    n2 = n2_in
    n3 = n3_in
    o2 = o2_in
    d2 = d2_in
    allocate(mask(n1*n2*n3))
    mask = mask_in
  end subroutine

  function avomodel_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 avomodel_lop2(adj,add,xx,yy )
    stat=0
  end function 

  subroutine avomodel_lop2(adj,add,xx,yy)
    logical,intent(in)   :: adj,add 
    real, dimension (:)  :: xx 
    real, dimension (:)  :: yy 
    real pi
    integer i, j       
    real theta
   
    do j = 1, n2
	theta=o2 + (j-1)*d2
	do i = 1, n1
	if ( adj) then
		yy((j-1)*n1+i)=yy((j-1)*n1+i)*mask((j-1)*n1+i)
		xx(i)=xx(i)+yy((j-1)*n1+i)
		xx(n1+i)=xx(n1+i)+yy((j-1)*n1+i)*sind(theta)**2
		xx(2*n1+i)=xx(2*n1+i)+yy((j-1)*n1+i)*sind(theta)**2*tand(theta)**2
	else
		yy((j-1)*n1+i)=yy((j-1)*n1+i)+xx(i)+xx(n1+i)*sind(theta)**2+xx(2*n1+i)*sind(theta)**2*tand(theta)**2
		yy((j-1)*n1+i)=yy((j-1)*n1+i)*mask((j-1)*n1+i)
	end if
	end do
    end do
  end subroutine 

  subroutine avomodel_close()
    deallocate(mask)
  end subroutine avomodel_close

end module avomodel
