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