adi.f

      SUBROUTINE ADI(A,B,C,D,E,F,G,U,JMAX,K,ALPHA,BETA,EPS)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(JJ=100,KK=6,NRR=32,MAXITS=100,ZERO=0.D0,TWO=2.D0,HALF=.5
     *D0)
      DIMENSION A(JMAX,JMAX),B(JMAX,JMAX),C(JMAX,JMAX),D(JMAX,JMAX),
     *    E(JMAX,JMAX),F(JMAX,JMAX),G(JMAX,JMAX),U(JMAX,JMAX),
     *    AA(JJ),BB(JJ),CC(JJ),RR(JJ),UU(JJ),PSI(JJ,JJ),
     *    ALPH(KK),BET(KK),R(NRR),S(NRR,KK)
      IF(JMAX.GT.JJ)PAUSE 'Increase JJ'
      IF(K.GT.KK-1)PAUSE 'Increase KK'
      K1=K+1
      NR=2**K
      ALPH(1)=ALPHA
      BET(1)=BETA
      DO 11 J=1,K
        ALPH(J+1)=SQRT(ALPH(J)*BET(J))
        BET(J+1)=HALF*(ALPH(J)+BET(J))
11    CONTINUE
      S(1,1)=SQRT(ALPH(K1)*BET(K1))
      DO 13 J=1,K
        AB=ALPH(K1-J)*BET(K1-J)
        DO 12 N=1,2**(J-1)
          DISC=SQRT(S(N,J)**2-AB)
          S(2*N,J+1)=S(N,J)+DISC
          S(2*N-1,J+1)=AB/S(2*N,J+1)
12      CONTINUE
13    CONTINUE
      DO 14 N=1,NR
        R(N)=S(N,K1)
14    CONTINUE
      ANORMG=ZERO
      DO 16 J=2,JMAX-1
        DO 15 L=2,JMAX-1
          ANORMG=ANORMG+ABS(G(J,L))
          PSI(J,L)=-D(J,L)*U(J,L-1)+(R(1)-E(J,L))*U(J,L)
     *        -F(J,L)*U(J,L+1)
15      CONTINUE
16    CONTINUE
      NITS=MAXITS/NR
      DO 27 KITS=1,NITS
        DO 24 N=1,NR
          IF(N.EQ.NR)THEN
            NEXT=1
          ELSE
            NEXT=N+1
          ENDIF
          RFACT=R(N)+R(NEXT)
          DO 19 L=2,JMAX-1
            DO 17 J=2,JMAX-1
              AA(J-1)=A(J,L)
              BB(J-1)=B(J,L)+R(N)
              CC(J-1)=C(J,L)
              RR(J-1)=PSI(J,L)-G(J,L)
17          CONTINUE
            CALL TRIDAG(AA,BB,CC,RR,UU,JMAX-2)
            DO 18 J=2,JMAX-1
              PSI(J,L)=-PSI(J,L)+TWO*R(N)*UU(J-1)
18          CONTINUE
19        CONTINUE
          DO 23 J=2,JMAX-1
            DO 21 L=2,JMAX-1
              AA(L-1)=D(J,L)
              BB(L-1)=E(J,L)+R(N)
              CC(L-1)=F(J,L)
              RR(L-1)=PSI(J,L)
21          CONTINUE
            CALL TRIDAG(AA,BB,CC,RR,UU,JMAX-2)
            DO 22 L=2,JMAX-1
              U(J,L)=UU(L-1)
              PSI(J,L)=-PSI(J,L)+RFACT*UU(L-1)
22          CONTINUE
23        CONTINUE
24      CONTINUE
        ANORM=ZERO
        DO 26 J=2,JMAX-1
          DO 25 L=2,JMAX-1
            RESID=A(J,L)*U(J-1,L)+(B(J,L)+E(J,L))*U(J,L)
     *          +C(J,L)*U(J+1,L)+D(J,L)*U(J,L-1)
     *          +F(J,L)*U(J,L+1)+G(J,L)
            ANORM=ANORM+ABS(RESID)
25        CONTINUE
26      CONTINUE
        IF(ANORM.LT.EPS*ANORMG)RETURN
27    CONTINUE
      PAUSE 'MAXITS exceeded'
      END