anchor.diff:  C   CODED PER B SCHNABEL OCT. 1971
   C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   C
  -      DIMENSION X(1)
  +      dimension X(1)
         XMAX = 0.
  -      DO 1 I = 1,MX
  +      do 1 I = 1,MX
         XA = ABS(X(I))
  -      IF (XMAX.GT.XA) GO TO 1
  +      if (XMAX .gt. XA) go to 1
         NXMAX = I
         XMAX = XA
  -    1 CONTINUE
  -      RETURN
  -      END
  +    1 continue
  +      return
  +      end
  
positive.code:       RD(1) = zmax*GGT*T(1)*T(1)/(PIW*PIW)
        PRV(1) = zmax*GGT*T(1)/PIW
        PAA(1) = zmax
        WRITE(6,112) D
        N = 1
        YY = SQRT(1.-D*D)
        DO 200 LOOP = 2, M
        W = 6.283185307/T(N)
        WD = YY*W
        W2 = W*W
        W3 = W2*W
        CALL CMPMAX(KUG,T(N),W,W2,W3,WD,D,DT,ZD,ZV,ZA,A)
        SV(NN,N) = ZV
        SA(NN,N) = ZA/GGT
         RD(N) = ZD
        PRV(N) = W*ZD
        PAA(N) = W2*ZD/GGT
    200 N = N + 1
        WRITE(6,312) GGT, (ID(NN,I), I = 1,10),D
        SUMSV = 0.
        SUMSA = 0.
        SUMT  = 0.
  
negative.code:       SAMAX = 0.
        TT1 = .1
        TT2 = 0.
        do 320 N = 1, M
        FREKV = 1./T(N)
        if (T(N) .lt. .0999 .or. TT2.gt.2.4999)  go to 320
        TT2 = (T(N+1) + T(N))/2.
        if (TT2.gt.2.5)  TT2 = 2.5
        TT = TT2 - TT1
        SUMSA = SA(NN,N)*TT + SUMSA
        SUMSV = SV(NN,N)*TT + SUMSV
        SUMT = SUMT + TT
        TT1 = TT2
        if (SVMAX.lt.SV(NN,N))  SVMAX = SV(NN,N)
        if (SAMAX.lt.SA(NN,N))  SAMAX = SA(NN,N)
    320 write(6,322) N,T(N),RD(N),SV(NN,N),PRV(N),SA(NN,N),PAA(N),FREKV
        write(6,2002) SUMSA,SUMSV,SAMAX,SVMAX
        do 11 K = 1,KG
     11 A(K) = A(K)/GGT
        return
  C
    112 format(/5X,41HTIMES AT WHICH MAX. SPECTRAL VALUES OCCUR /
       1 10X,33HTD = TIME FOR MAX. RELATIVE DISP. /
       2 10X,33HTV = TIME FOR MAX. RELATIVE VEL.  /
       3 10X,33HTA = TIME FOR MAX. ABSOLUTE ACC.   /
       4 5X, 15HDAMPING RATIO =  F5.2)
    312 format(5X,' SPECTRAL VALUES --'/
       15X,' [Acceleration of gravity used =' F8.2 ']'/
       210A6,2X,15HDAMPING RATIO =
       3 F5.2/5X,3HNO.,4X,6HPERIOD,5X,10HREL. DISP.,6X,9HREL. VEL.,3X,
       4 12HPSU.REL.VEL.,6X,9HABS. ACC.,3X,12HPSU.ABS.ACC. 5X,5HFREQ.)
    322 format(I8,F10.2,5F15.5,F10.2)
    402 format(8F9.5)
    412 format(I5,25H ACC. RESPONSE VALUES FOR   , 8A6)
    413 format(I5,25H VEL. RESPONSE VALUES FOR   , 8A6)
   2002 format(10X,40HVALUES IN PERIOD RANGE .1 to 2.5 SEC.        /
       115X,35HAREA OF ACC. RESPONSE SPECTRUM   =   F10.3/
       215X,35HAREA OF VEL. RESPONSE SPECTRUM   =   F10.3/
       315X,35HMAX. ACCELERATION RESPONSE VALUE =   F10.3/
       415X,35HMAX. VELOCITY     RESPONSE VALUE =   F10.3)
        end
  C  ********************************************************************
        subroutine CMPMAX (KUG,PR,W,W2,W3,WD,D,DT,ZD,ZV,ZA, UG)
  C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  C
  C   THIS ROUTINE COMPUTES RESPONSE VALUES FOR ONE SINGLE DEGREE OF
  C   FREEDOM SYSTEM USING STEP BY STEP METHOD
  C
  C   EXPLANATIOS to PARAMETERS GIVEN IN DRCTSP
  C
  C   CODED BY I. M. IDRISS  1967
  C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  C
        dimension XD(2), XV(2), T(3)
        dimension UG(1)
  C
        ZA = 0.
        ZD = 0.
  