subroutine asymp(np,r,sound,totr,totm,akrdsig2,mu,mr,rho,
     &                  gnu0,per_0,tdyn)
*
* this subroutine computes the asymptotic frequency and period
* spacings for a given equilibrium model.
* Formulae taken from 5 horsemen.
* Added in computation for Pi_H and Pi_He for trapping cycles
*
      implicit double precision(a-h,o-z)

      common/freq/partt(650),piece(650),acous(650),bvfreq(650),
     1 bvfrq(650),tfreq(650)
      common/res/res(650)
      common/element/xhe(650)
      common/trap/per0,perH,perHe,alpha2,ane,l

      real*8 bvn(650),r(650),sound(650),h(650),gnu0,per0,per_0
      real*8 bvnh(650),bvnhe(650),mu(650),akrdsig2(650),torfreq0,
     1 akrdsig(650),avgne(650),rho(650),mr(650),dpxtal(650),mxtal(650)  
      data bvn,h,bvnh,bvnhe/2600*0.0d0/
*
* get sound speed and brunt-vaisala frequency
*
      nmax=0
      g=6.67259e-8 
      pi=3.14159265358979
      tdyn=dsqrt(g*totm/totr**3)
      do 10 i=2,np
*
* integrand for torsional frequency spacing (constant, like p-modes)
*
       if(akrdsig2(i).gt.0.0 .and. mr(i)/mr(np).lt..99)then
         akrdsig(i)=dsqrt(akrdsig2(i))
       else
         akrdsig(i)=0.0
       endif
*
* integrand for Pi_o
*
        if(bvfreq(i).gt.0.0)then
          bvn(i)=dsqrt(bvfreq(i))
          bvn(i)=bvn(i)/totr/r(i)
*
* integrand for Pi_H
*
          if(xhe(i).lt.0.97 .and. 
     1       xhe(i-1).gt.xhe(i))then
            bvnh(i)=dsqrt(bvfreq(i))
            bvnh(i)=bvnh(i)/totr/r(i)
	    if (ihflag.eq.0) ihyzone=i
	    if (i.lt.np) then
	      avgne(i)=(r(i+1)-1.)*log(rho(i+1)/rho(i))/(r(i+1)-r(i))
	    else
	      avgne(i)=0.0
	    endif
            ihflag = 1
           elseif(ihflag.eq.1)then
            bvnh(i)=dsqrt(bvfreq(i))
            bvnh(i)=bvnh(i)/totr/r(i)
	    if (i.lt.np) then
	      avgne(i)=(r(i+1)-1.)*log(rho(i+1)/rho(i))/(r(i+1)-r(i))
	    else
	      avgne(i)=0.0
	    endif
           else
            bvnh(i)=0.0
	    avgne(i)=0.0
          endif
*
* integrand for Pi_He
*
          if(xhe(i).gt.0.03)then
            bvnhe(i)=dsqrt(bvfreq(i))
            bvnhe(i)=bvnhe(i)/totr/r(i)
           elseif(xhe(i).lt.0.03 .and.
     1       xhe(i-1).gt.xhe(i))then
            bvnhe(i)=dsqrt(bvfreq(i))
            bvnhe(i)=bvnhe(i)/totr/r(i)
            iheflag = 1
           elseif(iheflag.eq.1)then
            bvnhe(i)=dsqrt(bvfreq(i))
            bvnhe(i)=bvnhe(i)/totr/r(i)
           else
            bvnhe(i)=0.0
          endif
*
* add zero to integrand in convection zone
*
         else
          bvn(i)=0.0
          bvnh(i)=0.0
          bvnhe(i)=0.0
	  avgne(i)=0.0
        endif
10    continue
      do 30 n=2,np
        sound(n)=1./sound(n)
        rs=r(n)*totr
        if(n.ne.1)then
          h(n)=rs-rm1
         else
          h(1)=r(1)*totr
        endif
        rm1=rs
30    continue
      npm1=np-1
      call dosim(rslt,akrdsig,h,npm1)
      torint=rslt
      torfreq0=(1./2.)/torint
      pertor=1./torfreq0
      call dosim(rslt,sound,h,npm1)
      cint=rslt
      call dosim(rslt,bvn,h,npm1)
      gint=rslt
      gnu0=1./(2.*cint)
      pgnu0=1./gnu0
      per0=2.*pi*pi/gint

      do i=3,npm1
        mxtal(i)=mr(i)/mr(np)
        if (mxtal(i).lt.0.98) then
          dpxtal(i)=2.*pi*pi/(res(npm1)-res(i))
        endif
1004  format(2(e12.5,2x))
      enddo

      if(ihflag .eq. 1)then
        call dosim(rslt,bvnh,h,npm1)
        ghint=rslt
        perH=2.*pi*pi/ghint
        call dosim(ane,avgne,h,npm1)
	dr=totr*(1.-r(ihyzone))
	ane=ane/dr
       else
        perH=0.0
      endif
      if(iheflag .eq. 1)then
        call dosim(rslt,bvnhe,h,npm1)
        gheint=rslt
        perHe=2.*pi*pi/gheint
       else
        perHe=0.0
      endif
      per_0 = per0
*
* call routine to do asymptotic period spacing with mode
* trapping using formulae in Brassard etal 1992, ApJS, 80, 369
*
      call modetrap

      return
      end

************************************************************************