C C2DP      SOURCE    CHAT      05/01/12    21:44:59     5004
      SUBROUTINE C2DP(SIG,DSIG,YOUN,ANU,RT1,RT2,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),TENS(3)
      ITER=0
      Y=YOUN/(1.D0-ANU*ANU)
      SIREF=1.E-6*YOUN
      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,'C2DP COUPLAGE COMP 2 DRUCKER',/)
 9998 FORMAT(1X,'C2DP ',I4,'ITERATIONS INTERNES',/)
C
C------------------------------------------------
C     COUPLAGE COMPRESSION 2 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
      HDP=0.D0
      ADP=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)=0.D0
      EC(2)=1.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(2)+X*DSFG(2))/SG(2)
      F=SSI(2)/SG(2)
      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
      DL2=-E-F*DL3
      DO 145 ITYP=1,3
  145 SGG(ITYP)=SFG(ITYP)+X*DSFG(ITYP)-DL2*SG(ITYP)-DL3*SSI(ITYP)
C
C----------------------------------------------
C     SI UN DES DLAMDA < 0  FAUX COUPLAGE
C----------------------------------------------
C
      IF(DL2.GT.0.D0) IDAM(2)=0
      IF(DL3.LT.0.D0) IDAM(3)=0
      IF(IDAM(2).EQ.0.OR.IDAM(3).EQ.0) THEN
      DL2=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((DL2+XLAM2).GE.0.D0) GOTO 200
C
C------------------------------------------------
C     CAS OU LA FISSURE EST COMPLETEMENT FERMEE
C------------------------------------------------
C
      IDAM(2)=0
      DL2=-XLAM2
      IF(DSFG(2).EQ.0.D0) GOTO 160
      E=(DL2*SG(2)-SFG(2))/DSFG(2)
      F=SSI(2)/DSFG(2)
      DO 150 ITYP=1,3
      A(ITYP)=SFG(ITYP)-DL2*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+F*DL31
      XIN2=E+F*DL32
      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 C2DP FERM ',/)
      DL3=DL31
      X=XIN1
      GOTO 200
  160 DL3=(SFG(2)-DL2*SG(2))/SSI(2)
      RD=RDP-HDP*DL3
      IF(RD.LT.0.D0) THEN
      RDP=0.D0
      GOTO 110
      ENDIF
      DO 170 ITYP=1,3
      A(ITYP)=SFG(ITYP)-DL2*SG(ITYP)-DL3*SSI(ITYP)
  170 B(ITYP)=DSFG(ITYP)
      CALL XDP(A,B,RD,ADP,X,ITEST)
  200 CONTINUE
      DO 205 ITYP=1,3
      DSGG(ITYP)=X*DSFG(ITYP)
  205 SGG(ITYP)=SFG(ITYP)+DSGG(ITYP)-DL2*SG(ITYP)-DL3*SSI(ITYP)
C
C-----------------------------------------------
C     ON REGARDE SI ON ENDOMMAGE PAS LES AUTRES
C     CRITERES
C-----------------------------------------------
C-------------------------------------
C     ON ENDOMMAGE LA TRACTION 1
C-------------------------------------
C
      CALL CTRAF(SGG(1),RT1,VCTR)
      IF(VCTR.LE.0.D0) GOTO 1000
      IDAM(1)=1
      IDAM(2)=-1
      IDAM(3)=1
      AT(1)=-DSFG(1)
      AT(4)=Y
      AT(2)=Y*ANU
      AT(3)=-DSFG(2)
      BT(1)=SFG(1)-RT1
      BT(2)=SFG(2)
      CALL SYLIN2(AT,BT,X1,DL21)
      BT(1)=-SSI(1)
      BT(2)=-SSI(2)
      CALL SYLIN2(AT,BT,X2,DL22)
      DO 230 ITYP=1,3
      A(ITYP)=SFG(ITYP)-DL21*SG(ITYP)+X1*DSFG(ITYP)
  230 B(ITYP)=-DL22*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=DL21+DL31*DL22
      DLIN2=DL21+DL32*DL22
      XIN1=X1+DL31*X2
      XIN2=X1+DL32*X2
      IF(XIN1.GT.-1.D-10.AND.DLIN1.LT.1.E-10.AND.XIN1.LE.X
     1 .AND.CO1.GT.0.9) THEN
      DL3=DL31
      X=XIN1
      DL2=DLIN1
      GOTO 3000
      ENDIF
      IF(XIN2.GT.-1.D-10.AND.DLIN2.LT.1.E-10.AND.XIN2.LE.X
     1 .AND.CO2.GT.0.9) THEN
      DL3=DL32
      X=XIN2
      DL2=DLIN2
      GOTO 3000
      ENDIF
      IF(IIMPI.EQ.9) WRITE(IOIMP,20202)
20202 FORMAT(1X,'ERREUR DANS C2DP TRAC ',/)
      DL3=DL31
      DL2=DLIN1
      X=XIN1
      GOTO 3000
C
C-------------------------------------
C     ON ENDOMMAGE LA COMPRESSION 1
C-------------------------------------
C
 1000 CALL CCOAF(SGG(1),XLAM1,VCCO1)
      IF(VCCO1.LE.0.D0) GOTO 2000
      IDAM(1)=-1
      IDAM(3)=1
      IDAM(2)=-1
      AT(1)=-DSFG(1)
      AT(4)=Y
      AT(2)=Y*ANU
      AT(3)=-DSFG(2)
      BT(1)=SFG(1)
      BT(2)=SFG(2)
      CALL SYLIN2(AT,BT,X1,DL21)
      BT(1)=-SSI(1)
      BT(2)=-SSI(2)
      CALL SYLIN2(AT,BT,X2,DL22)
      DO 1010 ITYP=1,3
      A(ITYP)=SFG(ITYP)-DL21*SG(ITYP)+X1*DSFG(ITYP)
 1010 B(ITYP)=-DL22*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   CORRECTION PV ???
*     DLIN2=DL11+DL32*DL12   SUITE ERREUR APPOLO
      DLIN1=DL21+DL31*DL22
      DLIN2=DL21+DL32*DL22
      XIN1=X1+DL31*X2
      XIN2=X1+DL32*X2
      IF(XIN1.GT.-1.D-10.AND.DLIN1.LT.1.E-10.AND.XIN1.LE.X
     1 .AND.CO1.GT.0.9) THEN
      DL3=DL31
      X=XIN1
      DL2=DLIN1
      GOTO 3000
      ENDIF
      IF(XIN2.GT.-1.D-10.AND.DLIN2.LT.1.E-10.AND.XIN2.LE.X
     1 .AND.CO2.GT.0.9) THEN
      DL3=DL32
      X=XIN2
      DL2=DLIN2
      GOTO 3000
      ENDIF
      IF(IIMPI.EQ.9) WRITE(IOIMP,30303)
30303 FORMAT(1X,'ERREUR DANS C2DP COMP ',/)
      DL3=DL31
      DL2=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
      XT=XT+X
      RDP=RDP-DL3*HDP
      RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP)))
      IF(SIREF.GT.RMAX) THEN
      RDP=0.D0
      GOTO 110
      ENDIF
      DO 2010 ITYP=1,3
      SFG(ITYP)=SFG(ITYP)+X*DSFG(ITYP)-DL2*SG(ITYP)-DL3*SSI(ITYP)
 2010 EPDP(ITYP)=DL3*EPC(ITYP)+EPDP(ITYP)
      XLAM2=XLAM2+DL2
      XLAM3=XLAM3+DL3
      IF(XLAM2.LT.1.E-10) XLAM2=0.D0
      IF(XLAM2.LT.1.E-10) IDAM(2)=0
      C=1.D0-1.D-10
      IF(XT.LT.C) GOTO 2020
      DO 2015 ITYP=1,3
      IDAM(ITYP)=0
 2015 DSIG(ITYP)=0.D0
      CALL CHREP(SFG,SIG,-ANG)
      CALL CHREP(EPDP,EPPLDP,-ANG)
      IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER
      RETURN
 2020 IF(IDAM(2).EQ.-1) GOTO 100
      XLAM2=0.D0
      GAMDP=10.D0
      GAMTR1=10.D0
      GAMTR2=10.D0
      GAMCO1=10.D0
      DO 2030 ITYP=1,3
      DSFG(ITYP)=DSFG(ITYP)*(1.D0-XT)
 2030 TENS(ITYP)=SFG(ITYP)+DSFG(ITYP)
      CALL CTRAF(TENS(1),RT1,VCTR1)
      CALL CTRAF(TENS(2),RT2,VCTR2)
      CALL CCOAF(TENS(1),XLAM1,VCCO1)
      CALL CDP(TENS,ADP,RDP,VCDP)
      IF(VCTR1.GT.0.D0) CALL GAMTAF(SFG(1),DSFG(1),RT1,GAMTR1)
      IF(VCTR2.GT.0.D0) CALL GAMTAF(SFG(2),DSFG(2),RT2,GAMTR2)
      IF(VCCO1.GT.0.D0) CALL GAMCAF(SFG(1),DSFG(1),GAMCO1)
      IF(VCDP.GT.0.D0) CALL GDP(SFG,DSFG,RDP,ADP,GAMDP)
      IDAM(1)=0
      IDAM(2)=0
      IDAM(3)=0
      GAM=MIN(GAMTR1,GAMTR2,GAMCO1,GAMDP)
      IF(GAM.GE.1.D0) THEN
      DO 2040 ITYP=1,3
      SFG(ITYP)=SFG(ITYP)+DSFG(ITYP)
 2040 DSIG(ITYP)=0.D0
      CALL CHREP(SFG,SIG,-ANG)
      RETURN
      ENDIF
      IF(ABS(GAM-GAMTR1).LE.1.E-10) IDAM(1)=1
      IF(ABS(GAM-GAMTR2).LE.1.E-10) IDAM(2)=1
      IF(ABS(GAM-GAMCO1).LE.1.E-10) IDAM(1)=-1
      IF(ABS(GAM-GAMDP).LE.1.E-10) IDAM(3)=1
      DO 2050 ITYP=1,3
      SFG(ITYP)=SFG(ITYP)-GAM*DSFG(ITYP)
 2050 DSFG(ITYP)=(1.D0-GAM)*DSFG(ITYP)
      CALL CHREP(SFG,SIG,-ANG)
      CALL CHREP(DSFG,DSIG,-ANG)
      CALL CHREP(EPDP,EPPLDP,-ANG)
      RETURN
 3000 RDP=RDP-DL3*HDP
      RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP)))
      IF(SIREF.GT.RMAX) THEN
      RDP=0.D0
      GOTO 110
      ENDIF
      XT=XT+X
      DO 3010 ITYP=1,3
      SFG(ITYP)=SFG(ITYP)+X*DSFG(ITYP)-DL2*SG(ITYP)-DL3*SSI(ITYP)
      DSFG(ITYP)=(1.D0-XT)*DSFG(ITYP)
 3010 EPDP(ITYP)=DL3*EPC(ITYP)+EPDP(ITYP)
      CALL CHREP(SFG,SIG,-ANG)
      CALL CHREP(DSFG,DSIG,-ANG)
      CALL CHREP(EPDP,EPPLDP,-ANG)
      XLAM2=XLAM2+DL2
      XLAM3=XLAM3+DL3
      IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER
      RETURN
      END