subroutine pulsate
************************************************************************
* This code is taken from PRPWDXDP.F the prep code for models produced *
* by the evolution code WDXDFIT.F last modified by Mike Montgomery. It *
* calls the subroutine pulse at the end which is taken from CJHANRO.F  *
* a pulsation code originally written by Carl Hansen.                  *
* Specialized for Metacomputer use by Travis Metcalfe, 1998/1999       *
************************************************************************
      implicit double precision(a-h,o-z)

      common/freq/partt(650),piece(650),acous(650),bvfreq(650),
     1 bvfrq(650),tfreq(650)
      common/element/xhe(650)
      COMMON/modelp1/gnu0,per0,tdyn,mstar,model
      COMMON/modelp2/age,llsun,rrsun,teff,np_tmp
      COMMON/tape28a/brad,amass,nmod,nsurf
      COMMON/tape28b/xi2,rn,grav,rho,mr2
      COMMON/tape29/y1,voga1,y3,u,y5

      real*8 r(650),lr(650),t(650),rho(650),p(650),xx(21,650),
     1  kap(650),cv(650),chr(650),cht(650),epsr(650),epst(650),
     2  kapt(650),del(650),delad(650),eps(650),kapr(650)
      real*8 xi(650),u(650),v(650),voga1(650),ra(650),ga1(650),
     1  sound(650),derro(650),chtor(650),ram(650)
      real*8 mr(650),mstar,llsun,lnlsun,rrsun,teff
      real*8 gnu0,per0,tdyn
      real*8 xi2(650),rn(650),grav(650),mr2(650)
      real*8 x,g,dens
      real*8 y1(650),y2,y3(650),y4,y5(650)

*** stuff for Wojtek's pulsation code ***
      real*8 mu(650),akrdsig2(650),tfreq
*
* read evolved model file
*
      open(35,file='evolved.mod',status='old')

      read(35,*) model,age,pcen,tcen,ucen,rstr,teff,llsun,lnlsun,xtal

      rrsun=10.**rstr/6.96e+10
      llsun=10.**llsun
      teff=10.**teff

      read(35,1000) np
1000  format (i5)
      np_tmp = np
*
* read in equilibrium quantities
*
      do 50 i=1,20
        in=i
        if(i.eq.8) in=9
        if(i.eq.9) in=10
        if(i.eq.10) in=11
        if(i.eq.11) in=13
        if(i.eq.13) in=14
        if(i.eq.14) in=15
        if(i.eq.15) in=16
        if(i.eq.16) in=17
        if(i.eq.17) in=18
        if(i.eq.18) in=8
        if(i.eq.19) in=20
        if(i.eq.20) in=21
        read(35,8001) (xx(in,n),n=1,np)
50    continue
8001  format(1p,4e22.15)
      close(35)

      mstar=xx(2,np)/1.989e+33
      xx(2,np-1)=xx(2,np-2)
      xx(2,np)=xx(2,np-1)
*
* take 19 slices through 2d array
*
      do n=1,np

            r(n)=xx( 1,n)
           mr(n)=xx( 2,n)
           lr(n)=xx( 3,n)
            t(n)=xx( 4,n)
          rho(n)=xx( 5,n)
            p(n)=xx( 6,n)
          eps(n)=xx( 7,n)
          kap(n)=xx( 8,n)
           cv(n)=xx( 9,n)
          chr(n)=xx(10,n)
          cht(n)=xx(11,n)
         epsr(n)=xx(12,n)
         epst(n)=xx(13,n)
         kapr(n)=xx(14,n)
         kapt(n)=xx(15,n)
          del(n)=xx(16,n)
        delad(n)=xx(17,n)
          xhe(n)=xx(18,n)
        derro(n)=xx(19,n)

      enddo
*
* start computing pulsation quantities
*
      g=6.67259e-8 
      pi=3.14159265358979
      pi4=4.d0*pi 

      np=np-1
      totr=r(np) 
      totm=mr(np)

      do 40 n=2,np
        xi(n)=dlog(r(n)/p(n))
        chtor(n)=cht(n)/chr(n)
        v(n)=g*rho(n)*mr(n)/p(n)/r(n)
        ga1(n)=chr(n)+p(n)*cht(n)**2/rho(n)/t(n)/cv(n)
        voga1(n)=v(n)/ga1(n)
        u(n)=pi4*rho(n)*r(n)**3/mr(n)
        ram(n)=v(n)*chtor(n)*(del(n)-delad(n)-xx(20,n))
        bvfrq(n)=-g*mr(n)*ram(n)/r(n)**3
*** modified Ledoux case ***
        ra(n)=ram(n)
        bvfreq(n)=bvfrq(n)
        r(n)=r(n)/totr
        mr(n)=mr(n)/totm 
40    continue

      call asymp(np,r,sound,totr,totm,akrdsig2,mu,mr,rho,gnu0,per0,tdyn)
*
* below is everything needed for the pulse code
*
      nsurf=np-1
      nmod=model
      amass=totm
      brad=dlog10(totr)

      do 101 i=2,np
        x=xi(i)
        xi2(i) = xi(i)
        rn(i)=totr*r(i)
        mr(i)=mr(i)*totm
        mr2(i) = mr(i)
        grav(i)=g*mr(i)/rn(i)**2
        dens=rho(i)
        y1(i)=grav(i)/rn(i)
        y2=voga1(i)
        y3(i)=-1.*ra(i)
        y4=u(i)
        y5(i)=1./(1.+v(i))
101   continue

      call pulse

      return
      end 

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