WRF运算结果画图不知道哪里出错？

tianmeng

不知道问题错在哪里？

风之牧语

tianmeng 发表于 2014-10-29 15:32
;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics ...

我看楼主图的标题写的是70m的风场，但是你的程序却是读取的70hPa的风场！我建议你去画一下中低层的风场来看看！

tianmeng

;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.
;   Interpolating to specified height levels

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
;
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
files=systemfunc("ls -1 wrfreckm3_201312*")
  f=addfiles(files,"r")
  ListSetType(f,"cat")


; We generate plots, but what kind do we prefer?
; type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"
  wks_type = "png"
  wks_type@wkWidth = 2500
  wks_type@wkHeight = 2500
  wks = gsn_open_wks(wks_type,"plt_HeightLevel")
gsn_define_colormap(wks,"ncl_default")

; Set some basic resources
  res = True

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


; The specific height levels that we want the data interpolated to.
  height_levels = (/ 70./)   ; height levels to plot - in meter
  nlevels       = dimsizes(height_levels)     ; number of height levels

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

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need
times = wrf_user_getvar(f,"times",-1)  ; get all times in the file
   ntimes = dimsizes(times)         ; number of times in the file
   sum_u  = wrf_user_getvar(f,"ua",0)        ; u averaged to mass points
   sum_v  = wrf_user_getvar(f,"va",0)        ; v averaged to mass points
   p  = wrf_user_getvar(f, "pressure",0)
   z  = wrf_user_getvar(f, "z",0)
   lon2d = wrf_user_getvar(f,"XLONG",0)
   lat2d = wrf_user_getvar(f,"XLAT",0)

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

   do it = 1,ntimes-1,1

    u  = wrf_user_getvar(f,"ua",it)        ; u averaged to mass points
    v  = wrf_user_getvar(f,"va",it)        ; v averaged to mass point
sum_u=sum_u+u
    sum_v=sum_v+v
    delete(u)
    delete(v)
   end do

   avg_u=sum_u/ntimes
   avg_v=sum_v/ntimes

      height = height_levels
      result_u = wrf_user_intrp3d( avg_u,z,"h",height,0.,False)
      result_v = wrf_user_intrp3d( avg_v,z,"h",height,0.,False)
      result_uv = sqrt(result_u*result_u+result_v*result_v)
delete(avg_u)
      delete(avg_v)
      delete(sum_u)
      delete(sum_v)
      delete(result_u)
      delete(result_v)




  result_uv@long_name="average wind speed at 70m in December"

  ; Set some Basic Plot options
   res=True
   res@gsnAddCyclic = False
res@gsnMaximize             = True
   ;res@gsnDraw                 = False
   ;res@gsnFrame                = False
   res@mpMinLatF               = 34.
   res@mpMaxLatF               = 39.
   res@mpMinLonF               = 114.
   res@mpMaxLonF               = 123.
   res@mpDataBaseVersion       = "MediumRes"
   res@mpDataSetName           = "Earth..4"
   res@mpOutlineOn           = True         ; Turn on map outlines
   res@mpOutlineSpecifiers   = (/"China","Taiwan"/)       ;China:states
   res@mpAreaMaskingOn = True   ;cover it
   res@mpMaskAreaSpecifiers = (/"China","Taiwan"/)   ;China:states
   res@mpLandFillColor         = "white"
   res@mpInlandWaterFillColor  = "white"
res@mpOceanFillColor        = "white"
   res@mpOutlineBoundarySets   = "NoBoundaries"
   res@mpFillDrawOrder = "PostDraw"
   res@cnFillOn = True
   res@cnLinesOn = False


   plot = gsn_csm_contour_map(wks,result_uv,res)


end

又是那隻貓

楼主提问时能否给出，你想要解决什么问题，哪里没有解决，哪里遇到困难，错误信息是什么之类的。

tianmeng

又是那隻貓 发表于 2014-10-29 19:16
楼主提问时能否给出，你想要解决什么问题，哪里没有解决，哪里遇到困难，错误信息是什么之类的。

没有错误，画出来的图就是这样子的，但是这个图明显就不对呀

hxyj

风之牧语 发表于 2014-10-29 20:12
我看楼主图的标题写的是70m的风场，但是你的程序却是读取的70hPa的风场！我建议你去画一下中低层的风场来 ...

观察仔细！！！

tianmeng

风之牧语 发表于 2014-10-29 20:12
我看楼主图的标题写的是70m的风场，但是你的程序却是读取的70hPa的风场！我建议你去画一下中低层的风场来 ...

result_u = wrf_user_intrp3d( avg_u,z,"h",height,0.,False)
      result_v = wrf_user_intrp3d( avg_v,z,"h",height,0.,False)
不是的，我插值的是z，而z的值是70m啊

tianmeng

我找出问题了，如果直接用高度插值是有问题的，因为地形的关系，好多地方的山都比70m高

tianmeng

所以要算局地的高度再加上70m我去

风之牧语

tianmeng 发表于 2014-10-29 21:23
result_u = wrf_user_intrp3d( avg_u,z,"h",height,0.,False)
      result_v = wrf_user_intrp3d( avg_ ...

我觉得你还是先去看看wrf_user_intrp3d的用法吧！你这样用是不对的！