FUNCTION RHOEOS(T,P)
C     ====================
C
C     equation of state - determining density from
C     temperature and pressure 
C
C     Input:   T - temperature (K)
C              P - total pressure (cgs)
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      parameter(wmol0=1.67333E-24/2.3)
      data inirea /0/
C
C      initialize EOS tables
C
       if(inirea.eq.0) then
          call settrm
          inirea=1
       end if
C
      AN=P/BOLK/T
      RHO=AN*wmol0
C
      niteos=0
10    niteos=niteos+1
      CALL PRSENT(RHO,T,P0,S0)
      CALL PRSENT(RHO*1.01,T,P1,S1)
      DPDR=(P1-P0)/(.01*RHO)
      DRXX=(P-P0)/DPDR/rho
      if(drxx.lt.-0.9) drxx=-0.9
      rho=rho*(un+drxx)
      IF(ABS(DRXX).GT.1.d-5.and.niteos.lt.20) GO TO 10
C
      rhoeos=rho
      return
      end