SUBROUTINE RTECMC
C     =================
C
C     Solution of the radiative transfer equation with Compton scattering
C
      INCLUDE 'IMPLIC.FOR'                                                      
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ITERAT.FOR'
      COMMON/AUXRTE/
     *          COMA(MDEPTH),COMB(MDEPTH),COMC(MDEPTH),VL(MDEPTH),
     *          COME(MDEPTH),U(MDEPTH),V(MDEPTH),BS(MDEPTH),
     *          AL(MDEPTH),BE(MDEPTH),GA(MDEPTH)
      common/comgfs/gfm(mfreq,mdeptc),gfp(mfreq,mdeptc)
      DIMENSION BB(MDEPTC,MDEPTC+1),AA(MDEPTC),CC(MDEPTC),
     *          Z(MFREQ,MDEPTC),D(MFREQ,MDEPTC,MDEPTC),
     *          FF(MDEPTC,MDEPTC),ZZ(MDEPTC),
     *          drad(mfreq,mdeptc)
c
      nsti=1
      if(icomst.gt.1) nsti=icomst
      do isti=1,nsti
      DO IJ=1,NFREQ
         IJO=ijorig(ij)
         FR=FREQ(IJO)
         CALL OPACF1(IJO)
         CALL RTECF0(IJO)
         do id=1,nd
            do id1=1,nd
               bb(id,id1)=0.
            end do
         end do
         id=1
         bb(id,id)=be(id)
         bb(id,id+1)=-ga(id)
         do id=2,nd-1
            bb(id,id)=be(id)
            bb(id,id-1)=-al(id)
            bb(id,id+1)=-ga(id)
         end do
         id=nd
         bb(id,id)=be(id)
         bb(id,id-1)=-al(id)
         do id=1,nd
            if(ichcoo.eq.0) then
               bb(id,id)=bb(id,id)+un-comb(id)-bs(id)
             else
               bb(id,id)=bb(id,id)+un-comb(id)
            end if
            aa(id)=coma(id)
            cc(id)=comc(id)
         end do
c
c        linearization matrices for stimulated emission
c
         if(isti.gt.1) then
            do id=1,nd
               vl(id)=vl(id)-bb(id,id)*rad(ij,id)
               bb(id,id)=bb(id,id)-come(id)*rad(ij,id)
               aa(id)=aa(id)+u(id)*rad(ij,id) 
               cc(id)=cc(id)+v(id)*rad(ij,id) 
            end do
            id=1
            vl(id)=vl(id)-bb(id,id+1)*rad(ij,id+1)
            do id=2,nd-1
               vl(id)=vl(id)-bb(id,id-1)*rad(ij,id-1)-
     *                bb(id,id+1)*rad(ij,id+1)
            end do
            id=nd
            vl(id)=vl(id)-bb(id,id-1)*rad(ij,id-1)
            if(ij.gt.1) then
               do id=1,nd
                  vl(id)=vl(id)+aa(id)*rad(ij-1,id)
               end do
            end if
            if(ij.lt.nfreq) then
               do id=1,nd
                  vl(id)=vl(id)+cc(id)*rad(ij+1,id)
               end do
            end if
         end if
c
c        forward sweep of the grand matrix
c
         if(ij.eq.1) then
            call matinv(bb,nd,mdepth)
            do id=1,nd
               sum=0.
               do id1=1,nd
                  d(ij,id,id1)=bb(id,id1)*cc(id1)
                  sum=sum+bb(id,id1)*vl(id1)
               end do
               z(ij,id)=sum
            end do
c
          else 
            do id=1,nd
               do id1=1,nd
                  ff(id,id1)=bb(id,id1)-aa(id)*d(ij-1,id,id1)
               end do
            end do
            call matinv(ff,nd,mdepth)
            do id=1,nd
               do id1=1,nd
                  d(ij,id,id1)=ff(id,id1)*cc(id1)
               end do
            end do
            do id=1,nd
               zz(id)=vl(id)+aa(id)*z(ij-1,id)
            end do
            do id=1,nd
               sum=0.
               do id1=1,nd
                  sum=sum+ff(id,id1)*zz(id1)
               end do
               z(ij,id)=sum
            end do
         end if
      END DO
c
c     ----------------------------------
c     backward sweep of the grand matrix
c     ----------------------------------
c
      if(isti.eq.1) then
         ij=nfreq
         do id=1,nd
            rad(ij,id)=z(ij,id)
         end do
c
         DO IJ=NFREQ-1,1,-1
            do id=1,nd
                sum=0.
                do id1=1,nd
                   sum=sum+d(ij,id,id1)*rad(ij+1,id1)
                end do
                rad(ij,id)=z(ij,id)+sum
            end do
         END DO
      end if
c   
      if(isti.gt.1) then
         ij=nfreq
         do id=1,nd
            drad(ij,id)=z(ij,id)
         end do
c
         DO IJ=NFREQ-1,1,-1
            do id=1,nd
               sum=0.
               do id1=1,nd
                  sum=sum+d(ij,id,id1)*drad(ij+1,id1)
               end do
               drad(ij,id)=z(ij,id)+sum
            end do
         END DO
c
         chmax=0.
         DO IJ=1,NFREQ
            dri=0.
            do id=1,nd
               if(rad(ij,id).gt.0.) dr=drad(ij,id)/rad(ij,id)
               if(abs(dr).gt.chmax) chmax=abs(dr)
               if(abs(dr).gt.dri) dri=abs(dr)
               if(dr.gt.9.) dr=9.
               if(dr.lt.-0.999) dr=-0.999
               rad(ij,id)=rad(ij,id)*(un+dr)
            end do
         END DO
      end if
c
      if(isti.gt.1.and.chmax.lt.1.e-3) go to 100
      end do
c 
  100 continue  
      return
      end