NCL逐小时台站数据转nc格式问题

kevinniu23

读取2010年-2020年5月-10月台站逐小时降水数据，然后报错：
fatal:New: The dimension size list contains missing values, can't determine size
fatal:["Execute.c":8637]:Execute: Error occurred at or near line 126 in file 1.ncl

请大佬指教应如何修改脚本？原脚本如下：
begin
; nrow = numAsciiRow("E:\研一\描红\数据\hebing\hebing\数据\station_nw.list")

FillValue = -999.9
  missingvalue = -999.9
  NS=3789
  nhour=24

  hour=ispan(0,23,1)
  hour!0="hour"
  hour@long_name = "Station Time Hour"
  hour@size = 24

  sn=ispan(1,3789,1)
  sn!0="sn"
  sn@long_name = "Station Sequence Number"
  sn@size = 3789
  printVarSummary(sn)
  data=asciiread("E:\practise\shuju\station_nw.list",-1,"string")

  id= str_get_field(data, 1, " ")
  id!0="sn"
  id&sn=sn
  id@units = "None"
  id@long_name = "Current ID of Corresponding Station"
  id@_FillValue = FillValue
  printVarSummary(id)

;   ; lat=int2flt(data1(:,1)/100)+int2flt(data1(:,1)%100)/60.
;   ; lat=int2flt(data1(:,1))
  lat1= str_get_field(data, 3, " ")
  lat=stringtofloat(lat1)
  lat!0="sn"
  lat&sn=sn
  lat@units = "degree"
  lat@long_name = "Current Latitude of Corresponding Station"
  lat@_FillValue = FillValue
  printVarSummary(lat)

  ; lon=new(NS,string,FillValue)
  lon1= str_get_field(data, 4, " ")
  lon=stringtofloat(lon1)
  lon!0="sn"
  lon&sn=sn
  lon@units = "degree"
  lon@long_name = "Current Longitude of Corresponding Station"
  lon@_FillValue = FillValue
  printVarSummary(lon)

  ; altitude=new(NS,string,FillValue)
  ; altitude=int2flt(data1(:,3))
  altitude1= str_get_field(data, 5, " ")
  altitude=stringtofloat(altitude1)
  altitude!0="sn"
  altitude&sn=sn
  altitude@units = "degree"
  altitude@long_name = "Current altitude of Corresponding Station"
  altitude@_FillValue = FillValue
  printVarSummary(altitude)

  delete(data)


year=ispan(2010,2020,11)
month=ispan(5,10,1)
count= 184
ntime=count
  time   = new(ntime,double)
  Year   = new(ntime,integer)
  Month  = new(ntime,integer)
  Day    = new(ntime,integer)
  Hour   = new(ntime,integer)
  Minute = new(ntime,integer)
  Second = new(ntime,integer)
  Hour = 0
  Minute = 0
  Second = 0
  count=0

  do i=0,dimsizes(year)-1          ;把183天的“年”“月”“日”依次赋值

       Year(count:count+183)=year(i)

       Month(count:count+30)=5
       Day(count:count+30)=ispan(1,31,1)

       Month(count+31:count+60)=6
       Day(count+31:count+60)=ispan(1,30,1)

       Month(count+61:count+91)=7
       Day(count+61:count+91)=ispan(1,31,1)

       Month(count+92:count+122)=8
       Day(count+92:count+122)=ispan(1,31,1)

       Month(count+123:count+152)=9
       Day(count+123:count+152)=ispan(1,30,1)

       Month(count+153:count+183)=10
       Day(count+153:count+183)=ispan(1,31,1)


  end do

  time1=Year*10000+Month*100+Day
  printVarSummary(time1)
  print(time1(:20))
  time_units = "days since 2010-5-1 7:00:0.0"
  time = ut_inv_calendar(Year,Month,Day,Hour,Minute,Second,time_units,0)
  time!0 = "time"
  time@long_name = "time"
  time@units = time_units
  time@actual_range = (/"20100501-20201031"/)



precipitation=new((/ntime,NS,nhour/),float,FillValue)

do i=0,dimsizes(year)-1
    do k=0,dimsizes(month)-1


    nrow = numAsciiRow("E:/practise/shuju/SURF_PRE_HOR_NW-"+year(i)+"0"+month(k)+".TXT")
    print(nrow)
    stn=new(nrow,string) ;报错位置
    yy=new(nrow,integer)
    mm=new(nrow,integer)
    dd=new(nrow,integer)
    prt=new((/nrow,nhour/),float,FillValue)

    data=asciiread("E:/practise/shuju/SURF_PRE_HOR_NW-"+year(i)+"0"+month(k)+".TXT",-1,"string")

     stn=str_get_field(data, 1, " ")
     printVarSummary(stn)
     print(stn(0:10))

     yy=stringtoint(str_get_field(data,5 , " "))
     mm=stringtoint(str_get_field(data,6 , " "))
     dd=stringtoint(str_get_field(data,7 , " "))

     do h=0,23
       prt(:,h)=stringtofloat(str_get_field(data,8+h, " "))
     end do
     printVarSummary(prt)

     delete(data)


    do j=0,NS-1
    ; do j=0,0
        IND=ind(stn .eq. id(j))   ;返回的是同一个站号 第一次出现 第2次出现 第x次出现的行号位置集合
        print(IND)
        if (.not. any(ismissing(IND))) then
            time2=yy(IND)*10000+mm(IND)*100+dd(IND)
            print(time2)
            index=get1Dindex(time1,time2)
            print(index)
            exit
            precipitation(index,j,:)=prt(IND,:)

            delete(time2)
            delete(index)
        end if
        delete(IND)

     end do
    delete(stn)
    delete(yy)
    delete(mm)
    delete(dd)
    delete(prt)   
    printVarSummary(precipitation)
;     ; write loop content
  end do

end do
  precipitation=where(precipitation.lt.0.0,FillValue,precipitation)
  precipitation!0="time"
  precipitation!1="sn"
  precipitation!2="hour"
  precipitation&time=time
  precipitation&sn=sn
  precipitation&hour=hour
  precipitation@long_name = "Hourly Precipitation Observation of 3789 Stations in China"
  precipitation@units = "1 mm/h"
  precipitation@_FillValue = FillValue
  precipitation@missing_value = missingvalue
  precipitation@statistic = "Hourly Mean"
  precipitation@level_desc = "Surface"
  precipitation@dataset = "China Meteorology Station Observation Data"
  precipitation@var_desc = "Hourly Precip. Obs"

  system("/bin/rm -f " + "E:/practise/shuju/prec2010.1h.3789_new.nc")
  fout = addfile("E:/practise/shuju/prec2010.1h.3789_new.nc", "c")  
  ; system("/bin/rm -f " + "E:\研一\描红\数据\hebing\hebing\数据\prec"+year(i)+".1h.3789_new.nc")
  ; fout = addfile("E:\研一\描红\数据\hebing\hebing\数据\prec"+year(i)+".1h.3789_new.nc", "c")

  setfileoption(fout,"DefineMode",True)
  fAtt = True
  fAtt@title = "Observed Precipitation by Rain Gauge"
  fAtt@source = "Meteorology Observation Network, CMA"
  fAtt@platform = "Observation"
  fileattdef(fout,fAtt)


  dimNames = (/"time","sn","hour"/)
  dimSizes = (/-1,NS,24/)
  dimUnlim = (/True,False,False/)
  filedimdef(fout,dimNames,dimSizes,dimUnlim)

  filevardef(fout,"time",typeof(time),getvardims(time))
  filevardef(fout,"sn",typeof(sn),getvardims(sn))
  filevardef(fout,"hour",typeof(hour),getvardims(hour))
  filevardef(fout,"lat",typeof(lat),getvardims(lat))
  filevardef(fout,"lon",typeof(lon),getvardims(lon))
  filevardef(fout,"altitude",typeof(altitude),getvardims(altitude))
  filevardef(fout,"precipitation",typeof(precipitation),getvardims(precipitation))

  filevarattdef(fout,"time",time)
  filevarattdef(fout,"sn",sn)
  filevarattdef(fout,"hour",hour)
  filevarattdef(fout,"lat",lat)
  filevarattdef(fout,"lon",lon)
  filevarattdef(fout,"altitude",altitude)
  filevarattdef(fout,"precipitation",precipitation)

  setfileoption(fout,"DefineMode",False)

  fout->time = (/time/)
  fout->sn = (/sn/)
  fout->hour = (/hour/)
  fout->lat = (/lat/)
  fout->lon = (/lon/)
  fout->altitude = (/altitude/)
  fout->precipitation = (/precipitation/)


end

station_nw.list是台站信息文件，数据格式：站号、纬度、经度、高度

降水数据每个月一个文件，数据格式：站号、纬度、经度、高度、年、月、日、07-08时降水、08-09时降水 ...... 06-07时降水（每小时1个数，北京时间）、质控码（每小时一个数，0为正确数据）