      SUBROUTINE jacobi(a,n,np,d,v,nrot)                                JACO001
c --- modified from numerical recipes ----------------------------------JACO002
c                                                                       JACO003
c     one needs to set the threshold for sm from sm.eq.0 to sm.lt.10E-10JACO004
c     (anything in this range would be ok) due to underflow errors on   JACO005
c     some computers/compilers.                                         JACO006
c ----------------------------------------------------------------------JACO007
      INTEGER n,np,nrot,NMAX                                            JACO008
      REAL a(np,np),d(np),v(np,np)                                      JACO009
      PARAMETER (NMAX=500)                                              JACO010
      INTEGER i,ip,iq,j                                                 JACO011
      REAL c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX)                 JACO012
      do 12 ip=1,n                                                      JACO013
        do 11 iq=1,n                                                    JACO014
          v(ip,iq)=0.                                                   JACO015
11      continue                                                        JACO016
        v(ip,ip)=1.                                                     JACO017
12    continue                                                          JACO018
      do 13 ip=1,n                                                      JACO019
        b(ip)=a(ip,ip)                                                  JACO020
        d(ip)=b(ip)                                                     JACO021
        z(ip)=0.                                                        JACO022
13    continue                                                          JACO023
      nrot=0                                                            JACO024
      do 24 i=1,50                                                      JACO025
        sm=0.                                                           JACO026
        do 15 ip=1,n-1                                                  JACO027
          do 14 iq=ip+1,n                                               JACO028
            sm=sm+abs(a(ip,iq))                                         JACO029
14        continue                                                      JACO030
15      continue                                                        JACO031
c ---   modified convergence threshold ---                              JACO032
        if(sm.lt.10E-10)return                                          JACO033
        if(i.lt.4)then                                                  JACO034
          tresh=0.2*sm/n**2                                             JACO035
        else                                                            JACO036
          tresh=0.                                                      JACO037
        endif                                                           JACO038
        do 22 ip=1,n-1                                                  JACO039
          do 21 iq=ip+1,n                                               JACO040
            g=100.*abs(a(ip,iq))                                        JACO041
            if((i.gt.4).and.(abs(d(ip))+                                JACO042
     &g.eq.abs(d(ip))).and.(abs(d(iq))+g.eq.abs(d(iq))))then            JACO043
              a(ip,iq)=0.                                               JACO044
            else if(abs(a(ip,iq)).gt.tresh)then                         JACO045
              h=d(iq)-d(ip)                                             JACO046
              if(abs(h)+g.eq.abs(h))then                                JACO047
                t=a(ip,iq)/h                                            JACO048
              else                                                      JACO049
                theta=0.5*h/a(ip,iq)                                    JACO050
                t=1./(abs(theta)+sqrt(1.+theta**2))                     JACO051
                if(theta.lt.0.)t=-t                                     JACO052
              endif                                                     JACO053
              c=1./sqrt(1+t**2)                                         JACO054
              s=t*c                                                     JACO055
              tau=s/(1.+c)                                              JACO056
              h=t*a(ip,iq)                                              JACO057
              z(ip)=z(ip)-h                                             JACO058
              z(iq)=z(iq)+h                                             JACO059
              d(ip)=d(ip)-h                                             JACO060
              d(iq)=d(iq)+h                                             JACO061
              a(ip,iq)=0.                                               JACO062
              do 16 j=1,ip-1                                            JACO063
                g=a(j,ip)                                               JACO064
                h=a(j,iq)                                               JACO065
                a(j,ip)=g-s*(h+g*tau)                                   JACO066
                a(j,iq)=h+s*(g-h*tau)                                   JACO067
16            continue                                                  JACO068
              do 17 j=ip+1,iq-1                                         JACO069
                g=a(ip,j)                                               JACO070
                h=a(j,iq)                                               JACO071
                a(ip,j)=g-s*(h+g*tau)                                   JACO072
                a(j,iq)=h+s*(g-h*tau)                                   JACO073
17            continue                                                  JACO074
              do 18 j=iq+1,n                                            JACO075
                g=a(ip,j)                                               JACO076
                h=a(iq,j)                                               JACO077
                a(ip,j)=g-s*(h+g*tau)                                   JACO078
                a(iq,j)=h+s*(g-h*tau)                                   JACO079
18            continue                                                  JACO080
              do 19 j=1,n                                               JACO081
                g=v(j,ip)                                               JACO082
                h=v(j,iq)                                               JACO083
                v(j,ip)=g-s*(h+g*tau)                                   JACO084
                v(j,iq)=h+s*(g-h*tau)                                   JACO085
19            continue                                                  JACO086
              nrot=nrot+1                                               JACO087
            endif                                                       JACO088
21        continue                                                      JACO089
22      continue                                                        JACO090
        do 23 ip=1,n                                                    JACO091
          b(ip)=b(ip)+z(ip)                                             JACO092
          d(ip)=b(ip)                                                   JACO093
          z(ip)=0.                                                      JACO094
23      continue                                                        JACO095
24    continue                                                          JACO096
      pause 'too many iterations in jacobi'                             JACO097
      return                                                            JACO098
      END                                                               JACO099
C  (C) Copr. 1986-92 Numerical Recipes Software A2.Q2$2500.             JACO100
