Télécharger fldo3d.eso

Retour à la liste

fldo3d
C FLDO3D    SOURCE    FD218221  24/02/07    21:15:11     11834          
       subroutine fldo3d(XMAT,NMAT,sig0,sigf,deps,nstrs,VAR0,VARF,
     # NVARI,teta1,teta2,dt,ierr1,iso,mfr,ifou,istep,epst0,epstf,
     # tetaref,NBELAS3D,NB_MAT_BASE,NB_MAT_SUPP,NB_VAR_BASE,NB_VAR_SUPP,
     # NB_DECAL,NBNMAX3D,NBNB3D,IDIMB3D,XE3D,NBVIA3D,INLVIA3D,
     # NMAT0,NMAT1,NMAT2,NMAT3,NVAR0,NVAR1,NVAR2,NVAR3,NB_PARA_FIBRE_U)
 
c     A Sellier oct 2022
c     NMAT0 fin des parametres elastiques
c     NMAT1 fin du modele de la matrice 
c     NMAT2 fin des parametres pour les modules complementaires de fluendo (exple rag Multon ou fibre Gontero)
c     NMAT3 fin des parametres pour les renforts repartis
c     NMAT  fin des parametres pour les champs de phase via Helmholtz
 
c     NVAR0 0
c     NVAR1 fin des variables internes pour la matrice
c     NVAR2 fin des variables internes pour les modules complementaires
c     NVAR3 fin des variables internes pour les renforts repartis
c     NVARI fin des variables internes pour les champs de phases via Helmholtz
 
c     calcul de l ecoulement visco-elasto-plastique: Sellier mars 2014
c     istep pour la phase non locale : juillet 2015
c     actualisation DTT : aout 2018
 
c     08/05/2020 : hypothese non locale 
c     istep=0 : calcul local
c     istep=1 : preparation du 1er passge non local
c     istep=2 : dernier passage non local effectue
c     istep=3 : passage non local intermediaire en non local iteratif
 
c     tables de dimension fixe pour resolution des sytemes lineaires 
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
 
c******** declaration des variables externes ***************************
      integer nmat,nstrs,NVARI,nbelas3d,ierr1,mfr
      real*8 xmat(nmat),sig0(nstrs),sigf(nstrs),deps(nstrs)
c      integer IVALMAT3D(nmat)
      real*8 epst0(nstrs),epstf(nstrs)
      real*8 var0(NVARI),varf(NVARI)
      real*8 dt,teta1,teta2,tetaref  
c     variable logique pour matrice init isotrope
      logical iso  
      integer ifou,istep
 
c     *** nombre de parametres materiaux *******************************
      integer NB_MAT_BASE,NB_MAT_SUPP
 
c     *** nombre de variables internes ********************************* 
c     nombre de variables internes
      integer NB_VAR_BASE,NB_VAR_SUPP
 
c     *** nombre de variables facultative entre les obligatoires et les helmholtz
      integer NB_DECAL  
 
c**************Coordonnes des noeuds de l element fini *****************        
c     ligne a inclure dans les modeles utilisant les noeuds     
      integer NBNMAX3D,NBNB3D,IDIMB3D
      real*8  XE3D(3,NBNMAX3D)
c***********************************************************************
 
 
c************* parametres pour les renforts ****************************
c     lignes a inclure dans les modeles utilisants les renforts
c      include './nombre_renforts.h'
-INC HNBRREN
c      include './declare_materiaux_renforts.h'
-INC HDECREN
c***********************************************************************
 
c************** parametres pour Helmholtz ******************************
c     lignes a inclure dans les modeles utilisant Helmholtz
c      include './nombre_helmholtz.h'
-INC HNBRHEL
c      include './declare_materiaux_helmholtz.h'
-INC HDECHEL
c***********************************************************************
 
 
c***********************************************************************
c     parametres obligatoires pour pointer correctement les renforts et 
c     les formulations d helmholtz
      integer NMAT0,NMAT1,NMAT2,NVAR0,NVAR1,NVAR2,NVAR3
c***********************************************************************
 
 
c     *** tableau pour la resolution des systemes lineaires ************      
      integer ngf
      parameter (ngf=65)
      real*8 XN(ngf),BN(ngf),ANN(ngf,(ngf+1))
      integer ipzero(ngf)
 
 
c     *** parametres materiaux *****************************************   
c     parametres materiaux donnes      
      real*8 hyd0,hyds
      real*8 young00,nu00,rt00,reft00,delta00,beta00,ept00   
      real*8 gft00,epsm00,psik,xflu,tauk00,taum00,taum0,nrjm,dt80,nrjw
      real*8 biotw00,xnsat00,poro2,vrag00,ekdg,gfr00,kshr 
      real*8 vw2,mvgn,epc0,ekdc
      real*8 tetar,tetas,dfmx
      real*8 taar,nrjg,srsrag,alatr,epssr
      real*8 mdtt00,wdtt,pdtt
      real*8 vref,vmax,cvrt00
      real*8 tdef,nrjp,srsdef,vdef00,cna,ttdd,ttrp,nrjd,tfid,tdid     
      real*8 exmd,exnd,cnab,cnak,ssad,ttdf,nrjf,nsul
      real*8 ekds,skdw,tkvg
      real*8 ttrg,tdyn 
      real*8 cshr,kdef,cdef,vvdef
      real*8 hplt00,hplw00,hplw,hplg00,hpls00 
      real*8 tetarw
 
c     variables pour la deformation thermique transitoire
      real*8 tdtt 
 
c     parametres materiaux deduits des donnes 
c     deformation elastique de reference pour le fluage    
      real*8 epser,epc1,epc00,eweb
 
c     **** variables locales a fldo3d *******************************
 
c     tenseur unitaires de reference
      real*8 vref33(3,3),vref33t(3,3) 
c     variables generales      
      integer i,j,k,l
      integer ncal,nmax,nmax1
c     nmax nombre maxi de sous iteration en retour radial   
c     nmax1 nombre maxi de sous iteration sans couplage des criteres diffus localises   
      parameter(nmax=200,nmax1=190)
 
      integer nc
c     nombre de criteres      
      parameter (nc=13)
      real*8 alphadp(3),betadp(3),factif(nc)
      logical actif(nc),complet(nc)
      real*8 vnorm
      logical ppas
      real*8 deps6(6),epstf6(6),epst06(6),epstf33(3,3),epst033(3,3)
      real*8 theta
      logical plast_seule,fl3d,end3d
      real*8 xwl2,tau_trac_wl2,xwl2max,cwrt,tau_trac_wl2_max
      real*8 sig13(3)
      real*8 young,nu,rt,rc,reft,rc1,delta,beta,gft,gfr
      real*8 ept,epsm,xnsat,biotw,Vvgel,kgel,cgel
      real*8 lambda,Mub
      real*8 raideur66(6,6),souplesse66(6,6)
      real*8 xmt
      logical dtiso
      real*8 nc33(3,3)
      real*8 dth00,dth0,poro1,vw1,dth1,CTHp1,CTHv1
      real*8 epsk006(6),epsm006(6),sigef06(6)
      real*8 sigke06(6),sigm06(6) 
      real*8 dtr00,dtr0,dfl00,dfl0
      real*8 bw0,pw0,bg0,pg0,bs0,ps0
      real*8 epspc600(6),epspc60(6),epspc6p(6),depspc6(6)
      real*8 epspt600(6),epspt60(6),epspt6p(6),depspt6(6),depspt6p(6)
      real*8 epsvt600(6),epsvt60(6),depsvt6(6) 
      real*8 epspg600(6),epspg60(6),epspg6p(6),depspg6(6)
      real*8 epsps600(6),epsps60(6),epsps6p(6),depsps6(6)
      real*8 epspc6(6),epspt6(6),epsvt6(6),epspt33p(3,3)
      real*8 epspt33(3,3),epspt3(3),vepspt33(3,3),vepspt33t(3,3)
      real*8 epspg6(6),epsps6(6)
      real*8 pw,bw,srw2,CWv2 
      real*8 epse06(6) 
      real*8 tauk2,taum2 
      real*8 dfin00,CMp00
      real*8 cc03(3)
      real*8 cke66(6,6),ckk66(6,6),ckm66(6,6)
      real*8 xkm 
      real*8 depse6(6),epse16(6)
      real*8 depsk6(6),epsk16(6)
      real*8 depsm6(6),epsm16(6) 
      real*8 sig16(6),sigke16(6)
      real*8 At,St,M1,E1,vdef0,E2,AtF,StF,M1F,E1F,vdefF,E2F
      integer err1
      real*8 hydr
      real*8 sig133(3,3),vsig133(3,3),vsig133t(3,3),sig16p(6)
      real*8 rts00,rts,rtw00,rtw,refw00,refw,rtg00,rtg
      real*8 hplt,hplg,hpls,mgel,mdef,bdef,dphidef  
      real*8 dth2,CTHp2,CTHv2,Cwp2,CTHp,CTHv
      real*8 srw1,CWp1,CWv1,CWv,CWp
      real*8 tauk,taum
      real*8 aar0,aar1,def0,def1,bg1,pg1,bs1,ps1
      real*8 epsk06(6),epsm06(6)
      integer na 
      real*8 depleqc,epleqc0,epleqc1
      real*8 dflu1 
      real*8 umdt,sigp,dti,dtr
      real*8 dt3(3),dgt3(3),dgc3(3),dwt3(3),dwc3(3),dst3(3),dsc3(3)
      real*8 sigf6(6),sigmf6(6)
      real*8 sigrm33(3,3),sigrf33(3,3),sigrf33p(3,3),sigrf6p(6)
      real*8 sigrm33p(3,3),sigrm6p(6),sigrm6(6)
      real*8 sigrh6p(6),sigrh33p(3,3),sigrh33(3,3),sigrh6(6),sigf6d(6)
      real*8 mdtt
      logical errr
      real*8 zero6(6),zero3(3)
      real*8 wplt06(6),wpltx06(6),wplt6(6),wpltx6(6)
      real*8 wpl3(3),vwpl33(3,3),vwpl33t(3,3)
      real*8 wplx3(3),vwplx33(3,3),vwplx33t(3,3)  
      real*8 dr3(3),nc3(3)
      real*8 errgf1,errgf0
      real*8 alphag,DCDW,alphas,rt_app,rc_app,tau_trac
      real*8 dtw0,dtg0,dts0,dcw0,dcg0,dcs0,kwrt,kwrc,href,gft_app
      real*8 dtl00,dtl0,dtl1 
      real*8 epspmf33(3,3),epspmf6(6)  
      real*8 bw1,pw1,mwat
      real*8 Kb        
c     eau dans les CSH et coeff de DTT
      real*8 wcsh0,wcshf,CDTT     
c     directions principales des increments de deformations
      real*8 deps33(3,3),deps3(3),vdeps33(3,3),vdeps33t(3,3)
      real*8 taumdtt
      real*8 coeff_grand,coeff_petit
      parameter (coeff_grand=1.0d20,coeff_petit=1.0d-20) 
c     poro meca / pressions
      real*8 bgel,dphigel,dphiw
c     increment des contraintes effectives      
      real*8 dsigef6(6) 
c     pression w csh
      real*8 pwcsh0,pwcshf,dpg,dps,dpw 
c     activite de la refermeture localisee 
      logical refermtl
c     endo pre pic compression
      real*8 dcpp,epc,xmc,dcpp00,dcpp0
      logical dciso
c     traiement anti battement
      logical actif0(nc),actif1(nc)
c     resistances effectives et endos aux pics
      real*8 rceff,epleqc00,rteff,dpict,dpicc,rteff_app,rceff_app
c     taux de franchissement des seuils par les contraintes capillaires
      real*8 fshr6(6),fshr33(3,3),fshr33p(3,3),fshr006(6),fshr06(6)
c     caracteristique ecoulement initial (liquide)
      real*8 taum_ja,youn_ja,nu_ja,ss_ja
c     contrainte  d endo hydrique
      real*8 skdw00,epse_tild6(6) 
c     exposant de couplage pour cwp
      real*8 T1,T2,Tref 
c     coeff de dilatation thermique des renforts
      real*8 dalpr
c     coeff drucker prager jeune age      
      real*8 dlja 
c     evolution jeune age Kelvin
      real*8 tauk_ja,tauk0,psi_ja,psi00  
c     resistances jeune age
      real*8 rcja,rtja
c     couplage des endommagements traction compression
      real*8 altc,dct3(3),dct0,dcc3(3),dcc0 
c     tau de cisaillement de Drucker Pragger
      real*8 taudp,taudpmax
c     indicateur de battement sur la resolution des criteres
      logical  battcr3d 
c     autorisation de resolution couplee en multi criteres      
      logical CPTCR(nc,nc)
c     numeros de la varible non locale a traiter
      integer nl,ivari
c     controle de l hypothese d ecoulement
      logical garder,travail_fissure,endo_interface_explicit
      parameter (endo_interface_explicit=.true.)   
c     calcul des deformations non locales
      real*8 eps_nl06(6),w_nl06(6),w_nlx06(6)      
      real*8 eps_nlf6(6),w_nlf6(6),w_nlxf6(6)
      real*8 deps_nl3(3),deps_nl6(6),deps_l6(6)
      real*8 deps_nl33(3,3),source_nl_pt3(3),source_nl_pt6(6)
      real*8 v_source_nl_pt33(3,3)
      real*8 source_nl_pt33(3,3),Depspl_nl6(6)
      real*8 eps_nl6(6),eps_nl33(3,3),eps_nl3(3)
      real*8 veps_nl33t(3,3),veps_nl33(3,3)
      real*8 eptmasq033(3,3),eptmasq33(3,3),eptmasq33p(3,3)
      real*8 eptmasq06(6),eptmasq6(6)
      real*8 alpha_nl
      real*8 xx33(3,3)
c     logique mise au point plastc
      logical log_plastc,log_plastt
 
c     romain declarations pour les fibres
c      include './declare_fibres.h'
-INC HDECFIB
 
c     variables supplementaires pour le traitement non local des increment de ept00
      real*8 depspt33(3,3),depspt3(3),vdepspt33(3,3),vdepspt33t(3,3)
 
c     variables pour les traitements non locaux specifiques a fldo3d
      logical log_H_TWL2
      integer Num_H_TWL2 
      logical log_H_DNL(6)
      integer Num_H_DNL(6)
c     autres variables pour non local
      logical ENDO_NL
      integer IVARNL
c     longueur caracteristique deduite de la dissusion et capaci te d Helm holtz (cf. charge_materiaux_helmholtz.h)
      real*8 LcH
      logical Method_H,Method_N
 
c***********************************************************************
c     initialisation des parametres 
 
c     nombre de renforts actifs
      NRENF00=int(max(xmat(NMAT1),0.d0))
 
c     ** indicateur de premier passage 
      if (var0(1).ne.1.) then
         ppas=.true.
      else
         ppas=.false.
      end if 
 
c     ******************************************************************
 
c     initialisation d une  matrice identité
       do i=1,3
         do j=1,3
            if(i.eq.j) then
                 vref33(i,j)=1.d0
                 vref33t(i,j)=1.d0
            else
                 vref33(i,j)=0.d0
                 vref33t(i,j)=0.d0 
            end if
        end do
      end do
c     initialisation de la variable de controle d erreur
      ierr1=0
c     coeff pour theta methode
      theta=0.5d0 
c     pour ne pas avoir de fluage mettre plast_seule à vrai
      plast_seule=.false. 
c     pour fluage 
      fl3d=.true.
c     pour endo      
      end3d=.true. 
c     detecteur de battement lors du retour plastique
      battcr3d=.false. 
 
 
c     *********** chargement des parametres materiaux depuis xmat ******
 
c     l hydratation est considéree en fin de pas pour ne pas avoir
c     a recuperer celle du pas precedent  les exposants de De Shutter
c     sont  dans hydramat, il faur recompiler pour les modifier
 
C       do i=1,nmat
C         print*,'xmat(',i,')=',xmat(i)
C       end do
C       read*
 
c     Module d young                            
      young00=xmat(nbelas3d+88)
      if( young00.eq.0.) then
c         .or.(IVALMAT3D(nbelas3d+88).eq.0)) 
c         Print*,'Donner Young YORF pour l hydratation de reference HREF'
          young00=xmat(1)
c          print*,'E(href)=',young00
c         err1=1
      end if      
c     coefficient de Poisson      
      nu00=xmat(nbelas3d+89)
      if( nu00.eq.0.) then
c      .or.(IVALMAT3D(nbelas3d+89).eq.0)) then
c         Print*,'Donner coefficient de Poisson NURF pour l hydratation'
c         print*,'de reference HREF'
c         err1=1
          nu00=xmat(2)
c          print*,'nu(href)=',nu00
      end if      
c      print*,'alpha fluendo3d=',alpha      
c     hydratation 
      hydr=xmat(nbelas3d+1)      
c     hydratation seuil
      hyds=xmat(nbelas3d+2)        
c     resistances et pression limite
      rt00=xmat(nbelas3d+3)
      if(rt00.eq.0.) then
       print*,'Il manque la resistance a la traction RT'
       err1=1
      end if
c     seuil pour la refermeture des fissures
      reft00=xmat(nbelas3d+4)
      if(reft00.lt.rt00/100.d0) then
         print*,'La contrainte de refermeture REF est trop faible'
         err1=1
      end if  
c     resistance en compression-par cisaillement
      rc00=xmat(nbelas3d+5)
      if(rc00.eq.0.) then
         print*,'Il manque la resistance en compression RC '
         err1=1
      end if       
c     coeff drucker prager
      delta00=xmat(nbelas3d+6)
      if(delta00.eq.0.) then
         print*,'Il manque DELT pour Drucker Prager  '
         err1=1
      end if          
c     coeff dilatance de cisaillement
      beta00=xmat(nbelas3d+7)
      if (beta00.eq.0.) then
         print*,'Il manque BETA pour Drucker Prager  '
         err1=1 
      end if           
c     deformation au pic de traction      
      ept00=xmat(nbelas3d+8)
      if (ept00.eq.0.) then
           ept00=rt00/young00 
      end if        
c     taux d ecrouissage des phases effectives / gel  
      hplg00=xmat(nbelas3d+9)
      if (hplg00.eq.0.) then
           hplg00=0.03*young00 
      end if        
c     coeff de concentration de contrainte capillaire pour le critere hydrique     
      cshr=xmat(nbelas3d+10)      
c     compressibilite du gel      
      kgel=xmat(nbelas3d+11)
c     energie de fissuration en traction directe
      gft00=xmat(nbelas3d+12)
      if (gft00.eq.0.) then
         print*,'Il manque l energie de fissuration GFT '
         err1=1
      end if 
c     deformation caracteristique du potentiel de fluage      
      epsm00=xmat(nbelas3d+13)
      if (epsm00.eq.0.) then
         print*,'Il manque le potentiel de fluage EPSM '
         err1=1
      end if       
c     raideur relative Kelvin / Young      
      psi00=xmat(nbelas3d+14)
      if (psi00.eq.0.) then
         print*,'Il manque le multiplicateur PSIK pour Kelvin '
         err1=1
      end if 
c     coeff amplicateur  à 2/3 rc pour le  fluage      
      xflu=xmat(nbelas3d+15)
      if (xflu.eq.0.) then
           xflu=1.0
      end if   
c     temps caracteristique pour Maxwell      
      taum00=xmat(nbelas3d+17)
      if(taum00.eq.0.) then
         print*,'Il manque le temps caracterisqtique de fluage TAUM'
         err1=1
      end if         
c     temps caracteristique pour Kelvin      
      tauk00=xmat(nbelas3d+16)
      if(tauk00.eq.0.) then
         print*,'Il manque le temps caracterisqtique de fluage TAUK'
         err1=1
      end if      
c     nrj activation du potentiel de fluage
      nrjm=xmat(nbelas3d+18)       
c     endommagement thermique a 80°C
      dt80=xmat(nbelas3d+19)  
c     porosite      
      poro2=xmat(nbelas3d+22)      
c     coeff de biot pour l eau
      biotw00=xmat(nbelas3d+20)     
c     module de biot pour le non sature
      xnsat00=xmat(nbelas3d+21) 
c     volume maximal de rgi
      vrag00=xmat(nbelas3d+23)     
c     volume d eau dans les pores
      vw2=xmat(nbelas3d+24) 
      if(vw2.eq.0.) then
        vw2=poro2
      end if        
c     coeff de couplage endo hydrique endo de compression
      DCDW=xmat(nbelas3d+25)
c     exposant de Van Genuchten pour la pression capillaire
      mvgn=xmat(nbelas3d+26)
      if(mvgn.eq.0.) then
         mvgn=0.5d0
      end if
c     deformation au pic de compression
      epc0=xmat(nbelas3d+27)
c     deformation cracateristique endo de compression
      ekdc=xmat(nbelas3d+28)
      if(ekdc.eq.0.) then
         ekdc=1.d0
      end if
c     deformation caracteristique pour l endo de rgi
      ekdg=xmat(nbelas3d+29) 
      if(ekdg.eq.0.) then
         ekdg=0.003d0
      end if
c     energie de fissuration pour l endo de traction
      gfr00=xmat(nbelas3d+30)  
      if(gfr00.eq.0.) then
         print*,'Il manque l energie de refermeture de fissure GFR'
         err1=1
      end if          
c     volume des vides accessibles pour les RGI
      vvgel=xmat(nbelas3d+31)
c     coeff de concentration de contrainte des RGI
      cgel=xmat(nbelas3d+32)
      if(cgel.eq.0.) then
c         par défaut on prend la solution de la sphère élastique
          cgel=2.0d0
      end if          
c     temperature de reference des parametres materiaux (celsius)      
      tetar=xmat(nbelas3d+33)         
c     temperature seuil pour l endo thermique (celsius)      
      tetas=xmat(nbelas3d+34) 
c     endommagement maximum par fluage
      dfmx=xmat(nbelas3d+35)
c     temps caracterisqtique de la rag à tref      
      taar=xmat(nbelas3d+36)
c     nrj d'activation de la rgi
      nrjg=xmat(nbelas3d+37)
c     seuil de saturation minimal pour avoir la rgi
      srsrag=xmat(nbelas3d+38)     
c     constante de calage de la defom therm transitoire
      mdtt00=xmat(nbelas3d+39)
c     volume de reference pour Weibull
      vref=xmat(nbelas3d+40)
c     volume d integration maxi pour Weibull
      vmax=xmat(nbelas3d+41)
      if(vmax.eq.0.) then
c         on met VMAX = VREF pour emepcher l effet d echelle probabiliste
          VMAX=VREF
      end if
c     coefficient de variation de la resistance a la traction
      cvrt00=xmat(nbelas3d+42)         
c     temps cracteristique pour la def      
      tdef=xmat(nbelas3d+43)
c     energie d activation de precipitation de la def
      nrjp=xmat(nbelas3d+44)
c     seuil de saturation pour declancher la def
      srsdef=xmat(nbelas3d+45)
c     quantite maximale de def pouvant etre realise      
      vdef00=xmat(nbelas3d+46)       
c     teneur en alcalin pour la def
      cna=xmat(nbelas3d+47)
c     temperature de reference pour la precipitation de la def (Celsius)
      ttrp=xmat(nbelas3d+57)
c     energie d activation des processus de dissolution des phases primaires
      nrjd=xmat(nbelas3d+56)     
c     temps caracteristique pour la fixation des aluminiums en temperature
      tfid=xmat(nbelas3d+55)
c     temps caracteristique pour la dissolution de l ettringite primaire
      tdid=xmat(nbelas3d+54)      
c     temperature de reference pour la dissolution de la def (Celsius)
      ttdd=xmat(nbelas3d+53)
c     exposant de la loi de couplage precipitation def alcalins     
      exmd=xmat(nbelas3d+52)
c     exposant de la loi de couplage temperature de dissolution alcalins      
      exnd=xmat(nbelas3d+51)
c     concentration en alcalin de blocage de la def      
      cnab=xmat(nbelas3d+50)
c     concentration caracteristique pour les lois de couplage de la def      
      cnak=xmat(nbelas3d+49)
c     rapport molaire S03 / Al2O3 du ciment      
      ssad=xmat(nbelas3d+48) 
c     defcaracteristique pour l endo du a la def
      ekds=xmat(nbelas3d+58)   
c     contrainte caracteristique pour l endo hydrique
      skdw00=xmat(nbelas3d+59) 
*  pourquoi skdw ne pourrait pas etre nul?
**    if(skdw00.eq.0.) then
**       print*,'Il manque la contrainte SKDW pour l endo capillaire '
**       err1=1
**    end if
c     tetak pour exposant de van genuchten en fonction de teta 
      tkvg=xmat(nbelas3d+60)
c     temperature minimale pour la fixation des alu en cas de def
      ttdf=xmat(nbelas3d+61)
c     energie activation pour la vitesse de fixation des alu en cas de def      
      nrjf=xmat(nbelas3d+62) 
c     nombre de moles de sulfates par unite de volume
      nsul=xmat(nbelas3d+63)
c     temps caracteristique pour l amplification dynamique
      tdyn=xmat(nbelas3d+64)
c     temperature caracteristique pour la rag (le gel)
      ttrg=xmat(nbelas3d+65)
c     module d ecrouissage pour la resistance a la pression de DEF      
      hpls00=xmat(nbelas3d+66)
c     module de compressibilite pour la def
      kdef=xmat(nbelas3d+67)
c     volume des vides pour la def      
      Vvdef=xmat(nbelas3d+68)
c     coeff de concentration de contrainte pour la def      
      cdef=xmat(nbelas3d+69)
c     module d crouissage pour la resistance effective a la pression d eau
      hplw00=xmat(nbelas3d+70)
c     temperature de reference pour les coeffs de l isotherme
      tetarw=xmat(nbelas3d+71)         
c     coeff de couplage  endo de rag endo de compression
      alphag=xmat(nbelas3d+72) 
c     coeff de couplage enddo de def endo de compression  
      alphas=xmat(nbelas3d+73)     
c     effet capillaire sur traction apparente
      kwrt=xmat(nbelas3d+74)
c     effet capillaire sur compression apparente
      kwrc=xmat(nbelas3d+75)     
c     hydratation de reference pour la donnee des parametres mecaiques
      href=xmat(nbelas3d+76)
      if(href.eq.0.) then
         print*,'le degre d hydratation de reference HREF est nul '
         err1=1
      end if          
c     temps caracteristique pour la deformation thermique transitoire
      tdtt=xmat(nbelas3d+77)  
c     eau pour la DTT
      wdtt=xmat(nbelas3d+78) 
c     pression d eau de reference dans les csh
      pdtt=xmat(nbelas3d+79)   
c     contrainte seuil pour initier le fluage de maxwell
      ss_ja=xmat(nbelas3d+80)
c     temps caracteristique etat liquide
      taum_ja=xmat(nbelas3d+81) 
c     young jeune age
      youn_ja=xmat(nbelas3d+82)  
c     poisson jeune age
      nu_ja=xmat(nbelas3d+83)
c     delta drucker prager jeune age
      dlja=xmat(nbelas3d+84) 
c     rc jeune age
      rcja=xmat(nbelas3d+85)
c     rt jeune age
      rtja=xmat(nbelas3d+86) 
c     supplement de dilatation thermique des acier (alphar-alpha)
      dalpr=xmat(nbelas3d+87)       
c     couplage traction cisaillement
      altc =xmat(nbelas3d+90) 
c     avancement latent rag
      alatr=xmat(nbelas3d+91)
c     seuil de deformation pour endo de rag
      epssr=xmat(nbelas3d+92)      
c     module d ecrouissage traction directe 
      hplt00=0.d0
c     nrj activation pour l eau
      nrjw=17000.d0 
c     amplitude reversible du fluage jeune age
      psi_ja=psi00*coeff_grand
c     temps de kelvin jeune age
      tauk_ja=taum_ja*coeff_grand
c     recuperation des parametres pour les fibres
c     include './charge_materiaux_fibres.h'
-INC HMATFIB
 
c     ***********test de conformite sur les donnees pour la matrice ****
 
      if(err1.eq.1) then
         print*,'Incoherence dans les donnees de fluendo3d'
         ierr1=1
         return
      end if
 
c     ** verif validite de la deformation au pic de traction
      ept00=max(rt00/young00,ept00)
 
c     ** verif coherence deformation pic de compression 
      epc1=rc00/young00
c     eps pic comp ne peut pas etre inferieur à rc/E      
      epc00=max(epc0,epc1)      
 
c     ** verif validite de la dilatance
      if(beta00.gt.sqrt(3.d0)) then
        print*,'Dilatance BETA excessive dans fluendo3d'
        ierr1=1
        return
      end if
 
c     ** verif presence d un ecrouissage pour la rag      
      if(hplg00.le.0.d0) then
           if(vrag00.gt.0.) then
             print*,'Taux d ecrouissage HRAG doit etre > 0'
             ierr1=1
             return
           end if
      end if 
 
c     ** verif presence d un ecrouissage pour la def      
      if(hpls00.le.0.d0) then
           if((vdef00.gt.0.).or.(nsul.gt.0.)) then
             print*,'Taux d ecrouissage HDEF doit etre > 0'
             ierr1=1
             return
           end if
      end if 
 
c     ** test nbr maxi de criteres actifs pour le couplage
      if(nc.ne.13) then
          err1=1   
          print*,'pb nbr max de criteres ds fldo3d .ne. 13'
          return
      end if 
 
c     *** initialisation des parametres calculees a partir des donnees *
 
c     ** stockage hydratation 
      if (ppas) then
c       au 1er passage on suppose hyd0=hydr
        hyd0=hydr
c       on initialise var03d en cas de sous incrementation par fluage        
        var0(68)=hyd0
      else
c       on recupere l hydratation du pas precedent      
        hyd0=var0(68)
      end if
      varf(68)=hydr  
 
c     ** initialisation des variables internes associee au volume d eau
 
c     eau et porosite (initialise debut=fin 1er pas)
      if (ppas) then
         var0(69)=vw2
         var0(70)=poro2
      end if
c     stockage pour le prochain pas 
      varf(69)=vw2       
      varf(70)=poro2
 
c     tables a zero
      do j=1,6
         zero6(j)=0.d0
      end do
      do j=1,3
         zero3(j)=0.d0
      end do
 
c****************** parametres pour Helmholtz **************************
c     lignes a inclure dans les modeles utilisant des formulations d helmholtz
c      include './charge_materiaux_helmholtz.h'
-INC HMATHEL    
      do NL=1,NB_HELM
c       lecture des variable internes et rangements dans Helm0(NL,ii)
c       et Helm1(NL,ii) avec ii=1,NB_VARI_PAR_HELM
c        include './lire_vari_helmholtz.h'
-INC HLVIHEL
      end do
c***********************************************************************
 
 
c*************** parametres pour les renforts **************************
c     lignes a inclure dans les modeles utilisant des renforts
c      include './charge_materiaux_renforts.h'
-INC HMATREN
c*********************************************************************** 
 
c********** traitement des autorisations d activation des criteres *****
      do i=1,nc
        do j=1,nc
c         initialisation de la compatibilite des criteres a faux
          CPTCR(i,j)=.false. 
          if(i.ne.j) then 
             if ((i.le.9).and.(j.le.9)) then
c                tous les criteres de pression sont inter compatibles             
                 CPTCR(i,j)=.true.
             else if ((i.le.12).and.(j.le.12).and.
     #                (i.ge.10).and.(j.ge.10)) then
c                les macro fissures sont compatibles entre elles
                 CPTCR(i,j)=.true.
             else if((i.eq.13).and.(j.ge.10).and.(j.le.12)) then
c                DP et macros fissures incompatibles             
                 CPTCR(i,j)=.false.
             else
c                autres couplages interdits              
                 CPTCR(i,j)=.false.            
             end if
          else
c            cas des criteres seuls  
             CPTCR(i,i)=.true.  
          end if
c         on desactive tous les criteres
c         CPTCR(i,j)=.false.          
        end do
      end do      
 
c     ******************************************************************
c     affichage des options du calcul
c     print*,'fluage',fl3d
c     print*,'endo',end3d    
c***********************************************************************      
 
 
c***********************************************************************
c     chargement increment de deformation imposee et deformation finale
c     remarque 4-5-6 sont des gama jusqu ici, ensuite des epsilons
      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 gama 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*************** traitement non local de TWL2(VARI numero 72) **********
      ivari_twl2=72
c     on teste l existence d un champ de phase pour cette variable et
c     on recupere son numero Num_H_TWL2
      call exhmtz(istep,NBVIA3D,INLVIA3D,NB_HELM,
     #ivari_twl2,log_H_TWL2,Num_H_twl2)
c     si le champ de phase existe      
      if(log_H_TWL2) then
c       calcul approché de l'exposant de Weibull  (m)    
        if(cvrt00.ne.0.) then
             if (cvrt00.gt.1.) then
                print*,'Dans fluendo3d, erreur de donnees car le '
                print*,'coeff de variation de Rt pour Weibull est >1 '
                ierr1=1
                return
             end if        
             eweb=1.2d0/cvrt00-0.2d0
c            print*,'eweb=',eweb 
        else
            print*,'Donner le coeff de variation de Rt pour Weibull'
            print*,'Mettre Vmax a zero pour annuler  l effet de WL2'
            print*,'ou bien Vmax a Vref pour utiliser WL2'
            ierr1=1
            return
        end if   
c       print*,'Etape 2-1 ds fldo3d Istep ',istep
c       istep=1 : 1er passage pour preparer le calcul non local
        if(istep.eq.1) then
c         avec cwrt du pas precedent et contraintes pas actuel
          if (ppas.or.(hydr.le.hyds)) then
             if ((vmax.ne.0.).and.(eweb.ne.0.).and.(vref.ne.0.)) then
                 xweb=1.d0/eweb
                 cwrt=(vref/vmax)**xweb
             else                
                 cwrt=1.d0
             end if
c            1er passage il faut initialiser cwrt             
             var0(76)=cwrt
          else
c            on prend le dernier cwrt estime          
             cwrt=var0(76)                    
          end if
c         on stocke cwrt dans varf pour les sous increments locaux
          varf(76)=cwrt
c         variable pour le stockage du max de xwl2 non local (mxwl ds idvar4)
c        (a mettre absolument a un car utilise avec cette hypothese dans unpas)          
          varf(73)=1.d0 
c         variable pour le stockage du taux de charge maxi sur la structure (mtwl ds idvar4)
c        (a mettre absolument a un car utilise avec cette hypothese dans wl2d3d)   
          varf(74)=1.d0     
        else 
c         nouveau passage non local             
c         *** recup variable non locales WL2 *************************
c         recup de xwl2=//(s/r)**m .psi .dv=Veq.(max(s)/rt)**m
          xwl2=var0(ivari_twl2)
c         lecture de vari supplementaires associees a helmholtz
          ierr1=1
          return 
          if(ierr1.eq.1) then
            print*,'Pb calcul Helmholtz twl2 ds fldo3d'
            return
          else
c           lecture des variables internes complementaires pour Helmholtz
            NL=Num_H_TWL2
c           IVARI=IVARI
c            include './lire_vari_helmholtz.h'
-INC HLVIHEL
c           les var0 de cette formulation sont alors dans HELM0(Nl,i)
          end if
c         recup du  maxi sur le modele de sigma/rt_app            
          tau_trac_wl2=var0(72) 
c         recup du maxi sur le modele xwl2 non local 
          xwl2max=var0(73)    
c         tau de charge maxi
          tau_trac_wl2_max=var0(74)
c         indicateur de localisation de la fissure
c         inutile de l affecter par l hydratation cela a ete fait lorsque istep=1
          dtl0=var0(75)
c         endommagement localise probabble fin de pas (estime pour istep=1)
          dtl1=var0(119)            
c         coeff de Weibull rt pas precedent
          cwrt=var0(76)         
c         pression capillaire initiale
          pw0=var0(81)
c         saturation en eau * biot            
          bw0=var0(80)
c         ** calcul de la nouvelle resistance a la traction **********
c         contrainte et resistance ne sont pas utilisees a cette etape
c         on les met  a zero
          do i=1,3
               sig13(i)=0.d0
          end do
          call wl2d3d(2,vmax,vref,sig13,eweb,0.d0,cwrt,
     #    xwl2,tau_trac_wl2,dtl0,dtl1,xwl2max,tau_trac_wl2_max,
     #    bw0,pw0,kwrt)
c         on stocke les variables non locales pour verif eventuelle
          varf(71)=xwl2
c         variable non locale : sigm_maxi  /rt          
          varf(72)=tau_trac_wl2
c         stockage de la valeur maxi de smax /rt
          varf(74)=tau_trac_wl2_max
c         endo fin de pas            
c         on copie endommagement fin de pas utilisee pour cwrt pour le
c         debut de pas suivant car dtl calcule en fin de programme 
c         n est pas celui utilise pour calculer cwrt
c         DTRA debut de pas
          varf(75)=dtl1 
c         DTL0
          varf(119)=varf(75)            
c         coeff d effet d echelle sur rt 
c         print*,'cwrt ds fldo3d apres istep=2',cwrt           
          varf(76)=cwrt  
c         ecriture des variables internes HELM1 associees a la formulation NL de Helmholtz
c          include './ecrire_vari_helmholtz.h' 
-INC HEVIHEL                              
        end if                   
      else
c       pas de calcul non local pour wl2
c       ** rt00 n est pas affectee par la methode WL2 **************          
c       istep=0 :le calcul sans effet d echelle Weibull
        cwrt=1.d0
c       liste des variables non locales : XWL2
        varf(71)=0.d0
c       variable non locale : sigma maxi /rt   TWL2        
        varf(72)=0.d0  
c       max de xwl2 non utilisse on y laisse 1. MXWL
        varf(73)=1.d0
c       contrainte relative maxi non locale non utiliseee MTWL
        varf(74)=0.d0 
c       endommagement maxi de traction DTMX
        varf(75)=var0(75)
c       modif de rt par effet d echelle inactif     CWRT        
        varf(76)=cwrt  
      end if            
 
c***********************************************************************            
c     influence de l hydratation sur les parametres materiau  
c     les resistances locales sont prise egales à la resistance 
c     macroscopique idem pour les contraintes de refermeture
      rts00=rt00     
      rtg00=rt00
      rtw00=rt00
      refw00=reft00
c     pas d effet d hydratation sur coeff de def therm transitoire
      mdtt=mdtt00 
c     controle calcul raideur iso      
      if(.not.iso) then
         print*,' cas imprevu dans fluendo3d '      
         print*,' matrice de raideur initiale non isotrope'
         ierr1=1
         return
      end if 
c     modification des parametres materiau en fonction de l hydratation              
      call hydm3d(hyd0,hydr,hyds,young00,young,nu00,nu,rt00,rt,
     #rteff,reft00,reft,rc00,rc,rceff,delta00,delta,beta00,beta,gft00,
     #gft,ept00,ept,epsm00,epsm,xnsat00,xnsat,biotw00,biotw,iso,lambda,
     #Mub,Kb,rtg00,rtg,raideur66,souplesse66,xmt,dtiso,cwrt,err1,teta2,
     #href,tetar,tkvg,gfr00,gfr,rts00,rts,rtw00,rtw,hplt00,hplt,hplg00,
     #hplg,hpls00,hpls,epc00,epc,xmc,dciso,dpict,dpicc,taum_ja,youn_ja,
     #nu_ja,ss_ja,ssfl,taum00,taum0,skdw00,skdw,dlja,tauk_ja,
     #tauk00,tauk0,psi_ja,psi00,psik,rcja,rtja)      
c     traitement de l erreur eventuelle 
      if(err1.eq.1) then
        print*,'Pb lors de la mise a jour avec hydramat ds as3d'
        ierr1=1
        return
      end if     
c      print*,'istep',istep,'rt',rt         
c      print*,'Etape 4 ds fldo3d apres hydramat'         
 
c***********************************************************************
c     Nombre maxi de fissures localisees en fonction de la direction
      if(.not.ppas) then
c        chargement des longueurs des elements dans les directions des
c        renforts   
         do i=1,NRENF00
            lcr(i)=var0(NVAR2+(i-1)*NB_VARI_PAR_RENF+13)
         end do
      end if
c     calcul du nb max de fissures sur la zone d integration     
      call nfin3d(ppas,lcr,nc33,NRENF00,vecr,NB_RENF,longr,lsr,
     # lfr,deqr,syr,taur,rhor,rt00,XE3D,NBNMAX3D,NBNB3D,IDIMB3D,
     # NB_HELM,TAILH,log_H_RENF,Num_H_RENF,err1)
      if(err1.eq.1) then
         print*,'Erreur lors de l appel a nfin3d dans fldo3d'
         print*,'pour istep=',istep
         ierr1=1
         return
      end if
c     sauvegarde des longueurs des elements ds les directions des renforts     
      do i=1,NRENF00
        varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+13)=longr(i)
        if(longr(i).eq.0.) then
             print*,'Ds fldo3d longr',i,'=', longr(i)
             print*,' istep',istep,'ppas',ppas
             ierr1=1
             return
        end if
      end do
c     print*,'Etape 5 ds fldo3d apres nfin3d' 
 
c***********************************************************************        
c     chargement des principales variables internes 
c     (etat du squelette solide)
       do i=1,6
          epsk006(i)=var0(i+7)
          epsm006(i)=var0(i+13)
          epsvt600(i)=var0(31+i)              
          sigef06(i)=var0(i+49)
          sigke06(i)=var0(i+55)
          sigm06(i)=var0(i+61)            
       end do 
c      endommagement effectif par fluage 
       dfl00=var0(78)        
c      endommagement thermique 
       dth00=var0(77)  
c      endommagement iso prepic de traction (DTPP)
       dtr00=var0(79)
c      endo iso pre pic de cisaillement (DCPP)
       dcpp00=var0(122)         
c      pression capillaire          
       bw0=var0(80)
       pw0=var0(81)
c      pression gel         
       bg0=var0(82)
       pg0=var0(83)
c      pression def
       bs0=var0(84)
       ps0=var0(85)         
c      tenseurs des deformations plastiques et criteres hydrique         
       do j=1,6 
c         deformations plastiques  
          epspc600(j)=var0(19+j)         
          epspt600(j)=var0(25+j)
          epspg600(j)=var0(37+j) 
          epsps600(j)=var0(43+j)
c         critere pression capillaire              
          fshr006(j)=var0(129+j)
       end do 
c      def plastique equivalente Drucker Prager         
       epleqc00=var0(94)
 
 
c***********************************************************************         
c      influence du degre d hydratation sur les variables internes         
       call hydv3d(ierr1,hyd0,hydr,hyds,dth00,dth0,dtr00,dtr0,
     # dfl00,dfl0,dcpp00,dcpp0,epspt600,epspt60,epspc600,epspc60,
     # epsvt600,epsvt60,epspg600,epspg60,epsps600,epsps60,epsk006,
     # epsk06,epsm006,epsm06,epleqc00,epleqc0,fshr006,fshr06)
c      print*,'Etape 6 ds fldo3d apres hydv3d' 
c      print*,'deformation totale initiale'
c      call afic6(epspt60)
 
c***********************************************************************
c      increments non locaux des deformations plastiques de traction 
c      (DPTi=VARI numero 140+i)
       ENDO_NL=.true.
       do i=1,6
c           numero de la variable interne a traiter par Helmholtz (a partie de idvar4)
            ivarinl=140+i
c           test du traitement non local + assigantion logique log_H_DNL(i)
c           et recuperation du numero de la formulation non locale associee
c           dans Num_H_DNL(i) 
            call exhmtz(istep,NBVIA3D,INLVIA3D,NB_HELM,
     #      ivarinl,log_H_DNL(i),Num_H_DNL(i))
c           on verifie que les 3 deformations plastiques principales de
c           traction soient traitees par la formulation d Helmholtz
            ENDO_NL=(ENDO_NL.and.log_H_DNL(i))
       end do
c      print*,'fluendo ENDO_NL etape 1',ENDO_NL
 
       if(ENDO_NL) then
c        longueur non locale moyenne pour l endommagement
         LcH=1.d0
         do i=1,6
             LcH=LcH*max(LcarH(Num_H_DNL(i),1),LcarH(Num_H_DNL(i),2),
     #       LcarH(Num_H_DNL(i),3))
         end do
         LcH=LcH**(1./6.)
c        indicateur de la methode de calcul des dimensions           
         Method_H=.true.
         Method_N=.false.
C        print*,'LcH',LcH,'Method_H',Method_H,'Method_N',Method_N
C        read*
c        on est dans la procedure incrementale non locale
         do i=1,6             
c           position de la variable non locale
            ivarinl=140+i  
c           recuperation de la valeur non locale
            eps_nl6(i)=var0(ivarinl)
         end do 
      else
         LcH=1.d0
         Method_H=.false.
         Method_N=.true.
      end if
 
c***********************************************************************
c     influence de la temperature sur les parametres  materiau et
c     actualisation de l endo thermique debut de pas 
      poro1=var0(70)    
      vw1=var0(69)
      call thrm3d(teta1,nrjm,tetas,tetar,DT80,dth0,
     # dth1,CTHp1,CTHv1,poro1,vw1,nrjw)
c     endommagement thermique en fin de pas     
      call thrm3d(teta2,nrjm,tetas,tetar,DT80,dth1,
     # dth2,CTHp2,CTHv2,poro2,vw2,nrjw)
      varf(77)=dth1 
c     calcul des coeffs moyens pour le potentiel de fluage
      CTHp=0.5d0*(CTHp1+CTHp2)
      CTHv=0.5d0*(CTHv1+CTHv2)
 
 
c***********************************************************************
c     calcul des avancements physico-chimiques
c***********************************************************************
 
c     **** calcul de la pression capillaire due a l eau en fin de pas **
 
      call watr3d(mfr,biotw,poro2,vw2,xnsat,mvgn,pw,bw,srw2,teta2,
     # tetarw,tkvg,pw0,Mwat,dphiw,bw0,bw1)
      varf(80)=bw1
 
c     modif eventuelle des viscosites de Kelvin en fonction de srw2
      if(poro1.ne.0.) then
            srw1=vw1/poro1
      else
            srw1=1.d0
      end if
      CWv1=1.d0/srw1       
      CWv2=1.d0/srw2
      CWv=0.5d0*(CWv1+CWv2)
c     modif potentiel de fluage de Maxwell en fonction de srw
      Cwp2=srw2
      CWp1=srw1
      CWp=0.5d0*(CWp1+Cwp2) 
C c     pas d effet de l eau sur le potentiel 
C       Cwp2=1.d0
C       CWp1=1.d0
C       Cwp=1.d0      
 
c     *** cas de l eau des CSH ****************************************
 
      if(ppas) then
c         initialisation de la quantite d'eau dans les CSH a la valeur 
c         de reference
          wcsh0=wdtt
          pwcsh0=0.d0
      else
          wcsh0=var0(120)
          pwcsh0=var0(121)
      end if
c     prise en compte de la pression dans les csh pour la DTT
c     cf. F. Manzoni, T. Vidal, A. Sellier, X. Bourbon, and G. Camps,
c     “On the origins of transient thermal deformation of concrete,” 
c     Cem. Concr. Compos., vol. 107, 2020, 
c     https://doi.org/10.1016/j.cemconcomp.2019.103508       
      call wdtt3d(wcsh0,wcshf,tdtt,mdtt,wdtt,pdtt,CTHv,
     # teta1,teta2,pw,dt,CDTT,pwcsh0,pwcshf,tetarw,err1)
      if(err1.eq.1) then
           print*,'Pb lor du calcul de la DTT dans fldo3d'
           ierr1=1
           return
      end if
      varf(120)=wcshf
      varf(121)=pwcshf 
      varf(129)=cdtt 
 
c     ** Modification eventuelle de la viscosite de Kelvin (T,Sr) ******
 
      tauk=tauk0*CWv/CTHv
c     print*,'fldo3d tauk0',tauk0,cwv,cthv
 
c     ** Modification de la viscosite de Maxwell (T,Sr) **************** 
 
c     pour Sr, deja integre dans CC via Cwp
      taum=taum0/CTHv   
 
c     ** Modification de la viscosite de Maxwell (T,Sr) pour la DTT ****
 
      if((CDTT.gt.coeff_petit).and.(hydr.gt.hyds)) then
c            dtt active       
c             print*,'Fluendo3d CDTT',CDTT
             taumdtt=taum/CDTT
      else
c            pas de dtt avent percolation       
             taumdtt=taum*coeff_grand
      end if       
c      print*,'ds fluando3d taum=',taum
c      print*,'ds fluando3d taumdtt=',taumdtt     
 
c     **** avancement de la rag ****************************************
 
      aar0=var0(86)
      bg0=var0(82)
      call rag3d(taar,nrjg,ttrg,aar0,srw2,srsrag,teta1,
     # teta2,dt,aar1,vrag00,Kb,Kgel,bgel,mgel,dphigel,bg0,bg1,
     # vvgel,cgel,Rt,alatr,err1)
      if(err1.eq.1) then
            print*,'Pb lors du calcul de rag dans fldo3d'
            ierr1=1
            return
      end if
c      print*,'fldo3d',dphigel,bgel,mgel
c      read*
      varf(86)=aar1
      varf(82)=bg1
 
 
c     **** avancement de la def  ***************************************
 
      bs0=var0(84)
      E1=var0(87)
      M1=var0(88)
      At=var0(89)
      St=var0(90)         
      E2=var0(91)
      def0=var0(92)
      vdef0=var0(93)
      if((vdef00.gt.0.).or.(nsul.gt.0.)) then
        call def3d(ppas,tdef,nrjd,def0,srw2,srsdef,teta1,teta2,dt,
     #  vdef00,def1,CNa,nrjp,ttrp,tfid,ttdd,tdid,exmd,exnd,
     #  cnab,cnak,ssad,At,St,M1,E1,vdef0,E2,AtF,StF,M1F,E1F,vdefF,E2F,
     #  ttdf,nrjf,nsul,ierr1,Kb,Kdef,Vvdef,Cdef,Rt,mdef,bdef,dphidef,
     #  bs0,bs1) 
        if(ierr1.eq.1) then
            print*,'pb lors du calcul de la def dans fldo3d'
            return
        end if
      else
        bs1=bs0
        E1F=E1
        M1F=M1
        AtF=At
        StF=St
        E2F=E2
        def1=def0
        vdefF=vdef0
        bdef=0.d0
        dphidef=0.d0
        mdef=0.d0
      end if
      varf(84)=bs1
      varf(87)=E1F
      varf(88)=M1F
      varf(89)=AtF
      varf(90)=StF
      varf(91)=E2F
      varf(92)=def1
      varf(93)=vdefF
 
c     ******************************************************************      
c     prise en compte de l effet de pression (coeffs kwrt et kwrc)
c     et de l effet d echelle probabiliste (coeff cwrt) sur la  
c     resistance apparente de traction et en compression         
      rt_app=(rteff-bw*pw*kwrt)*cwrt*(1.d0-dpict)
      rc_app=(rceff-bw*pw*kwrc)*(1.d0-dpicc)
c     deformation elastique de reference  pour le fluage prise a 1/3 
c     de rc sat effective
      epser=(rc00/young00)/3.d0      
 
c     ****************************************************************** 
c     reevaluation de la deformation compatible avec l etat
c     de contrainte debut de pas  
      call hydc3d(souplesse66,sigke06,epsk06,psik,
     # epse06,sigef06,fl3d,plast_seule)
c      print*,'Etape 7 ds fldo3d apres hydc3d'      
 
c***********************************************************************
c     tir visco elastique
c***********************************************************************
 
c     ** initialisation des deformations *******************************
      do i=1,6
c        increment deformation elastique            
         epse16(i)=epse06(i)      
c        increment deformation kelvin         
         epsk16(i)=epsk06(i)
c        deformation maxwell           
         epsm16(i)=epsm06(i)
c        deformation visqueuse transitoire
         epsvt6(i)=epsvt60(i)
c        deformations plastiques cisaillement
         epspc6(i)=epspc60(i)
c        deformations plastiques traction 
         epspt6(i)=epspt60(i)
c        deformations plastiques gel       
         epspg6(i)=epspg60(i)
c        deformations plastiques def
         epsps6(i)=epsps60(i)
      end do
 
 
c     ** directions principales des deformations imposees et matrice de passage associee
      call x6x33(deps6,deps33)
      call b3_v33(deps33,deps3,vdeps33)
c      do i=1,3
c         print*,'fldo3d deps3:',i,deps3(i)
c      end do
c     construction matrice de passage inverse         
      call traps1(vdeps33t,vdeps33,3)  
 
c     ** effet du chargement initial sur le potentiel de fluage *******
      call dffn3d(sigm06,delta,rc_app,rt_app,xflu,
     # ssfl,dfin00,CMp00,dfmx,hydr,hyds,err1)
      if(err1.eq.1)then
         print*,'Erreur lors du calcul de CMP dans dffn3d fldo3d'
         ierr1=1
         return
      end if
c      print*,'fldo3d cmp=',cmp00 
 
c     ** effet de l endommagement capillaire sur le potentiel de fluage
      call ampl3d(souplesse66,sigm06,fshr06,skdw,DCDW,epse06,sigef06,
     # epse_tild6,refw)     
 
c     ** actualisation des coeffs de consolidation avec chargement ****
 
C     cf. Alain Sellier, Stéphane Multon, Laurie Buffo-Lacarrière, Thierry Vidal, Xavier Bourbon, Guillaume Camps,
C     Concrete creep modelling for structural applications: non-linearity, multi-axiality, hydration, temperature and drying effects,
C     Cement and Concrete Research, Volume 79, 2016,Pages 301-315,ISSN 0008-8846,
C     https://doi.org/10.1016/j.cemconres.2015.10.001.
 
      call conso3d(xkm,epsm,epser,epsm06,epse_tild6,cc03,
     # vdeps33,CMp00,CTHp,CTHv,Cwp)
c      do i=1,3
c         print*,'fldo3d cc3:',i,cc03(i)
c      end do
 
c     ** endommagement par fluage ************************************     
      call dflu3d(xkm,cc03,dfl0,dflu1,dfin00)
c     stockage endo de fluage
      varf(78)=dflu1
 
c     ** initialisation des variables pour le tir visco elastique ****
c     pour le fluage
      do i=1,3
         alphadp(i)=0.d0
         betadp(i)=0.d0
      end do
c     pour les criteres
      do i=1,nc
         actif(i)=.false.
         complet(i)=.true.
         factif(i)=0.d0
      end do  
 
c     ** initialisation du nombre de criteres actifs (pas de criteres actifs lors du tir)      
      NA=0        
 
c     ** assemblage des matrices de couplage de fluage ds base fixe      
      call coupl3d(ann,xn,bn,ngf,err1,deps3,vdeps33,vdeps33t,dt,
     # cc03,taum,taumdtt,tauk,psik,epse06,epsk06,raideur66,Mgel,bgel,
     # dphigel,Cgel,Hplg,Mdef,bdef,dphidef,Cdef,Hpls,Mwat,bw,dphiw,
     # cshr,alphadp,betadp,actif,factif,nc,dsigef6,depse6,depsk6,
     # depsm6,depspc6,depspt6,depsvt6,depspg6,depsps6,.true.,ipzero,
     # dpg,dps,dpw,zero6,.false.,complet,NA,CMp00)
 
c     ** traitement de l erreur eventuelle lors du tir     
      if(err1.eq.1) then
          print*,'pb assemblage couplagf3d ds fldo3d'
          ierr1=1
          return
      end if 
c     print*,'Etape 8 ds fldo3d apres coupl3d'       
 
c     ** recuperation des increments et actualisation des deformations 
      do i=1,6
c        increment deformation elastique            
         epse16(i)=epse16(i)+depse6(i)       
c        increment deformation kelvin         
         epsk16(i)=epsk16(i)+depsk6(i) 
c        deformation maxwell           
         epsm16(i)=epsm16(i)+depsm6(i)
c        deformation visqueuse transitoire (dtt)
         epsvt6(i)=epsvt6(i)+depsvt6(i)
      end do
 
c     ******* recuperation deformation equivalente de compression ******          
      epleqc1=epleqc0
 
c     ***** actualisation des pressions ********************************
c     gel rag
      pg1=pg0+dpg
c     def        
      ps1=ps0+dps      
c     depression capillaire        
      pw1=pw0+dpw             
 
c     ***** etat des contraintes effectives apres tir visco elastique **
      do i=1,6
         sig16(i)=sigef06(i)+dsigef6(i)
         sigke16(i)=sigke06(i)
         do j=1,6
             sigke16(i)=sigke16(i)+raideur66(i,j)*depsk6(j)*psik
         end do
C          print*,'Fluendo3d etape 8-2 de tir visco elastique'
C          print*,'sig16(',i,')=',sig16(i)
C          print*,'sigke16(',i,')=',sigke16(i)
C          print*,'deps6(',i,')=',deps6(i)
C          print*,'depse6(',i,')=',depse6(i)
C          print*,'depsk6(',i,')=',depsk6(i)
C          print*,'depsm6(',i,')=',depsm6(i)     
      end do
 
 
c***********************************************************************
c     test des criteres de plasticité
c***********************************************************************
 
c     initialisation du compteur de sousiterations locales
      NCAL=0
c     initialisation indicateur de plasticite en cisaillement
      log_plastc=.false.       
 
c     ** increment des iteration locales        
20    NCAL=NCAL+1
c     print*,'******numero sous iter retour radial ds Fluendo3d:',ncal
      if(NCAL.eq.1) then
c           diagonalisation du vecteur des contraintes 
            call x6x33(sig16,sig133)        
c           print*,'fldo3d sig133'
c           call afic33(sig133)
            call b3_v33(sig133,sig13,vsig133)
c           construction matrice de passage inverse         
            call traps1(vsig133t,vsig133,3)
c           print*,'etape 2-0-1 matrice de passage ds fldo3d' 
c           call afic33(vsig133)             
      else
c           les directions principales ne changent pas lors du retour
            call chrep6(sig16,vsig133,.false.,sig16p)
            do j=1,3
                 sig13(j)=sig16p(j)
            end do
      end if            
 
c     ------------------------------------------------------------------          
 
c     passage de la loi de comportement dans la base principale
c     des contraintes
      if(.not.iso) then         
           print*,'programmer changement base loi de comp ds fldo3d'
           ierr1=1
           return
      end if 
 
c     passage des resistances et contraintes de refermeture
c     dans la base principale des contraintes
      if(.not.iso) then
          print*,'programmer changement base resistances ds fldo3d'
          read*
          ierr1=1
      end if 
 
c     passage dans la base principale du tir visco elastique
      call chrep6(epspc6,vsig133,.false.,epspc6p)        
      call chrep6(epspt6,vsig133,.false.,epspt6p)      
      call chrep6(epspg6,vsig133,.false.,epspg6p)
      call chrep6(epsps6,vsig133,.false.,epsps6p) 
 
c     conditions de positivite des valeurs normales  pour les fissures
      do j=1,3
          epspg6p(j)=max(epspg6p(j),0.d0)
          epsps6p(j)=max(epsps6p(j),0.d0)  
          epspt6p(j)=max(epspt6p(j),0.d0)          
      end do       
 
c     ****** traitement anti-battement *********************************
      if (ncal.ge.3) then
c           stockage des criteres actifs precedents        
            do i=1,nc
                 actif0(i)=actif1(i)
            end do 
      end if 
      if (ncal.ge.2) then
c           stockage des criteres actifs precedents        
            do i=1,nc
                 actif1(i)=actif(i)
            end do 
      end if               
 
c     ****** evaluation des criteres de plasticite *********************
      call crir3d(bgel,pg1,bdef,ps1,bw1,pw1,rteff,cwrt,kwrt,reft,
     # rceff,kwrc,delta,beta,rteff_app,rceff_app,actif,factif,nc,err1,
     # na,sig13,cgel,hplg,cdef,hpls,cshr,hplw,epspg6p,epsps6p,epspt6p,
     # epspc6p,alphadp,betadp,refermtl,taudp,taudpmax,complet,CPTCR)
 
c     traitement de l erreur eventuelle     
      if(err1.eq.1) then
             print*,'Pb crir3d dans Fluendo3d'
             ierr1=1
             return
      end if
c      print*,'Etape 9 ds fldo3d apres crir3d'         
 
c     ****** traitement anti-battement *********************************
      battcr3d=.false. 
      if(ncal.ge.3) then
            do i=1,nc
                 if(actif(i).and.actif0(i)) then
c                    print*,'Fluendo3d Battement critere',i ,'iter',ncal
                     battcr3d=.true. 
c                    on active la compatibilité de ce critere avec tous
c                    les autres pour la suite des calculs
                     do j=1,nc
                        CPTCR(i,j)=.true.
                        CPTCR(j,i)=.true.
                     end do                        
                 end if                     
            end do           
      end if    
 
      if(na.gt.0) then
C         print*,'Etat des criteres dans fldo3d: na=',na
C         do i=1,nc
C            print*,'actif(',i,'):',actif(i),'factif',factif(i)
C         end do 
          if(actif(13)) then
             log_plastc=.true.
          end if
      end if
c     read*             
 
c     ****** ecoulement plastique le cas échéant ***********************        
      if(na.gt.0) then
 
c         *calcul des increments dans la base principale des contraintes
c         *re-assemblage des matrices de fluage avec increment 
c         deformations totales nulles et initialisation de depsp
 
c         on résoud le retour radial sans reconsiderer les termes 
c         sources car ils ont deja été consideres lors du tir
          dphigel=0.d0
          dphidef=0.d0
          dphiw=0.d0
 
c         reevaluation des coeff de consolidation dans la direction
c         principale des contraintes effectives avec les deformations
c         initiales (pour conserver les CC du tir) cf.
 
          call conso3d(xkm,epsm,epser,epsm06,epse_tild6,cc03,
     #    vsig133,CMp00,CTHp,CTHv,Cwp)
c          print*,'fluendo3 etape2 cmp=',cmp00
 
c         assemblage du systeme et resolution sans increment de deformations (zero3)
c         et sans la source de fluage due aux def elastiques initiales (zero6)            
          call coupl3d(ann,xn,bn,ngf,err1,zero3,vsig133,vsig133t,dt,
     #    cc03,taum,taumdtt,tauk,psik,zero6,zero6,raideur66,Mgel,bgel,
     #    dphigel,Cgel,Hplg,Mdef,bdef,dphidef,Cdef,Hpls,Mwat,bw,dphiw,
     #    cshr,alphadp,betadp,actif,factif,nc,dsigef6,depse6,depsk6,
     #    depsm6,depspc6,depspt6,depsvt6,depspg6,depsps6,.false.,
     #    ipzero,dpg,dps,dpw,epspt6p,refermtl,complet,na,CMp00)
 
c         affichage def plastique de Drucker Prager
C         do i=1,6
C            print*,'depspc6(',i,')=',depspc6(i)
C         end do
 
c         traitement de l erreur eventuelle         
          if(err1.eq.1) then
             ierr1=1
             print*,'Pb lors du retour radial ds Fluendo3d'
             print*,'Etat des criteres lors du pb:'
             do i=1,nc
                print*,'critere:',actif(i),'valeur',factif(i)
             end do
             print*,'contrainte effective matrice'
             do i=1,3
                print*,'sig(',i,')=',sig13(i)
             end do
             return
          end if 
c          print*,'Etape 10-1 ds fldo3d apres coupl3d iter',ncal
 
      else 
 
c         pas d ecoulement plastique
          do i=1,6
             depse6(i)=0.d0          
             depsk6(i)=0.d0
             depsm6(i)=0.d0
             depsvt6(i)=0.d0 
             depspg6(i)=0.d0
             depsps6(i)=0.d0             
             depspc6(i)=0.d0
             depspt6(i)=0.d0
             dsigef6(i)=0.d0             
          end do
          dpg=0.d0
          dps=0.d0
          dpw=0.d0
c          print*,'Etape 10-2 ds fldo3d apres coupl3d iter',ncal          
 
      end if
 
 
c     ******************************************************************          
 
c     ***** actualisation des deformations (base fixe) *****************
 
      do i=1,6
c           elastique            
            epse16(i)=epse16(i)+depse6(i)
            varf(i+1)=epse16(i)        
c           kelvin         
            epsk16(i)=epsk16(i)+depsk6(i)
            varf(i+7)=epsk16(i)
c           stockage depsk pour mise au point         
            varf(122+i)=varf(i+7)-var0(i+7)  
            if((dt.eq.0.).and.(abs(varf(122+i)).gt.1.0d-10)) then
                 print*,'Pb dans fldo3d Etape 10-3 avec dt=0' 
                 print*,'depsk(',i,')=',varf(122+i)
                 print*,'devrait etre nul !'
                 ierr1=1
                 return
            end if                 
c           maxwell          
            epsm16(i)=epsm16(i)+depsm6(i)
            varf(i+13)=epsm16(i)
c           cas des deformations plastiques
c           cisaillement et dilatance
            epspc6(i)=epspc6(i)+depspc6(i)            
            varf(19+i)=epspc6(i)
c           print*,'actualisation varf epspc',i,varf(19+i)
c           traction base fixe
            epspt6(i)=epspt6(i)+depspt6(i)          
            varf(25+i)=epspt6(i)
c           dtt
            epsvt6(i)=epsvt6(i)+depsvt6(i)
            varf(31+i)=epsvt6(i)           
c           rag
            epspg6(i)=epspg6(i)+depspg6(i)
            varf(37+i)=epspg6(i)
c           rsi
            epsps6(i)=epsps6(i)+depsps6(i)
            varf(43+i)=epsps6(i)            
      end do
 
c     ****** actualisation deformation equivalente de compression ******         
      depleqc=0.d0
      do i=1,6
            depleqc=depleqc+(epspc6(i)-epspc60(i))**2  
      end do
      depleqc=sqrt(depleqc*2.d0/3.d0)         
c      print*,'depleqc',depleqc
c      actualisation et stockage  
      epleqc1=max(epleqc0+depleqc,epleqc0)       
      varf(94)=epleqc1
 
 
c     ****** actualisation des pressions *******************************
c     gel rag
      pg1=pg1+dpg
      varf(83)=pg1 
c     stockage pression maxi du gel
      varf(140)=max(var0(140),pg1)
c     def        
      ps1=ps1+dps
      varf(85)=ps1        
c     depression capillaire        
      pw1=pw1+dpw      
      varf(81)=pw1         
 
c     **** etat des contraintes effectives apres tir visco elastique 
      do i=1,6
         sig16(i)=sig16(i)+dsigef6(i)
         do j=1,6
             sigke16(i)=sigke16(i)+raideur66(i,j)*depsk6(j)*psik
         end do
c        stockage dans variables internes             
         varf(49+i)=sig16(i)
         varf(55+i)=sigke16(i)
C         print*,'apres test criter iter',ncal
C         print*,'sig16(',i,')=',sig16(i)
C         print*,'sigke16(',i,')=',sigke16(i)
      end do
 
c     ******************************************************************
c     sous iteration en cas de plasticite 
      if(na.ne.0) then
c        on re-teste les criteres apres cette iteration
         if(ncal.le.nmax) then
           if (ncal.ge.nmax1) then
c            on commence a afficher les pb eventuels                  
             print*,'Iter dans Fluendo3d:',ncal
c            print*,'Avant ecoulement :', ncal+1                   
             do i=1,nc
               if(actif(i)) then
                 print*,'critere',i,actif(i),factif(i)
                 if((i.ge.10).and.(i.le.12)) then
                    do j=1,3
                        print*,'defpt',j,epspt6p(j)
                    end do
                 end if
               end if
             end do 
             goto 30  
             do i=1,6
               print*,'sig(',i,')=',sig16(i)
             end do                       
             print*, 'Apres ecoulement:' , ncal                   
             do i=1,6
                print*,'epspt6(',i,')=',epspt6(i)
             end do
             do i=1,6
                print*,'epspc6(',i,')=',epspc6(i)
             end do       
             do i=1,6
                print*,'epsvt6(',i,')=',epsvt6(i)
             end do     
             do i=1,6
                print*,'epspg6(',i,')=',epspg6(i)
             end do  
             do i=1,6
                print*,'epsps6(',i,')=',epsps6(i)
             end do                       
             print*,'----------------------------'
c            read*   
30           continue               
           end if
c          on va faire une sous iteration plastique  
           goto 20
         else
           print*,'Approximation plasticite ds Fluendo3d'
           print*,'apres',ncal,' iterations'
           do i=1,nc
             print*,'critere',actif(i),factif(i)
           end do
           print*,'Nbr maxi de sous iter atteinte dans Fluendo3d'
           print*,'Augmenter NMAX ds Fluendo3d ou modifier le Pb'
           ierr1=1
           return
         end if
      end if
 
c***********************************************************************
c       fin des procedures visco-elasto-plastiques
c       print*,'Etape 11 ds fldo3d apres convergence',ncal
c***********************************************************************
C       print*,'dt ds fldo3d=',dt
C       print*,'Fluendo3d etape 10-2 a la convergence locale'             
C       do i=1,6
C        print*,'sig16(',i,')=',sig16(i)
C        print*,'sigke16(',i,')=',sigke16(i)
C        print*,'deps6(',i,')=',deps6(i)
C        print*,'depse6(',i,')=',depse6(i)
C        print*,'depsk6(',i,')=',depsk6(i)
C        print*,'depsm6(',i,')=',depsm6(i) 
C        print*,'depsvt6(',i,')=',depsvt6(i)
C        print*,'depspc6(',i,')=',depspc6(i)
C        print*,'depspt6(',i,')=',depspt6(i)
C        print*,'depspg6(',i,')=',depspg6(i) 
C        print*,'depsps6(',i,')=',depsps6(i)
C       end do 
C       read*
 
c***********************************************************************
c     actualisation des taux de franchissement par depression capillaire
 
c     passage en matrice 3x3
      call x6x33(fshr06,fshr33)
c     passage en base principale des contraintes
      call chre3(fshr33p,fshr33,vsig133)
c     actualisation des taux de passage des seuils en contrainte
      do i=1,3
           fshr33p(i,i)=max(fshr33p(i,i),factif(6+i),0.d0)
      end do
c     retour en base fixe et stockage        
      call chre3(fshr33,fshr33p,vsig133t)
c     retour en vecteur 6
      call x33x6(fshr33,fshr6)
c     stockage        
      do i=1,6
           varf(129+i)=fshr6(i)
      end do             
 
c***********************************************************************
c     actualisation des ouvertures de fissures actuelles et maxi de traction
 
c     recupereration des ouvertures initiales
      do i=1,6
c        ouverture actuelle
         wplt06(i)=var0(96+i)
c         print*,'av ma:',wplt06(i),var0(96+i),96+i
c        ouverture maxi sur le temps
         wpltx06(i)=var0(102+i)
c         print*,'av ma:',wpltx06(i),var0(102+i),102+i
      end do
c     mise a jour des ouvertures des 3 fissures localisees
      if(ENDO_NL) then
 
c        ************ deformation plastique non locale ***************** 
         alpha_nl=1.d0
c        rangement des déformation plastiques non locales en tableau 3*3
         call x6x33(eps_nl6,eps_nl33)
c        diagonalisation tenseur def non locales
         call b3_v33(eps_nl33,eps_nl3,veps_nl33)
c        creation de la matrice de passage inverse    
         call traps1(veps_nl33t,veps_nl33,3)
c        Passage de la déformation locale dans la base principale de la déformation non locale
c        Rangement des déformation plastiques locales en tableau 3*3
         call x6x33(epspt6,epspt33)
         call chre3(epspt33p,epspt33,veps_nl33)
c                Indicateur et deformation non-locale equivalente
c        Chargement du epsilon non-local equivalent initial
         do i=1,6
            eptmasq06(i)=var0(146+i)
         end do
         call x6x33(eptmasq06,eptmasq033)
         do i=1,3  
            do j=1,3
                eptmasq33p(i,j)=0.d0  
            end do
            if((epspt33p(i,i).gt.alpha_nl*eps_nl3(i)).and.(epspt33p(i,i)
     #          .gt.0.d0).and.(eps_nl3(i).gt.0.d0))then    
                eptmasq33p(i,i)=eps_nl3(i)
            end if
         end do
c        retour en base fixe  
         call chre3(eptmasq33,eptmasq33p,veps_nl33t)
         call x33x6(eptmasq33,eptmasq6)
         do i=1,6
            varf(146+i)=eptmasq6(i)
         end do
c        ************ fin deformation plastique non locale *************
 
c        print*,'fluendo endo_nl av majw3d', Method_N,LcH,Method_H
c        l endommagement et les ouvertures sont calculee avec la 
c        deformation plastique non locales         
         call majw3d(eptmasq06,eptmasq6,wplt06,wplt6,wpltx06,wpltx6,
     #   wpl3,vwpl33,vwpl33t,wplx3,vwplx33,vwplx33t,XE3D,NBNMAX3D,
     #   NBNB3D,IDIMB3D,Method_N,LcH,Method_H,err1)
c        print*,'fluendo endo_nl ap majw3d', Method_N,LcH,Method_H      
      else
         call majw3d(epspt60,epspt6,wplt06,wplt6,wpltx06,wpltx6,
     #   wpl3,vwpl33,vwpl33t,wplx3,vwplx33,vwplx33t,XE3D,NBNMAX3D,
     #   NBNB3D,IDIMB3D,Method_N,LcH,Method_H,err1)
      end if
      if(err1.eq.1)then
         print*,'Pb lors du calcul ds majw3d dans fldo3d'
         ierr1=1
         return
      end if 
      do i=1,6
         varf(96+i)=wplt6(i)
         varf(102+i)=wpltx6(i)
c        print*,'fldo3d ap maj finale',wpltx6(i)
      end do               
c     ouverture maxi actuelle
      varf(109)=max(wpl3(1),wpl3(2),wpl3(3))
c     indicateur de progression de la fissure
      if(abs(varf(109)-var0(109)).gt.coeff_petit) then
          travail_fissure=.true.
      else
          travail_fissure=.false.
      end if            
C       print*,'Etape 12 ds fldo3d apres majw3d',ncal 
C       read*       
 
c***********************************************************************
c     calcul des contraintes dans les fibres ( Romain Gontero)
c      include './contraintes_dans_les_fibres.h'
-INC HSIGFIB
 
c***********************************************************************
c     contraintes dans solide et rgi en fin de pas avec
c     prise en compte de l endo thermique et de fluage
      umdt=(1.d0-dth1)*(1.d0-dflu1)           
c     resultante des pressions intraporeuses RGI et 
c     Capillaire (depression)
      sigp=-bw1*pw1-bg1*pg1-bs1*ps1 
c     print*,bw1,pw1,bg1,pg1,bs1,ps1
c     print*,'fluendo sigp',sigp
c     effet sur la contrainte apparente en non sature        
      do i=1,6
            if(i.le.3) then
c               prise en compte de la pression rgi           
                sigf6(i)=(sig16(i)+sigp)*umdt
            else
                sigf6(i)=sig16(i)*umdt
            end if
c           print*,'sigf6',sigf6(i),' sigp',sigp,'umdt',umdt
      end do
 
c***********************************************************************
c     prise en compte de l'endommagement mécanique 
 
      if (end3d.and.(.not.plast_seule)) then
c        chargement de la variable de controle d erreur de gf
         errgf0=var0(110)
c        endo pre pic traction         
         dtr=dtr0
c        endo pre pic de compression
         dcpp=dcpp0         
c        actualisation de gft en fonction de la depression capillaire
         gft_app=gft*rt_app/rt  
c        print*,'fldo3d Method_H', Method_H
c        read*       
c        calcul des contraintes endommagee et des endommagements         
         call kend3d(wpl3,vwpl33,vwpl33t,wplx3,vwplx33,vwplx33t,
     #   gft_app,gfr,iso,sigf6,sigf6d,err1,reft,young,nu,souplesse66,
     #   xmt,dtiso,dtr,dt3,dr3,dct3,0.,dcc3,nc3,nc33,ept,errgf1,
     #   errgf0,epspc6,0.d0,ekdc,epspg6,ekdg,alphag,dgt3,dgc3,fshr6,
     #   dwt3,dwc3,skdw,0.d0,epsps6,dst3,dsc3,ekds,alphas,rteff_app,
     #   ann,bn,xn,ipzero,ngf,tau_trac,xmc,dciso,dcpp,taudp,taudpmax,
     #   XE3D,NBNMAX3D,NBNB3D,IDIMB3D,altc,epssr,Method_H,Method_N,LcH)
c        stockage variable de controle de l erreur sur Gft     
         varf(110)=errgf1
c        endo pre pic         
         varf(79)=dtr
c        endo pre pic de compression - cisaillement
         varf(122)=dcpp         
c        traitement des erreurs         
         if(err1.eq.1) then
             print*,'pb calcul des endommagements dans fluendo3d'
             ierr1=1
             return
         end if
c        contraintes totale dans la matrice apres endommagement 
         do i=1,6
             sigmf6(i)=sigf6d(i)
             varf(61+i)=sigmf6(i)
         end do        
       else
c        pas d endommagement
         dcpp=0.d0
         dtr=0.d0
         varf(79)=dtr
         varf(122)=dcpp
         do i=1,3
             dt3(i)=0.d0
             nc3(i)=1.d0
             dr3(i)=0.d0
             dwt3(i)=0.d0
             dwc3(i)=0.d0
             dgt3(i)=0.d0
             dgc3(i)=0.d0
             dst3(i)=0.d0
             dsc3(i)=0.d0
             dcc3(i)=0.d0
             dct3(i)=0.d0
         end do 
         do i=1,6
             sigmf6(i)=sigf6(i)
             varf(61+i)=sigmf6(i)
         end do
      end if  
 
c***********************************************************************
c     combinaison fibres matrice endommagee (Romain Gontero 2022)
c      include './combinaison_matrice_fibres.h'
-INC HCOMFIB
 
c***********************************************************************
c     contrainte dans la matrice (eventuellement fibree)
      do i=1,6
c        debut romain 
         if(ifibre.eq.1) then            
            sigmf6(i)=(1.d0-rhof)*sigmf6(i)+sfib(i)
         end if 
c        print*,"combinaison contrainte fibres"             
c        fin romain
         varf(61+i)=sigmf6(i)
      end do
 
c*********** variables initiales pour les renforts *********************
c     lignes a copier si le modele utilise des renforts 
 
      if (NRENF00.ne.0) then
         do i=1,NRENF00
c            include './lire_vari_renforts.h'
-INC HLVIREN
c           recuperation des non locales a utiliser pour les renforts
c           pour les renforts le non local est fait sur SER(i) de numero 
c           NVAR2+(i-1)*NB_VARI_PAR_RENF+26
c           on teste si ce renfort est traite par helmholtz
            if(log_H_RENF(i)) then 
c             recuperation des variables pour une traitement non local du renfort          
              ivarnl=NVAR2+(i-1)*NB_VARI_PAR_RENF+26 
c             numero de la formulation non locale
              nl=Num_H_RENF(i)
c             position des variables de Helmholtz associees a cette formulation
              if(istep.ne.0) then                
c                recup des non locales associees aux renforts
                 if(istep.ne.1) then   
c                    recuperation de la contrainte non locale intermediaire
                     sigr_nl(i)=var0(ivarnl)    
c                    recuperation du gradient de la contrainte sur Hi2
                     grad_sigr(i)=Helm0(nl,2)                  
c                    survie du renfort
                     coeff_nlr(i)=Helm0(nl,3) 
c                    survie de l interface                     
                     coeff_nli(i)=Helm0(nl,4)
c                    coeff de taux dans l erreur
                     Helm1(nl,6)=Helm0(nl,6)
c                    contrainte non locale de l iteration non locale precedente
                     sigr_nl0(i)=Helm0(nl,7)                     
c                    recuperation de la variable de controle des hypotheses
                     Pb_nl(i)=Helm0(nl,8)
c                    hypothese ecoulement local
                     hypo_nl(i)=Helm0(nl,9) 
c                    copie du coeff de normalisation de taux 
                     helm1(nl,10)=helm0(nl,10)               
                 else
c                    cas istep=1, on a pas encore fait d etape non locale 
c                    contrainte non locale de l iteration non locale precedente
                     sigr_nl0(i)=Helm0(nl,1)                 
c                    recuperation de la contrainte non locale convergee precedente                 
                     sigr_nl(i)=sigr_nl0(i)                         
c                    recuperation du gradient de la contrainte non locale convergee precedente sur YU2i
                     grad_sigr(i)=Helm0(nl,2)                   
c                    on recharge les hypotheses sur endo et refermeture
                     if(ppas) then
c                        survie du renfort
                         coeff_nlr(i)=1.0d0
c                        survie de l interface 
                         coeff_nli(i)=1.0d0
                     else
c                        survie du renfort
                         coeff_nlr(i)=Helm0(nl,3)
c                        survie de l interface 
                         coeff_nli(i)=Helm0(nl,4)                          
                     end if
c                    hypothese pas d ecoulement plastique
                     Pb_nl(i)=0.d0 
                     hypo_nl(i)=0.d0
                 end if   
              else
c                traitement local du renfort
                 grad_sigr(i)=0.d0               
                 coeff_nlr(i)=1.d0
                 coeff_nli(i)=1.d0
                 Pb_nl(i)=0.d0 
                 hypo_nl(i)=0.d0                   
              end if
            else
c                traitement local du renfort
                 grad_sigr(i)=0.d0
                 coeff_nlr(i)=1.d0
                 coeff_nli(i)=1.d0
                 Pb_nl(i)=0.d0 
                 hypo_nl(i)=0.d0                    
           end if
         end do
      end if
 
c     ****** variables initiales pour WL2 ******************************
 
      if (log_H_TWL2) then 
c        resultat non local      
         xwl2=var0(71)
c        max de alpha  
         xwl2max=var0(73)
c        max du taux de charge          
         tau_trac_wl2_max=var0(74)
c        endo localise          
         dtl00=var0(75)         
c        coeff d effet d echelle  precedent        
         cwrt=var0(76)
      end if
 
c************* Etapes non locales intermediaires ***********************
 
      if((istep.ne.0).and.(istep.ne.2)) then 
c         *** reinitialisation des variables locales pour le calcul reel
c         car les sous iterations non locales ne doivent pas conserver 
c         les variables locales non convergees lors de la copie de varf 
c         sur var0 lors du passage dans UNPAS en non local sauf si 
c         istep est egal a 2 (il  s agit alors du dernier passage non local)
          do i=1,NVARI
             varf(i)=var0(i)
          end do 
 
c         **** etape non locale intermediaire pour WL2 *****************
 
c         'TWL2' est la variable a moyenner pour la methode wl2          
          ivarnl=ivari_twl2
          call exhmtz(istep,NBVIA3D,INLVIA3D,NB_HELM,
     #    ivarnl,log_H_TWL2,Num_H_TWL2)
c         numero de la formulation associee
          NL=Num_H_TWL2
 
          if (log_H_TWL2) then
c            ** le calcul de xwl2 est a faire **************************                   
c            print*,'fldo3d, istep=1 avec cwrt=',cwrt
c            recuperation des contraintes dans la matrice          
             call x6x33(sigmf6,sig133)         
             call b3_v33(sig133,sig13,vsig133)
c            prise en compte de la dilution d endommagement lors de l hydratation
             call hydrvari3d(dtl00,dtl0,hyd0,hydr,hyds) 
c            recuperation de l indicateur de localisation fin du nouveau pas
             dtl1=max(dt3(1),dt3(2),dt3(3))           
c            taux de charge en traction actualise pour la nouvelle deformation         
             tau_trac_wl2=tau_trac
c            pression capillaire
             pw0=pw1
c            saturation en eau * biot            
             bw0=bw1        
c            1er passage dans wl2          
             call wl2d3d(1,vmax,vref,sig13,eweb,rt_app,cwrt,xwl2,
     #       tau_trac_wl2,dtl0,dtl1,xwl2max,tau_trac_wl2_max,bw0,
     #       pw0,kwrt)
c            *** on actualise des vari a passer dans non local pour Weibull
c            les variables de Weibull (alpha(xwl2),smax(tau_trac_wl2)) 
c            alpha
             varf(ivarinl)=xwl2 
c            sigma/rt maxi pour max(structure)   
             print*,'a refaire avec les helm1...'
             ierr1=1
             return      
             varf(72)=tau_trac_wl2
c            max xwl2 un pour multiplier par xwl2 non local maxi 
c            ds la procedure non locale          
             varf(73)=1.d0 
c            max twl2 un pour multiplier par twl2 non local maxi
c            ds la procedure non locale          
             varf(74)=1.d0   
c            reconduction de l endommagement localise debut de ne 
c            pas affecte par hydratation
             varf(75)=dtl0
c            reconduction de cwrt en cas de localisation
             varf(76)=cwrt 
c            pression capillaire fin de pas pour rt istep=2
             varf(81)=pw1
c            saturation en eau * biot fin de pas pour istep=2           
             varf(80)=bw1 
c            endommagement fin de pas
             varf(119)=dtl1     
          else
c            ** rt00 n est pas affectee par la methode WL2 *************           
c            pas de non local pour xwl2       
             cwrt=1.d0
c            liste des variables non locales : alpha
             varf(71)=0.d0
c            variable non locale : sigma maxi /rt           
             varf(72)=0.d0  
c            max de xwl2 non utilisse on y laisse 1.
             varf(73)=1.d0
c            contrainte relative maxi non locale non utiliseee
             varf(74)=0.d0 
c            pas d endo de traction utilisable pour xwl2
             varf(75)=0.d0             
c            modif de rt par effet d echelle inactif            
             varf(76)=cwrt          
          end if
 
c         *** etape non locale intermediaire pour les renforts *********
 
          if (NRENF00.ne.0) then
 
c            intialisation des tenseurs pour les renforts
             call prrf3d(dalpr,teta1,teta2,tetaref,epst06,epst033,
     #       epstf6,epstf33,epspmf6,epspt6,epspmf33)
 
 
c            calcul de la source non locale pour chaque renfort             
             do i=1,NRENF00
 
c             cas des renforts traites par helmholtz
              if (log_H_RENF(i)) then
c               numero de la formulation d Helmholtz pour le renfort i 
                nl=Num_H_RENF(i)
c               recuperation des non locales a utiliser pour les renforts
c               pour les renforts le non local est fait sur SERi de numero 
                ivarnl=NVAR2+(i-1)*NB_VARI_PAR_RENF+26
                if(rhor(i).gt.0.) then                                         
c                etape non locale de rnfr3d on calcule la source 
c                non locale pour le pas suivant           
                 call rnfr3d(istep,NB_RENF,i,epstf33,vecr,epsr0(i),
     #           epsrf(i),eplr0(i),eplrf(i),yor(i),syr(i),sigre0(i),
     #           sigref(i),hplr(i),tor(i),ekr(i),skr(i),ATRR(i),khir(i),
     #           gamr(i),sprec(i),teta1,teta2,dt,ppas,theta,eprm0(i),
     #           eprmf(i),ttaref(i),rhor(i),mu_r0(i),fl3d,errr,xnr(i),
     #           xmuthr(i),eprk0(i),eprkf(i),tokr(i),yksyr(i),
     #           plast_seule,ann,xn,bn,ngf,ipzero,epspmf33,spre0(i),
     #           spref(i),epst033,depsa(i),epleq0(i),epleqf(i),
     #           epsmf(i),epsm0(i),sigr_nl(i),coeff_nlr(i),eper0(i),
     #           epart(i),eparv(i),hypo_nl(i),plr_iso,log_H_RENF(i),
     #           garder,Pb_nl(i),travail_fissure)   
c                traitement de l erreur     
                 if(errr) then
                     print*,'erreur dans rnfr3d'
                     ierr1=1
                     return
                 end if
c                mise a jour de la def anelastique pour la source de
c                l iteration suivante             
                 if(garder) then
c                   deformation visqueuse 
                    eparv(i)=eprkf(i)+eprmf(i)
c                   deformation anelastique totale 
                    epart(i)=eparv(i)+eplrf(i)                 
                 else                 
c                   deformation visqueuse initiale
                    eparv(i)=eprk0(i)+eprm0(i)
c                   deformation anelastique totale initiale
                    epart(i)=eparv(i)+eplr0(i)
                 end if                 
c                stockage de la longueur de l element fini                
                 varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+13)=longr(i)                                                       
c                calcul des endommagements  et de la correction de source               
                 call endorf(sigref(i),eplrf(i),syr(i),sur(i),
     #           epsrf(i),epu(i),hplr(i),yor(i),wpr(i),longr(i),
     #           damrt0(i),damrtf(i),damrt1(i),damrt2(i),damrc0(i),
     #           damrcf(i),refrf(i),eplrt(i),eplrc(i),plr_iso,
     #           diffr0(i),diffrf(i),yor(i),deqr(i),Hintr(i),
     #           grad_sigr(i),tauie(i),taur(i),dami0(i),damif(i),
     #           epsmf(i),eplr0(i),sigre0(i),log_H_RENF(i),epsm0(i),
     #           epsr0(i),depsa(i),Et0(i),Etf(i),Hs0(i),Hsf(i),sigrf(i),
     #           spref(i),epart(i),eparv(i),source_nl(i),coeff_nlr(i),
     #           taui0(i),tauif(i),coeff_nli(i),eprkf(i),
     #           eprmf(i),hypo_nl(i),istep,ierr_renfort,garder,
     #           endo_interface_explicit)       
c                traitement erreur endo renfort     
c                 include './traitement_erreur_renforts.h'
-INC HERRREN                                    
c                sauvegarde de la contrainte  a moyenner 
                 varf(ivarnl)=source_nl(i)                 
c                survie du renfort : coeff_nlr(i)
                 Helm1(nl,3)=coeff_nlr(i)
c                survie de l interface  :   coeff_nli(i)  
                 Helm1(nl,4)=coeff_nli(i)
c                hypotheses d ecoulements realises :  hypo_nl(i)
                 Helm1(nl,9)= hypo_nl(i)   
c                endommagement de l interface pour le calcul implicite du gradient endommagement d interface :   damif(i) 
                 varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+20)=damif(i)
c                contrainte de cisaillement implicite
                 varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+19)=tauiF(i)                                              
c                sauvegarde de la contrainte
                 if(garder) then
c                 les hypotheses sont coherentes on garde la solution                 
c                 valeur de la contrainte non locale avant ecoulement (SERi)
                  Helm1(nl,1)=sigr_nl(i)
c                 valeur de la contrainte endommagee  (SNRi)                
                  varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+14)= sigrf(i)               
c                 actualisation de la diffusion non locale (DFRi)
                  varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+18)=diffrf(i)
c                 stockage du module tangent du renfort ETRi
                  varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+21)=Etf(i)
c                 stockage module secant interface  HSRi               
                  varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+22)=Hsf(i)
c                 deformation anelastique EARi
                  varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+24)=epart(i)
c                 stockage de la deformation visqueuse du renfort EVRi
                  varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+25)=eparv(i)
                 else
c                 les hypotheses etaient incoherentes, on reprend la 
c                 valeur de la contrainte  locale avant l etape non locale
                  if(istep.ne.1) then
c                    les hypotheses etaient incoherente on reprend la 
c                    valeur de la contrainte  locale avant l etape non locale                  
                     Helm1(nl,1)=sigr_nl0(i)
                  else
c                    comme on a pas encore fait d etape non locale
c                    on prend la contrainte calculée localement                           
                     Helm1(nl,1)=sigrf(i)
                  endif
                 end if
                 if(endo_interface_explicit.and.(istep.eq.1)) then
                   varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+20)=dami0(i)
                   varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+19)=taui0(i)                   
                 end if              
                else
c                pas de renfort num i
c                rho(i)=0               
                 do j=1,9
c                 mise a zero des vari du renfort
                  varf(ivarnl)=0.d0
c                 valeur de la contrainte apres ecoulement 
                  sigrf(i)=0.d0                
                  varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+14)=sigrf(i)               
c                 actualisation de la diffusion non locale
                  diffrf(i)=0.d0 
                  varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+18)=diffrf(i)
c                 stockage du module du renfort
                  Etf(i)=0.d0
                  varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+21)=Etf(i)
c                 stockage module interface  
                  Hsf(i)=0.d0
                  varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+22)=Hsf(i)
c                 stockage de l hypothese de deformation anelastique
                  epart(i)=0.d0
                  varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+24)=epart(i)
c                 stockage de l autre hypothese
                  eparv(i)=0.d0
                  varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+25)=eparv(i)                   
                 end do
                end if
              end if                 
            end do              
          end if             
      end if
c******** fin l etape non locale intermidiaire *************************
 
 
c******** Derniere etape non locale ou calcul local ********************
 
c     *****  calcul des contraintes dans les renforts ******************
 
      if ((NRENF00.ne.0).and.((istep.eq.0).or.(istep.eq.2))) then
c          calcul du taux volumique d armatures
           rhov=0.d0
           do j=1,NRENF00
                rhov=rhov+rhor(j)
           end do
           if (rhov.gt.1.d0) then
               print*,'Pb car taux d armature > 1 dans fluendo3d ! :('
               ierr1=1
               return
           end if
c          initialisation des tenseurs de contraintes homogeneisee        
           do j=1,3
                do k=1,3
                    sigrm33(j,k)=0.d0
                    sigrf33(j,k)=0.d0
                end do
           end do 
c          intialisation des tenseurs pour les renforts
           call prrf3d(dalpr,teta1,teta2,tetaref,epst06,epst033,
     #     epstf6,epstf33,epspmf6,epspt6,epspmf33)
 
c          calcul des contraintes axiales dans chaque renfort             
           do i=1,NRENF00
c             recuperation des non locales a utiliser pour les renforts
c             pour les renforts le non local est fait sur SER de nuemero 
              ivarnl=NVAR2+(i-1)*NB_VARI_PAR_RENF+26
              if(rhor(i).gt.0.) then                
c                recuperation du tir elastique non local                    
                 if((istep.eq.2).and.log_H_RENF(i)) then
c                    derniere estimation pour les renforts non locaux                       
cc                    sauvegarde du tir non local precedent car istep=2
cc                   avec la precedente contrainte non locale
c                    numero de la formulation non locale
                     nl=Num_H_RENF(i)
c                    contrainte non locale
                     sigr_nl(i)=var0(ivarnl) 
c                    on utilse helm1(nl,1) pour stocker le resultat non local
                     Helm1(nl,1)=sigr_nl(i)
c                    recuperation du gradient
                     grad_sigr(i)=Helm0(nl,2)
c                    sauvegarde du gradient de sigr pour le pas suivant
                     helm1(nl,2)=grad_sigr(i)                                                        
c                    copie de la valeur initiale en cas de reprise
                     helm1(nl,7)=sigr_nl(i)
c                    coeff du gradient pour l erreur
                     Helm1(nl,10)=0.25d0*(deqr(i)/sur(i)) 
c                    recuperation des predicteurs pour le pas suivant 
                     Etf(i)=Et0(i)
                     Hsf(i)=Hs0(i)                     
                 else if (istep.eq.0) then
                    continue
                 else if ((istep.eq.2).and.(.not.log_H_RENF(i))) then
c                   certains aciers sont calcules en local                 
                    continue
                 else
c                   recuperation des predicteurs pour le pas suivant 
                    print*,'erreur d aiguillage dans fldo3d'
                    print*,'istep',istep,'Helmholtz(',nl,')=',
     #              log_H_RENF(i)
                    ierr1=1
                    return
                 end if
c                contrainte effective axiale dans les renforts                
                 call rnfr3d(istep,NB_RENF,i,epstf33,vecr,epsr0(i),
     #           epsrf(i),eplr0(i),eplrf(i),yor(i),syr(i),sigre0(i),
     #           sigref(i),hplr(i),tor(i),ekr(i),skr(i),ATRR(i),khir(i),
     #           gamr(i),sprec(i),teta1,teta2,dt,ppas,theta,eprm0(i),
     #           eprmf(i),ttaref(i),rhor(i),mu_r0(i),fl3d,errr,xnr(i),
     #           xmuthr(i),eprk0(i),eprkf(i),tokr(i),yksyr(i),
     #           plast_seule,ann,xn,bn,ngf,ipzero,epspmf33,spre0(i),
     #           spref(i),epst033,depsa(i),epleq0(i),epleqf(i),
     #           epsmf(i),epsm0(i),sigr_nl(i),coeff_nlr(i),eper0(i),
     #           epart(i),eparv(i),hypo_nl(i),plr_iso,log_H_RENF(i),
     #           garder,Pb_nl(i),travail_fissure)       
c                traitement erreur rnfr3d     
                 if(errr) then
                     print*,'erreur dans rnfr3d'
                     ierr1=1
                     return
                 end if 
c                action suivant conformite hypothese                 
                 if(.not.garder) then
                      print*,'Pb dans fldo3d'
                      print*,'pour istep=2 on devrait garder'
                      print*,'la solution du sous pas non local'
                      ierr1=1
                      return
                 end if                
c                deformation visqueuse 
                 eparv(i)=eprkf(i)+eprmf(i)
c                deformation anelastique totale 
                 epart(i)=eparv(i)+eplrf(i)
c                nouvelle deformation elastique 
                 eperf(i)=epsrf(i)-epart(i)                
c                print*,'fldo3d:sigref(',i,')=',sigref(i)
c                calcul des endommagements et actualisation lnl, dsource
                 call endorf(sigref(i),eplrf(i),syr(i),sur(i),
     #           epsrf(i),epu(i),hplr(i),yor(i),wpr(i),longr(i),
     #           damrt0(i),damrtf(i),damrt1(i),damrt2(i),damrc0(i),
     #           damrcf(i),refrf(i),eplrt(i),eplrc(i),plr_iso,
     #           diffr0(i),diffrf(i),yor(i),deqr(i),Hintr(i),
     #           grad_sigr(i),tauie(i),taur(i),dami0(i),damif(i),
     #           epsmf(i),eplr0(i),sigre0(i),log_H_RENF(i),epsm0(i),
     #           epsr0(i),depsa(i),Et0(i),Etf(i),Hs0(i),Hsf(i),sigrf(i),
     #           spref(i),epart(i),eparv(i),source_nl(i),coeff_nlr(i),
     #           taui0(i),tauif(i),coeff_nli(i),eprkf(i),
     #           eprmf(i),hypo_nl(i),istep,err1,garder,
     #           endo_interface_explicit)
C                  if (sigrf(i).ne. sigref(i)) then
C                    print*,'Ds fluendo3d sigrf(i).ne. sigref(i)'
C                    print*,i,sigrf(i),sigref(i)
C                    ierr1=1
C                    return
C                  end if   
c                print*,'Dans fldo3d,sigrf(',i,')=',sigrf(i)   
c                traitement erreur endommagement     
                 if(err1.eq.1) then
                     ierr1=1
                     return
                 end if
                 if((istep.eq.2).and.log_H_RENF(i)) then
c                    rappel numero de la formulation non locale
                     nl=Num_H_RENF(i)
c                    contrainte non locale derniere iteration
                     Helm1(nl,1)=sigr_nl(i)
c                    gradient de la contrainte
                     Helm1(nl,2)=grad_sigr(i)
c                    survie du renfort
                     Helm1(nl,3)=coeff_nlr(i)
c                    survie de l interface                     
                     Helm1(nl,4)=coeff_nli(i) 
c                    coeff de taux dans l erreur
                     Helm1(nl,6)=1.d0/sur(i)
c                    valeur de la contrainte non locale utilisee
                     Helm1(nl,7)=sigr_nl(i)
c                    numero sous iteration non locale
                     Helm1(nl,8)=Pb_nl(i)
c                    moyenne des hypotheses d ecoulements realises                    
                     Helm1(nl,9)=hypo_nl(i)                                      
c                   print*,'fluendo istep',istep,'errsigr',Helm1(nl,4)                      
                 end if 
              else
                 epsrf(i)=0.d0
                 eperf(i)=0.d0             
                 eplrf(i)=0.d0
                 sigref(i)=0.d0
                 eprmf(i)=0.d0
                 mu_r0(i)=0.d0
                 eprkf(i)=0.d0
                 spref(i) =0.d0
                 eplrc(i) =0.d0
                 damrtf(i)=0.d0
                 damrt1(i)=0.d0
                 damrt2(i)=0.d0                 
                 sigrf(i)=0.d0
                 epleqf(i) =0.d0
                 damrcf(i) =0.d0
                 eplrt(i)=0.d0
                 diffrf(i)=0.d0
                 tauif(i)=0.d0
                 damif(i)=0.d0          
                 Etf(i)=0.d0
                 Hsf(i)=0.d0
                 refrf(i)=0.d0
                 epart(i)=0.d0
                 eparv(i)=0.d0          
              end if               
c             print*,'fldo3d Etape 14 istep2 eper:',eperf(i)                    
c             stockage des variables internes des renforts
c              include './ecrire_vari_renforts.h'  
-INC HEVIREN            
c             print*,'Etape 13-3 istep 2 ds fluendp3d'   
c             print*,'vari sigrf',i,varf(nvar2+(i-1)*NB_VARI_RENF+14),sigrf(i)  
 
c             **** tenseur des contraintes dues aux renforts ************ 
              do j=1,3
                    do k=1,3
                        if(rhor(i).gt.0.) then
                            sigrm33(j,k)=sigrm33(j,k)+rhor(i)*sigrf(i)*
     #                      vecr(i,j)*vecr(i,k)
                        end if
                    end do                
              end do
c             fin de la boucle sur le renfort i              
           end do 
c          mise a zero des vari des renforts non actifs            
           if (NRENF00.lt.NB_RENF) then
                do i=(NRENF00+1),NB_RENF
                    do j=1,NB_VARI_PAR_RENF             
                        varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+j)=0.d0
                    end do
                end do
           end if 
      else
c            pas de renforts ou etape non locale intermediaire
             rhov=0.d0
c            **** tenseur des contraintes dues aux renforts ************ 
             do j=1,3
                    do k=1,3
                            sigrm33(j,k)=0.d0
                    end do                
             end do
c            variables internes des renforts
             if(NRENF00.eq.0) then
c               mise a zero des vari des renforts             
                do i=1,NB_RENF  
                    do j=1,NB_VARI_PAR_RENF             
                        varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+j)=0.d0
                    end do
                end do
             end if               
      end if
c     combinaison des contraintes dans les renforts et la matrice
c     print*,'contrainte dans les aciers'
      call x33x6(sigrm33,sigrm6)  
c     call afic6(sigrm6)
c     print*,'contraintes dans la matrice '
c     call afic6(sigmf6)                
      do i=1,6
         sigf6d(i)=(1.d0-rhov)*sigmf6(i)+sigrm6(i)
      end do
c     print*,'contrainte moyenne BA'
c     call afic6(sigf6d)
 
 
c***********************************************************************
c     calcul des endommagements equivalents pour affichage
c***********************************************************************
 
      if(end3d.and.(.not.plast_seule)) then
c            endommagement equivalent de traction 
             dtl1=max(dt3(1),dt3(2),dt3(3))
             if(istep.eq.0) then
                varf(75)=var0(119)
                varf(119)=dtl1                
             end if
             dtw0=max(dwt3(1),dwt3(2),dwt3(3))
             dtg0=max(dgt3(1),dgt3(2),dgt3(3))
             dts0=max(dst3(1),dst3(2),dst3(3))             
             dtra=1.d0-(1.d0-dtr)*(1.d0-dflu1)*(1.d0-dth1)*
     #                            (1.d0-dtl1)*
     #                            (1.d0-dtw0)*    
     #                            (1.d0-dtg0)*
     #                            (1.d0-dts0)
c            endommagement equivalent de compression
             dcw0=max(dwc3(1),dwc3(2),dwc3(3))
             dcg0=max(dgc3(1),dgc3(2),dgc3(3))
             dcs0=max(dsc3(1),dsc3(2),dsc3(3))
             dct0=max(dct3(1),dct3(2),dct3(3))
             dcc0=max(dcc3(1),dcc3(2),dcc3(3))
             dcom=1.d0-(1.d0-dflu1)*(1.d0-dth1)*(1.d0-dcpp)*
     #                            (1.d0-dcw0)*     
     #                            (1.d0-dcg0)*
     #                            (1.d0-dcs0)*
     #                            (1.d0-dct0)*
     #                            (1.d0-dcc0)
      else
c           cas sans endommagement         
            dtra=0.d0
            dcom=0.d0
      end if 
      varf(95)=dtra
      varf(96)=dcom 
 
c     endo meca traction     
      varf(111)=1.d0-(1.d0-dtl1)*(1.d0-dtr)
c     endo hydrique traction
      varf(112)=dtw0
c     endo rag traction
      varf(113)=dtg0
c     end do def traction
      varf(114)=dts0
 
c     endo meca compression 
      varf(115)=1.d0-(1.d0-dct0)*(1.d0-dcc0)*(1-dcpp)
c     endo hydrique compression
      varf(116)=dcw0
c     endo rag compression
      varf(117)=dcg0
c     endo def compression
      varf(118)=dcs0
c     endo meca prepic
      varf(122)=dcpp 
c     endo meca compression par couplage avec traction (buckling)
      varf(136)=dct0
c     contrainte equivalente de Drucker Prager 
      varf(137)=taudp
c     endommagement de compression par cisaillement  
      varf(138)=dcc0         
 
c***********************************************************************        
c     affectation dans le tableau de sortie des contraintes
c     et prise  encompte des contribution des renforts
      do i=1,nstrs
             sigf(i)=sigf6d(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***********************************************************************
c     VARIABLES INTERNES A ACTUALISER POUR LES FORMULATIONS DE HELMHOLTZ
C***********************************************************************
c     sauvegarde des vari de helmholtz HELM1(NL,ii=1..NBR_VARI_HELM) 
      do NL=1,NB_HELM
c        include './ecrire_vari_helmholtz.h'
-INC HEVIHEL   
      end do
c     autre variables internes a sauver en cas d etape non locale
      if(ENDO_NL) then
c       residu de deformation plastique locale a diffuser
        if((istep.eq.1).or.(istep.eq.3)) then
            do i=1,6
c               source locale
                varf(140+i)=epspt6(i)
c               -epspt60(i)
            end do
        else
            do i=1,6
c               stockage valeur non locale 
                varf(140+i)=eps_nl6(i)
            end do          
        end if
      end if
 
C***********************************************************************
c     stockage du pas de temps pour verif
      varf(139)=dt   
c***********************************************************************        
c     traitement erreur residuelle
      if(err1.eq.1) then        
             ierr1=1
             return 
      end if 
 
      return
c     fin         
c*********************************************************************** 
 
 
 
 
c***********************************************************************       
c     affichages pour etude des cas ou dt = 0 (dt est stocke dans vari
c     (139) si dt = 0, dans castem les vari non convergees de la 
c     dernieresi dt = 0, dans castem les vari non convergees de 
c     la derniere sous iteration avant la convergence forcee conservee        
      do i=1,6
          if((varf(139).eq.0.).and.(.not.ppas)) then
c            le pas de temps actuel est nul on affiche les var0 du pas
c            precedent converge          
             print*,'Fluendo3d Etape 14-1 avec dt=0 au pas actuel' 
             print*,'depsk0(',i,')=',var0(122+i)
             print*,'depskf(',i,')=',varf(122+i)             
             print*,'var0(',i,')=',var0(7+i)
             print*,'varf(',i,')=',varf(7+i)
c             read*
c             err1=1
          end if 
      end do
      do i=1,6
          if((var0(139).eq.0.).and.(.not.ppas).and.
     #      (abs(varf(7+i)-var0(7+i)).gt.1.0d-10)) then
c            on teste l evolution des vari visqueuses sur le pas
c            suivant la non covergence forcee, si on constate une
c            evolution des vari au pas actuel on l affiche     
             print*,'Fluendo3d Etape 14-2 avec dt=0 au pas precedent' 
             print*,'depsk0(',i,')=',var0(122+i)
             print*,'depskf(',i,')=',varf(122+i)             
             print*,'var0(',i,')=',var0(7+i)
             print*,'varf(',i,')=',varf(7+i)
c             err1=1
          end if 
      end do
 
c     dernier traitement des erreurs        
      if(err1.eq.1) then        
             ierr1=1
             return 
      end if 
 
c     cas ou on souhaite stopper le calcul si convergence forcee        
      if((varf(139).eq.0.).and.(.not.ppas)) then
             print*,'on sort de fldo3d avec le pas de temps nul' 
             print*,'---------------------------------------------'              
c             ierr1=1
      end if
 
c     variables internes fin de passage
C      do i=1,6
C            write(*,'(I4,1Xe10.3,1xe10.3)') 19+i,var0(19+i),varf(19+i)
C      end do
C      if(log_plastc) then
C             print*,'---on vient de plastifier en compression------'
C c            read*
C      end if
 
       return
      end
c***********************************************************************