SUBROUTINE RTECD
C     ================
C
C     solution of the radiative transfer equation by Feautrier method
C     for two continuum points 
C     used when one employs RTEDFE, ie. the DFE method for the
C     transfer equation for the inner frequency points
C
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      INCLUDE 'SYNTHP.FOR'
      DIMENSION D(3,3,MDEPTH),ANU(3,MDEPTH),AANU(MDEPTH),DDD(MDEPTH),
     *       AA(3,3),BB(3,3),CC(3,3),VL(3),AMU(3),WTMU(3),
     *       DT(MDEPTH),TAU(MDEPTH),
     *       RDD(MDEPTH),FKK(MDEPTH),ST0(MDEPTH),SS0(MDEPTH),
     *       RINT(MDEPTH,MMU)
      COMMON/RTEOPA/CH(MFREQ,MDEPTH),ET(MFREQ,MDEPTH),
     *              SC(MFREQ,MDEPTH)
      COMMON/EMFLUX/FLUX(MFREQ),FLUXC(MFREQC)
      COMMON/CONSCA/SCC1(mdepth),SCC2(MDEPTH)
      PARAMETER (UN=1.D0, HALF=0.5D0)
      PARAMETER (THIRD=UN/3., QUART=UN/4., SIXTH=UN/6.D0)
      PARAMETER (TAUREF = 0.6666666666667)
      DATA AMU/.887298334620742D0,.5D0,.112701665379258D0/,
     1     WTMU/.277777777777778D0,.444444444444444D0,.277777777777778D0
     1         /
C
      NMU=3
      ND1=ND-1
C
C     loop over two continuum frequencies 
C
      DO 100 IJ=1,2       
      TAUMIN=CH(IJ,1)/DENS(1)*DM(1)*HALF
      TAU(1)=TAUMIN
      DO I=1,ND1
         DT(I)=(DM(I+1)-DM(I))*(CH(IJ,I+1)/DENS(I+1)+CH(IJ,I)/DENS(I))*
     *         HALF
         ST0(I)=ET(IJ,I)/CH(IJ,I)
         SS0(I)=-SC(IJ,I)/CH(IJ,I)
         TAU(I+1)=TAU(I)+DT(I)
         IF(TAU(I).LE.TAUREF.AND.TAU(I+1).GT.TAUREF) IREF=I
      END DO
      ST0(ND)=ET(IJ,ND)/CH(IJ,ND)
      SS0(ND)=-SC(IJ,ND)/CH(IJ,ND)
      FR=FREQ(IJ)
      BNU=BN*(FR*1.E-15)**3
      PLAND=BNU/(EXP(HK*FR/TEMP(ND  ))-UN)
      DPLAN=BNU/(EXP(HK*FR/TEMP(ND-1))-UN)
      DPLAN=(PLAND-DPLAN)/DT(ND1)
C
C   +++++++++++++++++++++++++++++++++++++++++
C   FIRST PART  -  VARIABLE EDDINGTON FACTORS
C   +++++++++++++++++++++++++++++++++++++++++
C
C   Allowance for wind blanketing
C
      ALB1=0.
      DO I=1,NMU
C
C   ************************
C   UPPER BOUNDARY CONDITION
C   ************************
C
         ID=1
         DTP1=DT(1)
         Q0=0.
         P0=0.
C
C        allowance for non-zero optical depth at the first depth point
C
         TAMM=TAUMIN/AMU(I)
         IF(TAMM.GT.0.01) THEN
            P0=UN-EXP(-TAMM)
          ELSE
            P0=TAMM*(UN-HALF*TAMM*(UN-TAMM*THIRD*(UN-QUART*TAMM)))
         END IF
         EX=UN-P0
         Q0=Q0+P0*AMU(I)*WTMU(I)
C
         DIV=DTP1/AMU(I)*THIRD
         VL(I)=DIV*(ST0(ID)+HALF*ST0(ID+1))+ST0(ID)*P0
         DO J=1,NMU
            BB(I,J)=SS0(ID)*WTMU(J)*(DIV+P0)-ALB1*WTMU(J)
            CC(I,J)=-HALF*DIV*SS0(ID+1)*WTMU(J)
         END DO
         BB(I,I)=BB(I,I)+AMU(I)/DTP1+UN+DIV
         CC(I,I)=CC(I,I)+AMU(I)/DTP1-HALF*DIV
         ANU(I,ID)=0.
      END DO
C
C     Matrix inversion: instead of calling MATINV, a very fast inlined 
C     routine MINV3 for a specific 3 x 3 matrix inversion
C
C     CALL MATINV(BB,NMU,3)
C
C     ******************************
      BB(2,1)=BB(2,1)/BB(1,1)
      BB(2,2)=BB(2,2)-BB(2,1)*BB(1,2)
      BB(2,3)=BB(2,3)-BB(2,1)*BB(1,3)
      BB(3,1)=BB(3,1)/BB(1,1)
      BB(3,2)=(BB(3,2)-BB(3,1)*BB(1,2))/BB(2,2)
      BB(3,3)=BB(3,3)-BB(3,1)*BB(1,3)-BB(3,2)*BB(2,3)
C
      BB(3,2)=-BB(3,2)
      BB(3,1)=-BB(3,1)-BB(3,2)*BB(2,1)
      BB(2,1)=-BB(2,1)
C
      BB(3,3)=UN/BB(3,3)
      BB(2,3)=-BB(2,3)*BB(3,3)/BB(2,2)
      BB(2,2)=UN/BB(2,2)
      BB(1,3)=-(BB(1,2)*BB(2,3)+BB(1,3)*BB(3,3))/BB(1,1)
      BB(1,2)=-BB(1,2)*BB(2,2)/BB(1,1)
      BB(1,1)=UN/BB(1,1)
C
      BB(1,1)=BB(1,1)+BB(1,2)*BB(2,1)+BB(1,3)*BB(3,1)
      BB(1,2)=BB(1,2)+BB(1,3)*BB(3,2)
      BB(2,1)=BB(2,2)*BB(2,1)+BB(2,3)*BB(3,1)
      BB(2,2)=BB(2,2)+BB(2,3)*BB(3,2)
      BB(3,1)=BB(3,3)*BB(3,1)
      BB(3,2)=BB(3,3)*BB(3,2)
C     ******************************
C
      DO I=1,NMU
         DO J=1,NMU
            S=0.
            DO K=1,NMU
               S=S+BB(I,K)*CC(K,J)
            END DO
            D(I,J,ID)=S
            ANU(I,1)=ANU(I,1)+BB(I,J)*VL(J)
         END DO
      END DO
C
C   *******************
C   NORMAL DEPTH POINTS
C   *******************
C
      DO ID=2,ND1
         DTM1=DTP1
         DTP1=DT(ID)
         DT0=HALF*(DTM1+DTP1)
         AL=UN/DTM1/DT0
         GA=UN/DTP1/DT0
         BE=AL+GA
         A=(UN-HALF*AL*DTP1*DTP1)*SIXTH
         C=(UN-HALF*GA*DTM1*DTM1)*SIXTH
         B=UN-A-C
         VL0=A*ST0(ID-1)+B*ST0(ID)+C*ST0(ID+1)
         DO I=1,NMU
            DO J=1,NMU
               AA(I,J)=-A*SS0(ID-1)*WTMU(J)
               CC(I,J)=-C*SS0(ID+1)*WTMU(J)
               BB(I,J)=B*SS0(ID)*WTMU(J)
            END DO
         END DO
         DO I=1,NMU
            DIV=AMU(I)**2
            VL(I)=VL0
            AA(I,I)=AA(I,I)+DIV*AL-A
            CC(I,I)=CC(I,I)+DIV*GA-C
            BB(I,I)=BB(I,I)+DIV*BE+B
         END DO
         DO I=1,NMU
            S1=0.
            DO J=1,NMU
               S=0.
               S1=S1+AA(I,J)*ANU(J,ID-1)
               DO K=1,NMU
                  S=S+AA(I,K)*D(K,J,ID-1)
               END DO
               BB(I,J)=BB(I,J)-S
            END DO
            VL(I)=VL(I)+S1
         END DO
C
C     Matrix inversion: instead of calling MATINV, a very fast inlined 
C     routine MINV3 for a specific 3 x 3 matrix inversion
C
C     CALL MATINV(BB,NMU,3)
C
C     ******************************
      BB(2,1)=BB(2,1)/BB(1,1)
      BB(2,2)=BB(2,2)-BB(2,1)*BB(1,2)
      BB(2,3)=BB(2,3)-BB(2,1)*BB(1,3)
      BB(3,1)=BB(3,1)/BB(1,1)
      BB(3,2)=(BB(3,2)-BB(3,1)*BB(1,2))/BB(2,2)
      BB(3,3)=BB(3,3)-BB(3,1)*BB(1,3)-BB(3,2)*BB(2,3)
C
      BB(3,2)=-BB(3,2)
      BB(3,1)=-BB(3,1)-BB(3,2)*BB(2,1)
      BB(2,1)=-BB(2,1)
C
      BB(3,3)=UN/BB(3,3)
      BB(2,3)=-BB(2,3)*BB(3,3)/BB(2,2)
      BB(2,2)=UN/BB(2,2)
      BB(1,3)=-(BB(1,2)*BB(2,3)+BB(1,3)*BB(3,3))/BB(1,1)
      BB(1,2)=-BB(1,2)*BB(2,2)/BB(1,1)
      BB(1,1)=UN/BB(1,1)
C
      BB(1,1)=BB(1,1)+BB(1,2)*BB(2,1)+BB(1,3)*BB(3,1)
      BB(1,2)=BB(1,2)+BB(1,3)*BB(3,2)
      BB(2,1)=BB(2,2)*BB(2,1)+BB(2,3)*BB(3,1)
      BB(2,2)=BB(2,2)+BB(2,3)*BB(3,2)
      BB(3,1)=BB(3,3)*BB(3,1)
      BB(3,2)=BB(3,3)*BB(3,2)
C     ******************************
C
         DO I=1,NMU
            ANU(I,ID)=0.
            DO J=1,NMU
               S=0.
               DO K=1,NMU
                  S=S+BB(I,K)*CC(K,J)
               END DO
               D(I,J,ID)=S
               ANU(I,ID)=ANU(I,ID)+BB(I,J)*VL(J)
            END DO
         END DO
      END DO
C
C   ************
C   LOWER BOUNDARY CONDITION
C   ************
C
      ID=ND
      DO I=1,NMU
         AA(I,I)=AMU(I)/DTP1
         VL(I)=PLAND+AMU(I)*DPLAN+AA(I,I)*ANU(I,ID-1)
         DO J=1,NMU
            BB(I,J)=-AA(I,I)*D(I,J,ID-1)
         END DO
         BB(I,I)=BB(I,I)+AA(I,I)+UN
      END DO
C
C     Matrix inversion: instead of calling MATINV, a very fast inlined 
C     routine MINV3 for a specific 3 x 3 matrix inversion
C
C     CALL MATINV(BB,NMU,3)
C
C     ******************************
      BB(2,1)=BB(2,1)/BB(1,1)
      BB(2,2)=BB(2,2)-BB(2,1)*BB(1,2)
      BB(2,3)=BB(2,3)-BB(2,1)*BB(1,3)
      BB(3,1)=BB(3,1)/BB(1,1)
      BB(3,2)=(BB(3,2)-BB(3,1)*BB(1,2))/BB(2,2)
      BB(3,3)=BB(3,3)-BB(3,1)*BB(1,3)-BB(3,2)*BB(2,3)
C
      BB(3,2)=-BB(3,2)
      BB(3,1)=-BB(3,1)-BB(3,2)*BB(2,1)
      BB(2,1)=-BB(2,1)
C
      BB(3,3)=UN/BB(3,3)
      BB(2,3)=-BB(2,3)*BB(3,3)/BB(2,2)
      BB(2,2)=UN/BB(2,2)
      BB(1,3)=-(BB(1,2)*BB(2,3)+BB(1,3)*BB(3,3))/BB(1,1)
      BB(1,2)=-BB(1,2)*BB(2,2)/BB(1,1)
      BB(1,1)=UN/BB(1,1)
C
      BB(1,1)=BB(1,1)+BB(1,2)*BB(2,1)+BB(1,3)*BB(3,1)
      BB(1,2)=BB(1,2)+BB(1,3)*BB(3,2)
      BB(2,1)=BB(2,2)*BB(2,1)+BB(2,3)*BB(3,1)
      BB(2,2)=BB(2,2)+BB(2,3)*BB(3,2)
      BB(3,1)=BB(3,3)*BB(3,1)
      BB(3,2)=BB(3,3)*BB(3,2)
C     ******************************
C
      DO I=1,NMU
         ANU(I,ID)=0.
         DO J=1,NMU
            D(I,J,ID)=0.
            ANU(I,ID)=ANU(I,ID)+BB(I,J)*VL(J)
         END DO
      END DO
C
C   ************
C   BACKSOLUTION
C   ************
C
      DO ID=ND-1,1,-1
         DO I=1,NMU
            DO J=1,NMU
               ANU(I,ID)=ANU(I,ID)+D(I,J,ID)*ANU(J,ID+1)
            END DO
         END DO
         AJ=0.
         AK=0.
         DO I=1,NMU
            DIV=WTMU(I)*ANU(I,ID)
            AJ=AJ+DIV
            AK=AK+DIV*AMU(I)**2
         END DO
         FKK(ID)=AK/AJ
      END DO
C
C     surface Eddington actor
C
      AH=0.
      DO I=1,NMU
         AH=AH+WTMU(I)*AMU(I)*ANU(I,1)
      END DO
      FH=AH/AJ-HALF*ALB1
C
      FKK(ND)=THIRD
C
C
C   +++++++++++++++++++++++++++++++++++++++++
C   SECOND PART  -  DETERMINATION OF THE MEAN INTENSITIES
C   RECALCULATION OF THE TRANSFER EQUATION WITH GIVEN EDDINGTON FACTORS
C   +++++++++++++++++++++++++++++++++++++++++
C
      DTP1=DT(1)
      DIV=DTP1*THIRD
      BBB=FKK(1)/DTP1+FH+DIV+SS0(1)*(DIV+Q0)
      CCC=FKK(2)/DTP1-HALF*DIV*(UN+SS0(2))
      VLL=DIV*(ST0(1)+HALF*ST0(2))+ST0(1)*Q0
      AANU(1)=VLL/BBB
      DDD(1)=CCC/BBB
      DO ID=2,ND1
         DTM1=DTP1
         DTP1=DT(ID)
         DT0=HALF*(DTP1+DTM1)
         AL=UN/DTM1/DT0
         GA=UN/DTP1/DT0
         A=(UN-HALF*DTP1*DTP1*AL)*SIXTH
         C=(UN-HALF*DTM1*DTM1*GA)*SIXTH
         AAA=AL*FKK(ID-1)-A*(UN+SS0(ID-1))
         CCC=GA*FKK(ID+1)-C*(UN+SS0(ID+1))
         BBB=(AL+GA)*FKK(ID)+(UN-A-C)*(UN+SS0(ID))
         VLL=A*ST0(ID-1)+C*ST0(ID+1)+(UN-A-C)*ST0(ID)
         BBB=BBB-AAA*DDD(ID-1)
         DDD(ID)=CCC/BBB
         AANU(ID)=(VLL+AAA*AANU(ID-1))/BBB
      END DO
      BBB=FKK(ND)/DTP1+HALF
      AAA=FKK(ND1)/DTP1
      BBB=BBB-AAA*DDD(ND1)
      VLL=HALF*PLAND+DPLAN*THIRD
      RDD(ND)=(VLL+AAA*AANU(ND1))/BBB
      DO IID=1,ND1
         ID=ND-IID
         RDD(ID)=AANU(ID)+DDD(ID)*RDD(ID+1)
      END DO
      FLUX(IJ)=FH*RDD(1)
C
      if(ij.eq.1) then
         do id=1,nd
            scc1(id)=-rdd(id)*ss0(id)*ch(1,id)
         end do
       else
         do id=1,nd
            scc2(id)=-rdd(id)*ss0(id)*ch(2,id)
         end do
      end if
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)
      SCREF=EXP(LOG(-SS0(IREF))*X0+LOG(-SS0(IREF+1))*X1)
      SSREF=EXP(LOG(-SS0(IREF)*RDD(IREF))*X0+
     *           LOG(-SS0(IREF+1)*RDD(IREF+1))*X1)
      SREF=STREF+SSREF
      ALM=2.997925E18/FREQ(IJ)
      WRITE(96,636) IJ,ALM,IREF,DMREF,TREF,SCREF,STREF,SSREF,SREF
  636 FORMAT(1H ,I3,F10.3,I4,1PE10.3,0PF10.1,1X,1P3E10.3,E11.3)
      END IF
c
C   ********************************************************************
C
C   THIRD PART  -  DETERMINATION OF THE SPECIFIC INTENSITIES
C   RECALCULATION OF THE TRANSFER EQUATION WITH GIVEN SOURCE FUNCTION
C
      if(iflux.eq.0) go to 100
      DO IMU=1,NMU0
         ANX=ANGL(IMU)
         DTP1=DT(1)
         DIV=DTP1*THIRD/ANX
C
         TAMM=TAUMIN/ANX
         IF(TAMM.LT.0.01) THEN
            P0=TAMM*(UN-HALF*TAMM*(UN-TAMM*THIRD*(UN-QUART*TAMM)))
          ELSE
            P0=UN-EXP(-TAMM)
         END IF
C
         BBB=ANX/DTP1+UN+DIV
         CCC=ANX/DTP1-HALF*DIV
         VLL=(DIV+P0)*(ST0(1)-SS0(1)*RDD(1))
     *       +HALF*DIV*(ST0(2)-SS0(2)*RDD(2))
         AANU(1)=VLL/BBB
         DDD(1)=CCC/BBB
         DIV=ANX*ANX
         DO ID=2,ND1
            DTM1=DT(ID-1)
            DTP1=DT(ID)
            DT0=HALF*(DTP1+DTM1)
            AL=UN/DTM1/DT0
            GA=UN/DTP1/DT0
            A=(UN-HALF*DTP1*DTP1*AL)*SIXTH
            C=(UN-HALF*DTM1*DTM1*GA)*SIXTH
            AAA=DIV*AL-A
            CCC=DIV*GA-C
            BBB=DIV*(AL+GA)+UN-A-C
            VLL=A*(ST0(ID-1)-SS0(ID-1)*RDD(ID-1))+
     *          C*(ST0(ID+1)-SS0(ID+1)*RDD(ID+1))+
     *          (UN-A-C)*(ST0(ID)-SS0(ID)*RDD(ID))
            BBB=BBB-AAA*DDD(ID-1)
            DDD(ID)=CCC/BBB
            AANU(ID)=(VLL+AAA*AANU(ID-1))/BBB
         END DO
C
C        Lower boundary condition
C
         AAA=ANX/DTP1
         BBB=AAA+UN
         VLL=PLAND+ANX*DPLAN
C
         RINT(ND,IMU)=(VLL+AAA*AANU(ND1))/(BBB-AAA*DDD(ND1))
         DO IID=1,ND1
            ID=ND-IID
            RINT(ID,IMU)=AANU(ID)+DDD(ID)*RINT(ID+1,IMU)
         END DO
      END DO
c
      FLX=0.
      DO IMU=1,NMU0
         RINT(1,IMU)=RINT(1,IMU)/HALF
         FLX=FLX+ANGL(IMU)*WANGL(IMU)*RINT(1,IMU)
      END DO
      FLX=FLX*HALF
C
C     output of emergent specific intensities in continuum to Unit 18 
C
      if(iflux.ge.1) then
      WRITE(18,641) WLAM(IJ),FLX,(RINT(1,IMU),IMU=1,NMU0)
      end if
  100 CONTINUE
  641 FORMAT(1H ,f10.3,1pe15.5/(1P5E15.5))
c  
c     call rtedfe for the internal points
c
      CALL RTEDFE
c
      RETURN
      END