ncep日资料合成

晓彤

应老板要求，最近参与用5月apo指数预测6月长江中下游降水的工作，然后一个尴尬的问题来了，五月份还没有过完肿么办？！！
我用日资料计算5月1日至？日（以22日为例）的均值替代月平均，重复相关工作。由于ncep日资料为一年一个文件，计算起来相当不方便，就想将每年的5月1日至22日（地面温度场）均值提取出来放到一个文件里，参考了家园里各位前辈的思路，用grads来处理。虽然比较简单，但中间还是绕了点弯路，这里分享下我的脚本和注意事项，权做参考 。
脚本如下：

'reinit'
'set gxout fwrite'
'set fwrite D:\NCC\air_500_1may_22may_1948_2019.grd'    生成一个文件时要放在循环外
i=1948
while(i<=2019)
'sdfopen D:\NCC\data\air.'i'.nc'     原始ncep日资料（144*73*17*366/365）
'set x 1 144'
'set y 1 73'
'set z 6'
'set t 1'
'define aa'i'=ave(air,time=1may,time=22may)'    用“time=1may,time=22may”，避免了还要另外循环讨论平闰年的问题
'd aa'i''
'close 1'    关键！不可省略，也不可用'reinit'替换，否则生成的新文件里所有时次的值均被第一个时次覆盖或文件里只有一个时次的值。
i=i+1
endwhile
'disable fwrite'
'reinit'
;
然后就可以计算apo指数啦！
这里我是分别计算500、400、300、250和200hPa的apo指数后求各层平均的，尝试过用grads直接生成包含各个高度层的多年的均值场文件，失败，用fortran把各高度层的多年均值场合并到一个文件，失败，所以只好用这个笨方法。据说python里可以按固定维（e.g. 第三维）直接合并，我不太熟悉python也没有试过。如果有碰到过类似问题的朋友，很想知道你们是怎么解决的，期待交流指正！