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

  implicit none
  real, private		 	:: o2, d2
  integer, private, save 	:: mube,mufe,mubi,mufi
  integer, private		:: n2
  contains

  subroutine avomodel_init(n2_in,o2_in,d2_in,&
    & mube_in,mufe_in,mubi_in,mufi_in)
    optional 	:: mube_in,mufe_in,mubi_in,mufi_in
    real 	:: o2_in, d2_in
    integer 	:: mube_in,mufe_in,mubi_in,mufi_in, n2_in
    n2 = n2_in
    o2 = o2_in
    d2 = d2_in
    mube=0;mufe=9999;mubi=9999;mufi=0;
    if (mube_in.ne.0) mube=mube_in
    if (mufe_in.ne.0) mufe=mufe_in
    if (mubi_in.ne.0) mubi=mubi_in
    if (mufi_in.ne.0) mufi=mufi_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       
    real theta
    pi = 4*atan(1.)
    do i = 1, n2
	theta=pi*(o2 + (i-1)*d2)/180.
	if ( adj) then
		if ((i<mube).or.(i>mufe)) then
			yy(i)=0.
		elseif ((i>=mubi).and.(i<=mufi)) then
			yy(i)=0.
		end if
		xx(1)=xx(3)+yy(i)
		xx(2)=xx(2)+yy(i)*sin(theta)**2
		xx(3)=xx(3)+yy(i)*sin(theta)**2*tan(theta)**2
	else
		yy(i)=yy(i)+xx(1)+xx(2)*sin(theta)**2+xx(3)*sin(theta)**2*tan(theta)**2
		if ((i<mube).or.(i>mufe)) then
			yy(i)=0.
		elseif ((i>=mubi).and.(i<=mufi)) then
			yy(i)=0.
		end if
	endif
    end do
  end subroutine 
  subroutine avomodel_close()
  end subroutine 
end module 
