C CFLUE2    SOURCE    CB215821  25/04/08    21:15:08     12227          
      SUBROUTINE CFLUE2(wrk52,wrk53,wrk54,wrk2,wrk3,
     &  IB,IGAU,NBPGAU,NBGMAT,NELMAT,IFOURB)
*
* modele `lemaitre fluage` Lemaitre et Chaboche pp.289,411
* voir aussi syntheses RUPTHER SEMT/LISN/RT/99-124/A
* remarque ep et sigma' sont colineaires donc ep et (Sigma0 + K depst)' le sont
* Runge_Kutta d ordre 2 : Sigma1 = 1/2 YK0 + 1/2 YK1
* YK0 = seq0 * (Sigma0 + K depst) / seqtot
* etablit increment de r et D
* YK1 = seq1 * (YK0 + K depst) / seqtot
* etablit increment de r et D
* puis moyenne
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)

-INC PPARAM
-INC CCOPTIO
-INC DECHE
*
*
*
      SEGMENT WRK2
        REAL*8 TRAC(LTRAC)
      ENDSEGMENT
*
      SEGMENT WRK3
        REAL*8 WORK(LW),WORK2(LW2)
      ENDSEGMENT
*
      dimension YK0(8),YK1(8),spri0(8),delep0(8),spri1(8),delep1(8),
     &DIV(8),CRIGI(12)


      d0 = var0(3)
      r0 = var0(2)
      
      XUNTIER=1.D0/3.D0

* cas isotrope
      ip1 = 4
*
      youn0 = xmat0(1)
      sigy0 = xmat0(ip1+1)
      xn0 = xmat0(ip1+2)
      xm0 = xmat0(ip1+3)
      gk0 = xmat0(ip1+4)
      alp0 = xmat0(ip1+5)
      blp0 = xmat0(ip1+6)
      rini0 = xmat0(ip1+7)
      a0 = xmat0(ip1+8)
      pk0 = xmat0(ip1+9)
        x2mu0= xmat0(1)/(1.D0+xmat0(2))

         if (d0.le.0.) d0  = 0.D0
         if (1.D0 - d0.le.0.) then
           do ic = 1,nstrs
             sigf(ic) =  0.D0
           enddo

           depse = 0.d0
           do is = 1, nstrs
             defp(is) = depst(is)
             depse = depse + defp(is)*defp(is)
           enddo

           varf(1) = var0(1) + SQRT(depse*0.5D0)
           varf(2) = var0(2)
           varf(3) = 1.D0
           goto 1002
         endif

      if (r0.lt.0.D0) then
c        write(6,*) 'erreur valeur r ', r0
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE2-3'
        call erreur(943)
        return
      endif
      if (r0.eq.0.) r0 = 1.D-10

       if(ib.eq.1 .and. igau.eq. 1) then
*        write(6,*) 't0' , youn0, sigy0, xn0 ,xm0 ,gk0,pk0
       endif

      youn1 = xmat(1)
      sigy1 = xmat(ip1+1)
      xn1 = xmat(ip1+2)
      xm1 = xmat(ip1+3)
      gk1 = xmat(ip1+4)
      alp1 = xmat(ip1+5)
      blp1 = xmat(ip1+6)
      rini1 = xmat(ip1+7)
      a1 = xmat(ip1+8)
      pk1 = xmat(ip1+9)
        x2mu1= xmat(1)/(1.D0+xmat(2))

       if(ib.eq.1 .and. igau.eq. 1) then
*        write(6,*) 't1' , youn1, sigy1, xn1 ,xm1 ,gk1,pk1
       endif
*
      delt = tempf - temp0
* gerer la variable cachee avec les variations de temperature

C---------CARACTERISTIQUES GEOMETRIQUES---------------------------------
C
C         COQUES
C
      ALFAH=1.D0
      IF(MFR.EQ.3.OR.MFR.EQ.5.OR.MFR.EQ.9) THEN
        EP1=work(1)
        IF(MFR.NE.5) ALFAH=work(2)**2
      ENDIF
C---------COQUES AVEC CT------------------------------------------------
C         ON NE TRAVAILLE QUE SUR LES 6 PREMIERES COMPOSANTES

      IF(MFR.EQ.9) THEN
        NCONT=6
      ELSE
        NCONT=NSTRS
      ENDIF

* calcul YK0
* determine direction et sens de eplas
* calcul increments de contrainte / utilise xmatf
* remarque on utilise les caracteristiques elastiques a la date finale

      do ic = 1,valmat(/1)
       xmatf(ic) = valmat(ic)
* completer pour le cas anisotrope
       if (ic.eq.1) xmatf(ic) = xmatf(ic) * (1.D0 - d0)
      enddo

      CALL CALSIG(depst,DDAUX,NSTAUX,CMATE,xmatf,VALCAR,
     &        N2EL,N2PTEL,MFR,IFOURB,IB,IGAU,EPAIST,
     &      NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,
     &    XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,dsigt,IRTD)

      IF(IRTD.NE.1) THEN
         KERRE=69
         GOTO 1010
      ENDIF

      DO I=1,NSTRS
        DSIGT(I)=SIG0(I) + dsigt(i)
      ENDDO

C---------CAS DES POUTRES-----------------------------------------------

      IF(MFR.EQ.7) THEN
        DIV(1)=1.D0/work(4)
        DIV(2)=1.D0
        DIV(3)=1.D0
        DIV(4)=work(10)/work(1)
        DIV(5)=work(11)/work(2)
        DIV(6)=work(12)/work(3)
        IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(work(1)*work(4))
        IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(work(2)*work(4))
        IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(work(3)*work(4))
        DO I=1,NCONT
           DSIGT(I)= DSIGT(I)*DIV(I)
        ENDDO
      ENDIF

C-----------------------------------------------------------------------
C         CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ
C-----------------------------------------------------------------------
      seqtot=VONMIS0(dsigt,NSTRS,MFR,IFOURB,EP1,ALFAH)

      if (seqtot.gt.0.D0) then
        seq0 = seqtot
      else
*  contrainte spherique
           do ic = 1,nstrs
             YK0(ic) =  dsigt(ic)
           enddo

           do is = 1, nstrs
             delep0(is) = 0.D0
           enddo

           delr0 = 0.D0
           deld0 = 0.D0
           goto 500
        endif

* point fixe pour determiner le multiplicateur de sigtot'/seqtot
      icaz = 1
      do ipfx = 1,10
         if(ib.eq.1.and.igau.eq.1) then
c            write(6,*) 'pt fixe' , ipfx, seq0,seqtot,icaz
         endif
         goto (70,80) icaz
 70    continue
* fonction
       if (gk1.gt.0.D0 .and.xm1.gt.0.D0 .and.r0.gt.0.D0
     &       .and.(1.D0-d0).gt. 0.D0) then
        delr0 =  seq0/(1.D0-d0)/gk1/(r0 ** (1/xm1))
       else
        delr0 = 0.D0
       endif
       delr0 = max(delr0,0.d0)
c cas delr0 nul ?
       if (xn1.gt.0.D0 .and. x2mu1.gt.0.D0) then
             delr0 = delr0 ** xn1
       else
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE2-1'
         call erreur(943)
          return
       endif
       xmult = 1.5D0 * delr0 * x2mu1  * delt / (1.D0 - d0) / (1.D0 - d0)
        if(ib.eq.1.and.igau.eq.1) then
c         write(6,*) 'delr0' , delr0
        endif
        seq01 =  seqtot - xmult
       goto 90

* fonction associee
  80   continue
         xmult = seqtot - seq0
c en principe 1. - d0 > 0
       if (1. - d0.lt.0..or.delt.lt.0..or.x2mu1.le.0.)  then
*         write(6,*) 'erreur parametres 1' , delt, x2mu1, d0
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE2-2'
         call erreur(943)
          return
       endif
        delr0 = xmult / delt * (1.D0- d0) * (1.D0 - d0) /x2mu1 / 1.5D0
       if (xn1.gt.0) then
         delr0 = delr0 ** (1/xn1)
*          write(6,*) ' delr0 ', delr0
       else
c         write(6,*) 'erreur parametres 2' , xn1
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE2-4'
         call erreur(943)
          return
       endif
       if (gk1.gt.0..and.xm1.gt.0.) then
        seq01 = delr0 *  (r0 ** (1/xm1)) * gk1 * (1. - d0)
       else
c         write(6,*) 'erreur parametres 3' , gk1,xm1
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE2-5'
         call erreur(943)
          return
       endif

 90    continue
       if (ipfx.eq.1) then
         if (seq01.lt.0) then
           icaz = 2
           seq01 = seqtot/2.
         else if (seq01.eq.0.) then
           icaz =1
           seq01 = seqtot/2.
         endif
       endif
        varseq = abs((seq01 - seq0) / seq0)
           if (ib.eq.1.and.igau.eq.1) then
c        write(6,*) 'variation relative', varseq
           endif
        seq0 = seq01
         if (seq0 .lt.0.) then
c          write(6,*) 'erreur point fixe', seqtot, seq0
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE2-6'
          call erreur(943)
          return
         endif

      enddo

  100    if (seqtot - seq0 .lt.0.) then
c          write(6,*) 'erreur point fixe', seqtot, seq0
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE2-7'
          call erreur(943)
          return
         endif

*a revoir pour alp0 et beta0 .
      xki0 = seq0
      if ((xki0*a1).gt.0.and.(1. - d0).gt.0..
     &and.pk1.gt.0.and.a1.ne.0..and.rini1.gt.0.)then
        deld0 = ((xki0 /a1) ** rini1) * ((1. - d0) ** (pk1 * (-1.)))
      else
        deld0 = 0.
      endif

        if(ib.eq.1.and.igau.eq.1) then
*          write(6,*) 'deld0' , deld0
        endif

      trsig0 = (dsigt(1) + dsigt(2) + dsigt(3)) * XUNTIER
      spri0(1) = dsigt(1) - trsig0
      spri0(2) = dsigt(2) - trsig0
      spri0(3) = dsigt(3) - trsig0
      do is = 4,nstrs
        spri0(is) = dsigt(is)
      enddo

      if (deld0*delt.lt.1.) then

       if (seq0.gt.0.and.(1. - d0).gt.0..and.seqtot - seq0.gt.0.) then
           coep0 = (seqtot - seq0) * (1. - d0) / x2mu1/ seqtot
        else
         coep0 = 0.d0
       endif
       do is = 1, nstrs
         delep0(is) = spri0(is) * coep0
       enddo

       do ic = 1,nstrs
        YK0(ic) = spri0(ic) * seq0 / seqtot
       if (ic.le.3) YK0(ic) = YK0(ic) + trsig0 * seq0/seqtot
       enddo

      else

       do is = 1, nstrs
         delep0(is) = depst(is)
       enddo

       do ic = 1,nstrs
        YK0(ic) = 0.
       enddo
      endif

 500  continue
* calcul YK1
* gerer la variable cachee avec les variations de temperature
         r1 = r0 + delr0 * delt
      if (r1.lt.0.) then
c        write(6,*) 'erreur valeur r ', r1
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE2-8'
        call erreur(943)
        return
      endif
       if (r1.eq.0.) r1 = 1.e-10

      d1 = var0(3) + deld0*delt
      if (d1.le.0.) d1  = 0.D0
      if (d1.ge.1.) then
        do ic = 1,nstrs
          YK1(ic) = 0.
        enddo
        deld1 = 0.
        delr1 = 0.
        do is = 1, nstrs
          delep1(is)= depst(is)
        enddo
        goto 1000
      endif

      do ic = 1,valmat(/1)
       xmatf(ic) = valmat(ic)
* completer pour le cas anisotrope
       if (ic.eq.1) xmatf(ic) = xmatf(ic) * (1. - d1)
      enddo

      CALL CALSIG(depst,DDAUX,NSTAUX,CMATE,xmatf,VALCAR,
     &        N2EL,N2PTEL,MFR,IFOURB,IB,IGAU,EPAIST,
     &      NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,
     &    XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,dsigt,IRTD)

      IF(IRTD.NE.1) THEN
         KERRE=69
         GOTO 1010
      ENDIF

      DO I=1,NSTRS
        DSIGT(I)=SIG0(I) + dsigt(i)
      ENDDO
C---------CAS DES POUTRES-----------------------------------------------

      IF(MFR.EQ.7) THEN
        DIV(1)=1.D0/work(4)
        DIV(2)=1.D0
        DIV(3)=1.D0
        DIV(4)=work(10)/work(1)
        DIV(5)=work(11)/work(2)
        DIV(6)=work(12)/work(3)
        IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(work(1)*work(4))
        IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(work(2)*work(4))
        IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(work(3)*work(4))
        DO I=1,NCONT
           DSIGT(I)= DSIGT(I)*DIV(I)
        ENDDO
      ENDIF

C-----------------------------------------------------------------------
C         CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ
C-----------------------------------------------------------------------
      seqtot=VONMIS0(dsigt,NSTRS,MFR,IFOURB,EP1,ALFAH)


      if (seqtot.gt.0.) then
        seq1 = seqtot
      else
*  contrainte spherique
           do ic = 1,nstrs
             YK1(ic) =  dsigt(ic)
           enddo

           do is = 1, nstrs
             delep1(is) = 0.
           enddo

           delr1 = 0.
           deld1 = 0.
           goto 1000
      endif

* point fixe pour determiner le multiplicateur de sigtot'/seqtot
      icaz = 1
      do ipfx = 1,10
         if(ib.eq.1.and.igau.eq.1) then
c            write(6,*) 'pt fixe' , ipfx, seq1,seqtot,icaz
         endif
         goto (570,580) icaz

 570    continue
* fonction
       if (gk1.gt.0..and.xm1.gt.0..and.r1.gt.0..and.(1.-d1).gt.0.) then
        delr1 =  seq1/(1.-d1)/gk1/(r1 ** (1/xm1))
       else
        delr1 = 0.
       endif
       delr1 = max(delr1,0.d0)
c cas delr1 nul ?
       if (xn1.gt.0.) then
             delr1 = delr1 ** xn1
       else
c         write(6,*) 'erreur parametres 2' , xn1
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE2-9'
         call erreur(943)
          return
       endif
        if(ib.eq.1.and.igau.eq.1) then
*         write(6,*) 'delr1' , delr1
        endif
       xmult = 1.5 * delr1 * x2mu1 / (1.- d1) / (1. - d1) * delt
        seq11 =  seqtot - xmult
       goto 590

* fonction associee
 580   continue
         xmult = seqtot - seq1
c en principe 1. - d1 > 0
       if (1. - d1.lt.0..or.delt.lt.0..or.x2mu1.le.0.)  then
c         write(6,*) 'erreur parametres 1' , delt, x2mu1, d1
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE210'
         call erreur(943)
          return
       endif
        delr1 = xmult / delt * (1.- d1) * (1. - d1) /x2mu1 / 1.5
       if (xn1.gt.0) then
         delr1 = delr1 ** (1/xn1)
*          write(6,*) ' delr1 ', delr1
       else
c         write(6,*) 'erreur parametres 2' , xn1
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE211'
         call erreur(943)
          return
       endif
       if (gk1.gt.0..and.xm1.gt.0.) then
        seq11 = delr1 *  (r1 ** (1/xm1)) * gk1 * (1. - d1)
       else
c         write(6,*) 'erreur parametres 3' , gk1,xm1
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE212'
         call erreur(943)
          return
       endif

 590    continue
       if (ipfx.eq.1) then
         if (seq11.lt.0) then
           icaz = 2
           seq11 = seqtot/2.
         else if (seq11.eq.0.) then
           icaz =1
           seq11 = seqtot/2.
         endif
       endif
        varseq = abs((seq11 - seq1) / seq1)
         if(ib.eq.1.and.igau.eq.1) then
c        write(6,*) 'variation relative', varseq
         endif
        seq1 = seq11
         if (seq1 .lt.0.) then
c          write(6,*) 'erreur point fixe', seqtot, seq1
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE213'
          call erreur(943)
          return
         endif


      enddo


         if (seqtot - seq1 .lt.0.) then
c          write(6,*) 'erreur point fixe', seqtot, seq1
         moterr(1:16) = conm
         moterr(17:24) = 'CFLUE214'
          call erreur(943)
          return
         endif


*a revoir pour alp1 et bta1 .
      xki1 = seq1
      if ((xki1*a1).gt.0.and.(1. - d1).gt.0..
     &and.pk1.gt.0.) then
          deld1 = ((xki1 /a1) ** rini1) * ((1. - d1) ** (pk1 * (-1.)))
      else
        deld1 = 0.
      endif
        if(ib.eq.1.and.igau.eq.1) then
*        write(6,*) 'deld1' , deld1
       endif

      trsig1 = (dsigt(1) + dsigt(2) + dsigt(3)) * XUNTIER
      spri1(1) = dsigt(1) - trsig1
      spri1(2) = dsigt(2) - trsig1
      spri1(3) = dsigt(3) - trsig1
      do is = 4,nstrs
        spri1(is) = dsigt(is)
      enddo

      if (deld1*delt.lt.1) then
       if (seq1.gt.0.and.(1. - d1).gt.0..and.seqtot - seq1.gt.0.) then
           coep1 = (seqtot - seq1) * (1. - d1) / x2mu1/ seqtot
       else
         coep1 = 0.
       endif
       do is = 1, nstrs
         delep1(is) = spri1(is) * coep1
       enddo

       do ic = 1,nstrs
         YK1(ic) = spri1(ic) * seq1 / seqtot
        if (ic.le.3) YK1(ic) = YK1(ic) + trsig1 * seq1 / seqtot
       enddo

      else

       do is = 1, nstrs
         delep1(is) = depst(is)
       enddo

       do ic = 1,nstrs
         YK1(ic) = 0.
       enddo

      endif

c au final
 1000  continue
      do ic = 1,nstrs
        sigf(ic) = 0.5*YK0(ic) + 0.5*YK1(ic)
      enddo

      depse = 0.d0
      do is = 1, nstrs
        defp(is) = 0.5*(delep0(is) + delep1(is))
        depse = depse + defp(is)*defp(is)
      enddo

      varf(1) = var0(1) + SQRT(depse*0.5)
      vtd = var0(3) + (deld0 + deld1)*0.5*delt
      if (vtd.le.1.) then
       varf(3) = vtd
      else
       varf(3) = 1.
      endif
      varf(2) = (r0 + r1)*0.5


 1002  continue
       if(ib.eq.1.and.igau.eq.1) then
*        write(6,*) 't0' ,sigf(3),yk0(3),yk1(3),varf(3),depst(3)
       endif

 1010  continue
      RETURN
      END