fortran进行sst的eof分析程序无法执行

lily222

楼主需要做印度洋的海温eof分析，拿到的是1891年至今的全球海表温度月平均数据（1*1），用grads抠出了印度洋地区的海温，然后打算用fortran程序进行eof分析，没有报错但是一运行就会显示程序停止，显示界面如下图。
楼主以为是数据点太多了，就把时间点由1000多个缩减到前600个，空间点6161个（-30S—30N，30E—130E），然后运行结果还是一样。
望各位大神帮忙看看是否程序有问题，还是别的原因。非常感谢！
gs程序和fortran程序如下：
'reinit'
'sdfopen d:\sst.mon.mean.nc'
'set fwrite d:\sst.grd'
'set gxout fwrite'
'set x 30 130'
'set y -30 30'
'set z 1'
'set t 1 600'
'd sst'
'disable fwrite'
;

program EOF
implicit none
integer m,n,kg,ks
parameter (m=6161,n=600,kg=2,ks=2)
!   m空间点数，n时间点数
!   kg所要看的前个kg特征向量、特征值. kg<=min(m,n)
!   ks读入数据的处理：0、不处理；1、中心化；2、标准化
interface
    subroutine  jcb(m,A,V)
       integer m
    real A(m,m),v(m,m)
    endsubroutine jcb
    subroutine juping(m,n,X)
       integer m,n
    real X(m,n)
    endsubroutine juping
    subroutine jpbzh(m,n,X)
       integer m,n
    real X(m,n)
    endsubroutine jpbzh
    subroutine cheng(A,m1,n1,B,m2,n2,C)
       integer m1,n1,m2,n2
    real A(m1,n1),B(m2,n2),C(m1,n2)
    endsubroutine cheng
endinterface
!   子程序jcb：将实对称矩阵A处理成对角矩阵A（也就是求出了A其所有特
!   征值），同时求出矩阵V（空间函数）
!   子程序juping：将矩阵X中心化
!   子程序jpbzh：将中心化了的矩阵X再标准化
!   子程序cheng：两矩阵相乘，A(m1,n1)*B(m2,n2)=C(m1,n2)   
real X(m,n),A(min(m,n),min(m,n)),V(m,min(m,n)),Z(min(m,n),n),L(min(m,n))
real Ztemp(n),temp,sum,Vt(min(m,n),m),U(n,n),Xt(n,m),v1(m,min(m,n))
integer i,j,k
!   X(m,n)气象资料；V空间函数；Z时间函数；L特征值

!   读入数据－－注意 ：X(m,n)空间m、时间n ，要与气象资料相对应－－－－－－－
    open(22,file='d:\sst.grd',form='binary',recl=m*n*4)
    read(22,rec=1) ((X(i,j),i=1,m),j=1,n)    !!!!! 20-160E,40-80N (5.0×5.0)
close(22)
!   －－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

!   读入数据的处理－－ks：0、不处理；1、中心化；2、标准化－－－－－－－－－
    if (ks>0) then
       call juping(m,n,X)
    if (ks>1)  call jpbzh(m,n,X)
endif
!   －－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

!   EOF主程序－－－－－－－－－－-－－－－－－－－－－－－－－－－－－－－
    if (m<=n) then
    Xt=transpose(X)
       call cheng(X,m,n,Xt,n,m,A)
    call jcb(m,A,V)  
    do i=1,m-1
    do j=i+1,m
       if (A(i,i)<A(j,j)) then
          temp=A(i,i)
       A(i,i)=A(j,j)
             A(j,j)=temp
       do k=1,m
                Ztemp(k)=V(k,i)
       V(k,i)=V(k,j)
       V(k,j)=Ztemp(k)
       enddo
       endif
    enddo
    enddo  
    Vt=transpose(V)
    call cheng(Vt,m,m,X,m,n,Z)
else
    Xt=transpose(X)
       call cheng(Xt,n,m,X,m,n,A)
       call jcb(n,A,U)
    do i=1,n-1
    do j=i+1,n
       if (A(i,i)<A(j,j)) then
          temp=A(i,i)
       A(i,i)=A(j,j)
             A(j,j)=temp
       do k=1,n
                Ztemp(k)=U(k,i)
       U(k,i)=U(k,j)
       U(k,j)=Ztemp(k)
       enddo
       endif
    enddo
    enddo
       call cheng(X,m,n,U,n,n,V)
    do i=1,m
    do j=1,n
          if (A(j,j)>((1e-6)*A(1,1)) ) then
          V(i,j)=V(i,j)/sqrt(A(j,j))
    else
       V(i,j)=0
    endif
    enddo
    enddo
       Vt=transpose(V)
    call cheng(Vt,n,m,X,m,n,Z)
endif
!   －－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

!   输出结果－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
    open(22,file='c:\users\xin\desktop\result.txt',form='formatted')
sum=0
do i=1,min(m,n)
       L(i)=A(i,i)
    sum=sum+L(i)
enddo
write(22,'(a2,i2,a10)') "前",kg,"个特征值:"
write(22,'(<kg>g10.3)') (L(i),i=1,kg)
do i=1,min(m,n)
    L(i)=L(i)/sum
enddo
    write(22,'(a2,i2,a23)') "前",kg,"个特征相量 每个的贡献％"
write(22,'(<kg>f10.3)') (L(i)*100,i=1,kg)
do i=1,min(m,n)
    if (i==1) then
       L(1)=L(1)
    else
       L(i)=L(i-1)+L(i)
    endif
enddo
    write(22,'(a4,i2,a22)') "至前",kg,"个特征相量的累积贡献％"
write(22,'(<kg>f10.3)') (L(i)*100,i=1,kg)
write(22,'(a2,i2,a28,i2,a2)') "前",kg,"个特征相量所对应的时间序列：",kg,"列"
do j=1,n
    write(22,'(<kg>f10.3)') (Z(i,j),i=1,kg)
    enddo
write(22,'(a2,i2,a12,i2,a2)') "前",kg,"个特征相量:",kg,"列"
do i=1,m
    write(22,'(<kg>f10.4)') (V(i,j),j=1,kg)
    enddo
close(22)
!   －－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
      
!   输出数据－"前kg个特征相量"－－GRADS画图使用－－－－－－－－－－－－－－－
    open(22,file='d:\eof3v.grd',form='binary')
    write(22) ((V(i,j),i=1,m),j=1,kg)   
close(22)
open(22,file='d:\time1.grd',form='binary')
write(22) Z(1,1:n)
close(22)
open(22,file='d:\time2.grd',form='binary')
write(22) Z(2,1:n)
close(22)
end program

subroutine  jcb(m,A,V)
    implicit none
    integer m
    real A(m,m),v(m,m)
    integer i,j,p,q,state
    integer*4 count
    real sx,cx,a1,a2,a3,max
    real sum,eps,guan,Atemp(m,m)
    write(*,*) "jcb"
       sum=0
    do i=1,m
    do j=1,m
       if (i==j) then
       V(i,j)=1.0
    else
       V(i,j)=0.
    endif
    if (i>j)  sum=sum+abs(A(i,j))
    enddo
    enddo
       eps=(1e-5)*sum/(m*m)
       write(*,*) "精度控制eps=",eps
    state=1
       do while(state==1)
          state=0
    max=eps
    do i=1,m
    do j=1,m
       if (j>i) then
       if (abs(A(i,j))>max) then
       state=1
       max=abs(A(i,j))
       p=i
       q=j
    endif   
    endif
    enddo
    enddo
       if (state==1) then
      count=count+1
    if (mod(count,500)==0)  write(*,*) count,max
             a1=A(p,q)
    a2=(A(p,p)-A(q,q))/2
    if (abs(a2)<1e-9) then
    a3=1.0
    else
    a3=sign(1.0,a2)*a1/sqrt(a1*a1+a2*a2)
    endif
    sx=a3/sqrt(2*(1+sqrt(1-a3*a3)))
    cx=sqrt(1-sx*sx)        
    do j=1,m
    if ((j/=p).and.(j/=q)) then
     Atemp(p,j)=A(p,j)*cx+A(q,j)*sx
     Atemp(q,j)=(-1.0)*A(p,j)*sx+A(q,j)*cx
    endif
             enddo
    do i=1,m
    if ((i/=p).and.(i/=q)) then
       Atemp(i,p)=A(i,p)*cx+A(i,q)*sx
       Atemp(i,q)=-1.0*A(i,p)*sx+A(i,q)*cx
    endif
             enddo
    do i=1,m
    do j=1,m
    if ((i/=p).and.(i/=q).and.(j/=p).and.(j/=q)) then
        Atemp(i,j)=a(i,j)
       endif
    enddo
    enddo
       Atemp(p,p)=A(p,p)*cx*cx+A(q,q)*sx*sx+2*A(p,q)*sx*cx
    Atemp(q,q)=A(p,p)*sx*sx+A(q,q)*cx*cx-2*A(p,q)*sx*cx
    Atemp(p,q)=(A(q,q)-A(p,p))*sx*cx+A(p,q)*(cx*cx-sx*sx)
    Atemp(q,p)=Atemp(p,q)
       do i=1,m
       sum=V(i,p)*cx+V(i,q)*sx
       V(i,q)=-1.0*V(i,p)*sx+V(i,q)*cx
       V(i,p)=sum
    enddo
    do i=1,m
    do j=1,m
    A(i,j)=Atemp(i,j)
    enddo
    enddo                       
          endif
    end do
endsubroutine jcb

subroutine juping(m,n,X)
    integer m,n
    real X(m,n)
    real temp(m)
    integer i,j
    do i=1,m
       temp(i)=0
       do j=1,n
       temp(i)=temp(i)+X(i,j)
          enddo
    temp(i)=temp(i)/n
    enddo
    do i=1,m
    do j=1,n
       X(i,j)=X(i,j)-temp(i)
    enddo
    enddo
endsubroutine juping

subroutine jpbzh(m,n,X)
    integer m,n
    real X(m,n)
    real temp(m)
    integer i,j
    do i=1,m
       temp(i)=0
    do j=1,n
       temp(i)=temp(i)+X(i,j)**2
    enddo
    temp(i)=sqrt(temp(i)/n)
    enddo
    do i=1,m
    do j=1,n
       X(i,j)=X(i,j)/temp(i)
    enddo
    enddo  
endsubroutine jpbzh

subroutine cheng(A,m1,n1,B,m2,n2,C)
    integer m1,n1,m2,n2
    real A(m1,n1),B(m2,n2),C(m1,n2)
    integer i,j,k
    if (n1/=m2) then
       write(*,*) "wrong -cheng "
    stop 1111
    endif
    do i=1,m1
    do j=1,n2
       C(i,j)=0
    do k=1,n1
       C(i,j)=C(i,j)+A(i,k)*b(k,j)
    enddo
    enddo
    enddo
endsubroutine cheng

霹雳电光鼠

这种问题以后还是找身边的人问，没人愿意看这么长的程序！

lily222

open(22,file='c:\users\xin\desktop\result.txt',form='formatted')
这个文件名楼主没改过来，应该open(22,file='d:\result.txt',form='formatted')，文件全部保存在d盘的，因为楼主去借别人电脑运行程序发帖的时候没全部改过来...试别的电脑还是一样的结果，望大神们指点！

lily222

霹雳电光鼠 发表于 2017-3-23 14:36
这种问题以后还是找身边的人问，没人愿意看这么长的程序！

问了老师，老师说是数据点太多，但是改了结果没变化...

小小书童

堆栈，应改是数组太大的缘故。

lily222

小小书童 发表于 2017-3-23 17:33
堆栈，应改是数组太大的缘故。

是m和n太大了吗?但是印度洋格点数确实已经缩到最小了，而且我尝试缩小时间点从1200缩到600，还是停止运行。

lily222

补充一下，楼主在自己电脑上点击构建按钮后会出现一个warning如下：
Debug/EOF.exe : warning LNK4084: total image size 305029120 exceeds max (268435456); image may not run

EOF.exe - 0 error(s), 1 warning(s)
然后不管楼主怎么缩小空间点和时间点，image size 305029120 这个数字都不变，这是不是说明不是因为空间点和时间点太多而是别的东西超出容量了呢？
借别人内存大一点的电脑就没有warning了，但是结果是一样的停止运行。

lily222

楼主发现gs文件有错误，x、y应该用lon、lat代替才对。然而修改后eof程序运行结果还是溢出...附上修改后的gs文件和eof程序。

janicewang

楼主！请问问题解决了吗？我遇到同样问题了