where（）遇到的问题

追风的阳光

下面是我的脚本：请大家帮我看看，谢谢！
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/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
begin


  in  = addfile("/home/zssapr/WRF3.5.1/20110710/WRFV3/test/em_real/wrfout_d02_2011-07-05_00:00:00","r")
  FNM =in->FNM
  U =in->U
  ua = wrf_user_unstagger(U,U@stagger)
  V = in->V
  va = wrf_user_unstagger(V,V@stagger)
  lat = in->XLAT(0,:,0)
  lon = in->XLONG(0,0,:)


  scale=1.e05

  dv = new ( dimsizes(ua), typeof(ua), ua@_FillValue )
  dv@long_name  = "divengence"
  dv@units      = "scaled"
  dv = uv2dv_cfd (ua,va,lat,lon, 2)

  ud    = new ( dimsizes(ua), typeof(ua), ua@_FillValue )
  vd    = new ( dimsizes(va), typeof(va), va@_FillValue )

  copy_VarCoords(ua, ud )
  copy_VarCoords(ua, vd )

dv2uvg(dv,ud,vd) ; div  ==> divergent  wind components
  ud@long_name  = "Zonal Divergent Wind"
  ud@units      = ua@units
  vd@long_name  = "Meridional Divergent Wind"
  vd@units      = va@units
  dv2 = uv2dv_cfd (ud,vd,lat,lon, 2)


  vr = new ( dimsizes(ua), typeof(ua), ua@_FillValue )
  vr@long_name  = "vorticity"
  vr@units      = "scaled"
  vr = uv2vr_cfd (ua,va,lat,lon, 2) * scale

  uort    = new ( dimsizes(ua), typeof(ua), ua@_FillValue )
  vort    = new ( dimsizes(va), typeof(va), va@_FillValue )


  copy_VarCoords(ua, uort )
  copy_VarCoords(ua, vort )
  vr2uvg(vr,uort,vort) ; div  ==> rotational  wind components
  uort@long_name  = "Zonal Divergent Wind"
  uort@units      = ua@units
  vort@long_name  = "Meridional Divergent Wind"
  vort@units      = va@units
  vr2 = uv2vr_cfd (uort,vort,lat,lon, 2)


  svbit = new ( dimsizes(ua), typeof(ua), ua@_FillValue )

  svbit =where(vr2.le.0.and.dv2.le.0,9.e99,dv2/vr2)
  printVarSummary(svbit)

  copy_VarCoords(ua, svbit)
；copy_VarCoords(ua, vr)
；Svbit = dim_avg_n_Wrap(svbit(:, :, 18 :36 , 18 :45 ), (/2, 3/))
；log_Svbit=new ( dimsizes(Svbit), typeof(Svbit) )
；copy_VarCoords(Svbit,log_Svbit)
；log_Svbit =log10( abs(Svbit))


  wks  = gsn_open_wks("ps","svbit_udvr")
  gsn_define_colormap(wks,"BlueWhiteOrangeRed")
  res = True
  res@cnFillOn = False
  res@cnLinesOn       = True
  res@vpWidthF            = 0.5           ; change aspect ratio of plot
  res@vpHeightF           = 0.3
  res@cnSmoothingOn        = True
  res@cnConstFEnableFil =True
  res@cnLevelSelectionMode ="ExplicitLevels"
  res@cnLevels                =(/-8,-7,-6,-5,-4,-3,-2,-1,1,2,3,4,5,6,7,8/)

  plot = gsn_csm_contour(wks,Svbit(bottom_top|0:20,Time|:),res)
  end

longlivehj

追风的阳光 发表于 2014-5-15 15:28
谢谢，可以了，我还想问一下，俩个四维数组可以直接相除吗？比如脚本中的svbit=dv2/vr2

大小一致可以直接相除。
注意一下，where首先会计算参数里面的所有表达式，所以，如果除数包含0，就会出错。所以，最好在做除法之前把除数中的0设为缺测。

追风的阳光

不好意思。。又发grads版来了，编辑好问题后，还要再点击一下ncl版。。。。。

longlivehj

svbit =where(vr2.le.0.and.dv2.le.0,9.e99,dv2/vr2)
9.e99已经超过float的范围了，就是说溢出了，inf。改为下面的试试：
svbit =where(vr2.le.0.and.dv2.le.0,ua@_FillValue,dv2/vr2)

追风的阳光

longlivehj 发表于 2014-5-15 13:13
svbit =where(vr2.le.0.and.dv2.le.0,9.e99,dv2/vr2)
9.e99已经超过float的范围了，就是说溢出了，inf。改 ...

谢谢，可以了，我还想问一下，俩个四维数组可以直接相除吗？比如脚本中的svbit=dv2/vr2

追风的阳光

longlivehj 发表于 2014-5-15 15:34
大小一致可以直接相除。
注意一下，where首先会计算参数里面的所有表达式，所以，如果除数包含0，就会出 ...

好的，谢谢！