Télécharger iniree.eso

Retour à la liste

C INIREE    SOURCE    PV        16/12/08    21:15:05     9245           
      SUBROUTINE INIREE
c*************************************************************
c  evalue xpetit et xgrand pour la machine
c  on veut pouvoir prendre l 'inverse de xpetit et xgrand
c*************************************************************
      IMPLICIT INTEGER(I-N)
      real*8 xre1,xre2, xzauxi,xii
      real*4 xsre1,xsre2,xsre3,xsauxi
      common/cinire/xre1,xre2,xzauxi,xsre1,xsre2,xsre3,xsauxi
-INC CCOPTIO
-INC CCREEL
-INC CCTRACE
c calcul xpetit
      xre1 = 1.D0
*  nexpo ajuste pour ne pas produire de floating point exception sur x86_64
      do nexpo = 1,307
        xre2 = xre1/10.D0
        if (xre2 .eq. 0.D0) then
        if (xre2.gt.xre1) call gibtem(xii)
          goto 110
        endif
        xre1 = xre2
      enddo
**    moterr(1:6) = 'XPETIT'
**    call erreur(907)
*      write(6,*) ' pas marche ', nexpo,xre1
**    return
 
 110  continue
*     write(6,*) 'nexpo' , nexpo, xre1
      xpetit = xre1 * 1.D3
*     write(6,*) 'xpetit', xpetit
 
c calcul xgrand
      xre1 = 1.D0
      do nexpo = 1,307
        xre2=xre1*10.D0
        if (.NOT.(xre2*10.D0 .gt. xre2)) then
        if (xre2.lt.xre1) call gibtem(xii)
        goto 120
        endif
        xre1 = xre2
      enddo
*      write(6,*) ' pas marche ', nexpo,xre1
**    moterr(1:6) = 'XGRAND'
**    call erreur(907)
**    return
 
 120  continue
*      write(6,*) 'nexpo' , nexpo, xre1
      xgrand = xre1 / 1.D3
*      write(6,*) 'xgrand', xgrand
c
c  pour pouvoir inverser
      xzauxi = min(abs(log(xgrand)),abs(log(xpetit)))
*      write(6,*) xzauxi, abs(log(xgrand)) , abs(log(xpetit))
      if (abs(xzauxi - abs(log(xgrand))) .lt. 0.5D0) then
        xpetit = 1.D0/xgrand
*        write(6,*) 'xpetit', xpetit, 1/xpetit, 1/xgrand
      else
        xgrand = 1./xpetit
*        write(6,*) 'xgrand ' , xgrand, 1/xpetit, 1/xgrand
      endif
 
c calcul xspeti
      xsre1 = 1.E0
*  nexpo ajuste pour ne pas produire de floating point exception sur x86_64
      do nexpo = 1,36
          xsre2 = xsre1/10.E0
          xsre3 = xsre2/10.E0
*         if (nexpo .gt. 40) write(6,*) (b/a)
         if (xsre2/xsre3 .gt. 9.999D0) then
          xsre1=xsre2
         else
         if (xsre2.gt.xsre1) call gibtem(xii)
           goto 130
         endif
      enddo
       goto 130
*      write(6,*) ' pas marche ', nexpo,xsre1
**    moterr(1:6) = 'XSPETI'
**    call erreur(907)
**    return
 
 130  continue
*      write(6,*) 'nexpo' , nexpo, xsre1
      xspeti = xsre1 * 1E2
*      write(6,*) 'xspeti', xspeti
 
c calcul xsgran
      xsre1 = 1.E0
      do nexpo = 1,37
          xsre2 = xsre1*10.E0
          xsre3 = xsre2*10.E0
*          if (nexpo .gt. 30) write(6,*) a, (b.gt.a)
        if (xsre3 .gt. xsre2) then
        xsre1 = xsre2
          else
         if (xsre2.lt.xsre1) call gibtem(xii)
        goto 140
       endif
      enddo
      goto 140
**    moterr(1:6) = 'XSGRAN'
**    call erreur(907)
*      write(6,*) ' pas marche ', nexpo,xsre1
**    return
 
 140  continue
*       write(6,*) 'nexpo' , nexpo, xsre1
       xsgran = xsre1/1.E3
*       write(6,*) 'xsgran', xsgran
       amplit=xsgran
c  pour pouvoir inverser
      xsauxi = min(abs(log(xsgran)),abs(log(xspeti)))
*       write(6,*) xsauxi, abs(log(xsgran)) , abs(log(xspeti))
      if (abs(xsauxi - abs(log(xsgran))).lt.0.5) then
        xspeti = 1.E0/xsgran
*         write(6,*) 'xspeti', xspeti, 1/xspeti, 1/xsgran
      else
        xsgran = 1.E0/xspeti
*        write(6,*) 'xsgran ' , xsgran, 1/xspeti, 1/xsgran
      endif
 
c calcul xzprec ,
      xre1 = 1.D0
      do nexpol = 1,10000
        xre2=1.d0+xre1/10.d0
        if (xre2.le.1.d0) then
        if (xre2.lt.xre1) call gibtem(xii)
        goto 150
        endif
        xre1 = xre2-1.d0
      enddo
      moterr(1:6) = 'XZPREC'
      call erreur(907)
*      write(6,*) ' pas marche ', nexpo,xre1
      return
 
 150  continue
*      write(6,*) 'nexpo' , nexpo, xre1
      xzprec = xre1
*      write(6,*) 'xzprec', xzprec
 
c calcul xszpre ,
      xsre1 = 1.E0
      do nexpo = 1,10000
        xsre2=1.E0+xsre1/10.E0
        if (xsre2.le.1.e0) then
        if (xsre2.lt.xsre1) call gibtem(xii)
        goto 160
        endif
        xsre1 = xsre2-1E0
      enddo
      moterr(1:6) = 'XSZPRE'
      call erreur(907)
*      write(6,*) ' pas marche ', nexpo,xsre1
      return
 
 160  continue
*      write(6,*) 'nexpo' , nexpo, xsre1
      xszpre = xsre1
*      write(6,*) 'xszpre', xszpre
      return
 
      END
 
 
 
 

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