关于fortron的eof程序的一个问题

成信小学生啊

我用的是李建平老师的程序改了好多遍结果站点数换了一个数和文件就出不来结果了。。。 EOF.f90 (11.06 KB, 下载次数: 3)

2018-4-12 19:17 上传

点击文件名下载附件

成信小学生啊

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