Esope source files

   1 : C MOM       SOURCE    CHAT      05/01/13    01:51:10     5004           
   2 :       SUBROUTINE MOM(G,D,A,FAULT)
   3 :       IMPLICIT INTEGER(I-N)
   4 :       IMPLICIT REAL*8 (A-H,O-Z)
   5 :       REAL*8 A(6), B(6), C(6), G, D, ZZ, VV, Rdeux, 
   6 :      $ RIpi, W, E, R, H, T,U, Y, X, V, F, Z, S, P, Q, AA, AB, 
   7 :      $ EXPA, EXPB
   8 :       LOGICAL L, FAULT
   9 :       DATA ZZ, VV, LIMIT / 1.0E-5, 1.0E-8, 5000 /
  10 : C
  11 : C   Rdeux = sqrt(2)
  12 : C   RIpi = 1 / sqrt(pi)
  13 : C   EXPA = valeur t.q. EXP(EXPA) ne cause pas un overflow
  14 : C   EXPB = valeur t.q. 1 - EXP(EXPB) peut etre prise pour 1
  15 : C
  16 :       DATA Rdeux, RIpi, EXPA, EXPB
  17 :      $ / 1.414213562, 0.5641895835, 80.0, 23.7 /
  18 : C
  19 : C
  20 :       FAULT = .FALSE.
  21 :       DO 10 I=1,6
  22 :   10    C(I) = 0.D0
  23 :       W = G / D
  24 : C
  25 : C   Essai de valeur pour H
  26 : C      
  27 :       IF (W.GT.EXPA) GOTO 140
  28 :       E = DEXP(W) + 1.D0
  29 :       R = Rdeux / D
  30 :       H = 0.75D0
  31 :       IF (D.LT.3.D0) H = 0.25D0 * D
  32 :       K = 1
  33 :       GOTO 40
  34 : C
  35 : C  Debut de la boucle exterieur
  36 : C
  37 :   20  K = K + 1
  38 :       IF (K.GT.LIMIT) GOTO 140
  39 :       DO 30 I=1,6
  40 :   30  C(I) = A(I)
  41 : C
  42 : C  Pas de convergence - essai d'un H plus petit
  43 : C
  44 :       H = 0.5D0 * H
  45 :   40  T = W
  46 :       U = T
  47 :       Y = H * H
  48 :       X = 2.D0 * Y
  49 :       A(1) = 1.D0 / E
  50 :       DO 50 I=2,6
  51 :   50  A(I) = A(I-1) / E
  52 :       V = Y
  53 :       F = R * H
  54 :       M = 0
  55 : C
  56 : C Debut de la boucle interieur
  57 : C
  58 :   60  M = M + 1
  59 :       IF (M.GT.LIMIT) GOTO 140
  60 :       DO 70 I=1,6
  61 :   70  B(I) = A(I)
  62 :       U = U - F
  63 :       Z = 1.D0
  64 :       IF (U.GT.-EXPB) Z = DEXP(U) + Z
  65 :       T = T + F
  66 :       L = T.GT.EXPB
  67 :       IF (.NOT.L) S = DEXP(T) + 1.D0
  68 :       P = DEXP(-V) 
  69 :       Q = P
  70 :       DO 90 I=1,6
  71 :         AA = A(I)
  72 :         P = P / Z
  73 :         AB = AA
  74 :         AA = AA + P
  75 :         IF (AA.EQ.AB) GOTO 100
  76 :         IF (L) GOTO 80
  77 :         Q = Q / S
  78 :         AB = AA
  79 :         AA = AA + Q
  80 :         L = AA.EQ.AB
  81 :   80    A(I) = AA
  82 :   90  CONTINUE
  83 :  100  Y = Y + X
  84 :       V = V + Y
  85 :       DO 110 I=1,6
  86 :         IF (A(I).EQ.0.D0) GOTO 140
  87 :         IF (DABS((A(I) - B(I)) / A(I)).GT.VV) GOTO 60
  88 :  110  CONTINUE
  89 : C
  90 : C Fin de la boucle interieur
  91 : C
  92 :       V = RIpi * H
  93 :       DO 120 I=1,6
  94 :  120  A(I) = V * A(I)
  95 :       DO 130 I=1,6
  96 :         IF (A(I).EQ.0.D0) GOTO 140
  97 :         IF (DABS((A(I) - C(I)) / A(I)).GT.ZZ) GOTO 20
  98 :  130  CONTINUE
  99 : C
 100 : C Fin de la boucle exterieur
 101 : C
 102 :       RETURN
 103 :  140  FAULT = .TRUE.
 104 :       RETURN
 105 :       END
 106 :  
 107 :  
 108 :