***
      SUBROUTINE mrenv(kw,mass,mt,mc,lum,rad,rc,aj,tm,ltms,lbgb,lhei,
     &                 rzams,rtms,rg,menv,renv,k2e)
      implicit none
      integer kw
      real*8 mass,mt,mc,lum,rad,rc,aj,tm
      real*8 k2e,menv,menvg,menvt,menvz,renv,renvg,renvt,renvz,rzr
      real*8 A,B,C,D,E,F,x,y
      real*8 k2bgb,k2g,k2z,logm,logmt,lbgb,ltms,lhei,rg,rtms,rzams
      real*8 teff,tebgb,tetms,tau,tauenv,tautms
*
* A function to estimate the mass and radius of the convective envelope,
* as well as the gyration radius of the envelope.
* N.B. Valid only for Z=0.02!
*
* The following input is needed from HRDIAG:
*   kw = stellar type
*   mass = zero-age stellar mass
*   mt = actual mass
*   mc = core mass (not really needed, can also be done outside subroutine)
*   lum = luminosity
*   rad = radius
*   rc = core radius (not really needed...)
*   aj = age
*   tm = main-sequence lifetime
*   ltms = luminosity at TMS, lums(2)
*   lbgb = luminosity at BGB, lums(3)
*   lhei = luminosity at He ignition, lums(4)
*   rzams = radius at ZAMS
*   rtms = radius at TMS
*   rg = giant branch or Hayashi track radius, approporaite for the type. 
*        For kw=1 or 2 this is radius at BGB, and for kw=4 either GB or 
*        AGB radius at present luminosity.
*
      logm = log10(mass)
      A = MIN(0.81d0,MAX(0.68d0,0.68d0+0.4d0*logm))
      C = MAX(-2.5d0,MIN(-1.5d0,-2.5d0+5.d0*logm))
      D = -0.1d0
      E = 0.025d0
*
* Zero-age and BGB values of k^2.
*
      k2z = MIN(0.21d0,MAX(0.09d0-0.27d0*logm,0.037d0+0.033d0*logm))
      if(logm.gt.1.3d0) k2z = k2z - 0.055d0*(logm-1.3d0)**2
      k2bgb = MIN(0.15d0,MIN(0.147d0+0.03d0*logm,0.162d0-0.04d0*logm))
*
      if(kw.ge.3.and.kw.le.6)then
*
* Envelope k^2 for giant-like stars; this will be modified for non-giant
* CHeB stars or small envelope mass below.
* Formula is fairly accurate for both FGB and AGB stars if M <= 10, and
* gives reasonable values for higher masses. Mass dependence is on actual
* rather than ZA mass, expected to work for mass-losing stars (but not
* tested!). The slightly complex appearance is to insure continuity at 
* the BGB, which depends on the ZA mass.
*
         logmt = log10(mt)
         F = 0.208d0 + 0.125d0*logmt - 0.035d0*logmt**2
         B = 1.0d+04*mt**(3.d0/2.d0)/(1.d0+0.1d0*mt**(3.d0/2.d0))
         x = ((lum-lbgb)/B)**2
         y = (F - 0.033d0*log10(lbgb))/k2bgb - 1.d0
         k2g = (F - 0.033d0*log10(lum) + 0.4d0*x)/(1.d0+y*(lbgb/lum)+x)
      elseif(kw.eq.9)then
*
* Rough fit for for HeGB stars...
*
         B = 3.0d+04*mt**(3.d0/2.d0)
         x = (MAX(0.d0,lum/B-0.5d0))**2
         k2g = (k2bgb + 0.4d0*x)/(1.d0 + 0.4d0*x)
      else
         k2g = k2bgb
      endif
*
      if(kw.le.2)then
         menvg = 0.5d0
         renvg = 0.65d0
      elseif(kw.eq.3.and.lum.lt.3.d0*lbgb)then
*
* FGB stars still close to the BGB do not yet have a fully developed CE.
*
         x = MIN(3.d0,lhei/lbgb)
         tau = MAX(0.d0,MIN(1.d0,(x-lum/lbgb)/(x-1.d0)))
         menvg = 1.d0 - 0.5d0*tau**2
         renvg = 1.d0 - 0.35d0*tau**2
      else
         menvg = 1.d0
         renvg = 1.d0
      endif
*
      if(rad.lt.rg)then
*
* Stars not on the Hayashi track: MS and HG stars, non-giant CHeB stars,
* HeMS and HeHG stars, as well as giants with very small envelope mass.
*
         
         if(kw.le.6)then
*
* Envelope k^2 fitted for MS and HG stars.
* Again, pretty accurate for M <= 10 but less so for larger masses.
* [Note that this represents the whole star on the MS, so there is a 
* discontinuity in stellar k^2 between MS and HG - okay for stars with a 
* MS hook but low-mass stars should preferably be continous...]
*
* For other types of star not on the Hayashi track we use the same fit as 
* for HG stars, this is not very accurate but has the correct qualitative 
* behaviour. For CheB stars this is an overestimate because they appear
* to have a more centrally concentrated envelope than HG stars.
* NOTE: Max on radius ratio added to stop large k2e (5/5/2009, JH)
*
            rzr = MAX(1.d0,rad/rzams)
            k2e = (k2z-E)*rzr**C + E*rzr**D
         elseif(kw.eq.7)then
* Rough fit for naked He MS stars.
            tau = aj/tm
            k2e = 0.08d0 - 0.03d0*tau
         elseif(kw.le.9)then
* Rough fit for HeHG stars.
            k2e = 0.08d0*rzams/rad
         endif
*
* tauenv measures proximity to the Hayashi track in terms of Teff.
* If tauenv>0 then an appreciable convective envelope is present, and
* k^2 needs to be modified.
*
         if(kw.le.2)then
            teff = sqrt(sqrt(lum)/rad)
            tebgb = sqrt(sqrt(lbgb)/rg)
            tauenv = MAX(0.d0,MIN(1.d0,(tebgb/teff-A)/(1.d0-A)))
         else
            tauenv = MAX(0.d0,MIN(1.d0,(sqrt(rad/rg)-A)/(1.d0-A)))
         endif
*
         if(tauenv.gt.0.d0)then
            menv = menvg*tauenv**5
            renv = renvg*tauenv**(5.d0/4.d0)
            if(kw.le.1)then
* Zero-age values for CE mass and radius.
               x = MAX(0.d0,MIN(1.d0,(0.1d0-logm)/0.55d0))
               menvz = 0.18d0*x + 0.82d0*x**5
               renvz = 0.4d0*x**(1.d0/4.d0) + 0.6d0*x**10
               y = 2.d0 + 8.d0*x
* Values for CE mass and radius at start of the HG.
               tetms = sqrt(sqrt(ltms)/rtms)
               tautms = MAX(0.d0,MIN(1.d0,(tebgb/tetms-A)/(1.d0-A)))
               menvt = menvg*tautms**5
               renvt = renvg*tautms**(5.d0/4.d0)
* Modified expressions during MS evolution.
               tau = aj/tm
               if(tautms.gt.0.d0)then
                  menv = menvz + tau**y*menv*(menvt - menvz)/menvt
                  renv = renvz + tau**y*renv*(renvt - renvz)/renvt
               else
                  menv = 0.d0
                  renv = 0.d0
               endif
               k2e = k2e + tau**y*tauenv**3*(k2g - k2e)
            else
               k2e = k2e + tauenv**3*(k2g - k2e)
            endif
         else
            menv = 0.d0
            renv = 0.d0
         endif
      else
*
* All other stars should be true giants.
*
         menv = menvg
         renv = renvg
         k2e = k2g
      endif
*
      menv = menv*(mt - mc)
      renv = renv*(rad - rc)
      menv = MAX(menv,1.0d-10)
      renv = MAX(renv,1.0d-10)
*
      return
      end
***