SUBROUTINE SOLVES
C     =================
C
C     Same routine than SOLVE, but used for small systems of
C     equations (NN<=MSMX) to keep matrices in memory and to avoid I/O
C
C     Driving procedure for complete linearization
C     Solution of the system
C
C      A * del(PSI{ID-1}) + B * del(PSI{ID}) + C * del(PSI{ID+1}) = VECL
C
C     where PSI{ID} means vector PSI at depth ID,
C           del(PSI{ID}) linearization corrections to PSI;
C     Vector PSI is a vector composed of all unknown model parameters,
C     the choice and order of which is given by input parameters INHE,
C     INRE,INPC,INSE, and INMP (see START).
C     A,B,C are the so-called matrices of complete linearization
C     and VECL is the corressponding rhs vector, pertinent to depth ID.
C
C     The above block tridiagonal system is solved by the standard
C     Gaussian elimination, ie.
C
C       del(PSI{ID}) = BET{ID} - ALF{ID} * del{PSI{ID+1}),
C
C     (the so-called back-solution),
C     where the auxiliary matrix ALF and vector BET, at each depth,
C     are given by (the so-called forward elimination)
C
C       ALF{ID} = (B - A * ALF{ID-1})**-1 * C ,
C     and
C       BET{ID} = (B - A * ALF{ID-1})**-1 * (VECL - A * BET{ID-1})
C
C     Programming notes:
C     Although matrices A,B,C have size (NN x NN), only matrix B is
C     represented as array (MTOT x MTOT), MTOT being the maximum value
C     of NN; other matrices have large parts void, so that
C     A - is represented as array (MAROWS x MTOT), the maximum number of
C         rows of A is MFREX+2 (because only the transfer equations, the
C         hydrostatic equilibrium, and optionally radiative equilibrium
C         equations couple the depth point ID and ID-1)
C     C - is represented as array (MFREX x MCCOLS), ie.
C         the maximum number of rows is MFREX (only the transfer
C         equations couple the depth point ID and ID+1);
C         and since the upper square block (NFREQE x NFREQE), corresponding
C         to the radiative transfer equations, is diagonal, all the
C         transfer equations may be represemted as one column, thus
C         reducing the number of columns of C to number of constraint
C         equations + 1, ie. the maximum number to MCCOLS=MLEVEL+5.
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ITERAT.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ARRAY1.FOR'
      INCLUDE 'ALIPAR.FOR'
      DIMENSION ALF(MTOT,MTOT),BET(MTOT,MDEPTH),DPSI(MTOT)
      COMMON/STOMAT/STOA(MSMX,MSMX,MDEPTH),STOB(MSMX,MSMX,MDEPTH),
     *              STOALF(MSMX,MSMX,MDEPTH)
      COMMON/CMATZD/CZZ,CZN,CZE,CZM
      EQUIVALENCE (DPSI(1),Y2(1))
C
      N=NN
      M=NFREQE
      IF(IFALI.LE.5) THEN
         M1=NFREQE+2
         IF(ICONV.GT.0.OR.NRETC.NE.0) M1=NFREQE+3
       ELSE
         M1=N
      END IF
C
C     --------------------------------
C     First part - forward elimination
C     --------------------------------
C
      LMKA=.FALSE.
      IF(ITER.LT.NITER) LMKA=KANT(ITER+1).EQ.1
      LASO=KANT(ITER).EQ.1
      DO ID=1,ND
C
C        evaluate matrices A,B,C, and rhs vector VECL, corresponding
C        to depth ID
C
         CALL WNSTOR(ID)
         IF(.NOT.LASO) THEN
            CALL MATGEN(ID)
            IF(LMKA) THEN
               DO J=1,N
                  DO I=1,M1
                     STOA(I,J,ID)=A(I,J)
                  END DO
               END DO
            END IF
          ELSE
            CALL RHSGEN(ID)
            DO J=1,N
               DO I=1,M1
                  A(I,J)=STOA(I,J,ID)
               END DO
            END DO
         END IF
         IF(ID.GT.1) THEN
C
C           expression   VECL-A*BET
C
            DO 20 I=1,M1
               SUM=0.
               DO 10 J=1,N
                  SUM=SUM+A(I,J)*BET(J,ID-1)
   10          CONTINUE
               VECL(I)=VECL(I)-SUM
   20       CONTINUE
C
C           expression   B-A*ALF,  stored in array B
C
            IF(.NOT.LASO) THEN
               DO I=1,M1
                  DO J=1,N
                     SUM=0.
                     DO K=1,N
                        SUM=SUM+A(I,K)*ALF(K,J)
                     END DO
                     B(I,J)=B(I,J)-SUM
                  END DO
               END DO
            END IF
         END IF
C
C        expression  (B-A*ALF)**-1,  stored in array B
C
         IF(LASO) THEN
            DO J=1,N
               DO I=1,N
                  B(I,J)=STOB(I,J,ID)
               END DO
            END DO
         ELSE
            CALL MATINV(B,N,MTOT)
            IF(LMKA) THEN
               DO J=1,N
                  DO I=1,N
                     STOB(I,J,ID)=B(I,J)
                  END DO
               END DO
            END IF
         END IF
C
C        auxiliary vector BET = (B-A*ALF)**-1 * (L-A*BET)
C
         DO I=1,N
            SUM=0.
            DO J=1,N
               SUM=SUM+B(I,J)*VECL(J)
            END DO
            BET(I,ID)=SUM
         END DO
C
         IF(ID.LT.ND) THEN
C
C        auxiliary matrix   ALF = (B-A*ALF)**-1 * C
C
            IF(LASO) THEN
               DO J=1,N
                  DO I=1,N
                     ALF(I,J)=STOALF(I,J,ID)
                  END DO
               END DO
            ELSE
               DO I=1,N
                  DO J=1,M
                     ALF(I,J)=B(I,J)*C(J,J)
                  END DO
               END DO
               DO I=1,N
                  DO J=M+1,N
                     SUM=0.
                     DO K=1,M
                        SUM=SUM+B(I,K)*C(K,J)
                     END DO
                     IF(IFALI.GT.5.OR.ICONV.LE.2) THEN
                        DO K=M+1,N
                           SUM=SUM+B(I,K)*C(K,J)
                        END DO
                     END IF
                     ALF(I,J)=SUM
                  END DO
C
C     taking into account the separate elements of matrix C
C
                  BZ=0.
                  IF(INZD.GT.0) BZ=B(I,M+INZD)
                  IF(INZD.GT.0) ALF(I,M+INZD)=ALF(I,M+INZD)+BZ*CZZ
                  IF(INHE.GT.0) ALF(I,M+INHE)=ALF(I,M+INHE)+BZ*CZN
                  IF(INPC.GT.0) ALF(I,M+INPC)=ALF(I,M+INPC)+BZ*CZE
                  IF(INMP.GT.0) ALF(I,M+INMP)=ALF(I,M+INMP)+BZ*CZM
               END DO
C
C           store the auxiliary matrix ALF
C
               DO J=1,N
                  DO I=1,N
                     STOALF(I,J,ID)=ALF(I,J)
                  END DO
               END DO
            END IF
         END IF
      END DO
C
C     --------------------------
C     Second part - backsolution
C     --------------------------
C
      DO I=1,N
         VECL(I)=0.
      END DO
C
      DO IID=1,ND
         ID=ND-IID+1
         IF(ID.LT.ND) THEN
C
C           read old  PSI  
C
            DO I=1,N
               PSI0(I)=PSY0(I,ID)
            END DO
C
C           read auxiliary matrix ALF and vector BET
C
            DO J=1,N
               DO I=1,N
                  ALF(I,J)=STOALF(I,J,ID)
               END DO
            END DO
C
C           expression  ALF * delta(PSI{previous depth}), stored 
C           in array VECL
C
            DO I=1,N
               SUM=0.
               DO J=1,N
                  SUM=SUM+ALF(I,J)*DPSI(J)
               END DO
               VECL(I)=SUM
            END DO
         END IF
C
C        finally, evaluate delta(PSI) and corresponding relative changes
C        delta(PSI)/PSI
C
         DO I=1,N
            DPSI(I)=BET(I,ID)-VECL(I)
            CHAN=0.
            IF(PSI0(I).GT.0.) CHAN=DPSI(I)/PSI0(I)
            BET(I,ID)=CHAN
C
C           over-relaxation
C
            IF(I.GE.NFREQE+INSE) CHAN=ORELAX*CHAN 
C
C           To avoid instabilities, relative changes of all quantities
C           are artificially limited not to exceed certain predefined values
C
            IF(CHAN.LE.UN/DPSILG-UN) CHAN=UN/DPSILG-UN
            IF(CHAN.GE.DPSILG-UN) CHAN=DPSILG-UN
            IF(INRE.GT.0) THEN
               DPLP=DPSILT-UN
               DPLM=UN/DPSILT-UN
               IF(CHAN.LE.DPLM.AND.I.EQ.NFREQE+INRE) CHAN=DPLM
               IF(CHAN.GT.DPLP.AND.I.EQ.NFREQE+INRE) CHAN=DPLP
            END IF
            IF(INHE.GT.0) THEN
               DPLP=DPSILN-UN
               DPLM=UN/DPSILN-UN
               IF(CHAN.LE.DPLM.AND.I.EQ.NFREQE+INHE) CHAN=DPLM
               IF(CHAN.GT.DPLP.AND.I.EQ.NFREQE+INHE) CHAN=DPLP
            END IF
            IF(INPC.GT.0) THEN
               DPLP=DPSILN-UN
               DPLM=UN/DPSILN-UN
               IF(CHAN.LE.DPLM.AND.I.EQ.NFREQE+INPC) CHAN=DPLM
               IF(CHAN.GT.DPLP.AND.I.EQ.NFREQE+INPC) CHAN=DPLP
            END IF
            IF(INDL.GT.0) THEN
               DPLP=DPSILD-UN
               DPLM=UN/DPSILD-UN
               IF(CHAN.LE.DPLM.AND.I.EQ.NFREQE+INDL) CHAN=DPLM
               IF(CHAN.GT.DPLP.AND.I.EQ.NFREQE+INDL) CHAN=DPLP
            END IF
C
C           new vector PSI
C
            PSI0(I)=PSI0(I)*(CHAN+UN)
         END DO
         IF(INRE.GT.0) CHANT(ID)=BET(NFREQE+INRE,ID)
C
C        new vector PSI is stored
C
         DO I=1,N
            PSY0(I,ID)=PSI0(I)
         END DO
      END DO
C
C     print out the relative changes of vector PSI
C
      CALL PRCHAN(BET,CHMX,CHMT)
C
C     STOP if changes become too large
C
      IF(ITER.NE.1 .AND. ABS(CHMX).GT.1.D16) THEN
         WRITE(6,610) ITER,CHMX
         WRITE(10,610) ITER,CHMX
         STOP
      END IF
      
C     Reset iron lines cross-sections if changes are still large

      if(iter.le.1) LIROST=.FALSE.
      LITEK=.FALSE.
      if(iter.eq.7 .or. iter.eq.11 .or. iter.eq.15) LITEK=.TRUE.
      if(chmt.gt.chmaxt .and. ispodf.ge.1) LIROST=.TRUE.
      if(LITEK.and.LIROST) then
         call iroset
         LIROST=.FALSE.
      end if

C     Set up Kantorovich method
 
      IF(LASO) THEN
         WRITE(6,600) ITER
         WRITE(10,600) ITER
      END IF
C
C     Finally, set up quantity LFIN that indicates whether or not
C     this iteration of complete linearization is the last one
C
      LFIN=ABS(CHMX).LE.CHMAX.OR.ITER.GE.NITER
C
  600 FORMAT(' **** KANTOROVICH acceleration: ITER',I4)
  610 FORMAT(' **** STOP in SOLVE after ITER',I4,/,
     * ' Max change:',1PE12.2)
      RETURN
      END