Sources Esope | Cast3M

fluage
C FLUAGE    SOURCE    CHAT      05/01/13    00:04:07     5004
      SUBROUTINE FLUAGE(ITER,EPSTAR,S0,EPST,SSTAR,SI,DPSM1,
     1 DPSM2,IBOU,DEPS,ICONV,TIMEXI,TIMEXF,YUNG,KERRE,ecou,necou)
C
C            CE SOUS PROGRAMME CALCULE LA DEFORMATION DE FLUAGE
C
C      SORTIES AVEC PROBLEMES
C
C            KERRE = 3  DIMINUER LE PAS (TREANOR)
C            KERRE = 4  DEPS EST NEGATIF
C            KERRE = 5  PAS DE CONVERGENCE INTERNE
C                         ( ANCIENNEMENT NSTOP=0 )
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DSIGK1(6),DSIGK2(6),DSIGK3(6),DSIGK4(6)
C
       SEGMENT ECOU
***     COMMON/ECOU/TEST,ALFAH,
       REAL*8 TEST, ALFAH,
     1 HPAS,TEMPS,COVNMS(6),VECPRO(9),VALPRO(6),
     2 CVNMSD(12),STOT(6),SIGEL(6),DSIGP(6),SIGT(6),SIG5(6),SIGINT(6),
     1 DALPHA(6),DSIGO(6),E(12),XINV(3),
     2 SIPLAD(6),DSIGP0(6),TET,TETI
      ENDSEGMENT
 
C      COMMON/ECOU/TEST,ALFAH,
C     1 HPAS,TEMPS,COVNMS(6),VECPRO(9),VALPRO(6),
C     2 CVNMSD(12),STOT(6),SIGEL(6),DSIGP(6),SIGT(6),SIG5(6),SIGINT(6),
C     1 DALPHA(6),DSIGO(6),E(12),XINV(3),
C     2 SIPLAD(6),DSIGP0(6),TET,TETI
C
       SEGMENT NECOU
*      COMMON/NECOU/NCOURB,IPLAST,IT,IMAPLA,ISOTRO,
       INTEGER NCOURB,IPLAST,IT,IMAPLA,ISOTRO,
     . ITYP,IFOUR,IFLUAG,
     . ICINE,ITHER,IFLUPL,ICYCL,IBI,
     . JFLUAG,KFLUAG,LFLUAG,
     . IRELAX,JNTRIN,MFLUAG,JSOUFL,JGRDEF
      ENDSEGMENT
C      COMMON/NECOU/NCOURB,IPLAST,IT,IMAPLA,ISOTRO,
C     . ITYP,IFOUR,IFLUAG,
C     . ICINE,ITHER,IFLUPL,ICYCL,IBI,
C     . JFLUAG,KFLUAG,LFLUAG,
C     . IRELAX,JNTRIN,MFLUAG,JSOUFL,JGRDEF
C
      DATA ITMAX/15/
      ITER=1
      ICLFLU=0
      IF(IRELAX.EQ.0) GO TO 414
C
C ESSAI DE CALCUL AUTOMATIQUE DU PAS
C
C  IRELAX = 1         ALGORITHME LADOUCEUR
C  IRELAX = 2         ALGORITHME TREANOR
C
C
C  TESTS DE PRECISION POUR METHODE TREANOR
C
C  EPREC=INCREMENT DE DEFORMATION INELASTIQUE MAXIMUM (1ERE ESTIMATION)
C  EPREC2=INCREMENT DE DEFORMATION INELASTIQUE MAXIMUM (VALEUR FINALE)
C  PRERUN  TEST DE VALIDITE DU PAS DE TEMPS
C  PREDEF  VALEUR DE REFERENCE POUR LES DEFORMATIONS TROP PETITES
C  PRESIG  VALEUR DE REFERENCE POUR LES CONTRAINTES TROP PETITES
C  PREPH   TEST POUR PHMAX
C  PREBTA = BORNE SUP POUR LES ITERATIONS INTERNES DE 2-EME NIVEAU
C  RPREC  PRECISION SUR L ESTIMATION DE DELTA-EPS CUMULE ITER INT 2 NIV
C  PREZER  VALEUR POUR INTERPOLATION EN DELTA-EPSILON
C  PREMIN  VALEUR POUR ESTIMATION AUTOMATIQUE DU PAS EN TREANOR
C
      EPREC =1.D-3
      EPREC2=1.D-2
      PRECIS=TEST
      PRERUN=PRECIS
      PREDEF=1.D-5
      PRESIG=YUNG*PREDEF
      PREPH=1.D-3
      PREBTA=2.D00
      RPREC=PRECIS
      PREZER=1.D-7
      PREMIN=3.D-1*PRERUN
      IF(IT.EQ.1) GO TO 4430
      IF(DPSM1.LT.PREZER) GO TO 4430
      EPREC=MIN(EPREC,DPSM1)
      EPREC2=MIN(EPREC2,DPSM1)
4430  CONTINUE
      TIMO=TEMPS-HPAS
      TIMEXF=TIMEXI+HPAS
      DTET=TET-TETI
      DTETDT=0.D00
      IF(HPAS.NE.0.) DTETDT=DTET/HPAS
  200 CONTINUE
C  NITER= NOMBRE D'ITERATIONS SUR LE PAS
      NITER=3
C
C  BTAEPS=ITERPOLATION DU DSIGP-DSIGP0 SUIVANT LA VARIATION DE
C  EPSILON ---- L'INITIALISATION DE LA PREMIERE ITERATION EST BIDON
      BTAEPS=1.
      BTAA=1.
      IF(IT.EQ.1) GO TO 430
      BTAA=DPSM1
      IF(BTAA.GT.PREZER) BTAEPS=1./BTAA
C
C  IDECOU=NOMBRE DE DECOUPAGES POSSIBLE DU TEMPS
430   IDECOU=30
      IF(JSOUFL.NE.0) IDECOU=JSOUFL
C
C  ICONV=1 SI CONVERGENCE INTERNE NON EFFECTUEE
      ICONV=0
C
C  ICOMPT=NOMBRE D'ITERATIONS POUR AVOIR UNE BONNE ESTIMATION DE
C  LA VALEUR DE DELTA EPSILON ----- ICOMPT MAXIMUM=15
C
      ICOMPT=1
C
C  DSIGP0(IB)=PARTIE ELASTIQUE DU DSIGP(IB)--- INTERPOLEE EN TEMPS
C  DSIGP-DSIGP0=PARTIE INELASTIQUE DU DSIGP(IB)--- INTERPOLEE EN EPSILON
C  TEMINT=TEMPS INTERMEDIAIRE AU DEBUT DU PAS
C  TEMIN2=TEMPS INTERMEDIAIRE A LA FIN DU PAS
C  TRINS1  TEMPS INTRINSEQUE INTERMEDIAIRE DEBUT DU SS-PAS
C  TRINS2  TEMPS INTRINSEQUE INTERMEDIAIRE FIN   DU SS-PAS
C  EPSINT=DEFORMATION DE FLUAGE AU DEBUT DES TEMPS INTERMEDIAIRES
C  EPSIN2=DEFORMATION DE FLUAGE A LA FIN DES TEMPS INTERMEDIAIRES
C  SINT=CONTRAINTE CALCULEE AU DEBUT DU PAS INTERMEDIAIRE
C  SINT2=CONTRAINTE CALCULEE A LA FIN DU PAS INTERMEDIAIRE
C  SIGINT(IB)=TENSEUR DE CONTRAINTE AU DEBUT DES PAS INTERMEDIAIRES
C  TETINI = TEMPERATURE INTERMEDIAIRE AU DEBUT DU SS-PAS
C  TETINF = TEMPERATURE INTERMEDIAIRE A LA FIN DU SS-PAS
C  HPAS2 = PAS DE TEMPS ITERMEDIAIRES
C-----------------------------------------------------------------------
C                      ON DEMARRE
C-----------------------------------------------------------------------
499   TEMINT=TIMO
      TRINS1=TIMEXI
      EPSINT=EPSTAR
      SINT=S0
      DO 413 IB=1,IBOU
413   SIGINT(IB)=SIGEL(IB)
C-----------------------------------------------------------------------
C                    BOUCLE SUR LES SOUS-PAS
C-----------------------------------------------------------------------
      DO 416 I=1,IDECOU
C    ON INTERPOLE LINEAIREMENT
C
      TETINI=TETI+DTETDT*(TEMINT-TIMO)
      DECATE=TEMPS-TEMINT
      HPAS2=DECATE
      IF(IRELAX.NE.2) GO TO 8370
C
C            MODIFS POUR LA METHODE DE TREANOR
C  POUR LA PREMIERE ESTIMATION DU PAS DANS LES ITERATIONS EN DEPS,
C  ON UTILISE LA VALEUR HFIRST TROUVEE A L ITERATION PRECEDENTE
C  POUR LE DECOUPAGE EN SOUS-PAS,ON AJUSTE AUTOMATIQUEMENT LA TAILLE
C  DU PAS .
C
      IF(I.GT.1) GO TO 8371
      IF(ICOMPT.EQ.1) GO TO 8370
      HPAS2=HFIRST
      GO TO 8370
8371  HPAS2=HPADEC
8370  BETA=HPAS2/DECATE
      IF(BETA.GT.1.) BETA=1.
      CALL CRPLAW(VI,EPSINT,SINT,TETINI,TRINS1,HPAS2,ICLFLU,NCOURB)
      EPSPAS=HPAS2*VI
C  EPSPAS=INCREMENT D'ALLONGEMENT INELASTIQUE CALCULE A LA FIN
C  DU SOUS PAS INTERMEDIAIRE
C  BETA=COEFFICIENT CORRESPONDANT A LA FRACTION DU PAS CALCULE
C
      IF(EPSPAS.GT.EPREC) BETA=BETA*EPREC/EPSPAS
C
C-----------------------------------------------------------------------
C        DEBUT DES CALCULS POUR UN SOUS-PAS DE TAILLE DONNEE
C        ON RECOMMENCE ICI EN CAS DE REDUCTION DU PAS
C-----------------------------------------------------------------------
C
 424  HPAS2=DECATE*BETA
      BTAPAS=HPAS2/HPAS
      TEMIN2=TEMINT+HPAS2
      TRINS2=TRINS1+HPAS2
      TETINF=TETI+DTETDT*(TEMIN2-TIMO)
      TIMO5=TEMINT+0.5*HPAS2
      TRINS5=TRINS1+0.5*HPAS2
      TETIN5=TETI+DTETDT*(TIMO5-TIMO)
      IF(IRELAX.EQ.2) GO TO 1101
C
C  CALCUL PAR LA METHODE LADOUCEUR
C
      EPSPA2=EPSPAS*BETA
      EPSIN2=EPSINT+EPSPA2
      BTASIG=EPSPA2*BTAEPS
      IF(BTASIG.GT.PREBTA) BTASIG=PREBTA
C
C
C  IL FAUT TOUJOURS INTERPOLER EN DELTA-EPSILON, MEME SI DEPS EST PETIT,
C  POUR QUE LA RUSE DE LA 2-EME ITERATION MARCHE EN FLUAGE PUR
C
C
      IF(BTAA.LE.PREZER) BTASIG=BTAPAS
      ELTA=0.
      IF(SINT.NE.0.) ELTA=VI/SINT
      ELTF=HPAS2*ELTA
      DO 421 IB=1,IBOU
421   SIGT(IB)=SIGINT(IB)+DSIGP0(IB)*BTAPAS+
     1 (DSIGP(IB)-DSIGP0(IB))*BTASIG
     2 -E(IB)*SIGINT(IB)*ELTF
      SINT2=VNMISD(SIGT,ITYP,ALFAH,CVNMSD)
C
C DEBUT DES ITERATIONS SUR LA FRACTION DE PAS CONSIDERE
C
      DO 417 J=1,NITER
      CALL CRPLAW(VINT,EPSIN2,SINT2,TETINF,TRINS2,HPAS2,ICLFLU,NCOURB)
      ELTO=.25*HPAS2
      IF(J.EQ.1) EPSPA5=ELTO*(3*VI+VINT)*.5
      IF(J.NE.1.AND.(VI+VINT).NE.0.)
     1 EPSPA5=DEPS*(3*VI+VINT)/(VI+VINT)*.25
      EPSO5=EPSINT+EPSPA5
      BTASIG=EPSPA5*BTAEPS
      IF(BTASIG.GT.PREBTA) BTASIG=PREBTA
      IF(BTAA.LE.PREZER) BTASIG=BTAPAS*.5
      ELTB=0.
      IF(SINT2.NE.0.) ELTB=VINT/SINT2
      DO 422 IB=1,IBOU
      IF(J.EQ.1) SIG5(IB)=SIGINT(IB)+DSIGP0(IB)*BTAPAS*.5+
     1 (DSIGP(IB)-DSIGP0(IB))*BTASIG
     2 -ELTO*E(IB)*(3*SIGINT(IB)*ELTA+SIGT(IB)*ELTB)*.5
      IF(J.NE.1) SIG5(IB)=SIGINT(IB)+DSIGP0(IB)*BTAPAS*.5+
     1 (DSIGP(IB)-DSIGP0(IB))*BTASIG
     2 -E(IB)*ELTO*(SIGINT(IB)*ELTA+SIG5(IB)*ELTC*.25)
422   CONTINUE
      SO5=VNMISD(SIG5,ITYP,ALFAH,CVNMSD)
      CALL CRPLAW(VO5,EPSO5,SO5,TETIN5,TRINS5,HPAS2,ICLFLU,NCOURB)
      ELTO=.166666666667*HPAS2
      DEPS=ELTO*(VI+VINT+4*VO5)
      IF(DEPS.LT.EPREC2) GO TO 425
      BETA=BETA*.7
      GO TO 424
425   EPSIN2=EPSINT+DEPS
      BTASIG=DEPS*BTAEPS
      IF(BTASIG.GT.PREBTA) BTASIG=PREBTA
      IF(BTAA.LE.PREZER) BTASIG=BTAPAS
      ELTC=0.
      IF(SO5.NE.0.) ELTC=VO5/SO5*4.
      DO 423 IB=1,IBOU
423   SIGT(IB)=SIGINT(IB)+DSIGP0(IB)*BTAPAS+
     1 (DSIGP(IB)-DSIGP0(IB))*BTASIG
     2 -ELTO*E(IB)*(SIGINT(IB)*ELTA+SIGT(IB)*ELTB+SIG5(IB)*ELTC)
      SINT2=VNMISD(SIGT,ITYP,ALFAH,CVNMSD)
417   CONTINUE
C
      GO TO 1200
C
 1101 CONTINUE
C
C  CALCUL PAR LA METHODE DE TREANOR
C
      WEP=0.5D00
      CALL FLURUN(SIGINT,EPSINT,TETINI,TRINS1,HPAS2,ICLFLU,NCOURB,IBOU,
     . ITYP,BTAA,BETA,BTAEPS,BTAPAS,ICOD,DSIGK1,DEPSK1,SIGT,EPSINA,WEP,
     . EPSINT,EPREC2,ecou)
      IF(ICOD.EQ.1) GO TO 424
      CALL FLURUN(SIGT,EPSINA,TETIN5,TRINS5,HPAS2,ICLFLU,NCOURB,IBOU,
     . ITYP,BTAA,BETA,BTAEPS,BTAPAS,ICOD,DSIGK2,DEPSK2,SIGT,EPSINA,WEP,
     . EPSINT,EPREC2,ecou)
      IF(ICOD.EQ.1) GO TO 424
      WEP=1.D00
      CALL FLURUN(SIGT,EPSINA,TETIN5,TRINS5,HPAS2,ICLFLU,NCOURB,IBOU,
     . ITYP,BTAA,BETA,BTAEPS,BTAPAS,ICOD,DSIGK3,DEPSK3,SIGT,EPSINA,WEP,
     . EPSINT,EPREC2,ecou)
      IF(ICOD.EQ.1) GO TO 424
      PHMAX=0.
      IF(DEPSK2-DEPSK1.NE.0.) PHMAX=2.*(DEPSK2-DEPSK3)/(DEPSK2-DEPSK1)
      DO 4101 IB=1,IBOU
      IF(DSIGK2(IB)-DSIGK1(IB).EQ.0.) GO TO 4101
      PH=2.*(DSIGK2(IB)-DSIGK3(IB))/(DSIGK2(IB)-DSIGK1(IB))
      IF(PHMAX.LT.PH) PHMAX=PH
4101  CONTINUE
      IF(PHMAX.LE.0.) GO TO 4102
      IF(PHMAX.GE.PREPH) GO TO 4123
      XL1=1.-PHMAX*0.5
      XL2=0.5-PHMAX*0.166666666666666667
      XL3=0.166666666666666667-PHMAX*0.04166666666666666667
      XL4=-PHMAX*0.166666666666666667
      GO TO 4103
4123  CONTINUE
      XL1=(1.-EXP(-PHMAX))/PHMAX
      XL2=(1.-XL1)/PHMAX
      XL3=(0.5-XL2)/PHMAX
      XL4=XL1-2.*XL2
      GO TO 4103
4102  PHMAX=0.
      XL1=1.
      XL2=0.5
      XL3=0.166666666666666667
      XL4=0.
4103  CONTINUE
      XM1=-XL2+4.*XL3
      XM2=2.*(XL2-2.*XL3)
      XM3=4.*XL3-3.*XL2
C
C  NOUVELLE ESTIMATION DE LA SOLUTION
C
      EPSIN5=EPSINT+2.*XL2*DEPSK3+XL2*PHMAX*DEPSK2+XL4*DEPSK1
      DO 4104 IB=1,IBOU
4104  SIG5(IB)=SIGINT(IB)+2.*XL2*DSIGK3(IB)+XL2*PHMAX*DSIGK2(IB)
     . +XL4*DSIGK1(IB)
      WEP=1.D00
      CALL FLURUN(SIG5,EPSIN5,TETINF,TRINS2,HPAS2,ICLFLU,NCOURB,
     . IBOU,ITYP,BTAA,BETA,BTAEPS,BTAPAS,ICOD,DSIGK4,DEPSK4,SIGT,
     . EPSINA,WEP,EPSINT,EPREC2,ecou)
      IF(ICOD.EQ.1) GO TO  424
C
C  ESTIMATION FINALE DE LA SOLUTION
C
      FAC1=1.+PHMAX*XM3+2.*PHMAX*XM2
      FAC2=XL1+XM3+0.5*XM2*PHMAX
      FAC3=XM2*(1.+0.5*PHMAX)
      EPSIN2=EPSINT*FAC1+DEPSK1*FAC2+DEPSK2*FAC3+DEPSK3*XM2
     . +DEPSK4*XM1+EPSIN5*XM1*PHMAX
      DO 4105 IB=1,IBOU
4105  SIGT(IB)=SIGINT(IB)*FAC1+DSIGK1(IB)*FAC2+DSIGK2(IB)*FAC3
     . +DSIGK3(IB)*XM2+DSIGK4(IB)*XM1+SIG5(IB)*XM1*PHMAX
C
C  TEST DE PRECISION
C
      ERMAX=ABS(EPSIN2-EPSIN5)/MAX(PREDEF,ABS(EPSIN2-EPSINT))
      DO 4106 IB=1,IBOU
      ERM=ABS(SIGT(IB)-SIG5(IB))/MAX(PRESIG,
     .  ABS(SIGT(IB)-SIGINT(IB)))
      IF(ERMAX.LT.ERM) ERMAX=ERM
4106  CONTINUE
      IF(ERMAX.LE.PRERUN) GO TO 1240
      BETA=BETA*0.5
      GO TO 424
C
C-----------------------------------------------------------------------
C  METHODE DE TREANOR  - ESTIMATION AUTOMATIQUE DU PAS DE TEMPS
C-----------------------------------------------------------------------
C
1240  HPADEC=HPAS2
      IF(ERMAX.LE.PREMIN) HPADEC=2.D00*HPAS2
      SINT2=VNMISD(SIGT,ITYP,ALFAH,CVNMSD)
 1200 CONTINUE
C
C   CAS DU PREMIER DECOUPAGE : ON SAUVE LA TAILLE DU 1-ER SOUS-PAS
C
      IF(I.EQ.1.AND.ICOMPT.EQ.1) HFIRST=HPAS2
      DO 418 IB=1,IBOU
418   SIGINT(IB)=SIGT(IB)
C
C-----------------------------------------------------------------------
C        SORTIE DE LA BOUCLE DES SOUS-PAS SI BETA=1.
C        LE DECOUPAGE EST FINI
C-----------------------------------------------------------------------
C
      IF(ABS(BETA-1.D00).LE.PREZER) GO TO 419
      EPSINT=EPSIN2
C
C-----------------------------------------------------------------------
C  SI ON N'A PAS FINI,ON REGARDE SI ON PASSE LE TEST ANTI-DIVERGENCE
C  SI CA A TENDANCE A EXPLOSER,ON RECOMMENCE AVEC UN BTAA MODIFIE
C  SINON ON SORT EN DEMANDANT PLUS DE SOUS-PAS
C-----------------------------------------------------------------------
C
      IF(IT.LE.2) GO TO 502
      IF(BTAA.LE.PREDEF) GO TO 502
C
C  TEST ANTI-DIVERGENCE ITERATIONS INTERNES DE 2-EME NIVEAU
C
      IF(DEPS.LE.1.5*BTAA) GO TO 502
      BTAA=1.1*BTAA
      BTAEPS=1./BTAA
      GO TO 499
502   CONTINUE
      SINT=SINT2
      TEMINT=TEMIN2
      IF(JNTRIN.EQ.0) TRINS1=TRINS2
      IF(JNTRIN.NE.0.AND.DEPS.NE.0.) TRINS1=TRINS2
C
C-----------------------------------------------------------------------
C                    FIN DE LA BOUCLE SUR LES SOUS-PAS
C-----------------------------------------------------------------------
C
416   CONTINUE
      KERRE=3
      RETURN
C
C  FIN DU DECOUPAGE EN SOUS-PAS
C
419   SI=SINT2
      DEPS=EPSIN2-EPSTAR
      IF(DEPS.GE.0.) GO TO 3170
      KERRE=4
      DEPS=0.D00
3170  CONTINUE
C
C  ON TESTE LA VALEUR DE DEPS PAR RAPPORT A L'ESTIMATION FAITE
C  ON CONTINUE JUSQ'A CONVERGENGE A RPREC PRES
C
      IF(IT.EQ.1) GO TO 501
      IF(IT.EQ.2) GO TO 505
C
C  SI BTAA EST INFERIEUR A PREZER,ON NE FAIT PAS D'ITERATIONS
C  INTERNES DE 2-EME NIVEAU
C
      IF(BTAA.LE.PREZER) GO TO 501
      DIFFE=DEPS-BTAA
      RAPOR=ABS(DIFFE)/BTAA
      IF(RAPOR.LE.RPREC) GO TO 501
      IF(ICOMPT.LT.15) GO TO 415
      KERRE=5
      GO TO 505
415   CONTINUE
      ICOMPT=ICOMPT+1
      IF(ICOMPT.EQ.2) GO TO 503
      DESSU=BTAAZ*DEPS-BTAA*DEPSZ
      DESSOU=DIFFE-DIFFEZ
503   BTAAZ=BTAA
      DEPSZ=DEPS
      DIFFEZ=DIFFE
      IF(ICOMPT.EQ.2) GO TO 535
      IF(ABS(DESSOU).GT.1.D-20) GO TO 504
      BTAA=(BTAA+DEPS)*0.5
      GO TO 3390
504   BTAA=DESSU/DESSOU
      GO TO 3390
535   BTAA=DEPS
3390  IF(BTAA.GT.PREZER) BTAEPS=1./BTAA
      IF(BTAA.LT.0.) BTAA=0.5D00*BTAAZ
      GO TO 499
501   CONTINUE
      EPST=EPSIN2
      IF(JNTRIN.NE.0) TIMEXF=TRINS2
      RETURN
505   CONTINUE
      EPST=EPSIN2
      IF(JNTRIN.NE.0) TIMEXF=TRINS2
      IF(BTAA.LE.PREZER) RETURN
      ICONV=1
      RDEPS=DEPS/BTAA
      DO 506 IB=1,IBOU
506   DSIGP(IB)=SIGEL(IB)+DSIGP0(IB)+(DSIGP(IB)-DSIGP0(IB))*RDEPS
      RETURN
C
C-----------------------------------------------------------------------
C                   METHODE CLASSIQUE CASTEM
C-----------------------------------------------------------------------
C
414   TIMO=TEMPS-HPAS
      TIMEXF=TIMEXI+HPAS
      TIMEX5=0.5*(TIMEXI+TIMEXF)
      TET5=0.5*(TETI+TET)
      TIMO5=0.5*(TIMO+TEMPS)
      IF(IT.NE.1) GO TO 411
C
C  RUSE POUR LA 1-ERE ITERATION EN CAS THERMOPLASTIQUE
C
      CALL CRPLAW(VI,EPSTAR,S0,TETI,TIMEXI,HPAS,ICLFLU,NCOURB)
      DPSM1=VI
      EPST=EPSTAR+HPAS*VI
      CALL CRPLAW(VF,EPST,SSTAR,TET,TIMEXF,HPAS,ICLFLU,NCOURB)
      DPSM2=SSTAR
      SF=SSTAR
      ELTO=0.166666666667*HPAS
      SO5=0.5*(S0+SF)
      EPST=EPSTAR+0.25*HPAS*(VI+VF)
      CALL CRPLAW(VO5,EPST,SO5,TET5,TIMEX5,HPAS,ICLFLU,NCOURB)
      ELTA=0
      ELTB=0
      ELTC=0
      IF(S0.NE.0.) ELTA=1./S0
      IF(SF.NE.0.) ELTB=1./SF
      IF(SO5.NE.0.) ELTC=1./SO5
      ELTF=ELTO*(VF*ELTB+2.*VO5*ELTC)
      ELTO=ELTO*(VI*ELTA+2.*VO5*ELTC)
      DO 412 IB=1,IBOU
      ELTA=1.-E(IB)*ELTO
      ELTB=1.+E(IB)*ELTF
      SIGEL(IB)=SIGEL(IB)*ELTB/ELTA
  412 DSIGP(IB)=DSIGP(IB)*ELTB
  411 VI=DPSM1
      SF=DPSM2
      EPST=EPSTAR+HPAS*VI
      CALL CRPLAW(VF,EPST,SF,TET,TIMEXF,HPAS,ICLFLU,NCOURB)
      ELTO=0.166666666667*HPAS
      SO5=0.5*(S0+SF)
      EPST=0.25*HPAS*(VI+VF)+EPSTAR
      CALL CRPLAW(VO5,EPST,SO5,TET5,TIMEX5,HPAS,ICLFLU,NCOURB)
      ELTA=0.
      ELTB=0.
      ELTC=0.
      IF(S0.NE.0.) ELTA=1./S0
      IF(SF.NE.0.) ELTB=1./SF
      IF(SO5.NE.0.) ELTC=1./SO5
      ELTF=ELTO*(VF*ELTB+2.*VO5*ELTC)
      ELTO=ELTO*(VI*ELTA+2.*VO5*ELTC)
      DO 410 IB=1,IBOU
      ELTA=1.-E(IB)*ELTO
      ELTB=1.+E(IB)*ELTF
      SIGT(IB)=(DSIGP(IB)+SIGEL(IB)*ELTA)/ELTB
      IF(IT.GT.1) GO TO 410
      DSIGP(IB)=  SIGT(IB)*(1.+ELTB-ELTA)
  410 CONTINUE
      SI=VNMISD(SIGT,ITYP,ALFAH,CVNMSD)
      DPSM2=SI
      DEPS=0.166666666667*HPAS*(VI+VF+4.*VO5)
      EPST=EPSTAR+DEPS
      IF(JNTRIN.NE.0.AND.DEPS.EQ.0.) TIMEXF=TIMEXI
   57 CONTINUE
      RETURN
      END