C DDNSTH    SOURCE    CHAT      05/01/12    22:38:08     5004
C DDNSTH    SOURCE    AM        95/12/27    21:26:13     1962
      SUBROUTINE DDNSTH(IR,IC,DP,PN,U,E,E0,EYP,EYN,PCP,PCN,PYP,PYN,
     *                  UYP,UYN,CP,CN,PSP,PSN,BP,BN,UEP,PEP,UEN,
     *                  PEN,UCP,UCN,E1P,E1N,E2P,E2N,USP,USN,UTP,
     *                  PTP,UTN,PTN,EP,EN)
C=======================================================================
C     3-PARAMETER MODEL
C
C     THIS PROGRAM TRACES THE HYSTERETIC BEHAVIOR OF R/C COMPONENTS
C     - "DDNSTH" REPRESENTS DEGRADING, DETERIORATING,
C        NONSYMMETRIC, SLIP-TYPE & TRILINEAR HYSTERETIC MODEL.
C
C    *ESSENTIAL VARIABLES
C          IR= NUMBER OF BRANCH OF HYSTERETIC RULE
C          IC= NON-ZERO VALUE INDICATES CHANGE OF BRANCH
C          DP= INCREMENTAL STRESS
C          PN= NEW STRESS
C          U= NEW DEFORMATION
C          E= INSTANTANIOUS STIFFNESS
C    *PERMANENT VARIABLES FOR ENVELOPE CURVES (LAST LETTER "P" INDICATES
C     POSITIVE SIDE, AND "N" INDICATES NEGATIVE SIDE)
C          E0= INITIAL STIFFNESS
C          EYP,EYN= POST-YIELDING STIFFNESS
C          PCP,PCN= CRACKING STRESS
C          PYP,PYN= YIELDING STRESS
C          UYP,UYN= YIELDING DEFORMATION
C
C    *PERMANENT VARIABLES FOR DEGRADATION
C          CP,CN= COEFFICIENT FOR COMMON POINTS OF UNLOADING BRANCH
C
C    *PERMANENT VARIABLES FOR SLIPPAGE
C          PSP,PSN= SLIPPING YIELD STRESS (WHEN NO SLIP IS EXPECTED,
C                   SET PSP>PYP,PSN<PYN)
C
C    *PERMANENT VARIABLES FOR DETERIORATION
C          BP,BN= RATIO OF PLASTIC DEFORMATION INCREASE PER UNIT
C                 INCREASE OF ABSORBED ENERGY.
C
C    *TEMPORAL VARIABLES FOR DEGRADATION
C          UEP,PEP= EXPERIANCED MAXIMUM (POSITIVE) POINT
C          UEN,PEN= EXPERIANCED MINIMUM (NEGATIVE) POINT
C          UCP,UCN= UNLOADING (RELOADING) TURNING POINT ON X-AXIS
C          E1P,E1N= FIRST UNLOADING (RELOADING) STIFFNESS
C          E2P,E2N= SECOND RELOADING (UNLOADING)STIFFNESS
C
C    *TEMPORAL VARIABLES FOR SLIPPAGE
C          USP,USN= TARGET SLIP YIELD DEF. (CRACK CLOSING POINT)
C
C     *TEMPORAL VARIABLES FOR DETERIORATION
C          UPT,PTP= EXPECTED NEW MAXIMUN POINT
C          UTN,PTN= EXPECTED NEW MINIMUM POINT
C          EP     = INCREASE OF ABSORBED ENERGY TO CACULATE
C                   NEW MAXIMUM POINT
C          EN     = INCREASE OF ABSORBED ENERGY TO CACULATE
C                   NEW MINIMUM POINT
C
C    *NECESSARY INITIAL CONDITION
C          IR= 1 , PN= (LONG TERM LOADING) , PEN=PEP= 0
C
C===================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C
      IO=IR
      UO=U
      PO=PN
      PN=PO+DP
      GO TO(10,20,30,40,50,60,70,80,90,100,110,120,130,140,150),IR
C
C     -----
C     IR=1
C     -----
C
10    IF(PN.LT.PCP)GO TO 11
      IR=2
      PO=PCP
      UO=PCP/E0
      E=(PYP-PCP)/(UYP-UO)
      PN=PCP+(PN-PCP)*E/E0
      GO TO 22
11    IF(PN.LT.PCN)GO TO 12
      E=E0
      U=PN/E
      GO TO 200
12    IR=3
      PO=PCN
      UO=PCN/E0
      E=(PYN-PCN)/(UYN-UO)
      PN=PCN+(PN-PCN)*E/E0
      GO TO 32
C
C     -----
C     IR=2
C     -----
C
20    IF(DP.GE.0)GO TO 22
21    IR=6
      PEP=PO
      UEP=UO
      EP=0.
      E1P=(PO+CP*PYP)/(UO+CP*PYP/E0)
      UCP=U-PO/E1P
      USP=UCP
      PN=PO+DP*E1P/E
      E=E1P
      GO TO 63
22    IF(PN.GT.PYP)GO TO 23
      U=PCP/E0+(PN-PCP)/E
      EN=EN+0.5*(PN+PO)*(U-UO)
      GO TO 200
23    IR=4
      PN=PYP+(PN-PYP)*EYP/E
      U=UYP+(PN-PYP)/EYP
      EN=EN+0.5*(PYP+PO)*(UYP-UO)+0.5*(PN+PYP)*(U-UYP)
      E=EYP
      GO TO 200
C
C     -----
C     IR=3
C     -----
C
30    IF(DP.LE.0.)GO TO 32
31    IR=7
      PEN=PO
      UEN=UO
      EN=0.
      E1N=(PO+CN*PYN)/(UO+CN*PYN/E0)
      UCN=UO-PO/E1N
      USN=UCN
      PN=PO+DP*E1N/E
      E=E1N
      GO TO 73
32    IF(PN.LT.PYN)GO TO 33
      U=PCN/E0+(PN-PCN)/E
      EP=EP+0.5*(PN+PO)*(U-UO)
      GO TO 200
33    IR=5
      PN=PYN+(PN-PYN)*EYN/E
      U=UYN+(PN-PYN)/EYN
      EP=EP+0.5*(PYN+PO)*(UYN-UO)+0.5*(PN+PYN)*(U-UYN)
      E=EYN
      GO TO 200
C
C     -----
C     IR=4
C     -----
C
40    IF(DP.LT.0.)GO TO 21
41    U=UYP+(PN-PYP)/EYP
      E=EYP
      EN=EN+0.5*(PO+PN)*(U-UO)
      GO TO 200
C
C     -----
C     IR=5
C     -----
C
50    IF(DP.GT.0.)GO TO 31
51    U=UYN+(PN-PYN)/EYN
      E=EYN
      EP=EP+0.5*(PO+PN)*(U-UO)
      GO TO 200
C
C     -----
C     IR=6
C     -----
C
60    IF(PN.LE.PEP)GO TO 63
61    IF(PEP.GT.PYP)GO TO 62
      IR=2
      E2=(PYP-PCP)/(UYP-PCP/E0)
      PN=PEP+(PN-PEP)*E2/E
      EN=EN+0.5*(PEP+PO)*(UEP-UO)
      PO=PEP
      UO=UEP
      E=E2
      GO TO 22
62    IR=4
      PN=PEP+(PN-PEP)*EYP/E
      EN=EN+0.5*(PEP+PO)*(UEP-UO)
      PO=PEP
      UO=UEP
      GO TO 41
63    IF(PN.LT.0.)GO TO 64
      U=UCP+PN/E1P
      X=0.5*(PO+PN)*(U-UO)
      EP=EP+X
      EN=EN+X
      E=E1P
      GO TO 200
64    IR=8
      X=0.5*PO**2/E1P
      PO=0.
      UO=UCP
      EP=EP-X
      EN=EN-X
      IF(PEN.LE.PCN)GO TO 65
      E1N=E0
      E2N=-PCN/(UCP-PCN/E0)
      PN=PN*E2N/E1P
      E=E2N
      GO TO 84
CAVP2     65    IF(PSN.GT.PYN)GO TO 66
65    IF(PSN.GT.PYN.AND.PSN.GT.PEN)GO TO 66
      XE=0.5*PEN*(UEN-UCP)+EN
CAVP1
         IF(XE.LT.0.0)XE=0.0
CAVP1
      UTN=UEN-BN*XE
      E2=(PYN-PCN)/(UYN-PCN/E0)
      PTN=PCN+E2*(UTN-PCN/E0)
      IF(UTN.LT.UYN)PTN=PYN+EYN*(UTN-UYN)
      E1N=(CN*PYN+PTN)/(UTN+CN*PYN/E0)
      E2N=PTN/(UTN-UCP)
      PN=PN*E2N/E
      E=E2N
      GO TO 87
66    X=UEN-(PEN-PSN)/E1N-UCP
      PX=E0*(PSN*UCP/X)/(PSN/X-E0)
      XE=0.5*(UEN-UCP)*PX+0.5*PEN*(UEN-PX/E0)+EN
CAVP1
         IF(XE.LT.0.0)XE=0.0
CAVP1
      UTN=UEN-BN*XE
      E2=(PYN-PCN)/(UYN-PCN/E0)
      PTN=PCN+E2*(UTN-PCN/E0)
      IF(UTN.LT.UYN)PTN=PYN+EYN*(UTN-UYN)
      E1N=(CN*PYN+PTN)/(UTN+CN*PYN/E0)
      USN=UTN-PTN/E1N
      E2N=PSN/(UTN-(PTN-PSN)/E1N-UCP)
      PN=PN*E2N/E
      E=E2N
      GO TO 89
C
C     -----
C     IR=7
C     -----
C
70    IF(PN.GE.PEN)GO TO 73
71    IF(PEN.LT.PYN)GO TO 72
      IR=3
      E2=(PYN-PCN)/(UYN-PCN/E0)
      PN=PEN+(PN-PEN)*E2/E
      EP=EP+0.5*(PEN+PO)*(UEN-UO)
      PO=PEN
      UO=UEN
      E=E2
      GO TO 32
72    IR=5
      PN=PEN+(PN-PEN)*EYN/E
      EP=EP+0.5*(PEN+PO)*(UEN-UO)
      PO=PEN
      UO=UEN
      GO TO 51
73    IF(PN.GT.0.)GO TO 74
      U=UCN+PN/E1N
      X=0.5*(PO+PN)*(U-UO)
      EP=EP+X
      EN=EN+X
      E=E1N
      GO TO 200
74    IR=9
      X=0.5*PO**2/E1N
      PO=0.
      UO=UCN
      EP=EP-X
      EN=EN-X
      IF(PEP.GT.1.001*PCP)GO TO 75
      E1P=E0
      E2P=PCP/(PCP/E0-UCN)
      PN=PN*E2P/E1N
      E=E2P
      GO TO 94
CAVP1   75    IF(PSP.LT.PYP)GO TO 76
75    IF(PSP.LT.PYP.AND.PSP.LT.PEP)GO TO 76
      XE=0.5*PEP*(UEP-UCN)+EP
CAVP1
         IF(XE.LT.0.0)XE=0.0
CAVP1
      UTP=UEP+BP*XE
      E2=(PYP-PCP)/(UYP-PCP/E0)
      PTP=PCP+E2*(UTP-PCP/E0)
      IF(UTP.GT.UYP)PTP=PYP+EYP*(UTP-UYP)
      E1P=(CP*PYP+PTP)/(UTP+CP*PYP/E0)
      E2P=PTP/(UTP-UCN)
      PN=PN*E2P/E
      E=E2P
      GO TO 97
76    X=UEP-(PEP-PSP)/E1P-UCN
      PX=E0*(PSP*UCN/X)/(PSP/X-E0)
      XE=0.5*(UEP-UCN)*PX+0.5*PEP*(UEP-PX/E0)+EP
CAVP1
         IF(XE.LT.0.0)XE=0.0
CAVP1
      UTP=UEP+BP*XE
      E2=(PYP-PCP)/(UYP-PCP/E0)
      PTP=PCP+E2*(UTP-PCP/E0)
      IF(UTP.GT.UYP)PTP=PYP+EYP*(UTP-UYP)
      E1P=(CP*PYP+PTP)/(UTP+CP*PYP/E0)
      USP=UTP-PTP/E1P
      E2P=PSP/(UTP-(PTP-PSP)/E1P-UCN)
      PN=PN*E2P/E
      E=E2P
      GO TO 99
C
C     -----
C     IR=8
C     -----
C
80    IF(DP.LT.0.)GO TO 83
      IR=13
81    UCN=UO-PO/E1N
      PN=PO+DP*E1N/E2N
      E=E1N
      GO TO 73
83    IF(PEN.LT.1.001*PCN)GO TO 86
84    IF(PN.LT.PCN)GO TO 85
      E=E2N
      U=UCP+PN/E
      X=0.5*(PO+PN)*(U-UO)
      EP=EP+X
      GO TO 200
85    IR=3
      X=0.5*(PCN+PO)*(PCN/E0-UO)
      EP=EP+X
      E2=(PYN-PCN)/(UYN-PCN/E0)
      PN=PCN+E2*(PN-PCN)/E2N
      PO=PCN
      UO=PCN/E0
      E=E2
      GO TO 32
86    IF(PSN.GT.PYN)GO TO 89
87    IF(PN.GT.PTN)GO TO 88
      PEN=PTN
      UEN=UTN
      GO TO 71
88    E=E2N
      U=UO+(PN-PO)/E
      X=0.5*(PO+PN)*(U-UO)
      EP=EP+X
      EN=EN+X
      GO TO 200
CAVP89    PX=E0*E2N*UCP/(E2N-E0)
CPP189    PX=E2N*(USN-UCP)
CAVP89    PX=E2N*UCP*PTN/(E2N*UTN-PTN)
89    PX=E2N*(USN-UCP)
      IF(PN.GT.PX)GO TO 88
      IR=10
      X=0.5*(USN-UO)*(PX+PO)
      EP=EP+X
      EN=EN+X
      PO=PX
CAVP      UO=PX/E0
      UO=PX/E2N+UCP
      E2N=(PTN-PX)/(UTN-UO)
      PN=PX+(PN-PX)*E2N/E
      E=E2N
      GO TO 87
C
C     -----
C     IR=9
C     -----
C
90    IF(DP.GT.0.)GO TO 93
      IR=12
91    UCP=UO-PO/E1P
      PN=PO+DP*E1P/E2P
      E=E1P
      GO TO 63
93    IF(PEP.GT.1.001*PCP)GO TO 96
94    IF(PN.GT.PCP)GO TO 95
      E=E2P
      U=UCN+PN/E
      X=0.5*(PO+PN)*(U-UO)
      EN=EN+X
      GO TO 200
95    IR=2
      X=0.5*(PCP+PO)*(PCP/E0-UO)
      EN=EN+X
      E2=(PYP-PCP)/(UYP-PCP/E0)
      PN=PCP+E2*(PN-PCP)/E2P
      PO=PCP
      UO=PCP/E0
      E=E2
      GO TO 22
96    IF(PSP.LT.PYP)GO TO 99
97    IF(PN.LT.PTP)GO TO 98
      PEP=PTP
      UEP=UTP
      GO TO 61
98    E=E2P
      U=UO+(PN-PO)/E
      X=0.5*(PO+PN)*(U-UO)
      EP=EP+X
      EN=EN+X
      GO TO 200
CAVP99    PX=E0*EP2*UCN/(E2P-E0)
CAVP99    PX=E0*E2P*UCN/(E2P-E0)
CPP199    PX=E2P*(USP-UCN)
CAVP99    PX=E2P*UCN*PTP/(E2P*UTP-PTP)
99    PX=E2P*(USP-UCN)
      IF(PN.LT.PX)GO TO 98
      IR=11
      X=0.5*(USP-UO)*(PX+PO)
      EP=EP+X
      EN=EN+X
      PO=PX
CAVP      UO=PX/E0
      UO=PX/E2P+UCN
      E2P=(PTP-PX)/(UTP-UO)
      PN=PX+(PN-PX)*E2P/E
      E=E2P
      GO TO 97
C
C     -----
C     IR=10
C     -----
C
100   IF(DP.LT.0.)GO TO 87
      IR=15
      GO TO 81
C
C     -----
C     IR=11
C     -----
C
110   IF(DP.GT.0.)GO TO 97
      IR=14
      GO TO 91
C
C     -----
C     IR=12
C     -----
C
120   PX=(UCP-UCN)/(E1P-E2P)*E1P*E2P
      IF(PN.LT.PX)GO TO 63
      IR=9
      UX=UCP+PX/E1P
      X=0.5*(PX+PO)*(UX-UO)
      PN=PX+(PN-PX)*E2P/E1P
      E=E2P
      PO=PX
      UO=UX
      EP=EP+X
      EN=EN+X
      GO TO 93
C
C     -----
C     IR=13
C     -----
C
130   PX=(UCN-UCP)/(E1N-E2N)*E1N*E2N
      IF(PN.GT.PX)GO TO 73
      IR=8
      UX=UCN+PX/E1N
      X=0.5*(PX+PO)*(UX-UO)
      PN=PX+(PN-PX)*E2N/E1N
      E=E2N
      PO=PX
      UO=UX
      EP=EP+X
      EN=EN+X
      GO TO 83
C
C     -----
C     IR=14
C     -----
C
140   PX=(E1P*E2P*(UCP-UTP)+E1P*PTP)/(E1P-E2P)
      IF(PN.LT.PX)GO TO 63
      IR=11
      UX=UCP+PX/E1P
      X=0.5*(PX+PO)*(UX-UO)
      PN=PX+(PN-PX)*E2P/E1P
      E=E2P
      PO=PX
      UO=UX
      EP=EP+X
      EN=EN+X
      GO TO 97
C
C     -----
C     IR=15
C     -----
C
150   PX=(E1N*E2N*(UCN-UTN)+E1N*PTN)/(E1N-E2N)
      IF(PN.GT.PX)GO TO 73
      IR=10
      UX=UCN+PX/E1N
      X=0.5*(PX+PO)*(UX-UO)
      PN=PX+(PN-PX)*E2N/E1N
      E=E2N
      PO=PX
      UO=UX
      EP=EP+X
      EN=EN+X
      GO TO 87
C     ------------------------------------------------------------------
200   IC=0
      IF(IO.NE.IR)IC=IO
C=======================================================================
      RETURN
      END