      SUBROUTINE BMSOL(W,B,NBAND,NROW)      
*     Solution of band matrix equation W * X = B
*     with symmetric band matrix in special storage mode in array W.
*     Array B(NROW) is right-hand side, replaced by result.
      INTEGER I,J,N,NN,IMAX,JMAX,NDIM,NBAND,NROW
      PARAMETER (NDIM=100)
      REAL W(NBAND,NROW),B(NROW),DIAG(NDIM),RATIO
      DO N=1,NROW
       DIAG(N)=W(1,N)
      END DO
*     decomposition
      DO N=1,NROW
       IF(W(1,N)+DIAG(N).LE.DIAG(N)) THEN
*         singular case
          DO J=1,NBAND
           W(J,N)=0.0
          END DO
       ELSE
          W(1,N)=1.0/W(1,N)
          IMAX=MIN(NBAND-1,NROW-N)
          JMAX=IMAX
          DO I=1,IMAX
           RATIO=W(I+1,N)*W(1,N)
           DO J=1,JMAX
            W(J,N+I)=W(J,N+I)-W(J+I,N)*RATIO
           END DO
           JMAX=JMAX-1
           W(I+1,N)=RATIO
          END DO
       END IF
      END DO
*     solution by forward ...
      DO N=1,NROW
       DO J=1,MIN(NBAND-1,NROW-N)
        B(J+N)=B(J+N)-W(J+1,N)*B(N)
       END DO
      END DO
*     ... and backward substitution
      DO NN=1,NROW
       N=NROW+1-NN
       B(N)=B(N)*W(1,N)
       DO J=1,MIN(NBAND-1,NROW-N)
        B(N)=B(N)-W(J+1,N)*B(J+N)
       END DO
      END DO
      END
