SUBROUTINE BPOP(ID)
C     ===================
C
C     The part of matrix B corresponding to the statistical
C     equilibrium equations
C     i.e. the (NFREQE+INSE)-th thru (NFREQE+INSE+NLVEXP-1)-th rows;
C     and to the charge conservation equation, ie. the (NFREQE+INPC)-th
C     row
C
C     The formalism is similar to that described in Mihalas, Stellar
C     Atmospheres, 1978, pp. 143-145
C
C     Input: ID - depth index
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ARRAY1.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ODFPAR.FOR'
      INCLUDE 'ITERAT.FOR'
      dimension sbw(mlevel)
      dimension popp(mlevel)
C
      if(ioptab.lt.0) return
      DO I=1,NLVEXP
         ATT(I)=0.
         ANN(I)=0.
      END DO
C
      IF(.NOT.LTE .AND. IFPOPR.EQ.5.AND.IPSLTE.EQ.0) THEN
      CALL RATMAT(ID,IIFOR,0,ESEMAT,BESE)
      CALL LEVSOL(ESEMAT,BESE,POPP,IIFOR,NLVFOR,0)   
      DO I=1,NLEVEL
         II=IIEXP(I)
         IF(II.EQ.0.AND.IMODL(I).EQ.6) THEN
            III=ILTREF(I,ID)
            SBPSI(I,ID)=POPP(I)/POPP(III)
         END IF
      END DO
      DO ION=1,NION
         DO I=NFIRST(ION),NLAST(ION)
             SBW(I)=ELEC(ID)*SBF(I)*WOP(I,ID)
             IF(POPUL(NNEXT(ION),ID).GT.0..AND.IPZERO(I,ID).EQ.0)
     *       BFAC(I,ID)=POPP(I)/(POPP(NNEXT(ION))*SBW(I))
         END DO
      END DO
      END IF
C
      CALL LEVGRP(ID,IIEXP,0,POPP)
      CALL RATMAT(ID,IIEXP,0,ESEMAT,BESE)
C
      IF(IFPOPR.LE.3) CALL MATINV(ESEMAT,NLVEXP,MLEVEL)
C
C     Split BPOP in separate subroutines
C
      if(ipslte.eq.0) then
      IF(.NOT.LTE.AND.IBPOPE.GT.0.AND.ID.LT.IDLTE) THEN
         CALL BPOPE(ID)
         CALL BPOPF(ID)
      END IF
      end if
      CALL BPOPT(ID)
      IF(INPC.GT.0) CALL BPOPC(ID)
C
C     reset matrix elements for "small" populations
C
      DO I=1,NLVEXP
         IF(IGZERO(I,ID).NE.0) THEN
         DO J=1,NLVEXP
            B(NFREQE+INSE-1+I,NFREQE+INSE-1+J)=0.
         END DO
         B(NFREQE+INSE-1+I,NFREQE+INSE-1+I)=1.
         VECL(NFREQE+INSE-1+I)=0.
         END IF
      END DO
      RETURN
      END