C ======================================================================
      PROGRAM SAMPLE
C      Sample program to calculate density and sum of states of H2O
C ======================================================================
      PARAMETER (MAXVIV=20, MAXGRN=3000)
      INTEGER NVIB, NGRN, JFG(MAXVIV), IGR, IST
      DOUBLE PRECISION ENGMAX, SIZGRN, FRQNCY(MAXVIV),
     .                 RHO(MAXGRN), W(MAXGRN)
      DATA NVIB, ENGMAX, SIZGRN /3, 60000., 100./
      DATA (FRQNCY(I),I=1,3) /3657., 1595., 3756./
C
      CALL DCOUNT(NVIB, FRQNCY, SIZGRN, ENGMAX, RHO, W, NGRN, JFG, IST)
      IF (IST .NE. 0) WRITE(6, 610) IST
  610 FORMAT(' Error termination of DCOUNT with IST =', I3)
C
      WRITE(6, 620)
  620 FORMAT('   E[cm-1],     rho[cm],           W')
      DO IGR = 1, NGRN
        WRITE(6, 630) (IGR - 1) * SIZGRN, RHO(IGR), W(IGR)
  630   FORMAT(F10.2, ',', 1PE12.4, ',', 1PE12.4)
      END DO
      STOP
      END
C ----------------------------------------------------------------------
      SUBROUTINE DCOUNT(NV, FREQ, SIZG, EMAX, DENS, SUMS, NG, JFRQ, IST)
C      Direct count (Beyer-Swinehart algorithm)
C ----- input -----
C  NV     : number of vibrators
C  FREQ   : vibrational frequencies [cm-1]
C  SIZG   : grain size [cm-1]
C  EMAX   : maximum energy [cm-1]
C ----- output ----
C  DENS() : density of states (subscript 1 corresponds to E=0)
C  SUMS() : sum of states
C  NG     : number of grains
C  JFRQ() : frequencies in 'grain size' unit (rounded to integer)
C  IST    : status flag (0: normal, 1: error in freq,  2: error in emax)
C ----------------------------------------------------------------------
      INTEGER NV, NG, JFRQ(*), IST, IV, IG, JFR, NSTT, JG
      DOUBLE PRECISION FREQ(*), SIZG, EMAX, DENS(*), SUMS(*)
C     ----- Set frequencies and max energy in 'grain size' unit
      IST = 1
      IF (NV .LE. 0) RETURN
      DO IV = 1, NV
        JFRQ(IV) = IDNINT(FREQ(IV) / SIZG)
        IF (JFRQ(IV) .LE. 0) RETURN
      END DO
      NG = IDNINT(EMAX / SIZG) + 1
      IST = 2
      IF (NG .LE. 1) RETURN
      IST = 0
C     ----- Reset arrays
      DO IG = 1, NG
        DENS(IG) = 0.
        SUMS(IG) = 1.
      END DO
      DENS(1) = 1.
C     ----- Beyer-Swinehart direct count
      DO IV = 1, NV
        JFR = JFRQ(IV)
        NSTT = JFR + 1
        DO IG = NSTT, NG
          JG = IG - JFR
          DENS(IG) = DENS(IG) + DENS(JG)
          SUMS(IG) = SUMS(IG) + SUMS(JG)
        END DO
      END DO
C     ----- Convert unit of density of states from grain-1 to (cm-1)-1
      DO IG = 1, NG
        DENS(IG) = DENS(IG) / SIZG
      END DO
      RETURN
      END