Sources Esope | Cast3M

mom
C MOM       SOURCE    CB215821  16/04/21    21:17:49     8920
      SUBROUTINE MOM(G,D,A,FAULT)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 A(6), B(6), C(6), G, D, ZZ, VV, Rdeux,
     $ RIpi, W, E, R, H, T,U, Y, X, V, F, Z, S, P, Q, AA, AB,
     $ EXPA, EXPB
      LOGICAL L, FAULT
      DATA ZZ, VV, LIMIT / 1.0E-5, 1.0E-8, 5000 /
C
C   Rdeux = sqrt(2)
C   RIpi = 1 / sqrt(pi)
C   EXPA = valeur t.q. EXP(EXPA) ne cause pas un overflow
C   EXPB = valeur t.q. 1 - EXP(EXPB) peut etre prise pour 1
C
      DATA Rdeux, RIpi, EXPA, EXPB
     $ / 1.414213562, 0.5641895835, 80.0, 23.7 /
C
C
      FAULT = .FALSE.
      DO 10 I=1,6
  10    C(I) = 0.D0
      W = G / D
C
C   Essai de valeur pour H
C
      IF (W.GT.EXPA) GOTO 140
      E = EXP(W) + 1.D0
      R = Rdeux / D
      H = 0.75D0
      IF (D.LT.3.D0) H = 0.25D0 * D
      K = 1
      GOTO 40
C
C  Debut de la boucle exterieur
C
  20  K = K + 1
      IF (K.GT.LIMIT) GOTO 140
      DO 30 I=1,6
  30  C(I) = A(I)
C
C  Pas de convergence - essai d'un H plus petit
C
      H = 0.5D0 * H
  40  T = W
      U = T
      Y = H * H
      X = 2.D0 * Y
      A(1) = 1.D0 / E
      DO 50 I=2,6
  50  A(I) = A(I-1) / E
      V = Y
      F = R * H
      M = 0
C
C Debut de la boucle interieur
C
  60  M = M + 1
      IF (M.GT.LIMIT) GOTO 140
      DO 70 I=1,6
  70  B(I) = A(I)
      U = U - F
      Z = 1.D0
      IF (U.GT.-EXPB) Z = EXP(U) + Z
      T = T + F
      L = T.GT.EXPB
      IF (.NOT.L) S = EXP(T) + 1.D0
      P = EXP(-V)
      Q = P
      DO 90 I=1,6
        AA = A(I)
        P = P / Z
        AB = AA
        AA = AA + P
        IF (AA.EQ.AB) GOTO 100
        IF (L) GOTO 80
        Q = Q / S
        AB = AA
        AA = AA + Q
        L = AA.EQ.AB
  80    A(I) = AA
  90  CONTINUE
 100  Y = Y + X
      V = V + Y
      DO 110 I=1,6
        IF (A(I).EQ.0.D0) GOTO 140
        IF (ABS((A(I) - B(I)) / A(I)).GT.VV) GOTO 60
 110  CONTINUE
C
C Fin de la boucle interieur
C
      V = RIpi * H
      DO 120 I=1,6
 120  A(I) = V * A(I)
      DO 130 I=1,6
        IF (A(I).EQ.0.D0) GOTO 140
        IF (ABS((A(I) - C(I)) / A(I)).GT.ZZ) GOTO 20
 130  CONTINUE
C
C Fin de la boucle exterieur
C
      RETURN
 140  FAULT = .TRUE.
      RETURN
      END