program cno 
      implicit real*8(a-h,o-z)

c     ******************************************************
c      CNO Quadrup-cycle:
c        (with options for hydrostatic or explosive burning)
c     isotopes included are: p  4^he  12^c  13^c  14^n
c                            15^n  16^o  17^o  18^o  19^f
c                            13^n  15^o  17^f  18^f 
c
c     not really set up for beta-limited cno cycle, o14 not in network
c     *******************************************************


      parameter (nions=14)


c.... allocate space for nions isotopes

      double precision ai(nions),zi(nions),y(nions),dy(nions),x(nions)


c.... carry abundance arrays in common
c.... equivalence an array for indexing as a vector

      double precision rc12pg,r13nenu13c,rc13pg,rn14pginv,rn14pg,
     1                 r15oenu15n,rn15pa,rc12painv,rn15pg,ro16pginv,
     2                 ro16pg,r17fenu17o,ro17pa,rn14painv,ro17pg,
     3                 r18fenu18o,ro18pa,rn15painv,ro18pg,rf19pginv,
     4                 rf19pa,ro16painv,rn13pginv,ro15pginv,rf17pginv,
     5                 rf18pginv

      common  /rat/    rc12pg,r13nenu13c,rc13pg,rn14pginv,rn14pg,
     1                 r15oenu15n,rn15pa,rc12painv,rn15pg,ro16pginv,
     2                 ro16pg,r17fenu17o,ro17pa,rn14painv,ro17pg,
     3                 r18fenu18o,ro18pa,rn15painv,ro18pg,rf19pginv,
     4                 rf19pa,ro16painv,rn13pginv,ro15pginv,rf17pginv,
     5                 rf18pginv
     
      double precision ye,sumx,tot,enuc,snuc,time,tzero,rzero,
     1                 ndt,dth,t9,t,rho,delchi,chimin,fdtn,psi,ad

      double precision yp,y4he,y12c,y13c,y14n,y15n,y16o,y17o,y18o,y19f,
     1                 y13n,y15o,y17f,y18f
      common   /yy/    yp,y4he,y12c,y13c,y14n,y15n,y16o,y17o,y18o,y19f,
     1                 y13n,y15o,y17f,y18f
      equivalence (y(1),yp)

c.... initialize atomic weight number and nuclear charge arrays

      data  ai /1.d0,4.d0,12.d0,13.d0,14.d0,15.d0,16.d0,17.d0,18.d0,
     1          19.d0,13.d0,15.d0,17.d0,18.d0/
      data  zi /1.d0,2.d0,6.d0,6.d0,7.d0,7.d0,8.d0,8.d0,8.d0,9.d0,
     1          7.d0,8.d0,9.d0,9.d0/
   
c
c     read number of time steps, initial time step, temperature, 
c       density, and composition
c
      read (4,*)ndt
      write(6,*)ndt

      read (4,*)dth,t,rho
      write(6,*)dth,t,rho
       
      tzero=t/1.d9
      rzero=rho
      t9=tzero

c.... set those initial abundances not desired to be zero

c     input: mass fractions
      read (4,*)(x(ijk),ijk=1,nions)
      write(6,*)(x(ijk),ijk=1,nions)
      write(7,10) (zi(ijk),ijk=1,nions)
      write(7,10) (ai(ijk),ijk=1,nions)
   10 format (10x,14f10.1)

c     y = x / a
c     ye=sum(zi*yi)

      ye=0.d0
      do ijk = 1,nions
        y(ijk) = x(ijk)/ai(ijk)
        if(y(ijk).le.0.0)y(ijk) = 0.0d0
        ye=ye+zi(ijk)*y(ijk)
      end do

c     delchi=scaling for time step (usually 0.1, 0.05 is better)
c     chimin=minimum y to be considered in finding new time step (1.e-7?)
c     fdtn=maximum enlargement factor for new time step (usually 2.)
c
      read (4,*) delchi,chimin,fdtn
      write(6,*) delchi,chimin,fdtn

c     select hydrostatic (iburn=0), or explosive (iburn=1) burning
      read (4,*) iburn
      write(6,*) iburn

      if(iburn.eq.0) then
        t9=tzero
        rho=rzero
      endif

c
c     evolutionary loop
c

      time=0.d0
      do 180  m=1,ndt
      time = time + dth

c     ***************************************
      call burn (t9,rho,ye,dth,dy,enuc,snuc)
c     ***************************************


c     update abundances

      ye=0.d0
      do 80  i=1,nions
      y(i) = y(i) + dy(i)
  80  ye=ye+zi(i)*y(i)

c
c     update t9 and rho
c
       psi=1.0d0
       ad =3.0d0

c     ****************************************************
       call update (time,rzero,tzero,psi,ad,iburn,rho,t9)
c     ****************************************************

c
c     determine new time step
c
c     *****************************************
      call dtnuc (dth,chimin,delchi,fdtn,y,dy)
c     *****************************************

c
c     write abundance and energy generation history
c

      sumx = 0.d0
      do 170 i=1,nions
        x(i) = y(i)*ai(i)
  170 sumx = sumx + x(i) 
      tot = 1.d0 - sumx

c.. writing things out

      write (9,65) m, time, dth, t9, rho, enuc, snuc
   65 format (1x,i6,1pd12.4,5d12.2)
      write (7,70) time,(x(i), i=1,nions),tot
   70 format (1p16d10.2)
c.....reproducing clayton's fig. 5-15
c     write (8,75) time,(x(i)/(sumx-x(1)-x(2)), i=3,5),x(6)*1.d4/
c    1          (sumx-x(1)-x(2)),(x(i)/(sumx-x(1)-x(2)), i=7,nions),tot
c  75 format (1p14d10.2)

  180 continue


      stop
      end



cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

      subroutine dtnuc (dth,chimin,delchi,fdtn,y,dy)

      implicit real*8 (a-h,o-z)
      parameter (nions = 14)

      double precision y(nions),dy(nions),tau,fak,taug

      tau = fdtn*dth

      do 50  n=1,nions
      if (y(n).lt.(0.1d0*chimin)) go to 50
      fak = dy(n)
      if (fak)  20,50,20
   20 taug = dabs((y(n)+chimin)/fak)*delchi*dth
      if (taug) 50,50,30
   30 if (tau-taug) 50,50,40
   40 tau=taug
   50 continue
      dth=tau
      


      return
      end

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

      subroutine burn(t9,rho,ye,dth,dy,enuc,snuc)

      implicit real*8 (a-h,o-z)

c
c     this subroutine computes nucleosynthesis and energy generation

c.... required input

c        t9 = temperature in billions of degrees at new time point
c        rho =density in gm/cm**3 at new time point
c        dth = time step
c            as input variable h tells program over what interval
c               to evolve abundances.
c        y = the vector of old abundances.this will be returned unchanged.

c.... output quantities

c        dy = the vector of abundance changes that occurs during
c            the given time step dth.
c        enuc = nuclear energy generation in erg/gm/sec
c        snuc = net energy generation rate corrected for neutrino
c               losses (erg/g/sec) and weak interactions.
c
c        y is defined as mass fraction divided by nucleon number(moles/g).
c

      parameter (nions=14)
      double precision enuc,snuc,edot,sum,
     1                 n13neutl,o15neutl,f17neutl,f18neutl,
     2                 n13neu,o15neu,f17neu,f18neu
     
      common /abun/    aq(nions*nions),bq(nions)

c....matrices
      double precision a(nions,nions),b(nions)
      equivalence (aq,a)
      equivalence (bq,b)

c     reaction rate array

      double precision rc12pg,r13nenu13c,rc13pg,rn14pginv,rn14pg,
     1                 r15oenu15n,rn15pa,rc12painv,rn15pg,ro16pginv,
     2                 ro16pg,r17fenu17o,ro17pa,rn14painv,ro17pg,
     3                 r18fenu18o,ro18pa,rn15painv,ro18pg,rf19pginv,
     4                 rf19pa,ro16painv,rn13pginv,ro15pginv,rf17pginv,
     5                 rf18pginv

      common /rat/     rc12pg,r13nenu13c,rc13pg,rn14pginv,rn14pg,
     1                 r15oenu15n,rn15pa,rc12painv,rn15pg,ro16pginv,
     2                 ro16pg,r17fenu17o,ro17pa,rn14painv,ro17pg,
     3                 r18fenu18o,ro18pa,rn15painv,ro18pg,rf19pginv,
     4                 rf19pa,ro16painv,rn13pginv,ro15pginv,rf17pginv,
     5                 rf18pginv

      double precision yp,y4he,y12c,y13c,y14n,y15n,y16o,y17o,y18o,y19f,
     1                 y13n,y15o,y17f,y18f
      common /yy/      yp,y4he,y12c,y13c,y14n,y15n,y16o,y17o,y18o,y19f,
     1                 y13n,y15o,y17f,y18f

      double precision y(nions),dy(nions)
      equivalence (y(1),yp)

c     binding energy array for calculating nuclear energy generation rate
      double precision ai(nions),zi(nions),be(nions)

      data  ai   /1.d0,4.d0,12.d0,13.d0,14.d0,15.d0,16.d0,17.d0,
     1            18.d0,19.d0,13.d0,15.d0,17.d0,18.d0/
      data  zi   /1.d0,2.d0,6.d0,6.d0,7.d0,7.d0,8.d0,8.d0,8.d0,9.d0,
     1            7.d0,8.d0,9.d0,9.d0/
      data  be   /0.d0,28.2928d0,92.1624d0,97.1088d0,104.66d0,
     1            115.493d0,127.62d0,131.764d0,139.808d0,147.802d0,
     2            94.1064d0,111.956d0,128.221d0,137.371d0/


c
c     zeroing matrix
c

      do 10 izero=1,nions*nions
   10 aq(izero)=0.d0
      do 20 jzero=1,nions
   20 bq(jzero)=0.d0

c
c     subroutine rate generates thermonuclear reaction rates
c           ye is ne/(rho*na)
c
c    **********************
      call rate(t9,rho,ye)
c    **********************

c.... calculate the terms of the coupled, linearized rate equations
c.... for the abundances (a*dy=b). let the unknowns be the changes
c.... in y that occur in the time step h

      a(1,1) = 1.0d0/dth + y(3)*rc12pg + y(4)*rc13pg + y(5)*rn14pg 
     1       + y(6)*rn15pa + y(6)*rn15pg + y(7)*ro16pg + y(8)*ro17pa
     2       + y(8)*ro17pg + y(9)*ro18pa + y(9)*ro18pg + y(10)*rf19pa
      a(1,2) = -y(3)*rc12painv - y(5)*rn14painv - y(6)*rn15painv 
     1         - y(7)*ro16painv
      a(1,3) = y(1)*rc12pg - y(2)*rc12painv
      a(1,4) = y(1)*rc13pg
      a(1,5) = -rn14pginv + y(1)*rn14pg - y(2)*rn14painv
      a(1,6) = y(1)*rn15pa + y(1)*rn15pg - y(2)*rn15painv
      a(1,7) = -ro16pginv + y(1)*ro16pg - y(2)*ro16painv
      a(1,8) = y(1)*ro17pa + y(1)*ro17pg
      a(1,9) = y(1)*ro18pa + y(1)*ro18pg
      a(1,10)=- rf19pginv + y(1)*rf19pa 
      a(1,11)=- rn13pginv
      a(1,12)=- ro15pginv
      a(1,13)=- rf17pginv
      a(1,14)=- rf18pginv
      b(1)   =- y(1)*y(3)*rc12pg - y(1)*y(4)*rc13pg - y(1)*y(5)*rn14pg
     1        - y(1)*y(6)*rn15pa - y(1)*y(6)*rn15pg - y(1)*y(7)*ro16pg 
     2        - y(1)*y(8)*ro17pa - y(1)*y(8)*ro17pg - y(1)*y(9)*ro18pa 
     3        - y(1)*y(9)*ro18pg - y(1)*y(10)*rf19pa
     4        + y(5)*rn14pginv + y(2)*y(3)*rc12painv + y(7)*ro16pginv
     5        + y(2)*y(6)*rn15painv + y(2)*y(5)*rn14painv 
     6        + y(2)*y(7)*ro16painv + y(10)*rf19pginv
c.... reactions involving isotopes which beta decay
     7        + y(11)*rn13pginv + y(12)*ro15pginv + y(13)*rf17pginv
     8        + y(14)*rf18pginv

      a(2,1) =- y(6)*rn15pa - y(8)*ro17pa - y(9)*ro18pa - y(10)*rf19pa
      a(2,2) = 1.0d0/dth + y(3)*rc12painv + y(5)*rn14painv 
     1       + y(6)*rn15painv + y(7)*ro16painv
      a(2,3) = y(2)*rc12painv
      a(2,4) = 0.d0
      a(2,5) = y(2)*rn14painv
      a(2,6) =- y(1)*rn15pa + y(2)*rn15painv
      a(2,7) = y(2)*ro16painv
      a(2,8) =- y(1)*ro17pa
      a(2,9) =- y(1)*ro18pa
      a(2,10)=- y(1)*rf19pa
      a(2,11)= 0.0d0
      a(2,12)= 0.0d0
      a(2,13)= 0.0d0
      a(2,14)= 0.0d0
      b(2)   = y(1)*y(6)*rn15pa + y(1)*y(8)*ro17pa + y(1)*y(9)*ro18pa
     1       + y(1)*y(10)*rf19pa
     2       - y(2)*y(3)*rc12painv - y(2)*y(5)*rn14painv 
     3       - y(2)*y(6)*rn15painv - y(2)*y(7)*ro16painv

      a(3,1) = y(3)*rc12pg - y(6)*rn15pa 
      a(3,2) = y(3)*rc12painv
      a(3,3) = 1.d0/dth + y(1)*rc12pg + y(2)*rc12painv
      a(3,4) = 0.0d0
      a(3,5) = 0.0d0
      a(3,6) = - y(1)*rn15pa
      a(3,7) = 0.0d0
      a(3,8) = 0.0d0
      a(3,9) = 0.0d0
      a(3,10)= 0.0d0
      a(3,11)=- rn13pginv
      a(3,12)= 0.0d0
      a(3,13)= 0.0d0
      a(3,14)= 0.0d0
      b(3)   =- y(1)*y(3)*rc12pg + y(1)*y(6)*rn15pa 
     1        - y(2)*y(3)*rc12painv
c.... reactions involving isotopes which beta decay
     2        + y(11)*rn13pginv
 
      a(4,1) = y(4)*rc13pg 
      a(4,2) = 0.0d0
      a(4,3) = 0.0d0
      a(4,4) = 1.0d0/dth + y(1)*rc13pg
      a(4,5) = -rn14pginv
      a(4,6) = 0.0d0
      a(4,7) = 0.0d0
      a(4,8) = 0.0d0
      a(4,9) = 0.0d0
      a(4,10)= 0.0d0
      a(4,11)=- r13nenu13c
      a(4,12)= 0.0d0
      a(4,13)= 0.0d0
      a(4,14)= 0.0d0
      b(4)   =- y(1)*y(4)*rc13pg 
     1        + y(5)*rn14pginv
c.... reactions involving isotopes which beta decay
     2        + y(11)*r13nenu13c

      a(5,1) = y(5)*rn14pg - y(4)*rc13pg - y(8)*ro17pa
      a(5,2) = y(5)*rn14painv
      a(5,3) = 0.0d0
      a(5,4) =- y(1)*rc13pg
      a(5,5) = 1.0d0/dth + y(1)*rn14pg + rn14pginv + y(2)*rn14painv
      a(5,6) = 0.0d0
      a(5,7) = 0.0d0
      a(5,8) =- y(1)*ro17pa
      a(5,9) = 0.0d0
      a(5,10)= 0.0d0
      a(5,11)= 0.0d0
      a(5,12)=- ro15pginv
      a(5,13)= 0.0d0
      a(5,14)= 0.0d0
      b(5)   =- y(1)*y(5)*rn14pg + y(1)*y(4)*rc13pg + y(1)*y(8)*ro17pa
     1        - y(2)*y(5)*rn14painv - y(5)*rn14pginv
c.... reactions involving isotopes which beta decay
     2        + y(12)*ro15pginv

      a(6,1) =  y(6)*rn15pa + y(6)*rn15pg - y(9)*ro18pa
      a(6,2) = -y(3)*rc12painv + y(6)*rn15painv
      a(6,3) = -y(2)*rc12painv
      a(6,4) = 0.0d0
      a(6,5) = 0.0d0
      a(6,6) = 1.0d0/dth + y(1)*rn15pa + y(1)*rn15pg + y(2)*rn15painv
      a(6,7) = - ro16pginv
      a(6,8) = 0.0d0
      a(6,9) = -y(1)*ro18pa
      a(6,10)= 0.0d0
      a(6,11)= 0.0d0
      a(6,12)=- r15oenu15n
      a(6,13)= 0.0d0
      a(6,14)= 0.0d0
      b(6)   = - y(1)*y(6)*rn15pa - y(1)*y(6)*rn15pg
     1       + y(1)*y(9)*ro18pa
     2       + y(2)*y(3)*rc12painv - y(2)*y(6)*rn15painv 
     3       + y(7)*ro16pginv
c.... reactions involving isotopes which beta decay
     4       + y(12)*r15oenu15n

      a(7,1) =- y(6)*rn15pg + y(7)*ro16pg - y(10)*rf19pa
      a(7,2) = y(7)*ro16painv
      a(7,3) = 0.0d0
      a(7,4) = 0.0d0
      a(7,5) = 0.0d0
      a(7,6) = -y(1)*rn15pg
      a(7,7) = 1.0d0/dth + y(1)*ro16pg + ro16pginv + y(2)*ro16painv
      a(7,8) = 0.0d0
      a(7,9) = 0.0d0
      a(7,10)=- y(1)*rf19pa
      a(7,11)= 0.0d0
      a(7,12)= 0.0d0
      a(7,13)=- rf17pginv
      a(7,14)= 0.0d0
      b(7)   = y(1)*y(6)*rn15pg - y(1)*y(7)*ro16pg + y(1)*y(10)*rf19pa
     1       - y(2)*y(7)*ro16painv - y(7)*ro16pginv
c.... reactions involving isotopes which beta decay
     2       + y(13)*rf17pginv

      a(8,1) = y(8)*ro17pa + y(8)*ro17pg
      a(8,2) =- y(5)*rn14painv
      a(8,3) = 0.0d0
      a(8,4) = 0.0d0
      a(8,5) =- y(2)*rn14painv 
      a(8,6) = 0.0d0
      a(8,7) = 0.0d0
      a(8,8) = 1.0d0/dth + y(1)*ro17pa + y(1)*ro17pg
      a(8,9) = 0.0d0
      a(8,10)= 0.0d0
      a(8,11)= 0.0d0
      a(8,12)= 0.0d0
      a(8,13)=- r17fenu17o
      a(8,14)=- rf18pginv
      b(8)   = - y(1)*y(8)*ro17pa - y(1)*y(8)*ro17pg
     1       + y(2)*y(5)*rn14painv
c.... reactions involving isotopes which beta decay
     2       + y(13)*r17fenu17o + y(14)*rf18pginv

      a(9,1) = y(9)*ro18pa + y(9)*ro18pg
      a(9,2) =- y(6)*rn15painv
      a(9,3) = 0.0d0
      a(9,4) = 0.0d0
      a(9,5) = 0.0d0
      a(9,6) =- y(2)*rn15painv
      a(9,7) = 0.0d0
      a(9,8) = 0.0d0
      a(9,9) = 1.0d0/dth + y(1)*ro18pa + y(1)*ro18pg
      a(9,10)= - rf19pginv
      a(9,11) = 0.0d0
      a(9,12) = 0.0d0
      a(9,13) = 0.0d0
      a(9,14) =- r18fenu18o
      b(9)   = - y(1)*y(9)*ro18pa - y(1)*y(9)*ro18pg
     1       + y(2)*y(6)*rn15painv + y(10)*rf19pginv
c.... reactions involving isotopes which beta decay
     2       + y(14)*r18fenu18o 

      a(10,1) =- y(9)*ro18pg + y(10)*rf19pa
      a(10,2) =- y(7)*ro16painv
      a(10,3) = 0.0d0
      a(10,4) = 0.0d0
      a(10,5) = 0.0d0
      a(10,6) = 0.0d0
      a(10,7) =- y(2)*ro16painv 
      a(10,8) = 0.0d0
      a(10,9) =- y(1)*ro18pg
      a(10,10)= 1.0d0/dth + y(1)*rf19pa + rf19pginv
      a(10,11)= 0.0d0
      a(10,12)= 0.0d0
      a(10,13)= 0.0d0
      a(10,14)= 0.0d0
      b(10)   = y(1)*y(9)*ro18pg - y(1)*y(10)*rf19pa 
     1        + y(2)*y(7)*ro16painv - y(10)*rf19pginv

      a(11,1)= - y(3)*rc12pg
      a(11,2)= 0.0d0
      a(11,3)= - y(1)*rc12pg
      a(11,4)= 0.0d0
      a(11,5)= 0.0d0
      a(11,6)= 0.0d0
      a(11,7)= 0.0d0
      a(11,8)= 0.0d0
      a(11,9)= 0.0d0
      a(11,10)= 0.0d0
      a(11,11)= 1.0d0/dth + r13nenu13c + rn13pginv
      a(11,12)= 0.0d0
      a(11,13)= 0.0d0
      a(11,14)= 0.0d0
      b(11)   = y(1)*y(3)*rc12pg - y(11)*rn13pginv
     1        - y(11)*r13nenu13c 
    
      a(12,1)=- y(5)*rn14pg
      a(12,2)= 0.0d0
      a(12,3)= 0.0d0
      a(12,4)= 0.0d0
      a(12,5)=- y(1)*rn14pg
      a(12,6)= 0.0d0
      a(12,7)= 0.0d0
      a(12,8)= 0.0d0
      a(12,9)= 0.0d0
      a(12,10)= 0.0d0
      a(12,11)= 0.0d0
      a(12,12)= 1.0d0/dth + r15oenu15n + ro15pginv
      a(12,13)= 0.0d0
      a(12,14)= 0.0d0
      b(12)   = y(1)*y(5)*rn14pg - y(12)*ro15pginv 
     1        - y(12)*r15oenu15n

      a(13,1)=- y(7)*ro16pg 
      a(13,2)= 0.0d0
      a(13,3)= 0.0d0
      a(13,4)= 0.0d0
      a(13,5)= 0.0d0
      a(13,6)= 0.0d0
      a(13,7)=- y(1)*ro16pg 
      a(13,8)= 0.0d0
      a(13,9)= 0.0d0
      a(13,10)= 0.0d0
      a(13,11)= 0.0d0
      a(13,12)= 0.0d0
      a(13,13)= 1.0d0/dth + r17fenu17o + rf17pginv
      a(13,14)= 0.0d0
      b(13)   = y(1)*y(7)*ro16pg - y(13)*rf17pginv
     1        - y(13)*r17fenu17o

      a(14,1)=- y(8)*ro17pg 
      a(14,2)= 0.0d0
      a(14,3)= 0.0d0
      a(14,4)= 0.0d0
      a(14,5)= 0.0d0
      a(14,6)= 0.0d0
      a(14,7)= 0.0d0
      a(14,8)=- y(1)*ro17pg 
      a(14,9)= 0.0d0
      a(14,10)= 0.0d0
      a(14,11)= 0.0d0
      a(14,12)= 0.0d0
      a(14,13)= 0.0d0
      a(14,14)= 1.0d0/dth + r18fenu18o + rf18pginv
      b(14)   = y(1)*y(8)*ro17pg - y(14)*rf18pginv
     1        - y(14)*r18fenu18o

c
c     now invert the matrix to obtain the vector of new abundances b.

      nn=nions
      
c     *************
      call leqs(nn)
c     *************

c.... compute energy generation 
c.... b(i) is now delta y's :

50020 do 50024 i1=1,nions
50024 dy(i1) = b(i1)

c     energy generation here ignores neutrino losses

      edot=0.d0
      do 30  i1=1,nions
      edot = edot + dy(i1)*be(i1)
   30 continue
      enuc = edot*1.602d-6*6.02259d23/dth

c.
c..compute the average energy of the neutrinos from the expression
c..on page 142 of the ay289j notes
c.
c..for 13n(e+nu)13c
      sum = be(4) - be(11) -0.7826 - 1.022
      sum = 1.0 + sum/0.511
      n13neutl = 0.5 * sum * 0.511 * (1.0 - 1.0/(sum*sum)) 
     1          * (1.0 - 1.0/(4.0*sum) - 1.0/(9.0*sum*sum))
c..
c..for 15o(e+nu)15n
      sum = be(6) - be(12) -0.7826- 1.022
      sum = 1.0 + sum/0.511
      o15neutl = 0.5 * sum * 0.511 * (1.0 - 1.0/(sum*sum)) 
     1          * (1.0 - 1.0/(4.0*sum) - 1.0/(9.0*sum*sum))
c..
c..for 17f(e+nu)17o
      sum = be(8) - be(13) -0.7826 - 1.022
      sum = 1.0 + sum/0.511
      f17neutl = 0.5 * sum * 0.511 * (1.0 - 1.0/(sum*sum)) 
     1          * (1.0 - 1.0/(4.0*sum) - 1.0/(9.0*sum*sum))
c..
c..for 18f(e+nu)18o
      sum = be(9) - be(14) -0.7826 - 1.022
      sum = 1.0 + sum/0.511
      f18neutl = 0.5 * sum * 0.511 * (1.0 - 1.0/(sum*sum)) 
     1          * (1.0 - 1.0/(4.0*sum) - 1.0/(9.0*sum*sum))
c.

c....include neutrino losses 9.648e+17 is 6.02e+23*1.602e-6
c....compute neutrino energy loss rates
c.
       n13neu  = y(11)*r13nenu13c*9.647d+17*n13neutl
       o15neu  = y(12)*r15oenu15n*9.647d+17*o15neutl
       f17neu  = y(13)*r17fenu17o*9.647d+17*f17neutl
       f18neu  = y(14)*r18fenu18o*9.647d+17*f18neutl
c.

      snuc = enuc - n13neu - o15neu - f17neu - f18neu
      return
      end

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


      subroutine rate(tt9,rrho,ye)

      implicit real*8 (a-h,o-z)

c
c     this subroutine generates nuclear reaction rate factors lambda
c     using analytical fitting functions to experimental (cf88 & fczii)

      double precision rc12pg,r13nenu13c,rc13pg,rn14pginv,rn14pg,
     1                 r15oenu15n,rn15pa,rc12painv,rn15pg,ro16pginv,
     2                 ro16pg,r17fenu17o,ro17pa,rn14painv,ro17pg,
     3                 r18fenu18o,ro18pa,rn15painv,ro18pg,rf19pginv,
     4                 rf19pa,ro16painv,rn13pginv,ro15pginv,rf17pginv,
     5                 rf18pginv

      common /rat/     rc12pg,r13nenu13c,rc13pg,rn14pginv,rn14pg,
     1                 r15oenu15n,rn15pa,rc12painv,rn15pg,ro16pginv,
     2                 ro16pg,r17fenu17o,ro17pa,rn14painv,ro17pg,
     3                 r18fenu18o,ro18pa,rn15painv,ro18pg,rf19pginv,
     4                 rf19pa,ro16painv,rn13pginv,ro15pginv,rf17pginv,
     5                 rf18pginv
c.
c......calculate rate factors
c.

c..declare
      double precision t9,rho,
     1                 a,b,c,d,e,facta,factb,factd,facte,factf,
     2                 t9a13,t9a56,t912,t932,t913,t923,t943,t953,t965,
     3                 factor,theta
      parameter        (factor=0.1, theta=0.1)
c..
c..define some functions
      facta(a,b,c)=a/t923*exp(-b/t913-(t9/c)**2)
      factb(a,b,c,d,e)=1.d0+a*t913+b*t923+c*t9+d*t943+e*t953
      factd(a,b)=a*dexp(-b/t9)
      facte(a,b)=a*t932*exp(-b/t9)
      factf(a,b)=a/t932*exp(-b/t9)
      t9a13(a) = (t9/(1.d0+a*t9))**(1./3.)
      t9a56(a) = (t9/(1.d0+a*t9))**(5./6.)
c..
c..initialize
      t9=tt9
      rho=rrho
      t912=t9**(0.5)
      t932=t912*t9
      t913=t9**(1./3.)
      t923=t913**2
      t943=t923**2
      t953=t913**5 
      t965=t9**(6.0/5.0)
c..
c..10 cno cycle I reactions:
c.
c..12c(p,g)13n
c.
      rc12pg = rho * (facta(2.04d+07,1.369d+01,1.5d+00)
     1            * factb(3.d-02,1.19d0,2.54d-01,2.06d0,1.12d0)
     2            + factf(1.08d+05,4.925d0)
     3            + factf(2.15d+05,1.8179d+01))
c..inverse of 'rc12pg' is 'rn13pginv'
      rn13pginv = rc12pg/rho * facte(8.84d+09,2.2554d+01)
c.
c..13n(e+nu)13c  mean lifetime is 870 sec  (from clayton)
c.
      r13nenu13c = 1.0/870.0
c.
c..13c(p,g)14n  cf88 updated by caughlan & fowler tnrr5 1988 q = 7.551
c.
      rc13pg = rho * (8.01e+07/t923 * exp(-13.717/t913-(t9/2.0)**2)
     1                * (1.+0.030*t913+0.958*t923
     2                   + 0.204*t9+1.39*t943+0.753*t953)
     3                  + 1.21e+06/t965 * exp(-5.701/t9))
c..inverse of 'rc13pg' is 'rn14pginv'
      rn14pginv = rc13pg/rho * 1.19e+10*t932*exp(-87.621/t9)

c.
c..14n(p,g)15o  cf88 updated by caughlan & fowler tnrr5 1988 q= 7.297
c.
      rn14pg = rho * (4.90e+07/t923*exp(-15.228/t913-(t9/3.294)**2)
     1                * (1.+0.027*t913-0.778*t923
     2                   -0.149*t9+0.261*t943+0.127*t953)
     3                + 2.37e+03/t932*exp(-3.011/t9)
     4                + 2.19e+04*exp(-12.530/t9))
c..inverse of 'rn14pg' is 'ro15pginv'
      ro15pginv = rn14pg/rho * 2.70e+10*t932*exp(-84.678/t9)
c.
c..15o(e+nu)15n  mean lifetime is 178 sec  (from clayton)
c.
      r15oenu15n = 1.0/178.0
c.
c..15n(p,a)12c fcz3
c.
      rn15pa = rho * (1.08d+12/t923*exp(-15.251/t913-(t9/0.522)**2)
     1                * (1.+0.027*t913+2.62*t923+0.501*t9+5.36*t943
     2                   + 2.60*t953)
     3                + 1.19d+08/t932*exp(-3.676/t9) 
     4                + 5.41d+08/t912 * exp(-8.926/t9)
     5                + theta*(4.72d+08/t932*exp(-7.721/t9) +
     6                         2.20d+09/t932*exp(-11.418/t9)))
c..inverse of 'rn15pa' is 'rc12painv'
      rc12painv = rn15pa*7.06d-01*exp(-57.625/t9)
c..
c.. complete cycle I
c..
c..7 cno cycle II reactions:
c.
c..15n(p,g)16o
c.
      rn15pg = rho * ((facta(9.78d+08,1.5251d+01,4.5d-01)
     1                * factb(2.7d-02,2.19d-01,4.2d-02,6.83d0,3.32d0)
     2                + factf(1.114d+04,3.328d0)
     3                + factf(1.49d+04,4.665d0)
     3                + factf(3.80d+06,1.1048d+01)))
c..inverse of 'rn15pg' is 'ro16pginv'
      ro16pginv = rn15pg/rho * facte(3.62d+10,1.40741d+02)
c.
c..16o(p,g)17f
c.
      ro16pg = rho*(1.50d+08/(t923*(1.+2.13*(1.-exp(-0.728
     1            *t923))))*exp(-16.692/t913))
c..inverse of 'ro16pg' is 'rf17pginv'
      rf17pginv = ro16pg*facte(3.03d+09,6.970d0)/rho
c.
c..17f(e+nu)17o  mean lifetime is 95 sec (from clayton)
c.
      r17fenu17o = 1.0/95.0
c.
c..17o(p,a)14n
c.
      ro17pa = rho * (facta(1.53d+07,1.6712d+01,5.65d-01)
     1                * factb(2.5d-02,5.39d0,9.4d-01,1.35d+01,5.98d0)
     2                + factd(2.92d+06,4.247d0)*t9 + factor
     3                * facta(4.81d+10,1.6712d+01,4.d-02)*t953 + factor
     4                * factf(5.05d-05,7.23d-01)
     4                + factor*factf(1.31d+01,1.961d0))
c..inverse of 'ro17pa' is 'rn14painv'
      rn14painv = ro17pa * factd(6.76d-01,1.3845d+01)
c..
c.. complete cycle II
c..
c..5 cno cycle III reactions:
c.
c..17o(p,g)18f
c.
      ro17pg = rho * (7.97e+07*t9a56(2.69d0)/t932*exp(-16.712/
     1                  t9a13(2.69d0))
     2                + 1.51e+08/t923*exp(-16.712/t913)
     3                * factb(2.5d-02,-5.1d-02,-8.82d-3,0.d0,0.d0)
     4                + factd(1.56d+05,6.272d0)/t9
     5                + factor * factf(1.31d+01,1.961d0))
c..inverse of 'ro17pg' is 'rf18pginv'
      rf18pginv = ro17pg * facte(3.66d+10,6.5092d+01)/rho
c.
c..18f(e+nu)18o  mean lifetime is 9504 sec (from nuclear wallet card)
c.
      r18fenu18o = 1.0/9504.0
c.
c..18o(p,a)15n  cf88 updated by caughlan & fowler tnrr5 1988  q = 3.980
c.
      ro18pa = rho * (3.63e+11/t923*exp(-16.729/t913-(t9/1.361)**2)
     1                * (1.+0.025*t913+1.88*t923+0.327*t9
     2                   +4.66*t943+2.06*t953)
     3                + 9.90e-14/t932*exp(-0.231/t9)
     4                + 2.66e+04/t932*exp(-1.670/t9)
     5                + 2.41e+09/t932*exp(-7.638/t9)
     6                + 1.46e+09/t9*exp(-8.310/t9))
c..inverse of 'ro18pa' is 'rn15painv'
      rn15painv = ro18pa *1.66e-01*exp(-46.191/t9)
c..
c.. complete cycle III
c..
c..4 cno cycle IV reactions:
c.
c..18o(p,g)19f  cf88 updated by caughlan & fowler tnrr5 1988 q = 7.994
c.
      ro18pg =  rho * (3.45e+08/t923*exp(-16.729/t913-(t9/0.139)**2)
     1                * (1.+0.025*t913+2.26*t923
     2                    + 0.394*t9+30.56*t943+13.55*t953)
     3                + 1.25e-15/t932*exp(-0.231/t9) 
     4                + 1.64e+02/t932*exp(-1.670/t9)
     5                + 1.28e+04*t912*exp(-5.098/t9))
c..inverse of 'ro18pg' is 'rf19pginv'
      rf19pginv = ro18pg/rho * 9.20e+09*t932*exp(-92.769/t9) 
c.
c..19f(p,a)16o fcz3
c.
      rf19pa = rho * (3.55d+11/t923*exp(-18.113/t913-(t9/0.845)**2)
     1                * (1.+0.023*t913+1.96*t923
     2                    + 0.316*t9+2.86*t943+1.17*t953)
     3                + 3.67d+06/t932*exp(-3.752/t9) 
     4                + 3.07d+08*exp(-6.019/t9)) /
     5                (1.+4.*exp(-2.090/t9)+7.*exp(-16.440/t9))
c..inverse of 'rf19pa' is 'ro16painv'
      ro16painv =  rf19pa * factd(6.54d-01,9.4159d+01) 


      return
      end

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


      subroutine leqs(n)

      implicit real*8 (a-h,o-z)
      parameter (nions=14)

c....
c.....leqs performs matrix inversion
c....

      double precision r,c

      common  /abun/   aq(nions*nions),bq(nions)
      double precision a(nions,nions),b(nions)
      equivalence (aq,a)
      equivalence (bq,b)

      n1 = n - 1

c     find maximum element in each row, and divide
      do 19 i=1,n
      r=abs(a(i,1))

      do 16 j=2,n
      c=abs(a(i,j))
      if (r.lt.c) r=c
   16 continue
      do 17 j=1,n
   17 a(i,j)=a(i,j)/r
      b(i)=b(i)/r
   19 continue
      do 9 j=1,n1
      l=j+1
c.....no pivoting
      do 9 i=l,n
      r=-a(i,j)
      if (r.eq.0.0) go to 9
      r=r/a(j,j)
      do 7 k=l,n
    7 a(i,k)=a(i,k)+r*a(j,k)
      b(i)=b(i)+r*b(j)
    9 continue
c....
c     the matrix is now in triangular form
c     start the back substitution and find solution
c....
      b(n)=b(n)/a(n,n)
      do 13 l=1,n1
      i=n-l
      r=0.0
      imax=i+1
      do 12 j=imax,n
      jj=i+n+1-j
   12 r=r+a(i,jj)*b(jj)
   13 b(i)=(b(i)-r)/a(i,i)
      return
      end


cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

      subroutine update (tiem,rzero,tzero,psi,ad,iburn,rho,t9)
c
      implicit real*8 (a-h,o-z)
      double precision tiem,rzero,tzero,psi,ad,rho,t9,vzero,v,taud,taut

      if(iburn.eq.1) then
       vzero=1.0/rzero
       taud=psi*446.*sqrt(vzero)
c....
       taut=ad*taud
       v=vzero*exp(tiem/taud)
       rho=1.0/v
       t9=tzero/exp(tiem/taut)
      else
c....
       rho=rzero
       t9 = tzero
      endif
c....
      return
      end