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 ************************************************************************