Fortran编程求助

yh1223

我用的是1973年-2013年，每年6月至8月的500hPa温度场的资料，已经用GrADS提取成了grd格式的。用Fortran计算APO指时，编译和连接时都没有错误，可以运行完后数据好像不太对，我的程序如下：PROGRAM IAPO
implicit none
dimension tt(144,73,41),tp(73,41),tr(144,73,41),apo(41)
      dimension ap(41),bz(41),tr1(41),tr2(41)
real fc,a,ave
      open(10,file='d:\bylw\data\air500.grd',form='binary',status='old')

      !读数据
do it=1,41
do j=1,73
do i=1,144
read(10)tt(i,j,it)
enddo;enddo;enddo
      write(*,*)'read data ok'

do j=1,73
      do it=1,41
      tp(j,it)=0.0
enddo
enddo
      do it=1,41
      do j=1,73
      do i=1,144
      tp(j,it)=tp(j,it)+tt(i,j,it)/144
      enddo;enddo;enddo
print*,tp

      do it=1,41
do j=1,73
do i=1,144
tr(i,j,it)=tt(i,j,it)-tp(j,it)
enddo;enddo;enddo

       write(*,*)'select region ok111'

do it=1,41
tr1(it)=0
enddo
      do it=1,41
do j=43,57
do i=25,49
        tr1(it)=tr1(it)+tr(i,j,it)/375
      enddo
enddo
enddo
do it=1,41
tr2(it)=0
enddo
      do it=1,41
do j=43,57
do i=73,97
        tr2(it)=tr2(it)+tr(i,j,it)/375
      enddo
enddo
enddo

      do it=1,41
   apo(it)=tr1(it)-tr2(it)
enddo
print*,apo
      write(*,*)'select region ok222'
!标准化处理
ave=0.0
do it=1,41
   ave=ave+apo(it)/41
enddo
do it=1,41
   ap(it)=apo(it)-ave
enddo
print*,ap
a=0.0
do it=1,41
  a=a+ap(it)**2
enddo
fc=sqrt(a/41)
print*,fc
      write(*,*)'select region ok333'
do it=1,41
bz(it)=ap(it)/fc
enddo
print*,bz
      write(*,*)'select region ok444'
      open(20,file='d:\bylw\data\apo.grd',form='binary')
write(20) (apo(it),it=1,41)
open(30,file='d:\bylw\data\bz.grd',form='binary')
write(30) (bz(it),it=1,41)
close(20)
close(30)
end

这里是用亚洲和太平洋（我选的区域是15-50N,180-120W,60E-120E）区域平均的扰动温度（T'）之差作为APO（亚洲-太平洋涛动）指数，T'=T-T的纬向平均。
是我程序有问题还是说本身提取的grd数据就有问题呢？麻烦各位大神指点下。

蓝莲花

数据可能有缺测的，你试试先把它挑出来

lqouc

要人家帮你看的话就尽量把代码的注释写清楚。你的每一步都是干什么，省的人家猜。
看起来程序没什么大问题，尽量检查一下原始数据和导出数据。
至于你做的这些其实再grads里面更好实现。完全没必要用fortran。
ps：不知道你有没有看到，科学网上已经有人对apo进行打假了，理由是在每周也有一块异常，但是这个异常被忽略了，这是极其不严谨的，这个apo也就不成立了吧。所以楼主做之前再仔细思考下哈。

yh1223

lqouc 发表于 2014-4-22 16:02
要人家帮你看的话就尽量把代码的注释写清楚。你的每一步都是干什么，省的人家猜。
看起来程序没什么大问题 ...

额。。。我的毕业论文定的是这个。。。继续写这个会不会对毕业答辩有影响啊？

lqouc

yh1223 发表于 2014-4-22 16:22
额。。。我的毕业论文定的是这个。。。继续写这个会不会对毕业答辩有影响啊？

本科毕业论文哈，那就没关系了，你自己研究的话还是仔细点好。

yh1223

lqouc 发表于 2014-4-22 16:35
本科毕业论文哈，那就没关系了，你自己研究的话还是仔细点好。

哦哦，那就好，我就是写个论文，吓死我了，还以为要重新准备呢。。。谢谢！

yh1223

lqouc 发表于 2014-4-22 16:35
本科毕业论文哈，那就没关系了，你自己研究的话还是仔细点好。

学长，我想问下我要是求500-200hPa平均的扰动温度，是不是要把500-200hPa的温度资料放到一个文件里？还是说每一层分开提取好点？

river

先看看你提取的资料是不是对的吧，你验证了你提取的二进制文件是对的吗？如果没有验证，那就是一个大漏洞····还有就是这个apo指数的定义还算比较简单，在grads的脚本里面可以直接计算

lqouc

yh1223 发表于 2014-4-22 16:42
学长，我想问下我要是求500-200hPa平均的扰动温度，是不是要把500-200hPa的温度资料放到一个文件里？还是 ...

直接求平均不就好了？这些grads都可以算，如果你就喜欢fortran的话也可以自己编程。

yh1223

river 发表于 2014-4-22 16:47
先看看你提取的资料是不是对的吧，你验证了你提取的二进制文件是对的吗？如果没有验证，那就是一个大漏洞· ...

我用GrADS提取你nc数据的gs文件是
'sdfopen d:\bylw\ncdata\air.mon.mean.nc'
'set gxout fwrite'
'set fwrite d:\bylw\data\air200.grd'
'set x 1 144'
'set y 1 73'
'set z 10'
t=306
while(t<=792)
tt=0
while(tt<=2)
'set t 't +tt
'd air'
tt=tt+1
endwhile
t=t+12
endwhile
'disable fwrite'
'reinit'
；
自己写的ctl文件：
dset d:\bylw\data\air200.grd
title NMC reanalysis atlas
undef -9.99E+33
xdef 144 linear     0.0      2.5
ydef 73  linear   -90.0      2.5
zdef 1  levels     200
tdef 123 linear 00Z01JUN1973     1mo
vars 1
air=>air 1 t,z,y,x Monthly Mean of Air temperature
endvars
这个有什么问题吗？