SUBROUTINE RATMAT(ID,IICAL,IMODE,A,B)
C     =====================================
C
C     Evaluation of the rate matrix
C
C     More specifically, the set of statistical equilibrium + constraint
C     equations (see further) is represented as
C
C         A  *  vector of populations  =  B
C
C     (see Mihalas, 1978) pp.138-139 for details
C
C     Constraint equations:
C        the rows corresponding to the last explicit level of each atom
C        (ie. a one-level ion representing the highest ionization degree
C        of the given atom) are the abundance definition equations
C
C     Input: ID - depth index
C     Input transmitted by COMMON blocks:
C            RRU - array of upward radiative rates in all transitions
C            RRD - array of downward radiative rates in all transitions
C            COLRAT - array of collisional rates (only upward) in all
C                     transitions
C            SBF - array of Saha-Boltzmann factors
C            USUM - array of upper sums
C
C     Output: A - rate matrix
C             B - the right-hand-side vector
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      DIMENSION A(MLEVEL,MLEVEL),B(MLEVEL)
      DIMENSION AIJ(MTRANS),AJI(MTRANS)
      DIMENSION SBW(MLEVEL),IICAL(MLEVEL)
      DIMENSION LLTE(MLEVEL)
c     DATA ICOMP /0/
c
      if(ioptab.lt.0) return
c
      if(ipslte.ne.0) then
         do itr=1,ntrans
            rru(itr,id)=0.
            rrd(itr,id)=0.
            drdt(itr,id)=0.
         end do
      end if
c
      T=TEMP(ID)
      ANE=ELEC(ID)
      HKT=HK/T
      TK=HKT/H
      DO I=1,NLEVEL
         B(I)=0.
         SBW(I)=ANE*SBF(I)*WOP(I,ID)
         ILT=ILTION(IEL(I))
         LLT=ILT.EQ.1.AND.IMODE.EQ.0
         LLTE(I)=LLT.OR.LTE.OR.ILTLEV(I).GE.1.OR.ILT.GE.2
         LLTE(I)=LLTE(I).OR.ID.GE.IDLTE
         DO J=1,NLEVEL
            A(J,I)=0.
         END DO
      END DO
C
C     determine reference levels
C
      CALL REFLEV(ID,IABS(IMODE))
C
C ******* First part - the rows corresponding to true statistical
C                      equilibrium equations
C                      i.e. NLTE rate equations
C
C     1) simple expression in the case of LTE
C
      DO IAT=1,NATOM
         IF(IIFIX(IAT).NE.1.OR.IMODE.LT.0) THEN
            NREFI=NREFS(IAT,ID)
            DO I=N0A(IAT),NKA(IAT)
               II=IABS(IICAL(I))
               IF(I.NE.NREFI.AND.II.NE.0.AND.LLTE(I)) THEN
                  A(II,II)=UN
                  N=IABS(IICAL(ILTERF(I,ID)))
                  A(II,N)=A(II,N)-SBLPSI(I,ID)
               END IF
            END DO
         END IF
      END DO
C
C     2) NLTE rate equation
C
       IF(.NOT.LTE) THEN
          DO 40 ITR=1,NTRANS
             I=ILOW(ITR)
             IF(IIFIX(IATM(I)).EQ.1) GO TO 40
             J=IUP(ITR)
             NKE=NNEXT(IEL(I))
C
C            upward total rate
C
             AIJ(ITR)=(COLRAT(ITR,ID)+RRU(ITR,ID))*WOP(J,ID)
C
C            downward total rate
C
             IF(LINE(ITR)) THEN
C
C               bound-bound 
C
                AJI(ITR)=(COLTAR(ITR,ID)+RRD(ITR,ID)*
     *               G(I)/G(J)*EXP(HKT*FR0(ITR)))*WOP(I,ID)
             ELSE
C
C               Bound-free 
C               Quantity CORR is a factor which allows the bound-free
C               transition to end at an excited state
C
                CORR=UN
                IF(NKE.NE.J) CORR=G(NKE)/G(J)*
     *               EXP((ENION(NKE)-ENION(J))*TK)
                AJI(ITR)=COLTAR(ITR,ID)*WOP(I,ID)+
     *                   RRD(ITR,ID)*SBW(I)*CORR
             END IF
             IF(IICAL(I).LT.0.AND.IICAL(J).LT.0) THEN
                AIJ(ITR)=AIJ(ITR)*SBPSI(I,ID)
                AJI(ITR)=AJI(ITR)*SBPSI(J,ID)
             END IF
 40       CONTINUE
C
C      Elements of the rate matrix
C
         DO 41 ITR=1,NTRANS
            I=ILOW(ITR)
            IF(IIFIX(IATM(I)).EQ.1) GO TO 41
            NREFI=NREFS(IATM(I),ID)
            J=IUP(ITR)
            II=IABS(IICAL(I))
            JJ=IABS(IICAL(J))
            IF(IPZERO(I,ID).NE.0.OR.IPZERO(J,ID).NE.0) GO TO 41
            IF(I.NE.NREFI.AND.II.GT.0.AND..NOT.LLTE(I)) THEN
               A(II,II)=A(II,II)+AIJ(ITR)
               IF(JJ.GT.0) THEN
                  A(II,JJ)=A(II,JJ)-AJI(ITR)
                ELSE
                  JJJ=IICAL(ILTREF(J,ID))
                  A(II,JJJ)=A(II,JJJ)-AJI(ITR)*SBPSI(J,ID)
               END IF
            END IF
            IF(J.NE.NREFI.AND.JJ.GT.0.AND..NOT.LLTE(J)) THEN
               A(JJ,JJ)=A(JJ,JJ)+AJI(ITR)
               IF(II.GT.0) THEN
                  A(JJ,II)=A(JJ,II)-AIJ(ITR)
                ELSE
                  III=IICAL(ILTREF(I,ID))
                  A(JJ,III)=A(JJ,III)-AIJ(ITR)*SBPSI(I,ID)
               END IF
            END IF
   41    CONTINUE
      END IF
C
C        reset the rate matrix elements for "small" populations
C
         DO I=1,NLEVEL
            II=IICAL(I)
            IF(II.GT.0) THEN
               IF(IPZERO(I,ID).GT.0) THEN
                  DO J=1,NLEVEL
                     A(II,J)=0.
                 END DO
                 A(II,II)=1.
               END IF
             ELSE IF(II.LT.0) THEN
               IF(IGZERO(-II,ID).GT.0) THEN
                  DO J=1,NLEVEL
                     A(-II,J)=0.
                  END DO
                  A(-II,-II)=1.
               END IF
            END IF
         END DO
C
C ******** Second part - abundance definition equations
C
      DO 100 IAT=1,NATOM
         IF(IIFIX(IAT).EQ.1.AND.IMODE.GE.0) GO TO 100
         NREFII=IABS(IICAL(NREFS(IAT,ID)))
         DO I=N0A(IAT),NKA(IAT)
            IL=ILK(I)
            II=IICAL(I)
            IF(II.GT.0) THEN
               A(NREFII,II)=A(NREFII,II)+UN
               IF(IL.NE.0) A(NREFII,II)=A(NREFII,II)+ANE*USUM(IL)
             ELSE IF(II.LT.0) THEN
               A(NREFII,-II)=A(NREFII,-II)+SBPSI(I,ID)
               IF(IL.NE.0) A(NREFII,-II)=A(NREFII,-II)+ANE*USUM(IL)*
     *                                  SBPSI(I,ID)
             ELSE
               III=IICAL(ILTREF(I,ID))
               A(NREFII,III)=A(NREFII,III)+SBPSI(I,ID)
            END IF
         END DO
         B(NREFII)=B(NREFII)+
     *             DENS(ID)/WMM(ID)/YTOT(ID)*ABUND(IAT,ID)
  100 CONTINUE
C
      RETURN
      END