subroutine intg
*
*     Runge-Kutta-Gill method of integration.  See Ralston & Wiff
*     Math Methods for Digital Computers  p.110  (Wiley,1960)
*
      implicit double precision(a-h,o-z)
      double precision k(4,5),m 

      dimension y(4,5),q(4,5),d(4,5)
      dimension ax(5),bx(5),cx(5)

      common/tenv/rhok(3,60,35),eta(3,60,35),pp(3,60,35),datg(3,60,35),
     1 od(3,18,35),tk(3,60),uu(3,60,35),xtt(3,60,35),xrt(3,60,35)
      common/shells/r(800),m(800),p(800),t(800),el1(800),ol1(800),
     1 rl1(800),atg(800),rtg(800),phs(800)
      common/c/az(3),z(3),taux,ts,ps,sm,rs,err,te,gs,rho,ol1x,rl1x,el1x
     1 ,atgx,rtgx,tl,pl,xmass,acu,am,bk,an0,pi
      common/d/ ml(3,60),nmax,kind,lmax,nsuff,nlim,nlim1,j,iter,lmam, 
     1 nlum,ko,nhp
      common/slr/sl(11)
      common/je/je
      common/comp/amhyhe,amheca,x(4)
      common/grad/ drdtp(400),cp(400),ttg(400),drdtpx,cpx,pecx,ttgx,
     1 deptx,dpdr,dtdr,dmdr,dr,xr,xt,cv,drdptx
      common/thresh/thresh
      common/mixl/aml,mls

      data ax/0.e0,5.e-1,2.928932190e-1,1.7071067810e0,1.6666666667e-1/ 
      data bx/0.e0,2.e0,1.e0,1.e0,2.e0/
      data cx/0.e0,5.e-1,2.928932190e-1,1.7071067810e0,5.e-1/

      dept=deptx
      if ( ko .le. 1 ) then
         q(1,1)=0.
         q(2,1)=0.
         q(3,1)=0.
         q(4,1)=0.
      else
         q(1,1)=q(1,5)
         q(2,1)=q(2,5)
         q(3,1)=q(3,5)
         q(4,1)=q(4,5)
      endif
 
      h=dr
      n=ko
      y(1,1)=r(n)
      y(2,1)=p(n)
      y(3,1)=t(n)
      y(4,1)=m(n)
      d(1,2)=1.
      d(2,2) =0.4342945*dpdr*10.**(-y(2,1))
      d(3,2)=0.4342945*dtdr*10.**(-y(3,1))
      d(4,2)=0.4342945*dmdr/(1.989e+33*10.**y(4,1))
 
      do 1 jx=2,4
         mx=jx-1
         kx=jx+1
         do 2 ix=1,4
            k(ix,jx)=h*d(ix,jx)
            y(ix,jx)=y(ix,mx)+ax(jx)*(k(ix,jx)-bx(jx)*q(ix,mx))
            q(ix,jx)=q(ix,mx)+3.*ax(jx)*(k(ix,jx)-bx(jx)*q(ix,mx))
     1               -cx(jx)*k(ix,jx)
2        continue
         pl=y(2,jx) 
         tl=y(3,jx) 
         xmass=10.**y(4,jx)

         call ccomp 
         call estatem
         ol1x = opac3(rl1x,tl)
         rtgx = 7.732749e+9*10.**(ol1x+pl-4.*tl)*sl(j)/xmass
         d(1,kx)=1. 
         d(2,kx)=-0.4342945*27398.*10.**(y(4,jx)-y(2,jx))*rho*(6.96e+10/
     1            y(1,jx))**2 
*
*     the mixing length is evaluated locally (shell n, halfway, shell n+1)
*
         if ( rtgx .ge. atgx ) then
            phsx=dabs(0.4342945/d(2,kx)) 
            dept=dept+dabs(y(1,mx)-y(1,jx))
            if ( phsx .gt. dept ) phsx = dept
            pecx = 1.601743e+07 * cpx * 
     1             10.**(ol1x+2.5*rl1x+y(4,jx)-3.*y(3,jx)
     2             -0.5*y(2,jx)) * dsqrt(dabs(drdtpx))
     3             *phsx**2*(6.96e+10/y(1,jx))**2 
            pecx=9./8./pecx
            pecx=pecx/aml/aml 
*
* pick between ml1 and ml3 versions of mixing-length theory 
*
            if(mls.eq.0)then
*
* ml1
*
               call gradt
             else
*
* ml3
*
               pecx=(8./9.)*3./4.242440687*pecx
               call gradt2
             endif
         else
            ttgx=rtgx
         endif
         d(3,kx)=ttgx*d(2,kx) 
         d(4,kx)=0.4342945*4.*pi*rho*y(1,jx)**2/(1.989e+33*10.**y(4,jx))
1     continue
 
      do 7 ix=1,4
         k(ix,5)=h*d(ix,5)
         y(ix,5)=y(ix,4)+ax(5)*(k(ix,5)-bx(5)*q(ix,4))
         q(ix,5)=q(ix,4)+3.*ax(5)*(k(ix,5)-bx(5)*q(ix,4))-cx(5)*k(ix,5)
7     continue
 
      r(n+1)=y(1,5) 
      p(n+1)=y(2,5) 
      t(n+1)=y(3,5) 
      m(n+1)=y(4,5) 
      return
      end 

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