{已解决}{求助}将海温日资料处理为侯距平异常

十字天阶

现有2001-2011年的海温日数据，维数是365x1440x720，精度为0.25x0.25。要求是把日平均数据先做侯平均，再做侯距平，最后挑选出冬季的侯距平进行分析。我已经想了好几个月了......用了无数种方法，方法如下：1.用grads批量读取再处理，每年72个侯，一共792个侯，能出图但是部分有等值线部分没有，且有图的图都一样。

程序为：
'reinit'
year=2001
while(year<=2011)
'sdfopen d:\sst\sst.day.mean.'%year%'.nc'
'set lat 20 55'
'set lon 120 225'
day=1
mon=1
while(mon<=12)
pen=1
while(pen<=6)
**判定
**1,3,5,7,8,10,12都是最后一侯6天
  if(mon=1|mon=3|mon=5|mon=7|mon=8|mon=10|mon=12)
    if(pen=6)
     'a=ave(sst,t='day',t='day'+5)'
     day=day+6
    endif
**非闰年2月28天,最后一侯3天；闰年第6侯4天
  if(mon=2)
    if(pen=6)
*能被4整除但不能被100整除或者能被400整除
     rc=math_fmod(year,4)
     rc1=math_fmod(year,100)
     rc2=math_fmod(year,400)
     if((rc=0&rc1!=0)|(rc2=0))
        'a=ave(sst,t='day',t='day'+3)'
         day=day+4
     else
        'a=ave(sst,t='day',t='day'+2)'
        day=day+3
     endif
    endif
  endif
  else
   'd ave(sst,t='day',t='day'+4)'
   day=day+5
  endif
'printim c:\download\errors\sst'pen'.'mon'.'year'.png white x1000 y600'
'c'
  pen=pen+1
endwhile
mon=mon+1
endwhile
year=year+1
endwhile
;

（以上程序有参考气象家园中的大佬写的）
2.用ncl批量读取再处理，结果是过一段时间运行一次报一个错误.....
上次是报超出变量数目，我delete之后还是，这次又是这个错误（学姐建议我用超算，于是申请了超算帐号，无奈数据太多，光是上传数据就用了好多天，现在还没传输完。用已经传输完的试做了下，还是报一样的错误。）
emmmm老师建议我一行一行检查输出，还是报这个错误....

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
diri="d:/sst/"
filepath=systemfunc("ls "+diri+"sst.day.mean.*")+".nc"
in=addfiles(filepath, "r")
ListSetType(in, "cat")
year=ispan(2001, 2011, 1)
it_s=20011201 ;起始年月
it_e=20110228 ;结束年月
time=in[:]->time
YYYYMMDD=cd_calendar(time, -2);转换成公历日期
rec_s=ind(it_s.eq.YYYYMMDD)
rec_e=ind(it_e.eq.YYYYMMDD)
sst=in[:]->sst(rec_s:rec_e,:,:)
sst_pentad=new((/180,720,1440/), float)
day=0
year=2001
do while (year.le.2008)
    mon=0
    do while (mon.le.2)
        pen=1
        do while(pen.le.6)
            ;12.1月最后一侯都为6天
            if ((mon.eq.0).or.(mon.eq.1))then
                if(pen.eq.6)then
                    day_e=day+5
                    sst_pentad=dim_avg_n_Wrap(sst(day:day_e,:,:),0)
                    sst_pen_ano=dim_rmvmean_n_Wrap(sst_pentad, 0)
                    day=day+6
                end if
                ;平年2月28天，最后一侯3天；闰年第6侯4天
            else if(mon.eq.2)then
                if(pen.eq.6)then
                    if(((year%4.eq.0).and.(year%100.ne.0)).or.(year%400.eq.0))then
                        day_e=day+3
                        sst_pentad=dim_avg_n(sst(day:day_e,:,:),0)
                        sst_pen_ano=dim_rmvmean_n_Wrap(sst_pentad, 0)
                        day=day+4
                    else
                        day_e=day+2
                        sst_pentad=dim_avg_n(sst(day:day_e,:,:),0)
                        sst_pen_ano=dim_rmvmean_n_Wrap(sst_pentad, 0)
                        day=day+3
                    end if
                end if
            else
                day_e=day+4
                sst_pentad=dim_avg_n(sst(day:day_e,:,:),0)
                sst_pen_ano=dim_rmvmean_n_Wrap(sst_pentad, 0)
                day=day+5
            end if
            pen=pen+1
        end do
        mon=mon+1
    end do
    year=year+1
end do
wks=gsn_open_wks("pdf","plot-sst-pentad-200901")
gsn_define_colormap(wks, "rainbow")
res        =True
res@gsnAddCycle  =True
res@gsnDraw =False
res@gsnFrame =False
res@gsnLeftString =""
res@gsnRightString =""
res@mpFillOn =False
res@mpCenterLonF =180
res@mpGeophysicalLineThicknessF =0.5
res@pmTickMarkDisplayMode ="Always"
res@mpGridAndLimbOn =True
res@mpGridLatSpacingF =15
res@mpGridLonSpacingF =15
res@mpGridLineDashPattern =2
res@mpGridLineThicknessF =0.5

res@cnFillOn =True
res@cnLinesOn =True
res@cnLineColor ="White"
res@cnLineThicknessF =0.3
res@cnLineLabelsOn =False
plot=gsn_csm_contour(wks, sst_pen_ano, res)
draw(plot)
end

3.我本来想用Fortran批量读取，结果一直安装不上netcdf....
4.最后想用grads一年一年处理，也是错误不断...
5.用grads把nc数据转为grd格式，再用Fortran一年一年读取，还是错误不断...
而这还是本科毕业论文的内容。。我是不是毕不了业了。。哭了，求大神帮忙吧

river

怎么解决的，也给大家分享下经验吧。版主什么的会来加分的哦