CSTART OF AS 13
      SUBROUTINE MSTREE(N,A,Dlarge,D,B,C,Ifault)
      IMPLICIT NONE
*
*     On input:
*        N = number of points
*        A is a logical array of dimension (2:N) for use by MSTREE
*        Dlarge is larger than the largest distance between points
*        D is a triangular matrix of distances between pairs of points
*          fill D by
*             nd = 0
*             DO j = 2, N
*               DO i = 1, j-1
*                 nd = nd + 1
*                 D(nd) = distance between points i and j
*               ENDDO
*             ENDDO
*        B is an integer array dimensioned (2:N). 
*          On return B(i) is the number of a point to which point i is
*          connected.
*        C is a real array dimensioned (2:N).
*          On return C(i) is the distance from point i to point B(i)
*        Ifault is returned non-zero for certain errors
*
*
*     modified from:                --  WJM  29 Nov 1998
*
*      SUBROUTINE PRTREE(N, M, A, DLARGE, D, B, C, IFAULT)
C
C        ALGORITHM AS 13  APPL. STATIST. (1969) VOL.18, P.103
C
C        COMPUTES THE MINIMUM SPANNING TREE OF A DISTANCE MATRIX D
C
      INTEGER N, B(2:N), Ifault
      LOGICAL A(2:N)
      REAL    Dlarge, D(*), C(2:N)
*      REAL D(n*(n-1)/2)
*      REAL D(M)
      INTEGER i, j, k, l, next
      REAL    am, dist
*
      Ifault = 1
      IF ( N.LT.2 ) RETURN
*      IFAULT = 2
*      IF (N * (N - 1) / 2 .NE. M) RETURN
      Ifault = 0
C
C        A(I) IS .FALSE. IF I IS ALREADY ASSIGNED TO
C        THE TREE (INITIALLY CONSISTING OF NO. 1 ONLY),
C        OR .TRUE. OTHERWISE
C
!      PRINT *, 'in MSTREE', N, Dlarge
      DO i = 2, N
         A(i) = .TRUE.
         B(i) = 0
         C(i) = Dlarge
      ENDDO
      j = 1
      DO i = 2, N
         am = Dlarge
         DO k = 2, N
            IF ( A(k) ) THEN
               IF ( j.GT.k ) l = (j-1)*(j-2)/2 + k
               IF ( j.LE.k ) l = (k-1)*(k-2)/2 + j
               dist = D(l)
               IF ( dist.LT.C(k) ) THEN
                  C(k) = dist
                  B(k) = j
               ENDIF
               IF ( am.GT.C(k) ) THEN
                  am = C(k)
                  next = k
               ENDIF
            ENDIF
         ENDDO
         j = next
         A(j) = .FALSE.
      ENDDO
      RETURN
      END
CEND OF AS 13
