请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 3044|回复: 1

[求助] 关于fortron的eof程序的一个问题

[复制链接]

新浪微博达人勋

发表于 2018-4-12 19:18:49 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
我用的是李建平老师的程序改了好多遍结果站点数换了一个数和文件就出不来结果了。。。 EOF.f90 (11.06 KB, 下载次数: 3)
LM_G`I1NQ$IMIDEKFFXHQBN.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-4-12 19:23:21 | 显示全部楼层
parameter(m=2467,nyear=84,n=nyear*12,mnl=n,ks=-1)
integer:: sta,year,k,i,j,q
real tm(m,n),t1(m,n+60),egvt(m,mnl),ecof(mnl,n),er(mnl,4),lat,lon,tmax,tmin,pr,afw,rha,w10m

open(20,File='E:\www\REGCM45.txt')
read(20,*)
do k = 1, m
do j = 1,n+60
        read(20,*) sta,lat,lon,year,q,t1(k,j),tmin,tmax,pr,afw,rha,w10m
        IF(J>60)THEN
           Tm(k,j-60)=T1(k,j)
        ENDIF
end do
end do
CLOSE(20)
print *,tm(m,n)
pause 111

call eof(m,n,mnl,tm,ks,er,egvt,ecof)

        open(2,file="E:\www\eof\er.txt")
        do i=1,mnl
        write(2,"(4f20.8)") (er(i,j),j=1,4)
        enddo
        close(2)
         

        open(3,file="E:\www\eof\ecof.txt")
        do j=1,n
        write(3,"(4f20.8)") (ecof(i,j),i=1,mnl)
        enddo
        close(3)
       

        open(4,file="E:\www\eof\egvt.txt")
        do i=1,m
        write(4,"(4f20.8)") (egvt(i,j),j=1,mnl)
        enddo
        close(4)

end

        !---------------------------------------------------------------
        subroutine eof(m,n,mnl,f,ks,er,egvt,ecof)

                implicit none
                integer*4                                ::        m,n,mnl,ks
                real*4                                        ::        f(m,n),er(mnl,4),egvt(m,mnl),ecof(mnl,n)
                real*4,allocatable                ::        cov(:,:),s(:,:),d(:) !---work array
        !--------------------preprocessing data

                call transf(m,n,f,ks)
                allocate(cov(mnl,mnl))
        !--------------------crossed product matrix of the data f(m,n)       

                call crossproduct(m,n,mnl,f,cov)
               
                allocate(s(mnl,mnl))
                allocate(d(mnl))
        !--------------------eigenvalues and eigenvectors by the jacobi method

                call jacobi(mnl,cov,s,d,0.00001)
               
        !--------------------specified eigenvalues
                call arrang(mnl,d,s,er)
        !--------------------normalized eignvectors and their time coefficients
                call tcoeff(m,n,mnl,f,s,er,egvt,ecof)
                return
        end
        !-----------------------------------------------------------
       
        !--------------------------------------------------------------
        !---preprocessing data to provide a field by ks.
        !---input: m,n,f
        !------m: number of grid-points
        !------n: lenth of time series
        !------f(m,n): the raw spatial-temporal seires
        !------ks: contral parameter:ks=-1: self; ks=0: depature; ks=1: normalized depature
        !---output: f
        !------f(m,n): output field based on the control parameter ks.
        !------work variables: fw(n)
        !----------------------------------------------------------------
        subroutine transf(m,n,f,ks)
                implicit none
                integer*4                                        ::        m,n,ks
                real*4                                                ::        f(m,n)
                real*4,allocatable                        ::        fw(:),wn(:)     !---work array
                integer*4                                        ::        i0,i,j
                real*4                                                ::        af,sf,vf
               
                allocate(fw(n))
                allocate(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 has invalid data.'
                        write(*,*)'There are totally ',i0,'gridpionts with invalid data.'
                        write(*,*)'The array wn(m) stores the positions of those invalid grid-points.'
                        write(*,*)'You must pick off those invalid data from the orignal field---$ '  
                        write(*,*)'$---and then reinput a new field to calculate its eofs.'
                        write(*,*)'************************  fault  ************************************'
                endif            
                if(ks.eq.-1)return
                if(ks.eq.0)then                !anomaly of f
                        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
                if(ks.eq.1)then                 !normalizing f
                        do i=1,m
                                do j=1,n
                                        fw(j)=f(i,j)
                                enddo
                                call meanvar(n,fw,af,sf,vf)
                                if(sf==0.0)then
                                        do j=1,n
                                                f(i,j)=0.0
                                        enddo
                                else
                                        do j=1,n
                                                f(i,j)=(f(i,j)-af)/sf
                                        enddo
                                endif
                        enddo
                endif
                return
        end
        !--------------------------------------------------------
       
        !----------------------------------------------------------
        !---crossed product martix of the data. it is n times of covariance matrix of the data if ks=0 (i.e. for anomaly).
        !---input: m,n,mnl,f
        !---m: number of grid-points
        !---n: lenth of time series
        !---mnl=min(m,n)
        !---f(m,n): the raw spatial-temporal seires
        !---output: cov(mnl,mnl)  
        !---cov(m,n)=f*f' or f'f dependes on m and n.
        !---it is a mnl*mnl real symmetric matrix.
        subroutine crossproduct(m,n,mnl,f,cov)
                implicit none
                integer*4                ::        m,n,mnl
                real*4                        ::        f(m,n),cov(mnl,mnl)
                integer*4                ::        i,j,is,js
               
                        if((n-m)<0)then
                                do i=1,mnl
                                        do 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)
                                        enddo
                                enddo
                        else
                                do i=1,mnl
                                        do 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)
                                        enddo
                                enddo
                        endif
                       
                return       
        end
        !----------------------------------------------------------
       
       
        !----------------------------------------------------------
        !---computing eigenvalues and eigenvectors of a real symmetric matrix
        !---a(m,m) by the jacobi method.
        !---input: m,a,s,d,eps
        !------m: order of matrix
        !------a(m,m): the covariance matrix
        !------eps: given precision
        !---output: s,d
        !------s(m,m): eigenvectors
        !------d(m): eigenvalues
        subroutine jacobi(m,a,s,d,eps)
                implicit none
                integer*4                ::        m
                real*4                        ::        a(m,m),s(m,m),d(m),eps
                integer*4                ::        i,j,i1,l,iq,iq1,ip
                real*4                        ::        g,s1,s2,s3,v1,v2,v3,u,st,ct
                       
                        s=0.0
                        do i=1,m
                                s(i,i)=1.0
                        enddo
                        g=0.
                        do i=2,m
                                i1=i-1
                                do j=1,i1
                                        g=g+2.0*a(i,j)*a(i,j)
                                enddo
                        enddo
                        s1=sqrt(g)
                        s2=eps/float(m)*s1
                        s3=s1
                        l=0
50                        s3=s3/float(m)

60                        do iq=2,m

                                iq1=iq-1                                       
                                do ip=1,iq1
                                        if(abs(a(ip,iq)).lt.s3) exit
                                        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,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
                                                s(i,ip)=g
                                        enddo

                                        do i=1,m
                                                a(ip,i)=a(i,ip)
                                                a(iq,i)=a(i,iq)
                                        enddo

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

                        enddo

                        if((l-1)==0)then

                                l=0
                                go to 60
                        else
                                go to 150
                        endif
150                        if(s3.gt.s2) then
                         goto 50
                        end if
                        do i=1,m
                                d(i)=a(i,i)
                        enddo

                return
        end
        !-----------------------------------------------------
       
        !-----------------------------------------------------
        !---provides a series of eigenvalues from maximuim to minimuim.
        !---input: mnl,d,s
        !------d(mnl): eigenvalues
        !------s(mnl,mnl): eigenvectors
        !---output: er
        !------er(mnl,1): lamda (eigenvalues), its equence is from big to small.
        !------er(mnl,2): accumulated eigenvalues from big to small
        !------er(mnl,3): explained variances (lamda/total explain) from big to small
        !------er(mnl,4): accumulated explaned variances from big to small
        subroutine arrang(mnl,d,s,er)
                implicit none
                integer*4                ::        mnl
                real*4                        ::        d(mnl),s(mnl,mnl),er(mnl,4)
                integer*4                ::        i,mnl1,k1,k2
                real*4                        ::        c,tr
               
                tr=0.0
                do i=1,mnl
                        tr=tr+d(i)
                        er(i,1)=d(i)
                enddo

                mnl1=mnl-1
                do k1=mnl1,1,-1
                        do 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 i=1,mnl
                                                c=s(i,k2+1)
                                                s(i,k2+1)=s(i,k2)
                                                s(i,k2)=c
                                        enddo
                                endif
                        enddo
                enddo

                er(1,2)=er(1,1)
               
                do i=2,mnl
                        er(i,2)=er(i-1,2)+er(i,1)
                enddo

                do i=1,mnl
                        er(i,3)=er(i,1)/tr
                        er(i,4)=er(i,2)/tr
                enddo

                return
        end
        !--------------------------------------------
       
        !--------------------------------------------
        !---provides standard eigenvectors and their time coefficients
        !---input: m,n,mnl,f,s,er
        !------m: number of grid-points
        !------n: lenth of time series
        !------mnl=min(m,n)
        !------f(m,n): the raw spatial-temporal seires
        !------s(mnl,mnl): eigenvectors
        !------er(mnl,1): lamda (eigenvalues), its equence is from big to small.
        !------er(mnl,2): accumulated eigenvalues from big to small
        !------er(mnl,3): explained variances (lamda/total explain) from big to small
        !------er(mnl,4): accumulated explaned variances from big to small
        !---output: egvt,ecof
        !------egvt(m,mnl): normalized eigenvectors
        !------ecof(mnl,n): time coefficients of egvt
        subroutine tcoeff(m,n,mnl,f,s,er,egvt,ecof)
                implicit none
                integer*4                                        ::        m,n,mnl
                real*4                                                ::        f(m,n),s(mnl,mnl),er(mnl,4),egvt(m,mnl),ecof(mnl,n)
                real*4,allocatable                        ::        v(:)  !---work array
                integer*4                                        ::        i,j,js,is
                real*4                                                ::        c
               
                allocate(v(mnl))
                do j=1,mnl
                        do i=1,m
                                egvt(i,j)=0.
                        enddo
                        do i=1,n
                                ecof(j,i)=0.
                        enddo
                enddo
        !---normalizing the input eignvectors s
                do 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
                enddo
        !-----------------------------------------------
                if(m.le.n) then
                        do js=1,mnl
                                do i=1,m
                                        egvt(i,js)=s(i,js)
                                enddo
                        enddo
                        do 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)*s(i,is)
                                        enddo
                                        enddo
                        enddo
                else
                        do i=1,m
                                do j=1,n
                                        v(j)=f(i,j)
                                enddo
                                do js=1,mnl
                                        do j=1,n
                                                egvt(i,js)=egvt(i,js)+v(j)*s(j,js)
                                        enddo
                                enddo
                        enddo
                        do js=1,mnl
                                do j=1,n
                                        ecof(js,j)=s(j,js)*sqrt(abs(er(js,1)))
                                enddo
                                do i=1,m
                                        egvt(i,js)=egvt(i,js)/sqrt(abs(er(js,1)))
                                enddo
                        enddo
                endif
               
                return       
        end
        !---------------------------------------------

        !---------------------------------------------
        !---computing the mean ax, standard deviation sx and variance vx of a series x(i) (i=1,...,n).
        !---input: n and x(n)
        !------n: number of raw series
        !------x(n): raw series
        !---output: ax, sx and vx
        !------ax: the mean value of x(n)
        !------sx: the standard deviation of x(n)
        !------vx: the variance of x(n)
        !------by dr. li jianping, may 6, 1998.
        subroutine meanvar(n,x,ax,sx,vx)
                implicit none
                integer*4                ::        n
                real*4                        ::        x(n)
                real*4                        ::        ax,vx,sx
                integer*4                ::        i
               
                ax=0.
                vx=0.
                sx=0.
                do i=1,n
                        ax=ax+x(i)
                enddo                       
                ax=ax/float(n)
                do i=1,n
                        vx=vx+(x(i)-ax)**2
                enddo
                vx=vx/float(n)
                sx=sqrt(vx)
               
                return
        end
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

快速回复 返回顶部 返回列表