C
      SUBROUTINE RUNGE(M,V,X,F,SIZE,SPLENG,VZERO,XZERO,RK,NROW,NCOL,
     1 NSPRNG,NUMASS,FINALK,PRIMEK)
      INTEGER*4  STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR
      COMMON T,TLIMIT,H,SUM(3)
      COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC,
     2 KOUNT
      DIMENSION M(7,10),V(10,3),X(10,3),F(10,3),SIZE(10)
      DIMENSION SPLENG(10),FINALK(10,3),PRIMEK(10,3)
      DIMENSION VZERO(10,3),XZERO(10,3),RK(10,3,4)
C<<< INTEGRATING LOOP STARTS HERE
      J=0
      KOUNT=KOUNT+1
      TZERO=T
      N=0
C<<< STORE VALUES AT TIME TZERO
      DO 900 I=STMASS,ENDMAS
      N=N+1
      DO 900 K=1,3
      XZERO(N,K)=X(I,K)
      VZERO(N,K)=H*V(I,K)
  900 CONTINUE
C<<< CALCULATE RK CONSTANTS
  100 J=J+1
      IF(J.EQ.1) GO TO 901
      IF((NVDR+NFOR).EQ.0) GO TO 201
      CALL DRIVER(V,X,F,NCOL)
 201  CONTINUE
      CALL VELPOS(V,X,M,NROW,NCOL)
      IF(NDAM.EQ.0) GO TO 203
      CALL DAMPER(V,X,F,SIZE,NROW,NCOL)
 203  CONTINUE
      IF(NSPR.EQ.0) GO TO 204
      CALL SPRING(X,F,SIZE,SPLENG,NCOL,NSPRNG)
 204  CONTINUE
  901 CONTINUE
      N=0
      ITRANS=0
      N1=NN
      N2=NB
      DO 903 I=STMASS,ENDMAS
      N=N+1
      IROW=I
      CALL SYSUM(F,M,IROW,N1,N2,ITRANS,NROW,NCOL)
      DO 902 K=1,3
      F(I,K)=SUM(K)
      ACCELN= (F(I,K)/SIZE(I))
      RK(N,K,J)=0.5*H*H*ACCELN
  902 CONTINUE
  903 CONTINUE
      IF(J-3) 101,102,103
  101 CONTINUE
C<<< CALCULATE RK(2) AND RK(3)
      T=TZERO+0.5*H
      N=0
      DO 904 I=STMASS,ENDMAS
      N=N+1
      DO 904 K=1,3
      X(I,K)=XZERO(N,K)+0.5*VZERO(N,K)+0.25*RK(N,K,1)
      V(I,K)=(VZERO(N,K)+RK(N,K,J))/H
  904 CONTINUE
      GO TO 100
  102 CONTINUE
C<<< CALCULATE RK(4)
      T=TZERO+H
      N=0
      DO 905  I=STMASS,ENDMAS
      N=N+1
      DO 905 K=1,3
      X(I,K)=XZERO(N,K)+VZERO(N,K)+RK(N,K,3)
      V(I,K)=(VZERO(N,K)+2.*RK(N,K,3))/H
  905 CONTINUE
      GO TO 100
  103 CONTINUE
C<<< CALCULATE FINAL VALUES OF V AND X
      N=0
      DO 906  I=STMASS,ENDMAS
      N=N+1
      DO 906 K=1,3
      FINALK(N,K)=(RK(N,K,1)+RK(N,K,2)+RK(N,K,3))/3.0
      PRIMEK(N,K)=(RK(N,K,1)+2.*RK(N,K,2)+2.*RK(N,K,3)+RK(N,K,4))/3.0
      X(I,K)=XZERO(N,K)+VZERO(N,K)+FINALK(N,K)
      V(I,K)=(VZERO(N,K)+PRIMEK(N,K))/H
  906 CONTINUE
      RETURN
      END
C
      SUBROUTINE SYSUM(V,M,IROW,N1,N2,ITRANS,NROW,NCOL)
      INTEGER*4  STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR
      COMMON T,TLIMIT,H,SUM(3)
      COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC,
     2 KOUNT
      DIMENSION V(10,3),M(7,10)
C
C     SUBROUTINE TO CALCULATE M*V OR -M(TRANSPOSE)*V
C     SUM =-M(TRANSPOSE)*V IF ITRANS=1
C     SUM = M*V  IF ITRANS=0
C
C
      SUM(1)=0.
      SUM(2)=0.
      SUM(3) =0.
      IF(N2-N1) 40,41,41
  41  CONTINUE
      IF (ITRANS) 1000,42,43
C*****CALCULATE -M(TRANSPOSE)*V
  43  CONTINUE
      DO 44 K=1,3
      DO 44 J=N1,N2
      IF(M(J,IROW)) 45,44,46
  45  SUM(K)=SUM(K)+V(J,K)
      GO TO 44
  46  SUM(K)=SUM(K)-V(J,K)
  44  CONTINUE
      GO TO 40
C*****CALCULATE M*V
  42  CONTINUE
      DO 47 K=1,3
      DO 47 J=N1,N2
      IF(M(IROW,J)) 48,47,49
  48  SUM(K)=SUM(K)-V(J,K)
      GO TO 47
  49  SUM(K)=SUM(K)+V(J,K)
  47  CONTINUE
 40   RETURN
1000  CONTINUE
C
      STOP
      END