Télécharger endo3d.eso

Retour à la liste

endo3d
C ENDO3D    SOURCE    FD218221  24/02/07    21:15:10     11834          
      subroutine endo3d(XMAT,NMAT,sig0,sigf,deps,nstrs,
     # VAR0,VARF,NVARI,nbelas3d,teta1,teta2,dt,ierr1,iso,mfr,ifou,
     # istep,epst0,epstf,NBRFLU3D,NBSUPP3D,
c     #NBRENF3D,NBPARR3D,
     # NMAT1,NMAT2,NVARFLU3D,NVARSUP3D,NBVSCAL,NBVTENS,NBVPTENS,
c     NVARENF3D,
     # NBNMAX3D,NBNB3D,IDIMB3D,XE3D,tetaref)    
C      # NBNLOC3D,INLVIA3D,NBVIA3D,NBNLISO,NBNLUNI,NBPARISO,
C      # NBPARUNI,NBVARNL3D,NBVARISO,NBVARUNI,METHODNL,DHI,DHU,VU,
C      # IVALMAT3D)
 
 
c     tables de dimension fixe pour resolution des sytemes lineaires 
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
 
c     *** declaration complementaires si renforts **********************
 
c     nombre maxi de renforts==NBRENF3D
c      integer NBRENF
c     nombre de parametres materiaux par renfort==NBPARR3D      
c      integer NBPARR
c     nombre de variables internes par renfort NVARENF==NVARENF3D
c      integer  NVARENF
c     dimension compatibles avec cfluendo et idvar4      
c      parameter (NBRENF=0,NBPARR=25,NVARENF=25)
c     si NBRENF >0 inclure le fichier de declare_rnfr3d.h      
 
c     *** declaration pour le non local ********************************
C      include './declare_non_local3d.h'
 
c     precision
      real*8 precision3d
      parameter(precision3d=1.0d-8)
c     nombre maxi de sous iterations
      integer itermax,iprint
      parameter(itermax=50000,iprint=itermax-100)
c     endommagement maxi
      real*8 dmaxi
      parameter (dmaxi=1.d0-precision3d)
c     methode de prise en compte de l endo post pic de cisaillement
      logical orthodp
      parameter (orthodp=.false.)      
 
c     variable logique pour matrice init isotrope
      logical iso
 
c     decalage Gauss pour la plasticite en cas de fluage  nf0
c     et decalage de la taille nmat1, ifou numero de la formulation    
      integer nf0,nmat1,nmat2,ifou,istep
 
c     *** nombre de parametres materiaux *******************************
      integer NBRFLU3D,NBSUPP3D,NBRENF3D
c     nombre de parametres materiau par renfort      
c      integer NBPARR3D
c     *** nombre de vari scalaire et pseudovecteur *********************
      integer NBVSCAL,NBVTENS,NBVPTENS         
c     compteur associees
      integer iscal,itens,icomp3d      
 
c     *** nombre de variables internes ********************************* 
c     nombre de variables internes
      integer NVARFLU3D,NVARSUP3D     
c     nombre de variables internes par renfort
c      integer NVARENF3D      
 
c     **** tableau des coordonnees des neouds de l element ************* 
      integer NBNMAX3D,NBNB3D,IDIMB3D
      real*8 xe3d(3,NBNMAX3D)
c     methode de calcul de la taille des elements
      logical Mjacob,Mnoeuds
      parameter (Mjacob=.false.,Mnoeuds=.true.)       
 
c     *** declaration des variables externes ***************************
      integer nmat,nstrs,NVARI,nbelas3d,ierr1,mfr
      real*8 xmat(nmat),sig0(nstrs),sigf(nstrs),deps(nstrs)
      integer IVALMAT3D(nmat)
      real*8 epst0(nstrs),epstf(nstrs)
      real*8 var0(NVARI),varf(NVARI)
      real*8 dt,teta1,teta2,tetaref      
 
c     *** tableau pour la resolution des systemes lineaires ************
c     taille maxi du systeme lineaire
      integer ngf
      parameter (ngf=65)
      real*8 XN(ngf),BN(ngf),ANN(ngf,(ngf+1))
      integer ipzero(ngf)
c     variable de controle d erreur dans methode de gauss
      integer errgauss
 
c     **** variables locales a plastendo3d *****************************
c     criteres
      integer ncr
      parameter (ncr=5)
      logical actif(ncr)
      real*8 fcr(ncr)
c     parametres materiaux donnes      
      real*8 young,nu,rt,ept,gft,ref,gfr,rc,epc,delta,beta,dcpp,dcpk
      real*8 ekdc,poro,pcc,mcc,ppcc,pfcc
c     coeff elastique
      real*8 MK,MG      
c     variables generales      
      integer i,j,k,l
c     indicateur premier passage
      logical ppas
c     coeff pour theta methode
      real*8 theta
      parameter (theta=0.5d0) 
C     cohesions de DP
      real*8 Cdp_min,Cdp_max,Rcmin
c     pression de consolidation de CC
      real*8 Ptmax,Pt,Pcmin,Pc
c     indicateur d endommagement prepic
      logical logdtpp,logdcpp
c     resistance effective et endommagement prepic
      real*8 dtppmax,dcppmax,rteff,rceff 
c     invariant de la deformation plastique de DP
      real*8 EPLDPV,EPLDPD 
c     invariant de la deformation plastique de CC
      real*8 EPLCCV,EPCCD       
c     deformations elastique et plastique de compression
      real*8 EEPC,EPLPIC,EDPPIC
c     module d ecrouissage
      real*8 Hcc,Hdp
c     vecteur de deformations mecaniques
      real*8 deps6(6),epst06(6),epstf6(6)
c     matrice deps
      real*8 deps33(3,3),deps3(3),vdeps33(3,3),vdeps33t(3,3)
c     contraintes effectives
      real*8 dsig3(3),dsig6(6),sigef33(3,3),sigef3(3),vsigef33(3,3)
c     contraintes bornées utilisees dans DP et CC      
      real*8 sigdp3(3)
c     tenseurs variables internes DP et CC debut de pas
      real*8 EPLDP6(6),EPLCC6(6)
c     fin de passage
      real*8 epldp6f(6),eplcc6f(6)
c     plasticite de traction
      real*8 eplr60(6),eplr6f(6),eplr33(3,3),eplr3(3)
c     gestion des ouvertures de fissures
      real*8 wplt06(6),wplt6(6),wpltx06(6)
      real*8 wpltx6(6),wpl3(3),vwpl33(3,3),vwpl33t(3,3)
      real*8 wplx3(3),vwplx33(3,3),vwplx33t(3,3)
c     contraintes finales
      real*8 sigef6(6),sigf6d(6),sigf6(6),sigt6(6),sigt6d(6)
c     pseudo tenseur auxiliaire
      real*8 xp6(6),xp33(3,3),vxp33(3,3),vxp33t(3,3),xp3(3)
      real*8 aux3(3) 
c     logique pour reduire l increment
      logical log_reduc,ecoul 
      real*8 reduc
      integer iter
c     endommagement
      real*8 umdt
c     contrainte hydrique
      real*8 sigp 
c     base ouverture des fissures
      real*8 depr3(3),vdepr33(3,3),vdepr33t(3,3),dir3(3),long3(3)
c     increment def plastique traction 
      real*8 deplr3(3),deplr6(6)     
c     increment def plastique dp
      real*8 depldp3(3),depldp6(6) 
c     increment def plastique cc
      real*8 deplcc3(3),deplcc6(6)
c     increment de def elastique pour retour radial
      real*8 depse3(3),depse6(6)
c     multiplicateurs plastiques
      real*8 lr1,lr2,lr3,ldp,lcc 
c     parametre de la loi de reduction thermique pour rc,rt,e 
c     ordre defini dans idvic et idvar4
      real*8 tt0(3),tt1(3),mtt(3),ptt(3),red(3)
      integer ntt
      parameter (ntt=4) 
c     prise en compte variation module due a theta
      real*8 dx,dMK,dMG 
c     compte sur la boucle de reduction d increment plastique
      integer jter
c     exposant endommagement prepic
      real*8 m1,dcpp0,dpic,spic,smin,s 
c     endommagements de traction-refermeture      
      real*8 dt3(3),ref3(3),wkt
c     calcul de la matrice d endo orthotrope      
      real*8 sn3(3),dt66(6,6)
c     partition des contraintes effectives
      real*8 sigf3(3),vsigf33(3,3),vsigf33t(3,3)
      real*8 sigft6(6),sigfc6(6),sigfc3(3),sigft3(3)
      real*8 sigft6p(6),sigft6d(6),sigft6dp(6)
c     refermeture des fissures
      real*8 dr66(6,6)
      real*8 sigfc6p(6),sigfc6dp(6),sigfc6d(6)
      real*8 sigf61(6),sigft61(6)
      real*8 sigf62(6),sigft62(6)
c     flambage des bielttes comprimees
      real*8 dct3(3)
c     coeff de couplage endo traction/flambage bielettes de compression
      real*8 altc,dctmax
c     endo ortho de DP
      real*8 dcc3(3),dcdp,dcom,pentendo,nendo  
c     indicateur de mise a zero du reducteur de retour plastique
      logical log_reduc_zero,log_err      
 
c     ***** declaration parametres materiaux des renforts repartis *****
c     include './declare_rnfr3d.h'
c     ! ne rien déclarer apres cet include qui contient des commandes !
 
c***********************************************************************
c     initialisation des parametres 
 
c     ** indicateur de premier passage 
      if (var0(1).ne.1.) then
         ppas=.true.
      else
         ppas=.false.
      end if 
 
c     initialisation de la variable de controle d erreur
      ierr1=0
 
 
c     *********** chargement des parametres materiaux depuis xmat ******
 
c     Module d young                            
      young=xmat(1)  
c     coefficient de Poisson      
      nu=xmat(2)
c     calcul des autres coeffs elastiques
      MK=young/(3.d0*(1.d0-2.0d0*nu))
      MG=young/(2.d0*(1.d0+nu))      
c     resistance a la traction      
      rt=xmat(nbelas3d+1)
c     deformation au pic de traction      
      ept=xmat(nbelas3d+2)
      if (ept.eq.0.) then
           ept=rt/young 
      end if 
c     energie de fissuration en traction directe
      gft=xmat(nbelas3d+3)
      if (gft.eq.0.) then
         print*,'Il manque l energie de fissuration GFT ds plasendo3d '
         ierr1=1
      end if       
c     seuil pour la refermeture des fissures
      ref=xmat(nbelas3d+4)
      if(ref.lt.(precision3d*rt)) then
         ref=rt*precision3d
      end if
c     energie de fissuration pour l endo de traction
      gfr=xmat(nbelas3d+5)  
      if(gfr.eq.0.) then
         print*,'Il manque l energie de refermeture GFR ds plastendo3d'
         ierr1=1
      end if 
c     resistance en compression-par cisaillement
      rc=xmat(nbelas3d+6)
c     deformation au pic de compression
      epc=xmat(nbelas3d+7)
c     endommagement au pic de compression
      dcpk=xmat(nbelas3d+8)      
c     coeff drucker prager
      delta=xmat(nbelas3d+9)         
c     coeff dilatance de cisaillement
      beta=xmat(nbelas3d+10)          
c     deformation caracteristique endo de compression
      ekdc=xmat(nbelas3d+11)
      if(ekdc.eq.0.) then
         print*,'il manque EKDC dans plastendo3d'
         ierr1=1
      end if
c     porosite initiale     
      poro=xmat(nbelas3d+12)       
c     Module d ecrouissage initial pour Cam Clay
      mcc=xmat(nbelas3d+13) 
 
c     recup des caracteristiques pour la reduc thermique
      do i=1,3
        tt0(i)=xmat(nbelas3d+14+(i-1)*ntt+1)      
        tt1(i)=xmat(nbelas3d+14+(i-1)*ntt+2)
        mtt(i)=xmat(nbelas3d+14+(i-1)*ntt+3) 
        ptt(i)=xmat(nbelas3d+14+(i-1)*ntt+4)
c       print*,tt0(i),tt1(i),mtt(i),ptt(i)          
      end do
c     coeff de couplage endo traction/endo compression     
      altc=xmat(nbelas3d+27)
c     pression de préconsolidation (initiale) pour camclay
      ppcc=xmat(nbelas3d+28) 
c     pour eviter une singulariute sur le gradient de CC iso avec pt=pc
c     il faut PPCC>RT
      if(ppcc.lt.rt*(1.d0+precision3d)) then
        print*,'Dans endo3d il faut PPCC>Rt'
        ierr1=1
      endif
c     pression de fin de consolidation pour CamClay 
      pfcc=xmat(nbelas3d+14)        
 
c***********************************************************************
c     initialisation des variables internes finales 
c***********************************************************************
 
      do iscal=1,NBVSCAL
         varf(iscal)=var0(iscal)
      end do
      do itens=1,NBVTENS
        do ICOMP3D=1,NBVPTENS
             nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D
             varf(nv3d)=var0(nv3d)
        end do
      end do
 
c***********************************************************************
c      prise en compte de la temperature 
c***********************************************************************
 
 
c     *** reduction des proprietes mecaniquess ************************* 
 
c     actualisation de la temperature maximale atteinte
      varf(15)=max(var0(15),teta2)
 
c     calcul de la reduction fin de pas  
      do i=1,3  
        call thed3d(teta2,tt0(i),tt1(i),mtt(i),ptt(i),red(i))
c       print*,'red(',i,')=',red(i)
c       le coeff de reductin thermique ne peut pas reaugmenter 
        if(.not.ppas) then       
           red(i)=min(red(i),var0(11+i))
        end if
        varf(11+i)=red(i)
      end do      
c     ppt meca effectees par la temperature
c     compression et parametres associes    
      rc=rc*red(1)
c     epc=epc/red(1)
c     traction et parametres associes      
      rt=rt*red(2)
      gft=gft*red(2)      
      ref=ref*red(2)
      gfr=gfr*red(2)
c     module d young et deformations associees      
      young=young*red(3) 
 
 
c     ** prise en compte de la variation du module *********************
 
c      variation du coeff de reduction du module
       dx=varf(14)-var0(14)
c      rajout de dE*epse sir le module varie       
       if((abs(dx).gt.precision3d).and.(.not.ppas)) then
c         le module d Young a varie
c         recup def elastique pas precedent
          ITENS=2        
          do ICOMP3D=1,6
            nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
            xp6(ICOMP3D)=var0(nv3d)
          end do
c         modification de la contrainte du a la variation du module
          dMK=dx*MK
          dMG=dx*MG         
c         diagonalisation du vecteur des deformations elastiques
          call x6x33(xp6,xp33)
          call b3_v33(xp33,xp3,vxp33)
c         increment de contrainte associees          
          call DSDE3D(dsig3,xp3,dMK,dMG)
c         passage en base fixe
          call CHRE36(dsig3,dsig6,vxp33)
c         mise a jour des contraintes
          ITENS=1 
          do ICOMP3D=1,6
             nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
             varf(nv3d)=varf(nv3d)+dsig6(ICOMP3D)
             sigef6(ICOMP3D)=varf(nv3d)
c            print*,'dans endo3d A sigef(',icomp3d,')=',sigef6(icomp3d)
          end do 
c          read*          
       end if    
 
c      actualisation des coeffs elastiques
       MK=young/(3.d0*(1.d0-2.0d0*nu))
       MG=young/(2.d0*(1.d0+nu))  
 
c     ***********test de coherence des donnees pour la matrice *********
 
c     ** existance d un endommagement pre pic en traction **************
 
      if (ept.gt.((rt/young)*(1.d0+precision3d)))then
c       il existe un endommagement pre pic en traction      
        logdtpp=.true.
c       endo pre pic maxi
        dtppmax=1.d0-rt/(young*ept)        
        rteff=rt/(1.d0-dtppmax)
      else
        logdtpp=.false.
        dtppmax=0.d0
        ept=rt/young
        rteff=rt
      end if
c     on continue avec rteff      
      rt=rteff
c     deformation elastique au pic de traction
      EEPT=rt/young 
c     pas de plasticite pre pic en traction
      EPPT=0.d0
c     deformation au pic de traction
      EPT=EEPT      
 
c     ** existance d un endommagement pre pic en compression ***********
 
c     ** verif validite de la dilatance
      if(beta.gt.(sqrt(3.d0)*(1.d0-precision3d))) then
        print*,'Dilatance beta excessive dans plastendo3d'
        ierr1=1
      end if
 
c     ** verif validite de l angle de frottement
      if(delta.gt.(sqrt(3.d0)*(1.d0-precision3d))) then
        print*,'Coeff de confinement excessif dans plastendo3d'
        ierr1=1
      end if  
 
      if(dcpk.gt.precision3d) then
c       on veut un endo de compression pre pic      
        if(epc.gt.((rc/young)*(1.d0+precision3d))) then
c         endo pre pic compression possible
          logdcpp=.true.        
          if(dcpk.le.(1.d0-rc/(young*epc)))then
c           l endommagement prepic de compression est conserve
            dcppmax=dcpk            
          else
c           l endommagement pre pic est trop grand , on le limite
            dcppmax=1.d0-rc/(young*epc) 
          end if
          rceff=rc/(1.d0-dcppmax)
        else
c         endo compression impossible
c          print*,'dcpp,Rc, et EPC sont incompatibles'
c          ierr1=1
          logdcpp=.false.
          dcppmax=0.d0
          rceff=rc        
        end if
      else
c       pas d endo de compression
        logdcpp=.false.
        dcppmax=0.d0
        rceff=rc
      end if 
c     deformation elastique au pic de compression
      EEPC=rceff/young
c     deformation plastique au pic de compression
      EPLPIC=EPC-EEPC 
c     deformation plastique equivalente au pic de compression
      EDPPIC= 0.3D1 * eplpic / (sqrt(0.3D1) - beta)
c     on continue avec rceff
      rc=rceff      
 
c     *** initialisation des parametres calculees a partir des donnees *
 
c     cohesion finale de drucker prager pour atteindre Rc
      Cdp_max=rc*(1.d0/sqrt(3.d0)-delta/3.d0)
c     cohesion initiale de Drucker prager pour activer Rankine en 
c     avant d'activer DP    
      Cdp_min = Rt * delta     
c     Cdp_minmanuel=0.8*rc*(1.d0/sqrt(3.d0)-delta/3.d0)
c     Cdp_min = min(max(Rt * delta, Rt*(1.d0/sqrt(3.d0)
c     # +(2.d0*delta/3.d0)),Cdp_minmanuel),0.9*Cdp_max*(1.0-dcppmax))   
      Cdp_min = min(max(Rt * delta, Rt*(1.d0/sqrt(3.d0)
     # +(2.d0*delta/3.d0))),0.9*Cdp_max*(1.0-dcppmax))  
c     coherence Rt, Rc pour que les Rankine et DP puissent etres actifs     
      if((Cdp_min*(1.d0+precision3d)).gt.Cdp_max) then
        print*,'Rt/Rc trop grand dans plastendo3d pour activer DP'
        print*,'il faut:  Rt/Rc< (sqrt(3)-delta)/(sqrt(3)+2*delta):'        
        ierr1=1
      end if 
 
c     *** initialisation des pressions de consolidation ****************
 
      if(ppcc.lt.precision3d*rt) then
        print*,'dans ENDO3D'
        print*,'La pression de preconsolidation ne peut pas etre nulle'
        print*,'il faut PPCC=',PPCC,'>> ',precision3d*rt
        ierr1=1
        return
      else if (pfcc.lt.ppcc*(1.d0+precision3d)) then
        print*,'dans ENDO3D'
        print*,'La pression de fin de consolidation ne peut pas '
        print*,'etre inferieure a la pression de preconsolidation'
        print*,'il faut PFCC=',PFCC,'>> ',PPCC
        ierr1=1
        return       
      end if
 
 
c     ********** parametres pour le non local **************************
 
C      include './initialise_non_local3d.h'      
 
c     *** bilan analyse des donnees ************************************
      if(ierr1.eq.1) then
         print*,'Incoherence dans les donnees de plastendo3d '
         return
      end if 
 
c***********************************************************************
c     chargement increment de deformation imposee et deformation finale
c***********************************************************************
      if(nstrs.lt.6) then
         do i=1,nstrs
             deps6(i)=deps(i)
             epstf6(i)=epstf(i)
             epst06(i)=epst0(i)             
         end do
         do i=nstrs+1,6
             deps6(i)=0.d0
             epstf6(i)=0.d0
             epst06(i)=0.d0
         end do
      else
         do i=1,6
             deps6(i)=deps(i)
             epstf6(i)=epstf(i)
             epst06(i)=epst0(i)
         end do      
      end if
c     passage en epsilon
      do i=4,6
         deps6(i)=0.5d0*deps6(i)
         epstf6(i)=0.5d0*epstf6(i)
         epst06(i)=0.5d0*epst06(i)         
      end do
c     mise a jour deformation elastique dans les varf
      ITENS=2
      do ICOMP3D=1,6
          nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
          varf(nv3d)=varf(nv3d)+deps6(ICOMP3D)
      end do    
 
 
c***********************************************************************
c      tir elastique
c***********************************************************************
 
c      ** directions principales de l increment de deformation *********
c      diagonalisation du vecteur des d eps 
       call x6x33(deps6,deps33)       
       call b3_v33(deps33,deps3,vdeps33)     
c      loi elastique 
       call DSDE3D(dsig3,deps3,MK,MG)
       if(dsig3(1).ne.dsig3(1)) then
          print*,'Pb2-1 lors du tir elastique dans endo3d'
          ierr1=1
          return
       end if          
c      passage en base fixe
       call CHRE36(dsig3,dsig6,vdeps33)
c      mise a jour des contraintes
       ITENS=1 
       do ICOMP3D=1,6
          nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
          varf(nv3d)=varf(nv3d)+dsig6(ICOMP3D)
          sigef6(ICOMP3D)=varf(nv3d)
c         print*,'dans endo3d A sigef(',icomp3d,')=',sigef6(icomp3d)
       end do               
c***********************************************************************
c      criteres de plasticité et ecoulements 
c***********************************************************************
 
 
c     ******************************************************************
c     recuperation deformation plastique debut et fin de de pas pour
c     controle de refermeture 
      itens=3 
      do ICOMP3D=1,6
         nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
         eplr60(ICOMP3D)=var0(nv3d)
         eplr6f(ICOMP3D)=eplr60(ICOMP3D)
      end do 
c     recuperation des deformation debut et fin de pas pour dp
      itens=4 
      do ICOMP3D=1,6
         nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
         epldp6(ICOMP3D)=var0(nv3d)
         epldp6f(ICOMP3D)=epldp6(ICOMP3D)
      end do       
 
c      **********debut des sous iterations locales *********************
       iter=0
c      mise a jour num iteration plastique       
10     iter=iter+1
c      reinitialisation du compteur sur la boucle de reduction de la
c      taille de l increment plastique
       jter=0
 
c     ******************************************************************
c     actualisation des parametres ecrouissables des criteres
c     ****************************************************************** 
 
c     ***** cohesion de DP (norme du deviateur) ************************
 
c     deformation plastique equivalente de DP debut de sous iteration     
      EPLDPD=varf(9)
c     actualisation de la cohesion de DP
      if(EDPPIC.gt.precision3d) then
c       il existe un ecrouissage pre pic en compression
        call ecrp3d(epldpd,edppic,Cdp_min,Cdp_max,Cdp,Hdp)
      else
c       cohesion max directement      
        Cdp=Cdp_max        
c       ecrouissage pre pic neglige            
        Hdp=0.d0
      endif     
c     stockage de la nouvelle cohesion
      varf(2)=Cdp
c     actualisation de la resistance a la compression pour DP
      Rc=Cdp/(1.d0/sqrt(3.d0)-delta/3.d0)
c      if(ppas) then
c         print*,'EPLDPD',EPLDPD
c         print*,'rcmin',Cdp_min/(1.d0/sqrt(3.d0)-delta/3.d0)       
c         print*,'rcmax',Cdp_max/(1.d0/sqrt(3.d0)-delta/3.d0)
c         print*,'Hdp',hdp
c         read*
c      end if  
 
c     **** pressions limites pour camclay ******************************
 
c     prise en compte de la deformation volumique de CamClay
      EPLCCV=0.d0
      ITENS=5
      do ICOMP3D=1,6
         nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D
         EPLCC6(ICOMP3D)=varf(nv3d)
         if(ICOMP3D.le.3) then
            EPLCCV=EPLCCV+EPLCC6(ICOMP3D)
         end if
      end do 
c      print*,'endo3d pbcc', eplccv,iter      
c     actualisation pc et pt  en fonction de de la variation de volume
      call PCCD3D(precision3d,poro,EPLCCV,ppcc,pfcc,pc,rt,pt)
      if(pc.ne.pc.or.pt.ne.pt) then
         print*,'Pb4 dans endo3: calcul de pc et pt de CamClay'
         print*,'precision3d,poro,EPLCCV,ppcc,pfcc,pc,rt,pt'
         print*,precision3d,poro,EPLCCV,ppcc,pfcc,pc,rt,pt
         ierr1=1
         return
      else
        varf(3)=pc
        varf(4)=pt
      end if
c     actualisation du module tangent de CamClay
c     (si EPLCCV mis a zero : on garde un ecrouissage constant) 
      call HCCD3D(pfcc,ppcc,poro,EPLCCV,Hcc)
c     **** tir elastique ***********************************************
 
c      diagonalisation du vecteur des contraintes effectives
       call x6x33(sigef6,sigef33)
       if(sigef33(1,1).ne.sigef33(1,1)) then
            print*,'Pb3 dans endo3d sigef6',sigef6
            print*,'iter',iter,' jter',jter
            do i=1,ncr        
                write(*,100) 'critere:',i,'activite:',
     #          actif(i),'fcr(',i,'):',fcr(i)   
            end do
            do itens=1,7
                do ICOMP3D=1,6
                  nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
                  print*,'Itenseur',itens,'v(',icomp3d,')=',varf(nv3d)
                end do  
            end do
            ierr1=1
            return
       end if 
       call b3_v33(sigef33,sigef3,vsigef33)
c      deformation plastique normale de Rankine dans la base principales
       call chrep6(eplr6f,vsigef33,.false.,xp6) 
       do i=1,3
         eplr3(i)=xp6(i)
       end do
c      ***** test activite des criteres ********************************
 
       call CREN3D(sigef3,precision3d,Rt,Rc,Ref,eplr3,
     # sigdp3,Cdp,delta,Mcc,pt,pc,ncr,actif,fcr,ppcc,pfcc,poro,eplccv)
       if(ierr1.eq.1) return
c      affichage des causes de sous iterations    
100    format (a9,i2,1x,a9,l1,1x,a4,i1,a2,e10.3) 
       if(iter.ge.iprint) then       
         print*,'Endo3d sous iteration locale:',iter
         do i=1,ncr
             write(*,100) 'critere:',i,'activite:',
     #       actif(i),'fcr(',i,'):',fcr(i)
         end do      
         if (iter.eq.(itermax+1)) then
             print*,'Nbr maxi de sous iterations locales'
             print*,'atteint dans endo3d: ', iter-1
             ierr1=1
             return
         end if
       end if     
c      test necessite de l ecoulement     
20     ecoul=.false.
       do i=1,ncr
            ecoul=ecoul.or.actif(i)
       end do
       if ((jter.ge.iprint).and.ecoul) then
            print*,'Criteres residuels pour iteration ',iter             
            do i=1,ncr
                if ((iter.ge.iprint).or.(jter.ge.iprint)) then
                    print*,'dans endo3d actif(',i,')=',actif(i)            
                    write(*,100) 'critere:',i,'activite:',
     #              actif(i),'fcr(',i,'):',fcr(i)               
                end if
            end do
       end if
c      print*,'dans endo3d ecoul:',ecoul
c      reinitialisation du log_reduc
       log_reduc=.false.
c      retour le cas echeant
       if(ecoul) then  
c           initialisation des increments
            do i=1,3
                deplr3(i)=0.d0
                depldp3(i)=0.d0
                deplcc3(i)=0.d0 
            end do                
c           choix de l ecoulement en fonction des criteres actifs
            if(actif(1).and.actif(2).and.actif(3)) then 
c               tri-traction-refermeture
                call EPLR3D(MG,MK,sigef3(1),sigef3(2),sigef3(3),
     #          fcr(1),fcr(2),fcr(3),deplr3,3)
c                print*,'endo3d: R1 R2 R3', deplr3(1),deplr3(2),deplr3(2)           
c               verif eplr3 final > 0
                call TSTD3D(fcr,eplr3,deplr3,actif,ncr,
     #          precision3d,log_reduc,reduc)
c               pas  cc ni dp
                lcc=0.d0
                ldp=0.d0
            else if (actif(1).and.actif(2).and.actif(4)) then
c               2 Rankine et et un DP
c               choix des contraintes limite pour DP
                if(fcr(1).gt.0.d0) then
                   R1=Rt
                else
                   R1=-Ref
                end if
                if(fcr(2).gt.0.d0) then
                   R2=Rt
                else
                   R2=-Ref
                end if                
c               calcul des multiplicateurs
                call LR12LDP3D(MK,MG,beta,delta,sigef3(1),sigef3(2),
     #          sigef3(3),Hdp,fcr(1),R1,fcr(2),R2,fcr(4),lr1,lr2,ldp,
     #          precision3d,log_err)   
                if(log_err) then
c                   on resoud par familles de criteres
                    if(0.5d0*(abs(fcr(1))+abs(fcr(2))).gt.fcr(4)) then
                        actif(1)=.true.
                        actif(2)=.true.
                        actif(4)=.false.
                    else
                        actif(1)=.false.
                        actif(2)=.false.                 
                        actif(4)=.true.
                    end if
c                   on recommence                 
                    goto 20                 
                end if                 
c               print*,'endo3d: R1 R2 et DP ',lr1, lr2 ,ldp      
c               calcul de l ecoulement de DP                   
                call EPLD3D(R1,R2,sigef3(3),beta,ldp,depldp3)
c               ecoulement R1 et R2                
                deplr3(1)=lr1
                deplr3(2)=lr2
                deplr3(3)=0.d0
                call TSTD3D(fcr,eplr3,deplr3,actif,ncr,
     #          precision3d,log_reduc,reduc) 
c               pas  cc
                lr3=0.d0
                lcc=0.d0
            else if (actif(2).and.actif(3).and.actif(4)) then
c               2 Rankine et  un DP
c               choix des contraintes limite pour DP
                if(fcr(2).gt.0.) then
                   R2=Rt
                else
                   R2=-Ref
                end if
                if(fcr(3).gt.0.) then
                   R3=Rt
                else
                   R3=-Ref
                end if                
c               calcul des multiplicateurs
                call LR12LDP3D(MK,MG,beta,delta,sigef3(2),sigef3(3),
     #          sigef3(1),Hdp,fcr(2),R2,fcr(3),R3,fcr(4),lr2,lr3,ldp,
     #          precision3d,log_err)
                if(log_err) then
c                   on resoud par familles de criteres
                    if(0.5d0*(abs(fcr(2))+abs(fcr(3))).GT.fcr(4)) then
                        actif(2)=.true.
                        actif(3)=.true.
                        actif(4)=.false.
                    else
                        actif(2)=.false.
                        actif(3)=.false.                 
                        actif(4)=.true.
                    end if
c                   on recommence                 
                    goto 20
                end if      
c                print*,'endo3d: R2 R3 et DP ',lr2, lr3 ,ldp     
c               calcul de l ecoulement de DP
                call EPLD3D(R2,R3,sigef3(1),beta,ldp,aux3)
                depldp3(1)=aux3(3)
                depldp3(2)=aux3(1)
                depldp3(3)=aux3(2)
c               ecoulement R1 et R2                
                deplr3(1)=0.d0
                deplr3(2)=lr2
                deplr3(3)=lr3                                
                call TSTD3D(fcr,eplr3,deplr3,actif,ncr,
     #          precision3d,log_reduc,reduc) 
c               pas  cc
                lcc=0.d0
                lr1=0.d0                   
            else if (actif(1).and.actif(3).and.actif(4)) then
c               2 Rankine et et un DP
c               choix des contraintes limite pour DP
                if(fcr(1).gt.0.) then
                   R1=Rt
                else
                   R1=-Ref
                end if
                if(fcr(3).gt.0.) then
                   R3=Rt
                else
                   R3=-Ref
                end if                
c               calcul des multiplicateurs
                call LR12LDP3D(MK,MG,beta,delta,sigef3(1),sigef3(3),
     #          sigef3(2),Hdp,fcr(1),R1,fcr(3),R3,fcr(4),lr1,lr3,ldp,
     #          precision3d,log_err)
                if(log_err) then
c                   on resoud par familles de criteres
                    if(0.5d0*(abs(fcr(1))+abs(fcr(3))).GT.fcr(4)) then
                        actif(1)=.true.
                        actif(3)=.true.
                        actif(4)=.false.
                    else
                        actif(1)=.false.
                        actif(3)=.false.                 
                        actif(4)=.true.
                    end if
c                   on recommance                 
                    goto 20                 
                end if     
c                print*,'endo3d: R1 R3 et DP ',lr1, lr3 ,ldp       
c               calcul de l ecoulement de DP
                call EPLD3D(R1,R3,sigef3(2),beta,ldp,aux3)
                depldp3(1)=aux3(1)
                depldp3(2)=aux3(3)
                depldp3(3)=aux3(2)
c               ecoulement R1 et R3                
                deplr3(1)=lr1
                deplr3(2)=0.d0
                deplr3(3)=lr3                
                call TSTD3D(fcr,eplr3,deplr3,actif,ncr,
     #          precision3d,log_reduc,reduc)
c               pas  cc
                lcc=0.d0
                lr2=0.d0                 
            else if (actif(1).and.actif(2)) then
c               bi-traction-refermeture       
                call EPLR3D(MG,MK,sigef3(1),sigef3(2),sigef3(3),
     #          fcr(1),fcr(2),fcr(3),deplr3,2)
c                print*,'endo3d: R1 et R2  ', deplr3(1),deplr3(2)      
c               verif eplr3 final > 0
                call TSTD3D(fcr,eplr3,deplr3,actif,ncr,
     #          precision3d,log_reduc,reduc)
c               pas  cc ni dp
                lcc=0.d0
                ldp=0.d0     
            else if (actif(1).and.actif(3)) then
c               bi-traction-refermeture       
                call EPLR3D(MG,MK,sigef3(1),sigef3(3),sigef3(2),
     #          fcr(1),fcr(3),fcr(2),aux3,2)
c                print*,'endo3d: R1 et R3  ', aux3(1),aux3(2) 
                deplr3(1)=aux3(1)     
                deplr3(2)=aux3(3)            
                deplr3(3)=aux3(2)
c               verif eplr3 final > 0
                call TSTD3D(fcr,eplr3,deplr3,actif,ncr,
     #          precision3d,log_reduc,reduc)
c               pas  cc ni dp
                lcc=0.d0
                ldp=0.d0     
            else if (actif(2).and.actif(3)) then
c               bi-traction-refermeture       
                call EPLR3D(MG,MK,sigef3(2),sigef3(3),sigef3(1),
     #          fcr(2),fcr(3),fcr(1),aux3,2)
c                print*,'endo3d: R2 et R3  ', aux3(1),aux3(2)      
                deplr3(1)=aux3(3)
                deplr3(2)=aux3(1)        
                deplr3(3)=aux3(2)
c               verif eplr3 final > 0
                call TSTD3D(fcr,eplr3,deplr3,actif,ncr,
     #          precision3d,log_reduc,reduc)
c               pas  cc ni dp
                lcc=0.d0
                ldp=0.d0    
            else if (actif(1).and.actif(4)) then 
c               1 rankine et dp actifs
c               choix de la contrainte limite pour DP
                if(fcr(1).gt.0.) then
                   R1=Rt
                else
                   R1=-Ref
                end if
c               multiplicateurs plastiques                
                call LR1LDP3D(MK,MG,beta,delta,sigef3(1),sigef3(2),
     #          sigef3(3),Hdp,fcr(1),R1,fcr(4),lr1,ldp,precision3d,
     #          log_err)
                if(log_err) then
                    if(abs(fcr(1)).gt.fcr(4)) then
                        actif(1)=.true.
                        actif(4)=.false.
                    else
                        actif(1)=.false.
                        actif(4)=.true.
                    end if
                    goto 20
                end if     
c                print*,'endo3d: R1 et DP ',lr1 ,ldp
c               calcul de l ecoulement de DP
                call EPLD3D(R1,sigef3(2),sigef3(3),beta,ldp,aux3)
                depldp3(1)=aux3(1)
                depldp3(2)=aux3(2)
                depldp3(3)=aux3(3)
c               ecoulement R1 et R2                
                deplr3(1)=lr1
                deplr3(2)=0.d0
                deplr3(3)=0.d0                 
                call TSTD3D(fcr,eplr3,deplr3,actif,ncr,
     #          precision3d,log_reduc,reduc) 
c               pas  cc
                lcc=0.d0     
            else if (actif(2).and.actif(4)) then
c               1 rankine et dp actifs
c               choix de la contrainte limite pour DP
                if(fcr(2).gt.0.) then
                   R2=Rt
                else
                   R2=-Ref
                end if
c               multiplicateurs plastiques                
                call LR1LDP3D(MK,MG,beta,delta,sigef3(2),sigef3(1),
     #          sigef3(3),Hdp,fcr(2),R2,fcr(4),lr2,ldp,precision3d,
     #          log_err)
                if(log_err) then
                    if(abs(fcr(2)).gt.fcr(4)) then
                        actif(2)=.true.
                        actif(4)=.false.
                    else
                        actif(2)=.false.
                        actif(4)=.true.
                    end if
                    goto 20
                end if
c               print*,'endo3d: R4 et DP ',lr2 ,ldp
c               calcul de l ecoulement de DP
                call EPLD3D(R2,sigef3(1),sigef3(3),beta,ldp,aux3)
                depldp3(1)=aux3(2)
                depldp3(2)=aux3(1)                
                depldp3(3)=aux3(3)
c               ecoulement R1 et R2                
                deplr3(1)=0.d0
                deplr3(2)=lr2
                deplr3(3)=0.d0                 
                call TSTD3D(fcr,eplr3,deplr3,actif,ncr,
     #          precision3d,log_reduc,reduc)
c               pas  cc
                lcc=0.d0    
            else if (actif(3).and.actif(4)) then
c               1 rankine et dp actifs
c               choix de la contrainte limite pour DP
                if(fcr(3).gt.0.) then
                   R3=Rt
                else
                   R3=-Ref
                end if
c               multiplicateurs plastiques                
                call LR1LDP3D(MK,MG,beta,delta,sigef3(3),sigef3(2),
     #          sigef3(1),Hdp,fcr(3),R3,fcr(4),lr3,ldp,precision3d,
     #          log_err)
                if(log_err) then
                    if(abs(fcr(3)).gt.fcr(4)) then
                        actif(3)=.true.
                        actif(4)=.false.
                    else
                        actif(3)=.false.
                        actif(4)=.true.
                    end if
                    goto 20
                end if
c               print*,'endo3d: R3 et DP ',lr3 ,ldp
c               calcul de l ecoulement de DP
                call EPLD3D(R3,sigef3(2),sigef3(1),beta,ldp,aux3)
                depldp3(1)=aux3(3) 
                depldp3(2)=aux3(2)                
                depldp3(3)=aux3(1)
c               ecoulement R1 et R2                
                deplr3(1)=0.d0
                deplr3(2)=0.d0
                deplr3(3)=lr3               
                call TSTD3D(fcr,eplr3,deplr3,actif,ncr,
     #          precision3d,log_reduc,reduc)
c               pas  cc
                lcc=0.d0      
            else if (actif(4).and.actif(5)) then
c               DP et CC actifs 
c               calcul de derivee de dcc=d(fcc) / d(lcc)
                call DCCD3D(pt,pc,Mcc,Hcc,sigef3(1),sigef3(2),
     #          sigef3(3),dcc)
c               calcul des multiplicateurs   
                call LCCL3D(Mcc,MK,MG,beta,delta,pt,pc,sigef3(1),
     #          sigef3(2),sigef3(3),dcc,Hdp,lcc,ldp,fcr(5),fcr(4),
     #          precision3d,log_err)
                if(log_err) then
                   if(fcr(5).gt.fcr(4)) then
                      actif(5)=.true.
                      actif(4)=.false.
                   else
                      actif(5)=.false.
                      actif(4)=.true.
                   end if
                   goto 20
                end if
c               ecoulement de DP
                call EPLD3D(sigef3(1),sigef3(2),sigef3(3),beta,
     #          ldp,depldp3) 
c               ecoulement de CC
                call EPLC3D(pt,pc,Mcc,sigef3(1),sigef3(2),sigef3(3),
     #          lcc,deplcc3)
                print*,'endo3d: DP et CC actifs ',ldp ,lcc
c               test de proprosité non negative                   
                call TSTC3D(poro,EPLCCV,deplcc3,precision3d,
     #          log_reduc,reduc)  
                if(reduc.lt.precision3d) then
                    actif(5)=.false.
                    fcr(5)=0.d0
                    lcc=0.d0
                    log_reduc=.true.
                    do i=1,3
                       deplcc3(i)=0.d0
                    end do
                end if      
            else if (actif(1)) then
                call EPLR3D(MG,MK,sigef3(1),sigef3(2),sigef3(3),
     #          fcr(1),fcr(2),fcr(3),aux3,1)
c                print*,'endo3d: R1 seul ', aux3(1)      
                deplr3(1)=aux3(1)
                deplr3(2)=aux3(2)        
                deplr3(3)=aux3(3)
c               verif eplr3 final > 0
                call TSTD3D(fcr,eplr3,deplr3,actif,ncr,
     #          precision3d,log_reduc,reduc)
c               pas  cc ni dp
                lcc=0.d0
                ldp=0.d0     
            else if (actif(2)) then
                call EPLR3D(MG,MK,sigef3(2),sigef3(1),sigef3(3),
     #          fcr(2),fcr(1),fcr(3),aux3,1)
c                print*,'endo3d: R2 seul ', aux3(1)  
                deplr3(1)=aux3(2)
                deplr3(2)=aux3(1)        
                deplr3(3)=aux3(3)
c               verif eplr3 final > 0
                call TSTD3D(fcr,eplr3,deplr3,actif,ncr,
     #          precision3d,log_reduc,reduc)
c               pas  cc ni dp
                lcc=0.d0
                ldp=0.d0
            else if (actif(3)) then
                call EPLR3D(MG,MK,sigef3(3),sigef3(1),sigef3(2),
     #          fcr(3),fcr(1),fcr(2),aux3,1)
c                print*,'endo3d: R3 seul ', aux3(1)  
                deplr3(1)=aux3(2)
                deplr3(2)=aux3(3)        
                deplr3(3)=aux3(1)
c               verif eplr3 final > 0
                call TSTD3D(fcr,eplr3,deplr3,actif,ncr,
     #          precision3d,log_reduc,reduc)
c               pas  cc ni dp
                lcc=0.d0
                ldp=0.d0                
            else if (actif(4)) then
c               Drucker Prager seul
                call LDPD3D(fcr(4),MK,MG,beta,delta,Hdp,ldp)
c                print*,'endo3d: DP seul ',ldp                
c               ecoulement de DP
                call EPLD3D(sigef3(1),sigef3(2),sigef3(3),beta,
     #          ldp,depldp3)
c               pas de rankine ni de cc
200             format(a8,i1,a2,e10.3)
                lcc=0.d0    
            else if (actif(5)) then
c               CamClay seul
c               calcul de derivee de dcc=d(fcc) / d(lcc)
                call DCCD3D(pt,pc,Mcc,Hcc,sigef3(1),sigef3(2),
     #          sigef3(3),dcc)
c               multiplicateur plastique     
                call LCCD3D(pt,pc,Mcc,sigef3(1),sigef3(2),sigef3(3),MK,
     #           MG,dcc,fcr(5),lcc,precision3d)     
c                print*,'endo3d: CC seul lcc=',lcc 
c                print*,pt,pc,Mcc,sigef3(1),sigef3(2),sigef3(3),MK,
c     #           MG,dcc,fcr(5),lcc
c               ecoulement cc                
                call EPLC3D(pt,pc,Mcc,sigef3(1),sigef3(2),sigef3(3),
     #          lcc,deplcc3)
c                do i=1,3
c                  print*,'deplcc3(',i,')=',deplcc3(i)
c                end do                
c                test de porosité non negative                   
                 call TSTC3D(poro,EPLCCV,deplcc3,precision3d,
     #           log_reduc,reduc)
                 if(reduc.lt.precision3d) then
                    actif(5)=.false.
                    fcr(5)=0.d0
                    lcc=0.d0
                    log_reduc=.true.
                    do i=1,3
                       deplcc3(i)=0.d0
                    end do
                 end if                  
            else
                print*,'configuration imprevue dans endo3d'
                do i=1,5
                   print*,'critere:',i,actif(i),' fcr(',i,')=',fcr(i)
                end do
                ierr1=1
                return                
            end if
c           read*
c           test de dissipation sur les criteres CC et DP
            if(reduc.lt.precision3d) then
                log_reduc_zero=.true.                
            else
                log_reduc_zero=.false.
            end if
            if(((ldp.lt.-precision3d).and.actif(4)) .or.
     #         ((lcc.lt.-precision3d).and.actif(5))) then
c                print*,' Endo3d multiplicateurs  negatif',ldp,lcc
c               on indique qu il ne faut pas garder ce retour
                log_reduc=.true.
c               on met le retour a zero
                reduc=0.d0 
c               indicateur associe                
                log_reduc_zero=.true.                               
c               on desactive les critere fautifs
                if ((ldp.lt.-precision3d).and.actif(4)) then   
c                  print*,'ldp',ldp
                   actif(4)=.false.                 
                else if ((lcc.lt.-precision3d).and.actif(5)) then                      
c                  print*,'lcc',lcc
                   actif(5)=.false.
                else
                   print*,'pb6 dans endo3d'
                   ierr1=1
                   return                   
                end if                                   
            end if                 
c           modification de l ecoulement en cas de reduction de la refermeture
            if(log_reduc) then 
c              print*,'Cas log_reduc: reduc=',reduc
              do i=1,3
c                 print*,'avant reduc deplr3(',i,')=',deplr3(i)
                 deplr3(i)=reduc*deplr3(i)
                 depldp3(i)=reduc*depldp3(i)                     
                 deplcc3(i)=reduc*deplcc3(i)                  
              end do
c             reduction des multiplicaterus possiblement utilises apres              
              lcc=lcc*reduc              
              ldp=ldp*reduc 
c             traitement suivant la reduction                         
c             on recommance l ecoulement sans les criteres desactives
c             dans testref3d
              jter=jter+1
              if(jter.ge.iprint) then
                    print*,'sous iteration Endo3d'
                    print*,'reduction a la refermeture:',jter  
              end if                
              if(jter.eq.itermax) then
                    print*,'Itermax Reduction a la refartmeture atteint'
                    print*,'Criteres residuels pour itertation ',iter
                    do i=1,ncr                
                        write(*,100) 'critere:',i,'activite:',
     #                  actif(i),'fcr(',i,'):',fcr(i)
                    end do                     
                    ierr1=1
                    return
              else                  
c                   nouvelle iteration 
c                    print*,'log_reduc_zero=',log_reduc_zero
                    if(log_reduc_zero) then 
c                     pas de plasticité, on recommence en eliminant les
c                     criteres non dissipatifs  sans modifier la contrainte 
                      goto 20
                    end if
              end if
            end if                
c           mise a jour des deformations plastiques de traction
            call CHRE36(deplr3,deplr6,vsigef33)
            itens=3 
            do ICOMP3D=1,6
                nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
                varf(nv3d)=varf(nv3d)+deplr6(ICOMP3D)
                eplr6f(ICOMP3D)=varf(nv3d)
            end do             
c           mise a jour des deformations plastiques de dp
            call CHRE36(depldp3,depldp6,vsigef33)
            itens=4 
            do ICOMP3D=1,6
                nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
                varf(nv3d)=varf(nv3d)+depldp6(ICOMP3D)
                epldp6f(ICOMP3D)=varf(nv3d)
            end do 
c           deformation equivalente de dp
            epldpd=epldpd+ldp
            varf(9)=epldpd           
c           mise a jour des deformations plastiques de cc
            call CHRE36(deplcc3,deplcc6,vsigef33)
            itens=5 
            do ICOMP3D=1,6
                nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
                varf(nv3d)=varf(nv3d)+deplcc6(ICOMP3D)
                eplcc6f(ICOMP3D)=varf(nv3d)
            end do             
c           mise a jour def elastique dans les varf
            log_err=.false.
            do i=1,3
                depse3(i)=-(deplr3(i)+depldp3(i)+deplcc3(i))
                if(depse3(i).ne.depse3(i)) then
                  log_err=.true.
                  print*,'Dans endo3d pb5-0'
                  print*,'deplr3(i),depldp3(i),deplcc3(i)'
                  print*, deplr3(i),depldp3(i),deplcc3(i)
                end if
            end do
            if (log_err) then
                print*,'Endo3d Pb5 a l iter:',iter
                do i=1,ncr
                    write(*,100) 'critere:',i,'activite:',
     #              actif(i),'fcr(',i,'):',fcr(i)
                end do
                ierr1=1
                return
            end if                     
            call chre36(depse3,depse6,vsigef33)
            ITENS=2
            do ICOMP3D=1,6
                nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
                varf(nv3d)=varf(nv3d)+depse6(ICOMP3D)
            end do            
c           mise a jour des contraintes
            call DSDE3D(dsig3,depse3,MK,MG) 
c           retour en base fixe
            call CHRE36(dsig3,dsig6,vsigef33)
c           mise a jour des contraintes
            ITENS=1 
            do ICOMP3D=1,6
                nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
                varf(nv3d)=varf(nv3d)+dsig6(ICOMP3D)
                sigef6(ICOMP3D)=varf(nv3d)
c               print*,'dans endo3d B sigef(',icomp3d,')=',sigef6(icomp3d)                
            end do            
c           on teste a nouveau les criteres       
            goto 10
       else
c           aucun des criteres n est actif
            continue
       end if
c      recuperation des contraintes effectives convergees
30     ITENS=1 
       do ICOMP3D=1,6
            nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
            sigef6(ICOMP3D)=varf(nv3d)
c           print*,'dans endo3d B sigef(',icomp3d,')=',sigef6(icomp3d)                
       end do              
c      fin des procedures elasto-plastiques                 
c***********************************************************************
 
 
c***********************************************************************
c       ouvertures de fissure debut de pas 
c***********************************************************************
 
        itens=6
        do ICOMP3D=1,6
          nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
          wplt06(ICOMP3D)=var0(nv3d)
        end do
c       ouvertures de fissure maxi debut de pas
        itens=7
        do ICOMP3D=1,6
          nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
          wpltx06(ICOMP3D)=var0(nv3d)
        end do 
c       base principale des ouvertures de fissures
        call bsfs3d(eplr6f,eplr60,depr3,vdepr33,vdepr33t)
c       longueurs de l element dans cette base
        do i=1,3
           do l=1,3
            dir3(l)=vdepr33(l,i)
           end do
           call LNUD3D(long3(i),xe3d,NBNMAX3D,NBNB3D,dir3)
        end do           
c       mise a jour des ouvertures de fissures
        call wfis3d(depr3,long3,vdepr33t,wplt06,wplt6,
     #  vwpl33,vwpl33t,wpl3,wpltx06,wpltx6,vwplx33,vwplx33t,wplx3)
        itens=6
        do ICOMP3D=1,6
          nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
          varf(nv3d)=wplt6(ICOMP3D)
        end do
c       ouvertures maxi de fissure fin de pas
        itens=7
        do ICOMP3D=1,6
          nv3d=NBVSCAL+(ITENS-1)*NBVPTENS+ICOMP3D      
          varf(nv3d)=wpltx6(ICOMP3D)
        end do      
c       ouverture maxi actuelle
        varf(5)=max(wpl3(1),wpl3(2),wpl3(3)) 
 
 
c***********************************************************************
c       endommagements
c***********************************************************************
 
c       **** partition du tenseur des contraintes **********************
 
        call prtt3d(sigef6,sigf3,vsigf33,vsigf33t,
     #  sigft6,sigfc6,sigfc3,sigft3) 
 
c       *** endommagement prepic de compression ************************
 
        if(logdcpp) then
c           endo pre pic de compression possible
            dcpp0=var0(8)      
            dpic=dcppmax
            spic=Cdp_max
            smin=Cdp_min        
c           contrainte equivalente de DP        
c            s=fcr(4)+Cdp 
            s=Cdp
            if((epldpd.lt.edppic).and.(edppic.ge.precision3d)) then
c              endommagement pre pic lineaire
c               dcpp=dpic * ((s - smin) / (spic - smin)) 
               if (s.gt.smin) then
                   pentendo=(Cdp_max*(1.d0-dpic)-Cdp_min)
     #             /sqrt(Cdp_max-Cdp_min)
                   if (Cdp.eq.0.) then
                       ierr1=1
                       print*,'il faut dcp_min>0 dans endo3d'
                   end if
                   dcpp=1.d0-(Cdp_min+pentendo*sqrt(Cdp-Cdp_min))/Cdp                 
c                   print*,dcpp
               else
                   dcpp=0.d0
               end if
            else
c              pas de plasticite pre pic, on utilise le critere DP en 
c              contrainte effective pour l endo pre pic
c              exposant de la loi d endommagement
               m1=(-0.1D1 + dpic) / spic / dpic * (-spic + smin)
               dcpp=dpic * ((s - smin) / (spic - smin)) ** m1
            end if
            dcpp=max(dcpp,dcpp0)
         else
            dcpp=0.d0
            dcpp0=0.d0            
         end if
         varf(8)=dcpp
 
c       *** endommagement prepic de traction  **************************
 
        if(logdtpp) then
c           endo pre pic de compression possible
            dtpp0=var0(7)      
            dpic=dtppmax
            spic=rteff
            smin=0.d0        
c           contrainte equivalente de DP        
            s=max(sigef3(1),sigef3(2),sigef3(3),0.d0)
c           exposant de la loi d endommagement
            m1=(-0.1D1 + dpic) / spic / dpic * (-spic + smin)
c           endommagement pre pic tension            
            dtpp=dpic * ((s - smin) / (spic - smin)) ** m1
            dtpp=max(dtpp,dtpp0)
         else
            dtpp=0.d0  
         end if
         varf(7)=dtpp  
 
c       *** endommagement complet de traction **************************
 
c       ouverture caracteristique pour l' endommagement de traction
        wkt=gft/rt
        dtra=0.d0
        do i=1,3
c           endo localise base principale de fissuration       
            dt3(i)=wplx3(i)*(wplx3(i)+2.d0*wkt)/(wkt+wplx3(i))**2
c           borne de dt a dmaxi
            dt3(i)=min(dmaxi,dt3(i))
c           endommagement global de traction pour affichage             
            dtra=max(dtra,dt3(i))
c           indice de fissure
            sn3(i)=1.d0/(1.d0-dt3(i))            
        end do
c       matrice d endommagement de traction incluant l endo pre pic
        call b3_d66(nu,sn3,dt66,.false.,.false.)
c       passage  des contraintes effectives dans la base prin des endo
        call chrep6(sigft6,vwplx33,.false.,sigft6p)  
c       application du tenseur d endommagement 
        do i=1,6
          sigft6dp(i)=sigft6p(i)
          do j=1,6
             sigft6dp(i)=sigft6dp(i)-dt66(i,j)*sigft6p(j)
          end do          
        end do
c       retour des contraintes positives en base fixe
        call chrep6(sigft6dp,vwplx33t,.false.,sigft6d)
c       reconstruction du tenseur des contraintes endommage 
        do i=1,6
c           prise en compte de l endo iso prepic de traction        
            sigf61(i)=sigfc6(i)+sigft6d(i)*(1.d0-dtpp)
        end do
        varf(11)=1.d0-(1.d0-dtra)*(1.d0-dtpp)        
c       Nouvelle partition du tenseur des contraintes endommagés 
        call prtt3d(sigf61,sigf3,vsigf33,vsigf33t,
     #  sigft61,sigfc6,sigfc3,sigft3)        
 
c       *** fonction de refermeture ************************************
 
c       ouverture caracteristique pour l' endommagement de traction
        wkr=gfr/ref
        do i=1,3
c           endo localise base principale de fissuration       
            ref3(i)=1.d0-wpl3(i)*(wpl3(i)+2.d0*wkr)/(wkr+wpl3(i))**2
c           ajustement borne  mini            
            ref3(i)=max(ref3(i),0.d0)
c           ajustement borne maxi            
            ref3(i)=min(ref3(i),1.d0)
c           indice de refermeture
            sn3(i)=1.d0/ref3(i)            
        end do         
c       matrice des fonctions de refermeture
        call b3_d66(nu,sn3,dr66,.false.,.false.)
c       passage  des contraintes effectives dans la base prin des endo
        call chrep6(sigfc6,vwplx33,.false.,sigfc6p)  
c       application du tenseur de refermeture aux contraintes negatives
        do i=1,6
          sigfc6dp(i)=sigfc6p(i)
          do j=1,6
             sigfc6dp(i)=sigfc6dp(i)-dr66(i,j)*sigfc6p(j)
          end do          
        end do
c       retour des contraintes positives en base fixe
        call chrep6(sigfc6dp,vwplx33t,.false.,sigfc6d)
c       reconstruction du tenseur des contraintes avec refermetures
        do i=1,6
            sigf62(i)=sigfc6d(i)+sigft61(i)
        end do
c       Nouvelle partition du tenseur des contraintes endommagés 
        call prtt3d(sigf62,sigf3,vsigf33,vsigf33t,
     #  sigft62,sigfc6,sigfc3,sigft3)         
 
c       *** endommagement de flambage des bielettes comprimees *********
 
        call buck3d(altc,dt3,dct3,vwplx33,vwplx33t,sigfc6)
        dctmax=max(dct3(1),dct3(2),dct3(3))
 
 
c       *** endommagement de cisaillement ******************************
 
c       dilatance transverse au pic de l essai uniaxial
        if(orthodp) then
            dvtp = eplpic *(0.2D1*beta+sqrt(0.3D1))/(sqrt(0.3D1)-beta)
        else
            dvtp = 0.3D1 * beta * eplpic / (sqrt(0.3D1) - beta)
        end if
c       endo de cisaillement initial
        dcdp=1.d0-(1.d0-var0(10))/(1.d0-dcpp0)
c       actualisation endo de cisaillement        
        call damdp3d(precision3d,ekdc,dvtp,epldp6f,dcc3,sigfc6,
     #  sigft62,dcdp,orthodp)
        if(orthodp) then
            dccmax=max(dct3(1),dct3(2),dct3(3))       
        else
            dccmax=dcdp
        end if
c       prise en compte endommagement isotrope prepic de compression 
        do i=1,6
c             sigf6(i)=(sigft62(i)+sigfc6(i))*(1.d0-dcpp)
             sigf6(i)=sigft62(i)+sigfc6(i)*(1.d0-dcpp)
        end do
c       actualisation endo de cisaillement complet        
        dcom=1.d0-(1.d0-dccmax)*(1.d0-dcpp)
        varf(10)=dcom
 
 
c***********************************************************************        
c       traitement des erreurs         
        if(ierr1.eq.1) then
             print*,'pb lors du calcul des endommagements dans endo3d'
             return
        end if
 
 
c***********************************************************************        
c       affectation dans le tableau de sortie des contraintes
c       et prise  encompte des contribution des renforts
        do i=1,nstrs
             sigf(i)=sigf6(i)
c            print*,'sigf(',i,')=',sigf(i)
        end do 
 
c       indicateur de fin de 1er pas 
        if(ppas.and.((istep.eq.2).or.(istep.eq.0))) then     
             varf(1)=1.d0
        else
            if(ppas) then
                 varf(1)=0.d0
            else
                 varf(1)=var0(1)
            end if
        end if
 
c ******traitement erreur residuelle************************************
        if(ierr1.eq.1) then        
             return 
        end if 
 
        return
 
       end
c***********************************************************************