Télécharger iniree.eso

Retour à la liste

Numérotation des lignes :

  1. C INIREE SOURCE PV 16/12/08 21:15:05 9245
  2. SUBROUTINE INIREE
  3. c*************************************************************
  4. c evalue xpetit et xgrand pour la machine
  5. c on veut pouvoir prendre l 'inverse de xpetit et xgrand
  6. c*************************************************************
  7. IMPLICIT INTEGER(I-N)
  8. real*8 xre1,xre2, xzauxi,xii
  9. real*4 xsre1,xsre2,xsre3,xsauxi
  10. common/cinire/xre1,xre2,xzauxi,xsre1,xsre2,xsre3,xsauxi
  11. -INC CCOPTIO
  12. -INC CCREEL
  13. -INC CCTRACE
  14. c calcul xpetit
  15. xre1 = 1.D0
  16. * nexpo ajuste pour ne pas produire de floating point exception sur x86_64
  17. do nexpo = 1,307
  18. xre2 = xre1/10.D0
  19. if (xre2 .eq. 0.D0) then
  20. if (xre2.gt.xre1) call gibtem(xii)
  21. goto 110
  22. endif
  23. xre1 = xre2
  24. enddo
  25. ** moterr(1:6) = 'XPETIT'
  26. ** call erreur(907)
  27. * write(6,*) ' pas marche ', nexpo,xre1
  28. ** return
  29.  
  30. 110 continue
  31. * write(6,*) 'nexpo' , nexpo, xre1
  32. xpetit = xre1 * 1.D3
  33. * write(6,*) 'xpetit', xpetit
  34.  
  35. c calcul xgrand
  36. xre1 = 1.D0
  37. do nexpo = 1,307
  38. xre2=xre1*10.D0
  39. if (.NOT.(xre2*10.D0 .gt. xre2)) then
  40. if (xre2.lt.xre1) call gibtem(xii)
  41. goto 120
  42. endif
  43. xre1 = xre2
  44. enddo
  45. * write(6,*) ' pas marche ', nexpo,xre1
  46. ** moterr(1:6) = 'XGRAND'
  47. ** call erreur(907)
  48. ** return
  49.  
  50. 120 continue
  51. * write(6,*) 'nexpo' , nexpo, xre1
  52. xgrand = xre1 / 1.D3
  53. * write(6,*) 'xgrand', xgrand
  54. c
  55. c pour pouvoir inverser
  56. xzauxi = min(abs(log(xgrand)),abs(log(xpetit)))
  57. * write(6,*) xzauxi, abs(log(xgrand)) , abs(log(xpetit))
  58. if (abs(xzauxi - abs(log(xgrand))) .lt. 0.5D0) then
  59. xpetit = 1.D0/xgrand
  60. * write(6,*) 'xpetit', xpetit, 1/xpetit, 1/xgrand
  61. else
  62. xgrand = 1./xpetit
  63. * write(6,*) 'xgrand ' , xgrand, 1/xpetit, 1/xgrand
  64. endif
  65.  
  66. c calcul xspeti
  67. xsre1 = 1.E0
  68. * nexpo ajuste pour ne pas produire de floating point exception sur x86_64
  69. do nexpo = 1,36
  70. xsre2 = xsre1/10.E0
  71. xsre3 = xsre2/10.E0
  72. * if (nexpo .gt. 40) write(6,*) (b/a)
  73. if (xsre2/xsre3 .gt. 9.999D0) then
  74. xsre1=xsre2
  75. else
  76. if (xsre2.gt.xsre1) call gibtem(xii)
  77. goto 130
  78. endif
  79. enddo
  80. goto 130
  81. * write(6,*) ' pas marche ', nexpo,xsre1
  82. ** moterr(1:6) = 'XSPETI'
  83. ** call erreur(907)
  84. ** return
  85.  
  86. 130 continue
  87. * write(6,*) 'nexpo' , nexpo, xsre1
  88. xspeti = xsre1 * 1E2
  89. * write(6,*) 'xspeti', xspeti
  90.  
  91. c calcul xsgran
  92. xsre1 = 1.E0
  93. do nexpo = 1,37
  94. xsre2 = xsre1*10.E0
  95. xsre3 = xsre2*10.E0
  96. * if (nexpo .gt. 30) write(6,*) a, (b.gt.a)
  97. if (xsre3 .gt. xsre2) then
  98. xsre1 = xsre2
  99. else
  100. if (xsre2.lt.xsre1) call gibtem(xii)
  101. goto 140
  102. endif
  103. enddo
  104. goto 140
  105. ** moterr(1:6) = 'XSGRAN'
  106. ** call erreur(907)
  107. * write(6,*) ' pas marche ', nexpo,xsre1
  108. ** return
  109.  
  110. 140 continue
  111. * write(6,*) 'nexpo' , nexpo, xsre1
  112. xsgran = xsre1/1.E3
  113. * write(6,*) 'xsgran', xsgran
  114. amplit=xsgran
  115. c pour pouvoir inverser
  116. xsauxi = min(abs(log(xsgran)),abs(log(xspeti)))
  117. * write(6,*) xsauxi, abs(log(xsgran)) , abs(log(xspeti))
  118. if (abs(xsauxi - abs(log(xsgran))).lt.0.5) then
  119. xspeti = 1.E0/xsgran
  120. * write(6,*) 'xspeti', xspeti, 1/xspeti, 1/xsgran
  121. else
  122. xsgran = 1.E0/xspeti
  123. * write(6,*) 'xsgran ' , xsgran, 1/xspeti, 1/xsgran
  124. endif
  125.  
  126. c calcul xzprec ,
  127. xre1 = 1.D0
  128. do nexpol = 1,10000
  129. xre2=1.d0+xre1/10.d0
  130. if (xre2.le.1.d0) then
  131. if (xre2.lt.xre1) call gibtem(xii)
  132. goto 150
  133. endif
  134. xre1 = xre2-1.d0
  135. enddo
  136. moterr(1:6) = 'XZPREC'
  137. call erreur(907)
  138. * write(6,*) ' pas marche ', nexpo,xre1
  139. return
  140.  
  141. 150 continue
  142. * write(6,*) 'nexpo' , nexpo, xre1
  143. xzprec = xre1
  144. * write(6,*) 'xzprec', xzprec
  145.  
  146. c calcul xszpre ,
  147. xsre1 = 1.E0
  148. do nexpo = 1,10000
  149. xsre2=1.E0+xsre1/10.E0
  150. if (xsre2.le.1.e0) then
  151. if (xsre2.lt.xsre1) call gibtem(xii)
  152. goto 160
  153. endif
  154. xsre1 = xsre2-1E0
  155. enddo
  156. moterr(1:6) = 'XSZPRE'
  157. call erreur(907)
  158. * write(6,*) ' pas marche ', nexpo,xsre1
  159. return
  160.  
  161. 160 continue
  162. * write(6,*) 'nexpo' , nexpo, xsre1
  163. xszpre = xsre1
  164. * write(6,*) 'xszpre', xszpre
  165. return
  166.  
  167. END
  168.  
  169.  
  170.  
  171.  

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