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 

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