SUBROUTINE RYBHEQ
C     =================
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      common/grdpra/GRD(MDEPTH),pra(mdepth),pgs0(mdepth),ANTP(MDEPTH)
      common/rybpgs/CS(MDEPTH),PRAD2D(MDEPTH),F1HE
      dimension pradfc(mdepth)
c
      DO ID=1,ND
         GRD(ID)=0.
         pra(id)=0.
         pradt(id)=0.
         pradfc(id)=0.
      END DO
      PRD0=0.
c
      if(ifprad.gt.0) then
         CALL OPAINI(1)
         DO IJ=1,NFREQ
            CALL OPACF1(IJ)
            CALL RTEFR1(IJ)
            if(.not.lskip(1,ij))
     *      GRD(1)=GRD(1)+W(IJ)*ABSO1(1)*FH(IJ)*RAD1(1)
            DO ID=2,ND
               if(.not.lskip(id,ij)) then
               GRD(ID)=GRD(ID)+(RAD1(ID)*FAK1(ID)-
     *                 RAD1(ID-1)*FAK1(ID-1))*W(IJ)
               pra(id)=pra(id)+RAD1(ID)*FAK1(ID)*W(IJ)
               end if
            END DO
            pra(1)=pra(1)+RAD1(1)*FAK1(1)*W(IJ)-ABSO1(1)*W(IJ)*
     *             (RAD1(1)*FH(IJ)-HEXTRD(IJ))
         END DO
         if(idisk.eq.0) then
            GRD(1)=PCK*GRD(1)/DENS(1)
          else
            grd(1)=pck*grd(2)
         end if
         PRD0=PRD0*DENS(1)*DM(1)*PCK
c
         do id=1,nd
            pra(id)=pra(id)*pck
            pradt(id)=pra(id)
            pradfc(id)=pra(id)/(2.5213e-15*temp(id)**4)
         end do
C
         do id=2,nd-1
            dmm=un/(dm(id)-dm(id-1))
            dmp=un/(dm(id+1)-dm(id))
            dm0=two/(dm(id+1)-dm(id-1))
            qq=((pra(id+1)-pra(id))*dmp-(pra(id)-pra(id-1))*dmm)*dm0
            prad2d(id)=qq
         end do
         prad2d(1)=prad2d(2)
         prad2d(nd)=0.
      end if
c
      if(idisk.eq.0) then
         PGS0(1)=DM(1)*(GRAV-GRD(1))
         DO ID=2,ND
            PGS0(ID)=PGS0(ID-1)-PCK*GRD(ID)+GRAV*(DM(ID)-DM(ID-1))
         END DO
       else
         pgs0(1)=pgs(1)
         do id=2,nd
            grv=(dm(id)-dm(id-1))*qgrav*(zd(id)+zd(id-1))*half
            pgs0(id)=pgs0(id-1)-pck*grd(id)+grv
         end do
      end if
      if(iprybh.gt.0) write(6,603) iter
  603 format(/' rybheq in iter',i4,
     *     '  dm,dens,pgs0,pra,pradfc,grd,grad,ggrav'/)
c
      if(iprybh.gt.0) then
         do id=1,nd
            if(id.eq.1) then
               gra=grd(id)
               ggr=qgrav*zd(id)
             else
               gra=pck*grd(id)/(dm(id)-dm(id-1))
               ggr=qgrav*(zd(id)+zd(id-1))*half
            end if
            if(idisk.eq.0) then
               write(6,604) id,dm(id),dens(id),pgs0(id),pra(id),
     *         pradfc(id),grd(id),gra
             else
               write(6,604) id,dm(id),dens(id),pgs0(id),pra(id),
     *         pradfc(id),grd(id),gra,ggr,zd(id)
            end if
         end do
      end if
  604 format(i4,1p9e13.4)
C
      if(idisk.eq.0) then
         DO ID=1,ND
            T=TEMP(ID)
            pgs(id)=pgs0(id)
            AN=PGS(ID)/BOLK/T
            CALL ELDENS(ID,T,AN,ANE,ENRG,ENTT,WM,1)
            RHO=WMM(ID)*(AN-ANE)
            DENS(ID)=RHO
            ELEC(ID)=ANE
            CALL WNSTOR(ID)
            CALL STEQEQ(ID,POP,1)
c            IF(.NOT.LCHC.and.iter.lt.ielcor) CALL ELCOR(ID)
         END DO
       else
c
c        --------------------------
c        the rest is only for disks
c        --------------------------
c
         cprad=2.5213e-15
         do id=1,nd
            cs(id)=pgs0(id)/dens(id)/temp(id)
c           pi(id)=cprad*pradfc(id)
         end do
c
         hr1=grd(1)/qgrav
         hg1=sqrt(two*cs(1)*temp(1)/qgrav)
         x=(zd(1)-hr1)/hg1
         IF(X.LT.3.) THEN
            IF(X.LT.0.) X=0.
            F1HE=8.86226925D-1*EXP(X*X)*ERFCX(X)
          ELSE
            F1HE=HALF*(UN-HALF/X/X)/X
         END IF
c
         ntemp=2
         call pgset(ntemp)
c
         do id=1,nd
            pgs(id)=pgs0(id)
            ptotal(id)=pgs0(id)+pra(id)
        end do
c
      END IF
C
      RETURN
      END