C---    DPA BY MARK S. JOHNSON AND LES CADWALLADER
C---        MAR. 1985
C---        JAN  2015 REV FOR WINDOWS 7
C---      March 2023 REV FOR WINDOWS 10 BY GTW      
C---   THIS PROGRAM GENERATES P-Y CURVES AND DATA FILE  TAPE5 (?)
      PROGRAM DPA
      USE Fortran_WinPrint
      DIMENSION SAVE(200), PSF(200), IT(15)
      CHARACTER IANS*5, IDT*8, TITLE*60,FLNAME*20
      LOGICAL EXX
      INTEGER TYPE(20), PYSTA(31)
      REAL MT, MOM
      COMMON /BK1/ TITLE
      COMMON /RANGE/ IRANGE,INCR,THICKNESS,E,EI,IPIPE
      COMMON /LOADS/ ILOADS,KD
      COMMON /BK3/ MT, HT, VT,IPROB, DEPTH, DIAM, CAP, Z(20), C(20),
     1DENSE(20), E50(20), PHI(20), H2O, NUMST, IU8, TYPE,IDT
      COMMON /BK4/ IEND
100   FORMAT (' PROGRAM TERMINATED - DEPTH OR DIAMETER = 0')
102   FORMAT ('00001      DPA - BEAMCOL 28 VERSION - IBM 2.4.1',/
     1,'00002  ',A60, /
     1,'00003   LATERALLY LOADED DRILLED SHAFT FOUNDATIONS')
103   FORMAT(I2.2, '000  1  H=', F6.1, ' D=', F4.1, ' CAP=', F4.1,1P,
     1' MT=', E10.3, ' HT=', E10.3, ' VT=', E10.3)
104   FORMAT(I2.2, '200', I5, T21, F4.0, T31, '100.0', T41,F5.3)
105   FORMAT (I2.2, '401',  T31, 1P, E10.3, T51, E10.3)
106   FORMAT (I2.2, '402', T20, '1', 1PE10.3, T71, 1PE10.3)
107   FORMAT (I2.2, '403', T11, I5, T21, 1PE10.3)
108   FORMAT ('00001      DPA - BEAMCOL 28 VERSION - IBM 2.4.1',/
     1,'00002  ',A60, /
     1,'00003   STEEL PIPE PILE FOUNDATION, WALL THICKNESS = ',
     1   F8.5,' IN.')
120   FORMAT (//////,15X,35('*'),//,19X,'DPA BY MARK S. JOHNSON',
     1 /,17X,'IBM VERSION 2.4.1 - JAN 1987'
     1//,15X,(35('*')),///,
     1T20,'SOIL TYPES    DEFAULTS = *', //
     1T20,'TYPE 0 - DEFAULT', /
     1T20,'TYPE 1 - SOFT CLAY',   /
     1T20,'TYPE 2 - STIFF CLAY',   /
     1T18,'* TYPE 3 - SAND BELOW THE WATER TABLE',   /
     1T20,'TYPE 4 - STIFF CLAY BELOW THE WATER TABLE',      /
     1T20,'TYPE 5 - PADLL',/,
     1T18,'* TYPE 6 - CLAY82 - OCT. 1982 ASCE',/,
     1T20,'TYPE 7 - CLAY79 - AUG. 1979 ASCE',///)
121   FORMAT (A5,A60)
122   FORMAT (I2.2, '100', 15I5)
      IDT = '          '
C     OPEN (1,FILE='CON')
123   FORMAT (1X,4A)
      WRITE (*,120)
C---INITIAILIZE VALUES
C
C     INPUT FILE
C
29    WRITE (*,'(A)')' ENTER: INPUT FILE NAME ---->'
C      READ (*,'(A)') FLNAME
      FLNAME='D:INFILE.DSP'
  
      INQUIRE (FILE=FLNAME,EXIST=EXX)

      IF (EXX) THEN
                OPEN (8,FILE=FLNAME,STATUS='OLD')
      ELSE
                WRITE (*,'(3A)') ' FILE --', FLNAME,
     1    '-- NOT FOUND --->  PLEASE TRY AGAIN '
                GO TO 29
      END IF
C
C     OUTPUT FILE
C
      WRITE (*,'(A)')' ENTER: OUTPUT FILE NAME (PRN)            ---->'
      READ (*,'(A)') FLNAME
      
      IF (FLNAME.EQ.'                    ')
     1       FLNAME='PRN                 '

      INQUIRE (FILE=FLNAME,EXIST=EXX)
      IF (EXX) THEN

           OPEN (9,FILE=FLNAME,STATUS='OLD')
           REWIND 9
      ELSE
           OPEN (9, FILE=FLNAME,STATUS='NEW')
      END IF
C
C     PYCURVES FILE
C
      WRITE (*,'(A)')'+ENTER: BEAMCOL INPUT FILE NAME (PYCURVES.TXT)-->'
      READ (*,'(A)') FLNAME

      IF (FLNAME.EQ.'                    ')
     1       FLNAME='PYCURVES.TXT'

      INQUIRE(FILE=FLNAME,EXIST=EXX)
      IF (EXX) THEN

         OPEN(5,FILE=FLNAME,STATUS='OLD')
         REWIND 5
      ELSE
         OPEN (5,FILE=FLNAME,STATUS='NEW')
      END IF
C
C     SCRATCH FILE (UNIT 7)
C
      FLNAME='TAPE7.TXT'
     
      INQUIRE(FILE=FLNAME,EXIST=EXX)
      IF (EXX) THEN
         OPEN(7,FILE=FLNAME,STATUS='OLD')
         REWIND 7
      ELSE
         OPEN (7,FILE=FLNAME,STATUS='NEW')
      END IF

      IU8 = 8
      IU9 = 9
      IPROB = 1
C
C     'E' CAN BE MADE A DEFAULT FOR ALL DRILLIT PROBLEMS
C
      WRITE (*,'(A)')
     1          '+ENTER: MODULUS OF ELASTICITY (E) (3122018.6)->'
      READ(*,'(F20.0)') E
      IF (E.LE.0) E = 3122018.6

      WRITE (*,'(A)')
     1'+ENTER: INCREMENT LENGTH, (6) OR 12 INCHES --->'
      READ(*,'(I80)') INCR
      IF (INCR.LE.0) INCR = 6

      WRITE (*,'(A)')
     1          '+ENTER: DEFLECTION CLOSURE TOLERANCE (.001)--->'
      READ(*,'(F20.0)') DCTOL
      IF (DCTOL.LE.0.00001) DCTOL = .001
      IF (DCTOL.GT.0.2)     DCTOL = .001

1     CONTINUE
      READ(8,121,END=9991) IANS, TITLE
C---   SET VALUE OF NANS TO:
C---        0  FOR NORMAL
C---        1  FOR HOLD DEPTH AND VARY LOAD
C---        2  FOR HOLD LOAD AND VARY DEPTH
      NANS = 0
      IPRNT = 0
      IRANGE = 0
      IPIPE = 0
      THICKNESS=0
      KD = 1
      KLL = 1
      ILOADS = 0
      IF (IANS.EQ.'RPIPE') THEN
                IPIPE = 1
                IANS='RANGE'
      ENDIF
      IF (IANS.EQ.'RANGE') THEN
                NANS = 2
                IF (IPIPE.EQ.1) THEN
                        READ (8,*,END=9991) TO,THICKNESS
                        IF (THICKNESS.LE.0) IPIPE=0
                ELSE
                        READ (8,*,END=9991) TO
                ENDIF
                IRANGE = 2
      ELSE
                IF (IANS.EQ.'CURVE') IPRNT = 1
                IF (IANS.EQ.'LOADS') THEN
                        NANS = 1
                        KD=2
                        ILOADS=2
                ENDIF
                IF (IANS.EQ.'DEPTH') NANS = 2
      ENDIF

      IF (IPIPE.NE.1) THEN
                WRITE (5, 102) TITLE
      ELSE
                WRITE (5, 108) TITLE,THICKNESS
      ENDIF

      CALL INPUT1

      IF (IEND.EQ.1) GO TO 9991

      CALL SOILAVE (SAVE, PSF, K, INCR)

      CALL PYPNT (PYSTA, J, K, INCR)

      LL = IFIX(DEPTH*12./INCR)
      DO 302 L = 1, 5
302   IT(L) = 0
      IT(6) = 1
      IT(7) = 0
      IT(8) = 3
      IT(9) = J*3
      IT(10) = 25
      IT(11) = 0
310   CONTINUE
      IT(12) = LL/4
      IT(13) = LL/2
      IT(14) = IFIX(LL*3./4)
      IT(15) = LL
C---   WRITE ID CARD
      WRITE (5, 103) IPROB, DEPTH, DIAM, CAP, MT, HT, VT
C---   TABLE 1
      WRITE (5,122) IPROB, IT
      IF (IT(6).EQ.0) GO TO 304
C---    WRITE TABLE 2
      WRITE (5, 104) IPROB, LL,FLOAT(INCR), DCTOL
C---   WRITE TABLE 4
304   IF (IT(8).EQ.0) GO TO 305
      MOM = MT*12000.
      IF (ABS(CAP).GE..1) MOM = MOM + HT*CAP*12000.
      WRITE (5, 105) IPROB, HT*1000., MOM
      IF (IPIPE.EQ.1) THEN
                EI = E * 3.141592654 * ((DIAM*12)**4 -
     1                    (DIAM*12 - THICKNESS)**4) / 64
      ELSE
                EI = 0.049087 * (DIAM*12.)**4 * E
      ENDIF
      WRITE (5, 106)    IPROB, EI, VT*1000.
      WRITE (5,107)     IPROB, LL, EI
305   CONTINUE
      IF (IT(9).EQ.0)  GO TO 350
      IF (IPRNT.EQ.1) WRITE (9,'(1X,A1,15X,A,//,11X,A,//)') CHAR(12),
     1  ' P - Y CURVES DPA IBM  VERSION 2.4.1',TITLE
      CALL TABLE 5 (J, PYSTA, PSF, INCR, SAVE, IPRNT)
350   CONTINUE
      IPROB = IPROB + 1
      IF (IPROB.GT.99) IPROB = 0
      IF (NANS - 1) 400, 360, 370
C---   VARY LOAD
360   CONTINUE
C   .... NEW CARD ADDED GO TO 400
      GO TO 400
C
C     THE NEXT 3 LINES NEVER EXECUTE
C     AND THE VARIABLE ISIZEG IS USED
C     BUT NEVER DEFINED.
C
      KD = KD + 1
      IF (KD.GT.ISIZEG) GO TO 400
      GO TO 380
370   CONTINUE
      IF (IRANGE.GT.0) THEN
                DEPTH = DEPTH + INCR / 12.
                IF (DEPTH.GT.TO) GO TO 400
                IF (DEPTH.GT.100) GO TO 400
      ELSE
                KLL = KLL + 1
                IF (KLL.GE.10) GO TO 400
                DEPTH = DEPTH - INCR / 12.
      ENDIF
      LL = IFIX (DEPTH * 12. / INCR + .0001)
380   CONTINUE
      IT(5) = 1
      IT(7) = 0
      IT(9) = 0
      IF (NANS.EQ.1) IT(2) = 1
      IF (IT(2).EQ.1) IT(6) = 0
      IF (IT(4).EQ.1) IT(8) = 0
      GO TO 310
400   CONTINUE
      KD = 1
      IF (IPRNT.EQ.0) CALL BMCOL
      REWIND 7
      REWIND 5
      GO TO 1
9991  CONTINUE
      WRITE (9,'(1X,A1)') CHAR(12)
      END

C      $SUBTITLE: 'SUBROUTINE INPUT '


      SUBROUTINE INPUT1
      REAL MT
      INTEGER TYPE(20)
      COMMON/BK3/ MT,HT,VT,IPROB,H,D,CAP,Z(20),C(20),DENSE(20),
     1E50(20),PHI(20),H2O,NSTRA,IU8,TYPE,IDT
      COMMON /BK4/ IEND
      CHARACTER*8 IDT
      IEND = 0
      Z(1)=0.
      READ (IU8,*,END=250) MT, HT, VT
      READ (IU8,*,END=250) H, D, CAP
      READ (IU8,*,END=250) H2O,XNSTRA
      NSTRA = XNSTRA
      J=NSTRA+1
      DO 200 I=2,J
                 READ (IU8,*,END=250)  Z(I),C(I),DENSE(I),E50(I),TYPE(I)
150              PHI(I)=0.
                 IF (C(I).LT.10.) GO TO 200
                 PHI(I)=C(I)
                 C(I)=0.
200   CONTINUE
      RETURN
210   NSTRA=I-2
      RETURN
250   CONTINUE
      IEND = 1
      RETURN
      END

C     $SUBTITLE: 'SUBROUTINE PYPNT'

      SUBROUTINE PYPNT (PYSTA,  J, K, INCR)
      COMMON /BK3/ MT, HT, VT,IPROB, DEPTH, DIAM, CAP, Z(20), C(20),
     1DENSE(20), E50(20), PHI(20), H2O, NUMST, IU8, TYPE,IDT
      CHARACTER*8 IDT
      REAL MT
      INTEGER TYPE(20), PYSTA(31)
C---   CALCULATE P-Y STATIONS
      DO 260 I = 1, 31
260   PYSTA(I) = -1
      J=0
      II1 = NINT(2.*12./INCR)
      II2 = NINT(4.*12./INCR)
      IF (NUMST.LT.10) GO TO 265
      II1 = 500
      II2 = 1
265   CONTINUE
C---   SET UP STRATUMS
      DO 300 I = 2, K
C---            INCREMENT J COUNTER
                J = J+1
C---            TOP OF STRATUM
                PYSTA(J) = NINT(Z(I-1)*12./INCR)
                DZ = Z(I) - Z(I-1)
                IDZ = NINT(DZ*12./INCR)
                NPYSTA = (IDZ-II1)/II2
                IF (NPYSTA.LT.0) NPYSTA = 0
                IF (NPYSTA.EQ.0) GO TO 280
C---   NO INTERMEDIATE POINTS
C---   CALCULATE INTERMEDIATE POINTS
                IC = IDZ
                DO 270 L = 1, NPYSTA
                      J = J + 1
                      INC = IC/(NPYSTA - L + 2)
                      PYSTA(J) = PYSTA(J - 1) + INC
                      IC = IC - INC
270             CONTINUE
280             CONTINUE
C---   BOTTOM OF STRATUM
                J = J + 1
                PYSTA (J) = NINT(Z(I)*12./INCR)
                IF (PYSTA(J).LT.40)   GO TO 266
                II1 = 500
                II2 = 1
266             CONTINUE
300   CONTINUE
301   CONTINUE
      RETURN
      END

C     $SUBTITLE: 'SUBROUTINE TABLE 5 '


      SUBROUTINE TABLE 5 (J, PYSTA, PSF, INCR, SAVE, IPRNT)
      COMMON /BK3/ MT, HT, VT,IPROB, DEPTH, DIAM, CAP, Z(20), C(20),
     1DENSE(20), E50(20), PHI(20), H2O, NUMST, IU8, TYPE(20),IDT
      REAL MT
      CHARACTER IDT*8,CP(13)*5, CY(13)*5
      DIMENSION P(13), Y(13), PSF(200),SAVE(200)
      INTEGER TYPE, PYSTA(31)
108   FORMAT (I2.2, '5', I2.2, 3I5, 1P2E10.3, 2I5)
109   FORMAT (I2.2, '5', I2.2, 10X, 13A5)
112   FORMAT (//,' TYPE DEPTH DIAM   C   S AVE  PHI     WT',
     1'  E50 P-Y STA   J   I',/,
     17X,'(FT) (FT) (TSF) (KSF) (DEG)  (PCF)',/,
     1I5, F6.1, F5.1, F6.3, F6.3, F5.1, F7.1, F5.3,
     1I6, I6, I4,//)
113   FORMAT (/,7X,I5,12I7,/)
114   FORMAT (' P (LB/IN)', I3, 12I7)
115   FORMAT (' Y (IN)   ', F3.0, 13F7.3)
116   FORMAT (I2.2, '5', I2.2, 10X, 13A5)
117   FORMAT (2X,A5)
119   FORMAT (//)
120   FORMAT (G12.5)
C---   TABLE 5 PY CURVES
C---   INITIALIZE VALUES
      ISYM = 1
      NP = 0
      NCTB5 = 0
      KK = J
      I = 2
      J = 0
320   IF (J.GE.KK) GO TO 400
C---   DETERMINE DEFAULT P-Y TYPE
      J = J + 1
      AVERS = C(2) * 2.
      IF (PYSTA(J).GT.0) AVERS = SAVE(PYSTA(J))
      IF (J.LE.1)GO TO 325
      IF (PYSTA(J).EQ.PYSTA(J-1)) I = I + 1
325   X = PYSTA(J)*INCR/12.
      WT = 0.
      IF (X.GT.0.) WT = PSF(PYSTA(J))/X
      IF (TYPE(I).LT.0) TYPE(I) = 0
      IF (TYPE(I).GT.7) TYPE(I)=0
      IF (TYPE(I).NE.0) GO TO 350
      IF (C(I).LE.0.5) TYPE(I) = 6
      IF (C(I).GT.0.5) TYPE(I) = 6
      IF (PHI(I).GE.10.0) TYPE(I) = 3
350   CONTINUE
      CALL PY (TYPE(I), X, DIAM, C(I), AVERS, PHI(I), WT,
     1        E50(I), P, Y, NP, DEPTH,H2O)
      IF (IPRNT.EQ.1) THEN
                WRITE (9,112) TYPE(I),X,DIAM,C(I),AVERS,PHI(I),
     1                        WT, E50(I), PYSTA(J), J, I
                WRITE (9,113) (I2, I2 =1,NP)
                WRITE (9,114) (NINT(P(I2)), I2=1, NP)
                WRITE (9,115) (Y(I2), I2 = 1, NP)
                WRITE (9,119)
      ENDIF
      PM=1.
      YM = 1.
      YMAX = 0.
      PMAX = 0.
      DO 377 MM = 1, NP
      IF (P(MM).GT.PMAX) PMAX = P(MM)
377   IF (Y(MM).GT.YMAX) YMAX = Y(MM)
      YM = YMAX/10000.
      PM = PMAX/10000.
      IF (PM.LE.0.) PM = 1.
      IF (YM.LE.0.) YM = 1.
      NCTB5 = NCTB5 + 1
      ISTF = 0
      ISTT = 0
      ICONT = 0
C---   TOP
      IF (J.GE.KK) GO TO 379
      IF (PYSTA(J).NE.NINT(Z(I-1)*12./INCR)) GO TO 379
      ISTF = PYSTA(J)
      ICONT = 1
      GO TO 385
C--- BOTTOM
379   IF (PYSTA(J).NE.NINT(Z(I)*12./INCR)) GO TO 380
      ISTT = PYSTA(J)
      GO TO 385
C---   MIDDLE
380   ISTT = PYSTA(J)
      ICONT = 1
385   CONTINUE
      WRITE (5, 108) IPROB, NCTB5, ISTF, ISTT, ICONT, -1*PM*INCR,
     1YM, NP, ISYM
      NCTB5 = NCTB5 + 1
      REWIND 7
      WRITE (7, 120) (P(MM)/PM, MM = 1, NP), (Y(MM)/YM, MM = 1, NP)
      REWIND 7
      READ (7, 117) (CP(I2),I2=1,NP), (CY(I3),I3=1, NP)
C---   WRITE TABLE 5
      WRITE (5, 109) IPROB, NCTB5, (CP(MM), MM = 1, NP)
      NCTB5 = NCTB5 + 1
      WRITE (5, 116) IPROB, NCTB5, (CY(MM), MM = 1, NP)
      GO TO 320
400   CONTINUE
      RETURN
      END

C    $SUBTITLE: 'SUBROUTINE SOILAVE '


      SUBROUTINE SOILAVE (SAVE, PSF, K, INCR)
      COMMON /BK3/ MT, HT, VT,IPROB, DEPTH, DIAM, CAP, Z(20), C(20),
     1DENSE(20), E50(20), PHI(20), H2O, NUMST, IU8, TYPE,IDT
      REAL MT
      CHARACTER*8 IDT
      INTEGER TYPE(20)
      DIMENSION SAVE(200),PSF(200)
C---   INITIALIZE VALUES
      WATER = 0.
      SUM = 0.
      SUMC = 0.
      M = 0
      K = NUMST + 1
      H2O = FLOAT(NINT(H2O*12./INCR))/
     1  FLOAT(NINT(12./INCR))
C---   SOIL AVERAGING
C---   DOUBLE DO LOOPS 250 AND 240
      DO 250 I = 2, K
                Z(I) = FLOAT(NINT(Z(I)*12./INCR))/FLOAT(NINT(12./INCR))
                L = NINT(Z(I-1)*12./INCR) + 1
                M = NINT(Z(I)*12./INCR)
                TANGENT = TAN(PHI(I)/57.29577951)
                DO 240 N = L, M
                        IF (NINT(H2O*12./INCR).LT.N) WATER = 62.4
                        SUM = SUM + (DENSE(I) - WATER) / (12./INCR)
                        PSF(N) = SUM
                        SUMC = SUMC + C(I) + (PSF(N) -
     1                  (DENSE(I)-WATER))*INCR/24. * TANGENT/1000.
                        SAVE(N) = SUMC/(N/2.)
240             CONTINUE
250   CONTINUE
C---   END OF SOIL AVERAGING
      RETURN
      END

C     $SUBTITLE: 'SUBROUTINE PY '

      SUBROUTINE PY (ITYPE, H, DIAM, CC, SAVE, PHII, WTT, E50, P, Y, NP,
     1 D,H2O2)
C---REVISED 03-26-82 BY L.E. CADWALLADER
C---THIS PROGRAM COMPUTES SOIL RESISTANCE-DEFLECTION (P-Y) DATA POINTS
C---ACCORDING TO THE FOLLOWING PROCEDURES:
C---SOFT CLAY ABOVE OR BELOW WATER BY MATLOCK, 1970
C---STIFF CLAY ABOVE WATER BY REESE AND WELCH, 1975
C---STIFF CLAY BELOW WATER BY REESE, COX, AND KOOP, 1975
C---SAND ABOVE OR BELOW WATER BY REESE, COX, AND KOOP, 1974
      DIMENSION P(13), Y(13)
      REAL KO,KP,KA
      FX(A,B,C,D,X)=A*X*X*X+B*X*X+C*X+D
C---START OF SOFT CLAY ABOVE OR BELOW WATER AND STIFF CLAY ABOVE WATER
      GOTO (1, 2, 3, 4, 5, 6, 7) ITYPE



C     ENTRY CLAY79
7     XN=0.5
      XJ=2.0
      C1=2.0
      C2=0.5
      GO TO 210



C     ENTRY CLAY 82
6     XN = .5
      XJ = 2.83
      C1 = 2.0
      C2 = .5
      GO TO 210
C     ENTRY SOFT
1     XN = 1./3.
      GO TO 19



C     ENTRY STIFFA
2     XN = .25
19    CONTINUE
      XJ = .5
      C1 = 2.5
      C2 = 0.5
210   CONTINUE
      IF (CC.LT.0.001) GO TO 600
      A = (1. / C2)**(1. / XN)
      IF (E50.GT.0.0) GO TO 1001
      E50=.02
      IF (CC.GE..25) E50=.01
      IF (CC.GE..50) E50=.005
      IF (CC.GE.2.) E50=.004
1001  CONTINUE
C---CONVERT TO LBS. AND INCHES
      C=CC*2000./144.
      WT=WTT/1728.
      DIA=DIAM*12.
      DEPTH=H*12.
C---COMPUTE SOIL RESISTANCE P (LB/IN) AND DEFLECTION Y (IN)
      Q=3.+WT*DEPTH/C+XJ*DEPTH/DIA
      IF (Q.GT.9.) Q=9.
      PU=Q*C*DIA
      Y50= C1*E50*DIA
      YU=A*Y50
C---ESTABLISH 10 POINTS ON THE P-Y CURVE
      NP=10
      Y(1)=0.
      P(1)=0.
      Y(9)=YU
      P(9)=PU
      Y(10)=100.
      P(10)=PU
      DP=PU/8.
      DO 20 I=1,7
              P(I+1)=DP*I
              Y(I+1) = Y50 * ( P(I+1) / (PU*C2) )**(1. / XN)
20    CONTINUE
      RETURN



C---START OF STIFF CLAY BELOW WATER
C     ENTRY STIFFB
4     IF (E50.GT.0.) GO TO 1120
      IF (CC.LT.0.001) GO TO 600
      E50=.007
      IF (CC.GE.1.) E50=.005
      IF (CC.GE.2.) E50=.004
1120  K=500
      IF (CC.GE.1.) K=1000
      IF (CC.GE.2.) K=2000
C---COMPUTE AS
      XB=H/DIAM
      AS=FX (.015463838,-.126917197,.366979043,.2,XB)
      IF (XB.GE.3.64) AS=FX(0.,0.,.003349282,.587808612,XB)
C---CONVERT TO LBS. AND INCHES
      C=CC*2000./144.
      WT=WTT/1728.
      DIA=DIAM*12.
      DEPTH=H*12.
C---COMPUTE SOIL RESISTANCE P (LB/IN) AND DEFLECTION Y (IN)
      Q=2.+WT*DEPTH/C+2.83*DEPTH/DIA
      IF (Q.GT.11.) Q=11.
      PC=Q*C*DIA
      Y50=E50*DIA
      IF (DEPTH.EQ.0.) GO TO 29
      YK=(.5*PC/(SQRT(Y50)*K*DEPTH))**2
      PK=.5*PC*SQRT(YK/Y50)
29    CONTINUE
      YL=AS*Y50
      PL=.5*PC*SQRT(AS)
      W=3.*AS*Y50
      YP=W
      DO 23 I=1,4
      DP=.25*PC/SQRT(W*Y50)-.06875*PC*(((W/(AS*Y50))-1.)**.25)/(AS*Y50)
      DDP=-.125*PC/SQRT((W**3)*Y50)-.0171875*PC/(((W/(AS*Y50))-1.)**.75)
     1/(AS*Y50)**2
      W=W-DP/DDP
      IF(ABS((W-YP)/W).LE..001) GO TO 24
23    YP=W
24    YP=W
      PP=.5*PC*SQRT(YP/Y50)-.055*PC*((YP/(AS*Y50))-1.)**1.25
      YM=6.*AS*Y50
      PM=.5*PC*SQRT(6.*AS)-.411*PC
      YU=18.*AS*Y50
      PU=.5*PC*SQRT(6.*AS)-.411*PC-.75*PC*AS
      IF (PU.LT.0.) PU=.05*PP
C---ESTABLISH 10 POINTS ON THE P-Y CURVE
      NP=10
      Y(1)=0.
      P(1)=0.
      IF (DEPTH.EQ.0.) GO TO 30
      IF (YK.GE.YL) GO TO 30
      Y(2)=YK
      P(2)=PK
      Y(5)=YL
      P(5)=PL
      Y(6)=YL+(YP-YL)/3.
      P(6)=.5*PC*SQRT(Y(6)/Y50)-.055*PC*((Y(6)/(AS*Y50))-1.)**1.25
      Y(7)=YP
      P(7)=PP
      Y(8)=YM
      P(8)=PM
      Y(9)=YU
      P(9)=PU
      Y(10)=100.
      P(10)=PU
      DP=(PL-PK)/3.
      DO 22 I=1,2
              P(I+2)=PK+DP*I
              Y(I+2)=Y50*(P(I+2)*2./PC)**2
22    CONTINUE
      RETURN
30    CONTINUE
C---ESTABLISH 7 POINTS ON THE P-Y CURVE
      NP=7
      Y(2)=YL
      P(2)=PL
      Y(3)=YL+(YP-YL)/3.
      P(3)=.5*PC*SQRT(Y(3)/Y50)-.055*PC*((Y(3)/(AS*Y50))-1.)**1.25
      Y(4)=YP
      P(4)=PP
      Y(5)=YM
      P(5)=PM
      Y(6)=YU
      P(6)=PU
      Y(7)=100.
      P(7)=PU
      RETURN



C     ENTRY SAND
3     IF (PHII.LT.0.001) GO TO 600
      IF (ITYPE.EQ.3) GO TO 1005
C---START OF SAND ABOVE OR BELOW WATER

C---Estimate K
C--- Above the water table
      IF (H.GT.H2O2) GO TO 1005
      K=25
      IF (PHII.GE.30.) K=90
      IF (PHII.GE.35.) K=225
      GO TO 12

1005  CONTINUE

C--- Below the water table
      K=20
      IF (PHII.GE.30.) K=60
      IF (PHII.GE.35.) K=125



12    CONTINUE
      IF (H.LE.0.) GO TO 40
C---COMPUTE AS,BS
      XB=H/DIAM
      AS=FX(0.,.133333333,-1.033333333,3.,XB)
      IF (XB.GE.3.8) AS=FX(0.,0.,0.,.88,XB)
      BS=FX(-.066666667,.08,-.543333333,2.11,XB)
      IF (XB.GE.1.5) BS=FX(.066666667,-.4,.383333333,1.35,XB)
      IF (XB.GE.3.0) BS=FX(.006666667,.03,-.666666667,2.25,XB)
      IF (XB.GE.5.0) BS=FX(0.,0.,0.,.5,XB)
C---CONVERT TO RADIANS, LBS., AND INCHES
      PHI=PHII/57.29578
      WT=WTT/1728.
      DIA=DIAM*12.
      DEPTH=H*12.
C---COMPUTE FUNCTIONS OF PHI ANGLE
      TANA=TAN(PHI)
      TANB=TAN(PHI/2.)
      TANC=TAN(.7854+PHI/2.)
      TAND=TAN(.7854-PHI/2.)
      SINC=SIN(.7854+PHI/2.)
      COSB=COS(PHI/2.)
      KO=.4
      KP=TANC**2
      KA=TAND**2
C---COMPUTE SOIL RESISTANCE P (LB/IN) AND DEFLECTION Y (IN)
      PCS=WT*DEPTH*(KO*DEPTH*TANA*SINC/(TAND*COSB)+KP*(DIA+DEPTH*TANC
     1*TANB)+KO*DEPTH*TANC*(TANA*SINC-TANB)-KA*DIA)
      PCD=WT*DEPTH*DIA*(KA*((TANC**8)-1.)+KO*TANA*(TANC**4))
      PC=PCS
      IF (PC.GT.PCD) PC=PCD
      YU=3.*DIA/80.
      PU=AS*PC
      YM=DIA/60.
      PM=BS*PC
      S=(PU-PM)/(YU-YM)
      E=PM/(S*YM)
      C=PM/YM**(1./E)
      YK=(C/(K*DEPTH))**(E/(E-1.))
      IF (YK.GE.YM) YK=DIA/180.
      PK=C*YK**(1./E)
C---ESTABLISH 10 POINTS ON THE P-Y CURVE
      NP=10
      Y(1)=0.
      P(1)=0.
      Y(2)=YK
      P(2)=PK
      Y(8)=YM
      P(8)=PM
      Y(9)=YU
      P(9)=PU
      Y(10)=100.
      P(10)=PU
      DP=(PM-PK)/6.
      DO 27 I=1,5
      P(I+2)=PK+DP*I
      Y(I+2)=(P(I+2)/C)**E
27    CONTINUE
      RETURN
40    CONTINUE
600   CONTINUE
C---ESTABLISH 2 POINTS ON THE P-Y CURVE
      NP=2
      Y(1)=0.
      P(1)=0.
      Y(2)=100.
      P(2)=0.
      IF (ITYPE.EQ.5) GO TO 601
      RETURN



C     ENTRY PADLL
5     IF ((CC + PHII).LT.0.001) GO TO 600
      IF (E50.LT.0.001) GO TO 600
      C   =CC*2.
      PHI=PHII
      B=DIAM
      EP=E50
      ES=EP*144.
      PI=3.141592654
      RAD=PI/180.
      ZOB=H/B
      CONST = 4. / (3. * PI)
      QBAR = WTT * H / 1000.
      SRF = 1
      XKH=ES*5.7/(B*(D/B)**0.4)
C.....COMPUTE TRIGONOMETRIC FUNCTIONS
      SINPHI=SIN(PHI      *RAD)
      COSPHI=COS(PHI      *RAD)
      TANPHI=SINPHI/COSPHI
      IF(PHI      .GT.0.0) COTPHI=COSPHI/SINPHI
      SIN45P=SIN((45.0+PHI      /2.0)*RAD)
      COS45P=COS((45.0+PHI      /2.0)*RAD)
      TAN45P=SIN45P/COS45P
      SIN45M=SIN((45.0-PHI      /2.0)*RAD)
      COS45M=COS((45.0-PHI      /2.0)*RAD)
      TAN45M=SIN45M/COS45M
C
C.....COMPUTE COEFFICIENTS IN BRINCH HANSEN'S ULTIMATE CAPACITY FORMULA
C
      XK0=1.0-SINPHI
      DCINF=1.58+4.09*TANPHI**4
      XNC=5.14
      IF(PHI      .GT.0.0) XNC=(EXP(PI*TANPHI)*TAN45P**2-1.0)*COTPHI
      XKCINF=XNC*DCINF
      XKQINF=XKCINF*XK0*TANPHI
      A1=EXP((PI/2.0+PHI      *RAD)*TANPHI)*COSPHI*TAN45P
      XKC0=2.57
      IF(PHI      .GT.0.0) XKC0=(A1-1.0)*COTPHI
      XKQ0=A1-EXP(-(PI/2.0-PHI      *RAD)*TANPHI)*COSPHI*TAN45M
      AC=2.0*SIN45P*XKC0/(XKCINF-XKC0)
      AQ=0.0
      IF(PHI      .GT.0.0) AQ=XK0*SINPHI*XKQ0/((XKQINF-XKQ0)*SIN45P)
      XKC=(XKC0+XKCINF*AC*ZOB)/(1.0+AC*ZOB)
      XKQ=(XKQ0+XKQINF*AQ*ZOB)/(1.0+AQ*ZOB)
C
C.....APPLY STRENGTH REDUCTION FACTOR TO LAYER STRENGTH PROPERTIES
C
      CA=SRF      *C
      TANPHA=SRF      *TANPHI
C
C.....COMPUTE ULTIMATE LATERAL PRESSURE AND SIDE SHEAR MOMENT FOR EACH
C     ELEMENT IN SOIL-PIER MODEL
C
      FR=(C      *XKC+QBAR*XKQ-PI*CA/4.0)*B     /(1.0+CONST*TANPHA)
      VRTHET=CA*B     *PI/4.0+CONST*FR*TANPHA
      PULT      =(FR+VRTHET)/B
      IF (PULT.LT.0.0001) GO TO 600
      P(1)=0.
      Y(1)=0.
      DO 700 I=1,10
      P(I+1)=PULT*I*0.1*B/12.
      P(I+1)=P(I+1)*1000.
      Y(I+1)=PULT*(I*0.1)**2/(0.72*XKH)
      Y(I+1)=Y(I+1)*12.
 700  CONTINUE
      P(12)=P(11)
      Y(12)=100.
      NP=12
601   CONTINUE
C     WRITE (9,100) AC, AQ, XKC, XKQ, ZOB, XKC0, XKQ0, XKCINF, XKQINF,
C    1DCINF, XK0,XNC, QBAR, PULT
100   FORMAT (' AC', T11,F11.5, T30, 'AQ', T41, F11.5,/
     1' XKC   ', T11, F11.5, T30, 'XKQ   ', T41, F11.5,/
     1' ZOB ', T11, F11.5,/
     1' XKC0  ', T11, F11.5, T30, 'XKQ0  ', T41, F11.5,/
     1' XKCINF', T11, F11.5, T30, 'XKQINF', T41, F11.5,/
     1' DCINF ', T11, F11.5, T30, 'XK0   ', T41, F11.5,/
     1' XNC ', T11, F11.5, T30, 'QBAR  ', T41, F11.5,/
     1' PULT', T11, F11.5)
      RETURN
      END

C    $SUBTITLE: 'FUNCTION REDFAC '


      FUNCTION REDFAC (C, PHI)
      REDFAC=1.
      IF (C.GT.0.5) GO TO 700
      IF (PHI.GT.35.) REDFAC= 1. - 0.02 * (PHI - 35.)
      RETURN
 700  IF (C.GT.1.) GO TO 710
      REDFAC=1.25- 0.5*C
      RETURN
 710  IF (C.GT.2.) GO TO 720
      REDFAC=1.05-0.3*C
      RETURN
 720  IF (C.GT.2.5) GO TO 730
      REDFAC=0.65-0.1*C
      RETURN
 730  REDFAC=0.4
      RETURN
      END

C    $SUBTITLE:'SUBROUTINE BMCOL28'


      SUBROUTINE BMCOL
      COMMON /RANGE/ IRANGE,INCR,THICKNESS,E11,EI11,IPIPE
      COMMON /LOADS/ ILOADS,KD
      COMMON /BK3/ MT, HT, VT,IPROB, DEPTH, DIAM, CAP, Z(20), COH(20),
     1         DENSE(20), E50(20), PHI(20), H2O, NUMST, IU8, TYPE, IDT
      COMMON /BKBC/ F, Q, QF, S, SF, T, R, P, W, A, B, C, PETER, QFSAVE,
     1              PSAVE, TSAVE
      REAL MT, LDFAC, LDFAC1, LDFAC2
      SAVE /RANGE/
      DIMENSION    XXX(507), DEFLCT(507), ENMONT(507),
     1             F(507), Q(507), QF(507), S(507), SF(507), T(507),
     2             R(507), P(507), W(507), KEY(507), WS(30), DWS(30),
     3             KR(30), NPTS(30), JC1(30), JC2(30), QVN(15),
     4             WVN(15), QV(30,30), WV(30,30), IM(5), WM(5),
     5             A(507), B(507), C(507), ROTATE(507),
     6             QFSAVE(507), PSAVE(507), TSAVE(507)
      DOUBLE PRECISION  A, AA, AREV, ATEMP, B, BB, BREV, BTEMP, C,
     1                  CC, CREV, CTEMP, D, DBM, DD, DENOM, DREV, DTEMP,
     2                  DW, E, EE, F, FF, P, PART, Q, QF, Q1, Q2, R,
     3                  REACT, S, SF, SUM1, SUM2, S1, S2, T, W, WM, WP,
     4                  QFSAVE, PSAVE, TSAVE
      CHARACTER AN1(48)*5, AN2(14)*5,IDT*8,ERR*12
      INTEGER TYPE(20)
   12 FORMAT ( 16A5 )
   14 FORMAT ( I2, 5X, I1, 2X, 14A5 )
   20 FORMAT ( 5X, 15I5 )
   21 FORMAT ( 5X, I5, 10X,   3E10.3 )
   31 FORMAT ( 2(5X, I5),   2E10.3 )
   41 FORMAT ( 5X, 3I5,   6E10.3 )
   51 FORMAT ( 5X, 3I5,   2E10.3, 2I5 )
   52 FORMAT ( 15X, 13F5.0 )
9010  FORMAT (4X,60('*'),//,36X, 'OUTPUT', /,36X, '======', //,
     17X, 'DEFLECTION TOP =', F12.4, ' INCHES ',
     12X, 'BOTTOM = ', F12.4, ' INCHES', //,
     17X, 'ROTATION   TOP =', F12.4, ' DEGREES',
     12X, 'BOTTOM = ', F12.4, ' DEGREES', //,
     17X, 'MOMENT     TOP =', F12.4, ' IN-K   ',
     12X, 'BOTTOM = ', F12.4, ' IN-K', //,
     17X, 'MOMENT     MAX =', F12.4, ' IN-K.....',
     1    '......AT ', F12.4, ' FT', //,
     17X, 'MOMENT     MIN =', F12.4, ' IN-K.....',
     1    '......AT ', F12.4, ' FT', ///,
     125X, 'TOTAL APPLIED LOADS = ', 1PE12.5,//,
     125X, 'TOTAL ALL LOADS     = ', 1PE12.5)
      REWIND 5
      XMT = 0.
      XHT = 0.
C-----START EXECUTION OF PROGRAM - SEE GENERAL FLOW CHART, APPENDIX 1.1
1000  CONTINUE
C-----PROGRAM AND PROBLEM IDENTIFICATION
      READ (5,12) ( AN1(N), N = 1, 48 )
 1010 READ (5,14,END=9999)  NPROB, IPLOT,  ( AN2(N), N = 1, 14 )
      LDFAC = 1.
      LDFAC1 = 0.
      LDFAC2 = 1.
      LDLOCK = 0
C         IF ( NPROB ) 1020,  9999, 1020
1020  CONTINUE
C-----INPUT TABLE 1
 1100 READ (5,20) KEEPW, KEEP2, KEEP3, KEEP4, KEEP5, NCD2, NCD3, NCD4,
     1         NCD5, ITMAX, ( IM(N), N = 1, 5 )
          IF ( ITMAX ) 9999, 1200, 1110
1110  CONTINUE
C-----INPUT TABLE 2
1200  CONTINUE
          IF ( KEEP2 ) 9999, 1210, 1230
 1210 READ (5,21) M, H, WMAX, WTOL
          IF ( ITMAX ) 9999, 1240, 1215
1215  CONTINUE
          GO TO 1240
1230  CONTINUE
C-----COMPUTE CONSTANTS AND INDEXES
 1240          HT2 = H + H
               HE2 = H * H
               HE3 = H * HE2
               MP3 = M + 3
               MP4 = M + 4
               MP5 = M + 5
               MP6 = M + 6
               MP7 = M + 7
C-----INPUT TABLE 3
1300  CONTINUE
          IF ( KEEP3 ) 9999, 1310, 1305
1305  CONTINUE
          GO TO 1399
 1310     DO 1315 J = 1, MP7
               KEY(J) = 1
 1315     CONTINUE
          IF ( NCD3 ) 9999, 1320, 1325
1320  CONTINUE
          GO TO 1399
 1325          IF ( NCD3 - 20 ) 1327, 1327, 1326
1326  CONTINUE
          GO TO 9999
 1327          JS = 3
          IF (1.GT.NCD3) GO TO 1351
 1331     DO 1350 N = 1, NCD3
      READ (5,31) IN1, KASE, WS(N), DWS(N)
          IF ( IN1 + 4 - JS ) 1328, 1328, 1329
1328  CONTINUE
          GO TO 9999
 1329          JS = IN1 + 4
C-----SET INDEXES FOR FUTURE CONTROL OF SPECIFIED CONDITION ROUTINES
          IF (KASE-2) 1330,1335,1340
 1330          KEY(JS) = 2
          GO TO 1350
 1335          KEY(JS-1) = 3
               KEY(JS+1) = 5
          GO TO 1350
 1340          KEY(JS-1) = 3
               KEY(JS) = 4
               KEY(JS+1) = 5
 1350     CONTINUE
 1351     CONTINUE
 1399     CONTINUE
C-----INPUT TABLE 4
1400  CONTINUE
          IF ( KEEP4 ) 9999, 1401, 1403
1401      DO 1402 J = 1, MP7
               F(J) = 0.0
               QF(J) = 0.0
               SF(J) = 0.0
               T(J) = 0.0
               R(J) = 0.0
               P(J) = 0.0
1402     CONTINUE
          GO TO 1404
1403  CONTINUE
1404     IF ( NCD4 ) 9999, 1405, 1406
1405  CONTINUE
          GO TO 1499
C-----SEE FLOW CHART, APPENDIX 1.2 AND 1.3, TABLE 4 SEQUENCING PROCEDURE

1406          KR2 = 0
          DO 1480 N = 1, NCD4
               KR1 = KR2
      READ (5,41) IN1, IN2, KR2, FN2, QN2, SN2, TN2, RN2, PN2
      IF (QN2.GT.0.001) XHT = QN2
      IF (TN2.GT.0.001)  XMT = TN2
               JN = IN1 + 4
               J2 = IN2 + 4
               KSW = 1 + KR2 + 2 * KR1
          GO TO ( 1407, 1410, 1415, 1415 ), KSW
1407  CONTINUE
          GO TO 1420
1410  CONTINUE
          GO TO 1420
1415  CONTINUE
          GO TO 1435
 1420          J1 = JN
 1425          FN1 = FN2
               QN1 = QN2
               SN1 = SN2
               TN1 = TN2
               RN1 = RN2
               PN1 = PN2
          GO TO ( 1435, 1480, 9999, 1480 ), KSW
C-----SEE FLOW CHART, APPENDIX 1.3, TABLE 4 INTERPOL AND DISTRIB
 1435          JINCR = 1
               ESM = 1.0
          IF ( J2 - J1 ) 1437, 1450, 1440
1437  CONTINUE
          GO TO 9999
 1440          DENOM = J2 - J1
               ISW = 1
          GO TO 1455
 1450          DENOM = 1.0
               ISW = 0
 1455     DO 1460  J = J1, J2, JINCR
               DIFF = J - J1
               PART = DIFF / DENOM
               F(J) = F(J) + ( FN1 + PART * ( FN2 - FN1 ) ) * ESM
               QF(J) = QF(J) + ( QN1 + PART * ( QN2 - QN1 ) ) * ESM
               SF(J) = SF(J) + ( SN1 + PART * ( SN2 - SN1 ) ) * ESM
               T(J) = T(J) + ( TN1 + PART * ( TN2 - TN1 ) ) * ESM
               R(J) = R(J) + ( RN1 + PART * ( RN2 - RN1 ) ) * ESM
               P(J) = P(J) + ( PN1 + PART * ( PN2 - PN1 ) ) * ESM
 1460     CONTINUE

          IF ( ISW ) 9999, 1470, 1465
 1465          JINCR = J2 - J1
               ESM = - 0.5
               ISW = 0
          GO TO 1455
 1470     GO TO ( 1480, 9999, 1480, 1475 ), KSW
 1475          J1 = J2
          GO TO 1425
 1480     CONTINUE
          IF ( F(MP5) + F(MP6) + QF(MP5) + SF(MP5) + T(MP5) + T(MP6 ) +
     1      R(MP5) + R(MP6) + P(MP5) + P(MP6) ) 1485, 1499, 1485
1485  CONTINUE
          GO TO 9999
 1499     CONTINUE
C-----INPUT TABLE 5
1500  CONTINUE
C-----CLEAR OR HOLD PRIOR TABLE 5 DATA AND SET UP INPUT CONTROLS
          IF ( KEEP5 ) 9999, 1505, 1510
 1505          KCD5 = NCD5
               KCD = KCD5 / 3
               KST  = 1
          GO TO 1515
 1510  CONTINUE
               KST = KCD5 / 3 + 1
               KCD5 = KCD5 + NCD5
               KCD = KCD5 / 3
 1515     IF ( KCD - 30 ) 1524, 1524, 1520
 1520  CONTINUE
          GO TO 9999
 1524     IF ( NCD5 ) 9999, 1525, 1526
 1525  CONTINUE
          GO TO 1599
C-----INPUT Q-W CURVES
 1526          KR2 = 0
          IF (KST.GT.KCD) GO TO 1599
 1527     DO 1585 NC = KST, KCD
               KR1 = KR2
      READ (5,51) IN1, IN2, KR(NC), QBASE, WBASE, NPTS(NC), KSYM
               KR2 = KR(NC)
               NPT = NPTS(NC)
               KSW = 1 + KR2 + 2 * KR1
               JC1(NC) = IN1 + 4
               JC2(NC) = IN2 + 4
          GO TO ( 1530, 1535, 1540, 1540 ), KSW
1530  CONTINUE
          GO TO 1545
1535  CONTINUE
          GO TO 1545
1540  CONTINUE
 1545 READ (5,52) ( QVN(NP), NP = 1, NPT )
      READ (5,52) ( WVN(NP), NP = 1, NPT )
          DO 1550  NP = 1, NPT
               QV(NC,NP) = QVN(NP) * QBASE
               WV(NC,NP) = WVN(NP) * WBASE
 1550     CONTINUE
          IF ( NPT - 1 ) 1566, 1566, 1555
 1555     DO 1565 NP = 2, NPT
          IF ( WV(NC,NP) - WV(NC,NP-1) ) 1560, 1565, 1565
C-----INTERPOL ROUTINE AVOIDS DIV BY ZERO -- SEE APP 1.7
1560  CONTINUE
          GO TO 9999
 1565     CONTINUE
 1566     CONTINUE
          IF ( KSYM ) 9999, 1585, 1570
C-----APPLY SYMMETRY OPTION TO CONSTRUCT OTHER HALF OF Q-W CURVE
 1570          NP1 = NPT + 1
               NP2 = NPT * 2
          IF (NP1.GT.NP2) GO TO 1576
 1571     DO 1575  NP = NP1, NP2
               N = NP - NPT
               QV(NC,NP) = QV(NC,N)
               WV(NC,NP) = WV(NC,N)
 1575     CONTINUE
 1576     CONTINUE
          DO 1580  NP = 1, NPT
               N = NP2 + 1 - NP
               QV(NC,NP) = - QV(NC,N)
               WV(NC,NP) = - WV(NC,N)
 1580     CONTINUE
               NPTS(NC) = NP2
 1585     CONTINUE
 1599     CONTINUE
C-----CLEAR OR HOLD DEFLECTIONS W(J) FROM PRIOR PROBLEM
 2000     IF ( KEEPW ) 9999, 2010, 2025
 2010     DO 2020 J = 1, MP7
               W(J) = 0.0
 2020     CONTINUE
C-----PRINT MONITOR DATA HEADING ( TABLE 6 )
2025  CONTINUE
          IF ( ITMAX ) 9999, 2050, 2060
2050  CONTINUE
               NIT = 1
          GO TO 2200
2060  CONTINUE
               NIT = ITMAX
      CALL OUT
C-----START ITERATIVE SOLUTION
 2200     DO 7100 NITR = 1, NIT
C-----RESET FORCES Q AND SPRINGS S TO FIXED VALUES OF TABLE 4
          DO 2400  J = 3, MP5
               Q(J) = QF(J)
               S(J) = SF(J)
 2400     CONTINUE
C-----INTERPOLATE FOR Q AND S FROM Q-W CURVES
 3000          KOFFC = 0
          IF ( KCD5 ) 9999, 3999, 3001
C-----SEE FLOW CHARTS, APPENDIX 1.5, 1.6, 1.7
 3001          KR1 = 0
               KCD = KCD5 / 3
          DO 3095 NC = 1, KCD
          IF ( KR1 ) 9999, 3005, 3010
 3005          NC1 = NC
          IF ( KR(NC) ) 9999, 3010, 3090
C-----SEE FLOW CHARTS, APPENDIX 1.6 AND 1.7, INTERPOL BETWEEN CURVES
 3010          J1 = JC1(NC1)
               J2 = JC2(NC)
               DENOM = J2 - J1
               JEND = J2
               JINCR = 1
               ESM = 1.0
               ISW = 1
          IF ( DENOM )  3012, 3013, 3014
3012  CONTINUE
          GO TO  9999
 3013     IF ( MP4 - J1 ) 3080, 3020, 3020
 3014     IF ( MP4 - J1 ) 3080, 3080, 3015
 3015     IF ( MP4 - J2 ) 3016, 3030, 3030
 3016          JEND = MP4
          GO TO 3030
 3020          DENOM = 1.0
               ISW = 0
 3030     DO 3070  J = J1, JEND, JINCR
               L = 1
               N = NC1
C-----SEE FLOW CHART, APPENDIX 1.7, INTERPOL ALONG SINGLE CURVE
 3035          NPT = NPTS(N)
          IF (NPT.LE.2) GO TO 3042
 3037     DO 3040 NP = 2, NPT
          IF ( W(J) - WV(N,NP) ) 3045, 3055, 3040
 3040     CONTINUE
 3042          NP = NPTS(N)
               GO TO 3050
 3045     IF ( W(J) - WV(N,1) ) 3050, 3055, 3055
 3050          KOFFC = 1
 3055          NP = NP - 1
          IF ( WV(N,NP + 1) - WV(N,NP) ) 3057, 3058, 3057
 3057          S2 = - ( QV(N,NP+1) - QV(N,NP) ) / ( WV(N,NP+1) -
     1              WV(N,NP) )
               Q2 = QV(N,NP) + S2 * WV(N,NP)
          GO TO 3059
 3058          S2 = 0
               Q2 = 0
 3059     IF ( L ) 9999, 3065, 3060
 3060          Q1 = Q2
               S1 = S2
               L = 0
               N = NC
          GO TO 3035
 3065          DIFF = J - J1
               PART = DIFF / DENOM
               Q(J) = Q(J) + ( Q1 + PART * ( Q2 - Q1 ) ) * ESM
               S(J) = S(J) + ( S1 + PART * ( S2 - S1 ) ) * ESM
 3070     CONTINUE
          IF ( ISW ) 9999, 3080, 3075
 3075          JINCR = JEND - J1
               ESM = - 0.5
               ISW = 0
          GO TO 3030
 3080     IF ( KR(NC) ) 9999, 3090, 3085
 3085          JC1(NC) = JC2(NC)
               NC1 = NC
 3090          KR1 = KR(NC)
 3095     CONTINUE
 3999     CONTINUE
C-----START OF BEAM-COLUMN SOLUTION
 6000          NS = 1
               A(1) = 0
               A(2) = 0
               B(1) = 0
               B(2) = 0
               C(1) = 0
               C(2) = 0
          DO 6060 J = 3, MP5
C-----COMPUTE MATRIX COEFFS AT EACH STA J
           AA = F(J-1) - 0.25 * H * ( R(J-1) + H * P(J-1) )
           BB = - 2.0 * ( F(J-1) + F(J) )
           CC = F(J-1) + 4.0 * F(J) + F(J+1) + HE3 * S(J) +
     1          0.25 * H * ( ( R(J-1) + H * P(J-1) ) + ( R(J+1)
     2          + H * P(J+1) ) )
           DD = - 2.0 * ( F(J) + F(J+1) )
           EE = F(J+1) - 0.25 * H * ( R(J+1) + H * P(J+1) )
           FF = HE3 * Q(J) - 0.5 * HE2 * ( T(J-1) - T(J+1) )
C-----COMPUTE RECURSION OR CONTINUITY COEFFS AT EACH STA
               E = AA * B(J-2) + BB
               DENOM = E * B(J-1) + AA * C(J-2) + CC
          IF ( DENOM ) 6010, 6005, 6010
C-----NOTE IF DENOM IS ZERO, BEAM DOES NOT EXIST, D = 0 SETS DEFL = 0.
 6005          D = 0.0
          GO TO 6015
 6010          D = - 1.0 / DENOM
 6015          C(J) = D * EE
               B(J) = D * ( E * C(J-1) + DD )
               A(J) = D * ( E * A(J-1) + AA * A(J-2) - FF )
C-----CONTROL RESET ROUTINES FOR SPECIFIED CONDITIONS
               KEYJ = KEY(J)
          GO TO ( 6060, 6020, 6030, 6020, 6050 ), KEYJ
C-----RESET FOR SPECIFIED DEFLECTION
 6020          C(J) = 0.0
               B(J) = 0.0
               A(J) = WS(NS)
          IF ( KEYJ - 3 ) 6059, 6030, 6060
C-----RESET FOR SPECIFIED SLOPE AT NEXT STA
 6030          DTEMP = D
               CTEMP = C(J)
               BTEMP = B(J)
               ATEMP = A(J)
               C(J) = 1.0
               B(J) = 0.0
               A(J) = - HT2 * DWS(NS)
          GO TO 6060
C-----RESET FOR SPECIFIED SLOPE AT PRECEDING STATION
 6050          DREV = 1.0 / ( 1.0 - ( BTEMP * B(J-1) + CTEMP - 1.0 ) *
     1                D / DTEMP )
               CREV = DREV * C(J)
               BREV = DREV * ( B(J) + ( BTEMP * C(J-1) ) * D / DTEMP )
               AREV = DREV * ( A(J) + ( HT2 * DWS(NS) + ATEMP + BTEMP
     1                * A(J-1) ) * D / DTEMP )
               C(J) = CREV
               B(J) = BREV
               A(J) = AREV
 6059          NS = NS + 1
 6060     CONTINUE
C-----COMPUTE DEFLECTIONS  (NOTE THAT IF STIFFNESS VALUES  F, P, R  ARE
C         ZERO BEYOND END (M+4) OF REAL BMCOL,  B(MP5)=C(MP5)=C(MP4)=0.
C         AND VALUES OF  W(MP6)  AND  W(MP7)  ARE THEREFORE MOOT.
               W(MP6) = 0
               W(MP7) = 0
 6064          W(MP5) = A(MP5) + B(MP5) * W(MP6) + C(MP5) * W(MP7)
               KWCNT = 0
               KWMAX = 0
          DO 6100 L = 4, MP4
               J = M + 8 - L
               WP = W(J)
               W(J) = A(J) + B(J) * W(J+1) + C(J) * W(J+2)
C-----COUNT STAS WHERE DEFL IS EXCESSIVE
          IF ( WMAX - DABS( W(J) ) ) 6065, 6070, 6070
 6065          KWMAX = KWMAX + 1
C-----COUNT STAS NOT MEETING CLOSURE TOLERANCE
 6070     IF ( DABS( W(J) - WP ) - WTOL ) 6100, 6100, 6075
 6075          KWCNT = KWCNT + 1
 6100     CONTINUE
               W(3) = A(3) + B(3) * W(4) + C(3) * W(5)
          IF ( ITMAX ) 9999, 7200, 7005
C-----PRINT ITERATION MONITOR DATA
 7005     DO 7010 N = 1, 5
               JM = IM(N) + 4
               WM(N) = W(JM)
 7010     CONTINUE
          IF ( KOFFC ) 9999, 7015, 7020
7015  WRITE (*,611) NITR,KWCNT,(WM(N),N=1,5)
          GO TO 7025
7020  WRITE (*,612) NITR,KWCNT,(WM(N),N=1,5)
611   FORMAT ('+',2X,I4,7X,'NO',7X,I4,2X,5F10.6)
612   FORMAT ('+',2X,I4,7X,'YES',6X,I4,2X,5F10.3)
C-----CONTROL ITERATION PROCESS
 7025     IF ( KWMAX ) 9999, 7035, 7999
 7035     IF ( KWCNT ) 9999, 7999, 7100
 7100     CONTINUE
C-----7100 IS END OF MAIN ITERATION DO-LOOP
 7200          NITR = ITMAX
 7999     CONTINUE
C-----OUTPUT TABLE 7, RESULTS
         ERR=' '
         IF ( KOFFC ) 9999, 8010, 8005
8005  IF (IRANGE+ILOADS.EQ.0) THEN
          WRITE(9,900)
      ELSE
         ERR='OFF CURVES'
      ENDIF
900   FORMAT(/,'     ***** DEFLECTION OFF Q - W CURVES *****')
 8010     IF ( ITMAX ) 9999, 8100, 8015
 8015     IF ( KWMAX ) 9999, 8025, 8020
8020  IF (IRANGE+ILOADS.EQ.0) THEN
         WRITE (9,902) KWMAX
      ELSE
         ERR='DEF.EXCEEDED'
      ENDIF
902   FORMAT (/,'    ***** MAX ALLOWABLE DEFLECTION EXCEEDED AT ***',
     1 I4,'  *** STATIONS *****')
 8025     IF ( KWCNT ) 9999, 8100, 8030
8030  IF (IRANGE+ILOADS.EQ.0) THEN
         WRITE (9,901)
      ELSE
         ERR='NOT CLOSED'
      ENDIF
901   FORMAT (/,'    ***** SOLUTION NOT CLOSED WITHIN SPECIFIED',
     1  ' TOLERANCE *****')
C-----COMPUTE AND PRINT,RESULTS
8100  CONTINUE
               W(2) = 2.0 * W(3) - W(4)
               W(M+6) = 2.0 * W(MP5) - W(MP4)
               BMJ = 0.0
               BMK = 0.0
          DO 8200 J = 3, MP5
               ISTA = J - 4
               ZI = ISTA
               X = ZI * H
               DW = ( - W(J-1) + W(J+1) ) / HT2
               BMI = BMJ
               BMJ = BMK
               BMK = F(J+1) * ( W(J) - 2.0 * W(J+1) + W(J+2) ) / HE2
               DBM = ( - BMI + BMK ) / HT2
               REACT = ( BMI - 2.0 * BMJ + BMK ) / H
               AXIAL = P(J)
      ROTN = DW * 57.29577951
      DEF = W(J)
          IF ( J .EQ. 3 .OR. J .EQ. MP5 ) GO TO 8200
               DEFLCT(J-3) = W(J)
               ENMONT(J-3) = BMJ / 1000.0
               XXX(J-3)    = X / 12.0
      ROTATE(J-3)=ROTN
      LAST = J-3
 8200 CONTINUE
      SUM1 = 0.0
      SUM2 = 0.0
      DO 8350 J = 4, MP4
      SUM1 = SUM1 + QF(J)
      SUM2 = SUM2 + Q(J) - S(J) * W(J)
 8350 CONTINUE
      ENMONT(LAST)=ENMONT(LAST)*2.
      ENMONT(1)=ENMONT(1)*2.
      XMINM=ENMONT(1)
      XMAXM=ENMONT(1)
      LOCAT1=1
      LOCAT2=1
      DO 9012 J25=2,LAST
      IF(ENMONT(J25).LE.XMAXM) GO TO 9011
      XMAXM=ENMONT(J25)
      LOCAT1=J25
9011  IF(ENMONT(J25).GE.XMINM) GO TO 9012
      XMINM=ENMONT(J25)
      LOCAT2=J25
9012  CONTINUE

      IF (ILOADS.GE.2) THEN
          IF (LDLOCK.EQ.0) THEN
             DO 9014 J = 1, MP7
               QFSAVE(J) = QF(J)
               TSAVE(J)=T(J)
9014           PSAVE(J) = P(J)
             LDLOCK = 1
           ENDIF
          IF (KD.EQ.1) KD = 2
          WRITE(9,9130) LDFAC*100., DEFLCT(1),ROTATE(1),
     1           XMAXM, XXX(LOCAT1), SUM1/1000.,SUM2, ERR
          IF (ERR.NE.' ') THEN
             LDFAC2=LDFAC
             LDFAC=(LDFAC1+LDFAC2)/2.
          ELSE
             LDFAC1=LDFAC
             IF (LDFAC.GE.LDFAC2) THEN
                LDFAC=LDFAC*2.
                LDFAC2=LDFAC
             ELSE
                LDFAC=(LDFAC2+LDFAC1)/2.
             ENDIF
          ENDIF

          IF (ABS((LDFAC2-LDFAC1)/LDFAC2).LE.0.01) THEN
             WRITE (9,9133) LDFAC1*100.,LDFAC1*MT,LDFAC1*HT,LDFAC1*VT
9133         FORMAT (/,T20,'ULTIMATE FOUNDATION CAPACITY',F10.3,' %',
     1 /,'    MOMENT = ',F12.3,' FT-K   SHEAR = ', F12.3,' AXIAL = ',
     2 F12.3)
             RETURN
          ENDIF

          DO 9015 J = 1, MP7
             QF(J)=QFSAVE(J)*LDFAC
             T(J)=TSAVE(J)*LDFAC
9015         P(J)=PSAVE(J)*LDFAC

          GO TO 2000

       ELSE

      IF (IRANGE.EQ.0) WRITE(9,9010)DEFLCT(1),DEFLCT(LAST),ROTATE(1),
     1ROTATE(LAST),ENMONT(1),ENMONT(LAST),XMAXM,
     2XXX(LOCAT1),XMINM,XXX(LOCAT2),SUM1,SUM2

      IF (IRANGE.NE.0) WRITE(9,9120) M*INCR/12., DEFLCT(1),ROTATE(1),
     1XMAXM, XXX(LOCAT1), SUM1/1000.,SUM2, ERR

      ENDIF

9120  FORMAT (F10.1,2F10.3,F10.1,F7.3,F9.3,F9.5,1X,A)
9130  FORMAT (F10.2,2F10.3,F10.1,F7.3,F9.3,F9.5,1X,A)

C-----RETURN TO TAKE NEXT PROBLEM
          GO TO 1010
C-----COMPLETION OF RUN
 9999 CONTINUE
      RETURN
      END

C    $SUBTITLE: 'SUBROUTINE OUT '


      SUBROUTINE OUT
      COMMON /BK4/ IEND
      COMMON /RANGE/ IRANGE,INCR,THICKNESS,E,EI,IPIPE
      COMMON /LOADS/ ILOADS,KD
      COMMON /BK1/ TITLE
      COMMON /BK3/ MT, HT, VT,IPROB, DEPTH, DIAM, CAP, Z(20), C(20),
     1 DENSE(20), E50(20), PHI(20), H2O, NUMST, IU8, TYPE,IDT
      INTEGER TYPE(20)
      REAL MT
      CHARACTER IDT*8,TITLE*60
500   FORMAT( 14X,'DRILLED PIER LATERAL LOAD ANALYSIS IBM ',
     1 'VERSION 2.000009',
     1 /, 6X,'DATE',1X, A10   ,37X,'        RUN NO.',I3,
     1 //,6X,'BY ',5('_'),12X,'STRUCTURE NO. ',16('_'),
     2 //,6X,'CK ',5('_'),16X,'TYPE',5X,16('_'),
     3 ///,6X,13X,A60,///,6X,74('_')/)
505   FORMAT(
     1 //,6X,'DEPTH',F7.2,'   FT.',T45,'MOMENT',   F15.3,1X,'FT-KIP',
     1 //,6X,'DIAM.',F7.2,'   FT.',T45,'SHEAR',  F15.3,2X,'KIPS',
     2 //,6X,   'CAP  ',F7.2,'   FT.',T45,'AXIAL',F15.3,2X,
     2'KIPS (-C,+T)',
     2 //,1P,6X,'E ',G11.4,'  PSI',T45,'EI   ',G15.4,2X,'LB-IN^2  ',
     1 //,26X,0P,'SOIL PARAMETERS',
     1 /,26X,4('='),1X,10('='),
     1 //,22X,'WATER TABLE.....  ',F7.2,2X,'FT.',
     1 //, 7X, 'DEPTH',4X,'C OR PHI',5X,'DENSITY',11X,'E', 8X,'TYPE',
     1 /,8X,'(FT)',5X,'(TSF)', 9X,'(PCF)',/,
     1   4X,60('*'),/)
775   FORMAT (3F12.2,F12.3,I12)
810   FORMAT (4X,60('*'),///,
     1 '     DEPTH     DEFL.  ROTATION   MAX MOM     AT'
     1 '     SUM1     SUM2 ERROR STATUS',/,4X,76('*'))
820   FORMAT (4X,60('*'),///, T35,'ULTIMATE CAPACITY',//,
     1 '    % OF LOAD  DEFL.  ROTATION   MAX MOM     AT'
     1 '     SUM1     SUM2 ERROR STATUS',/,4X,76('*'))

      IF (ILOADS.EQ.3) RETURN
      IF (IRANGE.EQ.3) RETURN
         WRITE (9,500) IDT,IPROB,TITLE

         IF (IPIPE.NE.1) THEN
          WRITE (9,506) '                      DRILLED PIER FOUNDATION'
506         FORMAT (A)
         ELSE
          WRITE (9,507) '         STEEL PIPE PILE FOUNDATION',THICKNESS
507         FORMAT (A,' WALL THICKNESS = ',F8.5,' IN.')
         ENDIF

         WRITE (9,505) DEPTH,MT, DIAM,HT,CAP,VT,E,EI,H2O
         L1=NUMST+1
         DO 503 L2 = 2,L1

           IF (PHI(L2).GE.10.) THEN
             WRITE (9,775) Z(L2), PHI(L2), DENSE(L2),E50(L2),TYPE(L2)
           ELSE
             WRITE (9,775) Z(L2), C(L2), DENSE(L2), E50(L2),TYPE(L2)
           ENDIF

503      CONTINUE

         IF (ILOADS.EQ.2) THEN
             WRITE (9,820)
             ILOADS = 3
          ENDIF

         IF (IRANGE.EQ.2) THEN
             WRITE (9,810)
             IRANGE = 3
          ENDIF

999   RETURN
      END
