Télécharger la plaquette Cast3M


fr en

Télécharger prpfib.eso

Retour à la liste

Numérotation des lignes : 

prpfib

  1. C PRPFIB    SOURCE    PV090527  24/06/12    21:15:09     11940          
  2.       subroutine prpfib(p,Lf,mu,
  3.      # Td,Tmax,Ef,Hf,Df,sk,F0,alphad,M0,Rt,fu,
  4.      # L0,d,e,fy,mwor,nligne,niter,nangl,nrt,nlongfi,
  5.      # nangredu,phit,phicrit)    
  6.  
  7.  
  8. c      implicit none
  9.         implicit real*8 (a-h,o-z)
  10.         implicit integer (i-n)     
  11.  
  12.  
  13. c   ***************************************************************
  14. c   ************************ MODIFICATIONS ************************
  15. C   ***************************************************************
  16. c        variables pour la reduction de modele
  17.  
  18. c   ***************************************************************
  19. c   ***************** NOMS DES VARIABLES/PARAMETRES ***************
  20. c   ***************************************************************
  21.  
  22. c   i : boucle sur le nombre ligne
  23. c   j : boucle sur les angles
  24. c   k : boucle de sous-itération
  25. c   l : boucle sur les longueurs ancrées
  26. c   nligne : nombre de lignes pour 1 fibre
  27. c   niter : nombre d'itérations max sur k
  28. c   incremf : incrément de force élastique (MN) 
  29. c   incremld : incrément de longueur décollée (m)
  30. c   incremsp : incrément d'extraction de fibre (m)
  31. c   fel : force élastique imposée (MN)
  32. c   ld1 : longueur décollée coté 1 (m)
  33. c   ld2 : longueur décollée coté 2 (m)
  34. c   sp1 : extraction du coté le moins ancré de la fibre (m)
  35. c   fcrit1 : force de début de décollement coté moins ancré (MN)
  36. c   fcrit2 : force de début de décollement coté plus ancré (MN)
  37. c   L1 : longueur ancrée variable par ecaill jusqu'au pic de force coté le moins ancré (m)
  38. c   L2 : longueur ancrée variable par ecaill jusqu'au pic de force coté le plus ancré (m)
  39. c   dls1 : incrément d'écaillage coté le moins ancré  (m)
  40. c   dls2 : incrément d'écaillage coté le plus ancré  (m)
  41. c   ls1 : longueur d'écaillage coté le moins ancré (m)
  42. c   ls2 : longueur d'écaillage coté le plus ancré (m)  
  43. c   F : force pour une fibre droite (MN)
  44. c   w1 : glissement de décollement ou élastique lié au coté 1 (m) 
  45. c   w2 : glissement de décollement ou élastique lié au coté 2 (m) 
  46. c   Ef : module d'young de la fibre (MPa)
  47. c   Df : diamètre de la fibre (m)
  48. c   Af : section de la fibre (m^2)
  49. c   Hf : module de rigidité de l'interface (MPa/m)
  50. c   kc : coefficient de diffusion de l'equation differentielle des fibres (1/m)
  51. c   Lf : longueur totale de la fibre (m)
  52. c   Tmax : contrainte admissible maximale par l'interface fibre-matrice (MPa)
  53. c   Td : contrainte de frottement fibre-matrice après décohésion (MPa)
  54. c   sk : glissement caractéristique, pilote la pente de l'extraction (m)
  55. c   F0 : force d'about de fibre (MN)
  56. c   alpha : angle d'ouverture du tétraèdre d'écaillage et du cone de rupture du béton en rad
  57. c   Rt : résistance à la traction de la matrice (MPa)
  58. c   mu : coefficient de frottement fibre matrice pour fibres orientées
  59. c   phi : angle d'inclinaison des fibres // normale à la fissure en rad
  60. c   Fd1 : force de décollement côté le moins ancré (MN)
  61. c   sd1 : glissement de décollement côté le moins ancré (m)
  62. c   sd2 : glissement de décollement côté le plus ancré (m)
  63. c   La1 : ancrage initial côté le moins ancré (m)
  64. c   La2 : ancrage initial côté le plus ancré (m)
  65. c   lsf1 : longueur d'écaillage sauvegardé côté le moins ancré (m)
  66. c   lsf2 : longueur d'écaillage sauvegardé côté le plus ancré (m)
  67. c   sel1 : glissement élastique côté le moins ancré (m)
  68. c   sel2 : glissement élastique côté le plus ancré (m)
  69. c   ldind1 : sert à indiquer si un décollement s'est déjà produit sur le côté le moins ancré (m)
  70. c   ldind2 : sert à indiquer si un décollement s'est déjà produit sur le côté le plus ancré (m)
  71. c   phid : valeur angle inclinaison fibre // normale à la fissure en °
  72. c   alphad : valeur angle d'ouverture du tétraèdre d'écaillage et du cone de rupture du béton en °
  73. c   sdf1 : glissement de décollement sauvegardé côté le moins ancré (m)
  74. c   sdf2 : glissement de décollement sauvegardé côté le plus ancré (m)
  75. c   Ffo : force de la fibre orientée (MN)
  76. c   Lf1 : ancrage de la fibre sauvegardé côté le moins ancré (m)
  77. c   Lf2 : ancrage de la fibre sauvegardé côté le plus ancré (m)
  78. c   ldf1 : longueur décollée de la fibre sauvegardée côté le moins ancré (m)
  79. c   ldf2 : longueur décollée de la fibre sauvegardée côté le plus ancré (m)
  80. c   wf1 : ouverture de fissure totale sauvegardée (m)
  81. c   Ff : force sauvegardée avant orientation (MN)
  82. c   lsk1,lsk2 : paramètre servant à calculer la pente de la droite d'écaillage 
  83. c   Fco : force admissible sur la surface du cone de rupture final (MN)
  84. c   Fmax : force maximale admisssible par la fibre (MN)
  85. c   fy : limite élastique de la fibre (MPa)
  86. c   rupt : indicateur de rupture de la fibre ou du cone, ou d'extraction  totale passe à true si il y a rupture
  87. c   mi1,mi2 : positions des valeurs de la courbe f(w) d'une fibre encradrant la valeur de wlist lors du calcul de la force moy
  88. c   n : Numéro de ligne de la courbe force moyenne ouverture de fissure
  89. c   nangl +1 : nombre de discrétisation de l'angle
  90. c   nlongfi : nombre de discrétisation de la longueur ancrée des fibres
  91. c   nwlis : nombre d'ouvertures de fissures pour l'intrpoation pour calcul des forces moyennes
  92. c   nombre d'ouvertures de fissures lors du calul de Fmoy
  93. c   vark : valeurs de k lorsque l'on sort de la boucle d'équilibre d'écaillage
  94. c   nlign1 : valeur optimisée du nombre de lignes i pour le calcul des forces
  95. c   incremwlis : incrément d'ouverture de fissure pour le calcul des forces moyennes (m)
  96. c   precision1 : sert à ne pas avoir 0 là ou on veut pas
  97. c   npic : position de wlist pour laquelle on a la force moyenne au pic
  98. c   n0longfi : nombre de longueur de fibre pour calcul de la force moyenne ne prenant pas en compte les longueurs < 1mm
  99. c   M0 : module d'écrouissage liée à l'abrasion (MPa)
  100. c   L0 : longueur caractéristique prenant en compte l'influence de la longueur ancrée sur l'abrasion (m)
  101. c   somFL : somme des forces pour un angle et tous les ancrages pour calcul fmoy 
  102. c   ffo1 : tableau des forces orientées dépendant de i,phi,L (MN)
  103. c   ffo2 : tableau des forces orientées dépendant de n,phi,L (MN)
  104. c   FmoyL : Force moyenne pour une orientation (MN)
  105. c   wlist : liste d'ouverture de fissure de la force moyenne (m)
  106. c   Dwf : différence entre l'ouverture de wlist et de wf1 servant a intrpoer
  107. c   fu : contrainte ultime pour les fibres (MPa)
  108. c   wint0,wint1 : valeurs de wf1 qui encadrent la valeur de wlist pour intrpoation (m)
  109. c   fint0,fint1 : valeurs de ffo1 qui correspondent à wint0 et wint1 pour intrpoation (MN)
  110. c   fres : valeur intrpoée de la force pour calcul de la force moyenne (MN)
  111. c   k0 : rigidité initiale pour une fibre servant a calculer s0 (MN/m)
  112. c   fk0 : valeur de force au pas 2 servant a calculer k0 (MN)
  113. c   sk0 : glissement correspondant à fk0 (m)
  114. c   Fpic : force moyenne au pic pour chaque cas (MN)
  115. c   k0m : rigidité initiale moyenne au pic (MN/m) 
  116. c   fcrit0 : force de début de décollement pour la longueur ancrée initiale pour optimisation du nombre de pas phase élastique (MN)
  117. c   nloc : localisation de la force moyenne maximale
  118. c   Fmoylis : liste de force moyenne 1d servant à localiser le pic (MN)
  119. c   wpic : ouverture de fissure au pic de force moyenne (m)
  120. c   wfin : ouverture de fissure pour laquelle Fmoy vaut 0 pour la première fois après le 1er pas (m)
  121. c   spind : indicateur de mise en extraction
  122. c   rotule : mettre a false pour ne pas prendre en compte la rotule d'extraction des fibres orientées
  123. c   finf : indicateur marquant wfin
  124. c   p : numéro de la ligne des listes (Td,phid,point caractéristique des courbes)
  125. c   lisTd : Liste des Td pour lesquels les forces moyennes ont été calculées (MPa)
  126. c   lisphid : Liste des phid pour lesquels les forces moyennes ont été calculées (MPa)
  127. c   lisfp : Liste des forces moyennes au pic correspondant a Td phid (MN)
  128. c   lisk0m : Liste des rigidités initiales moyennes correspondant a Td phid (MN/m)
  129. c   liswf : Liste des ouvertures de fissures moyennes finales (Fmoy=0) correspondant a Td phid (m)
  130. c   lisw03 : Liste des ouvertures de fissures moyennes pour F=0,3*fpic correspondant a Td phid (m)
  131. c   f03 : 0.3*fpic (MN)
  132.  
  133.  
  134. c        parametres
  135.          integer i,j,k,l,t,mi1,mi2,n,nligne,niter,nangl,nlongfi,nwlis
  136.          integer vark,nlign1,npic,n0longfi,nrt,nangredu,iind
  137.          real*8 incremf,incremld,incremsp,Lf,incremwlis,pi,precision1 
  138.          real*8 phit,phicrit
  139.  
  140.          parameter (incremwlis=5.0d-6)
  141.          parameter (PI=3.141592653589793D0,precision1=1.0d-10)
  142.  
  143.            real*8 M0,L0
  144.            real*8 Ef,Df,Af,Hf,kc,Tmax,Td,sk,F0,alpha,Rt,mu,phi
  145.            real*8 La1,La2
  146.            real*8 phid,alphad
  147.            real*8 lsk1,lsk2
  148.            real*8 Fmax
  149.            real*8 somFL,fy,s0,nligner
  150.            real*8 Fr
  151.            real*8 fu,Fpasav
  152.            real*8 wint0,wint1,fint0,fint1,fres,k0,fk0,sk0,scrit
  153.            real*8 fcrit0,wpic       
  154.  
  155. C cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  156.            logical rupt,spind,rotule,finf,ldind1,ldind2,ecail1  
  157.            parameter(rotule=.true.,ecail1=.false.)
  158. C c        variables pour la reduction de modele
  159.            real*8 f03,w03,Fpic,k0m,wfin
  160.            integer p,d,e
  161. *
  162.        SEGMENT MWOR
  163.            real*8 fel(mligne),ld1(mligne,miter),ld2(mligne,miter)
  164.            real*8 sp1(mligne),fcrit1(mligne,miter),fcrit2(mligne,miter)        
  165.            real*8 L1(mligne,miter),L2(mligne,miter),dls1(mligne,miter)
  166.            real*8 dls2(mligne,miter),ls1(mligne,miter),ls2(mligne,miter)
  167.            real*8 F(mligne,miter),w1(mligne,miter),w2(mligne,miter)
  168.            real*8 Fd1(mligne,miter),sd1(mligne,miter),lsf1(mligne)
  169.            real*8 sel1(mligne,miter),sdf1(mligne)
  170.            real*8 Fo(mligne,miter),ldf2(mligne),Ffo(mligne)
  171.            real*8 lsf2(mligne),Lf2(mligne),sdf2(mligne),Fco(mligne)
  172.            real*8 sd2(mligne,miter),sel2(mligne,miter)
  173.            real*8 Dwf(mligne),Ff(mligne),Lf1(mligne),ldf1(mligne)   
  174.            real*8 ffo2(mligne,mangl1,mlongfi)
  175.            real*8 Ffo1(mligne,mangl1,mlongfi)
  176.            real*8 FmoyL(mligne,mangl1),wf1(mligne,mangl1,mlongfi)
  177.            real*8 Fmoylis(mligne),wlist(mligne)                                 
  178. *                                       
  179.           real*8 lisRt(mrtang1),lisphid(mrtang1)
  180.           real*8 lisfp(mrtang1),lisk0m(mrtang1)
  181.           real*8 liswf(mrtang1)
  182.           real*8 lisw03(mrtang1)                
  183. *               
  184.           real*8 lisrt1(mangnrt),lisphid1(mangnrt)
  185.           real*8 lisrt2(mangnan)
  186.           real*8 lisphid2(mangnan)  
  187.           real*8 liswp1(mangnrt),liswp2(mangnan)     
  188.           real*8 lisk01(mangnrt),lisk02(mangnan)                
  189.        ENDSEGMENT
  190. *
  191.  
  192. C          print*,Ef,"Ef"
  193. C          print*,Lf,"Lf"
  194. C          print*,Df,"df"
  195. C          print*,Hf,"hf"
  196. C          print*,Tmax,"tmax"
  197. C          print*,td,"td"
  198. C          print*,sk,"sk"
  199. C          print*,F0,"f0"
  200. C          print*,alphad,"alphad"
  201. C          print*,Rt,"Rt"
  202. C          print*,mu,"mu"
  203. C          print*,M0,"M0"
  204. C          print*,fu,"fu"
  205. C          print*,L0,"L0"
  206. C          print*,fy,"fy"
  207.  
  208.          nwlis=int(Lf/(incremwlis))      
  209.  
  210. c        initialisation de variables
  211.          ldf1(1)=0.d0
  212.          lsf1(1)=0.d0 
  213.          sp1(1)=0.d0
  214.          sdf1(1)=0.d0
  215.  
  216.  
  217.          ldf2(1)=0.d0
  218.          lsf2(1)=0.d0 
  219.          sdf2(1)=0.d0                
  220.  
  221.          ldind1=.false.
  222.          ldind2=.false.
  223.  
  224. c parametres increments 
  225.           incremld=Lf/200.d0         
  226.           incremsp=Lf/1000.d0                
  227.  
  228.  
  229. c section de la fibre
  230.         Af=(PI*Df**2)/4
  231. c coefficient de l exponentielle        
  232.         kc=sqrt(4.d0*Hf/(Ef*Df))
  233. c conversion de l'angle du cone en radian
  234.         alpha=PI*alphad/180
  235.  
  236. c     boucle sur les angles
  237.       do j=0,nangl
  238.       if(j.lt.nangredu) then
  239.       phid=phicrit*j/(nangredu-1)
  240.       else
  241. c      phid=phit*j/(nangl)
  242.       if(j.eq.nangredu) then
  243.       phid=phicrit
  244.       else
  245.       phid=phit*j/(nangl)
  246.       end if
  247.       end if
  248.  
  249.  
  250. C       if(j.eq.0) then
  251. C       ffo2=0.d0
  252. C       wlist=0.d0
  253. C       FmoyL=0.d0
  254. C       Fmoylis=0.d0
  255. C       Ffo1=0.d0
  256. C       wf1=0.d0
  257. C       wlist=0.d0
  258. C       end if
  259.  
  260.  
  261. c      print*,j,phid,nangredu
  262. c      read*
  263.  
  264. c      phid=15.d0*j
  265.       phi=PI*phid/180  
  266.       finf=.false.      
  267. c        longueurs des ancrages
  268. c       boucle sur les longueurs ancrées
  269. c on exclue les fibres d ancrage inf a x mm --> non
  270. c int(2.d0*x*nlongfi/Lf)+1
  271. c        n0longfi=int(2.d0*0.0d-3*nlongfi/Lf)+1
  272.         n0longfi=1        
  273. c        print*,n0longfi
  274.         do l=n0longfi,nlongfi
  275.         La1=(Lf/2.d0)*l/nlongfi
  276.         La2=Lf-La1
  277.         Lf1(1)=La1
  278.         Lf2(1)=La2                  
  279.         spind=.false.
  280.  
  281.         lsk1=La1/5.d0        
  282.         lsk2=La2/5.d0        
  283. c     ici on va optimiser le nombre de lignes pour chaque longueur ancrée
  284. c         calcul de la force critique pour L
  285.           fcrit0=Tmax*PI*Df*(exp(2.d0*kc*La1)-1.d0)/
  286.      #    (kc*(exp(2.d0*kc*La1)+1.d0))
  287.           incremf=fcrit0/10.d0
  288. c glissement de debut de decollement
  289.           scrit=Tmax/Hf
  290. c         calcul du nombre ligne nécessaires
  291.           nligner=fcrit0/incremf+la1/incremld+la1/incremsp
  292.           nlign1=int(nligner)
  293. c          print*,nlign1
  294.  
  295.           do i=1,nlign1
  296.  
  297.              if(i.eq.1) then
  298.              rupt=.false.
  299.              s0=0.d0
  300.              end if
  301.  
  302.              if(i.eq.1) then
  303.              fel(i)=0.d0
  304.              else
  305.              fel(i)=fel(i-1)+incremf
  306.              end if
  307.  
  308.              do k=1,niter
  309. c il suffit d'un seul ecaill plus grand que l'ancrage pour casser              
  310. c            longueur du coté le moins ancré             
  311. C                  if(i.eq.1) then
  312. C                      L1(i,k)=La1
  313. C                  else
  314. C                      if(k.eq.1) then
  315. C                      L1(i,k)=Lf1(i-1)
  316. C                      else
  317. C                        if ((Lf1(i-1).lt.dls1(i,k-1)).or.(Lf1(i-1) 
  318. C      #                 .eq.0.d0).or.(rupt)) then
  319. C                        L1(i,k)=0.d0
  320. C                        rupt=.true.
  321. C c                       exit
  322. C                        go to 10
  323. C                        else
  324. C                        L1(i,k)=Lf1(i-1)-dls1(i,k-1)
  325. C                        end if
  326. C                      end if
  327. C                  end if
  328. c            longueur du coté le mois ancré    
  329.  
  330.  
  331. c             print*,tmax ici ca marche          
  332.                  if(i.eq.1) then
  333.                      if(k.eq.1) then
  334.                      L1(i,k)=La1
  335.                      else
  336.                        if((La1.le.ls1(i,k-1)).or.(rupt)) then
  337.                        L1(i,k)=0.d0
  338.                        rupt=.true.
  339.                        vark=k
  340. c                       exit
  341.                        go to 10
  342.                        else
  343.                        L1(i,k)=La1-dls1(i,k-1)
  344.                        end if
  345.                      end if
  346.                  else
  347.                      if(k.eq.1) then
  348.                      L1(i,k)=Lf1(i-1)
  349.                      else
  350.                        if((Lf1(i-1).le.dls1(i,k-1)).or.(rupt)
  351.      #                 .or.(Lf1(i-1).eq.0.d0)) then
  352.                        L1(i,k)=0.d0
  353.                        rupt=.true.
  354.                        vark=k
  355. c                       exit
  356.                        go to 10
  357.                        else
  358.                        L1(i,k)=Lf1(i-1)-dls1(i,k-1)
  359.                        end if
  360.                      end if
  361.                  end if                 
  362. c            longueur du coté le plus ancré             
  363.                  if(i.eq.1) then
  364.                      if(k.eq.1) then
  365.                      L2(i,k)=La2
  366.                      else
  367.                        if((La2.le.ls2(i,k-1)).or.(rupt)) then
  368.                        L2(i,k)=0.d0
  369.                        rupt=.true.
  370.                        vark=k
  371. c                       exit
  372.                        go to 10
  373.                        else
  374.                        L2(i,k)=La2-dls2(i,k-1)
  375.                        end if
  376.                      end if
  377.                  else
  378.                      if(k.eq.1) then
  379.                      L2(i,k)=Lf2(i-1)
  380.                      else
  381.                        if((Lf2(i-1).le.dls2(i,k-1)).or.(rupt)
  382.      #                 .or.(Lf2(i-1).eq.0.d0)) then
  383.                        L2(i,k)=0.d0
  384.                        rupt=.true.
  385.                        vark=k
  386. c                       exit
  387.                        go to 10
  388.                        else
  389.                        L2(i,k)=Lf2(i-1)-dls2(i,k-1)
  390.                        end if
  391.                      end if
  392.                  end if
  393.  
  394. c               print*,tmax
  395.  
  396. c            calcul de la force critique coté moins ancré
  397.              fcrit1(i,k)=Tmax*PI*Df*(exp(2.d0*kc*L1(i,k))-1.d0)/
  398.      #       (kc*(exp(2.d0*kc*L1(i,k))+1.d0))
  399.  
  400.  
  401.  
  402.  
  403. c            calcul de la force critique coté plus ancré
  404.              fcrit2(i,k)=Tmax*PI*Df*(exp(2.d0*kc*L2(i,k))-1.d0)/
  405.      #       (kc*(exp(2.d0*kc*L2(i,k))+1.d0))
  406.  
  407.  
  408. c********************************************                 
  409. c            coté le moins ancré
  410. c********************************************
  411.              if(i.eq.1) then
  412.              ldind1=.false.
  413.              ldind2=.false.
  414.              end if
  415.  
  416.              if((fel(i).ge.fcrit1(i,k)).or.(ldind1)) then  
  417.              ldind1=.true.
  418. c            phase de décollement
  419. c            initialisation ld1
  420.                  if(i.eq.1) then 
  421.                  ld1(i,k)=0.d0
  422.                  else
  423.                      if(k.eq.1) then
  424.                      ld1(i,k)=ldf1(i-1)+incremld
  425.                      else
  426.                        if(ldf1(i-1)+incremld-dls1(i,k-1).le.0.d0) 
  427.      #                 then
  428. C                        print*,"decolllement  a 0",ldf1(i-1),dls1(i,k-1)
  429. C      #                  ,i,k
  430.                        ld1(i,k)=0.d0
  431.                        else
  432.                        ld1(i,k)=ldf1(i-1)+incremld-dls1(i,k-1)
  433.                        end if       
  434.                      end if                       
  435.                  end if
  436. c                 print*,ld1(i,k),i,k
  437. c            calcul de la force de decollement
  438.              Fd1(i,k)=PI*Df*Td*ld1(i,k)+Tmax*PI*Df*
  439.      #       (exp(2.d0*kc*L1(i,k))-exp(2.d0*kc*ld1(i,k)))/
  440.      #       (kc*(exp(2.d0*kc*L1(i,k))+exp(2.d0*kc*ld1(i,k))))
  441.  
  442. c            calcul du glissement de decollement
  443.              sd1(i,k)=Tmax/Hf+((PI*Df*Td*ld1(i,k)**2)/(2.d0*Ef*Af))+
  444.      #       (Tmax*PI*Df*ld1(i,k)/(kc*Af*Ef))*((exp(2.d0*kc*L1(i,k))-
  445.      #       exp(2.d0*kc*ld1(i,k)))/(exp(2.d0*kc*L1(i,k))+
  446.      #       exp(2.d0*kc*ld1(i,k))))
  447. c                      print*,sdf1(i-1),sd1(i,k),l1(i,k),ld1(i,k)     
  448.                  sp1(i)=0.d0  
  449. c                force glissement sauvegardés  
  450.                  if(i.eq.1) then 
  451.                  F(i,k)=0.d0
  452.                  w1(i,k)=0.d0
  453. c                initialisation de variables
  454.                  sdf1(i)=0.d0
  455.                  sd1(1,k)=0.d0 
  456.                  sp1(i)=0.d0    
  457.                  else 
  458. c                    extraction du coté le moins ancré
  459. c                    seulement si le glissement de décollement reduit ET Ld augmente
  460.                      if((Ld1(i,k).ge.L1(i,k)).or.
  461.      #               ((sd1(i,k)-sdf1(i-1).le.0.d0)
  462.      #                 .and.(ld1(i,k).gt.ldf1(i-1))).or.spind) then
  463.                      s0=f0/k0
  464.                      ld1(i,k)=ldf1(i-1)                     
  465.  
  466.                      if(.not.spind) then
  467.                      iind=i
  468.                      end if
  469.  
  470.                      if((spind).and.(i.gt.iind)) then
  471.                      sp1(i)=sp1(i-1)+incremsp
  472.                      else
  473.                      sp1(i)=0.d0
  474.                      end if
  475.                      sd1(i,k)=sdf1(i-1)
  476.                      w1(i,k)=sdf1(i-1) 
  477.  
  478. c                    si on commence a extraire une fois on extrait toujours apres                     
  479.                      spind=.true.
  480.  
  481.                        if(sp1(i).gt.L1(i,k)) then 
  482.                        F(i,k)=0.d0
  483.                        sp1(i)=sp1(i-1)
  484.                        ld1(i,k)=ldf1(i-1)
  485.                        rupt=.true.
  486.                        else
  487.                        F(i,k)=(Td*(sk/(sk+sp1(i)))+M0*(sp1(i))                     
  488.      #                 *L1(i,k)/(L1(i,k)+L0)**2)*PI*Df*(L1(i,k)-
  489.      #                 sp1(i))+F0     
  490. c                       print*,f(i,k),s0,k0,sp1(i),ld1(i,k),sd1(i,k)
  491.                        ld1(i,k)=ldf1(i-1)
  492.                        end if
  493.                      else 
  494. c                    decollement  du coté le moins ancré   
  495.                      sp1(i)=0.d0               
  496.                      F(i,k)=Fd1(i,k)
  497.                      w1(i,k)=sd1(i,k)
  498.  
  499.                      end if
  500.  
  501.                  end if
  502.  
  503.              else 
  504. c            glissement elastique 
  505.                if(L1(i,k).eq.0.d0) then
  506.                rupt=.true.
  507.                else
  508.                sel1(i,k)=fel(i)*kc*(exp(kc*2.d0*L1(i,k))+1.d0)/
  509.      #         (Hf*PI*Df*(exp(2.d0*kc*L1(i,k))-1.d0))
  510.                end if
  511. c            la longueur decollee vaut 0
  512.              Ld1(i,k)=0.d0
  513.  
  514.  
  515. c            force glissement sauvegardés 
  516.                 if(i.eq.1) then
  517.                 F(i,k)=0.d0
  518.                 w1(i,k)=0.d0
  519.                 else
  520.                 F(i,k)=fel(i)
  521.                 w1(i,k)=sel1(i,k)
  522.                 end if
  523. c            initialisation de variables     
  524.                 sdf1(i)=0.d0
  525.                 sd1(i,k)=0.d0 
  526.                 sp1(i)=0.d0  
  527.  
  528.  
  529.              end if
  530.  
  531. c**********************************************                 
  532. c            coté le plus ancré
  533. c**********************************************
  534.              if((F(i,k).ge.fcrit2(i,k)).or.(ldind2)) then  
  535.              ldind2=.true.
  536. c            phase de décollement
  537. c            initialisation ld2
  538.                  if(i.eq.1) then 
  539.                  ld2(i,k)=0.d0
  540.                  else
  541. c                si la force décroit la longueur décollée ld2 stagne
  542. c                sauf si les deux longueurs sont égales
  543. c                si la force croit à la fin du décollement le décollement s arrete quand meme
  544.                   if(((F(i,k)-Ff(i-1).gt.0.d0).or.(l.eq.nlongfi))
  545.      #            .and.(sp1(i).eq.0.d0)) then
  546.                      if(k.eq.1) then
  547.                      ld2(i,k)=ldf2(i-1)+incremld
  548.                      else
  549.                        if(ldf2(i-1)+incremld-dls2(i,k-1).le.0.d0) 
  550.      #                 then
  551.                        ld2(i,k)=0.d0
  552.                        else
  553.                        ld2(i,k)=ldf2(i-1)+incremld-dls2(i,k-1)
  554.                        end if       
  555.                      end if  
  556.                   else
  557.                   ld2(i,k)=ldf2(i-1)  
  558.                   end if                  
  559.                  end if
  560.  
  561. c            la force est toujours celle du coté 1
  562. c            calcul du glissement de decollement
  563.              sd2(i,k)=Tmax/Hf+((PI*Df*Td*ld2(i,k)**2.d0)/(2.d0*Ef*Af))+
  564.      #       (Tmax*PI*Df*ld2(i,k)/(kc*Af*Ef))*(exp(2.d0*kc*L2(i,k))-
  565.      #       exp(2.d0*kc*ld2(i,k)))/(exp(2.d0*kc*L2(i,k))+
  566.      #       exp(2.d0*kc*ld2(i,k)))
  567.  
  568. c                glissement sauvegardés  
  569.                  if(i.eq.1) then 
  570.                  w2(i,k)=0.d0
  571. c                initialisation de variables
  572.                  sdf2(i)=0.d0
  573.                  sd2(i,k)=0.d0                 
  574.                  else 
  575. c                decollement  du coté le plus ancré         
  576.                    if(spind) then
  577.                    sd2(i,k)=sdf2(i-1)
  578.                    w2(i,k)=sdf2(i-1)
  579.                    else                  
  580.                    w2(i,k)=sd2(i,k)
  581.                    end if
  582.                  end if
  583.  
  584.              else 
  585.  
  586. c            glissement elastique 
  587.              sel2(i,k)=F(i,k)*kc*(exp(kc*2.d0*L2(i,k))+1.d0)/
  588.      #       (Hf*PI*Df*(exp(2.d0*kc*L2(i,k))-1.d0))
  589. c            la longueur decollee vaut 0
  590.              Ld2(i,k)=0.d0
  591. c            force glissement sauvegardés 
  592.                   if(i.eq.1) then
  593.                   w2(i,k)=0.d0
  594.                   else
  595.                   w2(i,k)=sel2(i,k)
  596.                   end if
  597. c            initialisation de variables     
  598.              sdf2(i)=0.d0
  599.              sd2(i,k)=0.d0
  600.  
  601.              end if 
  602.  
  603.  
  604. ccccccccccccccccccccccccccccccccccccccccccccccc
  605. c rotule plastique pour les fibres orientées
  606. ccccccccccccccccccccccccccccccccccccccccccccc
  607. C             if(rotule.and.spind) then
  608. C               if(phi.gt.0.d0) then
  609. C             F(i,k)=F(i,k)+fy*Df**2/(6.d0*sin(phi))
  610. C               end if
  611. C             end if
  612.  
  613. C        if(rotule.and.spind) then
  614. C          if(i.gt.1) then
  615. C            if(k.gt.1) then
  616. C c           Fp=2.d0*fy*Df**3/(6.d0*(ls1(i,k-1)+w1(i,k)
  617. C c     #     +sp1(i)))
  618. C c       Fo(i,k)=(F(i,k)*(1.d0+mu*sin(phi/2.d0))+Fp)      
  619. C c     #  /(1.d0-mu*sin(phi/2.d0))  
  620. C         Mp=fy*Df**3/6.d0
  621. C         Fo(i,k)=(F(i,k)*L1(i,k)-6*Mp*mu)/((-3*sin(phi)*mu+cos(phi))
  622. C      #  *L1(i,k))
  623. C            else
  624. C c           Fp=2.d0*fy*Df**3/(6.d0*(lsf1(i-1)+w1(i,k)  
  625. C c     #     +sp1(i)))    
  626. C c       Fo(i,k)=(F(i,k)*(1.d0+mu*sin(phi/2.d0))+Fp)      
  627. C c     #  /(1.d0-mu*sin(phi/2.d0))  
  628. C         Mp=fy*Df**3/6.d0
  629. C         Fo(i,k)=(F(i,k)*L1(i,k)-6*Mp*mu)/((-3*sin(phi)*mu+cos(phi))
  630. C      #  *L1(i,k))     
  631. C            end if
  632. C          else 
  633. C          fo(i,k)=0.d0
  634. C          endif
  635. C        else       
  636. C c       Fo(i,k)=F(i,k)*(1.d0+mu*sin(phi/2.d0))/(1.d0-mu*sin(phi/2.d0)) 
  637. C         Mp=fy*Df**3/6.d0
  638. C         Fo(i,k)=(F(i,k)*L1(i,k)-6*Mp*mu)/((-3*sin(phi)*mu+cos(phi))
  639. C      #  *L1(i,k))         
  640. C        end if
  641. c       Fr=sin(phi/2.d0)*(F(i,k)+Fo(i,k))
  642.  
  643.             if(rotule) then
  644.               if(phid.gt.5.d0) then
  645.             if(w1(i,k).le.scrit) then
  646.             D1=w1(i,k)/scrit
  647.             else
  648.             D1=1.d0
  649.             end if
  650.  
  651.             if(l1(i,k).eq.0.d0) then 
  652.             F(i,k)=0.d0
  653.             else
  654.             F(i,k)=F(i,k)+D1*3.d0*fy*Df**3.d0/(6.d0*sin(phi)*l1(i,k))
  655.             end if
  656.               end if
  657.             end if
  658.  
  659.  
  660.          Fo(i,k)=F(i,k)/(1.d0-mu*sin(phi)) 
  661.          Fr=Fo(i,k)*sin(phi)
  662. ccccccccccccccccccccccccccccccccccccccccccccccc
  663. c ecaill ccccccccccccccccccccccccccccccccccc
  664. ccccccccccccccccccccccccccccccccccccccccccccccc             
  665. c dans l ecaill on rentre la force radiale 
  666. c on ne fait pas de calcul d ecaill si ce n est pas necessaire
  667.  
  668. c mise a zero de ff(i-1)       
  669.       if(i.eq.1) then
  670.       Fpasav=0.d0
  671.       else
  672.       Fpasav=Ff(i-1)
  673.       end if
  674.  
  675.       if((Fo(i,k)-Fpasav.gt.0.d0).and.(.not.rupt)) then
  676.  
  677.       if((L1(i,k).eq.0.d0).and.(i.gt.1)) then     
  678.       ls1(i,k)=lsf1(i-1)
  679.       else
  680.         if(L1(i,k).eq.0.d0) then
  681.         ls1(i,k)=0.d0
  682.         else
  683.         call ecaill(i,k,phi,lsk1,alpha,L1(i,k),ls1(i,k),Df,nlign1,
  684.      #  Rt,Fr,j) 
  685.         end if     
  686.       end if
  687.  
  688.       if((L2(i,k).eq.0.d0).and.(i.gt.1)) then 
  689.       ls2(i,k)=lsf2(i-1)
  690.       else
  691.         if((L2(i,k).eq.0.d0).or.(ecail1)) then
  692.         ls2(i,k)=0.d0
  693.         else
  694.         call ecaill(i,k,phi,lsk2,alpha,L2(i,k),ls2(i,k),Df,nlign1,
  695.      #  Rt,Fr,j) 
  696.         end if     
  697.       end if
  698.  
  699.       else   
  700.       if(i.eq.1) then
  701.         ls1(i,k)=0.d0
  702.         ls2(i,k)=0.d0
  703.         else
  704.         ls1(i,k)=lsf1(i-1)
  705.         ls2(i,k)=lsf2(i-1)
  706.         end if
  707.       end if
  708.  
  709. c        la matrice coté 1 ne s'écaille plus si la force décroit
  710.                  if(i.gt.1) then                 
  711.                    if((Fo(i,k)-Fpasav.lt.0.d0).or.(Fo(i,k).eq.0.d0)
  712.      #             .or.(ls1(i,k).lt.lsf1(i-1)).or.rupt)
  713.      #             then
  714.                    ls1(i,k)=lsf1(i-1)
  715.                    end if
  716.                  else
  717.                  ls1(i,k)=0.d0
  718.                  end if
  719.  
  720. c        la matrice coté 2 ne s'écaille plus si la force décroit
  721.                  if(i.gt.1) then                 
  722.                    if((Fo(i,k)-Fpasav.lt.0.d0).or.(Fo(i,k).eq.0.d0)
  723.      #             .or.(ls2(i,k).lt.lsf2(i-1)).or.rupt)
  724.      #             then
  725.                    ls2(i,k)=lsf2(i-1)
  726.                    end if
  727.                  else
  728.                  ls2(i,k)=0.d0
  729.                  end if                 
  730.  
  731.  
  732. c            coté le moins ancré
  733.              if(i.eq.1) then
  734.              dls1(i,k)=ls1(i,k)
  735.              else
  736.                  if(ls1(i,k)-lsf1(i-1).le.0.d0) then
  737.                  dls1(i,k)=0.d0
  738.                  else
  739.                  dls1(i,k)=ls1(i,k)-lsf1(i-1)
  740.                  end if
  741.              end if
  742.  
  743. c            coté le plus ancré 
  744.              if(i.eq.1) then
  745.              dls2(i,k)=ls2(i,k)
  746.              else
  747.                  if(ls2(i,k)-lsf2(i-1).le.0.d0) then
  748.                  dls2(i,k)=0.d0
  749.                  else
  750.                  dls2(i,k)=ls2(i,k)-lsf2(i-1)
  751.                  end if
  752.              end if
  753.  
  754. C                  if(i.gt.90) then
  755. C                  if((j.eq.7).and.(t.eq.5).and.(l.eq.50)) then
  756. C             print*,ls1(i,k),ls2(i,k),i,k,fo(i,k),Fpasav,fo(i,k)-Fpasav,
  757. C      #dls1(i,k),L1(i,k),Fr
  758. C                  read*
  759. C                  end if 
  760. C                  end if                 
  761.  
  762.  
  763. c         if(isnan(ld1(i,k))) then
  764. c         print*,"ld1",i,ldf1(i-1),ldind1,fel(i),L1(i,k)
  765. C          else if(isnan(ld2(i,k))) then
  766. C           print*,"ld2"        
  767. C          else if(isnan(fcrit1(i,k))) then
  768. C           print*,"fcrit1"          
  769. C          else if(isnan(fcrit2(i,k))) then
  770. C           print*,"fcrit2"          
  771. C          else if(isnan(L1(i,k))) then
  772. C           print*,"L1"  
  773. C          else if(isnan(L2(i,k))) then
  774. C           print*,"L2"  
  775. C          else if(isnan(dls1(i,k))) then
  776. C           print*,"dls1"  
  777. C          else if(isnan(dls2(i,k))) then
  778. C           print*,"dls2" 
  779. C          else if(isnan(ls1(i,k))) then
  780. C           print*,"ls1" 
  781. C          else if(isnan(ls2(i,k))) then
  782. C           print*,"ls2" 
  783. C          else if(isnan(F(i,k))) then
  784. C           print*,"F"
  785. C          else if(isnan(w1(i,k))) then
  786. C           print*,"w1"
  787. C          else if(isnan(w2(i,k))) then
  788. C           print*,"w2"          
  789. C          else if(isnan(Fd1(i,k))) then
  790. C           print*,"Fd1"  
  791. C          else if(isnan(sd1(i,k))) then
  792. C           print*,"sd1"  
  793. C          else if(isnan(sel1(i,k))) then
  794. C           print*,"sel1" 
  795. C          else if(isnan(Fo(i,k))) then
  796. C           print*,"Fo" 
  797. C          else if(isnan(sd2(i,k))) then
  798. C           print*,"sd2"          
  799. C          else if(isnan(sel2(i,k))) then
  800. C           print*,"sel2"
  801. c         end if          
  802.  
  803. C           print*,fel(i),ld1(i,k),fcrit1(i,k),fcrit2(i,k),L1(i,k),L2(i,k)
  804. C      #,sp1(i),dls1(i,k),dls2(i,k),ls1(i,k),ls2(i,k),sd1(i,k)
  805.  
  806. C           print*,fcrit1(i,k),L1(i,k),sp1(i),dls1(i,k),sd1(i,k),Ld1(i,k),
  807. C      # F(i,k),ldf1(i-1),i,k,rupt
  808.  
  809. c test iter k
  810. c posiibilité d enlever le if crit1=100?
  811.              if(k.gt.1) then
  812.              crit1=abs((dls1(i,k)-dls1(i,k-1)))/la1
  813.              else
  814.              crit1=100.d0
  815.              end if
  816. c             print*,dls1(i,k),i,k,j,l,crit1
  817. c             read*
  818.             if(((dls1(i,k).eq.0.d0).or.(crit1.lt.1.0d-6))
  819.      #       .and.(k.gt.1)) then   
  820.             vark=k            
  821. c            exit
  822.             go to 10
  823.             else
  824.             vark=niter
  825.             end if
  826. c fin boucle sur k        
  827.             end do 
  828.  
  829.    10   continue
  830. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  831. c           sauvegarde des variables après itération sur k 
  832. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  833. c           force finale orientée
  834.             if(rupt) then
  835.              Ff(i)=0.d0
  836.             ffo(i)=0.d0
  837.             else            
  838.             Ff(i)=Fo(i,vark)
  839.             ffo(i)=fo(i,vark)
  840.             end if
  841.  
  842. c ls est recalculé avec la force actuelle donc il faut prendre celui d'avant 
  843. c mais toutes les autres valeurs sont bonnes car elles sont actualisées avec 
  844. c la valeur de dls(i,k-1)     
  845.              if(((l1(i,vark).eq.0.d0).or.(Ffo(i).eq.0.d0)).and.(i.gt.1)) 
  846.      #       then
  847.              lsf1(i)=lsf1(i-1)
  848.              lsf2(i)=lsf2(i-1)
  849.              else
  850.              lsf1(i)=ls1(i,vark-1) 
  851.              lsf2(i)=ls2(i,vark-1)             
  852.              end if
  853.  
  854.              if(lsf1(i).gt.la1) then
  855.              lf1(i)=0.d0
  856.              rupt=.true.
  857.              else 
  858.              Lf1(i)=La1-lsf1(i)             
  859.              end if   
  860.              if(lsf2(i).gt.la2) then
  861.              lf2(i)=0.d0
  862.              rupt=.true.
  863.              else 
  864.              Lf2(i)=La2-lsf2(i)             
  865.              end if                
  866.              Lf2(i)=La2-lsf2(i)
  867.              ldf1(i)=ld1(i,vark)
  868.              ldf2(i)=ld2(i,vark)                                                                        
  869.              sdf1(i)=sd1(i,vark)
  870.              sdf2(i)=sd2(i,vark)
  871. C c            la matrice ne s'écaille plus si la force décroit (lsf est constant)
  872.  
  873. c****************************************************
  874. c*********** Rupture du cone de matrice**************
  875. c****************************************************
  876.  
  877. c    force admissible par le cone
  878.       alphaco=30.d0*pi/180.d0
  879. c on ne compte pas l ecaill dans le cone
  880.        Fco(i)=PI*(Lf1(i)+lsf1(i)-sp1(i))*tan(alphaco)*(Df+(Lf1(i)
  881.      #+lsf1(i)-sp1(i)))*Rt
  882.          if((Ffo(i).ge.Fco(i)).or.(rupt)) then
  883. c         print*,'cone'
  884.          Ffo(i)=0.d0
  885.          rupt=.true.
  886.          end if
  887.  
  888. c****************************************************
  889. c*********** Rupture de la fibre*********************
  890. c****************************************************
  891.  
  892. c    force admissible par la fibre
  893.          Fmax=(PI*Df**2/4.d0)*fu
  894.          if((Ffo(i).ge.Fmax).or.(rupt)) then
  895.          Ffo(i)=0.d0
  896.          rupt=.true.
  897.          end if          
  898.  
  899. c on ajoute la def des fibres sur la longueurs écaillée pour palier le fait que le glissement
  900. c de décollement diminue
  901.              if(i.gt.1) then
  902.                if(rupt) then
  903.                wf1(i,j+1,l)=wf1(i-1,j+1,l)+incremsp
  904.                else               
  905.              wf1(i,j+1,l)=w1(i,vark)+w2(i,vark)+(lsf1(i)+lsf2(i))*
  906.      #        (1.d0-cos(phi))+sp1(i)+s0
  907.      #       +ffo(i)*(lsf1(i)+lsf2(i))/(Af*Ef)
  908.                end if
  909.               if(wf1(i,j+1,l).lt.1.0d-10) then
  910.               wf1(i,j+1,l)=0.d0
  911.               end if
  912.              else 
  913.              wf1(i,j+1,l)=0.d0
  914.              end if
  915.  
  916. c sauvegarde de la rigidité initiale pour une fibre de longueur ancrée connue
  917.               if((i.eq.2).and.(wf1(i,j+1,l).gt.0.d0)) then
  918.               sk0=wf1(i,j+1,l)
  919.               fk0=fel(i)
  920.               k0=fk0/sk0
  921.  
  922.                 if(sk0.eq.0.d0) then
  923.                 print*,"rigidite initiale infinie",sk0,fk0  
  924.                 read*        
  925.                 end if              
  926.               end if
  927.  
  928. c           tau(i)= Ffo(i)/(PI*Df*(Lf1(i)-sp1(i)))              
  929. c          ici ffo et tau contiennent les valeurs pour un angle une longueur ancrée 
  930. c          et toutes les i ouvertures de fissures
  931. c          force unitaire pour chacun des parametres            
  932.            Ffo1(i,j+1,l)=ffo(i)
  933.  
  934.  
  935.  
  936. c         affichage des courbes pour 5 longueurs ancrée            
  937. C            if((j.eq.7).and.(t.eq.5)) then
  938. C           if(l.eq.15) then
  939. C           open(unit=1,file="ffol0.res")
  940. C           write(unit=1,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
  941. C           else if (l.eq.20) then
  942. C           open(unit=1,file="ffol1.res")
  943. C           write(unit=1,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
  944. C           else if (l.eq.30) then
  945. C           open(unit=1,file="ffol2.res")
  946. C           write(unit=1,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
  947. C           else if (l.eq.40) then
  948. C           open(unit=1,file="ffol3.res")
  949. C           write(unit=1,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))    
  950. C           else if (l.eq.50) then
  951. C           open(unit=1,file="ffol4.res")
  952. C           write(unit=1,fmt='(11g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l)),
  953. C      #w1(i,vark),w2(i,vark),lsf1(i),lsf2(i),sp1(i),s0,ffo(i),rupt
  954. C           end if
  955. C            endif
  956.  
  957. c plot 'ffol0.res','ffol1.res','ffol2.res','ffol3.res','ffol4.res'
  958.  
  959. c enregistrement des forces ouvertures de fissures pour plusieurs angles
  960. C           if(l.eq.nlongfi) then
  961. C           if(t.eq.nrt) then
  962. C           if((j.eq.0)) then
  963. C           open(unit=11,file="flee0.res")
  964. C           write(unit=11,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
  965. C           else if (j.eq.1) then
  966. C           open(unit=11,file="flee1.res")
  967. C           write(unit=11,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
  968. C           else if (j.eq.3) then
  969. C           open(unit=11,file="flee2.res")
  970. C           write(unit=11,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
  971. C           else if (j.eq.4) then
  972. C           open(unit=11,file="flee3.res")
  973. C           write(unit=11,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
  974. C           else if (j.eq.6) then
  975. C           open(unit=11,file="flee4.res")
  976. C           write(unit=11,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
  977. C           end if
  978. C           end if
  979. C           end if
  980.  
  981. c    plot 'flee0.res','flee1.res','flee2.res','flee3.res','flee4.res'
  982. c plot 'flee2.res','flee3.res'
  983.  
  984. C          if(isnan(ffo1(i,j+1,l))) then
  985. C          print*,"ffo1"
  986. C          else if(isnan(sp1(i))) then
  987. C          print*,"sp1"
  988. C         else if(isnan(lsf1(i))) then
  989. C         print*,"lsf1"
  990. C         else if(isnan(sdf1(i))) then
  991. C         print*,"sdf1"        
  992. C         else if(isnan(ldf2(i))) then
  993. C         print*,"ldf2"            
  994. C         else if(isnan(Ffo(i))) then
  995. C         print*,"Ffo"   
  996. C C         else if(isnan(lsf2(i))) then
  997. C C         print*,"lsf2"           
  998. C C         else if(isnan(Lf2(i))) then
  999. C C         print*,"Lf2"            
  1000. C C         else if(isnan(sdf2(i))) then
  1001. C C         print*,"sdf2"        
  1002. C C         else if(isnan(Fco(i))) then
  1003. C C         print*,"Fco"          
  1004. C C         else if(isnan(Ff(i))) then
  1005. C C         print*,"Ff"
  1006. C         else if(isnan(Lf1(i))) then
  1007. C         print*,"Lf1"
  1008. C C         else if(isnan(ldf1(i))) then
  1009. C C         print*,"ldf1"
  1010. C C         else if(isnan(ffo2(i,j+1,l))) then
  1011. C C         print*,"ffo2"
  1012. C C         else if(isnan(Ffo1(i,j+1,l))) then
  1013. C C         print*,"Ffo1"
  1014. C C         else if(isnan(wf1(i,j+1,l))) then
  1015. C C         print*,"wf1" 
  1016. C          end if
  1017.  
  1018.  
  1019.  
  1020.  
  1021.  
  1022.  
  1023. C c fin boucle sur i          
  1024.            end do     
  1025.       do n=1,nwlis 
  1026. c      construit une liste d'ouverture de fissure     
  1027.        if(n.eq.1) then
  1028.        wlist(1)=0.d0
  1029.        else 
  1030.        wlist(n)=wlist(n-1)+incremwlis
  1031.        end if   
  1032. c calcule le minimum du carré de la diff entre les listes d ouverture de fissure
  1033. c on localise la valeur la plus proche de wlist(n)      
  1034.  
  1035.          do i=1,nlign1
  1036.          Dwf(i)=wlist(n)-wf1(i,j+1,l)
  1037.          if(dwf(i).lt.0.d0) then
  1038.          mi2=i
  1039.          mi1=i-1
  1040. c         exit
  1041.          go to 20
  1042.          end if
  1043.          end do
  1044.    20   continue       
  1045. c      les valeurs de wf1 qui encadrent wlis sont
  1046.        wint0=wf1(mi1,j+1,l)
  1047.        wint1=wf1(mi2,j+1,l)
  1048. c      les valeurs de ffo1 qui correspondent a wint0 et wint1 sont 
  1049.        fint0=ffo1(mi1,j+1,l) 
  1050.        fint1=ffo1(mi2,j+1,l)   
  1051.  
  1052.        call intrpo(wint0,wint1,fint0,fint1,wlist(n),fres)    
  1053.  
  1054.  
  1055.  
  1056.        if((wint0.eq.wint1).or.(fint0.eq.fint1).or.
  1057.      #  (fres.lt.0.d0)) then
  1058.        ffo2(n,j+1,l)=0.d0
  1059.        else
  1060.        ffo2(n,j+1,l)=fres
  1061.        end if
  1062.  
  1063. C           if(t.eq.4) then
  1064. C           if(j.eq.0) then
  1065. C c          print*,n,fres,l,wint0,wint1,fint0,fint1
  1066. C            print*,n,ffo2(n,j+1,l)
  1067. C           read*
  1068. C           end if
  1069. C           end if   
  1070.  
  1071. c affichage de la force orientée pour 10 angles
  1072. C          if((l.eq.2).and.(t.eq.1)) then
  1073. C           if(j.eq.0) then
  1074. C           open(unit=22,file="fphi0.res")
  1075. C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1076. C           else if (j.eq.1) then
  1077. C           open(unit=22,file="fphi1.res")
  1078. C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1079. C           else if (j.eq.2) then
  1080. C           open(unit=22,file="fphi2.res")
  1081. C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1082. C           else if (j.eq.3) then
  1083. C           open(unit=22,file="fphi3.res")
  1084. C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1085. C           else if (j.eq.4) then
  1086. C           open(unit=22,file="fphi4.res")
  1087. C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1088. C           else if (j.eq.5) then
  1089. C           open(unit=22,file="fphi5.res")
  1090. C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1091. C           else if (j.eq.6) then
  1092. C           open(unit=22,file="fphi6.res")
  1093. C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1094. C           else if (j.eq.7) then
  1095. C           open(unit=22,file="fphi7.res")
  1096. C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l)) 
  1097. C           else if (j.eq.8) then
  1098. C           open(unit=22,file="fphi8.res")
  1099. C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1100. C           else if (j.eq.9) then
  1101. C           open(unit=22,file="fphi9.res")
  1102. C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1103. C           else if (j.eq.10) then
  1104. C           open(unit=22,file="fphi10.res")
  1105. C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1106. C           end if
  1107. C           end if
  1108.  
  1109. c         plot 'fphi0.res','fphi1.res','fphi2.res','fphi3.res','fphi4.res','fphi5.res','fphi6.res','fphi7.res','fphi8.res','fphi9.res','fphi10.res'
  1110. c plot 'fphi6.res'
  1111.  
  1112. C            if((j.eq.0).and.(t.eq.3)) then
  1113. C           if(l.eq.5) then
  1114. C           open(unit=17,file="fmoyphi0.res")
  1115. C           write(unit=17,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1116. C           else if (l.eq.10) then
  1117. C           open(unit=17,file="fmoyphi1.res")
  1118. C           write(unit=17,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1119. C           else if (l.eq.15) then
  1120. C           open(unit=17,file="fmoyphi2.res")
  1121. C           write(unit=17,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1122. C           else if (l.eq.30) then
  1123. C           open(unit=17,file="fmoyphi3.res")
  1124. C           write(unit=17,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
  1125. C           else if (l.eq.40) then
  1126. C           open(unit=17,file="fmoyphi4.res")
  1127. C           write(unit=17,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l)) 
  1128. C           end if        
  1129. C            end if
  1130. c   plot 'fmoyphi0.res','fmoyphi1.res','fmoyphi2.res','fmoyphi3.res','fmoyphi4.res'
  1131.  
  1132.        end do 
  1133. c       fin boucle n
  1134. c va chercher la force correspondant a la valeur de wlist
  1135. c  fin boucle sur l
  1136.         end do 
  1137.  
  1138. C         if(j.eq.5) then
  1139. C         print*,sum(test)/(nlongfi-n0longfi+1) ,phid,Lf/4*cos(phi)**2
  1140. C         read*
  1141. C         end if
  1142.  
  1143. c       calcul de la force moyenne ancrage variable pour chaque angle
  1144.           do n=1,nwlis
  1145.              do l=n0longfi,nlongfi
  1146.              if (l.eq.n0longfi) then
  1147.              somFL=ffo2(n,j+1,n0longfi)
  1148.              else 
  1149.              somFL=ffo2(n,j+1,l)+somFL
  1150.              endif
  1151.              end do 
  1152. c            fin de boucle sur l
  1153.  
  1154.           if(n.eq.1) then            
  1155.           FmoyL(1,j+1)=0.d0 
  1156.           else
  1157.           FmoyL(n,j+1)=somFL/(nlongfi-n0longfi+1) 
  1158.           end if
  1159.           Fmoylis(n)=FmoyL(n,j+1)
  1160.  
  1161. C           print*,Fmoylis(n),ffo2(n,j+1,l),wf1(n,j+1,l),n,l,t,nwlis,
  1162. C      #wlist(n)
  1163.  
  1164. c affichage de la force moyenne pour 10 angles
  1165. C           if(t.eq.4) then
  1166. C           if(j.eq.0) then
  1167. C           open(unit=19,file="fmoyphii0.res")
  1168. C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))
  1169. C           else if (j.eq.1) then
  1170. C           open(unit=19,file="fmoyphii1.res")
  1171. C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))
  1172. C           else if (j.eq.2) then
  1173. C           open(unit=19,file="fmoyphii2.res")
  1174. C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
  1175. C           else if (j.eq.3) then
  1176. C           open(unit=19,file="fmoyphii3.res")
  1177. C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))
  1178. C           else if (j.eq.4) then
  1179. C           open(unit=19,file="fmoyphii4.res")
  1180. C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
  1181. C           else if (j.eq.5) then
  1182. C           open(unit=19,file="fmoyphii5.res")
  1183. C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
  1184. C           else if (j.eq.6) then
  1185. C           open(unit=19,file="fmoyphii6.res")
  1186. C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
  1187. C           else if (j.eq.7) then
  1188. C           open(unit=19,file="fmoyphii7.res")
  1189. C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
  1190. C           else if (j.eq.8) then
  1191. C           open(unit=19,file="fmoyphii8.res")
  1192. C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))
  1193. C           else if (j.eq.9) then
  1194. C           open(unit=19,file="fmoyphii9.res")
  1195. C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
  1196. C           else if (j.eq.10) then
  1197. C           open(unit=19,file="fmoyphii10.res")
  1198. C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
  1199. C           end if
  1200. C           end if
  1201.  
  1202. c         plot 'fmoyphii0.res','fmoyphii1.res','fmoyphii2.res','fmoyphii3.res','fmoyphii4.res','fmoyphii5.res','fmoyphii6.res','fmoyphii7.res','fmoyphii8.res','fmoyphii9.res','fmoyphii10.res'
  1203.  
  1204.  
  1205. c  forces moyennes pour 6 angles et 5 contraintes laterales
  1206. C        if(t.eq.1) then
  1207. C          if(j.eq.0) then 
  1208. C          open(unit=19,file="fmoyt1j0.res")
  1209. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
  1210. C          elseif (j.eq.1) then
  1211. C          open(unit=19,file="fmoyt1j1.res")
  1212. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)          
  1213. C          elseif (j.eq.2) then
  1214. C          open(unit=19,file="fmoyt1j2.res")
  1215. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)   
  1216. C          elseif (j.eq.3) then
  1217. C          open(unit=19,file="fmoyt1j3.res")
  1218. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
  1219. C          elseif (j.eq.4) then
  1220. C          open(unit=19,file="fmoyt1j4.res")
  1221. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
  1222. C          elseif (j.eq.5) then
  1223. C          open(unit=19,file="fmoyt1j5.res")
  1224. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
  1225. C          end if
  1226. C        elseif(t.eq.2) then
  1227. C          if(j.eq.0) then 
  1228. C          open(unit=19,file="fmoyt2j0.res")
  1229. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
  1230. C          elseif (j.eq.1) then
  1231. C          open(unit=19,file="fmoyt2j1.res")
  1232. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)       
  1233. C          elseif (j.eq.2) then
  1234. C          open(unit=19,file="fmoyt2j2.res")
  1235. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
  1236. C          elseif (j.eq.3) then
  1237. C          open(unit=19,file="fmoyt2j3.res")
  1238. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)  
  1239. C          elseif (j.eq.4) then
  1240. C          open(unit=19,file="fmoyt2j4.res")
  1241. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
  1242. C          elseif (j.eq.5) then
  1243. C          open(unit=19,file="fmoyt2j5.res")
  1244. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
  1245. C          end if
  1246. C        elseif(t.eq.3) then
  1247. C          if(j.eq.0) then 
  1248. C          open(unit=19,file="fmoyt3j0.res")
  1249. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
  1250. C          elseif (j.eq.1) then
  1251. C          open(unit=19,file="fmoyt3j1.res")
  1252. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)         
  1253. C          elseif (j.eq.2) then
  1254. C          open(unit=19,file="fmoyt3j2.res")
  1255. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)  
  1256. C          elseif (j.eq.3) then
  1257. C          open(unit=19,file="fmoyt3j3.res")
  1258. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)  
  1259. C          elseif (j.eq.4) then
  1260. C          open(unit=19,file="fmoyt3j4.res")
  1261. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
  1262. C          elseif (j.eq.5) then
  1263. C          open(unit=19,file="fmoyt3j5.res")
  1264. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
  1265. C          end if
  1266. C        elseif(t.eq.4) then
  1267. C          if(j.eq.0) then 
  1268. C          open(unit=19,file="fmoyt4j0.res")
  1269. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
  1270. C          elseif (j.eq.1) then
  1271. C          open(unit=19,file="fmoyt4j1.res")
  1272. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)         
  1273. C          elseif (j.eq.2) then
  1274. C          open(unit=19,file="fmoyt4j2.res")
  1275. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)  
  1276. C          elseif (j.eq.3) then
  1277. C          open(unit=19,file="fmoyt4j3.res")
  1278. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
  1279. C          elseif (j.eq.4) then
  1280. C          open(unit=19,file="fmoyt4j4.res")
  1281. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
  1282. C          elseif (j.eq.5) then
  1283. C          open(unit=19,file="fmoyt4j5.res")
  1284. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
  1285. C          end if
  1286. C        elseif(t.eq.5) then
  1287. C          if(j.eq.0) then 
  1288. C          open(unit=19,file="fmoyt5j0.res")
  1289. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
  1290. C          elseif (j.eq.1) then
  1291. C          open(unit=19,file="fmoyt5j1.res")
  1292. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)          
  1293. C          elseif (j.eq.2) then
  1294. C          open(unit=19,file="fmoyt5j2.res")
  1295. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
  1296. C          elseif (j.eq.3) then
  1297. C          open(unit=19,file="fmoyt5j3.res")
  1298. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
  1299. C          elseif (j.eq.4) then
  1300. C          open(unit=19,file="fmoyt5j4.res")
  1301. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
  1302. C          elseif (j.eq.5) then
  1303. C          open(unit=19,file="fmoyt5j5.res")
  1304. C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
  1305. C          end if
  1306. C          end if   
  1307.  
  1308.  
  1309.  
  1310. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc              
  1311. c sauvergarde des paramètres pour la réduction du modèle
  1312. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  1313.  
  1314. c Rigidité initiale
  1315.        if(n.eq.2) then
  1316.        k0m=(FmoyL(2,j+1)-FmoyL(1,j+1))/(wlist(2)-wlist(1))
  1317.        end if
  1318.  
  1319. c force et ouverture fissure fin d extraction       
  1320.        if((FmoyL(n,j+1).eq.0.d0).and.(n.gt.1).and.(.not.finf)) then
  1321. c       Ffin=FmoyL(n-1,j+1)
  1322.        wfin=wlist(n-1)
  1323.        finf=.true.
  1324.        end if
  1325.  
  1326. C        open(unit=13,file="fmoyphi3d.res")
  1327. C        if(n.eq.1) then
  1328. C        write(unit=13,fmt='(11g13.5)')
  1329. C        end if
  1330. C        write(unit=13,fmt='(11g13.5)') wlist(n),phid,FmoyL(n,j+1)
  1331.  
  1332.  
  1333. c fin boucle sur n   
  1334.           end do
  1335.  
  1336. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc              
  1337. c sauvergarde des paramètres pour la réduction du modèle
  1338. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc   
  1339. c force et ouverture de fissure au pic
  1340. cc      picloc=maxloc(Fmoylis)      
  1341. cc       npic=int(picloc(1))
  1342.        call amxloc(Fmoylis,Fmoylis(/1),npic)
  1343.        Fpic=FmoyL(npic,j+1) 
  1344.        wpic=wlist(npic)
  1345.  
  1346.        F03=0.3d0*Fpic
  1347.  
  1348. c recherche de l'ouverture de fissure correspondant a f03 et f06       
  1349.        do n=npic,nwlis
  1350. C         if((fmoylis(n)-F06.le.0.d0).and.
  1351. C      #  (fmoylis(n-1)-F06.ge.0.d0)) then
  1352. C         w06=wlist(n)
  1353. C         end if
  1354.  
  1355.         if((fmoylis(n)-F03.le.0.d0).and.
  1356.      #  (fmoylis(n-1)-F03.ge.0.d0)) then
  1357.         w03=wlist(n)
  1358. c        exit
  1359.         go to 30
  1360.         end if
  1361.        end do
  1362.    30  continue       
  1363.  
  1364. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  1365. cccccccccccccccccc Déplacement total cccccccccccccccccccccc
  1366. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  1367. c on doit reconstruire la liste de déplacement après qu'on connaisse le pic de multifissuration
  1368. c pour gérer la localisation qui se fait au pic de force en fait donc est ce qu on a besoin de faire ca ici?
  1369. c et puis est ce qu on devrait pas faire la multifissuration qu a la fin?
  1370.  
  1371.  
  1372. c affichage de la force en 3d en fonction de w et angle
  1373. c         splot 'fmoyphi3d.res' u 1:2:3
  1374. C            open(unit=13,file="fmoyphi3d.res")
  1375. C            if(n.eq.0) then
  1376. C            write(unit=13,fmt='(11g13.5)')
  1377. C            end if
  1378. C            write(unit=13,fmt='(11g13.5)') wlist(n),phid,FmoyL(n,j+1)
  1379.  
  1380. c          flis(n)=FmoyL(n,j+1)
  1381. c          print*,flis            
  1382.  
  1383. c    2eme fin de boucle sur n pour calcul deplacement            
  1384. c            end do            
  1385.  
  1386. C            open(unit=7,file="kmoyphiredu3d.res")
  1387. C            if(j.eq.0) then
  1388. C            write(unit=7,fmt='(11g13.5)')
  1389. C            end if
  1390. C            write(unit=7,fmt='(11g13.5)') Rt,phid,k0m 
  1391. c liste de rigidité initiale
  1392.                 lisRt(p)=Rt 
  1393.                 lisphid(p)=phid 
  1394.                 lisk0m(p)=k0m               
  1395.  
  1396. C            open(unit=15,file="fpmoyphiredu3d.res")
  1397. C            if(j.eq.0) then
  1398. C            write(unit=15,fmt='(11g13.5)')
  1399. C            end if
  1400. C            write(unit=15,fmt='(11g13.5)') Rt,phid,Fpic 
  1401. c liste de force au pic
  1402.                 lisFp(p)= Fpic 
  1403.  
  1404. C            open(unit=18,file="wpmoyphiredu3d.res")
  1405. C            if(j.eq.0) then
  1406. C            write(unit=18,fmt='(11g13.5)')
  1407. C            end if
  1408. C            write(unit=18,fmt='(11g13.5)') Rt,phid,wpic
  1409. c liste d'ouverture au pic
  1410.                 if(j.lt.nangredu) then
  1411.                 lisRt1(d)=Rt 
  1412.                 lisphid1(d)=phid 
  1413.                 liswp1(d)= wpic 
  1414.  
  1415. C            open(unit=181,file="wp1moyphiredu3d.res")
  1416. C            if(j.eq.0) then
  1417. C            write(unit=181,fmt='(11g13.5)')
  1418. C            end if
  1419. C          write(unit=181,fmt='(11g13.5)') lisrt1(d),lisphid1(d),liswp1(d)    
  1420.                 else
  1421. c                listd2(1)=listd1(d-1)
  1422.                 lisRt2(e)=Rt 
  1423.                 lisphid2(e)=phid 
  1424. c                lisphid2(1)=lisphid1(d-1)
  1425.                 liswp2(e)=wpic
  1426. c                liswp2(1)=liswp1(d-1)
  1427.  
  1428. C            open(unit=182,file="wp2moyphiredu3d.res")
  1429. C            if(j.eq.nangredu) then
  1430. C            write(unit=182,fmt='(11g13.5)')
  1431. C            end if
  1432. C          write(unit=182,fmt='(11g13.5)') lisrt2(e),lisphid2(e),liswp2(e)                
  1433.                 end if          
  1434.  
  1435.                if(j.lt.nangredu) then
  1436.                lisk01(d)= k0m
  1437.                d=d+1 
  1438.                else
  1439.                lisk02(e)=k0m 
  1440.                 e=e+1                         
  1441.                end if 
  1442.  
  1443. C              print*, k0m,e,j,t,phid
  1444. C              read*             
  1445.  
  1446.  
  1447. C            open(unit=21,file="wfmoyphiredu3d.res")
  1448. C            if(j.eq.0) then
  1449. C            write(unit=21,fmt='(11g13.5)')
  1450. C            end if
  1451. C            write(unit=21,fmt='(11g13.5)') Rt,phid,wfin
  1452. c liste d'ouverture finale
  1453.                 liswf(p)= wfin
  1454.  
  1455. c liste d'ouverture a 30% de force 
  1456.  
  1457. C            open(unit=24,file="w03moyphiredu3d.res")
  1458. C            if(j.eq.0) then
  1459. C            write(unit=24,fmt='(11g13.5)')
  1460. C            end if
  1461. C            write(unit=24,fmt='(11g13.5)') Rt,phid,w03           
  1462.                 lisw03(p)=w03    
  1463.  
  1464.  
  1465. c p numéro de ligne pour la réduction   
  1466. c            print*,"Force moyenne ",p,"effectuee"       
  1467.             p=p+1            
  1468.  
  1469. c fin boucle sur j
  1470.         end do
  1471.  
  1472. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  1473. c fin partie prpfib      
  1474. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc            
  1475.  
  1476. C C         print*,ld2
  1477. C C         print*,ld1         
  1478.           return
  1479.  
  1480.  
  1481.        end 
  1482.  
  1483.  
  1484.  
  1485.  
  1486.  
  1487.  
  1488.  
  1489.  
  1490.  
  1491.  
  1492.  
  1493.  
  1494.  
  1495.  
  1496.  
  1497.  
  1498.  
  1499.  
  1500.  
  1501.  
  1502.  
  1503.  
  1504.  

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