SUBROUTINE RTEDF2(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     Analogous to RTEDF1, but using opacity and emissivity
c     instead of the source function
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'
      PARAMETER (SIXTH=UN/6.D0,
     *           THIRD=UN/3.D0,
     *           TWOTHR=TWO/3.D0,
     *           three=3.d0, 
     *           quart=0.25d0)
      DIMENSION DT(MDEPTH),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),
     *          chip0(mdepth),chim0(mdepth),ddm0(mdepth),
     *          chip(mdepth),chim(mdepth),
     *          etap(mdepth),etam(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)
         ddm0(id)=(dm(id+1)-dm(id))
         chip0(id)=(abso1(id)*dens1(id)+three*
     *              abso1(id+1)*dens1(id+1))*quart*ddm0(id)
         chim0(id)=(abso1(id)*dens1(id)*three+
     *              abso1(id+1)*dens1(id+1))*quart*ddm0(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.
      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)
            chip(id)=un+chip0(id)/amu(i)
            chim(id)=un+chim0(id)/amu(i)
            etap(id)=emis1(id+1)*dens1(id+1)/amu(i)*ddm0(id)
            etam(id)=emis1(id)*dens1(id)/amu(i)*ddm0(id)
         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)
         IF(IWINBL.EQ.0) THEN
            TAMM=TAUMIN/AMU(I)
            EX=EXP(-TAMM)
            P0=UN-EX
            QQ0=QQ0+P0*AMU(I)*WTMU(I)
            U0=U0+EX*WTMU(I)
            US0=US0+P0/TAMM*WTMU(I)
            rim(id)=st0(id)*p0
         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)
               aam=un/(un+chim(id)*chip(id))
               rim(id+1)=(two*rim(id)+etap(id)*chim(id)+etam(id))*aam
               rip(id)=(two*rim(id)*chim(id)+etam(id)*chip(id)-
     *                  etap(id))*aam
               aim(id+1)=bb*aa
               aip(id)=(cc+bb*aim(id))*aa
            end do
            do id=2,nd-1
               dtt=un/(dtau(id-1)+dtau(id))
               dtm=un/(ddm0(id-1)+ddm0(id))
               riin(id)=(rim(id)*ddm0(id)+rip(id)*ddm0(id-1))*dtm
               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               
            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)
               aam=un/(un+chim(id)*chip(id))
               rim(id)=(two*rim(id+1)+etam(id)*chip(id)+etap(id))*aam
               rip(id+1)=(two*rim(id+1)*chip(id)+etap(id)*chim(id)-
     *                    etam(id))*aam
               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))
               dtm=un/(ddm0(id-1)+ddm0(id))
               riup(id)=(rim(id)*ddm0(id-1)+rip(id)*ddm0(id))*dtm
               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)
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)=AH
      FH(IJ)=AH/RAD1(1)-HALF*ALB1
      FH0=FH(IJ)
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))
         END DO   
         IF(IBC.EQ.4) THEN
             ALI1(ND)= ALI1(ND)/(UN+SS0(ND))
         END IF
         ELSE IF(ILMCOR.EQ.3) THEN     
         DO ID=1,ND
             ALI1(ID)=ALI1(ID)/(UN+SS0(ID)*ALI1(ID))
         END DO   
         IF(IBC.EQ.4) THEN
             ALI1(ND)= ALI1(ND)*(UN+SS0(ND)*ALI1(ID))
         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