SUBROUTINE STEQEQ(ID,POP1,MODE)
C     ===============================
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ITERAT.FOR'
      COMMON/POPSTR/POPP(MLEVEL),
     *       POPP1(MLEVEL),POPP2(MLEVEL),POPP3(MLEVEL)  
      COMMON/PPAPAR/IPOPST(MATOM),NTERST,ITERST,
     *              IACPPP,IACPP0,IACPPD,LACPPP

      DIMENSION A(MLEVEL,MLEVEL),B(MLEVEL),POP0(MLEVEL)
      DIMENSION POP1(MLEVEL),SBW(MLEVEL)
C
      if(ioptab.lt.0) return
c 
      t=temp(id)
      aein=elec(id)
      an=dens(id)/wmm(id)+elec(id)
      if(ifmol.gt.0.and.t.lt.tmolim) then
         ipri=0
         call moleq(id,t,an,aein,ane,enrg,entt,wm,ipri)
C        don't change the electron density when
C        charge conservation is not solved
C        elec(id)=ane
         if(INPC.ne.0) elec(id)=ane
      end if
c
C     evaluation of the global rate matrix
C
      CALL SABOLF(ID)
      CALL RATMAT(ID,IIFOR,1,A,B)
C
C     new populations - solution of the rate equations
C
      CALL LEVSOL(A,B,POP0,IIFOR,NLVFOR,0) 
c
C     array of new populations
C
      DO I=1,NLEVEL
         SBW(I)=ELEC(ID)*SBF(I)*WOP(I,ID)
         II=IIFOR(I)
         IF(II.GT.0) THEN
            POP1(I)=POP0(II)
          ELSE IF(II.LT.0) THEN
            POP1(I)=POP0(-II)*SBPSI(I,ID)
          ELSE
            if(imodl(i).lt.0.or.iifix(iatm(i)).gt.0) then
               pop1(i)=popul(i,id)
             else
               III=IIFOR(ILTREF(I,ID))
               POP1(I)=SBPSI(I,ID)*POP0(III) 
            end if
         END IF
         if(iifix(iatm(i)).gt.0) pop1(i)=popul(i,id)
         if(ipzero(i,id).gt.0) pop1(i)=0.
      END DO
C
C     set up the parameter IPZERO indicating that a population is "small"
C     and will subsequently be set to zero
C
      IF(ITER.EQ.0) THEN
         LKIT=.TRUE.
      ELSE
         LKIT=KANT(ITER).EQ.0 .AND. ITER.LT.IACC
      ENDIF
      IF(LKIT) THEN
      DO IAT=1,NATOM
         POPM=DENS(ID)/WMM(ID)/YTOT(ID)*ABUND(IAT,ID)
         DO I=N0A(IAT),NKA(IAT)
            IF(POP1(I)/POPM.LT.POPZER) THEN
               POP1(I)=0.
               IPZERO(I,ID)=1
            END IF
         END DO
         if(nrefs(iat,id).gt.n0a(iat)) then 
            do i=nrefs(iat,id),n0a(iat),-1
               if(ipzero(i,id).gt.0.and.ilk(i).gt.0) then
                  do iii=nfirst(ilk(i)),nlast(ilk(i))
                     ipzero(iii,id)=1
                     pop1(iii)=0.
                  end do
               end if
            end do
         end if
      END DO
      END IF
C
C     if required (MODE=1), set up the global array POPUL
C          
      IF(MODE.NE.1) RETURN
      DO I=1,NLEVEL
         POPUL(I,ID)=POP1(I)
      END DO
C
C     array of b-factors
C
      DO I=1,NLEVEL
         BFAC(I,ID)=UN      
      END DO      
      IF(LTE.OR.IPSLTE.NE.0) RETURN        
      DO ION=1,NION
         DO I=NFIRST(ION),NLAST(ION)
            IF(POPUL(NNEXT(ION),ID).GT.0..AND.IPZERO(I,ID).EQ.0)
     *      BFAC(I,ID)=POPUL(I,ID)/(POPUL(NNEXT(ION),ID)*SBW(I))
         END DO
      END DO 
C
      RETURN
      END