SUBROUTINE RTEDFE
C     =================
C
C     Solution of the radiative transfer equation - frequency by
C     frequency - for the known source function.
C
C     The numerical method used:
c     Discontinuous Finite Element (DFE) method
c     Castor, Dykema, Klein, 1992, ApJ 387, 561.
C
C     Input through blank COMMON block:
C      CH     - two-dimensional array  absorption coefficient (frequency,
C               depth)
C      ET     - emission coefficient (frequency, depth)
C
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      INCLUDE 'SYNTHP.FOR'
      PARAMETER (ONE=1.,TWO=2.,HALF=0.5)
      PARAMETER (TAUREF = 0.6666666666667)
      DIMENSION DT(MDEPTH),ST0(MDEPTH),AB0(MDEPTH),DELDM(MDEPTH),
     *          dtau(mdepth),rip(mdepth),rim(mdepth),riup(mdepth),
     *          AMU(3),WTMU(3),RINT1(MMU),
     *          AMUI(MMU),AMUW(MMU),TAU(MDEPTH),SS0(MDEPTH)
      COMMON/RTEOPA/CH(MFREQ,MDEPTH),ET(MFREQ,MDEPTH),
     *              SC(MFREQ,MDEPTH)
      COMMON/EMFLUX/FLUX(MFREQ),FLUXC(MFREQC)
      COMMON/CONSCA/SCC1(mdepth),SCC2(MDEPTH)
      COMMON/REFDEP/IREFD(MFREQ)
C
C     angle points (AMU) and angular integration weights (WTMU)
C
      DATA AMU/.887298334620742D0,.5D0,.112701665379258D0/,
     * WTMU/.277777777777778D0,.444444444444444D0,.277777777777778D0/
C
      DO I=1,ND-1
         DELDM(I)=HALF*(DM(I+1)-DM(I))
      END DO
C
c     angle points
C
      IF(IFLUX.EQ.0) THEN
         NMUS=NMU
         do i=1,nmu
            amui(i)=amu(i)
            amuw(i)=amu(i)*wtmu(i)
         end do
       ELSE IF(IFLUX.EQ.1) THEN
         NMUS=NMU0
         do i=1,nmus
            amui(i)=angl(i)
            amuw(i)=angl(i)*wangl(i)
         end do
      END IF
C
C     overall loop over frequencies
C
      DO IJ=1,NFREQ
      FR=FREQ(IJ)
C
C     total source function
C  
      DO ID=1,ND
         AB0(ID)=CH(IJ,ID)
         SCT=FRX1(IJ)*SCC2(ID)+FRX2(IJ)*SCC1(ID)
         ST0(ID)=(ET(IJ,ID)+SCT)/AB0(ID)
         SS0(ID)=-SCT/AB0(ID)
      END DO
      AH=0.
C
C     optical depth scale
C
      TAU(1)=0.
      IREF=1
      DO ID=1,ND-1
         DT(ID)=DELDM(ID)*(AB0(ID+1)/DENS(ID+1)+AB0(ID)/DENS(ID))
         TAU(ID+1)=TAU(ID)+DT(ID)
         IF(TAU(ID).LE.TAUREF.AND.TAU(ID+1).GT.TAUREF) IREF=ID
      END DO
      IREFD(IJ)=IREF
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))-ONE)
      DPLAN=BNU/(EXP(HK*FR/TEMP(ND-1))-ONE)
      DPLAN=(PLAND-DPLAN)/DT(ND-1)
c
c     loop over angle poits
c
      DO I=1,NMUS
         do id=1,nd-1
            dtau(id)=dt(id)/amui(i)
         enddo
C
c           outgoing intensity
c 
            rip(nd)=PLAND+AMUI(I)*DPLAN
            id=nd-1
            dt0=dtau(id)
            dtaup1=dt0+one
            dtau2=dt0*dt0
            bb=two*dtaup1
            cc=dt0*dtaup1
            aa=dtau2+bb
            rim(id+1)=(aa*rip(id+1)-cc*st0(id+1)+dt0*st0(id))/bb
            do id=nd-1,1,-1
               dt0=dtau(id)
               dtaup1=dt0+one
               dtau2=dt0*dt0
               bb=two*dtaup1
               cc=dt0*dtaup1
               aa=one/(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
            enddo
            do id=2,nd-1
               riup(id)=(rim(id)*dtau(id-1)+rip(id)*dtau(id))/
     *                  (dtau(id-1)+dtau(id))
            enddo
            riup(1)=rim(1)
            riup(nd)=rip(nd)
c
         AH=AH+AMUW(I)*RIUP(1)
         RINT1(I)=RIUP(1)
         rint1(i)=max(rint1(i),1.e-40)
c
c     end of the loop over angle points
c
      END DO
c
      FLUX(IJ)=AH*HALF
      if(iflux.ge.1) then
C
C     output of emergent specific intensities to Unit 10 (line points)
C     or 18 (two continuum points)
C
      IF(IJ.GT.2) THEN
         WRITE(10,618) WLAM(IJ),FLUX(IJ),(RINT1(IMU),IMU=1,NMUS)
       ELSE
         WRITE(18,618) WLAM(IJ),FLUX(IJ),(RINT1(IMU),IMU=1,NMUS)
      END IF
      end if
  618 FORMAT(1H ,f10.3,1pe15.5/(1P5E15.5))
C
C     if needed (if iprin.ge.3), output of interesting physical
C     quantities at the monochromatic optical depth  tau(nu)=2/3
C
      IF(IPRIN.GE.3) THEN
         T0=LOG(TAU(IREF+1)/TAU(IREF))
         X0=LOG(TAU(IREF+1)/TAUREF)/T0
         X1=LOG(TAUREF/TAU(IREF))/T0
         DMREF=EXP(LOG(DM(IREF))*X0+LOG(DM(IREF+1))*X1)
         TREF=EXP(LOG(TEMP(IREF))*X0+LOG(TEMP(IREF+1))*X1)
         STREF=EXP(LOG(ST0(IREF))*X0+LOG(ST0(IREF+1))*X1)
         SSREF=EXP(LOG(-SS0(IREF))*X0+LOG(-SS0(IREF+1))*X1)
         SREF=STREF+SSREF
         ALM=2.997925E18/FREQ(IJ)
         WRITE(96,636) IJ,ALM,IREF,DMREF,TREF,STREF,SSREF,SREF
  636    FORMAT(1H ,I3,F10.3,I4,1PE10.3,0PF10.1,1X,1P3E10.3)
      END IF
C
C     end of the loop over frequencies
C
      END DO
      RETURN
      END