author hristov Tue, 21 Dec 2010 15:16:09 +0000 (15:16 +0000) committer hristov Tue, 21 Dec 2010 15:16:09 +0000 (15:16 +0000)

index b0c1bb9610babebbfaf4a048128d77b2cb1ae83d..521a846659eed98cec1ee493c1b3dd3d6bf5a4cf 100644 (file)
C...(Borrow slot N+2 for temporary direction.)
DO 850 I=IC1+1,IC2-1
IF((K(I,1).EQ.1.OR.K(I,1).EQ.2).AND.
&  KCHG(PYCOMP(K(I,2)),2).NE.0) THEN
-          FRAC1=FOUR(IC2,I)/(FOUR(IC1,I)+FOUR(IC2,I))
+          IF (ABS(FOUR(IC1,I)+FOUR(IC2,I)).GT.0.D0) THEN
+             FRAC1=FOUR(IC2,I)/(FOUR(IC1,I)+FOUR(IC2,I))
+          ELSE
+             FRAC1 = 1.D0
+          ENDIF
DO 840 J=1,4
P(N+2,J)=P(N+2,J)+FRAC1*P(I,J)
840     CONTINUE
C...Calculate Pg, a part of which will be added to Phad later. (EN)
ALPHA=1D0
BETA=1D0
ELSE
-          ALPHA=FOUR(NSAV+1,I2)/FOUR(I1,I2)
-          BETA=FOUR(NSAV+1,I1)/FOUR(I1,I2)
+           IF (ABS(FOUR(I1,I2)).GT.0.D0) THEN
+              ALPHA=FOUR(NSAV+1,I2)/FOUR(I1,I2)
+              BETA=FOUR(NSAV+1,I1)/FOUR(I1,I2)
+           ELSE
+              ALPHA=1D0
+              BETA=1D0
+           ENDIF
ENDIF
DO 1030 J=1,4
PG(J)=ALPHA*P(I1,J)+BETA*P(I2,J)