SUBROUTINE RTECOM
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/OPTDPT/DT(MDEPTH)
      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 aa(mdepth),bb(mdepth),cc(mdepth),
     *          d(mdepth),f(mdepth),z(mdepth),rd(mdepth)
c           
      PRD0=0.
      DO ID=1,ND
         PRADT(ID)=0.
      END DO
c           
c     ---------------------------------------------
c     "1st formal solution" to update Eddington factors 
c     ---------------------------------------------
c 
      if(ncfor1.gt.0) then
      do iform=1,ncfor1
         ij0=1
         if(icombc.gt.0) ij0=2
         DO IJ=ij0,NFREQ
            CALL OPACF1(IJ)
            CALL RTECF1(IJ)
         END DO
         if(icombc.gt.0) then
            ij=1
            iji=nfreq
            call rtecf0(ij)
            do id=1,nd
               rad(iji,id)=-rad(iji-1,id)*coma(id)/(comb(id)+bs(id))
            end do
         end if
      end do
      end if
c           
c     -----------------------------------------------
c     coupled solution for the frequency derivatives terms 
c     -----------------------------------------------
c
c     fully coupled treatment - traditional formulation
c
      if(ncfull.gt.0) then
      do icfull=1,ncfull
         call rtecmc 
c
c        iterative treatment of the derivative terms
c
         if(ncitot.gt.0) then
         do ictot=1,ncitot
c 
            if(nccoup.gt.0) then
            do iccoup=1,nccoup
               do ij=1,nfreq
                  ijo=ijorig(ij)
                  fr=freq(ijo)
                  call opacf1(ijo)
                  call rtecf0(ijo)
                  do id=1,nd
                     comb(id)=comb(id)+bs(id)
                     bb(id)=be(id)+un-comb(id)
                     aa(id)=al(id)
                     cc(id)=ga(id)
                     vl(id)=vl(id)+
     *                      (coma(id)*gfm(ij,id)+
     *                      comc(id)*gfp(ij,id))*rad(ij,id)
                  end do
c
c                 ----------------
c                 forward sweep
c                 ----------------
c
c                 i) upper boundary
c
                  f(1)=(bb(1)-cc(1))/cc(1)
                  d(1)=un/(un+f(1))
                  z(1)=vl(1)/bb(1)
c 
c                 ii) normal depth points
c
                  do id=2,nd-1
                     f(id)=(bb(id)-aa(id)-cc(id)+
     *                     aa(id)*f(id-1)*d(id-1))/cc(id)
                     d(id)=un/(un+f(id))
                     z(id)=(vl(id)+aa(id)*z(id-1))*d(id)/cc(id)
                  end do
c
c                 iii) lower boundary
c
                  id=nd
                  z(id)=(vl(id)+aa(id)*z(id-1))/(bb(id)-aa(id)*d(id-1))
c
c                 --------------------
c                 backward elimination 
c                 --------------------
c   
                  rd(nd)=z(nd)
                  do id=nd-1,1,-1
                     rd(id)=rd(id+1)*d(id)+z(id)
                  end do
c   
                  do id=1,nd
                     rad(ij,id)=rd(id)
                  end do
               end do
            end do
            end if
c           
c           ---------------------------------------------
c           "2nd formal solution" to update Eddington factors 
c           ---------------------------------------------
c 
            if(ncfor2.gt.0) then
            do iform=1,ncfor2
               ij0=nfreq
               if(icombc.gt.0) ij0=nfreq-1
               PRD0=0.
               DO ID=1,ND
                  PRADT(ID)=0.
               END DO
               DO IJ=1,ij0
                  ijo=ijorig(ij)
                  CALL OPACF1(IJo)
                  CALL RTECF1(IJo)
               END DO
               PRD0=PRD0*PCK
               DO ID=1,ND
                  PRADT(ID)=PRADT(ID)*PCK
               END DO
               if(icombc.gt.0) then
                  ij=1
                  iji=nfreq
                  call rtecf0(ij)
                  do id=1,nd
                     radcm(iji,id)=-radcm(iji-1,id)*coma(id)/
     *               (comb(id)+bs(id))
                  end do
                  flux(1)=radcm(iji,1)*fh(2)
               end if
c   
               do id=1,nd
                  DO IJ=1,NFREQ
                     rad(ij,id)=radcm(ij,id)
                  END DO
               end do
c
            end do
            end if
c 
         end do
         end if
c 
c     ---------------------------------------------
c     end of formal solutions
c     ---------------------------------------------
c
      end do
      end if
      return
      end