SUBROUTINE RTECF1(IJ)
C     =====================
C
C     Solution of the radiative transfer equation with Compton scattering
C     for one frequency (assuming the radiation intensity in i
C     other frequencies is given
C     solution is done for individual angles, and new Eddington factors 
C     are determined
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ITERAT.FOR'
      PARAMETER (SIXTH=UN/6.D0,
     *           THIRD=UN/3.D0,
     *           TWOTHR=TWO/3.D0)
      COMMON/OPTDPT/DT(MDEPTH)
      COMMON/SURFEX/EXTJ(MFREQ),EXTH(MFREQ)
      COMMON/EXTINT/WANGLE,EXTIN(MFREQ)
      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 RI(MDEPTH),RDH(MDEPTH),RDK(MDEPTH),RDN(MDEPTH),
     *          DTAU(MDEPTH),ST0(MDEPTH),RDWN(MMUC),
     *          ali(mdepth)
      DIMENSION AANU(MDEPTH),DDD(MDEPTH),FKK(MDEPTH),ali0(mdepth),
     *          SS0(MDEPTH),
     *          AAA(MDEPTH),BBB(MDEPTH),CCC(MDEPTH),EEE(MDEPTH),
     *          ZZZ(MDEPTH),ALRH(MDEPTH),ALRM(MDEPTH),ALRP(MDEPTH),
     *          ss0c(mdepth)
C
      IF(IJ.EQ.1) THEN
         if(icompt.gt.0.and.icombc.gt.0) then
            IJE=IJEX(IJ)
            DO ID=1,ND
               rad1(id)=rad(nfreq,id)
               fak1(id)=0.333333
               ali1(id)=0.
               if(ije.gt.0) then
                  RADEX(IJE,ID)=rad1(id)
                  FAKEX(IJE,ID)=fak1(id)
               END IF
            END DO
            return
         end if
      END IF
C
      WW=W(IJ)
      IJI=NFREQ-KIJ(IJ)+1
      FR=FREQ(IJ)
      CALL RTECF0(IJ)
c
      do id=1,nd
         rad1(id)=0.
         ali1(id)=0.
         rdh(id)=0.
         rdk(id)=0.
         rdn(id)=0.
         st0(id)=vl(id)+(comb(id)+bs(id))*rad(iji,id)
         ss0(id)=0.
      end do
      rdh1=0.
      rdhd=0.
c
      if(iji.gt.1) then
         do id=1,nd
            st0(id)=st0(id)+coma(id)*rad(iji-1,id)
         end do
      end if
      if(iji.lt.nfreq) then
         do id=1,nd
            st0(id)=st0(id)+comc(id)*rad(iji+1,id)
         end do
      end if
c
      if(idisk.eq.0.or.ifz0.lt.0) then
         FR15=FR*1.D-15
         BNU=BN*FR15*FR15*FR15
         PLAND=BNU/(EXP(HK*FR/TEMP(ND))-UN)*RRDIL
         DPLAN=BNU/(EXP(HK*FR/TEMP(ND-1))-UN)*RRDIL
         IF(TEMPBD.GT.0.) THEN
           PLAND=BNU/(EXP(HK*FR/TEMPBD)-UN)*RRDIL
           DPLAN=BNU/(EXP(HK*FR/TEMPBD)-UN)*RRDIL
         ENDIF
         DPLAN=(PLAND-DPLAN)/DT(ND-1)
      end if
c
      if(icomrt.eq.0) then
c
c     ========================================================
c     Formal angle-dependent solution done by Feautrier scheme
c     ========================================================
c
c     loop over angles points
c
      do i=1,nmu
         do id=1,nd-1
            dtau(id)=dt(id)/amu(i)
         end do
c
c        boundary conditions
c
         rup=0.
         rdown=0.
         rup=extint(ij,i)
         if(idisk.eq.0.or.ifz0.lt.0) rdown=pland+amu(i)*dplan
c
c        solution of the transfer equation
c
         call rtefe2(dtau,st0,rup,rdown,ri)
         ttau=0.
         do id=1,nd
            riid=wtmu(i)*ri(id)
            rad1(id)=rad1(id)+riid
            rdk(id)=rdk(id)+amu(i)*amu(i)*riid
         end do
         rdh1=rdh1+amu(i)*wtmu(i)*ri(1)
         rdhd=rdhd+amu(i)*wtmu(i)*ri(nd)
      end do
      rdh1=rdh1-half*hextrd(ij)
c
c     ----------------------
c     end of the loop over angle points
c
c     ===========================================
c     Formal angle-dependent solution done by DFE
c     ===========================================
c
      else
c
c     loop over angle points
c     ----------------------
c
      do i=1,nmuc
         do id=1,nd-1
            dtau(id)=dt(id)/abs(amuc(i))
         end do
c
c        boundary conditions
c
         rup=0.
         rdown=0.
         if(amuc(i).lt.0.) rup=extint(ij,i)
C
C        diffusion approximation for semi-infinite atmospheres
C
         if(idisk.eq.0.or.ifz0.lt.0) rdown=pland+amuc(i)*dplan
c
c        the case of finite slab - irradiation of the back side
c
         if(amuc(i).gt.0.) rdown=rdwn(nmuc-i+1)
c
c        solution of the transfer equation
c
         call rtesol(dtau,st0,rup,rdown,amuc(i),ri,ali)
         ttau=0.
         do id=1,nd
            riid=ri(id)*half
            rad1(id)=rad1(id)+wtmuc(i)*riid
            ali1(id)=ali1(id)+wtmuc(i)*ali(id)
            rdh(id)=rdh(id)+amuc1(i)*riid
            rdk(id)=rdk(id)+amuc2(i)*riid
            rdn(id)=rdn(id)+amuc3(i)*riid
         end do
         rdwn(i)=ri(nd)
         if(amuc(i).gt.0.) rdh1=rdh1+amuc1(i)*ri(1)*half
         rdhd=rdhd+abs(amuc1(i))*ri(nd)*half
      end do
c
c     ----------------------
c     end of the loop over angle points
c
      end if
c
      do id=1,nd
         fak1(id)=fak(ij,id)
         radk(ij,id)=rdk(id)
         if(icomve .gt. 0) then
            fkk(id)=rdk(id)/rad1(id)
          else
            fkk(id)=fak(ij,id)
         endif
         ss0(id)=0.
      end do
      if(icomve.gt.0) then
         do id=1,nd
            fak(ij,id)=rdk(id)/rad1(id)
            fak1(id)=fak(ij,id)
            fkk(id)=fak(ij,id)
         end do
      end if
      if(rad1(1).gt.0.) then
         flux(ij)=rdh1
         fhd(ij)=rdhd/rad1(nd)
      end if
c
      ah=rdh1
      if(iwinbl.lt.0) ah=ah+half*hextrd(ij)
      aj=rad1(1)
      fh(ij)=ah/aj
C
C     ********************
C
C     Again solution of the transfer equation, now with Eddington
C     FKK and FH determined above, to insure strict consistency of the
C     radiation field and Eddington factors
C
C     Upper boundary condition
C
      U0=0.
      QQ0=0.
      US0=0.
      TAUMIN=ABSO1(1)*DEDM1
      NMU=3
      DO I=1,NMU
         IF(IWINBL.EQ.0.AND.WANGLE.EQ.0.) THEN
C
C           allowance for non-zero optical depth at the first depth point
C
            TAMM=TAUMIN/AMU(I)
            EX=EXP(-TAMM)
            P0=UN-EX
            QQ0=QQ0+P0*AMU(I)*WTMU(I)
            U0=U0+EX*WTMU(I)
            if(tamm.gt.0.) US0=US0+P0/TAMM*WTMU(I)
         END IF
      END DO
      ID=1
      DTP1=DT(ID)
      IF(MOD(ISPLIN,3).EQ.0) THEN
         B=DTP1*HALF
         C=0.
       ELSE
         B=DTP1*THIRD
         C=B*HALF
      END IF
      BQ=UN/(B+QQ0)
      CQ=C*BQ
      BBB(ID)=(FKK(ID)/DTP1+FH(IJ)+B)*BQ+SS0(ID)
      CCC(ID)=(FKK(ID+1)/DTP1)*BQ-CQ*(UN+SS0(ID+1))
      ZZZ(ID)=UN/BBB(ID)
      VLL=ST0(ID)+CQ*ST0(ID+1)
c     IF(IWINBL.LT.0) VLL=VLL+TWO*HEXTRD(IJ)/DTP1
      AANU(ID)=VLL*ZZZ(ID)
      DDD(ID)=CCC(ID)*ZZZ(ID)
      IF(ISPLIN.GT.2) FFF=BBB(ID)/CCC(ID)-UN
C
C     Normal depth point
C
      DO ID=2,ND-1
         DTM1=DTP1
         DTP1=DT(ID)
         DT0=TWO/(DTP1+DTM1)
         ALP=UN/DTM1*DT0
         GAM=UN/DTP1*DT0
         IF(MOD(ISPLIN,3).EQ.0) THEN
            A=0.
            C=0.
          ELSE IF(ISPLIN.EQ.1) THEN
            A=DTM1*DT0*SIXTH
            C=DTP1*DT0*SIXTH
          ELSE
            A=(UN-HALF*DTP1*DTP1*ALP)*SIXTH
            C=(UN-HALF*DTM1*DTM1*GAM)*SIXTH
         END IF
         AAA(ID)=ALP*FKK(ID-1)-A*(UN+SS0(ID-1))
         CCC(ID)=GAM*FKK(ID+1)-C*(UN+SS0(ID+1))
         BBB(ID)=(ALP+GAM)*FKK(ID)+(UN-A-C)*(UN+SS0(ID))
         VLL=A*ST0(ID-1)+C*ST0(ID+1)+(UN-A-C)*ST0(ID)
         AANU(ID)=VLL+AAA(ID)*AANU(ID-1)
         IF(ISPLIN.LE.2) THEN
            ZZZ(ID)=UN/(BBB(ID)-AAA(ID)*DDD(ID-1))
            DDD(ID)=CCC(ID)*ZZZ(ID)
            AANU(ID)=AANU(ID)*ZZZ(ID)
          ELSE
            SUM=-AAA(ID)+BBB(ID)-CCC(ID)
            FFF=(SUM+AAA(ID)*FFF*DDD(ID-1))/CCC(ID)
            DDD(ID)=UN/(UN+FFF)
            AANU(ID)=AANU(ID)*DDD(ID)/CCC(ID)
         ENDIF
      END DO  
C
C     Lower boundary condition
C
      ID=ND
c
c     stellar atmospheric
c
      IF(IDISK.EQ.0.OR.IFZ0.LT.0) then
      IF(IBC.EQ.0) THEN
         BBB(ID)=FKK(ID)/DTP1+HALF
         AAA(ID)=FKK(ID-1)/DTP1
         VLL=HALF*PLAND+THIRD*DPLAN
       ELSE IF(IBC.LT.4) THEN
         B=UN/DTP1
         A=TWO*B*B
         BBB(ID)=UN+SS0(ID)+B*TWO*FHD(IJ)+A*FKK(ID)
         AAA(ID)=A*FKK(ID-1)
         VLL=ST0(ID)+B*(PLAND+TWOTHR*DPLAN)
       ELSE 
         B=UN/DTP1
         A=TWO*B*B
         BBB(ID)=B+A*FKK(ID)
         AAA(ID)=A*FKK(ID-1)
         VLL=B*(PLAND+TWOTHR*DPLAN)
      END IF
c
c     accretion disk - symmetric boundary
c
      ELSE
         B=TWO/DTP1
         BBB(ID)=FKK(ID)/DTP1*B+UN+SS0(ND)
         AAA(ID)=FKK(ID-1)/DTP1*B
         VLL=ST0(ID)
      END IF
C
      EEE(ND)=AAA(ID)/BBB(ID)
      ZZZ(ID)=UN/(BBB(ID)-AAA(ID)*DDD(ID-1))
      RAD1(ID)=(VLL+AAA(ID)*AANU(ID-1))*ZZZ(ID)
      FAK1(ID)=FKK(ND)
      ALRH(ID)=ZZZ(ID)
      frd=bbb(nd)*rad1(nd)-aaa(nd)*rad1(nd-1)
      frd1=(bbb(nd)-un)*rad1(nd)-aaa(nd)*rad1(nd-1)
C
C     Backsolution
C
      DO ID=ND-1,1,-1
         EEE(ID)=AAA(ID)/(BBB(ID)-CCC(ID)*EEE(ID+1))
         RAD1(ID)=AANU(ID)+DDD(ID)*RAD1(ID+1)
         FAK1(ID)=FKK(ID)
C        write(42,642),ij,id,rad1(id),st0(id),fak1(id)
         ALRH(ID)=ZZZ(ID)/(UN-DDD(ID)*EEE(ID+1))
         ALRM(ID)=0
         ALRP(ID)=0
      END DO
c
C        evaluate approximate Lambda operator
C
C        a) Rybicki-Hummer Lambda^star operator (diagonal)
C           (for JALI = 1)
C
      DO ID=1,ND
         ALIM1(ID)=0.
         ALIP1(ID)=0.
      END DO
      IF(JALI.EQ.1) THEN
         DO ID=1,ND 
            ALI1(ID)=ALRH(ID)
         END DO
c
         IF(IBC.EQ.0) THEN
            ali1(nd-1)=rad1(nd-1)/st0(nd-1)
            ali1(nd)=rad1(nd)/st0(nd)
         END IF
C
C        for IFALI>5:
C        tridiagonal Rybicki-Hummer operator (off-diagonal terms)
C
         IF(IFALI.GE.6) THEN
            ALIP1(1)=ALRH(2)*DDD(1)
            DO ID=2,ND-1 
               ALIM1(ID)=ALRH(ID-1)*EEE(ID)
               ALIP1(ID)=ALRH(ID+1)*DDD(ID)
            END DO
            ALIM1(ND)=ALRH(ND-1)*EEE(ND)
            IF(IBC.EQ.0) THEN
               ALIM1(nd)=0.
               ALIM1(nd-1)=0.
               ALIP1(nd)=0.
               ALIP1(nd-1)=0.
            END IF
         END IF
c
C        b) diagonal Olson-Kunasz Lambda^star operator, 
C           (for JALI = 2)
C
       ELSE IF(JALI.EQ.2) THEN
         DO ID=1,ND-1
            ALI0(ID)=0.
            DO I=1,NMU
               DIV=DT(ID)/AMU(I)
               ALI0(ID)=ALI0(ID)+(UN-EXP(-DIV))/DIV*WTMU(I)
            END DO
         END DO
         DO ID=2,ND-1
            ALI1(ID)=UN-HALF*(ALI0(ID)+ALI0(ID-1))
         END DO
         ALI1(1)=UN-HALF*(ALI0(1)+US0)
         ALI1(ND)=UN-ALI0(ND-1)
         ali1(nd-1)=rad1(nd-1)/st0(nd-1)
         ali1(nd)=rad1(nd)/st0(nd)
      END IF
C          
C      correction of Lambda^star for scattering
C  
      IF(ILMCOR.EQ.1) THEN     
         DO ID=1,ND
             ALI1(ID)=ALI1(ID)*(UN+SS0(ID))
             ALIM1(ID)=ALIM1(ID)*(UN+SS0(ID))
             ALIP1(ID)=ALIP1(ID)*(UN+SS0(ID))
         END DO   
       ELSE IF(ILMCOR.EQ.3) THEN     
         DO ID=1,ND
             ALI1(ID)=ALI1(ID)/(UN+SS0C(ID)*ALI1(ID))
             ALIM1(ID)=ALIM1(ID)/(UN+SS0C(ID)*ALIM1(ID))
             ALIP1(ID)=ALIP1(ID)/(UN+SS0C(ID)*ALIP1(ID))
         END DO   
      END IF
C
      DO ID=1,ND
         radcm(iji,id)=rad1(id)
      END DO
C
C     radiation pressure
C
      if(.not.lskip(1,IJ))
     *   PRD0=PRD0+ABSO1(1)*WW*(RAD1(1)*FH(IJ)-HEXTRD(IJ))
      DO ID=1,ND
         if(.not.lskip(ID,IJ))
     *   PRADT(ID)=PRADT(ID)+RAD1(ID)*FAK1(ID)*WW
         PRADA(ID)=PRADA(ID)+RAD1(ID)*FAK1(ID)*WW
      END DO
c
      if(chmax.ge.1.91e-3.and.chmax.le.2.03e-3) then
         tauij=taumin
         do id=1,nd
         if(id.gt.1) tauij=tauij+dt(id-1)
            write(97,697) ij,id,tauij,rad1(id),st0(id)/(un+ss0(id)),
     *      st0(id),un+ss0(id),ali1(id)
         end do
  697    format(2i4,1p6e12.4)
      end if
c
      do id=1,nd
         fak(ij,id)=fak1(id)
      end do
C
C     store quantities for explicit (linearized) frequencies
C
      IF(IJEX(IJ).LE.0) RETURN
      IJE=IJEX(IJ)
      DO ID=1,ND
         RADEX(IJE,ID)=RAD1(ID)
         FAKEX(IJE,ID)=FAK1(ID)
      END DO
c
      RETURN
      END