经验正交函数

  气象统计分析与预报方法    课程实验报告
实验名称
经验正交函数分解
 
大气科学woc
 
马兴邦
   
2012011009
 
大气121
实验地点
科教楼
实验日期
2014/10/21
 
指导老师
肖天贵
同组其他成员
第一组
一、 实验内容(含实验原理介绍):
(1) 对(- -850hPa高度场进行经验正交展开(EOF.FOR),输出分析主要参数指标;绘制环流型图和相应的时间系数序列图,并加以分析。
(2) 提供的资料为是20137100时到73118时逐6小时的850hPa高度场,时间序列长度为124次,网格点为1*1;经纬度范围:0N-90N,0E-180E
3EOF经验正交函数分解)是针对气象要素场进行的,其基本原理是把包含p个空间点
(变量)的场随时间变化进行分解。设抽取样本容量为n的资料.则场中任一空间点i和任一时间点j的距平观测值可看成由p个空间函数和时间函数(k=1,2,,p)的线性组合,表示成
EOF功能是从一个气象场多次观测资料中识别出主要空间型及其时间演变规律。
EOF展开就是将气象变量场分解为空间函数(V)和时间函数(T)两部分的乘积之和:
X=VT
二、实验目的:
经验正交函数分解(EOF)是统计天气分析中气象要素场最基础的研究模型,是必须理解和掌握的方法之一,是后续课程中许多气象要素场的计算结果的理解的基础理论,也是毕业设计和论文中的基本分析方法。该方法用个数较少的几个空间分布模态来描述环流形势,而且基本涵盖环流场的信息,既能作为天气分析模型,其方法的延拓又能作为天气预报模型,在实际工作中也有极强的实用意义。通过该实验,深刻理解气象要素场的统计模型的意义,掌握气象要素场分析的基本方法,为实际预报业务和科研工作打下一定的基础。
三、涉及实验的相关情况介绍(包含使用软件或实验设备等情况):
1)使用的软件是fortran5.0二甲基乙醇胺
四、实验结果(含程序、数据记录及分析和实验总结等,可附页):
一、程序
C**********************************************************************
C                                                                    *
C                        PROGRAM NOTES                                *
C                                                                    *
C          THIS PROGRAM USES EOF TO ANALYSIS TIME SERIES              *
C                    OF METEOROLOGICAL FIELD                          *
C                                                                    *
C**********************************************************************
C                                                                    *
C              ******** Parameter Table *********                    *
C                                                                    *
C    Mt===>LENTH OF TIME SERIES                                      *
C    N ===>NUMBER OF GRID-POINTS ( or STATIONS )                    *
C    KS=-1, SELF;  KS=0, DEPATURE;  KS=1, STANDERDLIZED DEPATURE  *
C    KV = NUMBER OF EIGENVALUES WILL BE OUTPUT                      *
C    KVT = NUMBER OF EIGENVECTORS AND TIME SERIES WILL BE OUTPUT    *
C    MNH = Minimum(Mt,N)                                            *
C    EGVT===>EIGENVECTORS, ECOF===>TIME COEFFICIENTS FOR EGVT        *
C    ER(KV,1)====>LAMDA;    LAMDA===>EIGENVALUE                      *
C    ER(KV,2)====>ACCUMULATE LAMDA                                  *
C    ER(KV,3)====>THE SUM OF COMPONENTS VECTORS PROJECTED ONTO      *
C                  EIGENVACTOR.                                      *
C    ER(KV,4)====>ACCUMULATE ER(KV,3)                                *
C                                                                    *
C**********************************************************************
      PARAMETER(N=??, MT=??, MNH=??)
      PARAMETER(KS=1, KV=5, KVT=5)   
      REAL F(N,MT),AVF(N),DF(N),ER(MNH,4)
      REAL A(MNH,MNH),S(MNH,MNH),V(MNH)
c**************************************************************************************
c    INFN-输入数据文件名;OUTERA-输出特征值及方差贡献、累积方差贡献的文件名(文本);
c    OUTTC1-输出时间系数文件(文本);OUTTC2-输出时间系数文件(二进制);
c    OUTTEVT-输出特征向量文件(二进制);
c**************************************************************************************   
    CHARACTER*50 INFN,OUTERA,OUTTC1,OUTTC2,OUTEVT
    DATA INFN/'??.dat'/
    DATA OUTERA/''/
    DATA OUTTC1/''/
    DATA OUTTC2/'hgt_OUTTC2.DAT'/
    DATA OUTEVT/'hgt_OUTTEVT.DAT'/细胞凋亡
   
C---------------- Read ORIGINAL DATA ----------------------------
      write(*,*)'Now is reading primative field !'
      OPEN (8,FILE=INFN,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=N)
      DO IT=1,MT
        READ (8,REC=IT)(F(IS,IT),IS=1,N)
      END DO
      pause
   
C************** START TO RUN EOF PROGRAM ******************************
      WRITE(*,*)
      write(*,*)' FIRST STEP'
      write(*,*)' forming the initial matrix (F) by using TRANSF !'
      CALL TRANSF(N,Mt,F,AVF,DF,KS)
      WRITE(*,*)
      write(*,*)' STEP 2'
      write(*,*)' achiving the covariance matrix by using the FORMA !'
      CALL FORMA(N,Mt,MNH,F,A)
      WRITE(*,*)
      write(*,*)' STEP 3 '
      write(*,*)' caculating the eigenvalue and eigenvectors '
      WRITE(*,*)' by using Jacob method !'
      CALL JCB(MNH,A,S,0.001)
      WRITE(*,*)
      write(*,*)' STEP 4'
      write(*,*)' arrange the eigenvalue and eigenvectors'
    WRITE(*,*)' by using ARRANG !'
      CALL ARRANG(MNH,A,ER,S)
      WRITE(*,*)
      write(*,*)' STEP 5'
      write(*,*)' the caculation of standard eigenvectors'
      WRITE(*,*)' by using TCOEFF !'
    CALL TCOEFF(KVT,N,Mt,MNH,S,F,V,ER)
      write(*,*)
      write(*,*)' STEP 6'
      write(*,*)' outputing eigenvalue and accumulation using OUTER !'
      CALL OUTER(MNH,ER,OUTERA)
      WRITE(*,*)
      WRITE(*,*)' STEP 7'
      write(*,*)' outputing the time coefficient of the eigenvecters !'
      CALL OUTVT1(KVT,N,Mt,MNH,S,F,OUTTC1,OUTTC2)
      WRITE(*,*)
      WRITE(*,*)' STEP 8'
      write(*,*)' outputing the eigenvecters !'
    CALL OUTVT3(KVT,N,Mt,MNH,S,F,OUTEVT)
      END
C *************** FINISH THE MAIN PROGRAM *****************************
C                                                                    *
C                    SUBROUTINE FUNCTION                            *
C                                                                    *
C    THIS SUBROUTINE PRINTS ARRAY ER                                *
C    ER(KV,1) FOR  SEQUENCE OF EIGENVALUE FROM BIG TO SMALL          *
C    ER(KV,2) FOR  EIGENVALUE FROM BIG TO SMALL                      *
C    ER(KV,3) FOR  SMALL LO=(LAMDA/TOTAL VARIANCE)                  *
C    ER(KV,4) FOR  BIG LO=SUM OF SMALL LO)                          *
C *********************************************************************
C -------- SAVING THE EIGENVALUE AND ERROR ---------------------------*
      SUBROUTINE OUTER(MNH,ER,OUTERA)
      DIMENSION ER(MNH,4)
    CHARACTER*50 OUTERA
核酸外切酶
      open (30,file=OUTERA)
      WRITE(30,510)
      WRITE(30,520)
      WRITE(30,530) (IS,(ER(IS,J),J=1,4),IS=1,MNH)
      CLOSE(30)
  510 FORMAT(25X,'EIGENVALUE AND ANALYSIS ERROR')
  520 FORMAT(5X,'N',8X,'LAMDA',10X,'SLAMDA',11X,'PH',12X,'SPH')
  530 FORMAT(I6,2E15.6,2F15.5)
    RETURN
      END
C**********************************************************************
C                        SUBROUTINE FUNCTION                          *
C                                                                    *
C              THIS SUBROUTINE PRINTS STANDARD EIGENVECTORS          *
C                AND ITS TIME-COEFFICENT SERIES                      *
C**********************************************************************
C ------------- save time-coeffivcent seried of S.E. ------------
      SUBROUTINE OUTVT1(KVT,N,M,MNH,S,F,OUTTC1,OUTTC2)
      DIMENSION F(N,M),S(MNH,MNH)
    CHARACTER*50 OUTTC1,OUTTC2
      OPEN(31,file=OUTTC1)
      OPEN(32,file=OUTTC2,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=KVT)
      WRITE(31,400)
      WRITE(31,200) (IS,IS=1,KVT)
      DO J=1,M
        IF(M.GE.N) THEN
          WRITE(31,300) J,(F(IS,J),IS=1,KVT)
          WRITE(32,REC=J) (F(IS,J),IS=1,KVT)
        ELSE
          WRITE(31,300) J,(S(J,IS),IS=1,KVT)
          WRITE(32,REC=J) (s(J,IS),IS=1,KVT)
        ENDIF
    END DO
      CLOSE(31)
  200 FORMAT(3X,10I15)
  300 FORMAT(I5,10E15.7)
  400 FORMAT(30X,'TIME-COEFFICENT SERIES OF S. E.')
      RETURN
      END
C --------- save standard eignvectors ------------------
      SUBROUTINE OUTVT3(KVT,N,M,MNH,S,F,OUTEVT)
      DIMENSION F(N,M),S(MNH,MNH)
    CHARACTER*50 OUTEVT
      OPEN(33,file=OUTEVT,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=N)
      DO JS=1,KVT
        IF(M.GE.N) THEN
          WRITE(33,REC=JS)(S(I,JS),I=1,N)
        ELSE
          WRITE(33,REC=JS)(F(I,JS),I=1,N)
        ENDIF
    END DO
      CLOSE(33)
    RETURN
    END
C***********************************************************************
C                    SUBROUTINE FUNCTION                              *
C                                                                      *
gis技术
C    THIS SUBROUTINE PROVIDES INITIAL F BY KS (optional parameter)    *
C      ks=-1, 0, or 1 according to primative field                    *
C***********************************************************************
      SUBROUTINE TRANSF(N,M,F,AVF,DF,KS)
      REAL F(N,M),AVF(N),DF(N)
      IF (KS)  30,10,10
10    DO I=1,N
        AVF(I)=0.0
        DF(I)=0.0
      END DO
      DO I=1,N
        DO J=1,M
        AVF(I)=AVF(I)+F(I,J)
        END DO
        AVF(I)=AVF(I)/M
        DO J=1,M
          F(I,J)=F(I,J)-AVF(I)
        END DO
      END DO
      IF (KS.EQ.1) THEN
        DO I=1,N
          DO J=1,M
            DF(I)=DF(I)+F(I,J)*F(I,J)
          END DO
          DF(I)=SQRT(DF(I)/M)
          DO J=1,M
            F(I,J)=F(I,J)/DF(I)
        END DO
        END DO
      END IF
30  CONTINUE
    RETURN
      END
C ----------------- FORMA -----------------------------
      SUBROUTINE FORMA(N,M,MNH,F,A)
      REAL F(N,M),A(MNH,MNH)
      IF (M-N) 40,50,50
40  DO I=1,MNH
        DO J=1,I
          A(I,J)=0.0
          DO IS=1,N
            A(I,J)=A(I,J)+F(IS,I)*F(IS,J)
          END DO
          A(J,I)=A(I,J)
      END DO
    END DO
      RETURN
50  DO I=1,MNH
        DO J=1,I
          A(I,J)=0.0
          DO JS=1,M
            A(I,J)=A(I,J)+F(I,JS)*F(J,JS)
          END DO
          A(J,I)=A(I,J)
      END DO
    END DO
    RETURN
      END
c***********************************************************************
C                    SUBROUTINE FUNCTION                              *
vip客户管理C                                                                      *
C    THIS SUBROUTINE COMPUTS EIGENVALUES AND EIGENVECTORS OF A        *
c***********************************************************************
      SUBROUTINE JCB(N,A,S,EPS)
      DIMENSION A(N,N),S(N,N)
      DO I=1,N
        DO J=1,N
          IF (I.EQ.J) THEN
            S(I,J)=1.0
        ELSE
          S(I,J)=0.0
        END IF
        END DO
      END DO
      G=0.0
      DO I=2,N
        I1=I-1
        DO J=1,I1
          G=G+2.*A(I,J)*A(I,J)
      END DO
    END DO
      S1=SQRT(G)
      S2=EPS/FLOAT(N)*S1
      S3=S1
      L=0
50  S3=S3/FLOAT(N)
60  DO 130 IQ=2,N
        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 I=1,N
            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
            S(I,IP)=G
          END DO
          DO I=1,N
            A(IP,I)=A(I,IP)
            A(IQ,I)=A(I,IQ)
        END DO
          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
      GOTO 60
150  IF (S3.GT.S2) GOTO 50
    RETURN
      END
c**********************************************************************
C                      SUBROUTINE FUNCTION                            *
C                                                                    *
C          THIS SUBROUTINE PROVIDES A SERIES OF EIGENVALUES          *
C                        FROM MAX TO MIN                              *
C**********************************************************************
      SUBROUTINE ARRANG(MNH,A,ER,S)
      DIMENSION A(MNH,MNH),ER(MNH,4),S(MNH,MNH)
      TR=0.0
      DO I=1,MNH
        TR=TR+A(I,I)
      ER(I,1)=A(I,I)
    END DO
      MNH1=MNH-1
      DO K1=MNH1,1,-1
        DO K2=K1,MNH1
          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 I=1,MNH
              C=S(I,K2+1)
              S(I,K2+1)=S(I,K2)
              S(I,K2)=C
          END DO
          END IF
      END DO
    END DO
      ER(1,2)=ER(1,1)
      DO I=2,MNH
        ER(I,2)=ER(I-1,2)+ER(I,1)
    END DO
      DO I=1,MNH
        ER(I,3)=ER(I,1)/TR
        ER(I,4)=ER(I,2)/TR
    END DO
    RETURN
      END
C**********************************************************************
C          THIS SUBROUTINE PROVIDES STANDARD EIGENVECTORS            *
C (M.GE.N, SAVED IN S;  M.LT.N, SAVED IN F) AND ITS TIME COEFFICENTS *
C        SERIES (M.GE.N, SAVED IN F;    M.LT.N, SAVED IN S)          *
C**********************************************************************
      SUBROUTINE TCOEFF(KVT,N,M,MNH,S,F,V,ER)
      DIMENSION S(MNH,MNH),F(N,M),V(MNH),ER(MNH,4)
      DO J=1,MNH
        C=0.0
        DO I=1,MNH
        C=C+S(I,J)*S(I,J)
      END DO
        C=SQRT(C)
        DO I=1,MNH
        S(I,J)=S(I,J)/C
      END DO
    END DO
      IF (M.GE.N) THEN
        DO J=1,M
          DO I=1,N
            V(I)=F(I,J)
            F(I,J)=0.0
        END DO
          DO IS=1,KVT
            DO I=1,N
            F(IS,J)=F(IS,J)+V(I)*S(I,IS)
          END DO
        END DO
      END DO
      ELSE
        DO I=1,N
          DO J=1,M
            V(J)=F(I,J)
            F(I,J)=0.0
        END DO
          DO JS=1,KVT
            DO J=1,M
              F(I,JS)=F(I,JS)+V(J)*S(J,JS)
          END DO
        END DO
      END DO
        DO JS=1,KVT
          DO J=1,M
            S(J,JS)=S(J,JS)*SQRT(ER(JS,1))
        END DO
          DO I=1,N
            F(I,JS)=F(I,JS)/SQRT(ER(JS,1))
        END DO
      END DO
      END IF
    RETURN
      END
二、 数据记录与图像输出
EIGENVALUE AND ANALYSIS ERROR
    N        LAMDA          SLAMDA          PH            SPH
    1  0.403676E+06  0.403676E+06        0.19765        0.19765
    2  0.265756E+06  0.669432E+06        0.13012        0.32777
    3  0.254969E+06  0.924401E+06        0.12484        0.45260
    4  0.172044E+06  0.109644E+07        0.08424        0.53684
    5  0.141533E+06  0.123798E+07        0.06930        0.60614
    6  0.134411E+06  0.137239E+07        0.06581        0.67195
    7  0.107424E+06  0.147981E+07        0.05260        0.72455
    8  0.868768E+05  0.156669E+07        0.04254        0.76708
    9  0.699882E+05  0.163668E+07        0.03427        0.80135
    10  0.551706E+05  0.169185E+07        0.02701        0.82836
    11  0.442305E+05  0.173608E+07        0.02166        0.85002
    12  0.433148E+05  0.177939E+07        0.02121        0.87123
三、 图像结果
   

本文发布于:2024-09-24 02:25:50,感谢您对本站的认可!

本文链接:https://www.17tex.com/xueshu/178496.html

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系,我们将在24小时内删除。

标签:实验   时间   气象要素   输出   分析   分解   方法   环流
留言与评论(共有 0 条评论)
   
验证码:
Copyright ©2019-2024 Comsenz Inc.Powered by © 易纺专利技术学习网 豫ICP备2022007602号 豫公网安备41160202000603 站长QQ:729038198 关于我们 投诉建议