subroutine eprep(nshint) * * EPREP -- white dwarf carbon stratified prep code subroutine * writes prepped models to tape50 * implicit double precision(a-h,o-z) double precision m common/shells/r(800),m(800),p(800),t(800),el1(800),ol1(800), 1rl1(800),atg(800),rtg(800),phs(800) common/xcomp/xcomp(400,3) common/je/je common/comp/amhyhe,amheca,x(4) common/bldx/bled common/terp/stpms,rat1 common/temp/s1,r1,s2,r2,b2,p2,t2,e2,xc2,xo2,fca2,f2,q2,w2,c 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,jjj,iter,lmam, 1 nlum,ko,nhp common/wprep/iprep common/dprep/aa(20,500), ecv(400), ext(400), exr(400) common/grad/ drdtp(400),cp(400),ttg(400),drdtpx,cpx,pecx,ttgx, 1 deptx,dpdr,dtdr,dmdr,dr,xr1,xt1,cv1,drdptx common/wprop/ bn2(400),ac2(400) common/kapder/rkapr(800),rkapt(800) common/dfin/rfin,blfin common/opacfin/ifinal common/tfit/tfit(100),told,ifit,itfit,nite1 data amsun,alsun/1.989e+33,3.90e+33/ 987 format(a80) 1002 format(1x,4e16.9) 1003 format(4e16.9) 1004 format(1x,16i4) 1005 format(16i4) icmp = 0 small = 1.e-4 nl=nshint xc2 = 1. xo2 = 0. * * call env for the converged model * call env ifin = 0 nel=nlim if( iprep .eq. 0 )then * * no reason to continue here * return endif nel = nel - 1 do 34 n=1,nel index = nel-n+1 * * these are placed in reverse order * rrr = r(index) amm = m(index) ttt = t(index) rho = rl1(index) ppp = p(index) aa(1,nl+n)=r(index) x(1) = xcomp(index,1) x(2) = xcomp(index,2) x(3) = xcomp(index,3) if ( n .eq. nel ) then x(1) = 0. x(2) = 0. x(3) = 1. endif * * mass (grams) * aa(2,nl+n)=10.**m(index)*amsun * * temperature (kelvins) * aa(4,nl+n)=10.**t(index) * * density (g/cm^3) * aa(5,nl+n)=10.**rl1(index) * * pressure (dynes/cm^2) * aa(6,nl+n)=10.**p(index) * * energy generation rate + derivatives (eps, epsr, and epst) * aa(7,nl+n)=0.0 aa(11,nl+n)=0.0 aa(12,nl+n)=0.0 * * use delrad as temp. gradient except in convection zone, then use * deltrue from gradt (or gradt2). * if(rtg(index).lt.atg(index))then tdel=rtg(index) else tdel=ttg(index) endif * * temperature gradient (del) * aa(15,nl+n)=tdel * * adiabatic temperature gradient (delad) * aa(16,nl+n)=atg(index) * * helium profile * aa(17,nl+n)=xcomp(index,2) * * radiative luminosity profile * aa(3,nl+n)=alsun*(10.0**(blfin))*(ttg(index)/rtg(index)) * * heat capacity Cv * aa(8,nl+n)=ecv(index) * * pressure gradients chi rho and chi T * aa(9,nl+n)=exr(index) aa(10,nl+n)=ext(index) * * opacity * tl=t(index) pl=p(index) rhol = rl1(index) oplo = opac3(rhol,tl) aa(18,nl+n)=10.**oplo ifinal = 1 * * start 5-point opacity derivative stuff here. * if(x(1).eq.1.0 .or. x(2).eq.1.0 .or. x(3).eq.1.0)then if(x(1).eq.1.0)then je=1 endif if(x(2).eq.1.0)then je=2 endif if(x(3).eq.1.0)then je=3 endif call opacder(rhol,tl,opcs,dlklr,dlklt) aa(13,nl+n)=dlklr aa(14,nl+n)=dlklt * * note: replace previous opacity value if we can get it here. * aa(18,nl+n)=10.**opcs else * * 5 point differencing * oplo5=opac3(rhol+0.02,tl) oplo1=opac3(rhol+0.01,tl) oplo2=opac3(rhol-0.01,tl) oplo6=opac3(rhol-0.02,tl) oplo7=opac3(rho,tl+0.02) oplo3=opac3(rhol,tl+0.01) oplo4=opac3(rhol,tl-0.01) oplo8=opac3(rho,tl-0.02) * * opacity derivatives (kapr(13) and kapt(14)) * aa(13,nl+n)=(oplo6-(8.*oplo2)+(8.*oplo1)-oplo5)/(12.*0.01) aa(14,nl+n)=(oplo8-(8.*oplo4)+(8.*oplo3)-oplo7)/(12.*0.01) endif * * if He fraction is constant, set Ledoux to zero. * if(aa(17,nl+n).eq.aa(17,nl+n-1))then aa(19,nl+n) = 0.0 go to 34 endif if(aa(17,nl+n-1).eq.0.0 .or. aa(17,nl+n).eq.0.0)then aa(19,nl+n) = 0.0 go to 34 endif if(index.eq.nel)then aa(19,nl+n) = 0.0 else * * call up modified ledoux term subroutine * call ledoux dly = ( dlog(aa(17,nl+n)) - dlog(aa(17,nl+n-1)) ) dlp = ( dlog(aa(6,nl+n)) - dlog(aa(6,nl+n-1)) ) bledoux=(1./aa(10,nl+n))*bled*(dly/dlp) * * Note: there is a minus sign difference between this ledoux term * in the core and the one in Brassard et al. (1991, ApJ, 367, 601). * This was "corrected" by the absolute value below, but I want to be * able to see if certain configurations are convectively unstable in the * core, so I am going to remove the dabs() (MHM April 1998). * aa(19,nl+n)=bledoux if(x(1).ne.0.0 .and. x(2).gt.0.99 .and. x(3).ne.0.0)then aa(19,nl+n)= 0.d0 endif endif 34 continue ifinal = 0 * * now write it out. * nshell=nl+nel if (tfit(ifit).lt.1.d-3) then write(50,1010) nshell, nl, nel endif 1010 format(i5,' = ',i4,' + ',i4) * * smooth out del and delad at core-envelope boundary * set mver for appropriate variable * 1. delad; 2. del; 3. chiT; 4. chiRho * do 33 i=1,nshell if(i.eq.nl)then * * delad * diff = aa(16,i+1) - aa(16,i) if(diff.gt.0.002)then mver = 1 call smooth(nl,mver) endif * * del * diff2 = aa(15,i+1) - aa(15,i) if(diff2.gt.0.002)then mver = 2 call smooth(nl,mver) endif * * chi T * diff3 = aa(10,i+1) - aa(10,i) if(diff3.gt.0.002)then mver = 3 call smooth(nl,mver) endif * * chi rho * diff2 = aa(9,i+1) - aa(9,i) if(diff4.gt.0.002)then mver = 4 call smooth(nl,mver) endif endif 33 continue * * Write out model to tape50 * if (tfit(ifit).lt.1.d-3) then do 2 i=1,20 write(50,1020) (aa(i,n), n=1,nshell) 2 continue endif 1020 format(1p,4e22.15) if (tfit(ifit).lt.1.d-3) then close(50) endif return end ************************************************************************