用wrfout数据画图，之前都可以画出来，往里面添加了一段程序画叠加图后提示error

黄小仙儿

错误和脚本如下，看了很久没有看出问题，请教各位大神给点提示，很小的提示也好，多谢
错误：
fatal:Subscript out of range, error in subscript #0
fatal:An error occurred reading times
fatal:["Execute.c":8567]:Execute: Error occurred at or near line 49 in file /home/Huanglei/wrf_Cloudgraup_add_L.ncl
脚本：
;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
  a = addfile("/home/Huanglei/data/wrfoutd032013080521"+".nc","r")            

; We generate plots, but what kind do we prefer?
; type = "eps"
type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"20130806plt_Cloudgraup1020")
  gsn_define_colormap(wks,"precip_11lev")
; Set some basic resources
  res = True
  res@MainTitle = "REAL-TIME WRF"
; res@gsnDraw      =  False                  
  ;res@gsnFrame     =  False
  mpres  = True  ; Map resources
  mpres@mpOutlineOn = False  ; Turn off map outlines
  mpres@mpFillOn    = False  ; Turn off map fill
  mpres@mpGridAndLimbOn = True
  pltres = True ; Plot resources
  pltres@PanelPlot  = True   ; Tells wrf_map_overlays not to remove overlays

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; What times and how many time steps are in the data set?
  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  ntimes = dimsizes(times)         ; number of times in the file
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;  do it =1, ntimes-1,1        ; TIME LOOP
   it = 136
    print("Working on time: " + times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots
   ; print(it + (/0,1,2,3,4,5,6,7,8,9,10,11,12,13,14/))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        
    if(isfilevar(a,"QGRAUP"))
   qi = wrf_user_getvar(a,"QGRAUP",it )
      qi = qi*1000.
      qi@units = "g/kg"   
    end if
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    do level = 13,15,2      ; LOOP OVER LEVELS
      display_level = level + 2
      opts = res
      opts@cnFillOn         = True
      opts@gsnSpreadColors  = False
      opts@ContourParameters       = (/ 1, 19, 2 /)
      opts@PlotLevelID      = "Eta Level  " + display_level
   opts@gsnDraw      =  False                  
      opts@gsnFrame     =  False

      if (isvar("qi"))
        qis  = qi(level + (/0,1,2,3,4,5,6/),:,:)
        qisum = dim_sum_n_Wrap(qis, 0)
        contour = wrf_contour(a,wks,qisum,opts)
        plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
        delete(contour)
      end if
;>============================================================<
;                      add China map
;>------------------------------------------------------------<
     
  shp_name1    = "/home/Huanglei/map/China/diquJie_polyline.shp"
  lnres                  = True
  lnres@gsLineColor      = "gray25"
  lnres@gsLineThicknessF = 0.5   
id = gsn_add_shapefile_polylines(wks,plot,shp_name1,lnres)
  shp_name2    = "/home/Huanglei/map/China/cnmap/cnhimap.shp"
  prres=True
  prres@gsLineThicknessF = 2.0      
  prres@gsLineColor = "black"
  plotcn3 = gsn_add_shapefile_polylines(wks,plot,shp_name2,prres)
   txres2  = True
   txres2@txFont  = 10
   txres2@txFontHeightF =0.01
   txres2@txFontColor = "Blue"
   txdum1 =gsn_add_text(wks, plot, "Chengdu", 104.06,30.67, txres2)
  draw(plot)       ; This will draw the map and the shapefile outlines.

   delete(opts)
ascii_filename = "/home/Huanglei/data/20130806/utc11.txt"
  seismic = asciiread(ascii_filename,(/183,3/),"float")   
  
    y = seismic(:,0)  ; Column 1 of file contains X values.
    x = seismic(:,1)  ; Column 2 of file contains Y values.
    z = seismic(:,2)  ; Column 3 of file contains Z values.
  txres2  = True
  txres2@txFont  = 0.01
  txres2@txFontHeightF =0.01
  txres2@txFontColor = "Red"
idx = ind(z .gt. 0)
print(idx)
if .not. all(ismissing(idx))
    str = new(dimsizes(idx), "string")
    str = "+"
    txdum1 = gsn_add_text(wks, plot, str, x(idx),y(idx), txres2)
end if
txres2@txFontColor = "Blue"
idx := ind(z .lt. 0)
if .not. all(ismissing(idx))
    str := new(dimsizes(idx), "string")
    str = "-"   
    txdum2 = gsn_add_text(wks, plot, str, x(idx),y(idx), txres2)
end if
   draw(plot)
  frame(wks)
    end do      ; END OF LEVEL LOOP
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;  end do        ; END OF TIME LOOP     
end

longlivehj

能把49行标红吗？

黄小仙儿

longlivehj 发表于 2015-1-8 20:55
能把49行标红吗？

谢谢longlive大神，我已经找到原因了，是我太粗心了