日志
Li Jianping老师的 REOF.f 子程序中错误已经修改好
热度 3 |
Li Jianping老师的REOF程序中错误已经修改好！
程序 REOF.f 如下，欢迎调用
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C ***********Li Jianping老师的 REOF 程序中错误已经修改好***********C
C A. STATISTICAL ANALYSIS C
C 8 EMPIRICAL ORTHOGONAL FUNCTIONS (EOF'S) AND C
C ROTATED EMPIRICAL ORTHOGONAL FUNCTIONS (REOF'S) C
C (1) REOF C
C (2) TRANSF C
C (3) CROSSPRODUCT C
C (4) JCB C
C (5) ARRANG C
C (6) TCOEFF C
C (7) ROTATOR C
C (8) ARRANGE2 C
C (9) MEANVAR C
C (10)INITIAL C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C BY DR. JIANPING LI AND DONG XIAO, JUL. 17, 2006.
C
!*****TIME-SPACE TRANSFERMATION IS USED FOR SAVING CALCULATIONAL TIME.
C THE MISSING VALUES (TERRAIN), ANNUAL CYCLE AND THE GRIDS WHOSE VARIANCE
C EQUAL TO ZERO ARE REMOVED (ADDED) BEFORE (AFTER) CALLING THE SUBROUTINE.
C
C REFERENCE: 《METHODS FOR DIAGNOSING AND FORECASTING CLIMATE VARIABILITY?
C BY HONGBAO WU AND LEI WU, PUBLISHED IN 2005. CHINA METOROLOGY PRESS.
C
C QUESTIONS, SUGGESTIONS OR GS FILE FOR DRAWING PICTURES? EMIAL US, PLEASE!
C E-MAIL: LJP@LASG.IAP.AC.CN, XIAODONG@MAIL.IAP.AC.CN.
C
C CITATION: LI, J., AND D. XIAO, 2006: REOF SUBROUTINE, HTTP://WEB.LASG.AC.CN/STAFF/LJP/SUBROUTINE/REOF.F.
C
C-----*----------------------------------------------------6---------7--
C
C WARNING: THE STACKS OF THIS SUBROUTINE MAY OVERFLOW IN PERSONAL COMPUTER (PC).
C IF YOU WANT TO COMPUTE IT ON PC, PLEASE DELETE THE EXPRESSIONS '*4'
C IN 68, 174, 181, 186 AND 193, AND GO ON.
C
C THE PARAMETERS YOU NEED TO CHANGE:
C N---LENGTH OF TIME SERIES, VIZ. SAMPLE SIZE;
C MX---GRIDS IN ROW;
C MY---GRIDS IN COLUMN;
C MNL---MIN(M,N)
C NP---NUM. OF EIGENVECTORS TO ROTATE;
C UNDEF---MISSING VALUE;
C MG1 AND MG2---THE EFFICIENT GRIDS (OUTPUT ON SCREEN) AND CANCEL THE
C "STOP" IN LINE 111, TRY AGAIN.
C KS=-1,0,1: SEE THE CAPTION OF THE REOF SUBROUTINE.
C KM=1: MONTHLY OR DAILY DATA. KM=0: DATA WITHOUT ANNUAL CYCLE.
C ND: THE NUMBER OF THE MONTHS OR DAYS IN A YEAR.
C
C-----*----------------------------------------------------6---------7--
PROGRAM MAIN
PARAMETER(N=722,MX=9,MY=9,M=MX*MY,MNL=M,NP=10)
PARAMETER(UNDEF=32767.0,MG1=81,MG2=81)
PARAMETER(KS=0,KM=1,ND=12,STD=0.0001)
!-----INPUT ARRAY
DIMENSION F0(N,M)
!-----WORK ARRAYS
DIMENSION F1(N,MG1),F2(N,MG1),G(MG1),H1(MG1,N),H2(MG1,N)
DIMENSION F(MG2,N),GVT(MG2,MNL),RGVT(MG2,NP),COF(MNL,N),RCOF(NP,N)
!-----OUTPUT ARRAYS
DIMENSION ER(MNL,4),EGVT(M,MNL),ECOF(MNL,N)
DIMENSION RER(NP,2),REGVT(M,NP),RECOF(NP,N)
C-----READ DATA.
OPEN(20,FILE='SST.DAT',STATUS='UNKNOWN'
&,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=M)
DO J=1,N
READ(20,REC=J)(F0(J,I),I=1,M)
END DO
CLOSE(20)
WRITE(*,*)'READ DATA OK'
C-----REMOVE THE TERRAIN OR MISSING VALUE.
L1=0
DO J=1,M
IF(F0(1,J).NE.UNDEF)THEN
L1=L1+1
DO K=1,N
F1(K,L1)=F0(K,J)
ENDDO
ENDIF
ENDDO
WRITE(*,*)'GRIDS WITHOUT TERRAIN'
WRITE(*,*)'MG1=',L1
C-----REMOVE ANNUAL CYCLE.
IF(KM.EQ.1)THEN
NY=N/ND
DO I=1,L1
CALL INITIAL(ND,NY,N,F1(1,I),F2(1,I))
END DO
DO I=1,L1
DO K=1,N
F1(K,I)=F2(K,I)
END DO
END DO
END IF
C-----REMOVE THE GRIDS WHOSE VARIANCE EQUAL TO ZERO.
L2=0
DO I=1,L1
CALL MEANVAR(N,F1(1,I),AX,G(I),VX)
IF(G(I).GT.STD)THEN
L2=L2+1
DO K=1,N
F(L2,K)=F1(K,I)
ENDDO
ENDIF
END DO
WRITE(*,*)'GRIDS WITHOUT TERRAIN AND CONSTANT VALUE'
WRITE(*,*)'MG2=',L2
STOP
C-----CALL THE SUBROUTINE.
CALL REOF(L2,N,MNL,NP,F,KS,ER,GVT,ECOF,RER,RGVT,RECOF)
WRITE(*,*)'REOF OK AND TRANSFORM TO THE ORIGINAL FORM IN THE NEXT'
C-----ADD THE GRIDS WHOSE VARIANCE EQUAL TO ZERO.
L3=0
DO I=1,MG2
IF(G(I).GT.STD)THEN
L3=L3+1
DO K=1,MNL
H1(I,K)=GVT(L3,K)
END DO
DO K=1,NP
H2(I,K)=RGVT(L3,K)
END DO
ELSE
DO K=1,MNL
H1(I,K)=UNDEF
END DO
DO K=1,NP
H2(I,K)=UNDEF
END DO
ENDIF
ENDDO
C-----ADD THE TERRAIN OR MISSING VALUE.
L4=0
DO I=1,M
IF(F0(1,I).NE.UNDEF)THEN
L4=L4+1
DO K=1,MNL
EGVT(I,K)=H1(L4,K)
END DO
DO K=1,NP
REGVT(I,K)=H2(L4,K)
END DO
ELSE
DO K=1,MNL
EGVT(I,K)=UNDEF
END DO
DO K=1,NP
REGVT(I,K)=UNDEF
END DO
ENDIF
ENDDO
C-----OUTPUT THE RESULT.
C-----OUTPUT THE ERROR.
WRITE(*,*)'OUTPUT THE RESULTS'
OPEN(10,FILE='ER.DAT',STATUS='UNKNOWN')
WRITE(10,*)'EOF LAMDA (EIGENVALUES) FROM BIG TO SMALL'
WRITE(10,*)(ER(I,1),I=1,MNL)
WRITE(10,*)'EOF ACCUMULATED EIGENVALUES FROM BIG TO SMALL'
WRITE(10,*)(ER(I,2),I=1,MNL)
WRITE(10,*)'EOF EXPLAINED VARIANCES'
WRITE(10,*)(ER(I,3),I=1,MNL)
WRITE(10,*)'EOF ACCUMULATED EXPLAINED VARIANCES'
WRITE(10,*)(ER(I,4),I=1,MNL)
WRITE(10,*)'REOF EXPLAINED VARIANCES'
WRITE(10,*)(RER(I,1)/N**2,I=1,NP)
WRITE(10,*)'REOF ACCUMULATED EXPLAINED VARIANCES'
WRITE(10,*)(RER(I,2)/N**2,I=1,NP)
CLOSE(10)
C-----OUTPUT EIGENVACTORS MATRIX OF EOF.
OPEN(11,FILE='EGVT.DAT',STATUS='UNKNOWN'
&,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=MX*MY)
DO J=1,MNL
WRITE(11,REC=J)(EGVT(I,J),I=1,M)
END DO
CLOSE(11)
C-----OUTPUT TIME COEFFICIENTS MATRIX OF EOF.
OPEN(12,FILE='ECOF.DAT',STATUS='UNKNOWN'
&,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=MNL)
DO I=1,N
WRITE(12,REC=I)(ECOF(J,I),J=1,MNL)
END DO
CLOSE(12)
C-----OUTPUT LOADING VECTORS OF REOF.
OPEN(13,FILE='REGVT.DAT',STATUS='UNKNOWN'
&,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=MX*MY)
DO J=1,NP
WRITE(13,REC=J)(REGVT(I,J),I=1,M)
END DO
CLOSE(13)
C-----OUTPUT TIME COEFFICIENTS MATRIX OF REOF.
OPEN(14,FILE='RECOF.DAT',STATUS='UNKNOWN'
&,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=NP)
DO J=1,N
WRITE(14,REC=J)(RECOF(I,J),I=1,NP)
END DO
CLOSE(14)
C-----
END
C-----*----------------------------------------------------6---------7--
C ROTATED EMPERICAL ORTHOGONAL FUNCTION (REOF)
C THIS SUBROUTINE APPLIES THE REOF APPROACH TO ANALYSIS
C TIME SERIES OF METEOROLOGICAL FIELD F(M,N).
C
C INPUT: M,N,MNL,NP,F(M,N),KS
C M: NUMBER OF GRID-POINTS
C N: LENGTH OF TIME SERIES, VIZ. SAMPLE SIZE
C MNL=MIN(M,N)
C NP: THE NUMBER OF THE EIGENVACTORS TO ROTATE
C (THE CORRESPONDING EOF ACCUMULATED VARIANCE IS ABOUT 60-80% OR THE
C CORRESPONDING EIGENVALUES OF EOF EXCEED 1.0 )
C F(M,N): THE RAW SPATIAL-TEMPORAL SEIRES
C KS: CONTRAL PARAMETER
C KS=-1: SELF; KS=0: DEPATURE; KS=1: STANDERLIZED DEPATURE
C SELF: GENERALLY NOT BE ADAPTED, THE FIRST MODES LIKE THE CLIMATE NORMAL.
C THE OTHERS ARE SAME AS THE MODES OF DEPARTURE FIELD.
C DEPARTURE: THE EIGENVECTORS STAND FOR THE AVERAGED ANOMALIES OF THE VARIABLE.
C THE LARGER THE ANOMALIES, THE MORE THE CONTRIBUTION.
C IT WOULD EMPHASIZE THE REGIONS WITH LARGER ANOMALIES.
C SUGGESTION: MAINLY USED IN LARGER SCALE FIELD.
C STANDARD: THE EIGENVECTORS STAND FOR THE CORRELATION COEFFICIENTS
C BETWEEN THE TIMESERIES OF THE GRIDS AND THE CORRESPONDING TIME COEFFICIENTS.
C STANDARD=DEPARTURE/(STANDARD DEVIATION). IT COMPUTE EVERY GRID IN SAME LEVEL.
C THERE IS NO REGIONAL DIFFERENCE IN THIS KIND FIELD.
C SUGGESTION: MAINLY USED IN THE CORRELATIONS OF DIFFERENT REGIONS.
C
C OUTPUT: EGVT,ECOF,ER,RER,REGVT,RECOF
C EGVT(M,MNL): ARRAY OF EIGENVACTORS
C ECOF(MNL,N): ARRAY OF TIME COEFFICIENTS FOR THE RESPECTIVE EIGENVECTORS
C ER(MNL,1): LAMDA (EIGENVALUES), ITS EQUENCE IS FROM BIG TO SMALL.
C ER(MNL,2): ACCUMULATED EIGENVALUES FROM BIG TO SMALL
C ER(MNL,3): EXPLAINED VARIANCES (LAMDA/TOTAL EXPLAIN) FROM BIG TO SMALL
C ER(MNL,4): ACCUMULATED EXPLAINED VARIANCES FROM BIG TO SMALL
C RER(MNL,1): EXPLAINED VARIANCES (LAMDA/TOTAL EXPLAIN)
C RER(MNL,2): ACCUMULATED EXPLAINED VARIANCES
C REGVT(M,NP): RODATED ARRAY OF LOADING VECTORS OF REOF
C RECOF(MNL,NP): RODATED ARRAY OF TIME COEFFICIENTS FOR THE RESPECTIVE EIGENVECTORS
C
C BY DONG XIAO AND DR. JIANPING LI, JUL. 17, 2006
C RECMOPLIED ON MAY 9, 2007
C
C REFERENCE: 《METHODS FOR DIAGNOSING AND FORECASTING CLIMATE VARIABILITY?
C BY HONGBAO WU AND LEI WU, PUBLISHED IN 2005. CHINA METOROLOGY PRESS.
C-----
SUBROUTINE REOF(M,N,MNL,NP,F,KS,ER,EGVT,ECOF,RER,REGVT,RECOF)
!-----INPUT ARRAY
DIMENSION F(M,N)
!-----WORK ARRAYS
DIMENSION COV(MNL,MNL),S(MNL,MNL),D(MNL)
!-----OUTPUT ARRAYS
DIMENSION ER(MNL,4),EGVT(M,MNL),ECOF(MNL,N)
DIMENSION RER(NP,2),REGVT(M,NP),RECOF(NP,N)
C---- PREPROCESSING DATA
CALL TRANSF(M,N,F,KS)
C---- CROSSED PRODUCT MATRIX OF THE DATA F
CALL CROSSPRODUCT(M,N,MNL,F,COV,SUM)
C---- EIGENVALUES AND EIGENVECTORS BY THE JACOBI METHOD
EPS=1E-7
CALL JCB(MNL,COV,S,D,EPS)
C---- SPECIFIED EIGENVALUES
CALL ARRANG(M,N,MNL,D,S,ER,TR)
WRITE(*,*)'ERROR=',SUM-TR
C---- NORMALIZED EIGNVECTORS AND THEIR TIME COEFFICIENTS
CALL TCOEFF(M,N,MNL,F,S,ER,EGVT,ECOF)
WRITE(*,*)'EOF ACCOMPLISHED'
C-----LINEAR TRANSPOSED (ROTATE)
CALL ROTAOR(M,N,MNL,NP,EGVT,ECOF,TR,RER,REGVT,RECOF)
WRITE(*,*)'REOF ACCOMPLISHED'
C---- SPECIFIED EIGNVECTORS AND TIME COEFFICIENTS WITH EXPLAINED VARIANCES OF REOF (IF YOU NEED)
CALL ARRANGE2(M,N,NP,REGVT,RECOF,RER)
C-----
RETURN
END
C-----*----------------------------------------------------6---------7--
C PREPROCESSING DATA TO PROVIDE A FIELD BY KS.
C
C INPUT: M,N,F,KS
C M: NUMBER OF GRID-POINTS
C N: LENGTH OF TIME SERIES, VIZ. SAMPLE SIZE
C F(M,N): THE RAW SPATIAL-TEMPORAL SEIRES
C KS: CONTRAL PARAMETER
C KS=-1: SELF; KS=0: DEPATURE; KS=1: STANDERDLIZED DEPATURE
C
C OUTPUT: F
C F(M,N): OUTPUT FIELD BASED ON THE CONTROL PARAMETER KS.
C
C BY DR. JIANPING LI.
C-----
SUBROUTINE TRANSF(M,N,F,KS)
!-----INPUT ARRAY
DIMENSION F(M,N)
!-----WORK ARRAYS
DIMENSION FW(N),WN(M)
I0=0
DO I=1,M
DO J=1,N
FW(J)=F(I,J)
ENDDO
CALL MEANVAR(N,FW,AF,SF,VF)
IF(SF.EQ.0.)THEN
I0=I0+1
WN(I0)=I
ENDIF
ENDDO
IF(I0.NE.0)THEN
WRITE(*,*)'**** FAULT ****'
WRITE(*,*)' THE PROGRAM CANNOT GO ON BECAUSE THE ORIGINAL FIELD'
WRITE(*,*)' HAS INVALID DATA WHOSE VARIANCE EQUAL ZERO.'
WRITE(*,*)
WRITE(*,*)' THERE ARE TOTALLY',I0,' GRIDPIONTS
*WITH INVALID DATA.'
WRITE(*,*)
WRITE(*,*)' THE ARRAY WN STORES THE POSITIONS OF THOSE INVALID '
WRITE(*,*)' GRID-POINTS. YOU MUST PICK OFF THOSE INVALID '
WRITE(*,*)' DATA FROM THE ORIGNAL FIELD AND THEN REINPUT '
WRITE(*,*)' A NEW FIELD TO CALCULATE ITS EOFS.'
WRITE(*,*)'**** FAULT ****'
STOP
ENDIF
IF(KS.EQ.-1)RETURN
C-----ANOMALY OF F
IF(KS.EQ.0)THEN
DO I=1,M
DO J=1,N
FW(J)=F(I,J)
ENDDO
CALL MEANVAR(N,FW,AF,SF,VF)
DO J=1,N
F(I,J)=F(I,J)-AF
ENDDO
ENDDO
RETURN
ENDIF
C-----NORMALIZING F
IF(KS.EQ.1)THEN
DO I=1,M
DO J=1,N
FW(J)=F(I,J)
ENDDO
CALL MEANVAR(N,FW,AF,SF,VF)
DO J=1,N
F(I,J)=(F(I,J)-AF)/SF
ENDDO
ENDDO
ENDIF
C-----
RETURN
END
C-----*----------------------------------------------------6---------7--
C CROSSED PRODUCT MARTIX OF THE DATA. IT IS N TIMES OF
C COVARIANCE MATRIX OF THE DATA IF KS=0 (I.E. FOR ANOMALY).
C
C INPUT: M,N,MNL,F
C M: NUMBER OF GRID-POINTS
C N: LENGTH OF TIME SERIES, VIZ. SAMPLE SIZE
C MNL=MIN(M,N)
C F(M,N): THE RAW SPATIAL-TEMPORAL SEIRES
C
C OUTPUT: COV(MNL,MNL),SUM
C COV(M,N)=(F*F') OR (F'F) DEPENDES ON M AND N.
C IT IS A MNL*MNL REAL SYMMETRIC MATRIX.
C SUM: THE TRA OF THE MATRIX (F*F')/N
C
C BY DR. JIANPING LI.
C RECOMPILED BY DONG XIAO, JUL. 17, 2006
C-----
SUBROUTINE CROSSPRODUCT(M,N,MNL,F,COV,SUM)
!-----INPUT ARRAY
DIMENSION F(M,N)
!-----WORK ARRAY
DIMENSION COV(MNL,MNL)
IF(N-M) 10,50,50
10 DO 20 I=1,MNL
DO 20 J=I,MNL
COV(I,J)=0.0
DO IS=1,M
COV(I,J)=COV(I,J)+F(IS,I)*F(IS,J)
ENDDO
COV(J,I)=COV(I,J)
20 CONTINUE
SUM=0.
DO I=1,MNL
SUM=SUM+COV(I,I)
END DO
SUM=SUM/(MNL*1.)
RETURN
50 DO 60 I=1,MNL
DO 60 J=I,MNL
COV(I,J)=0.0
DO JS=1,N
COV(I,J)=COV(I,J)+F(I,JS)*F(J,JS)
ENDDO
COV(J,I)=COV(I,J)
60 CONTINUE
SUM=0.
DO I=1,MNL
SUM=SUM+COV(I,I)
END DO
SUM=SUM/(MNL*1.)
C-----
RETURN
END
C-----*----------------------------------------------------6---------7--
C COMPUTING EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC MATRIX
C A(M,M) BY THE JACOBI METHOD
C
C INPUT: M,A,S,D,EPS
C M: ORDER OF MATRIX
C A(M,M): THE COVARIANCE MATRIX
C EPS: GIVEN PRECISION
C
C OUTPUT: S,D
C S(M,M): EIGENVECTORS
C D(M): EIGENVALUES
C
C BY DR. JIANPING LI.
C-----
SUBROUTINE JCB(M,A,S,D,EPS)
!-----INPUT ARRAY
DIMENSION A(M,M)
!-----OUTPUT ARRAYS
DIMENSION S(M,M),D(M)
DO 30 I=1,M
DO 30 J=1,I
IF(I-J) 20,10,20
10 S(I,J)=1.
GO TO 30
20 S(I,J)=0.
S(J,I)=0.
30 CONTINUE
G=0.
DO 40 I=2,M
I1=I-1
DO 40 J=1,I1
40 G=G+2.*A(I,J)*A(I,J)
S1=SQRT(G)
S2=EPS/FLOAT(M)*S1
S3=S1
L=0
50 S3=S3/FLOAT(M)
60 DO 130 IQ=2,M
IQ1=IQ-1
DO 130 IP=1,IQ1
IF(ABS(A(IP,IQ)).LT.S3) GOTO 130
L=1
V1=A(IP,IP)
V2=A(IP,IQ)
V3=A(IQ,IQ)
U=0.5*(V1-V3)
IF(U.EQ.0.0) G=1.
IF(ABS(U).GE.1E-10) G=-SIGN(1.,U)*V2/SQRT(V2*V2+U*U)
ST=G/SQRT(2.*(1.+SQRT(1.-G*G)))
CT=SQRT(1.-ST*ST)
DO 110 I=1,M
G=A(I,IP)*CT-A(I,IQ)*ST
A(I,IQ)=A(I,IP)*ST+A(I,IQ)*CT
A(I,IP)=G
G=S(I,IP)*CT-S(I,IQ)*ST
S(I,IQ)=S(I,IP)*ST+S(I,IQ)*CT
110 S(I,IP)=G
DO 120 I=1,M
A(IP,I)=A(I,IP)
120 A(IQ,I)=A(I,IQ)
G=2.*V2*ST*CT
A(IP,IP)=V1*CT*CT+V3*ST*ST-G
A(IQ,IQ)=V1*ST*ST+V3*CT*CT+G
A(IP,IQ)=(V1-V3)*ST*CT+V2*(CT*CT-ST*ST)
A(IQ,IP)=A(IP,IQ)
130 CONTINUE
IF(L-1) 150,140,150
140 L=0
GO TO 60
150 IF(S3.GT.S2) GOTO 50
DO 160 I=1,M
D(I)=A(I,I)
160 CONTINUE
C-----
RETURN
END
C-----*----------------------------------------------------6---------7--
C PROVIDES A SERIES OF EIGENVALUES FROM MAXIMUIM TO MINIMUIM.
C
C INPUT: M,N,MNL,D,S
C M: NUMBER OF GRID-POINTS
C N: LENGTH OF TIME SERIES, VIZ. SAMPLE SIZE
C D(MNL): EIGENVALUES
C S(MNL,MNL): EIGENVECTORS
C
C OUTPUT: ER,TR
C ER(MNL,1): LAMDA (EIGENVALUES), ITS EQUENCE IS FROM BIG TO SMALL.
C ER(MNL,2): ACCUMULATED EIGENVALUES FROM BIG TO SMALL
C ER(MNL,3): EXPLAINED VARIANCES (LAMDA/TOTAL EXPLAIN) FROM BIG TO SMALL
C ER(MNL,4): ACCUMULATED EXPLANED VARIANCES FROM BIG TO SMALL
C TR: THE TRA OF THE EIGENVECTORS MATRIX.
C
C BY DR. JIANPING LI.
C-----
SUBROUTINE ARRANG(M,N,MNL,D,S,ER,TR)
!-----INPUT ARRAYS
DIMENSION D(MNL),S(MNL,MNL)
!-----OUTPUT ARRAY
DIMENSION ER(MNL,4)
TR=0.0
DO 10 I=1,MNL
TR=TR+D(I)/(MNL*1.)
ER(I,1)=D(I)/(MNL*1.)
10 CONTINUE
MNL1=MNL-1
DO 20 K1=MNL1,1,-1
DO 20 K2=K1,MNL1
IF(ER(K2,1).LT.ER(K2+1,1)) THEN
C=ER(K2+1,1)
ER(K2+1,1)=ER(K2,1)
ER(K2,1)=C
DO 15 I=1,MNL
C=S(I,K2+1)
S(I,K2+1)=S(I,K2)
S(I,K2)=C
15 CONTINUE
ENDIF
20 CONTINUE
ER(1,2)=ER(1,1)
DO 30 I=2,MNL
ER(I,2)=ER(I-1,2)+ER(I,1)
30 CONTINUE
DO 40 I=1,MNL
ER(I,3)=ER(I,1)*100./TR
ER(I,4)=ER(I,2)*100./TR
40 CONTINUE
C-----
RETURN
END
C-----*----------------------------------------------------6---------7--
C PROVIDES EIGENVECTORS AND THEIR TIME COEFFICIENTS
C
C INPUT: M,N,MNL,F,S,ER
C M: NUMBER OF GRID-POINTS
C N: LENGTH OF TIME SERIES, VIZ. SAMPLE SIZE
C MNL=MIN(M,N)
C F(M,N): THE RAW SPATIAL-TEMPORAL SEIRES
C S(MNL,MNL): EIGENVECTORS
C ER(MNL,1): LAMDA (EIGENVALUES), ITS EQUENCE IS FROM BIG TO SMALL.
C ER(MNL,2): ACCUMULATED EIGENVALUES FROM BIG TO SMALL.
C ER(MNL,3): EXPLAINED VARIANCES (LAMDA/TOTAL EXPLAIN) FROM BIG TO SMALL.
C ER(MNL,4): ACCUMULATED EXPLANED VARIANCES FROM BIG TO SMALL.
C
C OUTPUT: EGVT,ECOF
C EGVT(M,MNL): NEW EIGENVECTORS
C DEPARTURE: THE EIGENVECTORS STAND FOR THE SWING OF THE VARIABLE.
C STANDARD: THE EIGENVECTORS STAND FOR THE CORRELATION COEFFICIENTS
C BETWEEN THE TIMESERIES AND THE CORRESPONDING TIME COEFFICIENTS.
C ECOF(MNL,N): TIME COEFFICIENTS OF EGVT
C
C BY DR. JIANPING.LI
C RECOMPILED BY DONG XIAO, JUL. 17, 2006
C-----
SUBROUTINE TCOEFF(M,N,MNL,F,S,ER,EGVT,ECOF)
!-----INPUT ARRAYS
DIMENSION F(M,N),S(MNL,MNL),ER(MNL,4)
!-----WORK ARRAYS
DIMENSION V(M),V1(N)
!-----OUTPUT ARRAYS
DIMENSION EGVT(M,MNL),ECOF(MNL,N)
DO J=1,MNL
DO I=1,M
EGVT(I,J)=0.
ENDDO
DO I=1,N
ECOF(J,I)=0.
ENDDO
ENDDO
C-----NORMALIZING THE INPUT EIGNVECTORS (C=SQRT(LAMDA))
DO 10 J=1,MNL
C=0.
DO I=1,MNL
C=C+S(I,J)*S(I,J)
ENDDO
C=SQRT(C)
DO I=1,MNL
S(I,J)=S(I,J)/C
ENDDO
10 CONTINUE
C-----(M<N)
IF(M.LE.N) THEN
DO JS=1,MNL
DO I=1,M
EGVT(I,JS)=S(I,JS)
ENDDO
ENDDO
DO 30 J=1,N
DO I=1,M
V(I)=F(I,J)
ENDDO
DO IS=1,MNL
DO I=1,M
ECOF(IS,J)=ECOF(IS,J)+V(I)*EGVT(I,IS)
ENDDO
ENDDO
30 CONTINUE
ELSE
C-----(M>N)
DO 40 I=1,M
DO J=1,N
V1(J)=F(I,J)
ENDDO
DO JS=1,MNL
DO J=1,N
EGVT(I,JS)=EGVT(I,JS)+V1(J)*S(J,JS)
ENDDO
ENDDO
40 CONTINUE
DO 60 J=1,N
DO I=1,M
V(I)=F(I,J)
ENDDO
!-----Z=V'(N*M)X(M*N)-------NOT BE ADOPTED, NEED TO BE NOMORALIZED AGAIN.
C DO IS=1,MNL
C DO I=1,N
C ECOF(IS,J)=ECOF(IS,J)+V(I)*EGVT(I,IS)
C ENDDO
C ENDDO
60 CONTINUE
C-----Z=V'SQRT(LAMDA)---(Z=V'X=>ZV=V'XV=V'(LAMDA*V)=(LAMDA*V')V=>Z=LAMDA*V' OR V=0)
C THE LAMDA HERE IS 1/N TIMES OF THAT IN THE REFERENCE.
DO 50 JS=1,MNL
DO J=1,N
ECOF(JS,J)=S(J,JS)*SQRT(ABS(ER(JS,1))/(MNL*1.0))
ENDDO
DO I=1,M
EGVT(I,JS)=EGVT(I,JS)/SQRT(ABS(ER(JS,1))/(MNL*1.0))
ENDDO
50 CONTINUE
C-----TRANSFER NORMALIZED EIGENVECTORS AND TIME COEFFICIENTS TO NEW ONES
C FOR NEW MEANS:
C DEPARTURE: THE EIGENVECTORS STAND FOR THE SWING OF THE VARIABLE.
C STANDARD: THE EIGENVECTORS STAND FOR THE CORRELATION COEFFICIENTS
C BETWEEN THE TIMESERIES AND THE CORRESPONDING TIME COEFFICIENTS.
DO JS=1,MNL
DO I=1,M
EGVT(I,JS)=EGVT(I,JS)*(SQRT(ABS(ER(JS,1))))
END DO
DO I=1,N
ECOF(JS,I)=(ECOF(JS,I))/(SQRT(ABS(ER(JS,1))))
END DO
END DO
END IF
C-----
RETURN
END
C-----*----------------------------------------------------6---------7--
C COMPUTING EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC MATRIX
C A(M,N) BY THE ROTATED METHOD
C
C INPUT: M,N,MNL,NP,EGVT,ECOF,TR
C M: GRIDS
C N: LENGTH OF TIME SERIES, VIZ. SAMPLE SIZE
C MNL=MIN(M,N)
C NP: NUM. OF EIGENVECTORS TO ROTATE
C EGVT(M,N): THE EIGENVACTORS MATRIX
C TR: THE RANK OF THE EIGENVECTORS MATRIX
C
C OUTPUT: REGVT,RECOF,RER
C REGVT(M,NP): THE ROTATED ARRAY OF EIGENVACTORS (LOADING VECTOR)
C RECOF(NP,N): ARRAY OF TIME COEFFICIENTS FOR THE RESPECTIVE EIGENVECTORS
C RER(MNL,1): EXPLAINED VARIANCES (LAMDA/TOTAL EXPLAIN)
C RER(MNL,2): ACCUMULATED EXPLAINED VARIANCES
C
C REFERENCE: 《METHODS FOR DIAGNOSING AND FORECASTING CLIMATE VARIABILITY?
C BY HONGBAO WU AND LEI WU, PUBLISHED IN 2005. CHINA METOROLOGY PRESS.
C
C BY PROF. HONGBAO WU
C BECOMPILED BY DONG XIAO, JULY 17, 2006
C-----
SUBROUTINE ROTAOR(M,N,MNL,NP,EGVT,ECOF,TR,RER,REGVT,RECOF)
!-----INPUT ARRAYS
DIMENSION EGVT(M,MNL),ECOF(MNL,N)
!-----WORK ARRAYS
DIMENSION H(M),ST(NP),VRLV(500),ER(MNL,4),S(M,NP)
!-----OUTPUT ARRAYS
DIMENSION REGVT(M,NP),RECOF(NP,N),RER(NP,2)
INTEGER T,K
WRITE(*,*)'ROTATE BEGINNING'
C-----TAKE OUT FORWARD NP EIGENVECTORS TO ROTATE
DO J=1,NP
DO I=1,M
REGVT(I,J)=EGVT(I,J)
END DO
END DO
C-----TAKE OUT FORWARD NP TIME COEFFICIENTS TO ROTATE
DO I=1,NP
DO J=1,N
RECOF(I,J)=ECOF(I,J)
END DO
END DO
C-----COMPUTE THE COMMON DEGREE
DO I=1,M
H(I)=0.0
DO J=1,NP
H(I)=H(I)+REGVT(I,J)**2
END DO
END DO
DO I=1,M
H(I)=SQRT(H(I))
END DO
C-----COMPUTE THE VARIANCE OF THE VARIANCE CONTRIBUTION BY PER COMMON DEGREE
C-----(EQUAL TO THE RATATED ONE)
VRLV0=0.0
DO K=1,NP
G1=0.0
G2=0.0
DO I=1,M
G1=G1+(REGVT(I,K)**2/H(I)**2)**2/REAL(M)
G2=G2+(REGVT(I,K)**2/H(I)**2)/REAL(M)
END DO
G2=G2**2
VRLV0=VRLV0+G1-G2
END DO
C-----ROTATED
LLL=0
WRITE(*,*)'ROTATE TIMES'
800 CONTINUE
DO 505 T=1,NP
DO 505 K=1,NP
IF(T.GE.K) GOTO 505
CALL ROT(REGVT,RECOF,H,T,K,M,N,NP)
505 CONTINUE
LLL=LLL+1
WRITE(*,*)'LLL=',LLL
VRLV(LLL)=0.0
DO K=1,NP
G1=0.0
G2=0.0
DO I=1,M
G1=G1+(REGVT(I,K)**2/H(I)**2)**2/REAL(M)
G2=G2+(REGVT(I,K)**2/H(I)**2)/REAL(M)
END DO
G2=G2**2
VRLV(LLL)=VRLV(LLL)+G1-G2
END DO
IF(LLL.LT.2) GOTO 800
CI=(VRLV(LLL)-VRLV(LLL-1))/VRLV0
IF(CI.LT.0.001) GOTO 600
IF(LLL.GE.100) GOTO 600
GOTO 800
600 CONTINUE
C-----COMPUTE THE EXPLAINED VARIANCES
DO J=1,NP
ST(J)=0.0
DO I=1,M
ST(J)=ST(J)+REGVT(I,J)**2
END DO
END DO
DO J=1,NP
RER(J,1)=ST(J)*100/TR
END DO
C-----COMPUTE THE ACCUMULATED EXPLAINED VARIANCES
DDD=0.0
DO K=1,NP
DDD=DDD+ST(K)/TR
RER(K,2)=DDD*100.
ENDDO
C-----
RETURN
END
C-----*----------------------------------------------------6---------7--
C COMPUTING EIGENVALUES AND EIGENVECTORS OF TWO COLUMNS OF
C A(M,T) AND A(M,K) BY THE ROTATED METHOD
C
C INPUT: A,B,H,T,K,MM,NN,NP
C A(M,NP): THE COVARIANCE MATRIX
C B(NP,N): THE TIME COEFFICIENT MATRIX
C H(M): COMMON DEGREE.
C T: COLUMN 1 (GRID 1)
C K: COLUMN 2 (GRID 2)
C M: GRIDS
C N: LENGTH OF TIME SERIES, VIZ. SAMPLE SIZE
C NP: NUM. OF EIGENVECTORS TO ROTATE
C
C OUTPUT: A,B
C A(M,NP): THE ROTATED COVARIANCE MATRIX
C B(NP,N): TIME COEFFICIENT MATRIX
C
C BY PROF. HONGBAO WU
C BECOMPILED BY DONG XIAO, JUL 17, 2006
C-----
SUBROUTINE ROT(A,B,H,T,K,M,N,NP)
!-----INPUT ARRAYS
DIMENSION A(M,NP),B(NP,N),H(M)
!-----WORK ARRAYS
DIMENSION U(M),VR(M),WT(M),WK(M),WBT(N),WBK(N)
!-----OUTPUT ARRAYS
C DIMENSION A(M,NP),B(NP,N)
INTEGER T,K
C-----COMPUTE FI
DO 25 I=1,M
U(I)=(A(I,T)/H(I))**2-(A(I,K)/H(I))**2
25 VR(I)=2.0*(A(I,T)/H(I))*(A(I,K)/H(I))
C=0.0
D=0.0
S=0.0
G=0.0
DO 30 I=1,M
C=C+(U(I)**2-VR(I)**2)
D=D+2.0*U(I)*VR(I)
S=S+U(I)
30 G=G+VR(I)
TG1=D-2.0*S*G/(M*1.)
TG2=C-(S**2-G**2)/(M*1.)
FI=ATAN2(TG1,TG2)/4.0
C-----COMPUTE NEW A AND B WITH FI
DO 75 I=1,M
WT(I)=A(I,T)
75 WK(I)=A(I,K)
DO 84 KK=1,N
WBT(KK)=B(T,KK)
84 WBK(KK)=B(K,KK)
DO 78 I=1,M
A(I,T)=WT(I)*COS(FI)+WK(I)*SIN(FI)
78 A(I,K)=WT(I)*(-SIN(FI))+WK(I)*COS(FI)
DO 89 IT=1,N
B(T,IT)=WBT(IT)*COS(FI)+WBK(IT)*SIN(FI)
89 B(K,IT)=WBT(IT)*(-SIN(FI))+WBK(IT)*COS(FI)
C-----
RETURN
END
C-----*----------------------------------------------------6--------7--
C SPECIFIED EIGNVECTORS AND TIME COEFFICIENTS WITH EXPLAINED VARIANCES
C OF REOF (IF YOU NEED). IF NOT, THE LOADING VECTORS OF REOF
C CORRESPONDS TO THE EOF MODES.
C
C INPUT: M,N,NP,REGVT,RECOF,TR
C M: GRIDS
C N: LENGTH OF TIME SERIES, VIZ. SAMPLE SIZE
C NP: NUM. OF EIGENVECTORS TO ROTATE
C REGVT(M,NP): THE ROTATED ARRAY OF EIGENVACTORS (LOADING VECTOR)
C RECOF(NP,N): ARRAY OF TIME COEFFICIENTS FOR THE RESPECTIVE EIGENVECTORS
C RER(MNL,1): EXPLAINED VARIANCES (LAMDA/TOTAL EXPLAIN)
C RER(MNL,2): ACCUMULATED EXPLAINED VARIANCES
C
C OUTPUT: REGVT,RECOF,RER
C REGVT(M,NP): THE ROTATED ARRAY OF EIGENVACTORS (LOADING VECTOR)
C RECOF(NP,N): ARRAY OF TIME COEFFICIENTS FOR THE RESPECTIVE EIGENVECTORS
C RER(MNL,1): EXPLAINED VARIANCES (LAMDA/TOTAL EXPLAIN)
C RER(MNL,2): ACCUMULATED EXPLAINED VARIANCES
C
C BY DONG XIAO, AUG 21, 2006
C-----
SUBROUTINE ARRANGE2(M,N,NP,REGVT,RECOF,RER)
!-----INPUT AND OUTPUT ARRAYS
DIMENSION REGVT(M,NP),RECOF(NP,N),RER(NP,2)
MNL1=NP-1
DO 20 K1=MNL1,1,-1
DO 20 K2=K1,MNL1
IF(RER(K2,1).LT.RER(K2+1,1)) THEN
C=RER(K2+1,1)
RER(K2+1,1)=RER(K2,1)
RER(K2,1)=C
DO I=1,M
D=REGVT(I,K2+1)
REGVT(I,K2+1)=REGVT(I,K2)
REGVT(I,K2)=D
END DO
DO I=1,N
E=RECOF(K2+1,I)
RECOF(K2+1,I)=RECOF(K2,I)
RECOF(K2,I)=E
END DO
ENDIF
20 CONTINUE
DO I=1,NP
RER(I,2)=0.0
END DO
RER(1,2)=RER(1,1)
DO 30 I=2,NP
RER(I,2)=RER(I-1,2)+RER(I,1)
30 CONTINUE
C-----
RETURN
END
C-----*----------------------------------------------------6---------7--
C COMPUTING THE MEAN AX, STANDARD DEVIATION SX
C AND VARIANCE VX OF A SERIES X(I) (I=1,...,N).
C
C INPUT: N AND X(N)
C N: LENGTH OF RAW SERIES
C X(N): RAW SERIES
C
C OUTPUT: AX, SX AND VX
C AX: THE MEAN VALUE OF X(N)
C SX: THE STANDARD DEVIATION OF X(N)
C VX: THE VARIANCE OF X(N)
C
C BY DR. JIANPING LI, MAY 6, 1999.
C-----
SUBROUTINE MEANVAR(N,X,AX,SX,VX)
!-----INPUT ARRAY
DIMENSION X(N)
AX=0.
VX=0.
SX=0.
DO 10 I=1,N
AX=AX+X(I)
10 CONTINUE
AX=AX/FLOAT(N)
DO 20 I=1,N
VX=VX+(X(I)-AX)**2
20 CONTINUE
VX=VX/FLOAT(N)
SX=SQRT(VX)
C-----
RETURN
END
C-----*----------------------------------------------------6---------7--
C THIS SUBROUTINE IS FOR REMOVING ANNUAL CYCLE (MONTHLY OF DAILY DATA)
C
C INPUT:
C ND: NUMBER OF DAYS OR MONTHES IN A YEAR.
C NY: NUMBER OF YEAR
C NDY=ND*NY
C F(NDY): THE DAILY OF MONTHLY DATA.
C
C OUTPUT:
C AF(NDY): THE DAILY OF MONTHLY DATA WITHOUT ANNUAL CYCLE
C
C BY DR JIANPING LI,
C RECOMPLIED BY DONG XIAO ON MAY 7, 2007
C-----
SUBROUTINE INITIAL(ND,NY,NDY,F,AF)
!-----INPUT ARRAY
REAL F(NDY)
!-----WORK ARRAYS
REAL FAVE(ND),FSTD(ND)
!-----OUTPUT ARRAY
REAL AF(NDY)
C-----
DO I=1,ND
FAVE(I)=0.
DO J=1,NY
FAVE(I)=FAVE(I)+F(I+(J-1)*ND)
ENDDO
FAVE(I)=FAVE(I)/NY
FSTD(I)=0.0
DO J=1,NY
FSTD(I)=FSTD(I)+(F(I+(J-1)*ND)-FAVE(I))**2
ENDDO
FSTD(I)=SQRT(FSTD(I)/FLOAT(NY))
DO J=1,NY
IJ=I+(J-1)*ND
C-----STANDARD DEPARTURE
C AF(IJ)=(F(IJ)-FAVE(I))/FSTD(I)
C-----DEPARTURE
AF(IJ)=(F(IJ)-FAVE(I))
ENDDO
ENDDO
C-----
RETURN
END
C-----*----------------------------------------------------6---------7--
不会这个，赶紧尝试！谢谢！