SUBROUTINE RTEINT
C     =================
C
C     Solution of the radiative transfer equation - for one frequency
C     - for the known source function.
C
C     Solution for specific intensities
C
C     The numerical method used:
C      for ISPLIN = 0  -  the ordinary Feautrier scheme
C                 = 1  -  the spline collocation method
C                 = 2  -  Hermitian fourth-order method
C                 = 3  -  improved Feautrier scheme
C                         (Rybicki & Hummer 1991, A&A 245, 171.)
C
C     In all cases, the overall matrix system is solved by the standard
C     Gaussian elimination, analogous to that described in SOLVE
C     (auxiliary matrix D is called ALF in SOLVE; auxiliary vector ANU
C     is called BET in SOLVE)
C
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)
      PARAMETER (MMA=20)
      DIMENSION ST0(MDEPTH),SS0(MDEPTH),AB0(MDEPTH),
     *          AA(MMA,MMA),BB(MMA,MMA),CC(MMA,MMA),VL(MMA),
     *          FFD(MMA,MMA),FF0D(MMA,MMA),FFPD(MMA,MMA),
     *          D(MMA,MMA,MDEPTH),ANU(MMA,MDEPTH)
      DIMENSION tau(mdepth),angl(mma),WANG(MMA)
      COMMON/OPTDPT/DT(MDEPTH)
C
      nmuf=nmu
      nmu=intens
      do imu=1,nmu
         angl(imu)=0.1+float(imu-1)*0.9/float(nmu-1)
         wang(imu)=0.9/float(nmu-1)
      end do
      wang(1)=wang(1)*0.5
      wang(nmu)=WANG(NMU)*0.5
C
      DO 500 IJO=1,NFREQ
      IJ=IJO
      IF(ispodf.eq.0) IJ=JIK(IJO)
      IF(IJX(IJ).EQ.-1) GO TO 500
      call opacf1(IJ)
c
      FR=FREQ(IJ)
C
C     total source function
C
      AH=0.
      DO ID=1,ND
         AB0(ID)=ABSO1(ID)
         ST0(ID)=EMIS1(ID)/AB0(ID)
         SS0(ID)=-SCAT1(ID)/AB0(ID)
         RAD1(ID)=0.
      END DO
C
C     optical depth scale
C
      tau(1)=absot(1)*dm(1)
      DO ID=1,ND-1
         DT(ID)=DELDMZ(ID)*(ABSOT(ID+1)+ABSOT(ID))
         tau(id+1)=tau(id)+dt(id)
      END DO
c
      U0=0.
      QQ0=0.
      US0=0.
      taumin=absot(1)*dm(1)/2.
C
      ALB1=0.
C
C ************** Forward elimination
C
C     Upper boundary condition
C
      ID=1
      DTP1=DT(1)
      P0=0.
      EX=UN
      IF(MOD(ISPLIN,3).EQ.0) THEN
         B=DTP1*HALF
         C=0.
       ELSE
         B=DTP1*THIRD
         C=B*HALF
      END IF
      QQ0=0.
      US0=0.
      DO I=1,NMU
         IF(IDISK.EQ.0) THEN
C
C        allowance for non-zero optical depth at the first depth point
C
             TAMM=TAUMIN/ANGL(I)
             EX=EXP(-TAMM)
             P0=UN-EX
             QQ0=QQ0+P0*ANGL(I)*WANG(I)
             U0=U0+EX*WANG(I)
             US0=US0+P0/TAMM*WANG(I)
          END IF
C
         BI=B/ANGL(I)
         CI=C/ANGL(I)
         VL(I)=(BI+P0)*ST0(ID)+CI*ST0(ID+1)
         IF(IWINBL.LT.0) VL(I)=VL(I)+EXTRAD(IJ)
         DO J=1,NMU
            BB(I,J)=SS0(ID)*WANG(J)*(BI+P0)-ALB1*WANG(J)
            CC(I,J)=-CI*SS0(ID+1)*WANG(J)
         END DO
         BB(I,I)=BB(I,I)+ANGL(I)/DTP1+UN+BI
         CC(I,I)=CC(I,I)+ANGL(I)/DTP1-CI
         ANU(I,ID)=0.
      END DO
C
      IF(ISPLIN.LE.2) THEN
         CALL MATINV(BB,NMU,MMA)
         DO I=1,NMU
            DO J=1,NMU
               D(I,J,ID)=0.
               DO K=1,NMU
                  D(I,J,ID)=D(I,J,ID)+BB(I,K)*CC(K,J)
               END DO
               ANU(I,1)=ANU(I,1)+BB(I,J)*VL(J)
            END DO
         END DO
      ELSE
         DO I=1,NMU
            DO J=1,NMU
               FF0D(I,J)=BB(I,J)/CC(I,I)
            END DO
            FF0D(I,I)=FF0D(I,I)-UN
         END DO
C
         CALL MATINV(BB,NMU,MMA)
         DO I=1,NMU
            ANU(I,ID)=0.
            DO J=1,NMU
               D(I,J,ID)=BB(I,J)*CC(J,J)
               ANU(I,ID)=ANU(I,ID)+BB(I,J)*VL(J)
            END DO
         END DO
      ENDIF
C
C     Normal depth points  1 < ID < ND
C
      DO ID=2,ND-1
         DTM1=DTP1
         DTP1=DT(ID)
         DT0=TWO/(DTM1+DTP1)
         AL=UN/DTM1*DT0
         GA=UN/DTP1*DT0
         BE=AL+GA
         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*AL*DTP1*DTP1)*SIXTH
            C=(UN-HALF*GA*DTM1*DTM1)*SIXTH
         END IF
         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)*WANG(J)
               CC(I,J)=-C*SS0(ID+1)*WANG(J)
               BB(I,J)=B*SS0(ID)*WANG(J)
            END DO
         END DO
         DO I=1,NMU
            VL(I)=VL0
            DIV=ANGL(I)*ANGL(I)
            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
            DO J=1,NMU
               VL(I)=VL(I)+AA(I,J)*ANU(J,ID-1)
            END DO
         END DO
         IF(ISPLIN.LE.2) THEN
            DO I=1,NMU
               DO J=1,NMU
                  S=0.
                  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
            END DO
C
            CALL MATINV(BB,NMU,MMA)
            DO I=1,NMU
               DO J=1,NMU
                  D(I,J,ID)=0.
                  DO K=1,NMU
                     D(I,J,ID)=D(I,J,ID)+BB(I,K)*CC(K,J)
                  END DO
               END DO
            END DO
         ELSE
            DO I=1,NMU
               BB(I,I)=-AA(I,I)+BB(I,I)-CC(I,I)
               DO J=1,NMU
                  FFPD(I,J)=AA(I,I)*FF0D(I,J)
               END DO
            END DO
            DO I=1,NMU
               DO J=1,NMU
                  S=0.
                  DO K=1,NMU
                     S=S+FFPD(I,K)*D(K,J,ID-1)
                  END DO
                  FFD(I,J)=(BB(I,J)+S)/CC(I,I)
               END DO
            END DO
            DO I=1,NMU
               DO J=1,NMU
                  FF0D(I,J)=FFD(I,J)
               END DO
               FFD(I,I)=FFD(I,I)+UN
            END DO
C
            CALL MATINV(FFD,NMU,MMA)
            DO I=1,NMU
               DO J=1,NMU
                  D(I,J,ID)=FFD(I,J)
                  BB(I,J)=FFD(I,J)/CC(J,J)
               END DO
            END DO
         END IF
         DO I=1,NMU
            ANU(I,ID)=0.
            DO J=1,NMU
               ANU(I,ID)=ANU(I,ID)+BB(I,J)*VL(J)
            END DO
         END DO
      END DO
C
C     Lower boundary condition
C
      ID=ND
C
C     First option:
C     b.c. is different from stellar atmospheres; expresses symmetry
C     at the central plane   I(taumax,-mu,nu)=I(taumax,+mu,nu)
C
      IF(IFZ0.GE.0.AND.IDISK.EQ.1) THEN
         B=DTP1*HALF
         A=0.
         DO I=1,NMU
            BI=B/ANGL(I)
            AI=A/ANGL(I)
            VL(I)=ST0(ID)*BI+ST0(ID-1)*AI
            DO J=1,NMU
               AA(I,J)=-AI*SS0(ID-1)*WANG(J)
               BB(I,J)=BI*SS0(ID)*WANG(J)
            END DO
            AA(I,I)=AA(I,I)+ANGL(I)/DTP1-AI
            BB(I,I)=BB(I,I)+ANGL(I)/DTP1+BI
         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     Second option:
C     b.c. is the same as in stellar atmospheres - the last depth point
C     is not at the central plane
C
       ELSE
         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
         END IF
         DPLAN=(PLAND-DPLAN)/DT(ND-1)
         IF(IBC.EQ.0.OR.IBC.EQ.4) THEN
            DO I=1,NMU
               AA(I,I)=ANGL(I)/DTP1
               VL(I)=PLAND+ANGL(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
          ELSE
            DO I=1,NMU
               A=ANGL(I)/DTP1
               B=HALF/A
               AA(I,I)=A
               VL(I)=B*ST0(ID)+PLAND+ANGL(I)*DPLAN+AA(I,I)*ANU(I,ID-1)
               DO J=1,NMU
                  BB(I,J)=B*SS0(ID)*WANG(J)-AA(I,I)*D(I,J,ID-1)
               END DO
               BB(I,I)=BB(I,I)+A+B+UN
            END DO
         END IF
      END IF
C
      CALL MATINV(BB,NMU,MMA)
C
      DO 230 I=1,NMU
         ANU(I,ID)=0.
         DO 230 J=1,NMU
            D(I,J,ID)=0.
            ANU(I,ID)=ANU(I,ID)+BB(I,J)*VL(J)
  230 CONTINUE
C
C ***************** Backsolution
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
      END DO
c 
      sum=0.
      sua=0.
      do imu=1,nmu
         sum=sum+anu(imu,1)*angl(imu)*wang(imu)
         sua=sua+angl(imu)*wang(imu)
      end do
C
      wlam=2.997925e18/freq(ij)
      WRITE(18,641) WLAM,flux(ij),sum,sua,(2.*ANU(IMU,1),IMU=1,NMU)
  641 format(f11.3,(1p13e11.3))
c
  500 continue
      nmu=nmuf
c
      return
      end