C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- PROGRAM ZEOF_PROGRAM PARAMETER(N=2,M=4) REAL F(M,N),FTRANS(N,M),R(M,M),D(M),V(M,M),DUM(M),C(M,N) REAL VTRANS(M,M),VNEW(M,M),FNEW(M,N) REAL T1(M),T2(M) DATA T1 /1, 2, 5, -1/ DATA T2 /1, 3, 6, -2/ ! !Program: This program shows a simple example of Empirical !Orthogonal Function (EOF) analysis applied to temperature !anomalies (differences from climatological averages) measured !in degrees Celsius over four defined regions (i.e., M = 4) for !two consecutive days (i.e., N = 2). The subroutine invoked is !called "ZEOF". In EOF analysis there will always be M !eigenvalues and M eigenvectors, regardless of the number of time !series values (N). In this example, temperature anomalies shown !below do not change much between day 1 and day 2, and the !leading eigenvector explains more than 99 percent of the !time-dependent spatial variability: ! ! Day 1 Day 2 ! +----------+----------+ +----------+----------+ ! | | | | | | ! | 1C | 2C | | 1C | 3C | ! | | | | | | ! +----------+----------+ +----------+----------+ ! | | | | | | ! | 5C | -1C | | 6C | -2C | ! | | | | | | ! +----------+----------+ +----------+----------+ ! !In EOF analysis, what is generally plotted as spatial maps are !not EOFs but instead derived Empirical Orthogonal eigenVectors !(EOVs). Function time dependence (the EOFs) follows only after !the EOVs are multiplied by derived projection time series from !the EOF analysis (i.e., in this program eigenvector matrix "V" !multiplied by projection matrix "C" yields the time dependent !reconstructed data matrix "F"). The eigenvectors are mutually !orthogonal (actually, orthonormal in EOF analysis) and therefore !the inner (i.e., scalar "dot") product between any two of them !is always zero. ! !The purpose of EOF analysis is to determine basic state spatial !eigenvectors which when only a few are retained and multiplied !by the projection time series matrix C yields as close as !possible, the original gridded time series values. In general !applications you may have M=5000 grid points, each with !perhaps N=2000 time series values. There will then be 5000 !eigenvalues and 5000 eigenvectors. If you're lucky, you can !reproduce closely the original data matrix F (i.e., F = VC) by !retaining only the first few leading eigenvectors in the !eigenvector matrix V. Iterative Jacobi rotation is invoked !(i.e., Numerical Recipes subroutine "jacobi") to solve for !eigenvalues and eigenvectors. Computation time will become !lengthy in cases of many grid cells M and long time !series length N. ! !Here is a reference for the EOF method: ! !Kutzbach, J. E., Empirical eigenvectors of sea-level pressure, !surface temperature and precipitation complexes over North !America, J. App. Meteorol., 6, 791-802, 1967. ! !Author: Dr. Jerry R. Ziemke !Define values for the "data" matrix "F": DO I=1,M F(I,1)=T1(I) F(I,2)=T2(I) ENDDO !Conduct EOF analysis: CALL ZEOF(F,M,N,FTRANS,R,D,V,VTRANS,DUM,C) !Print out eigenvalues and percent variance explained (Note: D(1) !is the leading eigenvalue, D(2) the second leading eigenvalue, !etc.): SUM=0. DO I=1,M SUM=SUM+D(I) ENDDO WRITE(*,*)' ' WRITE(*,*)'Eigenvalue, and percent of total variance explained:' DO I=1,M WRITE(*,*) D(I), 100.*D(I)/SUM ENDDO !Print out leading and second leading eigenvector (Note: The !eigenvectors are ordered within matrix V - V(*,1) is the leading !eigenvector, V(*,2) is the second leading eigenvector, etc.): WRITE(*,*)'Leading eigenvector:' DO J=1,M WRITE(*,*) V(J,1) ENDDO WRITE(*,*)'Second leading eigenvector:' DO J=1,M WRITE(*,*) V(J,2) ENDDO !Print out original data matrix F: WRITE(*,*)'Original data matrix F:' WRITE(*,*) F !Analysis check: Eigenvector matrix V multiplied by projection !time series matrix C should be the original data matrix F: CALL ZMULT(V,C,F,M,M,N) WRITE(*,*)'Check: Reconstructed matrix F (should be original F):' WRITE(*,*) F !Print out the 2 time series values (N=2) at the 4 gridpoints !(M=4) for the leading EOV (this is the reconstructed data !matrix F for the first eigenvector): DO K=1,M DO J=1,M VNEW(J,K)=V(J,K) IF (K.GT.1) VNEW(J,K)=0. ENDDO ENDDO CALL ZMULT(VNEW,C,FNEW,M,M,N) WRITE(*,*)'Reconstructed data matrix F for first eigenvector:' WRITE(*,*) FNEW !Print out the 2 time series values (N=2) at the 4 gridpoints !(M=4) for the combined leading and second leading EOVs (this is !the reconstructed data matrix F for the combined first and !second eigenvectors): DO K=1,M DO J=1,M VNEW(J,K)=V(J,K) IF (K.GT.2) VNEW(J,K)=0. ENDDO ENDDO CALL ZMULT(VNEW,C,FNEW,M,M,N) WRITE(*,*)'Reconstructed data matrix F for combined first and' WRITE(*,*)'second leading eigenvectors:' WRITE(*,*) FNEW STOP END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- SUBROUTINE ZEOF(F,M,N,FTRANS,R,D,V,VTRANS,DUM,C) C Subroutine: Computes empirical orthogonal vectors (EOVs) and the time C series projection matrix for determining the empirical orthogonal C functions (EOFs) for "M" locations, each of time series length "N". C C Parameters, etc: C F: 2D array holding the M input time series, each of length N C M,N: Dimensions of array F C FTRANS : Transpose of array F C R : Symmetric matrix C D : Vector holding "M" ordered eigenvalues (one for each EOV) C V : M by M array holding the ordered eigenvectors (EOVs); the first C column is the leading EOV having the largest eigenvalue, etc. C VTRANS : Transpose of V C DUM : Work array C C : Projection matrix holding time series values. C C Overview of method: (See Kutzbach, J. Appl. Meteorol., Oct, 1967) C C One seeks to maximize C 1 N C - SUM [f(n) dot e]^2/(e*e) (1) C N n=1 C where "f(n)" are N series vectors with length M, "dot" denotes dot C product, and "e" denotes the EOVs to be determined (these are C determined in this routine by invoking a singular value decomposition C (SVD) routine from Numerical Recipes). C C F (M by N array): C F = f(1),f(2),f(3),..,f(N) C where f(N) are each column vectors of length M. C C Example: Singapore monthly zonal winds. Here, there are 7 (M=7) C pressure levels and 168 months (N=168) for 1979-1992. C F (M by N matrix): C C f(1) f(2) f(3) f(168) C C 10hPa ==> -31.1 -33.6 -32.0 . . . . . 16.9 C 15hPa ==> -32.4 -33.5 -34.2 . . . . . 18.4 C 20hPa ==> -26.6 -26.3 -27.6 . . . . . 18.9 C 30hPa ==> 3.7 6.5 -7.5 . . . . . 17.6 C 40hPa ==> 11.0 11.0 10.5 . . . . . 13.7 C 50hPa ==> 11.6 11.2 9.9 . . . . . 4.4 C 70hPa ==> 2.1 3.8 4.4 . . . . . -10.5 C C R (M by M matrix): R = F*FTRANS/N (This subroutine ZEOF solves C R*e = e*lambda which follows from maximizing (1) above via Lagrange C multiplier method for eigenvectors "e" and eigenvalues lambda.) C C C (M by N matrix): C = VTRANS*F where VTRANS is the transpose of C array V defined as the eigenvector array with each column holding C an EOV. In principle, only leading EOVs are needed to explain most C of the variability in the original data matrix F: C C MSTAR C f(I)_trunc = SUM C(J,I)*e(J) , where MSTAR < M C J=1 C Author: Jerry R. Ziemke REAL F(M,N),FTRANS(N,M),R(M,M),d(M),v(M,M),DUM(M),C(M,N) REAL vtrans(M,M) CALL Z0JRH4(F,FTRANS,M,N) CALL ZMULT(F,FTRANS,R,M,N,M) C WRITE(*,*)'Final R matrix (see Kutzbach):' DO I=1,M DO J=1,M R(I,J)=R(I,J)/N ENDDO ENDDO DO I=1,M DO J=1,M DUM(J)=R(I,J) ENDDO C WRITE(*,'(1X,7F10.2)') DUM ENDDO CALL JRUO24(R,M,M,d,v,nrot) !Calculate e-values & e-vectors CALL E3TJ29(d,v,M,M) !Sort: e-values now largest to smallest C WRITE(*,*)'Number of "rotations" in Jacobi subroutine:',nrot CALL Z0JRH4(v,vtrans,M,M) !<==Compute transpose of matrix v CALL ZMULT(vtrans,F,C,M,M,N) !<==Compute projection matrix C RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- SUBROUTINE ZMULT(A,B,CMULT,L,M,N) REAL A(L,M),B(M,N),CMULT(L,N) !Subroutine: Multiplies matrices A and B with output CMULT=A*B DO K=1,N DO I=1,L CMULT(I,K)=0. DO J=1,M CMULT(I,K)=CMULT(I,K)+A(I,J)*B(J,K) ENDDO ENDDO ENDDO RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- SUBROUTINE Z0JRH4(F,FTRANS,M,N) REAL F(M,N),FTRANS(N,M) !Subroutine: Yields transpose matrix FTRANS of matrix F DO I=1,N DO J=1,M FTRANS(I,J)=F(J,I) ENDDO ENDDO RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- SUBROUTINE JRUO24(a,n,np,d,v,nrot) INTEGER n,np,nrot,NMAX REAL a(np,np),d(np),v(np,np) PARAMETER (NMAX=500) INTEGER i,ip,iq,j REAL c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX) 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 !pause 'too many iterations in JRUO24 (NR subroutine "jacobi")' return END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- SUBROUTINE E3TJ29(d,v,n,np) INTEGER n,np REAL d(np),v(np,np) INTEGER i,j,k REAL p do 13 i=1,n-1 k=i p=d(i) do 11 j=i+1,n if(d(j).ge.p)then k=j p=d(j) endif 11 continue if(k.ne.i)then d(k)=d(i) d(i)=p do 12 j=1,n p=v(j,i) v(j,i)=v(j,k) v(j,k)=p 12 continue endif 13 continue return END