SUBROUTINE RYBSOL
C     =================
C
C     driver of the complete-linearization solution in the Rybicki 
C     formalism
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ARRAY1.FOR'
      INCLUDE 'ITERAT.FOR'
      COMMON/RYBMTX/RA(MDEPTH),RB(MDEPTH),RC(MDEPTH),VR(MDEPTH),
     *              UA(MDEPTH),UB(MDEPTH),UC(MDEPTH),    
     *              VA(MDEPTH),VB(MDEPTH),VC(MDEPTH),WR(MDEPTH), 
     *              WM(MDEPTH,MDEPTH)
      common/imodlc/imodl0(mlevel)
      DIMENSION AL(MDEPTH),AU(MDEPTH),VAL(MDEPTH),UCOL(MDEPTH),
     *          VAU(MDEPTH,MDEPTH),CHANGT(MDEPTH)
      dimension pop1(mlevel),babs(mlevel,mdepth)
C
C     zeroing needed quantities 
C
      DO ID=1,ND
         REIT(ID)=0.
         REIN(ID)=0.
         AREIT(ID)=0.
         AREIN(ID)=0.
         CREIT(ID)=0.
         CREIN(ID)=0.
         REDT(ID)=0.
         REDTM(ID)=0.
         REDTP(ID)=0.
         REDN(ID)=0.
         REDNM(ID)=0.
         REDNP(ID)=0.
         FCOOLI(ID)=0.
         FLFIX(ID)=0.
         FLRD(ID)=0.
         ABROSD(ID)=0.
         SUMDPL(ID)=0.
      END DO
C
C     zero vectors and matrices of the Rybicki formalism
C
      DO ID=1,ND
         RA(ID)=0.
         RB(ID)=0.
         RC(ID)=0.
         UA(ID)=0.
         UB(ID)=0.
         UC(ID)=0.
         VA(ID)=0.
         VB(ID)=0.
         VC(ID)=0.
         VR(ID)=0.
         WR(ID)=0.
         AL(ID)=0.
         AU(ID)=0.
         VAL(ID)=0.
         DO ID1=1,ND
            WM(ID1,ID)=0.
            VAU(ID1,ID)=0.
         END DO
      END DO
C
      if(ioptab.lt.-1) call setdrt
c
      DO IJ=1,NFREQ
         FR=FREQ(IJ)
         W0=W0E(IJ)
c        IF(IOPTAB.GE.0) THEN
            CALL OPACTR(IJ)
c         ELSE
c           CALL OPACFD(IJ)
c        END IF
         CALL RTEFR1(IJ)
         CALL ALIFR1(IJ)
         CALL ROSSTD(IJ)
C
         CALL RYBMAT(IJ)
         CALL TRIDAG(RA,RB,RC,VR,AL,ND)
         ID=1
         VAL(ID)=VAL(ID)+VB(ID)*AL(ID)+VC(ID)*AL(ID+1)
         DO ID=2,ND-1
            VAL(ID)=VAL(ID)+VA(ID)*AL(ID-1)+VB(ID)*AL(ID)+
     *                      VC(ID)*AL(ID+1)
         END DO
         ID=ND
         VAL(ID)=VAL(ID)+VA(ID)*AL(ID-1)+VB(ID)*AL(ID)
C
         DO IDC=1,ND
            DO ID=1,ND
               UCOL(ID)=0.
            END DO
            UCOL(IDC)=UB(IDC)
            IF(IDC.GT.1) UCOL(IDC-1)=UC(IDC-1)
            IF(IDC.LT.ND) UCOL(IDC+1)=UA(IDC+1)
            CALL TRIDAG(RA,RB,RC,UCOL,AU,ND)
            ID=1
            VAU(ID,IDC)=VAU(ID,IDC)+VB(ID)*AU(ID)+VC(ID)*AU(ID+1)
            DO ID=2,ND-1
               VAU(ID,IDC)=VAU(ID,IDC)+
     *         VA(ID)*AU(ID-1)+VB(ID)*AU(ID)+VC(ID)*AU(ID+1)
            END DO
            ID=ND
            VAU(ID,IDC)=VAU(ID,IDC)+VA(ID)*AU(ID-1)+VB(ID)*AU(ID)
         END DO
      END DO
C
      DO ID=1,ND
         ABROSD(ID)= SUMDPL(ID)/ABROSD(ID)
         FCOOL(ID)=REINT(ID)*FCOOLI(ID)*DENS(ID)-REDIF(ID)*FLFIX(ID)
      END DO
      if(ioptab.lt.0) CALL ROSSTD(0)
C
C     final evaluation of matrices and the global inversion
C    
      CALL RYBENE
C
       do id=1,70,34
          write(6,603) id,wm(id,id),wm(id,id+1),wr(id)
       end do
  603  format('rybene ',8x,i4,1p4e11.3)
      DO ID=1,ND
         DO IDC=1,ND
            WM(ID,IDC)=WM(ID,IDC)-VAU(ID,IDC)
         END DO
         WR(ID)=WR(ID)-VAL(ID)
      END DO
      CALL LINEQS(WM,WR,CHANGT,ND,MDEPTH)
C
      IF(.NOT.LTE.and.ifryb.gt.1) THEN
         lte=.true.
         iflev0=iflev
         iflev=1
         do i=1,nlevel
            imodl(i)=imodl0(i)
         end do
         call levset
         do id=1,nd
            call steqeq(id,pop1,0)
            do i=1,nlevel
               babs(i,id)=un
               if(pop1(i).gt.0.) babs(i,id)=popul(i,id)/pop1(i)
            end do
         end do
         lte=.false.
         iflev=iflev0
         do i=1,nlevel
            imodl(i)=imodl0(i)
         end do
         call levset
      end if
C
      IF(.NOT.LTE.and.ifryb.gt.2) THEN
      DO IJ=1,NFREQ
         CALL OPACTR(IJ)
         CALL RTEFR1(IJ)
         CALL ALIFR1(IJ)
         CALL ROSSTD(IJ)
         CALL RYBMAT(IJ)
         ID=1
         WR(ID)=VR(ID)-UB(ID)*CHANGT(ID)-UC(ID)*CHANGT(ID+1)
         DO ID=2,ND-1
            WR(ID)=VR(ID)-UA(ID)*CHANGT(ID-1)-UB(ID)*CHANGT(ID)-
     *             UC(ID)*CHANGT(ID+1)
         END DO
         ID=ND
         WR(ID)=VR(ID)-UA(ID)*CHANGT(ID-1)-UB(ID)*CHANGT(ID)
         CALL TRIDAG(RA,RB,RC,WR,AL,ND)
         RAD(IJ,ID)=RAD1(ID)+AL(ID)
      END DO
      END IF
C
      CALL RYBCHN(CHANGT)
c
      if(.not.lte.and.ifryb.gt.1) then
         lte=.true.
         iflev0=iflev
         iflev=1
         do i=1,nlevel
            imodl(i)=imodl0(i)
         end do
         call levset
         do id=1,nd
            call steqeq(id,pop1,0)
            do i=1,nlevel
               popul(i,id)=pop1(i)*babs(i,id)
            end do
         end do
         lte=.false. 
         iflev=iflev0
         do i=1,nlevel
            imodl(i)=imodl0(i)
         end do
         call levset
      end if
c
      RETURN
      END