--- trunk/filtrez/jacobi.f	2013/11/15 18:45:49	76
+++ trunk/filtrez/jacobi.f90	2014/03/05 14:38:41	81
@@ -1,97 +1,99 @@
-!
-! $Header: /home/cvsroot/LMDZ4/libf/filtrez/jacobi.F,v 1.1.1.1 2004/05/19 12:53:09 lmdzadmin Exp $
-!
-      SUBROUTINE JACOBI(A,N,NP,D,V,NROT)
-      PARAMETER (NMAX=400)
-      DIMENSION A(NP,NP),D(NP),V(NP,NP),B(NMAX),Z(NMAX)
-      IF (n.gt.nmax) THEN
-         print*, 'n, nmax=', n, nmax
-         print*, 'Surdimensionnement insuffisant dans jacobi'
-         stop 1
-      ENDIF
-      DO 12 IP=1,N
-        DO 11 IQ=1,N
-          V(IP,IQ)=0.
-11      CONTINUE
-        V(IP,IP)=1.
-12    CONTINUE
-      DO 13 IP=1,N
-        B(IP)=A(IP,IP)
-        D(IP)=B(IP)
-        Z(IP)=0.
-13    CONTINUE
-      NROT=0
-      DO 24 I=1,50
-        SM=0.
-        DO 15 IP=1,N-1
-          DO 14 IQ=IP+1,N
-            SM=SM+ABS(A(IP,IQ))
-14        CONTINUE
-15      CONTINUE
-        IF(SM.EQ.0.)RETURN
-        IF(I.LT.4)THEN
-          TRESH=0.2*SM/N**2
-        ELSE
-          TRESH=0.
-        ENDIF
-        DO 22 IP=1,N-1
-          DO 21 IQ=IP+1,N
-            G=100.*ABS(A(IP,IQ))
-            IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP)))
-     *         .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN
-              A(IP,IQ)=0.
-            ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN
-              H=D(IQ)-D(IP)
-              IF(ABS(H)+G.EQ.ABS(H))THEN
-                T=A(IP,IQ)/H
-              ELSE
-                THETA=0.5*H/A(IP,IQ)
-                T=1./(ABS(THETA)+SQRT(1.+THETA**2))
-                IF(THETA.LT.0.)T=-T
-              ENDIF
-              C=1./SQRT(1+T**2)
-              S=T*C
-              TAU=S/(1.+C)
-              H=T*A(IP,IQ)
-              Z(IP)=Z(IP)-H
-              Z(IQ)=Z(IQ)+H
-              D(IP)=D(IP)-H
-              D(IQ)=D(IQ)+H
-              A(IP,IQ)=0.
-              DO 16 J=1,IP-1
-                G=A(J,IP)
-                H=A(J,IQ)
-                A(J,IP)=G-S*(H+G*TAU)
-                A(J,IQ)=H+S*(G-H*TAU)
-16            CONTINUE
-              DO 17 J=IP+1,IQ-1
-                G=A(IP,J)
-                H=A(J,IQ)
-                A(IP,J)=G-S*(H+G*TAU)
-                A(J,IQ)=H+S*(G-H*TAU)
-17            CONTINUE
-              DO 18 J=IQ+1,N
-                G=A(IP,J)
-                H=A(IQ,J)
-                A(IP,J)=G-S*(H+G*TAU)
-                A(IQ,J)=H+S*(G-H*TAU)
-18            CONTINUE
-              DO 19 J=1,N
-                G=V(J,IP)
-                H=V(J,IQ)
-                V(J,IP)=G-S*(H+G*TAU)
-                V(J,IQ)=H+S*(G-H*TAU)
-19            CONTINUE
-              NROT=NROT+1
-            ENDIF
-21        CONTINUE
-22      CONTINUE
-        DO 23 IP=1,N
-          B(IP)=B(IP)+Z(IP)
-          D(IP)=B(IP)
-          Z(IP)=0.
-23      CONTINUE
-24    CONTINUE
-      STOP '50 iterations should never happen'
-      RETURN
-      END
+
+! $Header: /home/cvsroot/LMDZ4/libf/filtrez/jacobi.F,v 1.1.1.1 2004/05/19
+! 12:53:09 lmdzadmin Exp $
+
+SUBROUTINE jacobi(a, n, np, d, v, nrot)
+  PARAMETER (nmax=400)
+  DIMENSION a(np, np), d(np), v(np, np), b(nmax), z(nmax)
+
+  IF (n>nmax) THEN
+    PRINT *, 'n, nmax=', n, nmax
+    PRINT *, 'Surdimensionnement insuffisant dans jacobi'
+    STOP 1
+  END IF
+  DO ip = 1, n
+    DO iq = 1, n
+      v(ip, iq) = 0.
+    END DO
+    v(ip, ip) = 1.
+  END DO
+  DO ip = 1, n
+    b(ip) = a(ip, ip)
+    d(ip) = b(ip)
+    z(ip) = 0.
+  END DO
+  nrot = 0
+  DO i = 1, 50
+    sm = 0.
+    DO ip = 1, n - 1
+      DO iq = ip + 1, n
+        sm = sm + abs(a(ip,iq))
+      END DO
+    END DO
+    IF (sm==0.) RETURN
+    IF (i<4) THEN
+      tresh = 0.2*sm/n**2
+    ELSE
+      tresh = 0.
+    END IF
+    DO ip = 1, n - 1
+      DO iq = ip + 1, n
+        g = 100.*abs(a(ip,iq))
+        IF ((i>4) .AND. (abs(d(ip))+g==abs(d(ip))) .AND. (abs(d(iq))+g==abs( &
+            d(iq)))) THEN
+          a(ip, iq) = 0.
+        ELSE IF (abs(a(ip,iq))>tresh) THEN
+          h = d(iq) - d(ip)
+          IF (abs(h)+g==abs(h)) THEN
+            t = a(ip, iq)/h
+          ELSE
+            theta = 0.5*h/a(ip, iq)
+            t = 1./(abs(theta)+sqrt(1.+theta**2))
+            IF (theta<0.) t = -t
+          END IF
+          c = 1./sqrt(1+t**2)
+          s = t*c
+          tau = s/(1.+c)
+          h = t*a(ip, iq)
+          z(ip) = z(ip) - h
+          z(iq) = z(iq) + h
+          d(ip) = d(ip) - h
+          d(iq) = d(iq) + h
+          a(ip, iq) = 0.
+          DO j = 1, ip - 1
+            g = a(j, ip)
+            h = a(j, iq)
+            a(j, ip) = g - s*(h+g*tau)
+            a(j, iq) = h + s*(g-h*tau)
+          END DO
+          DO j = ip + 1, iq - 1
+            g = a(ip, j)
+            h = a(j, iq)
+            a(ip, j) = g - s*(h+g*tau)
+            a(j, iq) = h + s*(g-h*tau)
+          END DO
+          DO j = iq + 1, n
+            g = a(ip, j)
+            h = a(iq, j)
+            a(ip, j) = g - s*(h+g*tau)
+            a(iq, j) = h + s*(g-h*tau)
+          END DO
+          DO j = 1, n
+            g = v(j, ip)
+            h = v(j, iq)
+            v(j, ip) = g - s*(h+g*tau)
+            v(j, iq) = h + s*(g-h*tau)
+          END DO
+          nrot = nrot + 1
+        END IF
+      END DO
+    END DO
+    DO ip = 1, n
+      b(ip) = b(ip) + z(ip)
+      d(ip) = b(ip)
+      z(ip) = 0.
+    END DO
+  END DO
+  STOP '50 iterations should never happen'
+  RETURN
+END SUBROUTINE jacobi