SUBROUTINE MATGEN(ID)
C     =====================
C
C     Auxiliary procedure for SOLVE
C     controls evaluation of matrices A,B, and C
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'
C
C
C     evaluation of the opacity, emissivity, scattering, and
C     their derivatives, at the current depth point ID;
C
         DO I=1,NN
            PSI0(I)=PSY0(I,ID)
         END DO
         IF(NFREQE.GT.0) THEN
         DO IJ=1,NFREQE
            IJT=IJFR(IJ)
            WDEP0(IJ)=W(IJT)
            RAD0(IJ)=RADEX(IJ,ID)
            FK0(IJ)=FAKEX(IJ,ID)
            ABSO0(IJ)=ABSOEX(IJ,ID)
            EMIS0(IJ)=EMISEX(IJ,ID)
            SCAT0(IJ)=SCATEX(IJ,ID)
            DABT0(IJ)=DABTEX(IJ,ID)
            DEMT0(IJ)=DEMTEX(IJ,ID)
            DABN0(IJ)=DABNEX(IJ,ID)
            DEMN0(IJ)=DEMNEX(IJ,ID)
            DABM0(IJ)=DABMEX(IJ,ID)
            DEMM0(IJ)=DEMMEX(IJ,ID)
            DO II=1,NLVEXP
               DRCH0(II,IJ)=DRCHEX(II,IJ,ID)
               DRET0(II,IJ)=DRETEX(II,IJ,ID)
            END DO
         END DO
         END IF
C
      IF(ID.GT.1) THEN
         DO I=1,NN
            PSIM(I)=PSY0(I,ID-1)
         END DO
         IF(NFREQE.GT.0) THEN
         DO IJ=1,NFREQE
            IJT=IJFR(IJ)
            RADM(IJ)=RADEX(IJ,ID-1)
            FKM(IJ)=FAKEX(IJ,ID-1)
            ABSOM(IJ)=ABSOEX(IJ,ID-1)
            EMISM(IJ)=EMISEX(IJ,ID-1)
            SCATM(IJ)=SCATEX(IJ,ID-1)
            DABTM(IJ)=DABTEX(IJ,ID-1)
            DEMTM(IJ)=DEMTEX(IJ,ID-1)
            DABNM(IJ)=DABNEX(IJ,ID-1)
            DEMNM(IJ)=DEMNEX(IJ,ID-1)
            DABMM(IJ)=DABMEX(IJ,ID-1)
            DEMMM(IJ)=DEMMEX(IJ,ID-1)
            DO II=1,NLVEXP
               DRCHM(II,IJ)=DRCHEX(II,IJ,ID-1)
               DRETM(II,IJ)=DRETEX(II,IJ,ID-1)
            END DO
         END DO
         END IF
      END IF
C
      IF(ID.LT.ND) THEN
         DO I=1,NN
            PSIP(I)=PSY0(I,ID+1)
         END DO
         IF(NFREQE.GT.0) THEN
         DO IJ=1,NFREQE
            IJT=IJFR(IJ)
            RADP(IJ)=RADEX(IJ,ID+1)
            FKP(IJ)=FAKEX(IJ,ID+1)
            ABSOP(IJ)=ABSOEX(IJ,ID+1)
            EMISP(IJ)=EMISEX(IJ,ID+1)
            SCATP(IJ)=SCATEX(IJ,ID+1)
            DABTP(IJ)=DABTEX(IJ,ID+1)
            DEMTP(IJ)=DEMTEX(IJ,ID+1)
            DABNP(IJ)=DABNEX(IJ,ID+1)
            DEMNP(IJ)=DEMNEX(IJ,ID+1)
            DABMP(IJ)=DABMEX(IJ,ID+1)
            DEMMP(IJ)=DEMMEX(IJ,ID+1)
            DO II=1,NLVEXP
               DRCHP(II,IJ)=DRCHEX(II,IJ,ID+1)
               DRETP(II,IJ)=DRETEX(II,IJ,ID+1)
            END DO
         END DO
         END IF
      END IF
C
C
C     ------------------------------------------------------------
C     Actual evaluation of matrices A,B,C, and the rhs vector VECL
C     ------------------------------------------------------------
C
C     First null arrays A,B,C,VECL
C
      DO I=1,NN0   
         VECL(I)=0.
         DO J=1,NN0   
            B(J,I)=0.
            A(J,I)=0.
            C(J,I)=0.
            E(J,I)=0.
         END DO
      END DO
C
C     for the sake of clarity, the matrices are evaluated by several
C     different subroutines
C
      IF(IDISK.EQ.0) THEN
         CALL BRTE(ID)
         IF(INHE.NE.0) CALL BHE(ID)
         IF(INRE.NE.0) CALL BRE(ID)
       ELSE
         IF(IZSCAL.EQ.0)  THEN
            CALL BRTE(ID)
            IF(INHE.NE.0.or.inzd.ne.0) CALL BHED(ID)
            IF(INRE.NE.0) CALL BRE(ID)
          ELSE
            CALL BRTEZ(ID)
            IF(INHE.NE.0.or.inzd.ne.0) CALL BHEZ(ID)
            IF(INRE.NE.0) CALL BREZ(ID)
         END IF
      END IF
C
C     contribution to matrix elements due to convection
C
      CALL MATCON(ID)
C
      IF(INSE.GT.0) THEN
      CALL SABOLF(ID)
      CALL BPOP(ID)
      CALL EMAT(ID)
C
C     skip rows corresponding to fully-zeroed populations
C
      NSE=NFREQE+INSE
      INONZ=NSE
      DO II=NSE,NN0
         IF(IGZERT(II-NSE+1).EQ.0) THEN
          IF(INONZ.NE.II) THEN
             DO JJ=1,NN0
                B(INONZ,JJ)=B(II,JJ)
                A(INONZ,JJ)=A(II,JJ)
                C(INONZ,JJ)=C(II,JJ)
             END DO
             VECL(INONZ)=VECL(II)
          END IF
          INONZ=INONZ+1
       END IF
      END DO
C
C     skip also corresponding columns
C
      INONZ=NSE
      DO II=NSE,NN0
         IF(IGZERT(II-NSE+1).EQ.0) THEN
          IF(INONZ.NE.II) THEN
             DO JJ=1,NN
                B(JJ,INONZ)=B(JJ,II)
                A(JJ,INONZ)=A(JJ,II)
                C(JJ,INONZ)=C(JJ,II)
             END DO
          END IF
          INONZ=INONZ+1
       END IF
      END DO
      END IF
C
      RETURN
      END