mcbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb
cbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb
cbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb
c
c  Begin of the file "h03.f"                        date: November 14, 1992

c ==========                                    ==========
c      <<<<<<    Stellar evolution package    >>>>>>
c ==========                                    ==========
c
c        This code has been originally developed by
c         Bohdan Paczynski  and  Maciej Kozlowski,
c          some parts of it have been written by
c                   Stephen J. Ratcliff
c
c            and this version has been prepared
c                   and distributed by:
c
c      Ryszard Sienkiewicz
c      Nicolaus Copernicus Astronomical Center
c      Bartycka 18, 00-716 Warsaw, POLAND
c
c      e-mail: rs@camk.edu.pl
c      FAX   : (intl) - 48 (Poland) - 22 (Warsaw) - 410046
c      telex : 81 3978  (zapan pl)
c
c ==========                                    ==========
c               ANY COMMENTS ARE MOST WELCOME
c ==========                                    ==========

c*************************************************************************
c   These programs and data come with absolutely no warranty of any kind.
c   The entire risk as to the quality and performance of this package is
c   with you.
c*************************************************************************
c
c   Attention:    read "readme.03" file, please.

c
c      This file contains "STANDARD" STELLAR EVOLUTION PROGRAM "h03".

c      This first version is able to follow evolution of hydrogen burning
c      (through p-p chain and CN cycle) models of single, non-rotating,
c      intermediate mass stars, up to the "helium flash" only. To follow
c      evolution of both low mass and massive stars the code has to be improved.
c      Future versions (in preparation) are expected to be more general purpose.
c      The main aim of this basic version is to provide an user with a simple,
c      well described code, parts of which can be easily changed/replaced
c      to the user's taste.

c      Some important limitations (features ?...):

c      - As concerns nuclear burning, the abundances of all elements - other
c        than hydrogen and helium-4 - are always assumed to have their local
c        equilibrium values, i.e. an evolution of these abundances due to
c        nuclear burning and/or convective mixing, is NOT followed explicitly.

c      - Only the CN cycle (instead of the CNO tri-cycle) is taken into account.
c        The "nitrogen-14(N14) + proton = oxygen-15 + gamma ray" reaction
c        represents the CN cycle, with the N14 content being the sum
c        of primordial carbon-12 and N14 abundances.

c      - This code does not mix matter within outer convective envelopes.
c        Improving this code, one has to remember as well, that:
c          - partial ionizations zones may, in some cases, extend very deep,
c            especially during a red giant phase,
c          - through large parts of the convective envelopes convection may be
c            extremely non-adiabatic.

c      - Any mass loss from stellar surface, as well as

c      - any mixing in semiconvective regions are ignored.

c      - An effect of convective overshooting is - only partly and indirectly -
c        taken into account by applying the Schwarzschild stability criterion
c        AT THE CONVECTIVE SIDE of a convective core boundary.
c        The convective core boundary is assumed to be at the mass point where
c        radiative temp. gradient /adiabatic temp. gradient is less than one
c        for the first time, going from a center.

c      - a term accounting for the gravitational energy release due to time
c        changes of a chemical composition is neglected.

c        The above-mentioned limitations have to be removed, if one wants
c        to follow evolution of some type of stars (e.g., the low-mass stars),
c        especially when precise determination of nuclear species abundances
c        is important and/or when the outer convective envelope penetrates
c        nuclear burning affected layers.
c        See also the files : "eostate.03", "nuclear.03", "opacity.03",
c        for other limitations concerning input physics.

c        When compiling, the program is automatically supplemented with
c        five extra disk files:
c         - "commonsb.03"
c         - "data.03"
c         - "opacity.03"
c         - "nuclear.03"
c         - "eostate.03"                (see following comments)
c
c ----------------------------------------------------------------------------
c Input data:
c    -) a disk file "mod####.h03" which contains:
c       -) a Zero-age Hydrogen Main Sequence static stellar model prepared
c          with a "Schwarzschild-type" program "s03" (in this case the file
c          name is "mod0000.h03"), or
c       -) a non-zero-age stellar model calculated with a "Henyey-type",
c          stellar evolution program "h03" (in this case a symbol "####"
c          stands for a sequential model number).
c       All these files have to be stored with the subroutine "writeh"
c       and are read with the subroutine "readh" (see the file "commonsb.03")
c
c    -) screen input (self-explanatory)
c
c    -) other input data files, as requested by included subroutines
c       (see "include" directives at the end of the code) i.e.:
c        -) opacity data files - "opint1",....,"opint7"
c        -) opacity data file  - "envopa" prepared with program "k02"
c            for given values of hydrogen and heavy elements content
c        -) equation-of-state data file - "electron.gas"
c
c Output data:
c    -) to a disk file "evol.h03" (one data line per model)
c    -) to a disk file "mod####.h03" (for every "nwrmod'th" model only)
c    -) screen output (see explanations inside the code).
c
c ----------------------------------------------------------------------------
c WARNING: rename, or remove existing "evol.h03" disk file
c          before running the program again !
c          Real*8 (double precision) has extended exponents.
c          On other machines there may be fatal underflows or zero divides.
c ----------------------------------------------------------------------------
c
c "h03.f"  is the main part of the code and includes the main program
c            and the six subroutines: "extrap", "profil", "hsolve",
c            "deriv", "rhsh", "addcor"  - all others subroutines are
c            located on separate disk files and are included through
c            "include" directives.
c
c   name:         subject:
c   -----         --------
c
c   main
c  program - SUBROUTINES CALLED:
c
c             - "readh"  - (included - file "commonsb.03")
c                          reads a model from the disk file "mod####.h03"
c             - "grid1"  - (included - file "commonsb.03")
c                          redistributes mass points
c             - "extrap" -
c             - "profil" -
c             - "hsolve" -
c             - "addcor" -
c             - "writeh" - (included - file "commonsb.03")
c                          writes a model on the disk file "mod####.h03"
c
c "extrap" - estimates new time step
c
c "profil" - estimates new hydrogen and helium profiles due to
c            nuclear burning ( + mixing, if convective core exists)
c            SUBROUTINES CALLED:
c             - "rhsh"
c             - "nburn"   (included - file "nuclear.03")
c
c "hsolve" - calculates time changes of all variables using a "Henyey-scheme"
c            SUBROUTINES CALLED:
c             - "env"     (included - file "commonsb.03")
c             - "deriv"
c             - "center"  (included - file "commonsb.03")
c
c "deriv"  - calculates the right hand sides of stellar structure equations
c            and their derivatives
c            SUBROUTINES CALLED:
c             - "rhsh"
c
c "rhsh"   - calculates the right hand sides of stellar structure equations
c            SUBROUTINES CALLED:
c             - "state"   (included - file "eostate.03")
c             - "opact"   (included - file "opacity.03")
c             - "nburn"   (included - file "nuclear.03")
c             - "mixlen"  (included - file "commonsb.03")
c "addcor" - adds the time changes to all variables.
c
c ----------------------------------------------------------------------------
c ----------------------------------------------------------------------------
c
c There are some separated files which provide the code with "standard"
c       physics and boundary conditions. They are listed below.
c T h e s e   f i l e s   m a y   b e   e a s i l y   r e p l a c e d, if ne-
c cessary. However, to avoid changes in the main part, new SUBROUTINES should:
c   -) provide the program with actually requested data
c   -) have the same names as the old ones
c   -) be connected with the code in the same manner (i.e., through proper
c      parameters and/or common blocks - see appropriate "call" statements).
c
c  disk file:             subject:
c  ----------             --------
c
c "data.03"     - model parameters and constants
c "eostate.03"  - equation of state includes contributions from:
c                  nondegenerate/degenerate electrons, nondegenerate ions
c                  and radiation, as well as from partially ionized gas
c                  (reads the disk file "electron.gas")
c "opacity.03"  - reads and interpolates coefficient of opacity from
c                  the Los Alamos radiative opacity tables with conductive
c                  opacities and Alexander's molecular opacities added
c                  (the disk files "opint1"..."opint7" are used in the
c                  internal part of a model and the disk file "envopa"
c                  is used in an envelope )
c "nuclear.03"  - calculates nuclear reaction rates for p-p chain and CN cycle
c                  + neutrino energy loss rate (if hydrogen content = 0)
c "commonsb.03" - includes subroutines common for the "s03" and "h03" codes.
c
c ----------------------------------------------------------------------------
c ----------------------------------------------------------------------------
c
c
c HEREAFTER:  "time difference" (or "time change") stands for the change
c             of a variable (the array "hx") during the current time step
c             (the array "dh").
c
c             "Gradient" alone stands for the temperature gradient.
c
c..............................................................................
c$DEBUG
      implicit double precision(a-h,o-z)
      include 'data.03'
      common/hen/ hx(m2,nt), dh(m2,nt), nh
      dimension outbur(20)
c                                     --- opening file for intermediate results
      iout2=2
      open(iout2,file='evol.h03',status='new')
c                                     --- Reading model. Asking about file name
c                                         Unit No 1 opened inside, then closed.
      call readh(iout2)
c                                     --- Screen input
  222 format(/,' SELECT the following parameters :',/,
     $       ' number of models to be calculated,         nofmod= ',$)
  223 format(' every nwrmod`th model data will be stored, nwrmod= ',$)
  224 format(/,' Thanks ! - end of screen input - now, wait...')
      write(*,222)
      read(*,*)nofmod
      write(*,223)
      read(*,*)nwrmod
      write(*,224)
c                                      --- end of input
c
c ..................................   --- begin of  N E W   M O D E L
      do 3 istep=1,nofmod
      model=model+1
      nh1=nh-1
      ifact=0

c ---------------------------- grid1:    "hx" used
      call grid1(2,nh1)
c ---------------------------- grid1:     NEW "nh", "hx" and "dh" redistributed

c ---------------------------- extrap:  "hx" and "dh" used
      call extrap
c ---------------------------- extrap:   NEW "dtime", "dh" set to 0.

c ---------------------------- profil:  "hx" used
    4 call profil(outbur)
c ---------------------------- profil:   NEW "dh" calculated (H and He-4, only)

      acc=1.0d-8
      iterat=0
      iterm=20
      icom(11)=0
c                              * * * * * *        --- begin of iteration loop
    1 continue
      iterat=iterat+1
      if(iterat.gt.iterm)go to 2

c ---------------------------- hsolve:  "hx" + "dh" used
      call hsolve(dismax,delmax)
c ---------------------------- hsolve:         "dh" corrected

      write(*    ,100)model,iterat,dismax,delmax
      if(dismax.gt.acc.or.delmax.gt.acc) go to 1

c                              * * * * * *        --- end of iteration loop

c ---------------------------- addcor:  "hx" = "hx" + "dh"
      call addcor
c ---------------------------- addcor:  time = time + dtime
c                              addcor: luminosity re-normalized

c  . . . . . . . . . . . . . .           begin of screen & disk  output:
      tilog=0.
      dtlog=0.
      if (time.ne.0.) tilog=dlog10(time/year)
      if (dtime.ne.0.)dtlog=dlog10(dtime/year)
      rhclg=hx(1,1)
      tclg=hx(2,1)
      telg=hx(1,nh)
      fllg=hx(2,nh)
      xc =outbur(1)
      fmc=outbur(2)*hx(mm,nh)
      if(outbur(3).ne.0.d0)fmc=outbur(3)*hx(mm,nh)
      gogcen=outbur(5)

      teflog = telg + dlog10(2.d0)/4
      do 60 k=nh-1,3,-1
      kk=k
      if(hx(2,k).lt.teflog.and.hx(2,k-1).gt.teflog)go to 61
   60 continue
      stop 'Teff not found'
   61 fact=(teflog-hx(2,kk))/(hx(2,kk-1)-hx(2,kk))
      refect=(1.d0-fact)*hx(3,kk)+fact*hx(3,kk-1)
      atm=dlog10(10.**hx(3,nh-1)/10.**refect-1.d0)
c ===
c model   = sequential model number
c iterat  = number of iterations
c nh      = number of mass points
c tilog   = log10 (time)                               [yr]
c dtlog   = log10 (time step)                          [yr]
c telg    = log10 ( surface temperature )              [K]
c fllg    = log10 ( surface luminosity / solar luminosity )
c tclg    = log10 ( central temperature )              [K]
c rhclg   = log10 ( central density )                  [g/cm**3]
c xc      = central hydrogen content
c fmc     = convective or helium core mass             [solar mass]
c gogcen  = central radiative/adiabatic gradient ratio
c atm     = log10 (relative geometrical thickness of optically thin layers)

      write(iout2,101)model,iterat,nh,tilog,dtlog,telg,fllg,tclg,rhclg,
     *            xc,fmc,gogcen,atm
  307 format(/,' mod it  nh   tilog dtlog  telg    fllg   tclg',
     *'  rhclg     xc     fmc    cen  atm')
      write(*,307)
      write(*    ,101)model,iterat,nh,tilog,dtlog,telg,fllg,tclg,rhclg,
     *            xc,fmc,gogcen,atm
      write(*,101)
c ===
      if (istep/nwrmod*nwrmod.eq.istep) call writeh
c  . . . . . . . . . . . . . .             end of screen & disk  output

c ..................................   --- end of  N E W   M O D E L
      go to 3

    2 continue
c ---------------------------
  102 FORMAT(1X,' Henyey iterations DO NOT CONVERGE:',
     *          ' DECREASE the TIME STEP,  IFACT =',I2)
      ifact=ifact+1
      write(*    ,102)ifact
      write(iout2,102)ifact
      if (ifact.gt.5)stop ': file "h03.f" - time step
     * decreased many times !...'
      dtime=dtime*0.5d0
c                           --- set all time changes = 0
      do 44 k=1,nh
      do 33 i=1,m2
         dh(i,k)=0.
   33 continue
   44 continue
      go to 4
c ---------------------------

    3 continue

  100 format('istep,iter,dismax,delmax=',2i5,2f20.15)
  101 format(i4,i3,i5,f8.4,f5.2,f7.4,f8.4,2f7.4,f9.6,f7.4,f5.2,f5.1)

c                                    --- end of main program
      stop '=> h03 - okay'
      end


c..............................................................................
c..............................................................................
c
      subroutine extrap
      implicit double precision (a-h,o-z)
      include 'data.03'
      common/hen/ hx(m2,nt), dh(m2,nt), nh
c
c subject:
c     establish new time step. Put time changes to the model, "dh" = 0.
c output:
c common /heninc/
c     dtime = "old" time step on input and "new" on output     [s]
c     model = sequential model no
c     fmtot = total stellar mass                               [g]
c     flunit= luminosity unit                                  [erg/s]
c common /hen/
c     nh - 2   = number of mass points
c     time step extrapolation is basing on time changes "dh (1..m2,2..nh-1)"
c     (excluding "dh(mm,2..nh-1)) to the model.
c     On output all the time changes are set to 0.
c predefined in 'data.03'
c     sunm
c     mm, m2
c     m2       is also dimension of the vector "st"
c     (st1..st8)*facst = limits for time differences
c............................................................................

      dimension st(m2)

      nh1=nh-1

      st(1)=st1 * facst
      st(2)=st2 * facst
      st(3)=st3 * facst
      st(4)=st4 * facst
      st(5)=st5 * facst
      st(6)=st6 * facst
      st(7)=st7 * facst
      st(8)=st8 * facst
c                           --- find maximum time difference in units of "st"
      ftinv=0.d0
      do 2 k=2,nh1
      do 1 i=1,m2
         if (i.ne.mm) ft=dabs(dh(i,k)/st(i))
         if (ft.gt.ftinv) then
            ftinv=ft
         end if
    1 continue
    2 continue
c                           --- if all time changes are 0 then time step will
c                               not be changed
      if (dabs(ftinv).lt.1.d-99) ftinv=1.
c                           --- if central hydrogen abundance > "endhra" then
c                               it should not fade faster than by half of value
c                               during single time step. For convective cores
c                               "endhra" may have to be replaced with "endhco".
      xcen=hx(mm+1,2)
      if(xcen.gt.endhra) then
         dxcen=dabs(dh(mm+1,2))
         ftlim=2.d0*dxcen/xcen
         if (ftinv.lt.ftlim) ftinv=ftlim
      end if

      facdt=1.d0/ftinv
c                           --- growth of time step will be limited to the
c                               factor 1.5 (but cutting down is off limits)
      if (facdt.gt.1.5) facdt=1.5d0
c                           --- estimate the new time step:
      dtime=dtime*facdt

c                           --- apply main sequence time scale of evolution
c                               if the time step not established yet
      if (dabs(dtime).lt.1.d-99.and.model.gt.1) then
         dtime=1.0e7/(fmtot/sunm)**2*year
      end if
c                           --- set all time changes = 0
c                               (they have been added earlier)
      do 4 k=1,nh
      do 3 i=1,m2
         dh(i,k)=0.d0
    3 continue
    4 continue
c                           --- end of "extrap"
      end
c
c............................................................................
c
      subroutine profil(out)
      implicit double precision (a-h,o-z)
      include 'data.03'
      dimension out(20)
      common/hen/ hx(m2,nt), dh(m2,nt), nh
c
c subject:
c     hydrogen burning.
c     hydrogen mixing in the convective core
c output:
c     out(1)= xcen   = new central hydrogen content
c     out(2)= fmccor = convective core mass / total stellar mass
c     out(3)= fmheco =     helium core mass / total stellar mass
c     out(5)= gogcen = central radiative/adiabatic gradient ratio
c----------------------------------------------------------------------------
c IMPORTANT: future ("new") convective core and chemical profiles
c            (for H and He-4 only) are predicted on the basis of actual model.
c            Boundary of the convective core is found by applying
c            the Schwarzschild stability criterion.
c----------------------------------------------------------------------------
c predefined in 'data.03'
c     mm, m2
c common /heninc/
c     dtime
c............................................................................

      dimension xvec(nx), dxtvec(nx), xx(m2), yy(m2), extra(20)
      dimension xvecme(nx), dxtmea(nx)

      nh1=nh-1
c                                    --- set some default values
      gogcen=0.d0
      fmccor=0.d0
      fmheco=0.d0
c                                    --- search for helium core boundary
      khecor=1
    1 khecor=khecor+1
      if (khecor.eq.nh1) stop ': file "h03.f" - subroutine "profil"
     * - pure helium star ...'
      if (hx(mm+1,khecor).lt.endhco) go to 1
      khecor=khecor-1
      if (khecor.gt.1) then
      fmheco=0.5d0*(dexp(hx(mm,khecor))/(1.d0+dexp(hx(mm,khecor)))
     $           +dexp(hx(mm,khecor+1))/(1.d0+dexp(hx(mm,khecor+1))))
      go to 10
      end if
c                                    --- no helium core

c                                    --- check for convection in central region
      kccore=2
      do 2 i=1,m2
         xx(i)=hx(i,2)
   2  continue
      call rhsh(xx,yy,0.d0,0.d0,0.d0,extra)
      gogcen=extra(20)
      if (gogcen.le.1.) go to 10
c                                    --- there is convection in central region.
c                                        Search for convective core boundary
c                                        by applying the Schwarzschild
c                                        stability criterion.
      do 4 k=2,nh1
         do 3 i=1,m2
            xx(i)=hx(i,k)
    3    continue
         call rhsh(xx,yy,0.d0,0.d0,0.d0,extra)
         kccore=k
         if (extra(20).lt.1.) go to 5
    4 continue
      if (kccore.eq.nh1) stop ': file "h03.f" - subroutine "profil"
     * - fully convective star ...'
c                                   --- convective core boundary is at the mass
c                                       point where radiative temp. gradient /
c                                       adiabatic temp. gradient is less than 1
c                                       for the first time.

    5 fmccor=dexp(hx(mm,kccore))/(1.d0+dexp(hx(mm,kccore)))
c                                   --- calculation of mean content of elements
c                                       in the convective core
      do 71 i=1,nx
      xvecme(i)=0.d0
   71 continue

      do 81 k=2,kccore

      do 72 i=1,nx
      xvec(i)=hx(mm+i,k)
      if(xvec(i).lt.0.d0)stop ': file "h03.f" - subroutine "profil"
     * - negative chemical abundance(s) ...'
   72 continue

      dm=0.5d0*(dexp(hx(mm,k+1))-dexp(hx(mm,k-1)))
      dm=dm/(1.d0+dexp(hx(mm,k+1)))/(1.d0+dexp(hx(mm,k-1)))
      if (k.eq.2) dm=0.5d0*(dexp(hx(mm,2))/(1.d0+dexp(hx(mm,2)))
     $                     +dexp(hx(mm,3))/(1.d0+dexp(hx(mm,3))))
      if(k.eq.kccore)dm=0.5d0*(dexp(hx(mm,kccore))-
     $   dexp(hx(mm,kccore-1)))
     $/(1.d0+exp(hx(mm,kccore)))/(1.d0+dexp(hx(mm,kccore-1)))

      do 73 i=1,nx
      xvecme(i) = xvecme(i) + xvec(i)*dm
   73 continue
   81 continue

      do 74 i=1,nx
      xvecme(i) = xvecme(i)/fmccor
   74 continue
      if(xvecme(1).lt.endhco)xvecme(1)=0.d0
c                                   --- calculation of mean burning rates
c                                       in the convective core
      do 75 i=1,nx
      dxtmea(i)=0.d0
   75 continue
      ppmea=0.d0
      cnmea=0.d0

      do 8 k=2,kccore
         rh=10.d0**hx(1,k)
         t =10.d0**hx(2,k)
         call nburn(rh,t,xvecme,epsx,dxtvec)

      dm=0.5d0*(dexp(hx(mm,k+1))-dexp(hx(mm,k-1)))
      dm=dm/(1.d0+dexp(hx(mm,k+1)))/(1.d0+dexp(hx(mm,k-1)))
      if (k.eq.2) dm=0.5d0*(dexp(hx(mm,2))/(1+dexp(hx(mm,2)))
     $                     +dexp(hx(mm,3))/(1+dexp(hx(mm,3))))
      if(k.eq.kccore)dm=0.5d0*(dexp(hx(mm,kccore))-
     $   dexp(hx(mm,kccore-1)))
     $/(1.d0+dexp(hx(mm,kccore)))/(1.d0+dexp(hx(mm,kccore-1)))

         do 76 i=1,nx
         dxtmea(i) = dxtmea(i) + dxtvec(i)*dm
   76    continue
         ppmea=ppmea+com(19)*dm
         cnmea=cnmea+com(20)*dm
    8 continue

      do 77 i=1,nx
      dxtmea(i) = dxtmea(i)/fmccor
   77 continue
      ppmea=ppmea/fmccor
      cnmea=cnmea/fmccor
c
c Notice: "ppmea" and "cnmea" are > 0., "cnmea" HAS ALWAYS to be > 0.
c
      if(xvecme(1).lt.endhco)then
        xnew = 0.d0
      else
      if(dabs(xvecme(1)/2.).gt.dabs(dxtmea(1)*dtime))then
        xnew = xvecme(1) + dxtmea(1)*dtime
      else
c------------------------------------------------
c  analytical integration of eq.: - dx/dt=ppmea*x**2+cnmea*x,(x=hydr. content),
c  ppmea, cnmea are assumed to be, approximately, independent on "x".

        xmea = xvecme(1)
        if(dtime*(ppmea*xmea+cnmea).lt.1.d-16)then
          xnew=xmea
        else
          if(cnmea*dtime.gt.1.d-8)then
            xnew=cnmea/((ppmea+cnmea/xmea)*dexp(cnmea*dtime)-ppmea)
          else
            xnew=xmea/(1.d0+dtime*(ppmea*xmea+cnmea))
          end if
        end if
c------------------------------------------------
      end if
      end if
      if (xnew.lt.endhco) xnew=0.d0

      do 9 k=2,kccore
         dh(mm+1,k) = xnew - hx(mm+1,k)
         dh(mm+2,k) = - dh(mm+1,k)
    9 continue
c                                   --- hydrogen mixing and burning in the
c                                       convective core has been calculated

c                                   --- calculate changes of chemical compo-
c                                       sition outside the convective or helium
c                                       core
   10 continue
      k1=2
      if (kccore.gt.2) k1=kccore+1
      if (khecor.ge.2) k1=khecor+1

      do 12 k=k1,nh1
         rh=10.**hx(1,k)
         t=10.**hx(2,k)
         do 11 i=1,nx
            xvec(i)=hx(mm+i,k)
   11    continue
         if(xvec(1).lt.endhco)xvec(1)=0.d0
         call nburn(rh,t,xvec,epsx,dxtvec)
         x=xvec(1)
         aq=com(19)
         bq=com(20)
c
c Notice: "aq" and "bq" are > 0., "bq" HAS ALWAYS to be > 0.
c
         if(x.lt.endhco)then
           dx=0.d0
         else
           if(dabs(dxtvec(1)*dtime).lt.dabs(x/2.))then
             dx=dxtvec(1)*dtime
           else
c------------------------------------------------
c  analytical integration of eq.: - dx/dt=aq*x**2+bq*x, (x=hydr. content)
c  aq, bq are assumed to be, approximately, independent on "x".

             if(bq+aq.lt.1.d-29.or.dtime*(aq*x+bq).lt.1.d-12)then
               dx=0.d0
             else
               if(bq*dtime.gt.1.d-8)then
                 dx=bq/((aq+bq/x)*dexp(bq*dtime)-aq)-x
               else
                 dx=x/(1.d0+dtime*(aq*x+bq))-x
               end if
             end if
c------------------------------------------------
           end if
           if(x+dx.lt.endhra)dx=-x
         end if
         if(dx.gt.0.d0)stop ': file "h03.f" - subroutine "profil"
     * - delta(X) > 0'

         dh(mm+1,k)=dx
         dh(mm+2,k)=-dx
   12 continue
c                                   --- chemical changes has been calculated

c                                   --- store boundary values
      dh(mm+1,1)=dh(mm+1,2)
      dh(mm+2,1)=dh(mm+2,2)
c                                   --- set outgoing data

      do 13 i=1,20
         out(i)=0.d0
   13 continue
      out(1) = hx(mm+1,1)+dh(mm+1,1)
      out(2) = fmccor
      out(3) = fmheco
      out(4) = 0.d0
      out(5) = gogcen
c                                   --- end of "profil"
      end
c
c.............................................................................
c
      subroutine hsolve(dismax,delmax)
      implicit double precision (a-h,o-z)
      include 'data.03'
      common/hen/ hx(m2,nt), dh(m2,nt), nh
c
c subject:
c     calculate the inner and the outer boundary conditions,
c     calculate all difference equations
c     calculate and adds the corrections "dx" to all time changes "dh"
c input:
c     the array of variables "hx" and time changes "dh"
c     at all mass points, stored in common/hen/
c output:
c     the array of variables "hx" and time changes "dh"
c     at all mass points, stored in common/hen/
c           with the corrections "dx" that have been added to "dh"
c     dismax = the largest discrepancy in the difference equations
c     delmax = the largest correction "dx"
c                     both in units of mass zone sizes, "ss"
c predefined in 'data.03'
c     nt, m2, nss, mm, nx
c............................................................................

c           --- below table declarations form a common block due to SVS/386 1.0
c               fortran compiler needs. This common is not used outside of this
c               subroutine and nothing must be saved from call to call.
c
      common /comhso/ ss(nss), s(4,9), cx(4,3,nt), dx(4,nt),
     *       fobc(4,3), fibc(4,3), sk(4,5), sk1(4,5), xvec(nx),
     *       resenv(20), resen1(20), xx(m2), xp(m2)


      dtinv=0.d0
      if (dtime.ne.0.)dtinv=1.0d0/dtime

c                                 === CONSTANTS ===
c                                 --- convert the constants:
      nh1=nh-1
c                     ---  ss(i), i=1,2,3,4 = the maximum step size
c                          in four variables. Defined in 'data.03'
      ss(1)=ss1
      ss(2)=ss2
      ss(3)=ss3
      ss(4)=ss4

c                                 === END OF CONSTANTS ===

c                      === BEGINNING OF OUTER BOUNDARY CONDITIONS ===
c                      --- calculate the array "fobc" defined as:
c                          fobc(i,1) = xx(i)-hx(i,nh1)-dh(i,nh1)
c                          fobc(i,2) = d ( xx(i) ) / dtelg
c                          fobc(i,3) = d ( xx(i) ) / dfllg
c                          where  i = 1, 2, 3, 4,    and:
c                          xx(i) = value of variable "i" at k=nh-1 as
c                                  calculated with subroutine "env"
c                          dtelg = perturbation of "telg", i.e. of the
c                                  log10 ( surface temperature)
c                          dfllg = perturbation of "fllg", i.e. of the
c                                  log10 ( surface luminosity )
      telg=hx(1,nh)+dh(1,nh)
      fllg=hx(2,nh)+dh(2,nh)
      starm=hx(mm,nh)
c                      --- prepare vector of chemical content at the outer edge
      do 40 i=1,nx
      xvec(i)=hx(mm+i,nh1) + dh(mm+i,nh1)
   40 continue

      dismax=0.d0
      kdisma=0
      idisma=0
      dlogt=0.d0
      dlogrh=0.d0

      call env(telg,fllg,starm,xvec,resenv)

      do 1 i=1,4
      fobc(i,1)=resenv(i)-hx(i,nh1)-dh(i,nh1)
      dis=abs(fobc(i,1))
      if (dismax.lt.dis) then
         dismax=dis
         kdisma=nh
         idisma=i
      end if
    1 continue

      del=1.0d-9

      call env(telg+del,fllg,starm,xvec,resen1)

      do 2 i=1,4
      fobc(i,2)=(resen1(i)-resenv(i))/del
    2 continue

      call env(telg,fllg+del,starm,xvec,resen1)

      do 3 i=1,4
      fobc(i,3)=(resen1(i)-resenv(i))/del
    3 continue

c                      --- the corrections to the surface parameters and the
c                          corrections to the variables at "nh-1" satify
c                          the following equations:
c  fobc(i,1)+fobc(i,2)*dtelg+fobc(i,3)*dfllg - dx(i,nh1) = 0  ,   i=1,2,3,4

c                     --- calculate the elements of the "cx" array defined with
c        dtelg=dx(1,nh)=cx(1,1,nh1)+cx(1,2,nh1)*dx(1,nh1)+cx(1,3,nh1)*dx(2,nh1)
c        dfllg=dx(2,nh)=cx(2,1,nh1)+cx(2,2,nh1)*dx(1,nh1)+cx(2,3,nh1)*dx(2,nh1)
c             dx(3,nh1)=cx(3,1,nh1)+cx(3,2,nh1)*dx(1,nh1)+cx(3,3,nh1)*dx(2,nh1)
c             dx(4,nh1)=cx(4,1,nh1)+cx(4,2,nh1)*dx(1,nh1)+cx(4,3,nh1)*dx(2,nh1)
c
c  the first two equations may be rewritten as:

c  [cx(1,1,nh1)+cx(1,2,nh1)*fobc(1,1)+cx(1,3,nh1)*fobc(2,1)]
c    + dtelg * [cx(1,2,nh1)*fobc(1,2)+cx(1,3,nh1)*fobc(2,2)-1]
c    + dfllg * [cx(1,2,nh1)*fobc(1,3)+cx(1,3,nh1)*fobc(2,3)]   = 0

c  [cx(2,1,nh1)+cx(2,2,nh1)*fobc(1,1)+cx(2,3,nh1)*fobc(2,1)]
c    + dtelg * [cx(2,2,nh1)*fobc(1,2)+cx(2,3,nh1)*fobc(2,2)]
c    + dfllg * [cx(2,2,nh1)*fobc(1,3)+cx(2,3,nh1)*fobc(2,3)-1] = 0

c  to fulfil this set of equations for each values of "dtelg" and "dfllg"
c  the expressions in square brackets have to be all equal zero.
c  These six equations allow us to estimate:

      det=fobc(1,2)*fobc(2,3)-fobc(2,2)*fobc(1,3)
      cx(1,1,nh1)=(fobc(2,1)*fobc(1,3)-fobc(1,1)*fobc(2,3))/det
      cx(1,2,nh1)=fobc(2,3)/det
      cx(1,3,nh1)=-fobc(1,3)/det
      cx(2,1,nh1)=(fobc(1,1)*fobc(2,2)-fobc(2,1)*fobc(1,2))/det
      cx(2,2,nh1)=-fobc(2,2)/det
      cx(2,3,nh1)=fobc(1,2)/det

      do 5 i=3,4
      do 4 j=1,3
      cx(i,j,nh1)=fobc(i,2)*cx(1,j,nh1)+fobc(i,3)*cx(2,j,nh1)
    4 continue
      cx(i,1,nh1)=cx(i,1,nh1)+fobc(i,1)
    5 continue
c                    --- the array "cx" at the point "nh-1" has been calculated
c                    === END OF OUTER BOUNDARY CONDITIONS ===

c                    --- calculate the right hand sides and their derivatives
c                        at k=nh-1:
      k=nh-1
      do 60 i=1,m2
      xx(i)=hx(i,k)+dh(i,k)
   60 continue
      dlogrh=dh(1,k)
      dlogt=dh(2,k)
      call deriv(xx,dtinv,dlogt,dlogrh,sk,plog,pt,pr)

c                     === BEGIN THE ZONES: k = nh-1, nh-2, ... 3, 2  ===
    6 continue

      do 62 i=1,4
      do 61 j=1,5
      sk1(i,j)=sk(i,j)
   61 continue
   62 continue
      plog1=plog
      pt1=pt
      pr1=pr

      k=k-1
c                     === CALCULATE THE "s" ARRAY FOR "k" , "k+1": ===
c
c                     --- the corrections "dx" have to satisfy the equations:
c                         (i=1,2,3,4)
c                         s(i,1) +
c                         + s(i,2)*dx(1,k) + s(i,3)*dx(2,k) + s(i,4)*dx(3,k) +
c                         + s(i,5)*dx(4,k) + s(i,6)*dx(1,k1) + s(i,7)*dx(2,k1)+
c                         + s(i,8)*dx(3,k1) + s(i,9)*dx(4,k1) = 0

      k1=k+1
      dm=hx(mm,k1)-hx(mm,k)
c                     --- the difference of values of an independent variable
c                         between "k+1" and "k" has been calculated

c                     --- calculate the right hand sides and their derivatives
      do 63 i=1,m2
      xx(i)=hx(i,k)+dh(i,k)
   63 continue
      dlogrh=dh(1,k)
      dlogt=dh(2,k)
      call deriv(xx,dtinv,dlogt,dlogrh,sk,plog,pt,pr)

c                     --- calculate the discrepancy between the two sides
c                         of the difference equations:

      s(1,1)=dm*(sk(1,1)+sk1(1,1))/2+plog-plog1
      dis=abs(s(1,1))
      if (dismax.lt.dis) then
         kdisma=k
         idisma=1
         dismax=dis
      end if

      do 11 i=2,4
      s(i,1)=dm*(sk(i,1)+sk1(i,1))/2.d0+hx(i,k)-hx(i,k1)+
     $   dh(i,k)-dh(i,k1)
      dis=dabs(s(i,1))
      if (dismax.lt.dis) then
         kdisma=k
         idisma=i
         dismax=dis
      end if
   11 continue

c                     --- "dismax" is the maximum discrepancy in units of "ss"

c                     --- calculate the derivatives of the discrepancies with
c                         respect to all variables in "k" and "k+1",
c                         and store them in the "s" array:
      do 13 j=1,4

      j1=j+1
      j5=j+5

      do 12 i=1,4
      s(i,j1)=dm*sk(i,j1)/2.d0
      s(i,j5)=dm*sk1(i,j1)/2.d0
   12 continue

      if(j.eq.1)go to 14
      s(j,j1)=s(j,j1)+1.d0
      s(j,j5)=s(j,j5)-1.d0

   14 continue
   13 continue

      s(1,2)=s(1,2)+pr
      s(1,3)=s(1,3)+pt
      s(1,6)=s(1,6)-pr1
      s(1,7)=s(1,7)-pt1

c                    === "s" ARRAY HAS BEEN CALCULATED ===

c                    === BEGIN CALCULATIONS OF THE "cx" ARRAY FOR "k","k+1" ===
c                    --- the arrray "s" has been calculated, it is  4 x 9 :
c                        equation #              <---- k --->  <--- k1 --->
c                                 1          y   +1  x  x  x   -1  x  x  x
c                                 2          y    x +1  x  x    x -1  x  x
c                                 3          y    0  x +1  0    0  x -1  0
c                                 4          y    x  x  0 +1    x  x  0 -1
c                        in this cartoon symbols have the following meaning:
c                        y = the "free" term, should be equal  0
c                        x = a non-zero element, proportional to zone thickness
c                       +1 = an element approximately equal +1
c                       -1 = an element approximately equal -1
c                        0 = an element that is identically equal  0

c                    --- the last array "cx" gives:
c                   dx(3,k1)=cx(3,1,k1)+cx(3,2,k1)*dx(1,k1)+cx(3,3,k1)*dx(2,k1)
c                   dx(4,k1)=cx(4,1,k1)+cx(4,2,k1)*dx(1,k1)+cx(4,3,k1)*dx(2,k1)

c                    --- eliminate "dx(3,k1)" and "dx(4,k1)" from the array "s"
      do 20 i=1,4
      s(i,1)=s(i,1)+s(i,8)*cx(3,1,k1)+s(i,9)*cx(4,1,k1)
      s(i,6)=s(i,6)+s(i,8)*cx(3,2,k1)+s(i,9)*cx(4,2,k1)
      s(i,7)=s(i,7)+s(i,8)*cx(3,3,k1)+s(i,9)*cx(4,3,k1)
   20 continue
c                    --- the array "s" is now  4 x 7 :
c                        y   +1  x  x  x   -1  x
c                        y    x +1  x  x    x -1
c                        y    0  x +1  0    x  x
c                        y    x  x  0 +1    x  x

      do 22 i=1,2
      fac3=s(i,4)/s(3,4)
      fac4=s(i,5)/s(4,5)
      do 21 j=1,7
      s(i,j)=s(i,j)-s(3,j)*fac3-s(4,j)*fac4
   21 continue
   22 continue

c                     --- the array "s" is now  4 x 7 :
c                         y   +1  x  0  0   -1  x
c                         y    x +1  0  0    x -1
c                         y    0  x +1  0    x  x
c                         y    x  x  0 +1    x  x

c                     --- solve the first two equations
c                         for "dx(1,k1)" and "dx(2,k1)":

      det=s(1,6)*s(2,7)-s(1,7)*s(2,6)

      do 25 m=1,2
      if(m.eq.1)fac1=-s(2,7)
      if(m.eq.1)fac2=s(1,7)
      if(m.eq.2)fac1=s(2,6)
      if(m.eq.2)fac2=-s(1,6)
      do 23 j=1,3
      cx(m,j,k)=(s(1,j)*fac1+s(2,j)*fac2)/det
   23 continue
   25 continue

c                     --- plug this solution into equations 3 and 4,
c                         and solve for "dx(3,k)" and "dx(4,k)"
      do 27 i=3,4
      do 26 j=1,3
      cx(i,j,k)=-(s(i,j)+cx(1,j,k)*s(i,6)+cx(2,j,k)*s(i,7))/s(i,i+1)
   26 continue
   27 continue

c                     === THE ARRAY "cx" FOR "k" , "k+1" HAS BEEN CALCULATED ==

      if(k.gt.2)go to 6

c                     === END OF THE ZONES: k = nh-1, nh-2, ... 3, 2 :

c                     --- we are at  k=2, the last array "cx" gives:
c                         dx(3,2)=cx(3,1,2)+cx(3,2,2)*dx(1,2)+cx(3,3,2)*dx(2,2)
c                         dx(4,2)=cx(4,1,2)+cx(4,2,2)*dx(1,2)+cx(4,3,2)*dx(2,2)
c                     --- the inner boundary will provide more equations

c                     === BEGINNING OF INNER BOUNDARY CONDITIONS: ===
c                     --- calculate the array "fibc" defined as:
c                         fibc(i,1) = xx(i)-hx(i,2)-dh(i,2)
c                         fibc(i,2) = d ( xx(i) ) / drhclg
c                         fibc(i,3) = d ( xx(i) ) / dtclg
c                         where  i = 1, 2, 3, 4,    and:
c                         xx(i) = value of variable "i" at   k=2
c                                 as calculated with subr. "center"
c                         dtclg = perturbation of "tclg", i.e. of the
c                                 log10 ( central temperature)
c                         drhclg= perturbation of ""rhclg, i.e. of the
c                                 log10 ( central density )
      rhclg=hx(1,1)+dh(1,1)
      tclg=hx(2,1)+dh(2,1)
      fmc=dexp(hx(mm,2))/(1.0d0+dexp(hx(mm,2)))
      dlogrh=dh(1,1)
      dlogt=dh(2,1)

      do 93 i=1,nx
      xvec(i)=hx(mm+i,2) + dh(mm+i,2)
   93 continue

      call center(rhclg,tclg,fmc,xvec,dtinv,dlogt,dlogrh,xx)

      do 28 i=1,4
      fibc(i,1)=xx(i)-hx(i,2)-dh(i,2)
      dis=dabs(fibc(i,1))
      if (dismax.lt.dis) then
         kdisma=1
         idisma=i
         dismax=dis
      end if
   28 continue

      del=1.0d-9

      call center(rhclg+del,tclg,fmc,xvec,dtinv,dlogt,dlogrh,xp)

      do 29 i=1,4
      fibc(i,2)=(xp(i)-xx(i))/del
   29 continue

      call center(rhclg,tclg+del,fmc,xvec,dtinv,dlogt,dlogrh,xp)

      do 30 i=1,4
      fibc(i,3)=(xp(i)-xx(i))/del
   30 continue


c                    --- the corrections to the central parameters and the
c                        corrections to the variables at  k=2  satify
c                       the following equations:
c                       fibc(i,1)+fibc(i,2)*drhclg+fibc(i,3)*dtclg - dx(i,2) =0
c                                                i=1,2,3,4

c                    --- calculate the elements of the "cx" array defined with:
c                  drhclg=dx(1,1)=cx(1,1,1)+cx(1,2,1)*dx(1,2)+cx(1,3,1)*dx(2,2)
c                   dtclg=dx(2,1)=cx(2,1,1)+cx(2,2,1)*dx(1,2)+cx(2,3,1)*dx(2,2)
c                         dx(3,2)=cx(3,1,1)+cx(3,2,1)*dx(1,2)+cx(3,3,1)*dx(2,2)
c                         dx(4,2)=cx(4,1,1)+cx(4,2,1)*dx(1,2)+cx(4,3,1)*dx(2,2)

      det=fibc(1,2)*fibc(2,3)-fibc(2,2)*fibc(1,3)
      cx(1,1,1)=(fibc(2,1)*fibc(1,3)-fibc(1,1)*fibc(2,3))/det
      cx(1,2,1)=fibc(2,3)/det
      cx(1,3,1)=-fibc(1,3)/det
      cx(2,1,1)=(fibc(1,1)*fibc(2,2)-fibc(2,1)*fibc(1,2))/det
      cx(2,2,1)=-fibc(2,2)/det
      cx(2,3,1)=fibc(1,2)/det

      do 32 i=3,4
      do 31 j=1,3
      cx(i,j,1)=fibc(i,2)*cx(1,j,1)+fibc(i,3)*cx(2,j,1)
   31 continue
      cx(i,1,1)=cx(i,1,1)+fibc(i,1)
   32 continue
c                    --- the array cx at the inner boundary has been calculated

c                    --- combine the two sets of equations expressing "dx(3,2)"
c                        and "dx(4,2)" in terms of "dx(1,2)" and "dx(2,2)":
      do 34 i=3,4
      do 33 j=1,3
      cx(i,j,1)=cx(i,j,1)-cx(i,j,2)
   33 continue
   34 continue

      det=cx(3,2,1)*cx(4,3,1)-cx(4,2,1)*cx(3,3,1)
      dx(1,2)=(cx(4,1,1)*cx(3,3,1)-cx(3,1,1)*cx(4,3,1))/det
      dx(2,2)=(cx(3,1,1)*cx(4,2,1)-cx(4,1,1)*cx(3,2,1))/det

c                     --- corrections: "dx(1,2)" and "dx(2,2)" have been found
      do 35 i=1,2
      dx(i,1)=cx(i,1,1)+cx(i,2,1)*dx(1,2)+cx(i,3,1)*dx(2,2)
   35 continue

c                     --- corr. to the central parameters have been calculated
c                     === END OF THE INNER BOUNDARY CONDITIONS ===

c                     === BEGIN CORRECTIONS: ===
c                     --- calculate all the corrections,
c                         look for the largest correction:
      delmax=dabs(dx(1,1)/ss(1))
      delm=dabs(dx(2,1)/ss(2))
      if(delmax.lt.delm)delmax=delm

      do 37 k=2,nh1
      k1=k+1
      dx(1,k1)=cx(1,1,k)+cx(1,2,k)*dx(1,k)+cx(1,3,k)*dx(2,k)
      dx(2,k1)=cx(2,1,k)+cx(2,2,k)*dx(1,k)+cx(2,3,k)*dx(2,k)
      dx(3,k) =cx(3,1,k)+cx(3,2,k)*dx(1,k)+cx(3,3,k)*dx(2,k)
      dx(4,k) =cx(4,1,k)+cx(4,2,k)*dx(1,k)+cx(4,3,k)*dx(2,k)

      do 36 i=1,4
      delm=dabs(dx(i,k)/ss(i))
      if(delmax.lt.delm)delmax=delm
   36 continue
   37 continue

      delm=dabs(dx(1,nh)/ss(2))
      if(delmax.lt.delm)delmax=delm
      delm=dabs(dx(2,nh)/ss(4))
      if(delmax.lt.delm)delmax=delm
c                     --- all the corrections have been calculated,
c                         the largest is "delmax"

c                     --- add corrections "dx" to all time changes "dh" :
      fac=1.0d0
      if((delmax*.3).gt.fac)fac=1.d0/(delmax*.3d0)

c                 ^                          ^
c            ---    Attention, please - relatively big corrections allowed.
c            ---    In the case of trubbles with convergence, you can try,
c            ---    as a first step, to increase the above indicated values.
      do 39 k=2,nh1

      do 38 i=1,4
      dh(i,k)=dh(i,k)+dx(i,k)*fac
   38 continue

   39 continue

c                     --- add corrections to the boundary parameters:
      dh(1,1)=dh(1,1)+dx(1,1)*fac
      dh(2,1)=dh(2,1)+dx(2,1)*fac
      dh(1,nh)=dh(1,nh)+dx(1,nh)*fac
      dh(2,nh)=dh(2,nh)+dx(2,nh)*fac

c                                           --- end of "hsolve"
      end
c
c..............................................................................
c
      subroutine deriv(xx,dtinv,dlogt,dlogrh,sk,plog,pt,pr)
      implicit double precision (a-h,o-z)
      include 'data.03'
      dimension xx(m2),sk(4,5)
c
c subject:
c     calculate right hand sides and the derivatives
c input:
c     xx(1)  = log10 (density)                      [g/cm**3]
c     xx(2)  = log10 (temperature)                  [K]
c     xx(3)  = log10 (radius)                       [cm]
c     xx(4)  = luminosity / luminosity unit
c     xx(5)  = log (mass / (total stellar mass - mass))
c     xx(6)  = x         hydrogen content                 [by mass fraction]
c     xx(7)  = y           helium content                 [by mass fraction]
c     xx(8)  = zn        nitrogen content for CNO burning [by mass fraction]
c     dtinv  = 1 / the current time step            [1/s]
c              ("dtime" of /heninc/ not used here)
c     dlogt  = change of log10 (temperature) during the current time step
c     dlogrh = change of log10 (density) during the current time step
c notice:
c     xx(i) already include the recent time change
c output:
c     sk(1,1)   = d plog / d xx(5)
c     sk(i,1)   = d xx(i) / d xx(5)    ,       i =    2, 3, 4
c     sk(i,j+1) = d sk(i,1) / d xx(j)  ,    i, j = 1, 2, 3, 4
c     plog      = log10 (pressure)                  [c.g.s.]
c     pt        = d ln p / d ln t     at constant density
c     pr        = d ln p / d ln rh    at constant temperature
c predefined in 'data.03'
c     m2  = number of variables in every mass point
c called:
c     rhsh     to calculate all the derivatives
c.............................................................................

      dimension yy(m2),xp(m2),yp(m2)
      dimension extra(20)

c                         --- calculate the right hand sides:

      call rhsh(xx,yy,dtinv,dlogt,dlogrh,extra)
      plog=extra(17)
      pt=extra(2)
      pr=extra(3)

      do 1 i=1,4
      sk(i,1)=yy(i)
    1 continue

c                         --- calculate the derivatives of the right hand sides
      del=1.0d-9
c                         --- "del" seems to be rather small...,  but it works!
      do 4 j=1,4

      do 2 i=1,m2
      xp(i)=xx(i)
    2 continue
      dlogtp=dlogt
      dlogrp=dlogrh

      xp(j)=xp(j)+del
      if(j.eq.1)dlogrp=dlogrp+del
      if(j.eq.2)dlogtp=dlogtp+del

      call rhsh(xp,yp,dtinv,dlogtp,dlogrp,extra)

      j1=j+1

      do 3 i=1,4
      sk(i,j1)=(yp(i)-yy(i))/del
    3 continue

    4 continue
c                         --- the array "sk" has been calculated
c                         --- end of "deriv"
      end
c
c..............................................................................
c
      subroutine rhsh(xx,yy,dtinv,dlogt,dlogrh,extra)
      implicit double precision(a-h,o-z)
      include 'data.03'
      dimension xx(m2),yy(m2),extra(20)
c
c subject:
c     calculate right hand sides of the stellar structure equations
c     including time-dependent terms.
c     However,  if dtinv=0. then time-dependent terms are neglected.
c input:
c     xx(1) = log10 (density)                      [g/cm**3]
c     xx(2) = log10 (temperature)                  [K]
c     xx(3) = log10 (radius)                       [cm]
c     xx(4) = luminosity / unit luminosity
c     xx(5) = log (mass / (total stellar mass - mass))
c     xx(6) = hydrogen content                     [by mass fraction]
c     xx(7) = helium content                       [by mass fraction]
c     xx(8..5+nx) = subsequent elements contents   [by mass fraction]
c     dtinv = 1 / the current time step            [1/s]
c     dlogt  = change of log10 (temperature) during the current time step
c     dlogrh = change of log10 (density)     during the current time step
c output
c     yy(1) = d log10 p / d xx(5)   attention: this is different than in "rhs"
c     yy(2) = d log10 t / d xx(5)
c     yy(3) = d log10 r / d xx(5)
c     yy(4) = d l/lunit / d xx(5)
c     yy(5) = d xx(5)   / d xx(5) = 1.0
c     yy(6)...yy(m2) = 0.
c     extra = number of additional functions
c     extra(1) = p     = pressure                  [c.g.s.]
c     extra(2) = pt    = d ln(p) / d ln(t)    (at constant density)
c     extra(3) = pr    = d ln(p) / d ln(ro)   (at constant temperature)
c     extra(4) = pgas  = gas pressure              [c.g.s.]
c     extra(5) = prad  = radiation pressure        [c.g.s.]
c     extra(6) = grad  = adiabatic temperature gradient
c     extra(7) = qt
c     extra(8) = qr
c     extra(9) = second helium ionization degree
c       (only if "termoh" is called, i.e. for t < te2, see "eostate.03")
c     extra(10)= (d u / d t)p   specific heat cap. at const. p
c       (only if "termoh" is called, i.e. for t < te2, see "eostate.03")
c     extra(11)= q =-(d ln ro / d ln t)p
c       (only if "termoh" is called, i.e. for t < te2, see "eostate.03")
c     extra(12)= not used
c     extra(13)= not used
c     extra(14)= not used
c     extra(15)= dxt   = d X / d t = hydrogen burning rate           [1/s]
c     extra(16)= drom  = d log10(rh) / d (mass fraction)
c     extra(17)= plog  = log10(p)
c     extra(18)= dpm   = d log10(p) / d (mass fraction)
c     extra(19)= grrad = radiative temperature gradient
c     extra(20)= grrad/grad = radiative/adiabatic gradient ratio
c called:
c     state    to calculate thermodynamic functions
c     opact    to calculate opacity (total: radiative and conductive)
c     nburn    to calculate energy generation rate in hydrogen burning
c              + neutrino energy loss rate (if hydrogen content = 0)
c     mixlen   to calculate effective temperature gradient
c              according to the "mixing length theory"
c common /heninc/
c     flunit  = luminosity unit, close to maximum or surface lumin.  [erg/s]
c     fmtot   = total stellar mass                                   [g]
c auxiliary variables:
c     p      = total pressure                                        [c.g.s.]
c     prad   = radiation pressure                                    [c.g.s.]
c     grad   = d ln t / d ln p     at constant entropy (adiabatic temperature
c                                  gradient)
c     grrad  = d ln t / d ln p     radiative temperature gradient
c     grt    = d ln t / d ln p     temperature gradient in the star
c     grrh   = d ln rh / d ln p     density gradient in the star
c     dpm    = d log10 p / d ( Mr / total mass ) pressure gradient in the star
c............................................................................

      dimension xvec(nx), dxtvec(nx)

      fcon1   = (flunit/fmtot)/(16.d0*pi*clight*gconst)
      fcon2   = gconst/4.d0/pi*fmtot**2/fln10
      fcon3   = fmtot/4.d0/pi/fln10
      fcon4   = fmtot/flunit

      rh=10.**xx(1)
      t=10.**xx(2)
      r=10.**xx(3)
      flr=xx(4)
      e5=dexp(xx(5))
      fmr=e5/(1+e5)

      do 1 i=1,nx
      xvec(i)=xx(mm+i)
    1 continue

      call opact(rh,t,xvec,fkap)
      call state(rh,t,xvec,extra)
      call nburn(rh,t,xvec,epsx,dxtvec)

      p    = extra(1)
      pt   = extra(2)
      pr   = extra(3)
      pgas = extra(4)
      prad = extra(5)
      grad = extra(6)
      qt   = extra(7)
      qr   = extra(8)

      plog=dlog10(p)
      grrad=fcon1*fkap*p/prad*flr/fmr
      grt=grad
      if(grt.gt.grrad)grt=grrad

c  ...........................................
      if (grrad.gt.grad) then
c                                 --- calculation of an effective temperature
c                                     gradient according to the "mixing length
c                                     theory" for t < te1.
c                                     Attention, please:
c                                     geometrical dillution NOT included
      if (t.lt.te1)
     $call mixlen(rh,t,r,fmr,p,extra(10),extra(11),fkap,grad,grrad,grt)

c                                 --- checking if covection is really adiabatic
c                                     in the intermediate zone, te1 < t < te2,
c                                     within relative accuracy of 0.001
      if (t.gt.te1.and.t.lt.te2) then
      grrad1=fcon1*com(14)*com(15)/prad*flr/fmr
      call mixlen(rh,t,r,fmr,com(15),com(16),com(17),com(14),com(18),
     $grrad1,grt12)
      if(dabs(grt12/com(18)-1.d0).gt.0.001)
     $stop 'nonadiabatic convection for t > te1, make "te1" higher'
      end if
      end if
c  ...........................................

      grrh=(1.d0-grt*pt)/pr
      dpm=-fcon2*fmr/r**4/p

      do 2 i=1,m2
      yy(i)=0.d0
    2 continue

      yy(1)=dpm
      yy(2)=grt*dpm
      yy(3)=fcon3/(rh*r**3)
      yy(4)=fcon4*epsx
      if (dtinv.ne.0.)
     *   yy(4)=yy(4)-fcon4*fln10*dtinv*(dlogt*qt-dlogrh*qr)

      do 3 i=1,m2
      yy(i)=yy(i)*e5/(1.d0+e5)**2
    3 continue

      yy(5)=1.d0

      drom=grrh*dpm
      extra(15) = dxtvec(1)
      extra(16) = drom
      extra(17) = plog
      extra(18) = dpm
      extra(19) = grrad
      extra(20) = grrad/grad
c                                           --- end of "rhsh"
      end
c
c..............................................................................
c
      subroutine addcor
      implicit double precision (a-h,o-z)
      include 'data.03'
      common/hen/ hx(m2,nt), dh(m2,nt), nh
c
c subject:
c     add time changes "dh" to the table "hx"
c     re-normalize luminosity (this is the only place where it is done)
c     add time step
c common /hen/
c     this is both input and output to this subroutine
c common /heninc/
c     flunit is re-normalized
c predefined in 'data.03'
c     m2,mm
c remark:
c     on output time changes are not set to 0. It is done  within
c     the subroutine "extrap".
c..............................................................................

      nh1=nh-1
c                                               --- find maximum luminosity
      flmax=0.d0
      do 1 k=2,nh1
         fl=hx(4,k) + dh(4,k)
         if (fl.gt.flmax) flmax=fl
    1 continue
      faclum=flmax
c                                               --- re-normalize luminosity
      do 2 k=2,nh1
         hx(4,k)=hx(4,k)/faclum
         dh(4,k)=dh(4,k)/faclum
    2 continue
      flunit=flunit*faclum
c                                               --- add time changes
      do 4 k=1,nh
      do 3 i=1,m2
         if (i.ne.mm) hx(i,k) = hx(i,k) + dh(i,k)
    3 continue
    4 continue
c                                               --- add time step
      time=time+dtime
c                                               --- end of "addcor"
      end
c
c..............................................................................

c this is a good place to include some extra code

      include 'eostate.03'
      include 'opacity.03'
      include 'nuclear.03'
      include 'commonsb.03'

      subroutine stupid
c This is stupid but NDP Fortran compiler needs here the "end" statement
      end
c..............................................................................
c        This file contains "STANDARD" STELLAR EVOLUTION PROGRAM "h03" -
c         main program and six subroutines called:
c             - "extrap" -
c             - "profil" -
c             - "hsolve" -
c             - "deriv"  -
c             - "rhsh"   -
c             - "addcor" -
c                                                    end of the file "h03.f"
ceeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
ceeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
ceeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee


