【已解决】关于GrADS求季节平均问题的阐述与讨论（求牛人来看）

爱知天气

       说一说昨天看到的关于用GrADS求时间平均的问题，（那个帖子的地址是http://bbs.06climate.com/forum.php?mod=viewthread&tid=22701&page=1#pid283450）我也有类似的问题想问一问，并且我用了论坛上[步绝尘_南山居士的方法试过了，并没有得到想要的结果（不知是不是我没有理解你的意思@步绝尘_南山居士一起来讨论），因此今天拿出来跟大家讨论一下。首先说明一下，因为本人比较愚钝，所以也没有弄清楚究竟是为什么，如果有牛人知道其中的道理，希望可以帮助大家解答一下，在此先谢过。帖内容可能很繁琐，比较长，耐心看完可能确实比较困难。言归正传，下面列出问题。

       我从网上下载了海温（sst）的nc资料，资料从1854年1月--2014年3月，共1923个时次，资料精度为2°*2°，现在我已经截取了1958年--1977年这20年的东印度洋的海温资料，并且差值成了2.5°*2.5°，经过画图检验是正确的（我挑一些个人觉得跟主题有关的ctl、gs以及fortran程序列出，之前数据处理的暂不列出，如有需要，我会在楼下补充)。

这里现列出1854--2014的东印度洋的海温ctl（sst.ctl）：

dset E:\BY\figure2\sst.grd

title EIOsst 1854-2014

undef -9.99e8

xdef 53 linear 70 2.5

ydef 21 linear -10 2.5

zdef 1 linear 0 1

tdef 1923 linear 00Z01JUN1854 1mo

vars 1

sst 0 0 sea surfcae tempreature

Endvars

然后我截取我要用的1958--1977年的资料，下面是gs（cutsst1958-1977.gs）：

'reinit'

'open E:\BY\figure2\sst.ctl'

'set fwrite E:\BY\figure2\sst1958-1977.grd'

'set gxout fwrite'

'set t 1249 1488'

'd sst'

'disable fwrite'

;

现在给出sst1958-1977.grd的ctl（sst1958-1977.ctl）：

dset E:\BY\figure2\sst1958-1977.grd

title EIOsst 1958-1977

undef -9.99e8

xdef 53 linear 70 2.5

ydef 21 linear -10 2.5

zdef 1 linear 0 1

tdef 240 linear 00Z01JAN1958 1mo        *时次已变为20yr*12mo=240

vars 1

sst 0 0 sea surfcae tempreature

endvars

       至此，已经得到需要的20年的资料。现在我要求这20年的每年夏季（6-10月）平均场。这里有两种办法：第一种就是我用的笨办法，分别把20年的夏季每一个月份资料分别提取出，生成5个grd文件，也就是20年的6月一个grd，20年的7月一个grd（以此类推），然后将这5个grd文件相加除以5，便得到每一年夏季平均。第二种便是在论坛里看到的，[步绝尘_南山居士说的方法，这是他给出的gs的模版：

tz=1173
while(tz<=1905)
'define s=ave(sst,t='tz',t='tz+2')'
'd s'
tz=tz+12
endwhile即用循环以及ave函数求每一年的夏季平均。但是这种方法求出的结果跟我的笨方法不一样，画出的图不同且随时次没有变化。

    现在分别给出两会总方法的具体步骤：

SⅠ：

我给出提取夏季每一个月的gs（笨办法，各位不要见笑）：

************collect the sst in specific month every year**************

'reinit'

'open E:\BY\figure2\sst1958-1977.ctl'

'set gxout fwrite'

'set fwrite e:\BY\figure2\20yrmean\sst1958-1977_6.grd'     *提取6月的资料

tt=6

while(tt<=240)

'set x 1 53'

'set y 1 21'

'set t 'tt''

'd sst'

tt=tt+12

endwhile

'disable fwrite'

;

'set fwrite e:\BY\figure2\20yrmean\sst1958-1977_7.grd'     *提取7月的资料

tt=7

while(tt<=240)

'set x 1 53'

'set y 1 21'

'set t 'tt''

'd sst'

tt=tt+12

endwhile

'disable fwrite'

;

'set fwrite e:\BY\figure2\20yrmean\sst1958-1977_8.grd'     *提取8月的资料

tt=8

while(tt<=240)

'set x 1 53'

'set y 1 21'

'set t 'tt''

'd sst'

tt=tt+12

endwhile

'disable fwrite'

;

'set fwrite e:\BY\figure2\20yrmean\sst1958-1977_9.grd'      *提取9月的资料

tt=9

while(tt<=240)

'set x 1 53'

'set y 1 21'

'set t 'tt''

'd sst'

tt=tt+12

endwhile

'disable fwrite'

;

'set fwrite e:\BY\figure2\20yrmean\sst1958-1977_10.grd'      *提取10月的资料

tt=10

while(tt<=240)

'set x 1 53'

'set y 1 21'

'set t 'tt''

'd sst'

tt=tt+12

endwhile

'disable fwrite'

;

这里我生成了5个grd，现在给出一个ctl（sst1958-1977_6.ctl），其他雷同：

dset E:\BY\figure2\20yrmean\sst1958-1977_6.grd

title sst 1958-1977 every JUNE

undef -9.99e08

xdef 53 linear 70 2.5

ydef 21 linear -10 2.5

zdef 1 levels 500

tdef 20 linear 00Z01JUN1958 1yr

vars 1

sst 1 0 Monthly T

Endvars

接下来，我就把每个文件都打开，求了平均，这里是gs（summerave_in20yr.gs）：

************calculate the average sst in summer every year**************

'reinit'

'open E:\BY\figure2\20yrmean\sst1958-1977_6.ctl'

'open E:\BY\figure2\20yrmean\sst1958-1977_7.ctl'

'open E:\BY\figure2\20yrmean\sst1958-1977_8.ctl'

'open E:\BY\figure2\20yrmean\sst1958-1977_9.ctl'

'open E:\BY\figure2\20yrmean\sst1958-1977_10.ctl'

'set gxout fwrite'

'set fwrite e:\BY\figure2\20yrmean\sst1958-1977_sum.grd'

'set t 1 20'

'd (sst.1+sst.2+sst.3+sst.4+sst.5)/5'

'disable fwrite'

;

至此，算出了夏季平均，给出ctl（sst1958-1977_sum.ctl）：

dset E:\BY\figure2\20yrmean\sst1958-1977_sum.grd

title sst 1958-1977 in summer

undef -9.99e08

xdef 53 linear 70 2.5

ydef 21 linear -10 2.5

zdef 1 levels 500

tdef 20 linear 00Z01JUN1958 1yr

vars 1

sst 1 0 Monthly T

Endvars

大功告成，现在可以画出20年的夏季平均图，每一年都不同。

接下来说论坛上的方法。

SⅡ：

对于这种方法，我先用fortran把每一年的夏季（6-10月）的资料提取出来，生成一个只有夏季资料的文件sst1958-1977JtoO.grd：

program main                  

implicit none

integer i,j,t,yr

integer,parameter::n=20     !20年时间，1958—1977年

real sst(53,21,240)

!!!把grd文件读入数组!!!

open(10,file='E:\BY\figure2\sst1958-1977.grd',form="binary")

read(10)(((sst(i,j,t),i=1,53),j=1,21),t=1,240)

close(10)

!!!提取所需数据重新写入新二进制文件，便于处理

open(11,file='E:\BY\figure2\data6to10\sst1958-1977JtoO.grd',form="binary")

do yr=1,n

      do t=(yr-1)*12+6,(yr-1)*12+10        ! 58年开始的6--10月

          do j=1,21

             do i=1,53

      if(sst(i,j,t)<-999.) then

      sst(i,j,t)=-9999.

  endif

                write(11)sst(i,j,t)

             end do

       end do

      end do

end do

close(11)

End

下面给出ctl文件（sst1958-1977JtoO.ctl）：

dset E:\BY\figure2\data6to10\sst1958-1977JtoO.grd

title sst 1958-1977

undef -9999

xdef 53 linear 70 2.5

ydef 21 linear -10 2.5

zdef 1 levels 500

tdef 100 linear 00Z01JUN1958 1mo           *每一年都是 6 7 8 9 10月

vars 1                                   *5*20=100

sst 0 0 sea surfcae tempreature

Endvars

现在给出用第二种方法编写的gs：

************calculate the average sst in summer every year**************

'reinit'

'open E:\BY\figure2\data6to10\sst1958-1977JtoO.ctl'

'set gxout fwrite'

'set fwrite e:\BY\figure2\20yrmean\sst1958-1977_sum.home.grd'

a1=1

while(a1<=100)

'define a2=ave(sst,t='a1',t='a1+4')'

'd a2'

a1=a1+5

endwhile

'disable fwrite'

;

这样生成grd后，给出ctl：

dset E:\BY\figure2\20yrmean\sst1958-1977_sum.home.grd

title sst 1958-1977 in summer

undef -9.99e08

xdef 53 linear 70 2.5

ydef 21 linear -10 2.5

zdef 1 levels 500

tdef 20 linear 00Z01JUN1958 1yr

vars 1

sst 1 0 Monthly T

Endvars

       至此，第二种方法完成，也可以画出图，但是20年的图完全一样，并且第一个时次与我的方法的第一个时次的图形也不一样。

第二种方法跟我之前想的另外一种方法有点像，暂时叫方法三（SⅢ）吧，我当时的结果也不对，也是20个时次的图相同，这里是那个gs：

'reinit'

'open E:\BY\figure2\data6to10\sst1958-1977JtoO.ctl'

'set fwrite e:\BY\figure2\20yrmean\sst1958-1977.grd'

'set gxout fwrite'

yr=1

while(yr<=20)

tt=(yr-1)*5+1

ttt=(yr-1)*5+5

'd ave(sst,t='tt',t='ttt')'

yr=yr+1

endwhile

'disable fwrite'

;

       三种方法都已经介绍完，我个人运行的结果，只有第一种方法可以把20年的图都输出，切每年不同；而后两种方法，20个时次的图都是一样的。我也很困惑，如果那位牛人知道其中道理，希望可以帮我们拨开迷雾。

刷牙

i=0
while(i<=20)
t1=i*12+6
t2=i*12+10
'display ave(sst,t='t1',t='t2')'
i=i+1
endwhile
'disable fwrite'
；
这样是不是就行了啊，细节你在捯饬倒持
帖子太详细了

才仔

了解的可以来回复一下

步绝尘_南山居士

我的这个语句是用来提取秋季平均海温的，用ctl打开之后，对于我自己的资料，61个时次，每年的图都是不一样的，你的图不随时次变化的原因，不是很清楚，你把画图的gs文件贴上来看下吧

刷牙

看的我胃疼，为啥非要用Fortran将数据提出来，麻烦不说，还容易出错，你直接用nc算不行吗
i=0
while(i<=50)
t1=i*12+1
t2=i*12+12
'display ave(hgt,t='t1',t='t2')'
i=i+1
endwhile
'disable fwrite'
我以前做过，这样做是没有错的，张张是不一样的，
你直接用nc试试吧

御风者

你自己再好好研究研究，帖子太长了。。。。

爱知天气

刷牙 发表于 2014-4-27 13:46
i=0
while(i

不行。。我帖子最后就是这个样子写的，我也不知道为什么

爱知天气

刷牙 发表于 2014-4-27 13:39
看的我胃疼，为啥非要用Fortran将数据提出来，麻烦不说，还容易出错，你直接用nc算不行吗
i=0
while(i

好。当时我想的是这样看着数据清楚，就这么做了

爱知天气

步绝尘_南山居士 发表于 2014-4-27 13:18
我的这个语句是用来提取秋季平均海温的，用ctl打开之后，对于我自己的资料，61个时次，每年的图都是不一样 ...

我也不知道为什么。不知道哪里出错了。感觉你说的就是对的，我以前也这么做，但不知道为什么还是出错。
open E:\BY\figure2\20yrmean\sst1958-1977_sum.home.ctl
set t 1 20
d sst.

是不是ave这个函数，t那里可以写表达式

刷牙

我fortran不太熟，所以没帮你看，你直接用插值后的数据算，别把6-10月的提取出来，再试试
能力有限，不好意思啊