subroutine interp
*
* interpolate model from old to new meshes
*
      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/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/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/terp/stpms,rat1
      common/contrl/ds,g,sm,wc,it,nite,ja,jb,j,k,l
      common/corats/ams(10),delmass(10),corat(10),dcorat(10),
     1              ncore, irdold, alph 

      if (nite.ge.20) niteo = 30
      s1=s2
      r1=r2
      b1=b2
      p1=p2
      t1=t2
      up1=up2
      ut1=ut2
      o1=o2
      ot1=ot2
      op1=op2
      e1=e2
      ep1=ep2
      et1=et2
      ea1=ea2
      q1=q2
      f1=f2
      w1=w2
      qe1=qe2
      qn1=qn2
      qnp1=qnp2
      qnt1=qnt2
*
* now that we've saved all the variables for the previous shell,
* let's go about stepping to the next shell.
* if irdold = 1 and we're near a mass where corat was defined on input, then
* reduce ds.
*
      fac = 1.
60    ds=c/dmax1(fac,q1,f1)
*
* 0.008 was the value I used before 2/28/91
* Matt uses .005 here. Gets lots more shells.
*
      dsx=dlog10(10.**s1+.008)-s1
*
* Matt's new value as of 2/28/91 
* want -1.0 instead of -4.0 to  prevent undesirable jumps in core
*
      if(s1.gt.-1.0.and.s1.lt.-.02)then
        ds=dsx
      endif
      if( nite .le. 20 .and. nite .ge. 10 
     1 .or.(nite.eq.0.and.niteo.le.20.and. niteo .ge. 10)) then
	 ds = .90*ds
	 niteo = nite
      elseif ( nite .lt. 10 .and. nite .ge. 2 
     1 .or.(nite.eq.0.and.niteo.lt.10.and. niteo .ge. 2)) then
	 ds = .80*ds
	 niteo = nite
      endif
      s2=s1+ds
*
* make sure that corat boundaries are mesh points 
*
      if ( irdold .eq. 1 ) then
         do 100 i = 2,ncore-1 
            qi = dlog10( ams(i) )
            if ( s1 .lt. qi  .and. s2 .ge. qi ) then
               s2 = qi
               ds = qi - s1
               goto 200
            endif
100      continue
      endif
 
200   if ( s2 .ge. stpms ) then
         it=1
         s2=stpms
         ds=stpms-s1
      endif
 
415   if ( j .le. jb ) then
         if ( s2 .le. s(j) ) then
            sfrac=(s2-s(j))/(s(j)-s(j-1))
            p2=p(j)+(p(j)-p(j-1))*sfrac 
            b2=b(j)+(b(j)-b(j-1))*sfrac 
            r2=r(j)+(r(j)-r(j-1))*sfrac 
            t2=t(j)+(t(j)-t(j-1))*sfrac 
            return
         endif
         j=j+1
         go to 415
      endif
 
      cs=dlog10(1.+ds/s1)
      p2=p1+cs
      t2=t1+w1*cs
      r2=r1+ds*q1
      b2=b1
 
      return
      end 

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