subroutine calc(line) * * set up variables for the step to the next shell. * 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/temp/s1,r1,s2,r2,b2,p2,t2,ea2,xc2,xo2,fca2,f2,q2,w2,c common/cal/ut1,ot1,et1,op1,o1,ep1,qnp1,qnt1,p1,q1,f1,t1,w1,b1,e1, 1 qe1,ea1,qn1,up1,qe2,qn2,qnp2,qnt2 common/xx/sin,sout,smid,grid1,grid2,ucent,cste,kon,kom common/contrl/ds,g,sm,wc,it,nite,ja,jb,j,k,l common/je/je common/comp/amhyhe,amheca,x(4) common/corats/ams(10),delmass(10),corat(10),dcorat(10), 1 ncore, irdold, alph common/wprep/iprep common/dprep/aa(20,500), ecv(400), ext(400), exr(400) common/bldxc/bledc common/crash/itcrashed,irestart logical converg data ck/1.5e0/ data alsun/3.90e+33/ line = 2 je = 3 k=k+1 if ( nite .eq. 1 ) then converg = .true. else converg = .false. endif * * if memory full, then resubmit with new time step * if ( k .ge. 400) then itcrashed = 1 return endif * * otherwise, go on about our business. * sk(k) = s2 rk(k) = r2 bk(k) = b2 pk(k) = p2 tk(k) = t2 kk=k 420 if ( l .le. ja ) then if ( s2 .le. sa(l) ) then sfrac = (s2-sa(l))/(sa(l)-sa(l-1)) ea2 = ea(l)+(ea(l)-ea(l-1))*sfrac xc2 = xca(l)+(xca(l)-xca(l-1))*sfrac xo2 = 1.-xc2 fca2 = fca(l)+(fca(l)-fca(l-1))*sfrac else l=l+1 go to 420 endif elseif ( w1 .lt. 0.4 ) then ea2=e1-1.914535e+08*dlog10(1.+ds/s1)*(1.-2.5*w1)/cei endif * * because i'm now specifying the c/o profile in the interior, i'm * going to interpolate seperately here if irdold = 1 * as opposed to interpolating in the 'alda' model when irdold = 0,2. * this should keep the profile * from getting spread out through quantization error. note also * that in interp, ds is decreased if near a corat boundary, and * that things are zoned so that a mass shell will be on all * corat boundaries. * if ( irdold .eq. 1 ) then do 500 iii = 2,ncore test = 10.**s2 if ( test .eq. ams(iii) ) then xc2 = corat(iii) xo2 = 1. - xc2 goto 600 else if ( test .lt. ams(iii) ) then xc2 = dcorat(iii) / delmass(iii) * (test - ams(iii-1)) 1 + corat(iii-1) xo2 = 1. - xc2 goto 600 endif 500 continue elseif ( irdold .eq. 2 ) then if ( s2 .lt. s(jb) ) then xo2 = .5 * ( 1. + alph ) * ( 1. - 10.**s2 )**alph xc2 = 1. - xo2 else xo2 = 0. xc2 = 1. endif endif 600 c=cste xc(k) = xc2 if ( s2 .gt. -sin .and. s2 .lt. -sout ) then c=c*(dabs(smid+s2)*grid2+grid1) endif cei=1./(.5+1./ci) if ( xc2 .gt. .000001 .and. xc2 .lt. .999999) then call istatco(p2,t2,0,converg) elseif ( xc2 .ge. .999999) then call istat1(p2,t2,0,12,converg) elseif ( xc2 .le. .000001) then call istat1(p2,t2,0,16,converg) endif e(k)=e2 fcc=(1.-w)*fca2+w*fcc if(k.le.1) ucent=u2 q2=10.**(-3.*r2+s2-u2-1.099210+sm) f2=10.**(-4.*r2-p2+2.*s2-8.275214+2.*sm) w2=o2*b2*10.**(p2-4.*t2-s2+43.18734-sm) wc=-ep2/et2 if (w2 .ge. wc) then w2 = wc if (kon .le. 0) then c = c/ck kon = k endif else if ( kon .gt. 0 ) then kom = kon kon = 0 c = c * ck endif endif qx=10.**(s2-33.22885+sm) qn2=(fcc-fn)*qx qnp2=-fp*qx qnt2=-qx*ft qe2=qx*10.**(t2-g) if ( k .gt. 1 ) line = 1 * * Now load up the variables in aa() if we're calculating * the converged model * if ( converg .and. iprep .eq. 1 ) then aa( 1,k) = 10.**r2 star=10.0**(sm) aa( 2,k) = 10.**s2 * star aa( 3,k) = b2 * alsun aa( 4,k) = 10.**t2 aa( 5,k) = 10.**u2 aa( 6,k) = 10.**p2 aa( 7,k) = fn cpx = ut2/wc*10.**(-u2+p2-t2) xt = -ut2/up2 cv = cpx*(1.e0-(xt*wc)) aa( 8,k) = -cv aa( 9,k) = 1./up2 aa(10,k) = xt aa(11,k) = ft aa(12,k) = fp/up2 * * Paul's opac business (Will have to go to opac3 when comp is mixed) * only calling opac for carbon here, because we lack O table. * Rad opac contribution is zilch here anyway. * We use 5pt differencing in the core. * cent = 0.01 dent = 0.02 oplo4=opac(u2+dent,t2) oplo1=opac(u2+cent,t2) oplo2=opac(u2-cent,t2) oplo3=opac(u2-dent,t2) * * kapr * aa(13,k)=(oplo3-(8.*oplo2)+(8.*oplo1)-oplo4)/(12.*0.01) oplo2=opac(u2,t2+dent) oplo3=opac(u2,t2+cent) oplo4=opac(u2,t2-cent) oplo1=opac(u2,t2-dent) * * kapt * aa(14,k)=(oplo1-(8.*oplo4)+(8.*oplo3)-oplo2)/(12.*0.01) aa(15,k) = w2 aa(16,k) = wc aa(17,k) = 0.0d0 aa(18,k) = o2 aa(20,k) = xo2 * * this is core ledoux term. (set to zero, then reset if need be) * aa(19,k) = 0.0 * * if at center, set Ledoux term to 0 by definition. * if(k.eq.1)then return endif if(xc2.gt.0.000001 .and. xc2.lt.0.999999)then * * if O fraction is constant, set Ledoux to zero. * if(aa(20,k).eq.aa(20,k-1))then return else * * otherwise, compute a Ledoux term. * call ledouxc if(aa(20,k).le.0.0 .or. aa(20,k-1).le.0.0)then return endif dlo = ( dlog(aa(20,k)) - dlog(aa(20,k-1)) ) dlp = ( dlog(aa(6,k)) - dlog(aa(6,k-1)) ) chtor = -aa(9,k)/aa(10,k) aa(19,k) = chtor*(dlo/dlp)*bledc * * Note: there is a minus sign difference between this ledoux term 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,k) = -aa(19,k) endif endif endif return end ************************************************************************