求hadley环流边界指数的程序，程序中有部分是错的，但是不知道具体怎么改，求讨论

1021582751

!求各年（1951-2010）全球纬向平均MSF值
!定义Hadley环流强度指数为：北半支最大MSF值为北半球强度指数（NHCI），南半支最小MSF值为南半球强度指数（SHCI）
!本程序一并求出NHCI和SHCI将60年的值分别写入txt文件，以备后期处理使用
!得到的所有值都是标准单位
!Copyright by Aires

program msf_year
implicit none

!声明程序所需变量并赋初值
real,parameter:: pi=3.14, a=6371000.0, g=9.8
real:: v(46,19,99), psi_up(46,19,99),psi_down(46,19,99),psi(46,19,99),w1(19),NHCE(99),SHCE(99)
real:: p(19)=(/101300,100000,92500,85000,70000,60000,&
               50000,40000,30000,25000,20000,15000,&
                           10000,7000,5000,3000,2000,1000,0/)
integer:: i,j,k,t,irec
data ((v(i,1,t),i=1,46),t=1,99) /4554*0/     !把第1层v设为0
data ((v(i,19,t),i=1,46),t=1,99) /4554*0/    !把第19层v设为0
data ((psi_up(i,1,t),i=1,46),t=1,99) /4554*0/
data ((psi_down(i,19,t),i=1,46),t=1,99) /4554*0/
real:: p_up(18)=(/-1300,-7500,-7500,-15000,-10000,-10000,-10000,-10000,-5000,&  !向上积分气压差值
                 -5000,-5000,-5000,-3000,-2000,-2000,-1000,-1000,-1000/)
real:: p_down(18)=(/1000,1000,1000,2000,2000,3000,5000,5000,5000,5000,10000,10000,&   !向下积分气压差值
                   10000,10000,15000,7500,7500,1300/)
real:: phi(73)=(/((i-37)*2.5*pi/180.0,i=1,73)/)      这个语句不太确定是什么意思，但是它是基于将纬度隔2.5共分为73个值，即90、87.5、85···············-87.5、-90，但是我的数据是隔4，即90、86、82··········-86、-90这样的，就是说我的纬度维数为46，不是73，不知道这句话该怎么设，导致后面写NHCE时就会出错，希望大家帮忙看看！
real ::  e=0.1e10
!real:: ph(46)=(/-90,-86,-82,-78,-74,-70,-66,-62,-58,-54,-50,-46,-42,-38,&
!          -34,-30,-26,-22,-18,-14,-10,-6,-2,2,6,10,14,18,22,26,30,&
!           34,38,42,46,50,54,58,62,66,70,74,78,82,86,90/)
!把二进制文件中的v读入数组中

open(10,file='d:\lunwen\2001-2099-01.grd',form='binary',status='old')
!irec=1
do t=1,99
  do j=2,18
    read(10) (v(i,j,t),i=1,46)
!        irec=irec+1
  end do
end do
close(10)

!求向上MSF

do t=1,99
  do j=2,19
    do i=1,46
      psi_up(i,j,t)=psi_up(i,j-1,t)+(2*pi*a*cos(phi(i))*(v(i,j,t)+v(i,j-1,t))*p_up(j-1))/(2*g)
    end do
  end do
end do

!同理，求向下MSF

do t=1,99
  k=1
  do j=19,2,-1
    do i=1,46
      psi_down(i,j-1,t)=psi_down(i,j,t)+(2*pi*a*cos(phi(i))*(v(i,j-1,t)+v(i,j,t))*p_down(k))/(2*g)
    end do
    k=k+1
  end do
end do

!求总MSF

do j=1,19
  w1(j)=p(j)/101300.0
end do

do t=1,99
  do j=1,19
    do i=1,46
      psi(i,j,t)=w1(j)*psi_up(i,j,t)+(1-w1(j))*psi_down(i,j,t)
    end do
  end do
end do

!加入边界条件

do t=1,99
  psi(:,1,t)=0
  psi(:,19,t)=0
  psi(1,:,t)=0
  psi(46,:,t)=0
end do
!open(200,file='d:\lunwen\msfX.grd',form='unformatted')
!write(200) ((psi(i,j,1),i=1,46),j=1,19)

!close(200)

!找出各年NHCE和SHCE分别写入两个txt文件
!NHCE、SHCE定义为向极地方向500hPa高度上msf为0所在纬度

open(20,file='d:\lunwen\NHCE.txt')
open(30,file='d:\lunwen\SHCE.txt')

do t=1,99
  NHCE(t)=0.0
  SHCE(t)=0.0
   do  i=2,72
if(abs(psi(i,7,t)-0.0)<=e)then
     if(phi(i)>NHCE(t))then
          NHCE(t)=phi(i)
         elseif(phi(i)<SHCE(t))then
          SHCE(t)=phi(i)
     end if
end if
   enddo
  write(20,*) NHCE(t)
  write(30,*) SHCE(t)
end do

close(20)
close(30)

end






  

lqouc

你要把想讨论的地方标红，或者突出出来，找你的问题都要找半天肯定没人搭理你。
还有把你尝试的结果也说一下，报错是什么？贴出来。
不管怎么说，先去看提问的智慧。

猫的他

可以参考，王盘兴1994年算大气经圈流函数的文章，这个算法应该是他提出来的