SUBROUTINE BREZ(ID)
C     ===================
C
C     The part of matrices A and B corresponding to the radiative
C     equilibrium equation
C     i.e. the (NFREQE+INRE)-th row
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'
      DIMENSION REXB(MLEVEL)
      EQUIVALENCE (REX(1),REXB(1))
C
      NRE=NFREQE+INRE
      NHE=NFREQE+INHE
      NPC=NFREQE+INPC
      NMP=NFREQE+INMP
      NSE=NFREQE+INSE-1
      IJ1=1
      if(icompt.gt.0.and.icombc.gt.0.and.ijex(1).gt.0) IJ1=2
C
      ittc=abs(nretc)/100
      if(iter.gt.ittc) then
      if(id.le.mod(abs(nretc),100)) then
         b(nre,nre)=1.
       if(nretc.lt.0) then 
          c(nre,nre)=-1.
          vecl(nre)=temp(id+1)-temp(id)
         end if
       return
      end if
      end if
C
C     the rhs vector accounts for total net cooling in ALI
C     transitions (FCOOL)
C 
      VECL(NRE)=FCOOL(ID)-reint(id)*TVISC(ID)
      if(reint(id).le.0) go to 100
C
C *********     integral equation part of the radiative
C     equilibrium equation
C
      BREPC=0.
      BREMP=0.
      DO I=1,NLVEXP
         REXB(I)=0.
      END DO
      IF(NFREQE.GT.0) THEN
      DO IJ=IJ1,NFREQE
         IJT=IJFR(IJ)
         BREPC=BREPC+((DABN0(IJ)-SIGEC(IJT))*RAD0(IJ)-
     *         DEMN0(IJ))*WDEP0(IJ)
         BREMP=BREMP+(DABM0(IJ)*RAD0(IJ)-DEMM0(IJ))*WDEP0(IJ)
         DO I=1,NLVEXP
            REXB(I)=REXB(I)+(DRCH0(I,IJ)*RAD0(IJ)-
     *             DRET0(I,IJ))*WDEP0(IJ)
         END DO
         B(NRE,NRE)=B(NRE,NRE)+(DABT0(IJ)*RAD0(IJ)-
     *              DEMT0(IJ))*WDEP0(IJ)*reint(id)
         HEAT=ABSO0(IJ)-SCAT0(IJ)
         B(NRE,IJ)=WDEP0(IJ)*HEAT*reint(id)
         VECL(NRE)=VECL(NRE)-(HEAT*RAD0(IJ)-EMIS0(IJ))*WDEP0(IJ)*
     *             reint(id)
c
c        additional terms for Compton scattering
c
         if(icompt.gt.5) then
          ijt=ijfr(ij)
          call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
          vecl(nre)=vecl(nre)+abso0(ij)*cms*wdep0(ij)*reint(id)
          if(icompt.gt.6) then
          if(icmdra.gt.0) then
            b(nre,ij)=b(nre,ij)-abso0(ij)*(cmb+cme)*wdep0(ij)*reint(id)
          else
          b(nre,ij)=b(nre,ij)-abso0(ij)*(cmb+cme)*reint(id)
          end if
            iji=nfreq-kij(ijt)+1
          if(iji.gt.1) then
             ijm=ijex(ijorig(iji-1))
             if(ijm.gt.0) then
             if(icmdra.gt.0) then
             b(nre,ijm)=b(nre,ijm)-abso0(ij)*cma*wdep0(ij)*reint(id)
             else
             b(nre,ijm)=b(nre,ijm)-abso0(ij)*cma*reint(id)
             end if
             end if
          end if
          if(iji.lt.nfreq) then
             ijp=ijex(ijorig(iji+1))
             if(ijp.gt.0) then
             if(icmdra.gt.0) then
             b(nre,ijp)=b(nre,ijp)-abso0(ij)*cmc*wdep0(ij)*reint(id)
             else
             b(nre,ijp)=b(nre,ijp)-abso0(ij)*cmc*reint(id)
             end if
             end if
          end if
          b(nre,nre)=b(nre,nre)-cmd*abso0(ij)*wdep0(ij)*reint(id)
          b(nre,npc)=b(nre,npc)-cms*abso0(ij)/elec(id)*wdep0(ij)*
     *                 reint(id)
          end if
       end if
C     
      END DO
      END IF
C
C     corrections for ALI frequency points
C
      B(NRE,NRE)=B(NRE,NRE)+REIT(ID)*reint(id)
      IF(INPC.GT.0) B(NRE,NPC)=B(NRE,NPC)+(BREPC+REIN(ID))*reint(id)
      IF(INPC.GT.0) B(NRE,NPC)=B(NRE,NPC)+(BREPC+REIN(ID))*reint(id)
      IF(INMP.GT.0) B(NRE,NMP)=B(NRE,NMP)+(BREMP+REIM(ID))*reint(id)
      IF(INHE.GT.0) B(NRE,NHE)=REIX(ID)*reint(id)
      A(NRE,NRE)=AREIT(ID)*reint(id)
      IF(INPC.GT.0) A(NRE,NPC)=AREIN(ID)*reint(id)
      C(NRE,NRE)=CREIT(ID)*reint(id)
      IF(INPC.GT.0) C(NRE,NPC)=CREIN(ID)*reint(id)
      IF(INMP.GT.0) C(NRE,NMP)=CREIM(ID)*reint(id)
      IF(INHE.GT.0) C(NRE,NHE)=CREIX(ID)*reint(id)
c     END IF
C
C     terms arising because of viscosity
C
      B(NRE,NRE)=B(NRE,NRE)+DTVIST(ID)*reint(id)
      IF(INPC.GT.0) B(NRE,NPC)=B(NRE,NPC)-DTVISR(ID)*reint(id)
      IF(INHE.GT.0) B(NRE,NFREQE+INHE)=
     *               (DTVISR(ID)+DTVISN(ID))*reint(id)
      IF(INMP.GT.0) B(NRE,NFREQE+INMP)=
     *                DTVISR(ID)*HMASS/WMM(ID)*reint(id)
C      
      DO II=1,NLVEXP
         B(NRE,NSE+II)=B(NRE,NSE+II)+(REXB(II)+REIP(II,ID))*reint(id)
      END DO
      IF(IFALI.GT.5.AND.ID.GT.1) THEN
      DO II=1,NLVEXP
         A(NRE,NSE+II)=A(NRE,NSE+II)+AREIP(II,ID)*reint(id)
      END DO
      END IF
      IF(IFALI.GT.5.AND.ID.LT.ND) THEN
      DO II=1,NLVEXP
         C(NRE,NSE+II)=C(NRE,NSE+II)+CREIP(II,ID)*reint(id)
      END DO
      END IF
C
C *********    differential equation part of the
C     radiative equilibrium equation
C
  100 CONTINUE
      if(redif(id).eq.0) return
C 
      TEFFD=TEFF**4*(UN-THETAV(ID))
      VECL(NRE)=VECL(NRE)+SIG4P*TEFFD*redif(id)
      if(id.eq.1) go to 200
C
      DDM=(ZD(ID-1)-ZD(ID))*HALF 
      AREN=0.
      BREN=0.
      AREPC=0.
      BREPC=0.
C
      GP=0.
      GN=UN
      IF(INMP.GT.0) THEN
         GP=UN
         GN=0.
      END IF
C
      DO I=1,NLVEXP
         REXB(I)=0.
         REXA(I)=0.
      END DO
C
      IF(NFREQE.GT.0) THEN
      DO IJ=1,NFREQE
         OMEG0=ABSO0(IJ)
         OMEGM=ABSOM(IJ)
         DTAUM=(OMEG0+OMEGM)*DDM
         FRD=FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ)
         GAMR=FRD/DTAUM
         A1=GAMR/(OMEG0+OMEGM)*WDEP0(IJ)
C
C     Corresponding elements of matrix A
C
         A(NRE,IJ)=-WDEP0(IJ)*FKM(IJ)/DTAUM*redif(id)
         AREPC=AREPC-A1*DABNM(IJ)
         A(NRE,NRE)=A(NRE,NRE)-A1*DABTM(IJ)*redif(id)
C
C     Corresponding elements of matrix B
C     Columns corresponding to mean intensities
C
         B(NRE,IJ)=B(NRE,IJ)+WDEP0(IJ)*FK0(IJ)/DTAUM*redif(id)
         BREPC=BREPC-A1*DABN0(IJ)
C
C     Column corresponding to temperature
C
         B(NRE,NRE)=B(NRE,NRE)-A1*DABT0(IJ)*redif(id)
C
C     auxiliary vectors for columns corresponding to populations
C
         DO I=1,NLVEXP
            REXA(I)=REXA(I)-A1*DRCHM(I,IJ)
            REXB(I)=REXB(I)-A1*DRCH0(I,IJ)
         END DO
C
C     The rhs vector
C
         VECL(NRE)=VECL(NRE)-WDEP0(IJ)*GAMR*redif(id)
      END DO
      END IF
C
C     Column corresponding to temperature
C
      A(NRE,NRE)=A(NRE,NRE)+REDTM(ID)*REDIF(ID)
      B(NRE,NRE)=B(NRE,NRE)+REDT(ID)*REDIF(ID)
      C(NRE,NRE)=C(NRE,NRE)+REDTP(ID)*REDIF(ID)
C
C     Column corresponding to electron density (for matrices A and B)
C
      IF(INPC.NE.0) THEN
         A(NRE,NPC)=A(NRE,NPC)+(AREPC+REDNM(ID))*redif(id)
         B(NRE,NPC)=B(NRE,NPC)+(BREPC+REDN(ID))*redif(id)
         C(NRE,NPC)=C(NRE,NPC)+REDNP(ID)*redif(id)
      END IF
C
C     Columns corresponding to populations (for matrices A and B)
C     Note that auxiliary arrays REXA and REX, which contain derivatives
C     of the absorption and emission coefficients wrt populations, have
C     been generated by BRTE
C
      DO II=1,NLVEXP
         A(NRE,NSE+II)=A(NRE,NSE+II)+(REXA(II)+REDPM(II,ID))*redif(id)
         B(NRE,NSE+II)=B(NRE,NSE+II)+(REXB(II)+REDP(II,ID))*redif(id)
         C(NRE,NSE+II)=C(NRE,NSE+II)+REDPP(II,ID)*redif(id)
      END DO
      RETURN
c     
C     upper boundary condition for the differential form
c
  200 CONTINUE
C
C     Columns corresponding to mean intensities; rhs vector
C
      IF(NFREQE.GT.0) THEN
      DO IJ=1,NFREQE
         IJT=IJFR(IJ)
         WF=WDEP0(IJ)*FH(IJT)*REDIF(ID)
         B(NRE,IJ)=B(NRE,IJ)+WF
         VECL(NRE)=VECL(NRE)-WF*RAD0(IJ)
      END DO
      END IF
C
C     Column corresponding to temperature
C
      B(NRE,NRE)=B(NRE,NRE)+REDT(ID)*REDIF(ID)
C
C     Column corresponding to electron density 
C
      IF(INPC.NE.0) THEN
         B(NRE,NPC)=B(NRE,NPC)+REDN(ID)*redif(id)
      END IF
C
C     Columns corresponding to populations
C
      DO II=1,NLVEXP
         B(NRE,NSE+II)=B(NRE,NSE+II)+REDP(II,ID)*redif(id)
      END DO
      RETURN
      END