SUBROUTINE RTEDF1(IJ)
C     =====================
C
C     Solution of the radiative transfer equation - for one frequency 
C     for the known source function.
C     Determination of the radiation field and variable Eddington
C     factors.
C
C     The numerical method used:
c     Discontinuous Finite Element method
c     Castor, Dykema, Klein, 1992, ApJ 387, 561.
C
C     different formulation of the boundary conditions
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ALIPAR.FOR'
      COMMON/OPTDPT/DT(MDEPTH)
      PARAMETER (SIXTH=UN/6.D0,
     *           THIRD=UN/3.D0,
     *           TWOTHR=TWO/3.D0)
      DIMENSION RDK(MDEPTH),FKK(MDEPTH),
     *          ST0(MDEPTH),SA0(MDEPTH),SS0(MDEPTH),
     *          dtau(mdepth),rip(mdepth),rim(mdepth),
     *          riin(mdepth),riup(mdepth),u(mdepth),
     *          aip(mdepth),aim(mdepth),al0(mdepth),
     *          aiin(mdepth),aiup(mdepth),
     *          ali0(mdepth),ss0c(mdepth),
     *          AAA(MDEPTH),BBB(MDEPTH),CCC(MDEPTH),EEE(MDEPTH),
     *          ZZZ(MDEPTH),ALRH(MDEPTH),ALRM(MDEPTH),ALRP(MDEPTH),
     *          DDD(MDEPTH),AANU(MDEPTH)
C
      FR=FREQ(IJ)
C
C     optical depth scale
C
      DO ID=1,ND-1
         DT(ID)=DELDMZ(ID)*(ABSOT(ID+1)+ABSOT(ID))
         SA0(ID)=EMIS1(ID)/ABSO1(ID)
         SS0(ID)=-SCAT1(ID)/ABSO1(ID)
      END DO
      SA0(ND)=EMIS1(ND)/ABSO1(ND)
      SS0(ND)=-SCAT1(ND)/ABSO1(ND)
C
      TAUMIN=ABSO1(1)*DEDM1
C
C     Allowance for wind blanketing
C
      ALB1=0.
      IF(IWINBL.GT.0) ALB1=TWO*ALBE(IJ)/(UN+ALBE(IJ))
C
C     quantities for the lower boundary condition
C
      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)
C
C     global ALI loop for treating electron scattering
C
      itrali=0
   10 itrali=itrali+1
C
C     total source function
C  
      DO ID=1,ND
         ST0(ID)=SA0(ID)-SS0(ID)*RAD(IJ,ID)
         RAD1(ID)=0.
         RDK(ID)=0.
         ALI1(ID)=0.
      END DO
      AH=0.
      ahout=0.
      ahd=0.
      U0=0.
      QQ0=0.
      US0=0.
c
c     loop over angle poits
c
      DO I=1,NMU
         AMU2=AMU(I)*AMU(I)*WTMU(I)
         do id=1,nd-1
            dtau(id)=dt(id)/amu(i)
         enddo
c               
c        incoming intensity 
c   
         ID=1
         P0=0.
         EX=UN
C
C        allowance for non-zero optical depth at the first depth point
C
c        rim(id)=EXTRAD(IJ)
         rim(id)=EXTINT(IJ,I)
c         IF(IWINBL.EQ.0) THEN
c            TAMM=TAUMIN/AMU(I)
c            EX=EXP(-TAMM)
c            P0=UN-EX
c            QQ0=QQ0+P0*AMU(I)*WTMU(I)
c            U0=U0+EX*WTMU(I)
c            US0=US0+P0/TAMM*WTMU(I)
c            rim(id)=st0(id)*p0
c         END IF
c               
c          incoming intensity 
c   
            aim(id)=0.
            do id=1,nd-1
               dt0=dtau(id)
               dtaup1=dt0+un
               dtau2=dt0*dt0
               bb=two*dtaup1
               cc=dt0*dtaup1
               aa=un/(dtau2+bb)
               rim(id+1)=(two*rim(id)+dt0*st0(id)+cc*st0(id+1))*aa
               rip(id)=(bb*rim(id)+cc*st0(id)-dt0*st0(id+1))*aa
               aim(id+1)=cc*aa
               aip(id)=(cc+bb*aim(id))*aa
            end do
            do id=2,nd-1
               dtt=un/(dtau(id-1)+dtau(id))
               riin(id)=(rim(id)*dtau(id)+rip(id)*dtau(id-1))*dtt
               aiin(id)=(aim(id)*dtau(id)+aip(id)*dtau(id-1))*dtt
            end do
            riin(1)=rim(1)
            riin(nd)=rim(nd)
            aiin(1)=aim(1)
            aiin(nd)=aim(nd)
C
c           outgoing intensity
c 
            if(idisk.eq.0) rim(nd)=PLAND+AMU(I)*DPLAN
            do id=nd-1,1,-1
               dt0=dtau(id)
               dtaup1=dt0+un
               dtau2=dt0*dt0
               bb=two*dtaup1
               cc=dt0*dtaup1
               aa=un/(dtau2+bb)
               rim(id)=(two*rim(id+1)+dt0*st0(id+1)+cc*st0(id))*aa
               rip(id+1)=(bb*rim(id+1)+cc*st0(id+1)-dt0*st0(id))*aa
               aim(id)=cc*aa
               aip(id+1)=(cc+bb*aim(id+1))*aa
            end do
            do id=2,nd-1
               dtt=un/(dtau(id-1)+dtau(id))
               riup(id)=(rim(id)*dtau(id-1)+rip(id)*dtau(id))*dtt
               aiup(id)=(aim(id)*dtau(id-1)+aip(id)*dtau(id))*dtt
            end do
            riup(1)=rim(1)
            riup(nd)=rim(nd)
            aiup(1)=aim(1)
            aiup(nd)=aim(nd)
c
c           final symmetrized (Feautrier) intensity -- (riin+riup)/2
c
            do id=1,nd
               u(id)=(riin(id)+riup(id))*half
               al0(id)=(aiin(id)+aiup(id))*half
            end do
c
         DO ID=1,ND
            RAD1(ID)=RAD1(ID)+WTMU(I)*U(ID)
            RDK(ID)=RDK(ID)+AMU2*U(ID)
            ALI1(ID)=ALI1(ID)+WTMU(I)*AL0(ID)
         END DO
         AH=AH+AMU(I)*WTMU(I)*U(1)
         ahd=ahd+amu(i)*wtmu(i)*u(nd)
         ahout=ahout+amu(i)*wtmu(i)*riup(1)
c
c     end of the loop over angle points
c
      END DO  
C
C     solution of the transfer equation
C     Variables:
C     RAD1    - mean intensity
C     FAK1    - Eddington factor f(K) = K/J
C     FH      - the "surface" Eddington factor
C     ALI1    - diagonal element of the lambda operator
C
      IF(IBC.EQ.0) THEN
         ALI1(ND)=RAD1(ND)/ST0(ND)
         ALI1(ND-1)=RAD1(ND-1)/ST0(ND-1)
      END IF
C
      DJTOT=0.
      DO ID=1,ND
         DELTAJ=(RAD1(ID)-RAD(IJ,ID))/(UN+SS0(ID)*ALI1(ID))
         RAD(IJ,ID)=RAD(IJ,ID)+DELTAJ
         DJTOT=MAX(DJTOT,ABS(DELTAJ/RAD(IJ,ID)))
      END DO
      IF(DJTOT.GT.DJMAX.AND.ITRALI.LE.NTRALI) GO TO 10
C
C     end of ALI loop for electron scattering
C
      DO ID=1,ND
         RAD1(ID)=RAD(IJ,ID)
         FAK1(ID)=RDK(ID)/RAD(IJ,ID)
         FKK(ID)=FAK1(ID)
      END DO
      FLUX(IJ)=AHout*half
      FH(IJ)=AH/RAD1(1)-HALF*ALB1
      FH0=FH(IJ)
      fhd(ij)=ahd/rad1(nd)
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
      if(ilmcor.eq.3) then
         do id=1,nd
            sa0(id)=st0(id)
            ss0c(id)=ss0(id)
            ss0(id)=0.
         end do
      end if
      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+FH0+B)*BQ+SS0(ID)
      CCC(ID)=(FKK(ID+1)/DTP1)*BQ-CQ*(UN+SS0(ID+1))
      ZZZ(ID)=UN/BBB(ID)
      VLL=SA0(ID)+CQ*SA0(ID+1)
      IF(IWINBL.LT.0) VLL=VLL+HEXTRD(IJ)*BQ
      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)
         AL=UN/DTM1*DT0
         GA=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*AL)*SIXTH
            C=(UN-HALF*DTM1*DTM1*GA)*SIXTH
         END IF
         AAA(ID)=AL*FKK(ID-1)-A*(UN+SS0(ID-1))
         CCC(ID)=GA*FKK(ID+1)-C*(UN+SS0(ID+1))
         BBB(ID)=(AL+GA)*FKK(ID)+(UN-A-C)*(UN+SS0(ID))
         VLL=A*SA0(ID-1)+C*SA0(ID+1)+(UN-A-C)*SA0(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
      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+A*FKK(ID)
         AAA(ID)=A*FKK(ID-1)
         VLL=SA0(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
      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)
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)
         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)/sa0(nd-1)
             ali1(nd)=rad1(nd)/sa0(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)/sa0(nd-1)
           ali1(nd)=rad1(nd)/sa0(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 
           IF(IBC.EQ.4) THEN
             ALI1(ND)= ALI1(ND)/(UN+SS0(ND))
             ALIM1(ND)= ALIM1(ND)/(UN+SS0(ND))
             ALIP1(ND)= ALIP1(ND)/(UN+SS0(ND))
           END IF
         END IF
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
      Q0(IJE)=QQ0
      UU0(IJE)=U0
      RETURN
      END