用fortran计算副高指数的时候越界是哪里出错了？

小亭子ng

     我用的是2013年ncep-1 的高度场资料（hgt.2013.nc）我用grads将其转换为grd.
    gs如下：
'reinit'
'sdfopen H:\data\fugao\ncep-1\hgt.2013.nc'
'set fwrite H:\data\fugao\ncep-1\hgt.grd'
'set gxout fwrite'
'set x 1 144'
'set y 1 73'
'set z 6'
t0=1
while (t0<=365)
'set t ' t0
'd hgt*0.1'
t0=t0+1
endwhile      
'disable fwrite'
;

然后我使用fortran计算副高指数时总是报错。
fortran程序如下：
!1 用NCEP  2.5分辨率逐日资料西太副高的5个特征指数，一天一个指数
!2 GM 面积指数，GQ 面积指数，GJ 西脊点指数，GX 脊线指数，GB 北界点指数
program fg_zhishu
implicit none
        integer,Parameter::x=144,y=73,mo=365
        integer,Parameter::lon1=37,lon2=73,lat10=41,lat80=69 ! lat10--北纬10度，lat80--北纬80度，LON1--东经90，LON2--东经180
        integer i,j,m,tem,tt,n,s,kr,kkr
        real h(x,y,mo),fgm(x,y,mo),fg(x,y,mo),n588(x,mo),s588(x,mo)
        real gm1(mo),gq2(mo),gj3(mo),gx4(mo),gb5(mo)

        real a,b,c,d,e,f,nn,ch,du
!读取多个文件
        character(len=1) :: index

        open(10,file='H:\data\fugao\ncep-1\grd\hgt.grd',form='binary')
             do m=1,mo
                read(10) ((h(i,j,m),i=1,x),j=1,y)
                enddo
        close(10)

!!!!!!!!!!!!!!!!!!!!!!!  判断副高的范围 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                do m=1,mo
                !588线南界
                        do i=lon1,lon2
                                s588(i,m)=0.0
                                do j=lat10,lat80
                                        if((h(i,j,m)/10.0>=588.0) .and. (h(i,j+1,m)/10.0>=588.0))then
                                                s588(i,m)=j               
                                                exit
                                        end if
                                end do
                        end do
                !588线北界
                        do i=lon1,lon2
                                n588(i,m)=0.0
                                do j=lat80,lat10,-1
                                        if((h(i,j,m)/10.0>=588.0) .and. (h(i,j-1,m)/10.0>=588.0))then
                                                n588(i,m)=j
                                                exit
                                        end if
                                end do
                        end do
                        do i=lon1,lon2
                !副高分布       
                                if(s588(i,m)/=0)then
                                        do j=s588(i,m),lat80
                                                fgm(i,j,m)=0.0;fg(i,j,m)=0.0
                                                if(h(i,j,m)/10.0>=588.0)then                                       
                                                        fgm(i,j,m)=1
                                                        fg(i,j,m)=h(i,j,m)
                                                end if                               
                                        end do
                                end if
                        end do
                end do


!!!!!!!!!!!!!!!!!!!!!!!! 1面积指数 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                do m=1,mo
                        gm1(m)=0.0
                        do i=lon1,lon2
                                if(s588(i,m)/=0)then
                                        do j=s588(i,m),n588(i,m)
                                                if((h(i,j,m)/10.0>=588.0))then
                                                        gm1(m)=gm1(m)+1
                                                end if
                                        end do
                                end if
                        end do
                end do
        print*, 'ok gm1'       

!!!!!!!!!!!!!!!!!!!!!!!! 2强度指数 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                do m=1,mo
                        gq2(m)=0.0
                        do i=lon1,lon2
                                if(s588(i,m)/=0)then
                                        do j=s588(i,m),n588(i,m)
                                                if((h(i,j,m)/10.0>=588.0))then
                                                        gq2(m)=gq2(m)+(h(i,j,m)/10.0-587.0)
                                                end if
                                        end do
                                end if
                        end do
                end do
        print*, 'ok gq2'
!!!!!!!!!!!!!!!!!!!!!!!!3 西脊点指数!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                do m=1,mo     
                        gj3(m)=0.0;tt=0
                        do i=lon1,lon2
                                do j=lat10,lat80
                                        if(h(i,j,m)/10.0>=588 .and. h(i+1,j,m)/10.0>=588)then
                                                tt=i
                                                exit                       
                                        end if
                                end do
                                if(tt/=0)then
                                        exit
                                end if       
                        end do
                        gj3(m)=(tt-1)*2.5
                        if(gj3(m)<=0)then
                                gj3(m)=180
                        end if
                end do
        print*, 'ok gj3'
!!!!!!!!!!!!!!!!!!!!!!!! 4 脊线位置指数 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                do m=1,mo
                        gx4(m)=0.0;tt=0
                        do i=lon1,lon2
                                if(n588(i,m)/=0 .and. s588(i,m)/=0)then
                                        tt=tt+1
                                        gx4(m)=gx4(m)+(n588(i,m)+s588(i,m))/2.0
                                end if
                        end do
                        if(gx4(m)/=0)then
                                gx4(m)=(gx4(m)-1)/float(tt)*2.5-90.0
                        end if
                end do       
        print*, 'ok gx4'
!!!!!!!!!!!!!!!!!!!!!!!! 5 北界位置指数 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       
                do m=1,mo
                        gb5(m)=0.0;tt=0
                        do i=lon1,lon2
                                if(n588(i,m)/=0.0)then
                                        tt=tt+1
                                        gb5(m)=gb5(m)+n588(i,m)                       
                                end if               
                        end do
                        if(gb5(m)/=0)then
                                gb5(m)=(gb5(m)-1)/float(tt)*2.5-90.0
                        end if
                end do

        print*, 'ok gb5'

!!!!!!!!!!!!!!!!!!!!!!!! 输出结果 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!循环输出文件    可根据需要输出数据
!                write( year,'(i4)' ) t+1977                !change 1977为起始年份-1
!                write(*,*) t+1977,year                     !!change
                open (9,file='2013.1.grd',form='binary',status='replace')
                        do m=1,mo
                                write(9) gm1(m)
                        end do
                close(9)
       
           open (9,file='2013.2.grd',form='binary',status='replace')
                        do m=1,mo               
                                write(9) gq2(m)
                        end do
                close(9)
                open (9,file='2013.3.grd',form='binary',status='replace')
                        do m=1,mo
                                write(9) gj3(m)
                        end do
                close(9)
                open (9,file='2013.4.grd',form='binary',status='replace')
                        do m=1,mo
                                write(9) gx4(m)
                        end do
                close(9)
                open (9,file='2013.5.grd',form='binary',status='replace')
                        do m=1,mo
                                write(9) gb5(m)
                        end do
                close(9)



end program

hhh走走看看

楼主，脊线位置和北界为9条经线交点的纬度的平均值，你那个程序计算了不止9条经线吧，这样也可以吗？

小亭子ng

楼主认怂了，已找到错误

nunu18

小亭子ng 发表于 2014-7-31 19:53
楼主认怂了，已找到错误

楼主知道原因了，不如和大家分享，这样好多人都可以受益

小亭子ng

nunu18 发表于 2014-7-31 20:28
楼主知道原因了，不如和大家分享，这样好多人都可以受益

这些个应该都没问题，只是我自己路径那边犯了个低级错误。。。

zhanghaiyan

对于菜鸟来说，这些程序很有借鉴意义，赞一个

hhh走走看看

{:eb502:}{:eb502:}

hhh走走看看

计算西太副高指数的范围定在北纬10度，北纬80度，东经90，东经180可以吗

katrina

我想问一下您用这个程序计算出的面积指数是否正确？

altao

小亭子ng 发表于 2014-7-31 21:15
这些个应该都没问题，只是我自己路径那边犯了个低级错误。。。

是路径哪里错了，我也是一样的问题，fortran不熟