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

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