subroutine rhs(sr,YY,dydx)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (nmax=4,jmax=1999)
c
c Computes right hand sides of equations of stellar structure, for
c use in shooting calculations as well as in Newton-Rapheson iterations
c
c Note: composition information must be shipped over in
c       common/torhs/cmp(15)
c   along with sp and st, the changes in lnP and lnT from last model to now
c
      DIMENSION YY(nmax),dydx(nmax)
c Common block CTRL contains most of the control parameters for passing
      common/CTRL/dllenv,dltenv,tcut,dtfac,sstep(7),tstep(7),
     1            nmod,nsav,nprint,nzp,itstep
      COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR
c Common block STRUC contains model structure
      common/STRUC/h(jmax,30),comp(jmax,15),srray(jmax),nj,modno
c Common block GLOB contains global model parameters: xat0 is the inital X
      common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl
c Common TORHS contains needed stuff imported into this subroutine
c  things imported are compostion variables cmp, sp and st are changes
c  in P and T from last model to this one
      common/TORHS/cmpin(15),sp,st,ishell
c Common block EQRATE contains info for equilibrium hydrogen burning
      common/EQRATE/ieqep,ieqecn
c Common block MLT contains choice of mixing length theory parameters
      common/MLT/infmlt
      dimension cmp(15),dcdt(15)
c
      fmr=fms-sr
      R=dexp(YY(1))
      P=dexp(YY(2))
      T=dexp(YY(3))
      tl=dlog10(t)
      Fl=YY(4)
      XC12=cmpin(8)
      XC13=cmpin(9)
      XN14=cmpin(10)
      XN15=cmpin(11)
      XO16=cmpin(12)
      do 10 iin=1,15
        cmp(iin)=cmpin(iin)
10    continue
c
c For E.O.S. and opacity purposes, call hydrogen and helium abundances
c equal to the abundances of all the tracked isotopes
c
      x=cmp(1)+cmp(2)
      y=cmp(4)+cmp(3)
c
c      tcut=5.5d0
c      tcut=6.0d0
c call EOS if ishell not set
      if (ishell.lt.1.or.ishell.gt.nj) then
        call EOS(P,t,x,y,XC12,tcut,q,cp,ga,d,xd,xt,u,cv,eta,fmu)
c otherwise use model quantities
      else
        q=h(ishell,9)
        ga=h(ishell,6)
        d=dexp(h(ishell,5))
        xd=h(ishell,11)
        xt=h(ishell,10)
        cp=p*q/(d*t*ga)
      endif
c
c If the first model in a new sequence, RHS is expected to
c compute the relevant quasi-equlibrium abundances:
c Abundances: use the routines PPEQ, CNEQ to compute equilibrium 
c abundances for He3, Li7, C12, C13, N14, N15.
c
c      if (time.le.0.and.cmp(1).gt.0.d0) then
c        call XMASS(dlog10(d),tl,fmr,10.d0*dtime,cmp)
c      endif
c
c call opacity and nuke stuff if ishell not set
c
c      if (ishell.lt.1.or.ishell.gt.nj) then
      if (ishell.eq.0) then
        call OPACITY(dlog10(d),dlog10(T),x,y,Z,XC12,XO16,
     1             fk,0,DLKLT,DLKLD)
        call BNUKE(dlog10(d),dlog10(t),cmp,dcdt,eps,ieqep,epptot,
     1            ieqecn,ecno,e3a,eac,enu,0,dlelt,dleld)
      else
        fk=h(ishell,12)
        eps=h(ishell,15)
      endif
c
c Reset ishell to zero to allow RHS to work from scratch
c
      ishell=0
c
c Continuity equation gives dlnR/dSr
c
      CR=CMS/(CP4*CRS**3)
      dydx(1)=-CR/R**3/D
c
c Hydrostatic equilibrium equation gives dlnP/dSr
c
      CHSE=CG*CMS**2/(CP4*CRS**4)
      DPM=CHSE*FMR/P/R**4
      dydx(2)=dpm
c
c Thermal equilibrium gives dlnT/dSr
c
      if (fl.eq.0.d0) fl=fmr*eps/cls*cms
c find the radiative gradient
      Crad=3.D0*CLS/(16.D0*CP4*CG*CSB*CMS)
      GR=Crad*fl*fK*P/fmr/T**4
c check for convective stability
      IF(GR.LE.GA) THEN
c   stable against convection; set gradient to radiative gradient
        GRAD=GR
      ELSE
c   Convectively unstable; use mixing length theory to evaluate
c   the temperature gradient... if T is less than 1e7 K
        if (t.lt.1.d7) then
          call TGRAD(fmr,r,t,p,d,alfa,infmlt,q,cp,fk,ga,gr,grad)
        else
          grad=ga
        endif
      endif
c with the proper gradient, find dlnT/dSr
      dydx(3)=CHSE*FMR*GRAD/P/R**4
c Energy conservation equation (dFl/dSr); 
      clum=-cms/cls
      if (st.eq.0.d0.and.sp.eq.0.d0.and.dtime.eq.0.d0) then
        sbar=0.d0
      else
        sbar=P*q/d/ga*(st-ga*sp)
      endif
      dydx(4)=clum*(eps-sbar/dtime)
c      print *,fmr,t,dydx(1),dydx(2),dydx(3),dydx(4)
      RETURN
      END
c
      SUBROUTINE TGRAD(fmr,r,t,p,d,alfa,infmlt,q,cp,fkap,dela,delr,del)
      implicit double precision (a-h,o-z)
c
c Compute true temperature gradient using standard mixing-length theory
c a-la Cox and Giuli.  r and fmr in solar units.
c
c Common block CONST contains fundamental physical constants
      COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR
c
      G=(CG*CMS/CRS**2)*FMR/R**2
      HP=P/D/G
      FLAM=ALFA*HP
c
c Standard mixing length theory a la Cox and Giuli, but using the 
c parameters definitions of Fontaine, Villeneuve, and 
c Wilson (Ap. J. 243, 550).
c
      af=0.125d0
      bf=0.500d0
      cf=24.00d0
c if infmlt is 1 then use the numbers corresponding to Gilles' ML2
      if (infmlt.eq.1) then
        af=1.000d0
        bf=2.000d0
        cf=16.00d0
      endif
c
      A=dsqrt(q)*cp*fkap*g*d**2.5d0*flam**2/csb/t**3/dsqrt(p)
      A=A*dsqrt(af)/cf
      a0=3.d0/16.d0*bf*cf
      B=(a**2/a0*(delr-dela))**(1.d0/3.d0)
c
      y=1.0d0
      do 43 it=1,10
        y1=y
        f = y1 + b*y1**2 + a0*b**2*y1**3 - a0*b**2
        dfdy = 1.d0 + 2.d0*b*y + 3.d0*a0*b**2*y**2
        dy=-f/dfdy
        y=y1+dy
        if (dabs(dy/y).le.1.d-4) go to 44
43    continue
44    y=y**3
      del=delr*(1.d0-y)+dela*y
c
      return
      end
c
      subroutine XMASS(dl,tl,fmr,dtime,cmp)
      implicit double precision (a-h,m-z)
c
c Subroutine which takes in the value of temperature
c and mass fraction and outputs the composition at that point, as
c determined using specified rules
c
c This version: Quasi-equilibrium abundances for hydrogen 
c               burning reactions
c
c Common block EQRATE contains info for equilibrium hydrogen burning
      common/EQRATE/ieqep,ieqecn
      dimension cmp(15)

      if (tl.gt.6.0d0.and.cmp(1).gt.0.d0) then
        call PPEQ(dl,tl,dtime,cmp)
c        call CNEQ(dl,tl,dtime,cmp)
      endif
c
      return
      end