C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     SUBROUTINE OFR SOIL RESPONSE BY DIRECT INTEG. -NONLINEAR
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C                                        CODED BY Y.OHSAKI
C
C     PURPOSE
C       TO COMPUTE, BY DIRECT INTEGRATION, RESPONSE ACCELERATION TIME-
C       HISTROIES AT THE TOPS OF LAYERS IN A HORIZONTALLY LAYERED NON-
C       LINEAR SOIL DEPOSIT, WHEN INPUT ACCELERATION TIME-HISTORY IS
C       GIVEN AT THE TOP OF BASE LAYER.  TO REPRESENT NONLINEARTITY OF
C       THE SOIL DEPOSIT, OHSAKI-HARA'S MODEL IS USED
C
C     USAGE
C       CALL SDIN(L,UW,TH,STF,G0,SU,ND1,DT,NN,DDY,ACC,ND2,ACMAX,SRMAX,
C                 STMAX,SR,IND)
C
C     DESCRIPTION OF PARAMETERS
C       L            - TOTAL NUMBER OF LAYERS, INCLUDING BASE LAYER
C                      L.LE.60
C       UW(ND1)      - UNIT WEIGHT IN TONS/CU.METER
C       TH(ND1)      - THICKNESS IN METERS
C       STF(ND1)     - SOIL TYPE FACTOR
C       G0(ND1)      - INITIAL SHEAR MODULUS IN TONS/SQ.METER
C       SU(ND1)      - SHEAR STRENGTH IN TONS/SQ.METER
C       ND1          - DIMENSION OF UW,TH,STF,G0,SU,ACC,ACMAX,SRMAX,
C                      STMAX, IN CALLING PROGRAM
C       DT           - TIME INCREMENT IN SEC
C       NN           - TOTAL NUMBER OF DATA IN TIME-HISTORIES
C       DDY(ND2)     - ACCELERATION TIME-HISTORY AT BASE LAYER IN GALS
C       ACC(ND2,ND1) - RESPONSE ACCELERATION TIME-HISTOIES IN GALS
C       ND2          - DIMENSION OF DDY,ACC IN CALLING PROGRAM
C       ACMAX(ND1)   - MAX. ACCELERATIONS IN GALS
C       SRMAX(ND1)   - MAX. SHEAR STRAINS IN DECIMAL FRACTION
C       STMAX(ND1)   - MAX. SHEAR STRESSES IN TONS/SQ.METER
C       IND          - IND.EQ.0 IF BASE LAYER IS INTERNAL
C                      IND.NE.0 IF BASE LAYER IS EXPOSED
C
C     REMARKS
C       (1) THICKNESS, SOIL TYPE FACTOR, AND SHEAR STRENGTH OF BASE
C           LAYER, TH(L),STF(L),SU(L), ARE NOT REQUIRED
C       (2) IF IND=0, UNIT WEIGHT AND SHEAR MODULUS OF BASE LAYER,
C           UW(L) AND G0(L), ARE ALSO NOTREQUIRED
C       (3) OTHER THAN DAMPING EFFECT BY NON-LINEALITY, MATERIAL
C           DAMPINGS OF RAYLEIGH TYPE ARE ASSUMED IN ALL LAYERS, WHERE
C           DAMPING FACTORS FOR THE FIRST AND SECOND MODES ARE 0.02
C       (4) FOR NUMERIACAL INTEGRATION OF EQUATIONS OF MOTION, WILSON'S
C           METHOD IS USED, IN WHICH THETA=1.4
C
C     SUBROUTINES AND FUCNTION SUBPROGRAMS REQUIRED
C       EIGN OHYS SKEL
C
      SUBROUTINE SDIN(L,UW,TH,STF,G0,SU,ND1,DT,NN,DDY,ACC,ND2,ACMAX,
     1                SRMAX,STMAX,SR,IND)
C
      DIMENSION UW(ND1), TH(ND1), STF(ND1), G0(ND1), SU(ND1), DDY(ND2),
     1          ACC(ND2,ND1), ACMAX(ND1), SRMAX(ND1), STMAX(ND1),
     2       SR(ND2,ND1)
      DIMENSION B(60), DDX(60), DX(60), X(60), U(60), DU(60), U1(60),
     1       DU1(60), V(60), G1(60), G01(60), SU1(60), W(2), EM(60,2),
     2       EK(60,2), EC(60,2), EMT(60,2), RM(60), RC(60), R(60),
     3       DELX(60)
      EQUIVALENCE (U(1),RM(1)),(DU(1),RC(1))
      DATA      THETA / 1.4 /, H0 / 0.02 /
C                    INITIALIZATION
      TAU = THETA*DT
      A1 = 6.0/TAU**2
      A2 = 3.0/TAU
      A3 = A2*2.0
      A4 = TAU/2.0
      A5 = A1/THETA
      A6 =-A3/THETA
      A7 = 1.0-3.0/THETA
      A8 = DT/2.0
      A9 = A8*DT/3.0
      UWL = UW(L)
      G0L = G0(L)
      UW(L)  = 0.0
      TH(L)  = 0.00001
      G0(L)  = 0.0
      STF(L) = 1.0
      DO 110  I = 1, L
      B(I) = 1.4+0.2*STF(I)
      DDX(I) = -DDY(1)
      DX(I)  = 0.0
      X(I)   = 0.0
      ACC(1,I) = 0.0
      U1(I)  = 0.0
      DU1(I) = 0.0
      G01(I) = G0(I)
      G1(I)  = G0(I)
      SU1(I) = SU(I)
      ACMAX(I) = 0.0
      SRMAX(I) = 0.0
      STMAX(I) = 0.0
 110  CONTINUE
C                    CONSISTENT MASS MATRIX AND INITIAL STIFFNESS MATRIX
      EM1 = UW(1)*TH(1)/58.8
      EK1 = G0(1)/TH(1)
      EM(1,1) = EM1*2.0
      EK(1,1) = EK1
      DO 120  I = 2, L
      EM2 = UW(I)*TH(I)/58.8
      EK2 = G0(I)/TH(I)
      EM(I-1,2) = EM1
      EK(I-1,2) =-EK1
      EM(I,1) = (EM1+EM2)*2.0
      EK(I,1) = EK1+EK2
      EM1 = EM2
      EK1 = EK2
 120  CONTINUE
      UW(L) = UWL
      G0(L) = G0L
C                    INITIAL DAMPING MATRIX
      N = L-1
      IF (IND.NE.0) N = L
      CALL EIGN(N,EK,EM,60,W,IND)
      W(1) = SQRT(W(1))
      W(2) = SQRT(W(2))
      WRITE(6,125) W
  125 FORMAT('0EIGENVALUE  W(1)=',1PE13.4,'  W(2)=',E13.4)
      B0 = 2.0*H0/(W(1)+W(2))
      A0 = B0*W(1)*W(2)
      DO 130  I = 1, L
      EC(I,1) = A0*EM(I,1)+B0*EK(I,1)
      EC(I,2) = A0*EM(I,2)+B0*EK(I,2)
 130  CONTINUE
      EC(L,1) = EC(L,1)+SQRT(UW(L)*G0(L)/9.8)
C                    TRIANGULARIZATION OF MASS MATRIX
      EMT(1,1) = SQRT(EM(1,1))
      DO 140  I = 2, L
      EMT(I-1,2) = EM(I-1,2)/EMT(I-1,1)
      EMT(I,1) = SQRT(EM(I,1)-EMT(I-1,2)**2)
 140  CONTINUE
C                    STEP-BY-STEP INTEGRATION
      DO 420  M = 2, NN
C                    EFFECTIVE STIFFNESS MATRIX
      EK1 = G1(1)/TH(1)
      EK(1,1) = EK1+A1*EM(1,1)+A2*EC(1,1)
      IF (N.EQ.1) GO TO 160
      DO 150  I = 2, N
      EK2 = G1(I)/TH(I)
      EK(I-1,2) = -EK1    +A1*EM(I-1,2)+A2*EC(I-1,2)
      EK(I,1)   =  EK1+EK2+A1*EM(I,1)  +A2*EC(I,1)
      EK1 = EK2
 150  CONTINUE
C                    TRIANGULARIZATION OF EFFECTIVE STIFFNESS MATRIX
 160  EK(1,1) = SQRT(EK(1,1))
      IF (N.EQ.1) GO TO 180
      DO 170  I = 2, N
      EK(I-1,2) = EK(I-1,2)/EK(I-1,1)
      EK(I,1) = SQRT(EK(I,1)-EK(I-1,2)**2)
 170  CONTINUE
C                    EFFECTIVE LOAD VECTOR
 180  R1 = -THETA*(DDY(M)-DDY(M-1))
      DO 190  I = 1, N
      RM(I) = R1+A3*DX(I)+3.0*DDX(I)
      RC(I) = 3.0*DX(I)+A4*DDX(I)
 190  CONTINUE
      IF (N.GE.2) GO TO 200
      R(1) = EM(1,1)*RM(1)+EC(1,1)*RC(1)
      GO TO 230
 200  R(1) = EM(1,1)*RM(1)+EM(1,2)*RM(2)+EC(1,1)*RC(1)+EC(1,2)*RC(2)
      IF (N.EQ.2) GO TO 220
      DO 210  I = 2, N-1
      R(I) = EM(I-1,2)*RM(I-1)+EM(I,1)*RM(I)+EM(I,2)*RM(I+1)
     1      +EC(I-1,2)*RC(I-1)+EC(I,1)*RC(I)+EC(I,2)*RC(I+1)
 210  CONTINUE
 220  R(N) = EM(N-1,2)*RM(N-1)+EM(N,1)*RM(N)
     1      +EC(N-1,2)*RC(N-1)+EC(N,1)*RC(N)
C                    SOLUTION FOR INCREMENTAL DISPLACEMENTS
 230  R(1) = R(1)/EK(1,1)
      IF (N.EQ.1) GO TO 250
      DO 240  I = 2, N
      R(I) = (R(I)-EK(I-1,2)*R(I-1))/EK(I,1)
 240  CONTINUE
 250  DELX(N) = R(N)/EK(N,1)
      IF (N.EQ.1) GO TO 270
      DO 260  I = 1, N-1
      NMI = N-I
      DELX(NMI) = (R(NMI)-EK(NMI,2)*DELX(NMI+1))/EK(NMI,1)
 260  CONTINUE
C                    RELATIVE ACCELERATIONS, SHEAR STRAINS AND STRESSES
 270  DO 280  I = 1, N
      DDXI = DDX(I)
      DXI = DX(I)
      DDX(I) = A5*DELX(I)+A6*DXI+A7*DDXI
      DX(I)= DXI+A8*(DDX(I)+DDXI)
      X(I) = X(I)+DT*DXI+A9*(DDX(I)+2.0*DDXI)
 280  CONTINUE
      IF (N.EQ.1) GO TO 300
      DO 290  I = 1, N-1
      U(I) = (X(I)-X(I+1))/TH(I)*0.01
      DU(I) = (DX(I)-DX(I+1))/TH(I)
 290  CONTINUE
 300  U(N) = X(N)/TH(N)*0.01
      DU(N) = DX(N)/TH(N)
      ICALL = M-1
C     WRITE(6,*) ' ** FROM SDIN ACC *** '
C     WRITE(6,888) M, (DDX(I),I=1,N)
C 888 FORMAT('0STEP=',I4,' ACCELERATION  DDX'/(1H ,10(1PE12.3)))
C      WRITE(6,887) (DU(I),I=1,N)
C  887 FORMAT(' STRAIN VELOCITY  DU'/(1H ,10(1PE12.3)))
C      WRITE(6,886) (V(I),I=1,N)
C  886 FORMAT(' STRESS  V'/(1H ,10(1PE12.3)))
C
      CALL OHYS(L-1,U,DU,U1,DU1,G01,SU1,B,V,G1,60,ICALL)
C
C      WRITE(6,885) (G1(I),I=1,N)
C  885 FORMAT(' SHEAR MODULUS  G1'/(1H ,10(1PE12.3)))
      DO 310  I = 1, L-1
      SRMAX(I) = AMAX1(SRMAX(I),ABS(U(I)))
      STMAX(I) = AMAX1(STMAX(I),ABS(V(I)))
      U1(I) = U(I)
      DU1(I) = DU(I)
      SR(M,I)=U(I)
 310  CONTINUE
C       IF (M.GT.50) GOTO 1111
C      WRITE(6,*)
C      WRITE(6,*) ' STEP= ',M
C      WRITE(6,*) ' DX(2,3,4) '
C      WRITE(6,*) (DX(I),I=2,4)
C      WRITE(6,*) ' V(2,3,4) '
C      WRITE(6,*) (V(I),I=2,4)
C      WRITE(6,*) ' EC(I,1) I=2,3,4 '
C      WRITE(6,*) (EC(I,1),I=2,4)
C      WRITE(6,*) ' EC(I,2) I=2,3,4 '
C      WRITE(6,*) (EC(I,2),I=2,4)
C      WRITE(6,*) ' R(I) I=2,3 '
C      WRITE(6,*) (R(I),I=2,3)
C      WRITE(6,*) ' EMT(I,1) I=2,3,4 '
C      WRITE(6,*) (EMT(I,1),I=2,4)
C      WRITE(6,*) ' EMT(I,2) I=2,3,4 '
C      WRITE(6,*) (EMT(I,2),I=2,4)
C      WRITE(6,*) ' DDX(2,3,4) '
C      WRITE(6,*) (DDX(I),I=2,4)
C                    RE-EVALUATION OF RELATIVE ACCELERATIONS
C111   IF (N.GE.2) GO TO 320
C      R(1) = EC(1,1)*DX(1)+V(1)*100.0
C      GO TO 350
C320   R(1) = EC(1,1)*DX(1)+EC(1,2)*DX(2)+V(1)*100.0
C      IF (N.EQ.2) GO TO 340
C      WRITE(6,*) ' DX(2,3,4) '
C      WRITE(6,*) (DX(I),I=2,4)
C      WRITE(6,*) ' V(2,3,4) '
C      WRITE(6,*) (V(I),I=2,4)
C      WRITE(6,*) ' EC(I,1) I=2,3,4 '
C      WRITE(6,*) (EC(I,1),I=2,4)
C      WRITE(6,*) ' EC(I,2) I=2,3,4 '
C      WRITE(6,*) (EC(I,2),I=2,4)
C      DO 330  I = 2, N-1
C      R1 = EC(I-1,2)*DX(I-1)
C      R2 = EC(I,1)*DX(I)
C      R3 = EC(I,2)*DX(I+1)
C      R4 = (V(I)-V(I-1))*100.0
C      R(I)=R1+R2+R3+R4
C        IF (I.EQ.3) WRITE(6,*) ' R1,2,3,4 ',R1,R2,R3,R4
C330   CONTINUE
C        IF (M.LE.50) WRITE(6,*) 'RE-R3 ',R(3)
C340   IF (IND.NE.0) V(N) = 0.0
C      R(N) = EC(N-1,2)*DX(N-1)+EC(N,1)*DX(N)+(V(N)-V(N-1))*100.0
C350   R(1) = R(1)/EMT(1,1)
C      IF (N.EQ.1) GO TO 370
C      DO 360  I = 2, N
C      R(I) = (R(I)-EMT(I-1,2)*R(I-1))/EMT(I,1)
C360   CONTINUE
C370   DDX(N) = R(N)/EMT(N,1)
C        IF (M.LE.50) WRITE(6,*) 'RE-R2,R3 ',R(2),R(3)
C      IF (N.EQ.1) GO TO 390
C      DO 380  I = 1, N-1
C      NMI = N-I
C      DDX(NMI) = (R(NMI)-EMT(NMI,2)*DDX(NMI+1))/EMT(NMI,1)
C380   CONTINUE
C390   DO 400  I = 1, N
C      DDX(I) = -DDY(M)-DDX(I)
C400   CONTINUE
C      WRITE(6,889) (DDX(I),I=1,N)
C 889  FORMAT(' RELATIVE ACC  DDX'/(1H ,10(1PE12.3)))
C                    ABSOLUTE ACCELERATIONS
      DO 410  I = 1, N
      ACC(M,I) = DDX(I)+DDY(M)
      ACMAX(I) = AMAX1(ACMAX(I),ABS(ACC(M,I)))
 410  CONTINUE
C      WRITE(6,*) 'NEXT STEP?'
C      READ(5,14) ANS
C   14 FORMAT(A1)
C      IF (ANS.EQ.'N') RETURN
 420  CONTINUE
      RETURN
      END
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     SUBROUTINE FOR EIGENVALUES OF TWO LOWEST MODES - CONTINUANT
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C                                        CODED BY Y.OHSAKI
C     PURPOSE
C       TO COMPUTE THE LOWEST AND THE SECOND LOWEST EIGENVALUES OF
C       A*U=EIG*B*U, WHERE A,B ARE SYMMETRIC CONTINUANT MATRICES
C
C     USAGE
C       CALL EIGN(N,A,B,ND,EIG,IND)
C
C     DESCRIPTION OF PARAMETERS
C       N       - ORDER OF MATRICES A,B  N.LE.60
C       A(ND,2) - DIAGONAL AND SUB-DIAGONAL ELEMENTS OF MATRIX A
C       B(ND,2) - DIAGONAL AND SUB-DIAGONAL ELEMENTS OF MATRIX B
C       ND      - DIMENSION OF A,B IN CALLING PROGRAM
C       EIG(2)  - THE LOWEST AND THE SECOND LOWEST EIGENVALUES
C       IND     - IND=0 IF MATRIX A IS REGULAR
C                 IND=1 IF MATRIX A IS SINGULAR
C
C     REMARKS
C       (1) A(N,2) AND B(N,2) ARE NOT REQUIRED
C       (2) IF MATRIX A IS SINGULAR, THE LOWEST AND THE SECOND LOWEST
C           EIGENVALUES ARE OBTAINED, EXCLUDING ZERO-EIGENVALUE
C
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C       NONE
C
      SUBROUTINE EIGN(N,A,B,ND,EIG,IND)
C
      DIMENSION A(ND,2),  B(ND,2),  EIG(2)
      DIMENSION AA(60,2), BB(60,60), U(60), AU(60), BU(60)
C                    INITIALIZATION
      IF (N.EQ.1) GO TO 360
      DO 120  I = 1, N
      AA(I,1) = A(I,1)
      AA(I,2) = A(I,2)
      DO 110  J = I, N
      BB(I,J) = 0.0
 110  CONTINUE
      BB(I,I) = B(I,1)
 120  CONTINUE
      DO 130  I = 1, N
      BB(I,I+1) = B(I,2)
      BB(I+1,I) = B(I,2)
 130  CONTINUE
      AA(1,1) = SQRT(AA(1,1))
      DO 140  I = 2, N
      AA(I-1,2) = AA(I-1,2)/AA(I-1,1)
      AA(I,1) = SQRT(AA(I,1)-AA(I-1,2)**2)
 140  CONTINUE
      NN = N
C                    REDUCTION FOR SINGULAR SYSTEM
      IF (IND.EQ.0) GO TO 190
      BU(1) = B(1,1)+B(1,2)
      SSB = BU(1)
      IF (N.EQ.2) GO TO 160
      DO 150  I = 2, N-1
      BU(I) = B(I-1,2)+B(I,1)+B(I,2)
      SSB = SSB+BU(I)
 150  CONTINUE
 160  BU(N) = B(N-1,2)+B(N,1)
      SSB = SSB+BU(N)
      DO 180  I = 1, N
      DO 170  J = I, N
      BB(I,J) = BB(I,J)-BU(I)*BU(J)/SSB
      BB(J,I) = BB(I,J)
 170  CONTINUE
 180  CONTINUE
      NN = N-1
C                    ITERATION
 190  DO 340  IEIG = 1, 2
      U(1) = 1.0
      DO 200  I = 2, NN
      U(I) = 0.0
 200  CONTINUE
      EIGEN = 0.0
      DO 260  ITER = 1, 10
      DO 220  I = 1, NN
      SB = 0.0
      DO 210  K = 1, NN
      SB = SB+BB(I,K)*U(K)
 210  CONTINUE
      BU(I) = SB
 220  CONTINUE
      BU(1) = BU(1)/AA(1,1)
      DO 230  I = 2, NN
      BU(I) = (BU(I)-AA(I-1,2)*BU(I-1))/AA(I,1)
 230  CONTINUE
      U(NN) = BU(NN)/AA(NN,1)
      S = U(NN)**2
      DO 240  I = 1, NN-1
      NNMI = NN-I
      U(NNMI) = (BU(NNMI)-AA(NNMI,2)*U(NNMI+1))/AA(NNMI,1)
      S = S+U(NNMI)**2
 240  CONTINUE
      S = SQRT(S)
      DO 250  I = 1, NN
      U(I) = U(I)/S
 250  CONTINUE
      IF (ABS(S-EIGEN)/S.LT.0.01) GO TO 270
      EIGEN = S
 260  CONTINUE
C                    RAYLEIGH'S QUOTIENT
 270  AU(1) = A(1,1)*U(1)+A(1,2)*U(2)
      SSA = AU(1)*U(1)
      IF (NN.EQ.2) GO TO 290
      DO 280  I = 2, NN-1
      AU(I) = A(I-1,2)*U(I-1)+A(I,1)*U(I)+A(I,2)*U(I+1)
      SSA = SSA+AU(I)*U(I)
 280  CONTINUE
 290  AU(NN) = A(NN-1,2)*U(NN-1)+A(NN,1)*U(NN)
      SSA = SSA+AU(NN)*U(NN)
      SSB = 0.0
      DO 310  I = 1, NN
      SB  = 0.0
      DO 300  K = 1, NN
      SB  = SB+BB(I,K)*U(K)
 300  CONTINUE
      SSB = SSB+U(I)*SB
 310  CONTINUE
      EIG(IEIG) = SSB/SSA
      IF (IEIG.EQ.2) GO TO 350
C                    DEFLATION FOR SECOND EIGENVALUE
      DO 330  I = 1, NN
      DO 320  J = I,NN
      BB(I,J) = BB(I,J)-EIG(1)*AU(I)*AU(J)/SSA
      BB(J,I) = BB(I,J)
 320  CONTINUE
 330  CONTINUE
 340  CONTINUE
 350  EIG(1) = 1.0/EIG(1)
      EIG(2) = 1.0/EIG(2)
      RETURN
 360  EIG(1) = A(1,1)/B(1,1)
      RETURN
      END
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     SUBROUTINE FOR OHSAKI-HARA'S HYSTERESIS MODEL
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C                                     CODED BY Y.OHSAKI
C
C     PURPOSE
C       TO COMPUTE SHEAR STRESSES AND TANGENTIAL SHEAR MODULI AT M-TH
C       STEP IN SOIL LAYERS CHARACTERIZED BY OHSAKI-HARA'S HYSTERESIS
C       MODEL, WHEN STRAINS AND STRAIN VELOCITIES AT M-TH AND (M-1)TH
C       STEP ARE GIVEN
C
C     USAGE
C       CALL OHYS(L,U,DU,U1,DU1,G0,SU,B,V,G,ND,ICALL)
C
C     DESCRIPTION OF PARAMETERS
C       L       - TOTAL NUMBER OF LAYERS  L.LE.60
C       U(ND)   - STRAINS AT M-TH STEP IN DECIMAL FRACTION
C       DU(ND)  - STRAIN VELOCITIES AT M-TH STEP
C       U1(ND)  - STRAINS AT (M-1)TH STEP IN DECIMAL FRACTION
C       DU1(ND) - STRAIN VELOCITIES AT (M-1)TH STEP
C       G0(ND)  - INITIAL SHEAR MODULI IN TONS/SQ.METER
C       SU(ND)  - SHEAR STRENGTHS IN TONS/SQ.METER
C       B(ND)   - CONSTANTS IN THE EQUATION OF SKELETON CURVE
C       V(ND)   - STRESSES AT M-TH STEP IN TONS/SQ.METER
C       G(ND)   - TANGENTIAL SHEAR MODULI IN TONS/SQ.METER
C       ND      - DIMENSION OF U,DU,U1,DU1,G0,SU,B,V,G IN CALLING
C                 PROGRAM
C       ICALL   - FOR THE FIRST CALL, ICALL.EQ.1  OTHERWIASE, ICALL.NE.1
C
C     SUBROUTINES AND FUCTION SUBPROGRAMS REQUIRED
C       SKEL
C
      SUBROUTINE OHYS(L,U,DU,U1,DU1,G0,SU,B,V,G,ND,ICALL)
C
      DIMENSION U(ND), DU(ND), U1(ND), DU1(ND), G0(ND), SU(ND), B(ND),
     1          V(ND), G(ND)
      DIMENSION A(60), C(60), K(60), D(60), U0(500,60), V0(500,60),
     1          VAPRX(60)
C                    INITIALIZATION
      IF (ICALL.NE.1) GO TO 120
      DO 110  I = 1, L
      C(I) = G0(I)/SU(I)
      A(I) = C(I)*0.01-1.0
C      IF (I.LE.14) R05=0.003
C      IF (I.LE.19.AND.I.GT.14) R05=0.0005
C      IF (I.GT.19) R05=0.0008
C     A(I) = (2/(R05*G0(I)))**1.3
      K(I) = 1
      D(I) = 1.0
      U0(1,I) = 0.0
      V0(1,I) = 0.0
      VAPRX(I) = 0.0
 110  CONTINUE
 120  DO 170  I = 1, L
C                    TURNING POINT
      IF (DU(I)*DU1(I).GE.0.0) GO TO 130
      UTN = (DU(I)*U1(I)-DU1(I)*U(I))/(DU(I)-DU1(I))
      UNORM = (UTN-U0(K(I),I))/D(I)
      CALL SKEL(UNORM,VNORM,GNORM,VAPRX(I),A(I),B(I),C(I))
      VTN = D(I)*VNORM+V0(K(I),I)
      K(I) = K(I)+1
      D(I) = 2.0
      U0(K(I),I) = UTN
      V0(K(I),I) = VTN
      GO TO 140
C                    TERMINAL POINT
 130  IF (K(I).EQ.1) GO TO 160
 140  IF (K(I).EQ.2) UEND = -U0(2,I)
      IF (K(I).GT.2) UEND =  U0(K(I)-1,I)
      IF ((UEND-U(I))*SIGN(1.0,DU(I)).GE.0.0) GO TO 160
      IF (K(I).GT.3) GO TO 150
      K(I) = 1
      D(I) = 1.0
      GO TO 160
 150  K(I) = K(I)-2
C                    SHEAR STRESS AND TANGENTIAL MODULUS
 160  UNORM = (U(I)-U0(K(I),I))/D(I)
      CALL SKEL(UNORM,VNORM,GNORM,VAPRX(I),A(I),B(I),C(I))
      V(I) = (D(I)*VNORM+V0(K(I),I))*SU(I)
      G(I) = GNORM*G0(I)
      VAPRX(I) = VNORM
 170  CONTINUE
      RETURN
      END
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     SUBROUTINE FOR STRESS ON SKELETON CURVE
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C                                        CODED BY Y.OHSAKI
C
C     PURPOSE
C       TO COMPUTE SHEAR STRESS ON OHSAKI-HARA'S SKELETON CURVE FOR
C       GIVEN STRAIN BY MEANS OF NEWTON-RAPHSON'S METHOD. TANGENTIAL
C       SHEAR MODULUS AT THE GIVEN STRAIN IS ALSO COMPUTED
C
C     USAGE
C       CALL SKEL(U,V,G,V1,A,B,C)
C
C     DESCRIPTION OF PARAMETERS
C
C       U     - GIVEN STRAIN
C       V     - STRESS NORMALIZED IN TERMS OF SHEAR STRENGTH
C       G     - TANGENTIAL MODULUS NORMALIZED IN TERMS OF INITIAL
C               MODULUS
C       V1    - APPROXIMATE VALUE OF V
C       A,B,C - PARAMETERS IN OHSAKI-HARA'S EQUATION OF SKELETON CURVE
C
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C       INCORPORATED
C
      SUBROUTINE SKEL(U,V,G,V1,A,B,C)
C
      DOUBLE PRECISION DV
C
      CU = C*U
      V = V1
      CALL FUNC(V,FV,DFV,A,B,CU)
      DO 140  I = 1, 20
      IF (ABS(FV).LE.0.0001) GO TO 150
      DV = DBLE(FV)/DBLE(DFV)
      V = SNGL(DBLE(V)-DV)
      IF (ABS(V).GT.0.01) GO TO 120
      H = 0.0
      CALL FUNC(H,FH,DFH,A,B,CU)
      IF (ABS(FH).GT.0.0001) GO TO 110
      V = 0.0
      G = 1.0/DFH
      RETURN
 110  IF (SNGL(DABS(DV)).LE.0.01) GO TO 150
      GO TO 130
 120  DV = DV/DBLE(V)
      IF (SNGL(DABS(DV)).GT.0.01) GO TO 130
      CALL FUNC(V,FV,DFV,A,B,CU)
      IF (ABS(FV).LE.0.0001) GO TO 150
 130  CALL FUNC(V,FV,DFV,A,B,CU)
 140  CONTINUE
 150  G = 1.0/DFV
      RETURN
      END
C
      SUBROUTINE FUNC(V,FV,DFV,A,B,CU)
C
      VB = ABS(V)**B
      FV = V*(1.0+A*VB)-CU
      DFV = 1.0+A*(1.0+B)*VB
      RETURN
      END
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     SUBROUTINE FOR SUBDIVISION OF SOIL LAYERS
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C                                        CODED BY Y.OHSAKI
C     PURPOSE
C       TO SUBDEVIDE EACH LAYER OF A SOIL DEPOSIT INTO SUBLAYERS WITH
C       EQUAL THICKNESS, THE NUMBER OF SUBDIVISION BEING DEPENDENT ON
C       FUNDAMENTAL PERIOD OF EACH LAYER
C
C     USAGE
C       CALL SUBD(L,UW,TH,G,DAMP,STF,ND)
C
C     DESCRIPTION OF PARAMETERS
C       L        - TOTAL NUMBER OF LAYERS, INCLUDING BASE LAYER
C                  L(AT RETURN).LE.60
C       UW(ND)   - UNIT WEIGHT
C       TH(ND)   - THICKNESS
C       G(ND)    - SHEAR MODULUS
C       DAMP(ND) - DAMPING FACTOR
C         STF(ND)  - SOIL TYPE FACTOR
C       ND       - DIMENSION OF UW,TH,G,DAMP IN CALLING PROGRAM
C
C     REMARKS
C       (1) L,UW,TH,G,DAMP ARE OF ORIGINAL/SUBDIVIDED LAYERS AT
C           CALL/RETURN
C       (2) BASE LAYER IS NO SUBDIVIDED
C
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C       NONE
C
      SUBROUTINE SUBD(L,UW,TH,G,DAMP,STF,ND)
C
      DIMENSION UW (ND), TH (ND), G (ND), DAMP (ND), STF(ND)
      DIMENSION UW1(60), TH1(60), G1(60), DAMP1(60), TDIV(6), STF1(60)
      DATA      TDIV / 0.033, 0.04, 0.045, 0.055, 0.09, 0.33 /
C
      L1 = 0
      DO 160  I = 1, L
      IF (I.EQ.L) GO TO 120
      T = 4.0*TH(I)*SQRT(UW(I)/G(I)/9.8)
      DO 110  K = 1, 6
      IF (T.LE.TDIV(K)) GO TO 130
 110  CONTINUE
      NDIV = 14
      GO TO 140
 120  NDIV = 1
      GO TO 140
 130  NDIV = K+1
 140  KS = L1+1
      KE = L1+NDIV
      DO 150  K = KS, KE
      UW1(K) = UW(I)
      TH1(K) = TH(I)/FLOAT(NDIV)
      G1(K)  = G(I)
      DAMP1(K) = DAMP(I)
      STF1(K) = STF(I)
      L1 = L1+1
 150  CONTINUE
 160  CONTINUE
      L = L1
      IF (L.GT.ND) GO TO 180
      DO 170  I = 1, L
      UW(I) = UW1(I)
      TH(I) = TH1(I)
      G(I)  = G1(I)
      DAMP(I) = DAMP1(I)
      STF(I) = STF1(I)
 170  CONTINUE
      RETURN
 180  WRITE(6,600) L
      STOP
C
 600  FORMAT(1H0,'ND MUST BE LARGER THAN',I3)
C
      END
