**==aa0001.spg  processed by SPAG 5.11R  at 11:51 on  8 Jan 1999
c
c program implementing Friedman-Rafsky two-sample tests.
c
c Friedman, J. H. and Rafsky, L. C. (1979). Multivariate analogs of
c    the Wald-Wolfowitz and Smirnov two-sample tests. Annals of
c    Statistics, 7, 697-717 (June).
c
c version of 10/17/79.
c
c coded and copyright (c) by:
c
c                        Jerome H. Friedman
c
c                     Department of Statistics
c                               and
c                  Stanford Linear Accelerator Center
c                        Stanford University
c
! modified structure and names somewhat  -  W.J.Metzger - 8 Jan 1999
c
! wjm c call mstcpr(p,n,n1,data,ntree,deg,mst)
! wjm c integer p,n,n1,ntree
c
! call mstclfrd(n1,n2,dlarge,dist,ntree,deg,mst,signif)
! integer p,n1,n2,ntree
! real dist(n*(n-1)/2),deg(n,2)
c
c integer mst(max0(5,ntree+3),n)
c
! wjm c   p = dimension.
! wjm c   n = cardinality of pooled sample.
c   n1 = cardinality of first sample.
c   n2 = cardinality of second sample.                ! WJM
c   ntree = number of mst's comprising graph.
! wjm c   data = pooled sample data matrix.
! wjm c          first n1 rows correspond to first sample.
! wjm c          rows n1+1 through n correspond to second sample.
! wjm  dist = triangular distance matrix:
! wjm           dist( (Max(i,j)-1)*(Max(i,j)-2)/2 + Min(i,j) )
! wjm         is the distance from point i to point j
! wjm  dlarge = a number larger than the largest distance, dist
c   deg = scratch storage for mstcpr. upon return
c         deg(1,1) - deg(n,1) contain degree of each node (point)
c         in the first mst.
c         deg(1,2) - deg(n,2) contain degree of each node in
c         the ntree-th mst.
c   mst = scratch storage for mstcpr. upon return
c         mst(k,1) - mst(k,n-1) = kth mst (1.le.k.le.ntree).
c
c labeled commons:
c
! wjm c common /answer/ signif(ltree+3)
c real signif(ntree+3)
c
c   signif(1) = significance "standard" multivariate
c                      smirnov test.
c   signif(2) = significance "radial" multivariate
c                      smirnov test.
c   signif(3) = significance sample 1 (degree 1) test.
c   signif(4) = significance runs test 1st mst.
c   signif(5) = significance runs test 1st and 2nd mst.
c   signif(k+3) = significance runs test union(1-k) mst.
c   signif(ntree+3) = significance runs test union(1-ntree) mst.
c
! wjm c   signif(1) = significance sample 1 (degree 1) test.
! wjm c   signif(2) = significance runs test 1st mst.
! wjm c   signif(3) = significance runs test 1st and 2nd mst.
! wjm c   signif(k+1) = significance runs test union(1-k) mst.
! wjm c   signif(ntree+1) = significance runs test union(1-ntree) mst.
! wjm c   signif(ntree+2) = significance "standard" multivariate
! wjm c                      smirnov test.
! wjm c   signif(ntree+3) = significance "radial" multivariate
! wjm c                      smirnov test.
c
c common /print/ iprint
c
c   iprint = print flag for program.
c          = 0 : no printed output.
c          .ne. 0: printed output.
c   default value = 1.
c
      BLOCKDATA 
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      INTEGER IPRint
C*** End of declarations inserted by SPAG
      COMMON /MSTCLFPR/ IPRint
      DATA IPRint/0/
      END
      SUBROUTINE MSTCLFRD(N1,N2,Dlarge,Dis,Lt,Deg,Istor,Signif)
! wjm SUBROUTINE MSTCPR(D,N,N1,X,Lt,Deg,Istor)
      IMPLICIT NONE
      REAL Dlarge, Signif(*)
! wjm REAL X(D,N), Deg(N,2)
      REAL Dis(*), Deg(N1+N2,2)
      INTEGER N2                     ! wjm
! wjm INTEGER D
      INTEGER N, N1, Istor(*)
      INTEGER IPRint
      COMMON /MSTCLFPR/ IPRint
C*** Start of declarations inserted by SPAG
      REAL e11, edges, ERF, eruns, p11, qshr, runs, std, std11, v11
      INTEGER ier, j, k, n1p1, nd1, nd11, nm1
C*** End of declarations inserted by SPAG
      INTEGER ft, Lt
      INTEGER ttn, ftn, node
      REAL*8 peqk, plek
      REAL*8 fn1, fn2, fn, p, var1, covar1, covar2, tqshr, var
      N = N1 + N2                    ! wjm
      ft = 1
      nm1 = N - 1
      ttn = MAX0(Lt*nm1,N+N)
      ftn = ttn + N + N
! wjm CALL PRIMK(D,N,X,ft,Lt,Istor,Deg,Istor(ttn+1),Istor(ftn+1))
      CALL PRIMK(N,Dis,Dlarge,ft,Lt,Istor,Deg,Istor(ttn+1),Istor(ftn+1))
      DO j = 1, nm1
         Deg(j,1) = 1.0
      ENDDO
      Deg(N,1) = 0.0
      DO j = 1, nm1
         Deg(Istor(Lt*(j-1)+1),1) = Deg(Istor(Lt*(j-1)+1),1) + 1.0
      ENDDO
      DO j = 1, N
         Deg(j,2) = Deg(j,1)
      ENDDO
      runs = 1.0
      nd1 = 0
      nd11 = nd1
      n1p1 = N1 + 1
      fn = N
      fn1 = N1
      fn2 = N - N1
      p = 2.0*fn1*fn2/(fn*(fn-1.0))
      var1 = p*(1.0-p)
      covar1 = 0.5*p - p**2
      covar2 = 2.0*p*(fn1-1.0)*(fn2-1.0)/((fn-2.0)*(fn-3.0)) - p**2
      DO k = 1, Lt
         IF ( k.EQ.1 ) THEN
            DO j = 1, N
               IF ( Deg(j,1).EQ.1 ) THEN
                  nd1 = nd1 + 1
                  IF ( j.LE.N1 ) nd11 = nd11 + 1
               ENDIF
            ENDDO
         ELSE
            DO j = 1, nm1
               Deg(j,2) = Deg(j,2) + 1.0
               Deg(Istor(Lt*(j-1)+k),2) = Deg(Istor(Lt*(j-1)+k),2) + 1.0
            ENDDO
         ENDIF
         DO j = 1, N1
            IF ( Istor(Lt*(j-1)+k).GT.N1 ) runs = runs + 1.0
         ENDDO
         DO j = n1p1, nm1
            IF ( Istor(Lt*(j-1)+k).LE.N1 ) runs = runs + 1.0
         ENDDO
         tqshr = 0.0
         DO j = 1, N
            IF ( Deg(j,2).GT.1.0 ) tqshr = tqshr + Deg(j,2)
     &           *(Deg(j,2)-1.0)
         ENDDO
         IF ( k.EQ.1 ) THEN
            p11 = fn1/fn
            e11 = p11*nd1
            v11 = e11*(1.0-p11)*(fn-nd1)/(fn-1.0)
            std11 = (nd11-e11)/SQRT(v11)
            CALL MDHYP(MIN0(nd11,nd1-nd11),nd1,N,N1,peqk,plek,ier)
! wjm       Signif(1) = 2.0*plek - peqk
            Signif(3) = 2.0*plek - peqk
            IF ( IPRint.NE.0 ) WRITE (6,99001) nd1, nd11, e11, std11, 
     &                                Signif(3)
! wjm&                                Signif(1)
         ENDIF
         IF ( k.GE.ft ) THEN
            edges = k*nm1
            eruns = p*edges + 1.0
            var = var1*edges + covar1*tqshr + 
     &            covar2*(edges*(edges-1.0)-tqshr)
            std = (runs-eruns)/DSQRT(var)
! wjm       Signif(k+1) = 0.5*(1.0+ERF(std/1.41421))
            Signif(k+3) = 0.5*(1.0+ERF(std/1.41421))
            IF ( IPRint.NE.0 ) THEN
               qshr = 0.5*tqshr/(N-2)
               WRITE (6,99002) k, qshr
               WRITE (6,99003) runs, eruns, std, Signif(k+3)
! wjm          WRITE (6,99003) runs, eruns, std, Signif(k+1)
            ENDIF
         ENDIF
      ENDDO
      IF ( ft.GT.1 ) RETURN
      Istor(ftn+1) = 0
      DO j = 2, N
         Istor(ftn+j) = Istor(ftn+j-1) + Deg(j-1,1)
      ENDDO
      DO j = 1, nm1
         Istor(ftn+j) = Istor(ftn+j) + 1
         Istor(ttn+Istor(ftn+j)) = Istor(Lt*(j-1)+1)
         node = Istor(Lt*(j-1)+1)
         Istor(ftn+node) = Istor(ftn+node) + 1
         Istor(ttn+Istor(ftn+node)) = j
      ENDDO
      Istor(ftn+1) = 0
      DO j = 2, N
         Istor(ftn+j) = Istor(ftn+j-1) + Deg(j-1,1)
      ENDDO
! wjm CALL SMRNOV(X,D,N,N1,Deg,Istor(ftn+1),Istor(ttn+1),Istor,Deg(1,2),
! wjm&            Lt+2)
      CALL SMRNOV(Dis,N,N1,Deg,Istor(ftn+1),Istor(ttn+1),Istor,
     &            Deg(1,2),Signif)
      RETURN
99001 FORMAT (/,' ',i5,'   deg=1 nodes.  class 1 =',i5,' expect =',f5.1,
     &        /,'    ',g13.5,' stds with hypergeometric significance ',
     &        g13.5)
99002 FORMAT (/,' ',i2,
     &        'th mst. (number of mst edges that share nodes)/(n-2) = ',
     &        g13.5)
99003 FORMAT (/,' runs =',f5.0,'  expect =',f5.0,/,'   ',g13.5,
     &        ' stds with normal significance ',g13.5)
      END
**==primk.spg  processed by SPAG 5.11R  at 11:51 on  8 Jan 1999
      SUBROUTINE PRIMK(N,Dis,Dlarge,Ftree,Ltree,Mst,Ui,Ji,Nit)
! wjm SUBROUTINE PRIMK(M,N,Points,Ftree,Ltree,Mst,Ui,Ji,Nit)
      IMPLICIT NONE
      INTEGER N, Ftree, Ltree, Mst(Ltree,N), Ji(N), Nit(N)
      REAL Dis(*), Dlarge, Ui(N)
C*** Start of declarations inserted by SPAG
      REAL d, uk
      INTEGER i, itree, k, kp, ktree, ktrem1, ni, nitp
C*** End of declarations inserted by SPAG
      DO ktree = Ftree, Ltree
         ktrem1 = ktree - 1
         nitp = N - 1
         kp = N
         DO i = 1, nitp
            Nit(i) = i
            Ui(i) = Dlarge
! wjm       Ui(i) = 1.0E73
         ENDDO
 50      DO i = 1, nitp
            ni = Nit(i)
            IF ( ktrem1.GT.0 ) THEN
               DO itree = 1, ktrem1
                  IF ( Mst(itree,ni).EQ.kp ) GOTO 60
                  IF ( Mst(itree,kp).EQ.ni ) GOTO 60
               ENDDO
            ENDIF
! wjm       d = 0.0
! wjm       DO i001 = 1, M
! wjm          d = d + (Points(i001,ni)-Points(i001,kp))**2
! wjm       ENDDO
            d = Dis((MAX(ni,kp)-1)*(MAX(ni,kp)-2)/2+MIN(ni,kp))     !wjm
            GOTO 80
! wjm 60         d = 1.0E73
 60         d = Dlarge
 80         IF ( Ui(i).GT.d ) THEN
               Ui(i) = d
               Ji(i) = kp
            ENDIF
         ENDDO
         uk = Ui(1)
         DO i = 1, nitp
            IF ( Ui(i).LE.uk ) THEN
               uk = Ui(i)
               k = i
            ENDIF
         ENDDO
         kp = Nit(k)
         Mst(ktree,kp) = Ji(k)
         Ui(k) = Ui(nitp)
         Nit(k) = Nit(nitp)
         Ji(k) = Ji(nitp)
         nitp = nitp - 1
         IF ( nitp.NE.0 ) GOTO 50
      ENDDO
      RETURN
      END
**==smrnov.spg  processed by SPAG 5.11R  at 11:51 on  8 Jan 1999
      SUBROUTINE SMRNOV(Dis,N,N1,Deg,Pntr,Auux,Stack,Depth,Sig)
! wjm SUBROUTINE SMRNOV(X,D,N,N1,Deg,Pntr,Auux,Stack,Depth,Place)
      IMPLICIT NONE
! wjm INTEGER D
      INTEGER N1, N, maxdph, enddia, center, DIAEND
      REAL Dis(*), Deg(*), Depth(*)
      REAL Sig(2)
      INTEGER Pntr(*), Auux(*), Stack(2,1)
      maxdph = N + 1
      CALL MAXTAG(1,Deg,Pntr,Auux,Stack,N,Depth)
      enddia = DIAEND(1,Deg,Pntr,Auux,Depth,maxdph)
      CALL MAXTAG(enddia,Deg,Pntr,Auux,Stack,N,Depth)
! wjm CALL STDSMR(enddia,Deg,Pntr,Auux,Depth,N1,N,Stack,Place)
      CALL STDSMR(enddia,Deg,Pntr,Auux,Depth,N1,N,Stack,Sig(1))
      maxdph = Depth(enddia)/2.0 + 1.0
      center = DIAEND(enddia,Deg,Pntr,Auux,Depth,maxdph)
! wjm CALL RADSMR(X,D,N,N1,center,Deg,Pntr,Auux,Stack,Depth,Place+1)
      CALL RADSMR(Dis,N,N1,center,Deg,Pntr,Auux,Stack,Depth,Sig(2))
      RETURN
      END
**==stdsmr.spg  processed by SPAG 5.11R  at 11:51 on  8 Jan 1999
      SUBROUTINE STDSMR(Root,Deg,Pntr,Auux,Depth,N1,N,Stack,Sig)
! wjm SUBROUTINE STDSMR(Root,Deg,Pntr,Auux,Depth,N1,N,Stack,Place)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL dphmin, f1, f2, Sig
      INTEGER i, ier, IPRint, jdeg
C*** End of declarations inserted by SPAG
      COMMON /MSTCLFPR/ IPRint
      INTEGER Root, N1, N, node, stkptr, parent, next, son
      REAL sum, Depth(1), Deg(1), distks
      INTEGER Pntr(1), Auux(1), Stack(1)
      f1 = -1.0/N1
      f2 = 1.0/(N-N1)
      node = Root
      parent = 0
      stkptr = parent
      Deg(node) = -Deg(node)
      sum = f1
      IF ( node.GT.N1 ) sum = f2
      distks = ABS(sum)
 100  jdeg = ABS(Deg(node))
      next = 0
      dphmin = 9.9E73
      DO i = 1, jdeg
         son = Auux(Pntr(node)+i)
         IF ( Deg(son).GT.0 ) THEN
            IF ( Depth(son).LT.dphmin ) THEN
               dphmin = Depth(son)
               next = son
            ENDIF
         ENDIF
      ENDDO
      IF ( next.NE.0 ) THEN
         stkptr = stkptr + 1
         parent = node
         Stack(stkptr) = parent
         node = next
         IF ( node.LE.N1 ) THEN
            sum = sum + f1
         ELSE
            sum = sum + f2
         ENDIF
         distks = AMAX1(distks,ABS(sum))
         Deg(node) = -Deg(node)
         GOTO 100
      ELSEIF ( stkptr.EQ.0 ) THEN
! wjm    CALL KSNULL(distks,MAX0(N1,N-N1),MIN0(N1,N-N1),p1side,p2side,
! wjm&               ier)
         CALL KSNULL(distks,MAX0(N1,N-N1),MIN0(N1,N-N1),Sig,ier)
         IF ( IPRint.NE.0 ) WRITE (6,99001) distks, Sig
      ELSE
         node = Stack(stkptr)
         parent = 0
         IF ( stkptr.GT.1 ) parent = Stack(stkptr-1)
         stkptr = stkptr - 1
         GOTO 100
      ENDIF
      RETURN
99001 FORMAT (/,' mst standard kolmogorov-smirnov distance = ',g13.5,/,
     &        ' with significance = ',g13.5)
      END
**==radsmr.spg  processed by SPAG 5.11R  at 11:51 on  8 Jan 1999
      SUBROUTINE RADSMR(Dis,N,N1,Start,Deg,Pntr,Auux,Stack,Depth,Sig)
! wjm SUBROUTINE RADSMR(X,D,N,N1,Start,Deg,Pntr,Auux,Stack,Depth,Place)
      IMPLICIT NONE
      REAL Dis(*), Deg(*), Depth(*), Sig
! wjm INTEGER D
      INTEGER Start, N1, N, Pntr(1), Auux(1), Stack(N,1)
      INTEGER IPRint
      COMMON /MSTCLFPR/ IPRint
C*** Start of declarations inserted by SPAG
      INTEGER i, ier, j, jdeg, kj, m, m0
C*** End of declarations inserted by SPAG
! wjm INTEGER node, stkptr, parent, next, istart, Place
      INTEGER node, stkptr, parent, next, istart
      REAL distks, f1, f2, sum
      node = Start
      parent = 0
      istart = 1
      stkptr = 0
 100  jdeg = ABS(Deg(node))
      i = istart
 200  IF ( (i).GT.(jdeg) ) THEN
         Depth(node) = stkptr
         IF ( stkptr.EQ.0 ) THEN
            DO i = 1, N
               Stack(i,1) = i
            ENDDO
! wjm       CALL SORT(Depth,Stack,1,N)
            CALL CLMSSORT(Depth,Stack,1,N)
            m = 0
         ELSE
            node = Stack(stkptr,1)
            parent = 0
            IF ( stkptr.GT.1 ) parent = Stack(stkptr-1,1)
            istart = Stack(stkptr,2)
            stkptr = stkptr - 1
            GOTO 100
         ENDIF
      ELSE
         next = Auux(Pntr(node)+i)
         IF ( next.NE.parent ) THEN
            stkptr = stkptr + 1
            Stack(stkptr,1) = node
            Stack(stkptr,2) = i + 1
            parent = node
            node = next
            istart = 1
            GOTO 100
         ELSE
            i = i + 1
            GOTO 200
         ENDIF
      ENDIF
 300  m = m + 1
      m0 = m
 400  IF ( Depth(m).EQ.Depth(m+1) ) THEN
         m = m + 1
         IF ( m.LT.N ) GOTO 400
      ENDIF
      IF ( m.GT.m0 ) THEN
         DO j = m0, m
            kj = Stack(j,1)
! wjm       Depth(j) = 0
! wjm       DO i = 1, D
! wjm          Depth(j) = Depth(j) + (X(i,kj)-X(i,Start))**2
! wjm       ENDDO
            Depth(j) = Dis((MAX(Start,kj)-1)*(MAX(Start,kj)-2)
     &                 /2+MIN(Start,kj))                          ! wjm
         ENDDO
! wjm    CALL SORT(Depth,Stack,m0,m)
         CALL CLMSSORT(Depth,Stack,m0,m)
      ENDIF
      IF ( m.LT.N ) GOTO 300
      distks = 0.0
      sum = distks
      f1 = -1.0/N1
      f2 = 1.0/(N-N1)
      DO j = 1, N
         IF ( Stack(j,1).GT.N1 ) THEN
            sum = sum + f2
         ELSE
            sum = sum + f1
         ENDIF
         distks = AMAX1(distks,ABS(sum))
      ENDDO
! wjm CALL KSNULL(distks,MAX0(N1,N-N1),MIN0(N1,N-N1),p1side,p2side,ier)
      CALL KSNULL(distks,MAX0(N1,N-N1),MIN0(N1,N-N1),Sig,ier)
      IF ( IPRint.NE.0 ) WRITE (6,99001) distks, Sig
! wjm IF ( IPRint.NE.0 ) WRITE (6,99001) distks, p2side
! wjm signif(Place) = p2side
      RETURN
99001 FORMAT (/,' mst radial  kolmogorov-smirnov distance = ',g13.5,/,
     &        ' with significance = ',g13.5)
      END
**==maxtag.spg  processed by SPAG 5.11R  at 11:51 on  8 Jan 1999
      SUBROUTINE MAXTAG(Start,Deg,Pntr,Auux,Stack,N,Count)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      INTEGER i, jdeg
C*** End of declarations inserted by SPAG
      INTEGER Start, N, node, stkptr, parent, next, istart
      REAL Count(1), Deg(1), maxcnt
      INTEGER Pntr(1), Auux(1), Stack(2,1)
      DO i = 1, N
         Count(i) = 0.0
      ENDDO
      node = Start
      parent = 0
      istart = 1
      stkptr = 0
 100  jdeg = ABS(Deg(node))
      i = istart
 200  IF ( (i).GT.(jdeg) ) THEN
         maxcnt = 0.0
         i = 1
 250     IF ( (i).GT.(jdeg) ) THEN
            Count(node) = maxcnt + 1.0
            IF ( stkptr.EQ.0 ) RETURN
            node = Stack(1,stkptr)
            parent = 0
            IF ( stkptr.GT.1 ) parent = Stack(1,stkptr-1)
            istart = Stack(2,stkptr)
            stkptr = stkptr - 1
            GOTO 100
         ELSE
            next = Auux(Pntr(node)+i)
            IF ( next.NE.parent ) maxcnt = AMAX1(maxcnt,Count(next))
            i = i + 1
            GOTO 250
         ENDIF
      ELSE
         next = Auux(Pntr(node)+i)
         IF ( next.NE.parent ) THEN
            stkptr = stkptr + 1
            Stack(1,stkptr) = node
            Stack(2,stkptr) = i + 1
            parent = node
            node = next
            istart = 1
            GOTO 100
         ELSE
            i = i + 1
            GOTO 200
         ENDIF
      ENDIF
      END
**==diaend.spg  processed by SPAG 5.11R  at 11:51 on  8 Jan 1999
      INTEGER FUNCTION DIAEND(Start,Deg,Pntr,Auux,Count,Maxdph)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      INTEGER i, jdeg
C*** End of declarations inserted by SPAG
      INTEGER Start, node, parent, next, nexts, depth, Maxdph
      REAL Count(1), Deg(1), maxcnt
      INTEGER Pntr(1), Auux(1)
      node = Start
      parent = 0
      depth = 0
 100  maxcnt = 0.0
      depth = depth + 1
      jdeg = ABS(Deg(node))
      i = 1
 200  IF ( (i).LE.(jdeg) ) THEN
         next = Auux(Pntr(node)+i)
         IF ( next.NE.parent ) THEN
            IF ( Count(next).GT.maxcnt ) THEN
               maxcnt = Count(next)
               nexts = next
            ENDIF
         ENDIF
         i = i + 1
         GOTO 200
      ELSEIF ( (maxcnt.NE.0.0) .AND. (depth.LT.Maxdph) ) THEN
         parent = node
         node = nexts
         GOTO 100
      ENDIF
      DIAEND = node
      RETURN
      END
**==ksnull.spg  processed by SPAG 5.11R  at 11:51 on  8 Jan 1999
      SUBROUTINE KSNULL(Pdif1,N,M,Pdif6,Ier)
! wjm SUBROUTINE KSNULL(Pdif1,N,M,Pdif5,Pdif6,Ier)
      IMPLICIT NONE
      REAL Pdif1, Pdif6
      INTEGER M, N, Ier
C*** Start of declarations inserted by SPAG
      REAL cc, fm, fmd, fn, fnd, pdif4
      INTEGER nl100, nl50
C*** End of declarations inserted by SPAG
      REAL PROBKS            ! use routine from Numerical Recipes - WJM
      Ier = 0
      fn = N
      fm = M
      fnd = 1./fn
      fmd = 1./fm
c                                  smirnov continuity correction
      cc = SQRT(fnd)*.6666667
      nl50 = 1
      nl100 = 1
c                                  set warnings, choose approximation
c                                  probabilities
      IF ( N.LT.50 ) nl50 = 0
      IF ( N.LT.100 ) nl100 = 0
      IF ( fm*fnd.GE..1 ) THEN
         IF ( N.EQ.M*(N/M) ) THEN
            IF ( nl100.EQ.0 ) Ier = 35
         ELSE
c                                  adjust smirnov continuity correction
            cc = .6*cc
            IF ( nl100.EQ.0 ) THEN
c                                  set error indicators
               Ier = 37
               IF ( nl50.EQ.0 ) Ier = 36
            ENDIF
         ENDIF
         pdif4 = Pdif1*SQRT((fn*fm)/(fn+fm)) + cc
      ELSEIF ( nl50.EQ.0 ) THEN
         Ier = 33
         pdif4 = Pdif1*SQRT(fm)
      ELSEIF ( M.LT.80 ) THEN
         Ier = 34
         pdif4 = Pdif1*SQRT(fm)
      ELSEIF ( nl100.NE.0 ) THEN
c                                  calculate the statistic z
         pdif4 = Pdif1*SQRT((fn*fm)/(fn+fm)) + cc
      ELSE
c                                  dm
         pdif4 = Pdif1*SQRT(fm)
      ENDIF
c                                  obtain probabilities
! wjm CALL MDSMR(pdif4,Pdif5,Pdif6)
! wjm Pdif5 = 1. - Pdif5
! wjm Pdif6 = 1. - Pdif6
      Pdif6 = PROBKS(pdif4)   ! use routine from Numerical Recipes - WJM
      RETURN
      END
**==mdsmr.spg  processed by SPAG 5.11R  at 11:51 on  8 Jan 1999
      SUBROUTINE MDSMR(X,P1,P2)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL P1, P2, X, x1, x8
C*** End of declarations inserted by SPAG
      x1 = -X*X
      P1 = EXP(x1+x1)
      IF ( X.GE..22 ) THEN
         x8 = 8.*x1
         IF ( X.LE.0.8 ) THEN
            P2 = (2.506628/X)*EXP(9.869604/x8)
            P1 = 1.0 - P1
         ELSEIF ( X.LE.3.15 ) THEN
            x8 = P1*P1
            x8 = x8*x8
            P2 = P1*x8*x8 - x8 + P1
            P2 = 1.0 - P2 - P2
            P1 = 1.0 - P1
         ELSE
            P2 = 1.0
            P1 = 1.0
         ENDIF
      ELSE
         P2 = 0.0
         P1 = 1.0 - P1
      ENDIF
      RETURN
      END
**==clmssort.spg  processed by SPAG 5.11R  at 11:51 on  8 Jan 1999
! wjm SUBROUTINE SORT(V,A,Ii,Jj)
      SUBROUTINE CLMSSORT(V,A,Ii,Jj)
      IMPLICIT NONE
      REAL V
      INTEGER A, Ii, ij
      DIMENSION A(Jj), V(*)
C*** Start of declarations inserted by SPAG
      INTEGER i, il, iu, j, Jj, k, l, m
      REAL vt, vtt
C*** End of declarations inserted by SPAG
      DIMENSION iu(20), il(20)
      INTEGER t, tt
c
c     puts into a the permutation vector which sorts v into
c     increasing order.  only elements from ii to jj are considered.
c     arrays iu(k) and il(k) permit sorting up to 2**(k+1)-1 elements
c
c     this is a modification of cacm algorithm #347 by r. c. singleton,
c     which is a modified hoare quicksort.
c
      m = 1
      i = Ii
      j = Jj
 100  IF ( i.GE.j ) GOTO 500
 200  k = i
      ij = (j+i)/2
      t = A(ij)
      vt = V(ij)
      IF ( V(i).GT.vt ) THEN
         A(ij) = A(i)
         A(i) = t
         t = A(ij)
         V(ij) = V(i)
         V(i) = vt
         vt = V(ij)
      ENDIF
      l = j
      IF ( V(j).LT.vt ) THEN
         A(ij) = A(j)
         A(j) = t
         t = A(ij)
         V(ij) = V(j)
         V(j) = vt
         vt = V(ij)
         IF ( V(i).GT.vt ) THEN
            A(ij) = A(i)
            A(i) = t
            t = A(ij)
            V(ij) = V(i)
            V(i) = vt
            vt = V(ij)
         ENDIF
      ENDIF
 300  l = l - 1
      IF ( V(l).GT.vt ) GOTO 300
      tt = A(l)
      vtt = V(l)
 400  k = k + 1
      IF ( V(k).LT.vt ) GOTO 400
      IF ( k.LE.l ) THEN
         A(l) = A(k)
         A(k) = tt
         V(l) = V(k)
         V(k) = vtt
         GOTO 300
      ELSE
         IF ( l-i.LE.j-k ) THEN
            il(m) = k
            iu(m) = j
            j = l
            m = m + 1
         ELSE
            il(m) = i
            iu(m) = l
            i = k
            m = m + 1
         ENDIF
         GOTO 600
      ENDIF
 500  m = m - 1
      IF ( m.EQ.0 ) RETURN
      i = il(m)
      j = iu(m)
 600  IF ( j-i.GE.11 ) GOTO 200
      IF ( i.EQ.Ii ) GOTO 100
      i = i - 1
 700  i = i + 1
      IF ( i.EQ.j ) GOTO 500
      t = A(i+1)
      vt = V(i+1)
      IF ( V(i).LE.vt ) GOTO 700
      k = i
 800  A(k+1) = A(k)
      V(k+1) = V(k)
      k = k - 1
      IF ( vt.LT.V(k) ) GOTO 800
      A(k+1) = t
      V(k+1) = vt
      GOTO 700
      END
**==mdhyp.spg  processed by SPAG 5.11R  at 11:51 on  8 Jan 1999
      SUBROUTINE MDHYP(K,N,L,Nd,Peqk,Plek,Ier)
      IMPLICIT NONE
      INTEGER K, N, L, Nd, Ier
      REAL*8 Peqk, Plek
C*** Start of declarations inserted by SPAG
      INTEGER icnt, iflag, k0, m, mm, mm1, nnd, nxd
C*** End of declarations inserted by SPAG
      REAL*8 eps, p, sp, u, a, b, a1, b1, aa, bb, annd, anxd, aj
      DATA eps/1.0D-78/
      IF ( (K.LT.0) .OR. (K.GT.N) ) THEN
         Ier = 129
c                                  check range of n
      ELSEIF ( (N.LT.1) .OR. (N.GT.L) ) THEN
         Ier = 130
c                                  check range of nd
      ELSEIF ( (Nd.LT.0) .OR. (Nd.GT.L) ) THEN
         Ier = 131
      ELSEIF ( N.LE.32767 ) THEN
c                                  initialization for recursive formula
         Ier = 0
         p = 0.0D0
         sp = 0.0D0
         IF ( K.GT.Nd ) THEN
            sp = 1.0D0 - sp
         ELSEIF ( (N-K).LE.(L-Nd) ) THEN
            p = 1.0D0
            nnd = MIN0(N,Nd)
            nxd = MAX0(N,Nd)
            b = L
            annd = nnd
            anxd = nxd
c                                  choose more efficient direction for
c                                  computation of probabilities
            IF ( (K*(L+2)).GT.((Nd+1)*(N+1)) ) THEN
c                                  initialization for calculation of
c                                  p(x .gt. k)
               iflag = 1
               aa = b - anxd - annd
               bb = anxd - annd
               a1 = annd
               b1 = 1.0D0
               mm1 = nnd - K
               mm = L - nxd
               a = b - annd
               IF ( nnd.LT.mm ) THEN
                  mm = nnd
                  a = anxd
               ENDIF
            ELSE
c                                  initialization for calculation of
c                                  p(x .le. k)
               iflag = 0
               aa = anxd - annd
               bb = b - anxd - annd
               k0 = N + Nd - L
               IF ( k0.LT.0 ) k0 = 0
               aj = k0
               a1 = annd - aj
               b1 = aj + 1.0D0
               mm1 = K - k0
               a = b - anxd + aj
               mm = nnd - k0
            ENDIF
            icnt = 0
            IF ( mm.NE.0 ) THEN
c                                  loop for peqk
               DO m = 1, mm
                  u = a/b
c                                  detect and correct for potential
c                                  underflow
                  IF ( u.LT.eps/p ) THEN
                     p = p/eps
                     icnt = icnt + 1
                  ENDIF
                  p = p*u
                  a = a - 1.0D0
                  b = b - 1.0D0
               ENDDO
            ENDIF
            IF ( mm1.NE.0 ) THEN
c                                  loop for plek
               DO m = 1, mm1
                  IF ( icnt.EQ.0 ) sp = sp + p
                  u = a1*(a1+aa)/(b1*(b1+bb))
                  p = u*p
                  IF ( p.GE.1.0D0 ) THEN
c                                  inverse underflow adjustment
                     p = p*eps
                     icnt = icnt - 1
                  ENDIF
                  a1 = a1 - 1.0D0
                  b1 = b1 + 1.0D0
               ENDDO
            ENDIF
            IF ( icnt.NE.0 ) p = 0.0D0
            IF ( iflag.NE.0 ) THEN
               sp = 1.0D0 - sp
            ELSE
c                                  p(x .le. k)=p(x .eq. k)+p(x .lt. k)
               sp = sp + p
            ENDIF
         ENDIF
         Peqk = p
         Plek = sp
      ELSE
         Ier = 132
      ENDIF
      RETURN
      END
