C T1DP      SOURCE    CHAT      05/01/13    03:29:29     5004
      SUBROUTINE T1DP(SIG,DSIG,YOUN,ANU,RT1,RT2,H1,RDP,ADP,HDP,
     1 XLAM1,XLAM2,XLAM3,IDAM,ANG,EPPLDP,ALPHA,KERRE)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

-INC PPARAM
-INC CCOPTIO
      DIMENSION SIG(3),DSIG(3),SGG(3),DSGG(3),SFG(3),DSFG(3),
     1 A(3),B(3),IDAM(3),EPC(3),SSI(3),SI(3),EPC1(3),
     1 EPDP(3),EC(3),EPPLDP(3),SG(3),AT(4),BT(2)
      ITER=0
      SIREF=1.E-6*YOUN
      Y=YOUN/(1.D0-ANU*ANU)
      DO 10 ITYP=1,3
      EPDP(ITYP)=0.D0
   10 EPPLDP(ITYP)=0.D0
      DL3=0.D0
      XT=0.D0
      IF(IIMPI.EQ.9) WRITE(IOIMP,9999)
9999  FORMAT(1X,'T1DP COUPLAGE TRACTION 1 DRUCKER',/)
 9998 FORMAT(1X,'T1DP ',I4,'ITERATIONS INTERNES',/)
C
C------------------------------------------------
C     COUPLAGE TRACTION DRUCKER
C------------------------------------------------
C----------------------------------------------
C     DEFINITION D UNE CONTRAINTE DE REFERENCE
C     SI RMAX DU DRUCKER INFERIEURE A CETTE VALEUR
C     ALORS SIGMA=0
C----------------------------------------------
C
      RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP)))
      IF(SIREF.GT.RMAX) THEN
      DO 20 ITYP=1,3
  20  SIG(ITYP)=SIG(ITYP)+DSIG(ITYP)
      RDP=0.D0
      ADP=0.D0
      HDP=0.D0
      CALL CPHOIN(SIG,EPPLDP,YOUN,ANU)
      CALL NORME(SIG,DL3)
      DL3=DL3/YOUN*10.D0
      XLAM3=XLAM3+DL3
      DO 30 ITYP=1,3
      IDAM(ITYP)=0
      DSIG(ITYP)=0.D0
  30  SIG(ITYP)=0.D0
      IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER
      RETURN
      ENDIF
C-----------------------------------------------
C     ON SEPLACE DANS LE REPERE DE FISSURATION
C-----------------------------------------------
C
      CALL CHREP(SIG,SFG,ANG)
      CALL CHREP(DSIG,DSFG,ANG)
C
C-----------------------------------------------------
C     ESTIMATION DU PAS D'INCREMENT DE CONTRAINTES
C-----------------------------------------------------
C
      CALL SIJ(SFG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC)
      CALL CPHOMO(EPC,SSI,YOUN,ANU,ALPHA)
      CALL SCAL2(EPC,SSI,VAL)
      VAL=VAL-HDP
      IF(VAL.LT.0.D0) THEN
      KERRE=459
      RETURN
      ENDIF
      CALL SCAL(SSI,DSFG,VAL)
      CALL NORME(DSFG,VA1)
      CALL NORME(SSI,VA2)
      IF(VA1.EQ.0.D0) THEN
      IDAM(1)=0
      IDAM(2)=0
      IDAM(3)=0
      RETURN
      ENDIF
      RMIN=MIN((RDP/1.73),(RDP/(1.D0+2.D0*ADP)))
      X=VAL/VA1/VA2
      IF(X.GE.1.D0) THEN
      X=1.D0
      GOTO 41
      ENDIF
      X=1.D0/SQRT(1.0001D0-X*X)*RMIN/VA1/8.D0
      IF(X.GE.1.D0) X=1.D0
C
C----------------------------------------------------------
C     ON ECOULE ET ON REGARDE LA VARIATION DE NORMALE
C----------------------------------------------------------
C
C
C---------------------------------------
C     RESOLUTION DU SYSTEME EN DL3
C---------------------------------------
   41 X=X/2.D0
      EC(1)=1.D0
      EC(2)=0.D0
      EC(3)=0.D0
      CALL CPHOOB(EC,SG,YOUN,ANU)
  100 X=X*2.D0
      CALL SIJ(SFG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC)
      CALL CPHOMO(EPC,SSI,YOUN,ANU,ALPHA)
      CALL SCAL2(EPC,SSI,VAL)
      VAL=VAL-HDP
      IF(VAL.LT.0.D0) THEN
      KERRE=459
      RETURN
      ENDIF
      IF((XT+X).GT.1.D0) X=1.D0-XT
  110 RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP)))
C
C --------------------------------------------------
C     CAS OU LE RAYON MAXIMAL DU DRUCKER-PRAGER
C     INFERIEUR A LA CONTRAINTE DE REFERENCE
C --------------------------------------------------
C
      IF(SIREF.LT.RMAX) GOTO 131
      DO 120 ITYP=1,3
  120 SFG(ITYP)=SFG(ITYP)+(1.D0-XT)*DSFG(ITYP)
      RDP=0.D0
      ADP=0.D0
      HDP=0.D0
      CALL NORME(SFG,DL3)
      DL3=DL3/YOUN*10.D0
      XLAM3=XLAM3+DL3
      CALL CPHOIN(SFG,EPPLDP,YOUN,ANU)
      DO 130 ITYP=1,3
      IDAM(ITYP)=0
      EPDP(ITYP)=EPDP(ITYP)+EPPLDP(ITYP)
      DSIG(ITYP)=0.D0
  130 SIG(ITYP)=0.D0
      CALL CHREP(EPDP,EPPLDP,-ANG)
      IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER
      RETURN
  131 ITER=ITER+1
      IF(ITER.GT.200) THEN
      KERRE=460
      RETURN
      ENDIF
C
C---------------------------------
C     ON REALISE LE COUPLAGE
C---------------------------------
C
      E=-(SFG(1)+X*DSFG(1)-RT1)/(SG(1)-H1)
      F=SSI(1)/(SG(1)-H1)
      DO 140 ITYP=1,3
      DSGG(ITYP)=X*DSFG(ITYP)
      A(ITYP)=SFG(ITYP)+DSGG(ITYP)+E*SG(ITYP)
  140 B(ITYP)=-SSI(ITYP)+F*SG(ITYP)
      CALL DLAMDP(A,B,DL3,RDP,ADP,HDP,ITEST)
      IF(ITEST.EQ.1) THEN
      X=X/2.D0
      GOTO 131
      ENDIF
      DL1=-E-F*DL3
      DO 145 ITYP=1,3
  145 SGG(ITYP)=SFG(ITYP)+DSGG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP)
C
C----------------------------------------------
C     SI UN DES DLAMDA < 0  FAUX COUPLAGE
C----------------------------------------------
C
      IF(DL1.LT.-1.D-10) IDAM(1)=0
      IF(DL3.LT.-1.D-10) IDAM(3)=0
      IF(IDAM(1).EQ.0.OR.IDAM(3).EQ.0) THEN
      DL1=0.D0
      DL3=0.D0
      X=0.D0
      GOTO 3000
      ENDIF
C
C-----------------------------------------------------
C     ON VERIFIE QUE LA NORMALE SUIVANT LE DRUCKER
C     PRAGER TOURNE PEU
C-----------------------------------------------------
C
      CALL SIJ(SGG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC1)
      CALL SCAL(EPC,EPC1,VAL)
      CALL NORME(EPC,VA1)
      CALL NORME(EPC1,VA2)
      CO=VAL/VA1/VA2
      IF(CO.LT.0.99) THEN
      X=X/2.D0
      GOTO 131
      ENDIF
      IF((RT1-H1*DL1).GE.-1.E-10) GOTO 200
C
C------------------------------------------------
C     CAS OU LA LIMITE EN TRACTION CHUTE EN
C     DESSOUS DE 0
C------------------------------------------------
C
      DL1=RT1/H1
      IF(DSFG(1).EQ.0.D0) GOTO 160
      E=(DL1*SG(1)-SFG(1))/DSFG(1)
      F=SSI(1)/DSFG(1)
      DO 150 ITYP=1,3
      A(ITYP)=SFG(ITYP)-DL1*SG(ITYP)+E*DSFG(ITYP)
  150 B(ITYP)=F*DSFG(ITYP)-SSI(ITYP)
      CALL DLAMD(A,B,DL31,DL32,RDP,ADP,HDP,ITEST)
      DO 151 ITYP=1,3
      SGG(ITYP)=A(ITYP)+DL31*B(ITYP)
  151 DSGG(ITYP)=A(ITYP)+DL32*B(ITYP)
      CALL SIJ(SGG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC1)
      CALL SCAL(EPC,EPC1,VAL)
      CALL NORME(EPC,VA1)
      CALL NORME(EPC1,VA2)
      CO1=VAL/VA1/VA2
      CALL SIJ(DSGG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC1)
      CALL SCAL(EPC,EPC1,VAL)
      CALL NORME(EPC,VA1)
      CALL NORME(EPC1,VA2)
      CO2=VAL/VA1/VA2
      XIN1=E+DL31*F
      XIN2=E+DL32*F
      IF(XIN1.GT.-1.D-10.AND.XIN1.LE.X.AND.CO1.GT.0.9) THEN
      DL3=DL31
      X=XIN1
      GOTO 200
      ENDIF
      IF(XIN2.GT.-1.D-10.AND.XIN2.LE.X.AND.CO2.GT.0.9) THEN
      DL3=DL32
      X=XIN2
      GOTO 200
      ENDIF
      IF(IIMPI.EQ.9) WRITE(IOIMP,10101)
10101 FORMAT(1X,'ERREUR DANS T1DP DEPA ',/)
      DL3=DL31
      X=XIN1
      GOTO 200
  160 DL3=(SFG(1)-DL1*SG(1))/SSI(1)
      RD=RDP-DL3*HDP
      IF(RD.LT.0.D0) THEN
      RDP=0.D0
      GOTO 110
      ENDIF
      DO 170 ITYP=1,3
      A(ITYP)=SFG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP)
  170 B(ITYP)=DSFG(ITYP)
      CALL XDP(A,B,RD,ADP,X,ITEST)
C
C-----------------------------------------------
C     ON REGARDE SI ON ENDOMMAGE PAS LES AUTRES
C     CRITERES
C-----------------------------------------------
C-------------------------------------
C     ON ENDOMMAGE LA TRACTION 2
C-------------------------------------
C
  200 CONTINUE
      DO 205 ITYP=1,3
      DSGG(ITYP)=X*DSFG(ITYP)
  205 SGG(ITYP)=SFG(ITYP)+DSGG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP)
      CALL CTRAF(SGG(2),RT2,VCTR)
      IF(VCTR.LE.0.D0) GOTO 1000
      IDAM(1)=1
      IDAM(2)=1
      IDAM(3)=1
      AT(1)=-DSFG(1)
      AT(2)=Y-H1
      AT(4)=Y*ANU
      AT(3)=-DSFG(2)
      BT(1)=SFG(1)-RT1
      BT(2)=SFG(2)-RT2
      CALL SYLIN2(AT,BT,X1,DL11)
      BT(1)=-SSI(1)
      BT(2)=-SSI(2)
      CALL SYLIN2(AT,BT,X2,DL12)
      DO 230 ITYP=1,3
      A(ITYP)=SFG(ITYP)-DL11*SG(ITYP)+X1*DSFG(ITYP)
  230 B(ITYP)=-DL12*SG(ITYP)+X2*DSFG(ITYP)-SSI(ITYP)
      CALL DLAMD(A,B,DL31,DL32,RDP,ADP,HDP,ITEST)
      DO 231 ITYP=1,3
      SGG(ITYP)=A(ITYP)+DL31*B(ITYP)
  231 DSGG(ITYP)=A(ITYP)+DL32*B(ITYP)
      CALL SIJ(SGG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC1)
      CALL SCAL(EPC,EPC1,VAL)
      CALL NORME(EPC,VA1)
      CALL NORME(EPC1,VA2)
      CO1=VAL/VA1/VA2
      CALL SIJ(DSGG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC1)
      CALL SCAL(EPC,EPC1,VAL)
      CALL NORME(EPC,VA1)
      CALL NORME(EPC1,VA2)
      CO2=VAL/VA1/VA2
      DLIN1=DL11+DL31*DL12
      DLIN2=DL11+DL32*DL12
      XIN1=X1+DL31*X2
      XIN2=X1+DL32*X2
      IF(XIN1.GT.-1.D-10.AND.DLIN1.GT.-1.E-10.AND.XIN1.LE.X
     1 .AND.CO1.GT.0.9) THEN
      DL3=DL31
      X=XIN1
      DL1=DLIN1
      GOTO 3000
      ENDIF
      IF(XIN2.GT.-1.D-10.AND.DLIN2.GT.-1.E-10.AND.XIN2.LE.X
     1 .AND.CO2.GT.0.9) THEN
      DL3=DL32
      X=XIN2
      DL1=DLIN2
      GOTO 3000
      ENDIF
      IF(IIMPI.EQ.9) WRITE(IOIMP,20202)
20202 FORMAT(1X,'ERREUR DANS T1DP TRAC ',/)
      DL3=DL31
      DL1=DLIN1
      X=XIN1
      GOTO 3000
C
C-------------------------------------
C     ON ENDOMMAGE LA COMPRESSION 2
C-------------------------------------
C
 1000 CALL CCOAF(SGG(2),XLAM2,VCCO2)
      IF(VCCO2.LE.0.D0) GOTO 2000
      IDAM(1)=1
      IDAM(3)=1
      IDAM(2)=-1
      AT(1)=-DSFG(1)
      AT(2)=Y-H1
      AT(4)=Y*ANU
      AT(3)=-DSFG(2)
      BT(1)=SFG(1)-RT1
      BT(2)=SFG(2)
      CALL SYLIN2(AT,BT,X1,DL11)
      BT(1)=-SSI(1)
      BT(2)=-SSI(2)
      CALL SYLIN2(AT,BT,X2,DL12)
      DO 1010 ITYP=1,3
      A(ITYP)=SFG(ITYP)-DL11*SG(ITYP)+X1*DSFG(ITYP)
 1010 B(ITYP)=-DL12*SG(ITYP)+X2*DSFG(ITYP)-SSI(ITYP)
      CALL DLAMD(A,B,DL31,DL32,RDP,ADP,HDP,ITEST)
      DO 1011 ITYP=1,3
      SGG(ITYP)=A(ITYP)+DL31*B(ITYP)
 1011 DSGG(ITYP)=A(ITYP)+DL32*B(ITYP)
      CALL SIJ(SGG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC1)
      CALL SCAL(EPC,EPC1,VAL)
      CALL NORME(EPC,VA1)
      CALL NORME(EPC1,VA2)
      CO1=VAL/VA1/VA2
      CALL SIJ(DSGG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC1)
      CALL SCAL(EPC,EPC1,VAL)
      CALL NORME(EPC,VA1)
      CALL NORME(EPC1,VA2)
      CO2=VAL/VA1/VA2
      DLIN1=DL11+DL31*DL12
      DLIN2=DL11+DL32*DL12
      XIN1=X1+DL31*X2
      XIN2=X1+DL32*X2
      IF(XIN1.GT.-1.D-10.AND.DLIN1.GT.-1.E-10.AND.XIN1.LE.X
     1 .AND.CO1.GT.0.9) THEN
      DL3=DL31
      X=XIN1
      DL1=DLIN1
      GOTO 3000
      ENDIF
      IF(XIN2.GT.-1.D-10.AND.DLIN2.GT.-1.E-10.AND.XIN2.LE.X
     1 .AND.CO2.GT.0.9) THEN
      DL3=DL32
      X=XIN2
      DL1=DLIN2
      GOTO 3000
      ENDIF
      IF(IIMPI.EQ.9) WRITE(IOIMP,30303)
30303 FORMAT(1X,'ERREUR DANS T1DP COMP ',/)
      DL3=DL31
      DL1=DLIN1
      X=XIN1
      GOTO 3000
 2000 CONTINUE
C
C------------------------------------------------
C     CAS OU AUCUN AUTRE CRITERE EST ENDOMMAGE
C------------------------------------------------
C-----------------------------------
C     REMISE A JOUR DES VARIABLES
C     PUIS NOUVEAU PAS SI TOUT
C     L INCREMENT N EST PAS ECOULE
C-----------------------------------
C
      RDP=RDP-DL3*HDP
      RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP)))
      IF(RMAX.LT.SIREF) THEN
      RDP=0.D0
      GOTO 110
      ENDIF
      RT1=RT1-H1*DL1
      IF(RT1.LE.1.E-10) THEN
      H1=0.D0
      RT1=0.D0
      ENDIF
      XT=XT+X
      DO 2010 ITYP=1,3
      SFG(ITYP)=SFG(ITYP)+X*DSFG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP)
 2010 EPDP(ITYP)=DL3*EPC(ITYP)+EPDP(ITYP)
      XLAM1=XLAM1+DL1
      XLAM3=XLAM3+DL3
      C=1.D0-1.D-10
      IF(XT.LT.C) GOTO 100
      DO 2020 ITYP=1,3
      IDAM(ITYP)=0
 2020 DSIG(ITYP)=0.D0
      CALL CHREP(SFG,SIG,-ANG)
      CALL CHREP(EPDP,EPPLDP,-ANG)
      IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER
      RETURN
 3000 RDP=RDP-DL3*HDP
      RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP)))
      IF(RMAX.LT.SIREF) THEN
      RDP=0.D0
      GOTO 110
      ENDIF
      RT1=RT1-H1*DL1
      IF(RT1.LE.1.E-10) THEN
      H1=0.D0
      RT1=0.D0
      ENDIF
      XT=XT+X
      DO 3010 ITYP=1,3
      SFG(ITYP)=SFG(ITYP)+X*DSFG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP)
      DSFG(ITYP)=DSFG(ITYP)*(1.D0-XT)
 3010 EPDP(ITYP)=DL3*EPC(ITYP)+EPDP(ITYP)
      CALL CHREP(SFG,SIG,-ANG)
      CALL CHREP(DSFG,DSIG,-ANG)
      CALL CHREP(EPDP,EPPLDP,-ANG)
      XLAM1=XLAM1+DL1
      XLAM3=XLAM3+DL3
      IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER
      RETURN
      END