用fortran计算南亚高压强度画图有问题

lucy_9208

资料用的是欧洲中心的资料，用matlab转成dat文件之后，用grads进行夏季（6-8月）平均，得到1958-1979年的夏季平均位势高度资料，范围是LAT 20°S-90°N，LON 0°-360°，高度五层。要计算南亚高压的强度指数，取100hpa上位势高度大于16600gpm的格点，其位势高度与16600之差的和为南亚高压的强度，范围选取在30°E-180°W，10°S-90°S。
1.这是计算夏季平均的gs
'reinit''open e:/graduate/hgt.1958-2002.ctl''set gxout fwrite''set fwrite e:/graduate/hgt.summer.dat''set x 1 144''set y 1 45'it=6while(it<=684)iz=1 while(iz<=5) 'set z 'iz'' 'd ave(z/9.8,t='it',t='it+2')' iz=iz+1 endwhileit=it+12endwhile'disable fwrite'*end

2.这是计算南亚高压强度的fortran程序
program SAHIntensity
implicit none
integer:: nx,ny,nz,nt,ix,iy,iz,it,r
parameter (nx=144,ny=45,nz=5,nt=57)
real:: z(nx,ny,nz,nt)
real:: i(nt)

open(11,file="e:\graduate\hgt.summer.dat",form="unformatted",recl=nx*ny*nz*nt*4)
read(11,rec=1) ((((z(ix,iy,iz,it),ix=1,nx),iy=1,ny),iz=1,nz),it=1,nt )
close(11)

do it=1,nt
iz=5
i(it)=0
  do ix=61,144
    do iy=13,29
          r=z(ix,iy,iz,it)-16600.0
          if(r .GE. 0.3E-10) then
            i(it)=i(it)+r
          end if
        end do
  end do
end do

open(12,file="e:\graduate\hgt.intensity.dat",form="unformatted",recl=nt*4)
write(12) (i(it),it=1,nt)
close(12)
end
不知道如何将图粘贴上来，就放在附件里了。
画出来的图在t=1的时候接近于零，t=2之后的时刻计算值就正常了。想请问各位是程序哪里出了问题呢？看了好久也没看出来。先谢谢各位了！

wohenkuang

的图在t=1的时候接近于零，t=2之后的时刻计算值就正常了

wohenkuang

t=2的时候接近于零，t=1之后的时刻计算值