      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 

************************************************************************

