从wrfout里面提取某个高度的风场

YJane

最近在处理固定高度处的风场，编了一个NCL脚本，放上来分享一下。主要思路是从wrfout里面提取某个高度的风场上，并另存为.nc文件，方便在matlab里面处理。（我的wrfout是逐小时输出的）






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

begin

  TC_name="Lekima"

; Make a list of all files we are interested in
  DATADir = "/home/data/user1data/Result/Lekima/1/"
  FILES = systemfunc (" ls -1 " + DATADir + "wrfout_d02* ")
  numFILES = dimsizes(FILES)
  print("numFILES = " + numFILES)
  print(FILES)
  print (" ")

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

  a = addfiles(FILES+".nc","r")

  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  ntimes = dimsizes(times)         ; number of times in the file
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;



  u  = wrf_user_getvar(a,"ua",-1)      ; 3D U at mass points
  v  = wrf_user_getvar(a,"va",-1)      ; 3D V at mass points
  w  = wrf_user_getvar(a,"wa",-1)       ; w averaged to mass points                  [m/s]
  z  = wrf_user_getvar(a, "z",-1)        ; grid point height

  u10 = wrf_user_getvar(a,"U10",-1)    ; u at 10 m, mass point
  v10 = wrf_user_getvar(a,"V10",-1)    ; v at 10 m, mass point

  lat= wrf_user_getvar(a,"lat",-1)       ; get latitude
  lon= wrf_user_getvar(a,"lon",-1)       ; get longitude
  slp= wrf_user_getvar(a,"slp",-1)       ; get sea level presser [hPa]

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

; The specific height levels that we want the data interpolated to.
; And interpolate to these levels

  height_levels = (/ 1000./)            ; height levels to plot - in meter
  nlevels       = dimsizes(height_levels)     ; number of height levels

  u_plane  = wrf_user_intrp3d( u,z,"h",height_levels,0.,False)
  v_plane  = wrf_user_intrp3d( v,z,"h",height_levels,0.,False)
  w_plane  = wrf_user_intrp3d( w,z,"h",height_levels,0.,False)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  variable_dims   = dimsizes(v_plane)
  m               = variable_dims(1)
  n               = variable_dims(2)

; print(variable_dims)
; print(getvardims(v_plane))
;===================================================================
  do it = 0,ntimes-1,1             ; TIME LOOP
    print("Working on time: " + times(it) )

   setfileoption("nc", "Format",  "NetCDF4")
   fn = height_levels + "_" + TC_name+".nc"
   system("/bin/rm -f " + fn) ; remove if exists
   f = addfile(fn, "c")

;===================================================================
;Define the file dimensions, NOTE that both dimensions are unlimited.
   dimNames = (/"south_north","west_east","Time","u_plane1","v_plane1","w_plane1","u10","v10","slp"/)
   dimSizes = (/     -1      ,    -1     ,  -1  ,   -1     ,   -1     ,   -1     , -1  ,  -1 ,  -1 /)
   dimUnlim = (/    True     ,   True    , True ,  True    ,  True    ,  True    ,True , True, True/)
   filedimdef(f, dimNames, dimSizes, dimUnlim)

;===================================================================
   u_height = new((/numFILES, m, n/), float)
   u_height@name = "u_height"
   u_height!0 = "Time"
   u_height!1 = "south_north"
   u_height!2 = "west_east"
   filevardef(f, "u_height", typeof(u_height), getvardims(u_height))
   filevarattdef(f,"u_height", u_height)


   v_height = new((/numFILES, m, n/), float)
   v_height@name = "v_height"
   v_height!0 = "Time"
   v_height!1 = "south_north"
   v_height!2 = "west_east"
   filevardef(f, "v_height", typeof(v_height), getvardims(v_height))
   filevarattdef(f,"v_height", v_height)



   w_height = new((/numFILES, m, n/), float)
   w_height@name = "w_plane1"
   w_height!0 = "Time"
   w_height!1 = "south_north"
   w_height!2 = "west_east"
   filevardef(f, "w_height", typeof(w_height), getvardims(w_height))
   filevarattdef(f,"w_height", w_height)


   u10_h = new((/numFILES, m, n/), float)
   u10_h@name = "u10_h"
   u10_h!0 = "Time"
   u10_h!1 = "south_north"
   u10_h!2 = "west_east"
   filevardef(f, "u10_h", typeof(u10_h), getvardims(u10_h))
   filevarattdef(f,"u10_h", u10_h)

   v10_h = new((/numFILES, m, n/), float)
   v10_h@name = "v10_h"
   v10_h!0 = "Time"
   v10_h!1 = "south_north"
   v10_h!2 = "west_east"
   filevardef(f, "v10_h", typeof(v10_h), getvardims(v10_h))
   filevarattdef(f,"v10_h", v10_h)

   slp_h = new((/numFILES, m, n/), float)
   slp_h@name = "slp_h"
   slp_h!0 = "Time"
   slp_h!1 = "south_north"
   slp_h!2 = "west_east"
   filevardef(f, "slp_h", typeof(slp_h), getvardims(slp_h))
   filevarattdef(f,"slp_h", slp_h)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    do k = 0, ntimes(0)-1
    do j = 0, m - 1
    do i = 0, n - 1
       u_height(k, i, j) = u_plane(k, i, j)
       v_height(k, i, j) = v_plane(k, i, j)
       w_height(k, i, j) = w_plane(k, i, j)
       u10_h(k, i, j) = u10(k, i, j)
       v10_h(k, i, j) = v10(k, i, j)
       slp_h(k, i, j) = slp(k, i, j)
    end do
    end do
    end do
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    f->u_height = (/u_height/)
    f->v_height = (/v_height/)
    f->w_height = (/w_height/)
    f->u10_h   = (/u10_h/)
    f->v10_h   = (/v10_h/)
    f->slp_h   = (/slp_h/)

  end do        ; END OF TIME LOOP



end

opeveu

收藏保存。楼主的分享对我很有帮助。谢谢

opeveu

请问一下楼主，这个高度怎么确定是海拔高度还是距离地面高度呢？

YJane

opeveu 发表于 2019-12-23 10:16
请问一下楼主，这个高度怎么确定是海拔高度还是距离地面高度呢？

这个我也说不清楚，我一般认为是离地高度。

opeveu

YJane 发表于 2019-12-23 15:42
这个我也说不清楚，我一般认为是离地高度。

好的。谢谢楼主

Liujiaxin

你好，高度应该不是对地高度吧，应该减去地形高度才对吧？