subroutine write4 * * print the structure of envelopes and convective transfer properties * both scaled down, also interpolate between shells to find exact * physical properties when eta=0,20, at the bottom of the convection * zone and when the mass of the envelope is=1e-6* stellar mass * implicit double precision(a-h,o-z) double precision m logical writeit common/shells/r(800),m(800),p(800),t(800),el1(800),ol1(800), 1 rl1(800),atg(800),rtg(800),phs(800) 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/je/je common/comp/amhyhe,amheca,x(4) common/xcomp/xcomp(400,3) common/slr/sl(11) 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/dirac/ psi(141),f12(141),f32(141) common/grad/ drdtp(400),cp(400),ttg(400),drdtpx,cpx,pecx,ttgx, 1 deptx,dpdr,dtdr,dmdr,dr,xr,xt,cv,drdptx common/surf/u(2,3),v(2,3),w(2,3),rm,bm,is,ks,ls,ms,kstart,npass common/terp/stpms,rat1 common/wprop/ bn2(400),ac2(400) common/writeit/writeit common/oldval/bold,teold common/mixl/aml,mls dimension s(400),d(400),nbn2(400),amfrac(400) * * Format statements * 12 format(/13h log(m/msun)=,f6.3,2x,13h log(r/rsun)=,f6.3,2x,11h log( 1grav)=,f6.3,2x,11h log(teff)=,f6.3,2x,13h log(l/lsun)=,f6.3,2x,12h 2 log(l*/m*)=,f6.3// 12h iterations=,i2, 359h nc=-1 if convv + nondeg, 0 if convv + deg, +1 if nonconvv/) 33 format(' too few shells') 37 format(4x,6hradius,5x,11hlog(1-r/rs),2x,11hlog(1-m/ms),2x,4hlogp, 1 4x,4hlogt,4x,6hlogrho,4x,7hopacity,5x,3heta,5x,3httg,8x,3hatg,8x, 2 3hrtg,6x,2hnc,3x,4hne/z) 900 format(/,4x,6hradius,5x,11hlog(1-r/rs),2x,11hlog(1-m/ms),10x,3hbn2 1,11x,3hac2,8x,2hje,4x,4hnbn2) 38 format(2x,1pe11.5,3x,0pf8.4,5x,f8.4,1x,2(1x,f7.3),2x,f7.3,3x, 1 1pe9.3,2x,0pf6.2,3(2x,e9.3),2x,i2,3x,f5.3) 902 format(2x,1pe11.5,3x,0pf8.4,5x,f8.4,2(2x,1pe16.9),0p,3f12.9,i5) 101 format(/1x,8henvelope,2x,9hstructure,16x,6hscale=,i2,5x, 1 ' total shells: ',i5) 103 format('log(1-r/rs)=',f8.4,2x,'log(1-m/ms)=',f8.4,2x,'logp=', 1 f7.3,2x,'logt=',f6.3,2x,'logrho=',f6.3,2x,'logkappa=',f6.3,2x, 2 'eta=',f6.2) 104 format(//57h physical conditions at the bottom of the convection z 1one/) 106 format(///////31h convective transfer properties,16x,6hscale=i2) 107 format(/4x,6hz=rs-r,5x,11hlog(1-r/rs),2x,11hlog(1-m/ms),3x, 1 7hlocal g,6x, 6hpeclet,6x,3httg,6x,5hsuper,5x,10hmix le 2ngth,2x,5hfc/ft,3x,9hconv velo,3x,5hv/vth) 108 format(2x,1pe11.5,3x,0pf8.4,5x,f8.4,4x,1pe9.3,4x,e9.3,3x,0pf5.3, 1 3x,1pe9.3,3x,e9.3,3x,0pf5.3,3x,1pe9.3,3x,0pf5.3) 109 format(//' physical conditions when m(env)/m(star)=',1pd14.6/) * * envelope structure * rads=r(1) nmod=1 alc4pig = -6.07666 - 2. * dlog10(gs) amass = 10.**sm * 1.989e33 do 4 n=1,nlim 987 format('m(n): ',f20.17) if( 1.-r(n)/rads .ge. 1.d-18 ) then d(n) = dlog10( 1.-r(n)/rads ) else d(n) = -10000. endif amfrac(n)=10.**m(n)/sm ame1=1.d0-amfrac(n) if (ame1 .ge. 1.d-18) then s(n)=dlog10(ame1) else s(n) = -10000. endif 4 continue xs=dlog10(sl(j)) xr=dlog10(rs) xg=dlog10(gs) xt=dlog10(te) xm=dlog10(sm) ratio=xs-xm ion=-1 do 11 n=1,nlim,nmod if ( atg(n) .gt. rtg(n) ) then nc=1 else if ( el1(n) .lt. 0. ) then nc=-1 else nc=0 endif endif opy=10.**ol1(n) * * compute degree of ionization xion * if ( ion .gt. 0. ) then xion=1.0 else zzz = xcomp(n,1)*z(1) + xcomp(n,2)*z(2) + xcomp(n,3)*z(3) xion=10.**(1.5*t(n)-rl1(n))*am/bk/zzz/1.329153 1 *f12(el1(n)) if ( xion .gt. 1.0) then xion=1.0 ion=1 endif endif 11 continue * * end of loop * * physical properties at places of interest * it=1 stop=dlog10(1.-10.**stpms) do 90 n=1,nlim l=n if ( s(n) .gt. stop) goto 91 90 continue 91 slope=(s(l-1)-stop)/(s(l-1)-s(l)) * * returns to 22 from below * 22 dfin=d(l-1)+(d(l)-d(l-1))*slope sfin=s(l-1)+(s(l)-s(l-1))*slope pfin=p(l-1)+(p(l)-p(l-1))*slope tfin=t(l-1)+(t(l)-t(l-1))*slope rfin=rl1(l-1)+(rl1(l)-rl1(l-1))*slope ofin=ol1(l-1)+(ol1(l)-ol1(l-1))*slope efin=el1(l-1)+(el1(l)-el1(l-1))*slope if ( it .le. 1 ) then stp = 10.**stop slum=sl(1) u(1,is)=pfin u(2,is)=tfin v(1,is)=dlog10(slum) v(2,is)=dlog10(1-10.**dfin)+w(2,is) it=it+1 * * skip over the next section with the following statement. * the next section computes the convective properties. * it's caused some problems in the hot models * if ( nlim1 .le. 0 ) go to 30 do 25 n=kind,nlim l=n if ( atg(n) .gt. rtg(n) ) goto 26 25 continue 26 slope=(atg(l-1)-rtg(l-1))/(atg(l-1)-rtg(l-1)-atg(l)+rtg(l)) go to 22 endif * * bottom of loop * * convection zone * nmod1=((nlim1-kind)/100)+1 * * compute geometrical depth zc, local gravity glocl, superadiabaticity * super, ratio of convective to total flux fcf, convective velocity * conv, ratio of convective to thermal velocities vcvt and peclet * number * do 52 n=kind,nlim1,nmod1 zc=rads-r(n) glocl=27398.*10.**m(n)*(6.96e+10/r(n))**2 super=ttg(n)-atg(n) fcf=1.-ttg(n)/rtg(n) if ( fcf .le. 0.) then fcf=0. conv=0. else conv=27398.*(te/5800.)**4*10.**(m(n)-t(n)-rl1(n)) 1 *3.90e+33*phs(n) 1 *drdtp(n)*fcf/(16.*pi*cp(n)*r(n)**2) conv=conv*aml * * choose mixing-length version here * (mls) =0 Bohm-Vitense or (mls=1) Bohm & Cassinelli * default is Bohm-Vitense version * if(mls.eq.1)then * * ml2 or ml3 version * conv=2.*dsqrt(2.d0)*conv endif if ( conv .gt. 0.0 ) conv = conv**0.333333 endif azave = xcomp(n,1)*az(1)+xcomp(n,2)*az(2)+xcomp(n,3)*az(3) vcvt=conv/(3.*bk*an0/azave*10.**t(n))**0.5 * * also need to adjust the peclet #. default is Bohm-Vitense. * peclet=cp(n)*phs(n)*conv/4.535688e-04*10.**(2.*rl1(n)+ol1(n)- 1 3.*t(n)) peclet=peclet*aml*aml * * Bohm & Cassinelli version * if(mls.eq.1)then peclet=peclet/dsqrt(2.d0) endif 52 continue 30 continue do 901 n=1,nlim,nmod nbn2(n)=bn2(n)/dabs(bn2(n)) if (bn2(n) .ge. 0.) then bn2(n)=dlog10(bn2(n)) else bn2(n) = -7. endif if (ac2(n) .ge. 0.) then ac2(n)=dlog10(ac2(n)) else ac2(n) = -7. endif 901 continue 1901 format(5(2x,f15.10)) return end ************************************************************************