Télécharger la plaquette Cast3M


fr en

Télécharger cflue2.eso

Retour à la liste

Numérotation des lignes : 

cflue2

  1. C CFLUE2    SOURCE    CB215821  25/04/08    21:15:08     12227          
  2.       SUBROUTINE CFLUE2(wrk52,wrk53,wrk54,wrk2,wrk3,
  3.      &  IB,IGAU,NBPGAU,NBGMAT,NELMAT,IFOURB)
  4. *
  5. * modele `lemaitre fluage` Lemaitre et Chaboche pp.289,411
  6. * voir aussi syntheses RUPTHER SEMT/LISN/RT/99-124/A
  7. * remarque ep et sigma' sont colineaires donc ep et (Sigma0 + K depst)' le sont
  8. * Runge_Kutta d ordre 2 : Sigma1 = 1/2 YK0 + 1/2 YK1
  9. * YK0 = seq0 * (Sigma0 + K depst) / seqtot
  10. * etablit increment de r et D
  11. * YK1 = seq1 * (YK0 + K depst) / seqtot
  12. * etablit increment de r et D
  13. * puis moyenne
  14.       IMPLICIT INTEGER(I-N)
  15.       IMPLICIT REAL*8 (A-H,O-Z)
  16.  
  17. -INC PPARAM
  18. -INC CCOPTIO
  19. -INC DECHE
  20. *
  21. *
  22. *
  23.       SEGMENT WRK2
  24.         REAL*8 TRAC(LTRAC)
  25.       ENDSEGMENT
  26. *
  27.       SEGMENT WRK3
  28.         REAL*8 WORK(LW),WORK2(LW2)
  29.       ENDSEGMENT
  30. *
  31.       dimension YK0(8),YK1(8),spri0(8),delep0(8),spri1(8),delep1(8),
  32.      &DIV(8),CRIGI(12)
  33.  
  34.  
  35.       d0 = var0(3)
  36.       r0 = var0(2)
  37.  
  38.       XUNTIER=1.D0/3.D0
  39.  
  40. * cas isotrope
  41.       ip1 = 4
  42. *
  43.       youn0 = xmat0(1)
  44.       sigy0 = xmat0(ip1+1)
  45.       xn0 = xmat0(ip1+2)
  46.       xm0 = xmat0(ip1+3)
  47.       gk0 = xmat0(ip1+4)
  48.       alp0 = xmat0(ip1+5)
  49.       blp0 = xmat0(ip1+6)
  50.       rini0 = xmat0(ip1+7)
  51.       a0 = xmat0(ip1+8)
  52.       pk0 = xmat0(ip1+9)
  53.         x2mu0= xmat0(1)/(1.D0+xmat0(2))
  54.  
  55.          if (d0.le.0.) d0  = 0.D0
  56.          if (1.D0 - d0.le.0.) then
  57.            do ic = 1,nstrs
  58.              sigf(ic) =  0.D0
  59.            enddo
  60.  
  61.            depse = 0.d0
  62.            do is = 1, nstrs
  63.              defp(is) = depst(is)
  64.              depse = depse + defp(is)*defp(is)
  65.            enddo
  66.  
  67.            varf(1) = var0(1) + SQRT(depse*0.5D0)
  68.            varf(2) = var0(2)
  69.            varf(3) = 1.D0
  70.            goto 1002
  71.          endif
  72.  
  73.       if (r0.lt.0.D0) then
  74. c        write(6,*) 'erreur valeur r ', r0
  75.          moterr(1:16) = conm
  76.          moterr(17:24) = 'CFLUE2-3'
  77.         call erreur(943)
  78.         return
  79.       endif
  80.       if (r0.eq.0.) r0 = 1.D-10
  81.  
  82.        if(ib.eq.1 .and. igau.eq. 1) then
  83. *        write(6,*) 't0' , youn0, sigy0, xn0 ,xm0 ,gk0,pk0
  84.        endif
  85.  
  86.       youn1 = xmat(1)
  87.       sigy1 = xmat(ip1+1)
  88.       xn1 = xmat(ip1+2)
  89.       xm1 = xmat(ip1+3)
  90.       gk1 = xmat(ip1+4)
  91.       alp1 = xmat(ip1+5)
  92.       blp1 = xmat(ip1+6)
  93.       rini1 = xmat(ip1+7)
  94.       a1 = xmat(ip1+8)
  95.       pk1 = xmat(ip1+9)
  96.         x2mu1= xmat(1)/(1.D0+xmat(2))
  97.  
  98.        if(ib.eq.1 .and. igau.eq. 1) then
  99. *        write(6,*) 't1' , youn1, sigy1, xn1 ,xm1 ,gk1,pk1
  100.        endif
  101. *
  102.       delt = tempf - temp0
  103. * gerer la variable cachee avec les variations de temperature
  104.  
  105. C---------CARACTERISTIQUES GEOMETRIQUES---------------------------------
  106. C
  107. C         COQUES
  108. C
  109.       ALFAH=1.D0
  110.       IF(MFR.EQ.3.OR.MFR.EQ.5.OR.MFR.EQ.9) THEN
  111.         EP1=work(1)
  112.         IF(MFR.NE.5) ALFAH=work(2)**2
  113.       ENDIF
  114. C---------COQUES AVEC CT------------------------------------------------
  115. C         ON NE TRAVAILLE QUE SUR LES 6 PREMIERES COMPOSANTES
  116.  
  117.       IF(MFR.EQ.9) THEN
  118.         NCONT=6
  119.       ELSE
  120.         NCONT=NSTRS
  121.       ENDIF
  122.  
  123. * calcul YK0
  124. * determine direction et sens de eplas
  125. * calcul increments de contrainte / utilise xmatf
  126. * remarque on utilise les caracteristiques elastiques a la date finale
  127.  
  128.       do ic = 1,valmat(/1)
  129.        xmatf(ic) = valmat(ic)
  130. * completer pour le cas anisotrope
  131.        if (ic.eq.1) xmatf(ic) = xmatf(ic) * (1.D0 - d0)
  132.       enddo
  133.  
  134.       CALL CALSIG(depst,DDAUX,NSTAUX,CMATE,xmatf,VALCAR,
  135.      &        N2EL,N2PTEL,MFR,IFOURB,IB,IGAU,EPAIST,
  136.      &      NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,
  137.      &    XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,dsigt,IRTD)
  138.  
  139.       IF(IRTD.NE.1) THEN
  140.          KERRE=69
  141.          GOTO 1010
  142.       ENDIF
  143.  
  144.       DO I=1,NSTRS
  145.         DSIGT(I)=SIG0(I) + dsigt(i)
  146.       ENDDO
  147.  
  148. C---------CAS DES POUTRES-----------------------------------------------
  149.  
  150.       IF(MFR.EQ.7) THEN
  151.         DIV(1)=1.D0/work(4)
  152.         DIV(2)=1.D0
  153.         DIV(3)=1.D0
  154.         DIV(4)=work(10)/work(1)
  155.         DIV(5)=work(11)/work(2)
  156.         DIV(6)=work(12)/work(3)
  157.         IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(work(1)*work(4))
  158.         IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(work(2)*work(4))
  159.         IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(work(3)*work(4))
  160.         DO I=1,NCONT
  161.            DSIGT(I)= DSIGT(I)*DIV(I)
  162.         ENDDO
  163.       ENDIF
  164.  
  165. C-----------------------------------------------------------------------
  166. C         CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ
  167. C-----------------------------------------------------------------------
  168.       seqtot=VONMIS0(dsigt,NSTRS,MFR,IFOURB,EP1,ALFAH)
  169.  
  170.       if (seqtot.gt.0.D0) then
  171.         seq0 = seqtot
  172.       else
  173. *  contrainte spherique
  174.            do ic = 1,nstrs
  175.              YK0(ic) =  dsigt(ic)
  176.            enddo
  177.  
  178.            do is = 1, nstrs
  179.              delep0(is) = 0.D0
  180.            enddo
  181.  
  182.            delr0 = 0.D0
  183.            deld0 = 0.D0
  184.            goto 500
  185.         endif
  186.  
  187. * point fixe pour determiner le multiplicateur de sigtot'/seqtot
  188.       icaz = 1
  189.       do ipfx = 1,10
  190.          if(ib.eq.1.and.igau.eq.1) then
  191. c            write(6,*) 'pt fixe' , ipfx, seq0,seqtot,icaz
  192.          endif
  193.          goto (70,80) icaz
  194.  70    continue
  195. * fonction
  196.        if (gk1.gt.0.D0 .and.xm1.gt.0.D0 .and.r0.gt.0.D0
  197.      &       .and.(1.D0-d0).gt. 0.D0) then
  198.         delr0 =  seq0/(1.D0-d0)/gk1/(r0 ** (1/xm1))
  199.        else
  200.         delr0 = 0.D0
  201.        endif
  202.        delr0 = max(delr0,0.d0)
  203. c cas delr0 nul ?
  204.        if (xn1.gt.0.D0 .and. x2mu1.gt.0.D0) then
  205.              delr0 = delr0 ** xn1
  206.        else
  207.          moterr(1:16) = conm
  208.          moterr(17:24) = 'CFLUE2-1'
  209.          call erreur(943)
  210.           return
  211.        endif
  212.        xmult = 1.5D0 * delr0 * x2mu1  * delt / (1.D0 - d0) / (1.D0 - d0)
  213.         if(ib.eq.1.and.igau.eq.1) then
  214. c         write(6,*) 'delr0' , delr0
  215.         endif
  216.         seq01 =  seqtot - xmult
  217.        goto 90
  218.  
  219. * fonction associee
  220.   80   continue
  221.          xmult = seqtot - seq0
  222. c en principe 1. - d0 > 0
  223.        if (1. - d0.lt.0..or.delt.lt.0..or.x2mu1.le.0.)  then
  224. *         write(6,*) 'erreur parametres 1' , delt, x2mu1, d0
  225.          moterr(1:16) = conm
  226.          moterr(17:24) = 'CFLUE2-2'
  227.          call erreur(943)
  228.           return
  229.        endif
  230.         delr0 = xmult / delt * (1.D0- d0) * (1.D0 - d0) /x2mu1 / 1.5D0
  231.        if (xn1.gt.0) then
  232.          delr0 = delr0 ** (1/xn1)
  233. *          write(6,*) ' delr0 ', delr0
  234.        else
  235. c         write(6,*) 'erreur parametres 2' , xn1
  236.          moterr(1:16) = conm
  237.          moterr(17:24) = 'CFLUE2-4'
  238.          call erreur(943)
  239.           return
  240.        endif
  241.        if (gk1.gt.0..and.xm1.gt.0.) then
  242.         seq01 = delr0 *  (r0 ** (1/xm1)) * gk1 * (1. - d0)
  243.        else
  244. c         write(6,*) 'erreur parametres 3' , gk1,xm1
  245.          moterr(1:16) = conm
  246.          moterr(17:24) = 'CFLUE2-5'
  247.          call erreur(943)
  248.           return
  249.        endif
  250.  
  251.  90    continue
  252.        if (ipfx.eq.1) then
  253.          if (seq01.lt.0) then
  254.            icaz = 2
  255.            seq01 = seqtot/2.
  256.          else if (seq01.eq.0.) then
  257.            icaz =1
  258.            seq01 = seqtot/2.
  259.          endif
  260.        endif
  261.         varseq = abs((seq01 - seq0) / seq0)
  262.            if (ib.eq.1.and.igau.eq.1) then
  263. c        write(6,*) 'variation relative', varseq
  264.            endif
  265.         seq0 = seq01
  266.          if (seq0 .lt.0.) then
  267. c          write(6,*) 'erreur point fixe', seqtot, seq0
  268.          moterr(1:16) = conm
  269.          moterr(17:24) = 'CFLUE2-6'
  270.           call erreur(943)
  271.           return
  272.          endif
  273.  
  274.       enddo
  275.  
  276.   100    if (seqtot - seq0 .lt.0.) then
  277. c          write(6,*) 'erreur point fixe', seqtot, seq0
  278.          moterr(1:16) = conm
  279.          moterr(17:24) = 'CFLUE2-7'
  280.           call erreur(943)
  281.           return
  282.          endif
  283.  
  284. *a revoir pour alp0 et beta0 .
  285.       xki0 = seq0
  286.       if ((xki0*a1).gt.0.and.(1. - d0).gt.0..
  287.      &and.pk1.gt.0.and.a1.ne.0..and.rini1.gt.0.)then
  288.         deld0 = ((xki0 /a1) ** rini1) * ((1. - d0) ** (pk1 * (-1.)))
  289.       else
  290.         deld0 = 0.
  291.       endif
  292.  
  293.         if(ib.eq.1.and.igau.eq.1) then
  294. *          write(6,*) 'deld0' , deld0
  295.         endif
  296.  
  297.       trsig0 = (dsigt(1) + dsigt(2) + dsigt(3)) * XUNTIER
  298.       spri0(1) = dsigt(1) - trsig0
  299.       spri0(2) = dsigt(2) - trsig0
  300.       spri0(3) = dsigt(3) - trsig0
  301.       do is = 4,nstrs
  302.         spri0(is) = dsigt(is)
  303.       enddo
  304.  
  305.       if (deld0*delt.lt.1.) then
  306.  
  307.        if (seq0.gt.0.and.(1. - d0).gt.0..and.seqtot - seq0.gt.0.) then
  308.            coep0 = (seqtot - seq0) * (1. - d0) / x2mu1/ seqtot
  309.         else
  310.          coep0 = 0.d0
  311.        endif
  312.        do is = 1, nstrs
  313.          delep0(is) = spri0(is) * coep0
  314.        enddo
  315.  
  316.        do ic = 1,nstrs
  317.         YK0(ic) = spri0(ic) * seq0 / seqtot
  318.        if (ic.le.3) YK0(ic) = YK0(ic) + trsig0 * seq0/seqtot
  319.        enddo
  320.  
  321.       else
  322.  
  323.        do is = 1, nstrs
  324.          delep0(is) = depst(is)
  325.        enddo
  326.  
  327.        do ic = 1,nstrs
  328.         YK0(ic) = 0.
  329.        enddo
  330.       endif
  331.  
  332.  500  continue
  333. * calcul YK1
  334. * gerer la variable cachee avec les variations de temperature
  335.          r1 = r0 + delr0 * delt
  336.       if (r1.lt.0.) then
  337. c        write(6,*) 'erreur valeur r ', r1
  338.          moterr(1:16) = conm
  339.          moterr(17:24) = 'CFLUE2-8'
  340.         call erreur(943)
  341.         return
  342.       endif
  343.        if (r1.eq.0.) r1 = 1.e-10
  344.  
  345.       d1 = var0(3) + deld0*delt
  346.       if (d1.le.0.) d1  = 0.D0
  347.       if (d1.ge.1.) then
  348.         do ic = 1,nstrs
  349.           YK1(ic) = 0.
  350.         enddo
  351.         deld1 = 0.
  352.         delr1 = 0.
  353.         do is = 1, nstrs
  354.           delep1(is)= depst(is)
  355.         enddo
  356.         goto 1000
  357.       endif
  358.  
  359.       do ic = 1,valmat(/1)
  360.        xmatf(ic) = valmat(ic)
  361. * completer pour le cas anisotrope
  362.        if (ic.eq.1) xmatf(ic) = xmatf(ic) * (1. - d1)
  363.       enddo
  364.  
  365.       CALL CALSIG(depst,DDAUX,NSTAUX,CMATE,xmatf,VALCAR,
  366.      &        N2EL,N2PTEL,MFR,IFOURB,IB,IGAU,EPAIST,
  367.      &      NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,
  368.      &    XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,dsigt,IRTD)
  369.  
  370.       IF(IRTD.NE.1) THEN
  371.          KERRE=69
  372.          GOTO 1010
  373.       ENDIF
  374.  
  375.       DO I=1,NSTRS
  376.         DSIGT(I)=SIG0(I) + dsigt(i)
  377.       ENDDO
  378. C---------CAS DES POUTRES-----------------------------------------------
  379.  
  380.       IF(MFR.EQ.7) THEN
  381.         DIV(1)=1.D0/work(4)
  382.         DIV(2)=1.D0
  383.         DIV(3)=1.D0
  384.         DIV(4)=work(10)/work(1)
  385.         DIV(5)=work(11)/work(2)
  386.         DIV(6)=work(12)/work(3)
  387.         IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(work(1)*work(4))
  388.         IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(work(2)*work(4))
  389.         IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(work(3)*work(4))
  390.         DO I=1,NCONT
  391.            DSIGT(I)= DSIGT(I)*DIV(I)
  392.         ENDDO
  393.       ENDIF
  394.  
  395. C-----------------------------------------------------------------------
  396. C         CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ
  397. C-----------------------------------------------------------------------
  398.       seqtot=VONMIS0(dsigt,NSTRS,MFR,IFOURB,EP1,ALFAH)
  399.  
  400.  
  401.       if (seqtot.gt.0.) then
  402.         seq1 = seqtot
  403.       else
  404. *  contrainte spherique
  405.            do ic = 1,nstrs
  406.              YK1(ic) =  dsigt(ic)
  407.            enddo
  408.  
  409.            do is = 1, nstrs
  410.              delep1(is) = 0.
  411.            enddo
  412.  
  413.            delr1 = 0.
  414.            deld1 = 0.
  415.            goto 1000
  416.       endif
  417.  
  418. * point fixe pour determiner le multiplicateur de sigtot'/seqtot
  419.       icaz = 1
  420.       do ipfx = 1,10
  421.          if(ib.eq.1.and.igau.eq.1) then
  422. c            write(6,*) 'pt fixe' , ipfx, seq1,seqtot,icaz
  423.          endif
  424.          goto (570,580) icaz
  425.  
  426.  570    continue
  427. * fonction
  428.        if (gk1.gt.0..and.xm1.gt.0..and.r1.gt.0..and.(1.-d1).gt.0.) then
  429.         delr1 =  seq1/(1.-d1)/gk1/(r1 ** (1/xm1))
  430.        else
  431.         delr1 = 0.
  432.        endif
  433.        delr1 = max(delr1,0.d0)
  434. c cas delr1 nul ?
  435.        if (xn1.gt.0.) then
  436.              delr1 = delr1 ** xn1
  437.        else
  438. c         write(6,*) 'erreur parametres 2' , xn1
  439.          moterr(1:16) = conm
  440.          moterr(17:24) = 'CFLUE2-9'
  441.          call erreur(943)
  442.           return
  443.        endif
  444.         if(ib.eq.1.and.igau.eq.1) then
  445. *         write(6,*) 'delr1' , delr1
  446.         endif
  447.        xmult = 1.5 * delr1 * x2mu1 / (1.- d1) / (1. - d1) * delt
  448.         seq11 =  seqtot - xmult
  449.        goto 590
  450.  
  451. * fonction associee
  452.  580   continue
  453.          xmult = seqtot - seq1
  454. c en principe 1. - d1 > 0
  455.        if (1. - d1.lt.0..or.delt.lt.0..or.x2mu1.le.0.)  then
  456. c         write(6,*) 'erreur parametres 1' , delt, x2mu1, d1
  457.          moterr(1:16) = conm
  458.          moterr(17:24) = 'CFLUE210'
  459.          call erreur(943)
  460.           return
  461.        endif
  462.         delr1 = xmult / delt * (1.- d1) * (1. - d1) /x2mu1 / 1.5
  463.        if (xn1.gt.0) then
  464.          delr1 = delr1 ** (1/xn1)
  465. *          write(6,*) ' delr1 ', delr1
  466.        else
  467. c         write(6,*) 'erreur parametres 2' , xn1
  468.          moterr(1:16) = conm
  469.          moterr(17:24) = 'CFLUE211'
  470.          call erreur(943)
  471.           return
  472.        endif
  473.        if (gk1.gt.0..and.xm1.gt.0.) then
  474.         seq11 = delr1 *  (r1 ** (1/xm1)) * gk1 * (1. - d1)
  475.        else
  476. c         write(6,*) 'erreur parametres 3' , gk1,xm1
  477.          moterr(1:16) = conm
  478.          moterr(17:24) = 'CFLUE212'
  479.          call erreur(943)
  480.           return
  481.        endif
  482.  
  483.  590    continue
  484.        if (ipfx.eq.1) then
  485.          if (seq11.lt.0) then
  486.            icaz = 2
  487.            seq11 = seqtot/2.
  488.          else if (seq11.eq.0.) then
  489.            icaz =1
  490.            seq11 = seqtot/2.
  491.          endif
  492.        endif
  493.         varseq = abs((seq11 - seq1) / seq1)
  494.          if(ib.eq.1.and.igau.eq.1) then
  495. c        write(6,*) 'variation relative', varseq
  496.          endif
  497.         seq1 = seq11
  498.          if (seq1 .lt.0.) then
  499. c          write(6,*) 'erreur point fixe', seqtot, seq1
  500.          moterr(1:16) = conm
  501.          moterr(17:24) = 'CFLUE213'
  502.           call erreur(943)
  503.           return
  504.          endif
  505.  
  506.  
  507.       enddo
  508.  
  509.  
  510.          if (seqtot - seq1 .lt.0.) then
  511. c          write(6,*) 'erreur point fixe', seqtot, seq1
  512.          moterr(1:16) = conm
  513.          moterr(17:24) = 'CFLUE214'
  514.           call erreur(943)
  515.           return
  516.          endif
  517.  
  518.  
  519. *a revoir pour alp1 et bta1 .
  520.       xki1 = seq1
  521.       if ((xki1*a1).gt.0.and.(1. - d1).gt.0..
  522.      &and.pk1.gt.0.) then
  523.           deld1 = ((xki1 /a1) ** rini1) * ((1. - d1) ** (pk1 * (-1.)))
  524.       else
  525.         deld1 = 0.
  526.       endif
  527.         if(ib.eq.1.and.igau.eq.1) then
  528. *        write(6,*) 'deld1' , deld1
  529.        endif
  530.  
  531.       trsig1 = (dsigt(1) + dsigt(2) + dsigt(3)) * XUNTIER
  532.       spri1(1) = dsigt(1) - trsig1
  533.       spri1(2) = dsigt(2) - trsig1
  534.       spri1(3) = dsigt(3) - trsig1
  535.       do is = 4,nstrs
  536.         spri1(is) = dsigt(is)
  537.       enddo
  538.  
  539.       if (deld1*delt.lt.1) then
  540.        if (seq1.gt.0.and.(1. - d1).gt.0..and.seqtot - seq1.gt.0.) then
  541.            coep1 = (seqtot - seq1) * (1. - d1) / x2mu1/ seqtot
  542.        else
  543.          coep1 = 0.
  544.        endif
  545.        do is = 1, nstrs
  546.          delep1(is) = spri1(is) * coep1
  547.        enddo
  548.  
  549.        do ic = 1,nstrs
  550.          YK1(ic) = spri1(ic) * seq1 / seqtot
  551.         if (ic.le.3) YK1(ic) = YK1(ic) + trsig1 * seq1 / seqtot
  552.        enddo
  553.  
  554.       else
  555.  
  556.        do is = 1, nstrs
  557.          delep1(is) = depst(is)
  558.        enddo
  559.  
  560.        do ic = 1,nstrs
  561.          YK1(ic) = 0.
  562.        enddo
  563.  
  564.       endif
  565.  
  566. c au final
  567.  1000  continue
  568.       do ic = 1,nstrs
  569.         sigf(ic) = 0.5*YK0(ic) + 0.5*YK1(ic)
  570.       enddo
  571.  
  572.       depse = 0.d0
  573.       do is = 1, nstrs
  574.         defp(is) = 0.5*(delep0(is) + delep1(is))
  575.         depse = depse + defp(is)*defp(is)
  576.       enddo
  577.  
  578.       varf(1) = var0(1) + SQRT(depse*0.5)
  579.       vtd = var0(3) + (deld0 + deld1)*0.5*delt
  580.       if (vtd.le.1.) then
  581.        varf(3) = vtd
  582.       else
  583.        varf(3) = 1.
  584.       endif
  585.       varf(2) = (r0 + r1)*0.5
  586.  
  587.  
  588.  1002  continue
  589.        if(ib.eq.1.and.igau.eq.1) then
  590. *        write(6,*) 't0' ,sigf(3),yk0(3),yk1(3),varf(3),depst(3)
  591.        endif
  592.  
  593.  1010  continue
  594.       RETURN
  595.       END
  596.  
  597.  
  598.  
  599.  
  600.  
  601.  
  602.  
  603.  
  604.  
  605.  
  606.  
  607.  
  608.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales