SUBROUTINE STATE(ID,TE,ANE,Q)
C
C     modified LTE Saha equations - possibly using 
C     radiation temperatures after 
C     Schaerer and Schmutz AA 288, 321, 1994
C
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'WINCOM.FOR'
      common/moltst/pfmol(600,mdepth),anmol(600,mdepth),
     *              pfato(100,mdepth),anato(100,mdepth),
     *              pfion(100,mdepth),anion(100,mdepth)
      common/ioniz2/anion2(30,mdepth)
      dimension FFI(MION0)
C
      Q=0.
      DO 50 I=1,NATOMS
         IF(LGR(I)) GO TO 50
         ION=IONIZ(I)
         RQ=0.
         RS=1.
         T=TRAD(INPOT(I,1),ID)
         if(t.le.0.) t=te
         X=SQRT(T/ANE)
         XMX=2.145E4*SQRT(X)
         CALL PARTF(I,1,T,ANE,XMX,UM)
         PFSTD(1,I)=UM
         JMAX=1
         DO J=2,ION
            J1=J-1
            T=TRAD(INPOT(I,J),ID)
            if(t.le.0.) t=te
            TLN=LOG(T)*1.5
            TK=BOLK*T
            THL=11605./T
            X=SQRT(T/ANE)
            XMX=2.145E4*SQRT(X)
            DCH=EH/XMX/XMX/TK
            DCHT=DCH*J1
            FI=36.113+TLN-THL*ENEV(I,J1)+DCHT
            X=J
            XMAX=XMX*SQRT(X)
            CALL PARTF(I,J,T,ANE,XMAX,U)
            PFSTD(J,I)=U
            FI=EXP(FI)*U/UM/ANE
            FFI(J)=FI
            IF(FFI(J).GT.1.) JMAX=J
            UM=U
         END DO
         IF(JMAX.LT.ION) THEN
            R=1.
            RQ=JMAX-1
            DO J=JMAX+1,ION
               R=R*FFI(J)
               RR(I,J)=R/PFSTD(J,I)
               RS=RS+R
               RQ=RQ+(J-1)*R
            END DO
         END IF
         IF(JMAX.GT.1) THEN
            R=1.
            DO JJ=1,JMAX-1
               J=JMAX-JJ
               R=R/FFI(J+1)
               RR(I,J)=R/PFSTD(J,I)
               RS=RS+R
               RQ=RQ+(J-1)*R
            END DO
         END IF
         ABND(I)=ABNDD(I,ID)
         RR(I,JMAX)=ABND(I)/RS
         DO J=1,ION
            IF(J.NE.JMAX) RR(I,J)=RR(I,J)*RR(I,JMAX)
            if(rr(i,j).lt.1.e-35) rr(i,j)=0.
         END DO
         RR(I,JMAX)=RR(I,JMAX)/PFSTD(JMAX,I)
         X=RQ/RS
c        IF(LRM(I)) GO TO 50
         if(i.gt.1) Q=X*ABND(I)+Q
         anato(i,id)=rr(i,1)*pfstd(1,i)
         pfato(i,id)=pfstd(1,i)
         anion(i,id)=rr(i,2)*pfstd(2,i)
         pfion(i,id)=pfstd(2,i)
   50 CONTINUE
c
      do i=2,30
         anion2(i,id)=rr(i,3)*pfstd(3,i)
      end do
c
      do imol=1,500
         anmol(imol,id)=0.
         pfmol(imol,id)=0.
      end do
c
      RETURN
      END