立即注册 登录
气象家园 返回首页

/xin村儿/的个人空间 http://bbs.06climate.com/?3576 [收藏] [复制] [分享] [RSS]

日志

Li Jianping老师的 REOF.f 子程序中错误已经修改好

热度 3已有 1800 次阅读2013-12-16 16:43 | FORTRAN子程序


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--


发表评论 评论 (3 个评论)

回复 hillside 2013-12-18 18:07
好!试试。
回复 yyy 2013-12-19 11:55
   不会这个,赶紧尝试!谢谢!
回复 豆角菲菲 2013-12-24 15:47
多谢分享

facelist doodle 涂鸦板

您需要登录后才可以评论 登录 | 立即注册

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

返回顶部