[ncl-talk] Problem in calculate a index in NCL

[grace]

Hi,all:    I am calculating a index and try to plot it.
     lindex=abs(sreh)*(qgsum+qisum+qssum)‍
    but after I write the script according the function above and try to run it,it just calculated the sreh and has no errors warning and wanted pictures appeared.
The bash window and script:



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/d032"+".nc","r")            




; We generate plots, but what kind do we prefer?
  type = "pdf"
; type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"plt_lightningindex")


  gsn_define_colormap(wks,"precip_11lev")
; Set some basic resources
  res = True
  res at MainTitle = "REAL-TIME WRF"
 ; res at gsnDraw      =  False                  
  ;res at gsnFrame     =  False


  mpres  = True  ; Map resources
  mpres at mpOutlineOn = False  ; Turn off map outlines
  mpres at mpFillOn    = False  ; Turn off map fill
  mpres at mpGridAndLimbOn = True
 ;res at mpProjection          = "Lambert"
  pltres = True ; Plot resources
  pltres at 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


    print("Working on time: " + times(it) )
    res at TimeLabel = times(it)   ; Set Valid time to use on plots
    sreh = wrf_user_getvar(a,(/"helicity","3000"/),it)    	; here 3km is specifically set - same as above


    cnres                      = res
    cnres at cnFillOn             = True
    cnres at cnSmoothingOn        = True
    cnres at cnSmoothingDistanceF = .005
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        
    if(isfilevar(a,"QGRAUP"))
	  qgraup = wrf_user_getvar(a,"QGRAUP",it )
      qgraup = qgraup*1000.
      qgraup at units = "g/kg"   
    end if


 if(isfilevar(a,"QICE"))


	  qice = wrf_user_getvar(a,"QICE",it )
      qice = qice*1000.
      qice at units = "g/kg"   
    end if


 if(isfilevar(a,"QSNOW"))
	  qsnow = wrf_user_getvar(a,"QSNOW",it )
      qsnow = qsnow*1000.
      qsnow at units = "g/kg"   
    end if
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


   
      opts = res
      opts at cnFillOn         = True
      opts at gsnSpreadColors  = False
      opts at ContourParameters       = (/ 1, 19, 2 /)
	  opts at gsnDraw      =  False                  
      opts at gsnFrame     =  False




      if (isvar("qgraup"))
        qgs  = qgraup(13 + (/0,1,2,3,4,5,6/),:,:)
        qgsum = dim_sum_n_Wrap(qgs, 0)
      end if


	   if (isvar("qice"))
        qis  = qice(13 + (/0,1,2,3,4,5,6/),:,:)
        qisum = dim_sum_n_Wrap(qis, 0)
      end if


	   if (isvar("qsnow"))
        qss  = qsnow(13 + (/0,1,2,3,4,5,6/),:,:)
        qssum = dim_sum_n_Wrap(qss, 0)
      end if


   lindex=abs(sreh)*(qgsum+qisum+qssum)
   contour = wrf_contour(a,wks,lindex,opts)
   plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
;>============================================================<
;                      add China map
;>------------------------------------------------------------<


     
  shp_name1    = "/home/Huanglei/map/China/diquJie_polyline.shp"


  lnres                  = True
  lnres at gsLineColor      = "gray25"
  lnres at 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 at gsLineThicknessF = 2.0       
  prres at gsLineColor = "black"
  plotcn3 = gsn_add_shapefile_polylines(wks,plot,shp_name2,prres)


   txres2  = True
   txres2 at txFont  = 10
   txres2 at txFontHeightF =0.01
   txres2 at 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.
  frame(wks)
   delete(opts)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


  end do        ; END OF TIME LOOP     


end‍



     Do you guys have some advice?
     any information will be appreciated

‍
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20140923/46a3fda1/attachment-0001.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/jpeg
Size: 75277 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20140923/46a3fda1/attachment-0001.jpe