subroutine fourier_matrix
      parameter(igit=64,nb=36,nfmax=igit/2+1,nfmaxl=10,nfleg=10)
      real fun(nfmax,igit)
      real theta(nb),pleg(nfmaxl,nb)
      common/legandr/pfleg(nb,nb)
      common/fourier/pf1(igit,igit),pf3(igit,igit),pf5 (igit,igit),
     $               pf7(igit,igit),pf9(igit,igit),pf11(igit,igit),
     $                              pf2(igit,igit),pf12(igit,igit)
      real a(nfmax,nfmax),am1(nfmax,nfmax),u(nfmax,nfmax),
     $     v(nfmax,nfmax),sigma(nfmax),wk1(nfmax),wk2(nfmax)
      real al(nfmaxl,nfmaxl),am1l(nfmaxl,nfmaxl),ul(nfmaxl,nfmaxl),
     $     vl(nfmaxl,nfmaxl),sigmal(nfmaxl),wk1l(nfmaxl),wk2l(nfmaxl)
      pi=2.*asin(1.)
      dlon=2.*pi/float(igit)
      do i=1,igit
        xlon=float(i-1)*dlon
        fun(1,i)=1.
        do if=1,igit/4
        iif=2*if
        fun(iif  ,i)=sin(float(if)*xlon)
        fun(iif+1,i)=cos(float(if)*xlon)
        end do
      end do
      DO IF=1,IGIT/2+1
       SNORM=0.
       DO I=1,IGIT
          SNORM=SNORM+FUN(IF,I)*FUN(IF,I)
       END DO
        DO I=1,IGIT
          FUN(IF,I)=FUN(IF,I)/SQRT(SNORM)
        END DO
      END DO
      relerr=1.e-4
      call filter(igit,igit,igit,nfmax, 3,fun,fun,a,am1,
     $            u,v,sigma,wk1,wk2, pf1,relerr)
      call filter(igit,igit,igit,nfmax, 5,fun,fun,a,am1,
     $            u,v,sigma,wk1,wk2, pf2,relerr)
      call filter(igit,igit,igit,nfmax, 7,fun,fun,a,am1,
     $            u,v,sigma,wk1,wk2, pf3,relerr)
      call filter(igit,igit,igit,nfmax,11,fun,fun,a,am1,
     $            u,v,sigma,wk1,wk2, pf5,relerr)
      call filter(igit,igit,igit,nfmax,15,fun,fun,a,am1,
     $            u,v,sigma,wk1,wk2, pf7,relerr)
      call filter(igit,igit,igit,nfmax,19,fun,fun,a,am1,
     $            u,v,sigma,wk1,wk2, pf9,relerr)
      call filter(igit,igit,igit,nfmax,23,fun,fun,a,am1,
     $            u,v,sigma,wk1,wk2,pf11,relerr)
      call filter(igit,igit,igit,nfmax,25,fun,fun,a,am1,
     $            u,v,sigma,wk1,wk2,pf12,relerr)
      do j=1,nb 
c------------------------------------------- co-latitude COMMA
        theta(j)=pi/180.*(2.5+ 5.*FLOAT(j-1))
      end do
      CALL POLLEG(nb,nfmaxl,nfleg,1,theta,pleg)
c      CALL PLLEG0(NJ,NFMAX,NF,ARG,PL0) 
      CALL FILTER(nb,nb,nb,nfmaxl,nfleg,pleg,pleg,
     $            Al,AM1l,Ul,Vl,SIGMAl,WK1l,WK2l,pfleg,RELERR)
      return
      end 
c	------------------------------------------------
      SUBROUTINE FILTER(NJMAX,NJS,NJD,NFMAX,NF,
     $ PLS,PLD,A,AM1,U,V,SIGMA,WK1,WK2,PF,RELERR)
      REAL PLD(NFMAX,NJD),PLS(NFMAX,NJS),
     $     A(NFMAX,NF),AM1(NFMAX,NF),U(NFMAX,NF),
     $     V(NFMAX,NF),SIGMA(NF),WK1(NF),WK2(NF),
     $     PF(NJMAX,NJD),WK(2*NF)
      DO 10 IF=1,NF
      DO 20 KF=1,NF
      S=0.
      DO 30 JD=1,NJD
   30 S=S+PLD(IF,JD)*PLD(KF,JD)
      A(KF,IF)=S
   20 CONTINUE
   10 CONTINUE
c      CALL INVSVD(NFMAX,NF,NF,A,AM1,U,V,SIGMA,WK1,RELERR)
      call lginf(a,nfmax,nf,nf,relerr,am1,nfmax,sigma,wk,ier)
      DO 40 JD=1,NJD
      DO 50 IF=1,NF
      S=0.
      DO 55 KF=1,NF
   55 S=S+AM1(KF,IF)*PLD(KF,JD)
      WK1(IF)=S
   50 CONTINUE
      DO 60 JS=1,NJS
      DO 70 KF=1,NF
c---- filter out the wave with zonal wave number 6, 1, 2
c      if(kf.eq.12.or.kf.eq.13) then
c      if(kf.eq.2.or.kf.eq.3) then
c      if(kf.eq.4.or.kf.eq.5) then
c        wk2(kf)=0.
c                             else
        WK2(KF)=PLS(KF,JS)
c      end if
   70 continue
      CALL PRVV(WK1,NF,WK2,RESULT)
      PF(JS,JD)=RESULT
   60 CONTINUE
   40 CONTINUE
      RETURN
      END
c------------------------------------------------------------
      SUBROUTINE INVSVD(NDIM,M,N,A,AM1,U,V,SIGMA,WORK,RELERR)
      DOUBLE PRECISION S
      REAL A(NDIM,N),U(NDIM,N),V(NDIM,N),
     *     SIGMA(N),WORK(N),AM1(NDIM,M)
      CALL SVD(NDIM,M,N,A,SIGMA,.TRUE.,U,.TRUE.,V,IERR,WORK)
      IF(IERR.NE.0) WRITE(99,1) IERR
   1  FORMAT(/'   /INVSVD/    TROUBLE.IERR = ',I4/)
      SIGMA1=0.
      DO 2 I=1,N
      IF(SIGMA(I).GT.SIGMA1) SIGMA1=SIGMA(I)
   2  CONTINUE
      TAU=RELERR*SIGMA1
      DO 4 I=1,N
      DO 4 K=1,M
      S=0.0D0
      DO 3 J=1,N
      IF(SIGMA(J).LE.TAU) GO TO 3
      S=S+DBLE(V(I,J)*U(K,J)/SIGMA(J))
   3  CONTINUE
      AM1(I,K)=S
   4  CONTINUE
      RETURN
      END
c------------------------------------------------------
      SUBROUTINE SVD(NM,M,N,A,W,MATU,U,MATV,V,IERR,RV1)
C
      REAL A(NM,N),W(N),U(NM,N),V(NM,N),RV1(N)
C
      LOGICAL MATU,MATV
C
                       IERR=0
      DO 100 I=1,M
      DO 100 J=1,N
  100 U(I,J)=A(I,J)
      G=0.
      SCALE=0.
      ANORM=0.
      DO 300 I=1,N
      L=I+1
      RV1(I)=SCALE*G
      G=0.
      S=0.
      SCALE=0.
      IF(I.GT.M) GO TO 210
      DO 120 K=I,M
      SCALE=SCALE+ABS(U(K,I))
  120 CONTINUE
      IF(SCALE.EQ.0.) GO TO 210
      DO 130 K=I,M
      U(K,I)=U(K,I)/SCALE
      S=S+U(K,I)**2
  130 CONTINUE
      F=U(I,I)
      G=-SIGN(SQRT(S),F)
      H=F*G-S
      U(I,I)=F-G
      IF(I.EQ.N) GO TO 190
      DO 150 J=L,N
      S=0.
      DO 140 K=I,M
      S=S+U(K,I)*U(K,J)
  140 CONTINUE
      F=S/H
      DO 150 K=I,M
      U(K,J)=U(K,J)+F*U(K,I)
  150 CONTINUE
  190 DO 200 K=I,M
  200 U(K,I)=SCALE*U(K,I)
  210 W(I)=SCALE*G
      G=0.
      S=0.
      IF(I.GT.M.OR.I.EQ.N) GO TO 290
      DO 220 K=L,N
  220 SCALE=SCALE+ABS(U(I,K))
C
      IF(SCALE.EQ.0.) GO TO 290
      DO 230 K=L,N
      U(I,K)=U(I,K)/SCALE
      S=S+U(I,K)**2
  230 CONTINUE
      F=U(I,L)
      G=-SIGN(SQRT(S),F)
      H=F*G-S
      U(I,L)=F-G
      DO 240 K=L,N
  240 RV1(K)=U(I,K)/H
      IF(I.EQ.M) GO TO 270
      DO 260 J=L,M
      S=0.
      DO 250 K=L,N
  250 S=S+U(J,K)*U(I,K)
      DO 260 K=L,N
  260 U(J,K)=U(J,K)+S*RV1(K)
  270 DO 280 K=L,N
  280 U(I,K)=SCALE*U(I,K)
  290 ANORM=AMAX1(ANORM,ABS(W(I))+ABS(RV1(I)))
  300 CONTINUE
      IF(.NOT.MATV) GO TO 410
      DO 400 II=1,N
      I=N+1-II
      IF(I.EQ.N ) GO TO 390
      IF(G.EQ.0.) GO TO 360
      DO 320 J=L,N
  320 V(J,I)=(U(I,J)/U(I,L))/G
      DO 350 J=L,N
      S=0.
      DO 340 K=L,N
  340 S=S+U(I,K)*V(K,J)
      DO 350 K=L,N
  350 V(K,J)=V(K,J)+S*V(K,I)
  360 DO 380 J=L,N
      V(I,J)=0.
      V(J,I)=0.
  380 CONTINUE
  390 V(I,I)=1.
      G=RV1(I)
      L=I
  400 CONTINUE
  410 IF(.NOT.MATU) GO TO 510
      MN=N
      IF(M.LT.N)MN=M
      DO 500 II=1,MN
      I=MN+1-II
      L=I+1
      G=W(I)
      IF(I.EQ.N)  GO TO 430
      DO 420 J=L,N
  420 U(I,J)=0.
  430 IF(G.EQ.0.) GO TO 475
      IF(I.EQ.MN) GO TO 460
      DO 450 J=L,N
      S=0.
      DO 440 K=L,M
  440 S=S+U(K,I)*U(K,J)
      F=(S/U(I,I))/G
      DO 450 K=I,M
      U(K,J)=U(K,J)+F*U(K,I)
  450 CONTINUE
  460 DO 470 J=I,M
  470 U(J,I)=U(J,I)/G
      GO TO 490
  475 DO 480 J=I,M
  480 U(J,I)=0.
  490 U(I,I)=U(I,I)+1.
  500 CONTINUE
  510 DO 700 KK=1,N
      K1=N-KK
      K=K1+1
      ITS=0
  520 DO 530 LL=1,K
      L1=K-LL
      L=L1+1
      IF(ABS(RV1(L))+ANORM.EQ.ANORM)  GO TO 565
      IF(ABS(W(L1)) +ANORM.EQ.ANORM)  GO TO 540
  530 CONTINUE
  540 C=0.
      S=1.
      DO 560 I=L,K
      F=S*RV1(I)
      RV1(I)=C*RV1(I)
      IF(ABS(F)+ANORM.EQ.ANORM) GO TO 565
      G=W(I)
      H=SQRT(F*F+G*G)
      W(I)=H
      C= G/H
      S=-F/H
      IF(.NOT.MATU) GO TO 560
      DO 550 J=1,M
      Y=U(J,L1)
      Z=U(J,I)
      U(J,L1)=Y*C+Z*S
      U(J,I)=-Y*S+Z*C
  550 CONTINUE
  560 CONTINUE
  565 Z=W(K)
      IF( L .EQ. K) GO TO  650
      IF(ITS.EQ.300) GO TO 1000
      ITS=ITS+1
      X=W(L)
      Y=W(K1)
      G=RV1(K1)
      H=RV1(K)
      F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.*H*Y)
      G=SQRT(F*F+1.)
      F=((X-Z)*(X+Z)+H*(Y/(F+SIGN(G,F))-H))/X
C
      C=1.
      S=1.
      DO 600 I1=L,K1
      I=I1+1
      G=RV1(I)
      Y=W(I)
      H=S*G
      G=C*G
      Z=SQRT(F*F+H*H)
      RV1(I1)=Z
      C=F/Z
      S=H/Z
      F= X*C+G*S
      G=-X*S+G*C
      H=Y*S
      Y=Y*C
      IF(.NOT.MATV) GO TO 575
      DO 570 J=1,N
      X=V(J,I1)
      Z=V(J,I)
      V(J,I1)=X*C+Z*S
      V(J,I)=-X*S+Z*C
  570 CONTINUE
  575 Z=SQRT(F*F+H*H)
      W(I1)=Z
      IF(Z.EQ.0.)   GO TO 580
      C=F/Z
      S=H/Z
  580 F=C*G+S*Y
      X=-S*G+C*Y
      IF(.NOT.MATU) GO TO 600
      DO 590 J=1,M
      Y=U(J,I1)
      Z=U(J,I)
      U(J,I1)=Y*C+Z*S
      U(J,I)=-Y*S+Z*C
  590 CONTINUE
  600 CONTINUE
      RV1(L)=0.
      RV1(K)=F
      W(K)  =X
      GO TO 520
  650 IF(Z.GE.0.)   GO TO 700
      W(K)=-Z
      IF(.NOT.MATV) GO TO 700
      DO 690 J=1,N
  690 V(J,K)=-V(J,K)
  700 CONTINUE
      GO TO 1001
 1000 IERR=K
C
 1001 RETURN
      END
C	-------------------<POLLEG>---------------------
      SUBROUTINE POLLEG(NJ,NFMAX,NF,M,ARG,PL)
      REAL ARG(NJ),PL(NFMAX,NJ)
      NFM1=NF-1
      DO 10 J=1,NJ
      X=ARG(J)
      X=COS(X)
      XX=X*X
      IF(M.GT.1) GO TO 20
      XSQ=SQRT(1.-XX)
      PL(1,J)=XSQ
      PL(2,J)=3.*X*XSQ
      GO TO 30
   20 XSQ=1.-XX
      PL(1,J)=3.*XSQ
      PL(2,J)=15.*X*XSQ
   30 CONTINUE
      DO 40 IF=2,NFM1
      KF=IF-M+1
   40 PL(IF+1,J)=((2.*IF+1.)*X*PL(IF,J)-(IF+M)*PL(IF-1,J))/KF
   10 CONTINUE
      RETURN 
      END
c---------------------------------------                  
      SUBROUTINE  PRMV(NDIM,NI,NJ,A,B,C)
      REAL  A(NDIM,NJ), B(NJ), C(NI)
      DO  2  I = 1, NI
          SUM  = 0.
      DO  1  J = 1, NJ
          SUM  = SUM + A(I,J) * B(J)
   1  CONTINUE
          C(I) = SUM
   2  CONTINUE
      RETURN
      END
c--------------------------------------
      SUBROUTINE  PRVV(A,NI,B,RESULT)
      REAL  A(NI), B(NI)
        RESULT = 0.
      DO  1  I = 1, NI
        RESULT = RESULT + A(I) * B(I)
   1  CONTINUE
      RETURN
      END