回 帖 发 新 帖 刷新版面

主题:[讨论]Fortran数组越界问题

Fortran数组越界问题
求助:有个关于最小二乘法拟合趋势的程序,编译没错,但运行是会出现数组越界问题,请教下高手是哪里的问题?
!        CC*CHANG 为拟合趋势程序,其中AA-是输入变量,151为资料长度,
!        CC       3,4表示拟合的阶数,3=2阶+1,4等于拟合阶数+2
!        CC           X, AT, ST均为工作数组
!        CC       调过此子程序后,AA的时间序列已去除了趋势
        PROGRAM MAIN
        PARAMETER(N=73,m=36*13)      
      DIMENSION X(2,N),AT(3),ST(3,4),AA(N),vwnd(m,n),vwnd2(m,n)
        !读入数据
        open(1,file='vwnd.1996.txt',form='formatted')
        do it=1,n
         read(1,*)(vwnd(i,it),i=1,m)        
        enddo
        close(1)           
        write(*,*)'input done'
        
         !线性拟合
        do i=1,m
          do it=1,n
           aa(it)=vwnd(i,it)
!           write(*,*) aa(it),it        
          enddo
        CALL CHANG(AA,n,3,4,X,AT,ST)
           do j=1,n
              vwnd2(i,j)=aa(j)                   !拟合后赋值输出
           enddo
        enddo
        write(*,*) 'call subpoutine done'                                                         
        open(1,file='vwndaa.1996.txt',form='formatted')
        
        write(1,*)((vwnd2(i,it),i=1,m),it=1,n)
        
        close(1)
        write(*,*) 'output done'
!        open(2,file='vwndaa.grd',form='formatted')
!        do i=1,m
!        write(2,rec=1)(vwnd2(i,it),it=1,n)
!        enddo
!        close(2)
          end
        
      SUBROUTINE CHANG(UU,N,M1,M2,X,A,S)                               
      DIMENSION X(2,N),A(M1),S(M1,M2),UU(N)                            
      DO 100 NT=1,N                               
100   X(1,NT)=FLOAT(NT)                           
      DO 300 NT=1,N                               
300   X(2,NT)=UU(Nt)                                               
      M=M1-1                                                      
      DO 5 I=1,M1                                                  
      DO 5 J=1,M2                                                  
5    S(I,J)=0.0                                                   
      S(1,1)=N                                                     
      DO 6 J=1,N                                                   
6    S(1,M2)=S(1,M2)+X(2,J)                                        
      DO 7 I=2,M1                                                   
      I1=I-1                                                        
      MI=M1+I-2                                                     
      DO 7 J=1,N                                                    
      S(I,1)=S(I,1)+X(1,J)**I1                                      
      S(M1,I)=S(M1,I)+X(1,J)**MI                                    
7    S(I,M2)=S(I,M2)+X(1,J)**I1*X(2,J)                             
      DO 8 J=2,M                                                    
      DO 8 I=J,M                                                    
8    S(I,J)=S(I+1,J-1)                                             
      DO 9 I=1,M                                                    
      I1=I+1                                                        
      DO 9 J=I1,M1                                                  
9    S(I,J)=S(J,I)                                                 
      CALL GS(M1,S,A,M2,FAIL)                                       
      DO 500 NT=1,N                                                 
500   UU(NT)=UU(NT)-A(1)-A(2)*FLOAT(NT)-A(3)*FLOAT(NT*NT)           
200   CONTINUE                                                      
      RETURN                                                        
      END                                                           
      SUBROUTINE GS(N,A,X,N1,FAIL)                                  
      DIMENSION A(N,N1),X(N)                                        
      DO 11 K=1,N                                                   
      DO 6 I=K,N                                                    
      IF(A(I,K).EQ.0.0) GO TO 6                                     
      I0=I                                                          
      GO TO 7                                                       
6    CONTINUE                                                      
      FAIL=-1.0                                                     
      RETURN                                                        
7    IF(I0.EQ.K) GO TO 9                                           
      DO 8 J=K,N1                                                   
      T=A(K,J)                                                      
      A(K,J)=A(I0,J)                                                
8    A(I0,J)=T                                                     
9    K1=K+1                                                        
      DO 10 J=K1,N1                                                 
10    A(K,J)=A(K,J)/A(K,K)                                          
      DO 11 I=K1,N                                                  
      DO 11 J=K1,N1                                                 
11    A(I,J)=A(I,J)-A(I,K)*A(K,J)                                   
      X(N)=A(N,N1)                                                  
      NI=N-1                                                        
      DO 12 K=1,N1                                                  
      K1=NI-K+1                                                     
      X(K1)=A(K1,N1)                                                
      K2=K1+1                                                       
      DO 12 J=K2,N                                                  
12   X(K1)=X(K1)-A(K1,J)*X(J)                                      
      FAIL=0.0                                                      
      RETURN                                                        
      END

回复列表 (共19个回复)

11 楼

谢谢回复!
  我用上面拟合的程序生成的txt文件转成grd文件,用下面的程序转的:
program ex
implicit none
integer i,j,k
integer,parameter :: x=144
integer,parameter :: y=73
integer,parameter :: t=365
real :: dat(x,y,t)

open(100,file='e:\ncep\uwndaa.txt',form='formatted')
 do k=1,t

     read(100,*) ((dat(i,j,k),i=1,x),j=1,y)

 
 print *,k
end do
close(100)

open(200,file='e:\ncep\uwndaaa.grd',form='binary')
do k=1,t
  do j=1,y
    do i=1,x
   write(200) dat(i,j,k)
    end do
  end do
 end do
close(200)

stop
end
结果只能读到364天的数据
显示错误如下:
forrt1: severe <24>: end-of file during read,unit 100,file e:\ncep\uwndaa.txt
若将程序中365改为364,程序就能运行。请问下是哪的问题呢?

12 楼

也有可能是txt转grd这个程序的问题,麻烦高手帮我看下吧?谢谢

13 楼

这个。。。
你输出的时候我也没看出来哪儿输出了365啊?

PS:NCEP?你是华北院的?

14 楼

哦  不好意思啊,我没讲清楚,拟合的程序后来我把格点改为了x方向144格点,y方向73格点,全球范围,t时间改为365天,资料是NCEP\NCAR再分析风场资料,
PARAMETER(N=73,m=36*13)
改为了:  
PARAMETER(N=365,m=144*73)
生成的txt转grd文件时提示错误,只能到364
我刚又数了下总数,那个拟合的程序应该没问题,可能是txt转grd文件时有问题,可不知道怎么改呢?
PS:我是南气院的 气象学专业

15 楼

拟合程序中的:
PARAMETER(N=73,m=36*13) 
这行中,M值你设置对了么?

16 楼

恩 对的 PARAMETER(N=365,m=144*73) 这样子的 
拟合后生成的txt文件中数据是这样的:
-0.3873872     -0.3964562     -0.4055951     -0.4161802     -0.4243862    
 -0.4266428     -0.4338557     -0.4344189     -0.4374759     -0.4480302    
 -0.4421558     -0.4449532     -0.4458366     -0.4488226     -0.4454437    
 -0.4480694     -0.4415421     -0.4411916     -0.4341181     -0.4280878    
 -0.4195198     -0.4159690     -0.4072542     -0.4031726     -0.3907706    
。。。。。。。。。

那个转格式的程序读这个txt读不进去,我有pause查了下程序,read那句有错误。


17 楼


冒昧的问下,您介意QQ交流下吗?我把元数据发你麻烦你帮我看下呢?我的QQ号:412974579
您介意的话也没关系 呵呵

18 楼

real :: dat(x*y)
    open(100,file='e:\ncep\uwndaa.txt',form='formatted')
    open(200,file='e:\ncep\uwndaaa.grd',form='binary')

    do k=1,t
        read(100,*)(dat(i),i=1,x*y)
        write(200) (dat(i),i=1,x*y)
    end do

    close(200)
    close(100)

19 楼

问题解决了,非常感谢!

我来回复

您尚未登录,请登录后再回复。点此登录或注册