SUBROUTINE GREYD
C     ================
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'ALIPAR.FOR'
      PARAMETER (ERRM0=1.E-3, NTRM=50)
C
      CHI0=20.
      IF(TEFF.GT.10000.) THEN
         XION=2.
       ELSE IF(TEFF.LT.6000.) THEN
         XION=1.
       ELSE
         XION=1.+(TEFF-6000.)/4000.
      END IF
      ANEREL=1.-1./XION
C
      DMP=0.
      ID=1
      C1=BOLK*XION/WMM(ID)
      C2=3.1415926/2./QGRAV
      C3=SQRT(C1*C2)
      C4=WBARM*OMEG32/ALPHAV
      C5=C4/C1
      C6=C1/WBARM/OMEG32
C
      ITRM=0
   10 CONTINUE
      ITRM=ITRM+1
      C1=BOLK*XION/WMM(ID)
      C3=SQRT(C1*C2)
      C5=C4/C1
      T=(0.375*TEFF**4*CHI0*C5)**0.2
      DM0=C5/T
      RHO=DM0/SQRT(T)/C3
      TEMP(ID)=T
      DENS(ID)=RHO
      CALL RHONEN(ID,T,RHO,AN,ANE)
      ELEC(ID)=ANE
      XION=AN/(AN-ANE)
      CALL WNSTOR(ID)
      CALL STEQEQ(ID,POP,1)
C
C     evaluation of the Rosseland and Planck mean opacities
C     for the new values of temperature, electron density, and
C     populations (OPROS - Rosseland opacity per 1 cm**3; OPPLA - Planck
C     mean opacity per 1 cm**3)
C
      CALL OPACF0(ID,NFREQ)
      CALL MEANOP(T,ABSO,SCAT,OPROS,OPPLA)
      ABROS=OPROS/DENS(ID)
      ABPLA=OPPLA/DENS(ID)
      if(abpla.lt.abpmin) abpla=abpmin
      CHI0=(ABROS+chi0)/2.
      WRITE(6,601) ITRM,T,DM0,RHO,CHI0,abros,CHI0*DM0,XION,ANE
  601 FORMAT(i5,1p8e9.2)
      ERRM=ABS(DM0-DMP)/DM0
      DMP=DM0
      IF(ERRM.GT.ERRM0.AND.ITRM.LT.NTRM) GO TO 10
C
      DMTOT=DM0
      RETURN
      END