subroutine read2 * * read model and some more control parameters * also contains write1, write2, and write3. * implicit double precision(a-h,o-z) common/shells/ sa(400),ra(400),ba(400),pa(400),ta(400), 1 ea(400),xca(400),fca(400),s(400),r(400),b(400), 2 p(400),t(400),e(400),xc(400),sk(400), 3 rk(400),bk(400),pk(400),tk(400) common/thermo/u2,up2,ut2,e2,ep2,et2,psi2,pg,o2,op2,ot2,fp,ft, 1 en2(5),fn,fcc,ce,ci,cif,w,nu common/surf/u(2,3),v(2,3),ww(2,3),rm,bm,is,ks,ls,ms,kstart,npass common/temp/s1,r1,s2,r2,b2,p2,t2,ea2,xc2,xo2,fca2,f2,q2,w2,c common/contrl/ds,g,sm,wc,it,nite,ja,jb,j,k,l common/xx/sin,sout,smid,grid1,grid2,ucent,cste,kon,kom common/times/tmax1,tmax2,tmax3,time,time1,f,dg,dg1,g1,ip5,ip6,ip8 common/terp/stpms,rat1 common/vca/modnr common/dvca/sgtmp common/xbnd/amxc,amxo,gold common/ixswch/ixswch common/corats/ams(10),delmass(10),corat(10),dcorat(10), 1 ncore, irdold, alph common/fin/ifin common/dfin/rfin,blfin common/tfit/tfit(100),told,ifit,itfit,nite1 common/crash/itcrashed,irestart dimension bn(5) dimension x(400,7),y(400,7),z(400,5) dimension cp(400), cv(400) equivalence (sa(1),x(1,1)),(s(1),y(1,1)),(sk(1),z(1,1)) * * Format Statements for READ2 * 1 format(12i4) 2 format(f6.3,f10.6,e14.7,e10.3,f5.2,2e9.2,i4,f9.6) 3 format(2e16.8,4f8.5,i4,f5.2) 4 format(7f7.5) 7 format(i10/(e15.8,f9.6,1pe11.4,0pf10.6,f9.6,e14.7,f8.6)) 71 format(i10/(e15.8,f9.6,1pe11.4,0pf10.6,f9.6,e14.7,f8.6,f8.6)) 8 format(4i10/2f10.6/(6d13.7)) 20 format(e15.8,f9.6,e15.8,f9.6,2e15.8) 21 format(2f9.3,f10.6) 103 format(///20h **** model nr. ,i4,7h *** ,f6.3,43h solar mass 1es (h - he - c/o )) *** age = ,1pd14.8,11h years ****/) 104 format(//6h ip5=i4,6h ip6=i4,6h ip7=i4,7h nmod=i4,8h nite1= 1i4,4h m=i4,5h md=i4,5h nn=i4,6h ip8=i4,6h w=f8.3) 114 format(60x,4hsin=f7.4,6h sout=f7.4,7h smid=f7.4,7h grid1=f7.4, 1 7h grid2=f7.4/) 105 format(4h dg=f10.6,4h g=f10.6,4h sg=e16.8,7h stpms=e11.4, 1 4h f=f10.6,8h tmax1=e10.4,8h tmax2=e10.4//) 106 format (5h sm=e16.8,6h rat1=e16.8,5h c=f6.3,6h xh=f9.5, 1 6h yh=f9.5,6h yyh=f9.5,6h ja=i4//) 153 format(e15.8,f9.6,1pe11.4,0pf10.6,f9.6,e14.7,2f8.6) 160 format(f6.3,f10.6,e14.7,e10.3,f5.2,2e9.2,i4,f9.6) c162 format(1h ,i2,f6.2,e13.6,6f8.3,7x,32hnn,sm,sg,p2,t2,ucent,rm,bl, c 1tel) 162 format(1h ,i2,f6.2,e13.6,6f8.3,7x,30hnn,sm,sg,p2,t2,ucent,rm,bl,tel) 163 format(1h ,'amx=',f6.3,' amxo=',f6.3) * * Format statements for WRITE2 * 309 format(42h **** convection in mass-shells j=i3,7h, s(j)= 1e16.8,12h, through j=i3,7h, s(j)=e16.8) 310 format(26h containing f11.5,25h solar masses with 1x12=f8.6,' x16=',f8.6) 319 format(/' s2=',1pe16.8,' r2=',0pf12.7,' b2=', 1pe12.4, 1 ' p2=', 0pf11.6,' t2=',f11.6,' e2=',1pe14.7, 2 ' xc2=',0pf8.6,' xo2=',f8.6) 321 format(' q2=', 1pe12.4, ' f2=', e12.4, ' w2=', e12.4,' ea2=', 1 e12.4,' psi=',e12.4,' cp= ',e12.4,' cv = ',e12.4) 323 format(4h u2= e12.4,5h o2= e12.4,4h b=f11.6,5h wc=f11.6, 1 5hvesc=e12.4,9h eratio=e12.4) 325 format(5h fte=e12.4,5h fcc=e12.4,5h fn=e12.4,23x,2hr=e12.6, 15h m=f9.6) 326 format(1h ,6h bnpl=f8.4,6h bnpa=f8.4,6h bnph=f8.4,6h bnrc=f8.4, 16h bnbr=f8.4) 327 format(5h epl=1pe12.4,5h eph=e12.4,5h epa=e12.4,5h erc=e12.4, 1 5h ebr=e12.4,6h etot=e12.4) * * Format statements for WRITE3 * 711 format(/7h fcc=e12.4,6h, s2=f8.3 ,7h, xdel=f12.6, 1 26h max.c12 rate, x12-change) 712 format(7h fn=e12.4,6h, s2=f8.3 ,23h max.neutrino rate) 713 format(7h bl=f12.4,6h, bcc=f12.4,6h, bn=f12.4,7h,bgrav=f12.4, 1 8h egrav=e14.5/) 724 format(2h f7.3,e14.5,3e14.4,23h m,r,v,a,shock at peak) 725 format(9h e14.5,2e14.4,35h r,v,a at last sh 1ell/) 333 format(8(1x,f8.4)) * * subroutine TAPE reads in carbon and oxygen interior eos table * call tape * * read in the model, and set some variables * note that the model is now tagged on the end of the input file * * * read from input file * read(7,*)ip5,ip6,ip7,nmod,nite1,m,md,nu,ip8,ip1,irdold,ip40 if ( irdold .eq. 1 .or. irdold .eq. 0) then iii = 0 888 iii = iii + 1 * * read from input file * read(7,*) ams(iii), corat(iii) if ( iii .gt. 1 ) then dcorat(iii) = corat(iii) - corat(iii - 1) delmass(iii) = ams(iii) - ams(iii-1) endif if ( ams(iii) .lt. 1.0 ) goto 888 ncore = iii elseif ( irdold .eq. 2 ) then * * read from input file * read(7,*) alph endif * * read from input file * read(7,3) sm,rat1,c,xh,yh,yyh,khomo,time read(7,4) ce,cif,sin,sout,smid,grid1,grid2 read(7,2) dg,g,sg,stpms,f,tmax1,tmax2,modnr,g3 * * hardwire time step * c dg = 13.8 * * change time step for restarted jobs * dg = dg - 0.1*float(irestart) tmax3=10.*tmax2 gold=g * * read from input file * read(7,20) speak,rpeak1,vpeak1,rlast,vexp1,egrav1 read(7,21) bgrav,w,g1 read(7,8) ks,ls,ms,kstart,rm,bm,u,v,ww * * here is where we read in the model. if we are reading in * an old model, then we also have to specify the run of the * desired c/o profile. (corat = 1 ==> pure c). * if irdold = 1, then do linear interpolation between the specified * mass points * if irdold = 2, then use the relation given in Barrat, Hansen, * & Mochkovitch 1988, A+A, 199, L15 * if ( irdold .eq. 1 ) then * * read from input file * read(7, 7) jb,((y(j,i),i=1,7),j=1,jb) do 1000 i = 1,jb amr = 10.**y(i,1) iii = 2 999 if ( amr .gt. ams(iii) ) then iii = iii + 1 goto 999 endif xc2 = dcorat(iii) / delmass(iii) * (amr - ams(iii-1)) 1 + corat(iii-1) xo2 = 1 - xc2 xc(i) = xc2 if ( xc2 .gt. .00001 .and. xc2 .lt. .99999)then call istatco(p(i),t(i),0,.false.) elseif ( xc2 .ge. .999999) then call istat1(p(i),t(i),0,12,.false.) elseif ( xc2 .le. .000001) then call istat1(p(i),t(i),0,16,.false.) endif e(i) = e2 1000 continue elseif ( irdold .eq. 2 ) then * * read from input file * read(7, 7) jb,((y(j,i),i=1,7),j=1,jb) do 1001 i = 1,jb amr = 10.**y(i,1) if ( i .ne. jb ) then xo2 = .5 * ( 1. + alph ) * ( 1. - amr )**alph else xo2 = 0. endif xc2 = 1 - xo2 xc(i) = xc2 if ( xc2 .gt. .000001 .and. xc2 .lt. .999999) then call istatco(p(i),t(i),0,.false.) elseif ( xc2 .ge. .999999) then call istat1(p(i),t(i),0,12,.false.) elseif ( xc2 .le. .000001) then call istat1(p(i),t(i),0,16,.false.) endif e(i) = e2 1001 continue corat(1) = xc(1) corat(2) = xc(jb) ncore = 2 else * * read from input file * read(7, 7) jb,((y(j,i),i=1,7),j=1,jb) do j=1,jb xc(j)=y(j,7) enddo endif close(7) smass=10.**(sm-33.298635) call armove(x,y,jb,7) * * perform a homology transform if requested, store new model in name.out * if(khomo.gt.0) then call homo(xh,yh,yyh,jb) endif ja=jb time1=time dg1=dg if ( w.lt.1. ) ip8=-1 npass=0 return *......................................................................* entry write1 * * write a (header) summary of a model, then set some * variables to start a new sequence. * nite=nite1 modnr=modnr+1 sg=sg+10.**(g-7.499095) smass=10.**(sm-33.298635) fcc1=0. xdel=0. sc1=0. fn1=0. sn1=0. k1=0 kon=0 kom=0 total=0. tmass=0. egrav=0. bcc=1e-10 do 207 i=1,5 bn(i)=1e-10 207 continue bnt=1e-10 cste=c g2=dg-g3 do 208 j=1,jb xc2 = xc(j) xo2 = 1 - xc2 if ( xc2 .gt. .000001 .and. xc2 .lt. .999999) then call istatco(p(j),t(j),0,.false.) elseif ( xc2 .ge. .999999) then call istat1(p(j),t(j),0,12,.false.) elseif ( xc2 .le. .000001) then call istat1(p(j),t(j),0,16,.false.) endif fca(j)=fcc/((1.-2.2892*fcc*10.**(g-18.)/0.999999)**2) 208 continue * * fca is fractional C abundance at time t + delta t * See eqn 3 of Kutter & Savedoff 1969, ApJ, 157, 1021. * return *......................................................................* entry write2 * * called by end, write2 writes a shell-by-shell summary of * the computed model * * compute specific heats cp and cv, and store in * cp() and cv(). write out both here and along with * the model (to tape9). * cp(k) = et2 cv(k) = et2 - ep2 * ut2 / up2 xmass=10.**s2-10.**s1 smr=10.**(sm-33.591065) xr=xmass*smr do 330 i=1,5 bn(i)=bn(i)+en2(i)*xr 330 continue bnt=bnt+fn*xr egrav=egrav-xmass*10.**(2.*sm-7.176004)*(10.**s2+10.**s1)/ 1 (10.**r2+10.**r1) * * compute radius, velocity, acceleration, and shock at location of max. * expansion. (At least that's what it used to do.) * Now it computes this for m/M* = 0.1 (fixed) * if ( speak .le. s2 .and. s1 .le. speak ) then rpeak=r1+(r2-r1)*(speak-s1)/(s2-s1) vpeak=(10.**rpeak-10.**rpeak1)*10.**(-g) apeak=2.*(vpeak-vpeak1)/(10.**g+10.**g1) shock=apeak*10.**(rpeak*2.-speak-sm+7.176004) endif * * mess with core convection zone, if one is present. * if(xc2.ne.0.) xctmp=xc2/(1.+2.2892*fcc*10.**(g2-18.)/xc2) if ( xc(k) .lt. 0.000001 ) xc(k)=0.000001 if ( kon .gt. 0 ) then k1=k tmass=tmass+xmass total=total+xctmp*xmass elseif ( k1 .gt. 0 ) then do 308 k3=kom,k1 xctmp=total/tmass 308 continue xotmp = xo2 tmass=tmass*10.**(sm-33.2986) k1=0 kom=0 total=0. tmass=0. endif * * carbon burning stuff * if ( fcc .gt. fcc1 ) then fcc1=fcc xdel=xc(k)-xc2 sc1=s2 endif if ( fn .gt. fn1 ) then fn1=fn sn1=s2 endif if(it.le.0 .and. mod(modnr,ip1).ne.0 1 .and. mod(k,10).ne.1)then if ( ip40 .eq. 1 ) then goto 1980 else return endif endif pz=10.**(pg-p2) eratio=10.**(u2+s2+sm-p2-r2-7.352095)/(2.-pz) vesc=10.**(.5*(s2+sm-r2-6.874974)) fte=(e2-ea2)*10.**(t2-g) r3=10.**r2 s3=10.**s2 theta = 7832.3 * dsqrt( 10.**u2 ) / 10.**t2 gammac = 10.**(u2/3. -t2) * 3.5772e6 gammao = gammac / 3.5772e6 * 5.7782e6 gamma = xc2 * gammac + xo2 * gammao * * before exiting, write to tape40 if ip40=1 * 1980 if ( ip40 .eq. 0 ) return bl=dlog10(b(jb)) tel=bl/4.-rm/2.+9.18458 if ( k .eq. 1 ) then endif return *......................................................................* entry write3 * * write header and model to tape9. * sn1=10.**sn1 sc1=10.**sc1 bl=dlog10(b(jb)) blfin=bl rfin=rm bcc=dlog10(bcc) do 331 i=1,5 bn(i)=dlog10(bn(i)) 331 continue bnt=dlog10(bnt) * * grav. potential luminosity * bgrav=dlog10(dabs(egrav-egrav1))-g-33.591065 * * old expansion quantities. velocity and acceleration * vexp=(10.**r(jb)-10.**rlast)*10.**(-g) aexp=2.*(vexp-vexp1)/(10.**g+10.**g1) spek=10.**speak rpek=10.**rpeak rlast=r(jb) rlas=10.**rlast * * dump out neutrino rates to separate file * rename variables to those of previous model. * egrav1=egrav vexp1=vexp vpeak1=vpeak rpeak1=rpeak g1=g g=g2 g3=.5*dmax1(bl,bcc) p2=p(1) t2=t(1) tel=bl/4.-rm/2.+9.18458 if ( ixswch .eq. 1 ) then pfreeze = phasec(t2) elseif ( ixswch .eq. 2 ) then pfreeze = phaseo(t2) elseif ( ixswch .eq. 3 ) then pfreeze = xc(1) * phasec(t2) + (1.-xc(1)) * phaseo(t2) elseif ( ixswch .eq. 4 ) then pfreeze = dmin1( phasec(t2), phaseo(t2) ) else stop endif xtest = y(1,4) - pfreeze amxc = -10. amxo = -10. if ( xtest .ge. 0. ) then call xbndry endif gold=g if ((ixswch.ne.4).and.(tfit(ifit).lt.1.d-3)) then open(50,file='evolved.mod',status='unknown') write(50,161) modnr,sg,p2,t2,ucent,rm,tel,bl,bnt,10.**amxc smx=10.**(sm-33.298635) endif 161 format(i4,1p,e12.4,0p,f7.3,2f6.3,f7.3,f9.6,2f9.4,2f7.3) 261 format(f4.1,2i5,1p,3e9.2,0p,f8.0,f9.3,f8.3) if ( ip7 .gt. 0 ) then m=m-1 endif nmod=nmod-1 if ( nmod.le.0 )then nshint = jb call eprep(nshint) stop endif return end ************************************************************************