SUBROUTINE RTEFR1(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      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      U0     - derivative of Q0 wrt taumin;
C               for "fixed" frequencies, U0 has the meaning of
C               absorption coefficient * second moment H
C               ( a quantity needed for lower boundary condition of the
C               hydrostatic equilibrium equation, specifically for
C               accounting for an effect of fixed-option transitions)
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)
      DIMENSION AANU(MDEPTH),DDD(MDEPTH),FKK(MDEPTH),
     *          RDD(MDEPTH),ST0(MDEPTH),SS0(MDEPTH),AB0(MDEPTH),
     *          AA(MMU,MMU),BB(MMU,MMU),CC(MMU,MMU),VL(MMU),
     *          FFD(MMU,MMU),FF0D(MMU,MMU),
     *          FFPD(MMU,MMU),ali0(mdepth),ss0c(mdepth),
     *          AAA(MDEPTH),BBB(MDEPTH),CCC(MDEPTH),EEE(MDEPTH),
     *          ZZZ(MDEPTH),ALRH(MDEPTH),ALRM(MDEPTH),ALRP(MDEPTH),
     *          D(MMU,MMU,MDEPTH),ANU(MMU,MDEPTH),scor(mdepth)
      DIMENSION rmmu(2*MMU),wmmu(2*MMU),rwmu(2*MMU),
     *          dtau(mdepth),ri(mdepth),ali(mdepth),alij1(mdepth)
      dimension tau(mdepth)
      COMMON/OPTDPT/DT(MDEPTH)
C
      WW=W(IJ)
      ISPL=ISPLIN
      IF(ISPLIN.GE.5) THEN
         ISPLIN=ISPL-5
         IF(IJALI(IJ).GT.0) THEN
            IF(IRTE.EQ.0) THEN
               CALL RTEDF1(IJ)
             ELSE 
               CALL RTEDF2(IJ)
            END IF
            ISPLIN=ISPL
            if(ifprad.eq.0) return
            DO ID=1,ND
               if(.not.lskip(ID,IJ))
     *         PRADT(ID)=PRADT(ID)+RAD1(ID)*FAK1(ID)*W(ij)
            END DO
            if(.not.lskip(1,IJ))
     *         PRD0=PRD0+ABSO1(1)*W(IJ)*(RAD1(1)*FH(IJ)-HEXTRD(IJ))
            DO ID=1,ND
               PRADA(ID)=PRADA(ID)+RAD1(ID)*FAK1(ID)*WW
            END DO
            RETURN
         END IF
      END IF
C
      if(icompt.gt.0.and.(iter.gt.1.or.ilam.gt.0)) then
         call rtecf1(ij)
         return
      end if
c
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)
         RAD1(ID)=0.
         ALI1(ID)=0.
      END DO
C
C     non-coherent electron scattering by lambda iteration
C  
      IF(NELSC.LE.0) THEN
         DO ID=1,ND
            SS0(ID)=-SCAT1(ID)/AB0(ID)
         END DO
       ELSE
         DO ID=1,ND
            ST0(ID)=ST0(ID)+SCAT1(ID)*EMEL1(ID)*RAD1(ID)/AB0(ID)
            SS0(ID)=0.
         END DO
      END IF
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.
c     TAUMIN=ABSO1(1)*DEDM1
      taumin=absot(1)*dm(1)/2.
      ALB1=0.
C
C     Allowance for wind blanketing
C
      IF(IWINBL.GT.0) ALB1=TWO*ALBE(IJ)/(UN+ALBE(IJ))
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/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)
         END IF
C
         BI=B/AMU(I)
         CI=C/AMU(I)
         VL(I)=(BI+P0)*ST0(ID)+CI*ST0(ID+1)
         IF(IWINBL.LT.0) VL(I)=VL(I)+EXTINT(IJ,I)
         DO J=1,NMU
            BB(I,J)=SS0(ID)*WTMU(J)*(BI+P0)-ALB1*WTMU(J)
            CC(I,J)=-CI*SS0(ID+1)*WTMU(J)
         END DO
         BB(I,I)=BB(I,I)+AMU(I)/DTP1+UN+BI
         CC(I,I)=CC(I,I)+AMU(I)/DTP1-CI
         ANU(I,ID)=0.
      END DO
C
      IF(ISPLIN.LE.2) THEN
         CALL MATINV(BB,NMU,MMU)
         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
c        CALL MINV3(BB)
         CALL MATINV(BB,NMU,MMU)
         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
      END IF
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)*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
            VL(I)=VL0
            DIV=AMU(I)*AMU(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,MMU)
            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,MMU)
            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/AMU(I)
            AI=A/AMU(I)
            VL(I)=ST0(ID)*BI+ST0(ID-1)*AI
            DO J=1,NMU
               AA(I,J)=-AI*SS0(ID-1)*WTMU(J)
               BB(I,J)=BI*SS0(ID)*WTMU(J)
            END DO
            AA(I,I)=AA(I,I)+AMU(I)/DTP1-AI
            BB(I,I)=BB(I,I)+AMU(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)=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
         ELSE
            DO I=1,NMU
               A=AMU(I)/DTP1
               B=HALF/A
               AA(I,I)=A
               VL(I)=B*ST0(ID)+PLAND+AMU(I)*DPLAN+AA(I,I)*ANU(I,ID-1)
               DO J=1,NMU
                  BB(I,J)=B*SS0(ID)*WTMU(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,MMU)
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 ***************** Backsolution
C
      FKK(ND)=THIRD
      AJ=0.
      AH=0.
      AK=0.
      DO I=1,NMU
         RMU=WTMU(I)*ANU(I,ID)
         AJ=AJ+RMU
         AH=AH+RMU*AMU(I)
         AK=AK+RMU*AMU(I)*AMU(I)
      END DO
      RDD(ID)=AJ
      IF(IBC.EQ.0) THEN
         FKK(ND)=THIRD
       ELSE
         FKK(ID)=AK/AJ
         FHD(IJ)=AH/AJ
      END IF
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
            RMU=WTMU(I)*ANU(I,ID)
            AJ=AJ+RMU
            AK=AK+RMU*AMU(I)*AMU(I)
         END DO
C
C         solution of the transfer equation
C         Variables:
C          ANU     - Feautrier intensity
C          RDD     - mean intensity
C          FKK     - Eddington factor f(K) = K/J
C
         FKK(ID)=AK/AJ
         RDD(ID)=AJ
      END DO
C
      if(idisk.eq.1) then
         do id=1,nd
            fak(ij,id)=fkk(id)
         end do
      end if
C
C     the "surface" Eddington factor fH
C
      AH=0.
      DO I=1,NMU
         AH=AH+WTMU(I)*AMU(I)*ANU(I,1)
      END DO
      FH0=AH/AJ-HALF*ALB1
      FH(IJ)=FH0
      q0(ij)=qq0
      uu0(ij)=u0

c     q0(ij)=0.
c     uu0(ij)=0.
  
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
      if(ilmcor.eq.2) then
         do id=1,nd
            scor(id)=un/(un+ss0(id))
         end do
      else if(ilmcor.eq.3) then
         do id=1,nd
            ss0c(id)=ss0(id)
            st0(id)=st0(id)-ss0(id)*rdd(id)
            ss0(id)=0.
         end do
      end if
C
C     Upper boundary condition
C
      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))
      VLL=ST0(ID)+CQ*ST0(ID+1)
      IF(IWINBL.LT.0) VLL=VLL+HEXTRD(IJ)*BQ
      if(ilmcor.eq.2) then
         bbb(id)=bbb(id)*scor(id)
         ccc(id)=ccc(id)*scor(id)
         vll=vll*scor(id)
      end if
      ZZZ(ID)=UN/BBB(ID)
      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*ST0(ID-1)+C*ST0(ID+1)+(UN-A-C)*ST0(ID)
         if(ilmcor.eq.2) then
            aaa(id)=aaa(id)*scor(id)
            bbb(id)=bbb(id)*scor(id)
            ccc(id)=ccc(id)*scor(id)
            vll=vll*scor(id)
         end if
         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)
         END IF
      END DO
C
C     Lower boundary condition
C
      ID=ND
c
c     stellar atmospheric
c
      IF(IDISK.EQ.0.OR.IFZ0.LT.0) then
      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*TWO*FHD(IJ)+A*FKK(ID)
         AAA(ID)=A*FKK(ID-1)
         VLL=ST0(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
c
c     accretion disk - symmetric boundary
c
      ELSE
         B=TWO/DTP1
         BBB(ID)=FKK(ID)/DTP1*B+UN+SS0(ND)
         AAA(ID)=FKK(ID-1)/DTP1*B
         VLL=ST0(ID)
      END IF
      if(ilmcor.eq.2) then
         aaa(id)=aaa(id)*scor(id)
         bbb(id)=bbb(id)*scor(id)
         vll=vll*scor(id)
      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
      FLUX(IJ)=FH(IJ)*RAD1(1)-half*hextrd(ij)
     *        -(st0(1)-ss0(1)*rad1(1))*q0(ij)

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)/st0(nd-1)
             ali1(nd)=rad1(nd)/st0(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)/st0(nd-1)
           ali1(nd)=rad1(nd)/st0(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
      if(ifalih.gt.0) then
c
c     solution for the individual angles - to get Lambda^star_H
C
      do id=1,nd
         alih1(id)=0.
         alij1(id)=0.
      end do
      nw=nmu
      do i=1,nw
        rmmu(i)=-amu(nw-i+1)
        rmmu(i+nw)=amu(i)
      wmmu(i)=wtmu(nw-i+1)
      wmmu(i+nw)=wtmu(i)
      end do
      do i=1,2*nw
         rwmu(i)=rmmu(i)*wmmu(i)*half
      end do
C
c     --------------------- loop over angles
c
      do i=1,2*nw
         do id=1,nd-1
            dtau(id)=dt(id)/abs(rmmu(i))
         end do
c
c        boundary conditions
c
c        rup=extrad(ij)
         rup=extint(ij,i)
C
C        diffusion approximation for semi-infinite atmospheres
C
         rdown=pland+rmmu(i)*dplan
c
c        solution of the transfer equation
c
         call rtesol(dtau,st0,rup,rdown,rmmu(i),ri,ali)
c
       DO ID=1,ND
            alih1(id)=alih1(id)+rwmu(i)*ali(id)
            alij1(id)=alij1(id)+wmmu(i)*ali(id)*half
       END DO
      end do
      end if
C
c     --------------------- end of loop over angles
c
      ISPLIN=ISPL
c
      if(idisk.ne.0) then
        iji=nfreq-kij(ij)+1
        DO ID=1,ND
           rad(iji,id)=rad1(id)
        END DO
      end if
C
C     radiation pressure
C
      if(ifprad.gt.0) then
      if(.not.lskip(1,IJ))
     *   PRD0=PRD0+ABSO1(1)*WW*(RAD1(1)*FH(IJ)-HEXTRD(IJ))
      DO ID=1,ND
         if(.not.lskip(ID,IJ))
     *      PRADT(ID)=PRADT(ID)+RAD1(ID)*FAK1(ID)*WW
         PRADA(ID)=PRADA(ID)+RAD1(ID)*FAK1(ID)*WW
      END DO
      end if
c
      if(chmax.ge.1.91e-3.and.chmax.le.2.03e-3) then
      tauij=taumin
      do id=1,nd
      if(ilmcor.eq.3) ss0(id)=ss0c(id)
      if(id.gt.1) tauij=tauij+dt(id-1)
      write(97,697) ij,id,tauij,rad1(id),st0(id)/(un+ss0(id)),st0(id),
     *              un+ss0(id),ali1(id)
      end do
  697 format(2i4,1p6e12.4)
      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
c     Q0(IJE)=QQ0
c     UU0(IJE)=U0
      RETURN
      END