求助！计算施能老师夏季风指数，结果与其相比差距较大，已贴gs & Fortran程序

timeover

在计算施能老师夏季风指数结果后与其文章里的结果对比差异较大，差异感觉并非单纯由资料的不同引起，主要在于指数的位相不同。（施能老师文章中使用的资料是英国气象局逐月海平面气压资料）

同时，将结果与国家气候中心监测网上的结果也不同。

施能老师夏季风指数定义如下：


在计算时，先计算海平面气压夏季平均，并转为dat格式，gs如下
------------------------------------------------------------------------------------------------------------------------
'reinit'
'sdfopen E:\data\mon\slp.mon.mean.nc'
'set gxout fwrite'
'set fwrite E:\workgroup\SN\slp_sum.dat'

time=66
i=1
while (i<=time)
'set lon 0 357.5'
'set lat -90 90'
aa=12*i-6
ab=12*i-4
'define u=ave(slp,t='aa',t='ab',1)'
say i
'd u'
i=i+1
endwhile
'disable'
'reinit'

----------------------------------------------------------------------------------------------------------
再计算夏季风指数，fortran程序如下：
----------------------------------------------------------------------------------------------------------
      PROGRAM INTRPREC
      PARAMETER NO=144,NA=73,NY1=66
      REAL DAT(NO,NA,NY1)
        real datb1(NO,NA,NY1),DA1(NA,NY1)
        real C(NY1),CC(NY1)
        real ave1,bzc1,std1,ave2,bzc2,std2
        integer i,j,k,l,t,n,m
*************************************************************
      OPEN(1,FILE='slp_sum.dat',form='binary')
        read(1)(((DAT(i,j,n),i=1,NO),j=1,NA),n=1,NY1)
        close(1)
       
*************************标准化******************************
        !第一次标准化
        do io=1,NO
                do ia=1,NA
                ave1=0.0
                        do iy=1,NY1
                        ave1=ave1+DAT(io,ia,iy)
                        enddo
                        ave1=ave1/66.0
                bzc1=0.0
                        do iy=1,NY1
                        bzc1=bzc1+(DAT(io,ia,iy)-ave1)**2
                        enddo
                        bzc1=sqrt(bzc1/66.0)
                        do iy=1,NY1
                        datb1(io,ia,iy)=(DAT(io,ia,iy)-ave1)/bzc1
                        enddo
                enddo
        enddo
       
*************************计算EASM******************************
       
        do ia=1,NA
                do iy=1,NY1
                DA1(ia,iy)=datb1(45,ia,iy)-datb1(65,ia,iy) !110-160
                enddo               
        enddo
             
                do iy=1,NY1
                        C(iy)=0.0
                        do ia=45,57,2      !累加
                        C(iy)=C(iy)+DA1(ia,iy)
                        enddo
                enddo
       
        call BZH(c,cc)
       
        open(999,file='sum_EASM.txt')
        do i=1,NY1
      write(999,*) 1947+i,CC(i),-1.0*CC(i)
        enddo
        close(999)

        END


*************************标准化***********************************
      subroutine BZH(c,cc)
      INTEGER,PARAMETER::N=66
          INTEGER i
      REAL C(n),ave,CC(N),std

      ave=0.0
      do i=1,N
          ave=ave+c(i)
        enddo
        ave=ave/N

        std=0.0
        do i=1,N
          std=std+(c(i)-ave)**2
        enddo
        std=sqrt(std/N)

        do i=1,N
                CC(i)=(c(i)-ave)/std
        enddo
             
        END

----------------------------------------------------------------------------------------
气候中心提供结果：

按以上程序计算结果（没有乘-1）


结果差异比较大，不知道有谁计算过，麻烦帮忙看看错误在哪里。。。不胜感激！

最终正确结果脚本及程序 将另发帖与大家共享

timeover

orz~~~过节前木有人逛论坛了么- -

赵尼玛

有没有施老师的学生，来帮帮忙

wc619lf

我想，问题可能出在你第一次标准化上，因为第一次标准化只是针各个纬度带的不均匀做的，不需要加时间的66年进去标准化，只做纬度循环。
只是个人一点想法，仅供参考

timeover

wc619lf 发表于 2013-9-18 11:05

                               
登录/注册后可看大图

我想，问题可能出在你第一次标准化上，因为第一次标准化只是针各个纬度带的不均匀做的，不需要加时间的66年 ...

你的意思是，第一次标准化不是对每个格点自身时间序列进行标准化，而是做纬度序列的标准化？

timeover

wc619lf 发表于 2013-9-18 11:05

                               
登录/注册后可看大图

我想，问题可能出在你第一次标准化上，因为第一次标准化只是针各个纬度带的不均匀做的，不需要加时间的66年 ...

我试下看看哈~~thanks for your help~~

wc619lf

timeover 发表于 2013-9-18 11:12

                               
登录/注册后可看大图

你的意思是，第一次标准化不是对每个格点自身时间序列进行标准化，而是做纬度序列的标准化？

嗯，从他的表述应该是这个意思，我也没计算过这个指数，但他说是第一次标准化是为了消除各个纬带的不均匀的

wc619lf

也就是说你可以把时间循环放在最外面

timeover

wc619lf 发表于 2013-9-18 11:19

                               
登录/注册后可看大图

也就是说你可以把时间循环放在最外面

C:\Documents and Settings\Harvey.Young\桌面\1.jpg


刚刚上传的有点问题，已修改

kongfeng0824

请问有现成的算好的比较公认的夏季风指数或者moonsoon的指数吗？或者每年对中国影响较大的副高指数这类数据。一年一个值的这种的。