	subroutine find_azm(i0,i1,npa,azx,aaox,azm,seazm,aaoim,seaaoim)
c
c ...... calculate az and aoi of phase onset
c ...... based on a Bob program
c
c
	real*4 azx(*),aaox(*)
	real*4 az(1000),aaoi(1000)
c
	npa= i1 - (i0 + 1) + 1
	npa2=npa/2
c
      do ii = 1, npa
        kk = i0 + 1 + ii - 1
        az(ii) = azx(kk)
        aaoi(ii) = aaox(kk)
      end do
c
c ...... find median and quartile limits for az & aaoi
c
	do 4 i=1,npa-1
	  do 3 j=i+1,npa
	    if(az(j).lt.az(i)) then
	      azt=az(i)
	      az(i)=az(j)
	      az(j)=azt
	    endif
	    if(aaoi(j).lt.aaoi(i)) then
	      aaoit=aaoi(i)
	      aaoi(i)=aaoi(j)
	      aaoi(j)=aaoit
	    endif
    3   continue
    4 continue
c
	npa2=npa/2
	if(npa2.lt.1) npa2=1
	azm=az(npa2)
	if(azm.lt.0.) then
	  azm=azm+360.
	elseif(azm.gt.360.) then
	  azm=azm-360.
	endif
	aaoim=aaoi(npa2)
	if(npa.ge.4) then
	  npa1=npa/4
	  npa3=(3*npa)/4
	  seazm=(az(npa3)-az(npa1))/2.
	  seazm=abs(seazm)
	  seaaoim=(aaoi(npa3)-aaoi(npa1))/2.
	  seaaoim=abs(seaaoim)
	endif
c
	return
	end
