      subroutine eprep(nshint)
*
* EPREP -- white dwarf carbon stratified prep code subroutine
*  writes prepped models to tape50
*
      implicit double precision(a-h,o-z)
      double precision m

      common/shells/r(800),m(800),p(800),t(800),el1(800),ol1(800),
     1rl1(800),atg(800),rtg(800),phs(800)
      common/xcomp/xcomp(400,3)
      common/je/je
      common/comp/amhyhe,amheca,x(4)
      common/bldx/bled
      common/terp/stpms,rat1
      common/temp/s1,r1,s2,r2,b2,p2,t2,e2,xc2,xo2,fca2,f2,q2,w2,c
      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,jjj,iter,lmam, 
     1 nlum,ko,nhp
      common/wprep/iprep
      common/dprep/aa(20,500), ecv(400), ext(400), exr(400)
      common/grad/ drdtp(400),cp(400),ttg(400),drdtpx,cpx,pecx,ttgx,
     1 deptx,dpdr,dtdr,dmdr,dr,xr1,xt1,cv1,drdptx
      common/wprop/ bn2(400),ac2(400)
      common/kapder/rkapr(800),rkapt(800)
      common/dfin/rfin,blfin
      common/opacfin/ifinal
      common/tfit/tfit(100),told,ifit,itfit,nite1

      data amsun,alsun/1.989e+33,3.90e+33/

987   format(a80)
1002  format(1x,4e16.9)
1003  format(4e16.9)
1004  format(1x,16i4)
1005  format(16i4)

      icmp = 0
      small = 1.e-4
      nl=nshint
      xc2 = 1.
      xo2 = 0.
*
*  call env for the converged model
*
      call env
      ifin = 0
      nel=nlim

      if( iprep .eq. 0 )then
*
* no reason to continue here
*
        return
      endif
      nel = nel - 1
      do 34 n=1,nel
	 index = nel-n+1
*
* these are placed in reverse order
*
         rrr = r(index)
         amm = m(index)
	 ttt = t(index)
         rho = rl1(index)
         ppp = p(index)

         aa(1,nl+n)=r(index)
	 x(1) = xcomp(index,1)
	 x(2) = xcomp(index,2)
	 x(3) = xcomp(index,3)
	 if ( n .eq. nel ) then
           x(1) = 0.
           x(2) = 0.
           x(3) = 1.
	 endif
*
* mass (grams)
*
         aa(2,nl+n)=10.**m(index)*amsun
*
* temperature (kelvins)
*
         aa(4,nl+n)=10.**t(index)
*
* density (g/cm^3)
*
         aa(5,nl+n)=10.**rl1(index)
*
* pressure (dynes/cm^2)
*
         aa(6,nl+n)=10.**p(index)
*
* energy generation rate + derivatives (eps, epsr, and epst)
*
         aa(7,nl+n)=0.0
         aa(11,nl+n)=0.0
         aa(12,nl+n)=0.0
*
* use delrad as temp. gradient except in convection zone, then use
* deltrue from gradt (or gradt2).
*
         if(rtg(index).lt.atg(index))then
            tdel=rtg(index)
         else
            tdel=ttg(index)
         endif
*
* temperature gradient (del)
*
         aa(15,nl+n)=tdel
*
* adiabatic temperature gradient (delad)
*
         aa(16,nl+n)=atg(index)
*
* helium profile
*
         aa(17,nl+n)=xcomp(index,2)
*
* radiative luminosity profile
*
         aa(3,nl+n)=alsun*(10.0**(blfin))*(ttg(index)/rtg(index))
*
* heat capacity Cv
*
         aa(8,nl+n)=ecv(index)
*
* pressure gradients chi rho and chi T
*
         aa(9,nl+n)=exr(index)
         aa(10,nl+n)=ext(index)
*
* opacity
*
         tl=t(index)
         pl=p(index)
         rhol = rl1(index)
	 oplo = opac3(rhol,tl)
         aa(18,nl+n)=10.**oplo
         ifinal = 1
*
* start 5-point opacity derivative stuff here.
*
         if(x(1).eq.1.0 .or. x(2).eq.1.0 .or. x(3).eq.1.0)then
           if(x(1).eq.1.0)then
             je=1
           endif
           if(x(2).eq.1.0)then
             je=2
           endif
           if(x(3).eq.1.0)then
             je=3
           endif
           call opacder(rhol,tl,opcs,dlklr,dlklt)
           aa(13,nl+n)=dlklr
           aa(14,nl+n)=dlklt
*
* note: replace previous opacity value if we can get it here.
*
           aa(18,nl+n)=10.**opcs
          else
*
* 5 point differencing
*
           oplo5=opac3(rhol+0.02,tl)
           oplo1=opac3(rhol+0.01,tl) 
           oplo2=opac3(rhol-0.01,tl) 
           oplo6=opac3(rhol-0.02,tl)
           oplo7=opac3(rho,tl+0.02)
           oplo3=opac3(rhol,tl+0.01)
           oplo4=opac3(rhol,tl-0.01)
           oplo8=opac3(rho,tl-0.02)
*
* opacity derivatives (kapr(13) and kapt(14))
*
         aa(13,nl+n)=(oplo6-(8.*oplo2)+(8.*oplo1)-oplo5)/(12.*0.01)
         aa(14,nl+n)=(oplo8-(8.*oplo4)+(8.*oplo3)-oplo7)/(12.*0.01)
         endif
*
* if He fraction is constant, set Ledoux to zero.
*
         if(aa(17,nl+n).eq.aa(17,nl+n-1))then
           aa(19,nl+n) = 0.0
           go to 34
         endif
         if(aa(17,nl+n-1).eq.0.0 .or. aa(17,nl+n).eq.0.0)then
           aa(19,nl+n) = 0.0
           go to 34
         endif
         if(index.eq.nel)then
           aa(19,nl+n) = 0.0
          else
*
* call up modified ledoux term subroutine
*
           call ledoux
           dly = ( dlog(aa(17,nl+n)) - dlog(aa(17,nl+n-1)) )
           dlp = ( dlog(aa(6,nl+n)) - dlog(aa(6,nl+n-1)) )
           bledoux=(1./aa(10,nl+n))*bled*(dly/dlp)
*
* Note: there is a minus sign difference between this ledoux term 
* in the core 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,nl+n)=bledoux
           if(x(1).ne.0.0 .and. x(2).gt.0.99 .and. x(3).ne.0.0)then
             aa(19,nl+n)= 0.d0
           endif
         endif
   34 continue
      ifinal = 0
*
* now write it out. 
*
      nshell=nl+nel

      if (tfit(ifit).lt.1.d-3) then
        write(50,1010) nshell, nl, nel
      endif

 1010 format(i5,'  =  ',i4,'  +  ',i4)
*
* smooth out del and delad at core-envelope boundary
* set mver for appropriate variable
* 1. delad; 2. del; 3. chiT; 4. chiRho
*
      do 33 i=1,nshell
        if(i.eq.nl)then
*
* delad
*
          diff = aa(16,i+1) - aa(16,i)
          if(diff.gt.0.002)then
            mver = 1
            call smooth(nl,mver)
          endif
*
* del 
*
          diff2 = aa(15,i+1) - aa(15,i)
          if(diff2.gt.0.002)then
            mver = 2
            call smooth(nl,mver)
          endif
*
* chi T
*
          diff3 = aa(10,i+1) - aa(10,i)
          if(diff3.gt.0.002)then
            mver = 3
            call smooth(nl,mver)
          endif
*
* chi rho
*
          diff2 = aa(9,i+1) - aa(9,i)
          if(diff4.gt.0.002)then
            mver = 4
            call smooth(nl,mver)
          endif
        endif
33    continue
*
* Write out model to tape50
*
      if (tfit(ifit).lt.1.d-3) then
        do 2 i=1,20
          write(50,1020) (aa(i,n), n=1,nshell)
2       continue
      endif

1020  format(1p,4e22.15)

      if (tfit(ifit).lt.1.d-3) then
        close(50)
      endif

      return
      end 

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

