用NCL提取wrfout中的SLP，出现很大震荡

Wanglethe · 发表于 2019-7-17 16:53:33

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

最近用NCL处理模拟出wrfout数据，画出的台风路径（用最小海平面气压slp定中心）尚可，画最大表面风速的时间序列也正常，但画台风最小中心气压slp随时间变化的序列图时，发现整体虽然有先减小后增大的趋势，但某些时刻slp有很大的震荡。于是画了每个时刻的SLP平面图，发现震荡时刻对应的平面图上，虽然还能显现出台风环流（这是为什么定中心还可以），但slp的值都偏大很多，且整体比较混乱。

又模拟了两次，每次出来的wrfout画slp都出现这种振荡错误，而且用ncview查看psfc这个变量随时间变化很连续，应该可以排除模拟的问题。因为我WRF模拟采用的是vortex following，出问题的也都是最里层移动网格的slp场，外面固定网格提取出的slp都没有震荡，现在考虑：是否NCL提取slp的命令slp = wrf_user_getvar(a,"slp",-1)，当对应移动网格的设置时计算出现问题？请问各位大神，还有没有另外计算slp的方法呢？

来自群组: 南京大学风云英华

Wanglethe · 发表于 2019-7-17 17:07:42

NCL提取slp的命令wrf_user_getvar定义在$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl中，里面显示的计算方法是：

if( variable .eq. "slp" ) then
if(isfilevar(nc_file,"T")) then
;; first compute theta - function wrf_tk needs theta and pressure (Pa) on input
;; THEN compute sea level pressure, from qv, p (Pa), tk, z
if ( time .eq. -1 ) then
if(ISFILE) then
T = nc_file->T
P = nc_file->P
PB = nc_file->PB
QVAPOR = nc_file->QVAPOR
PH = nc_file->PH
PHB = nc_file->PHB
else
T = file_handle[:]->T
P = file_handle[:]->P
PB = file_handle[:]->PB
QVAPOR = file_handle[:]->QVAPOR
PH = file_handle[:]->PH
PHB = file_handle[:]->PHB
end if
else
if(ISFILE) then
T = nc_file->T(time_in,:,:,:)
P = nc_file->P(time_in,:,:,:)
PB = nc_file->PB(time_in,:,:,:)
QVAPOR = nc_file->QVAPOR(time_in,:,:,:)
PH = nc_file->PH(time_in,:,:,:)
PHB = nc_file->PHB(time_in,:,:,:)
else
T = file_handle[:]->T(time_in,:,:,:)
P = file_handle[:]->P(time_in,:,:,:)
PB = file_handle[:]->PB(time_in,:,:,:)
QVAPOR = file_handle[:]->QVAPOR(time_in,:,:,:)
PH = file_handle[:]->PH(time_in,:,:,:)
PHB = file_handle[:]->PHB(time_in,:,:,:)
end if
end if
T = T + 300.
P = P + PB
QVAPOR = QVAPOR > 0.000
PH = ( PH + PHB ) / 9.81
z = wrf_user_unstagger(PH,PH@stagger)
tk = wrf_tk( P , T ) ; calculate TK
[b] slp = wrf_slp( z, tk, P, QVAPOR ) ; calculate slp[color=Red][/color][/b]
delete_VarAtts(T,(/"description","units"/))
copy_VarAtts(T,slp)

复制代码

所以其实ncl计算slp是通过提取wrfout中的z, tk, P, QVAPOR计算得出，tk是通过 P , T计算得出。

Wanglethe · 发表于 2019-7-17 17:11:24

附上我画时间序列代码中提取slp的部分，其实十分简单，画slp平面图也是这几步。

mslp = new(nt,float)
mwnd = new(nt,float)
do n=0,nt-1
slp = wrf_user_getvar(a1,"slp",n)
mslp(n) = min(slp)
u10 = wrf_user_getvar(a1,"U10",n) ; u at 10 m, mass point
v10 = wrf_user_getvar(a1,"V10",n) ; v at 10 m, mass point
wnd = sqrt(u10^2+v10^2)/0.5144 ;compute and convert to 'knot' units
mwnd(n) = max(wnd)
end do

复制代码

chensiqi · 发表于 2019-7-18 09:12:50

我之前用过一个脚本直接根据中心最低气压提取的台风路径比你这个震荡的还要厉害，当时看了整体气压的变化是连续的，但是在模拟台风后期，d01出现了又一个小低压，你是否方便把所有脚本都放上来，或许是，台风中心最低气压转为经纬度时出了错误

chensiqi · 发表于 2019-7-18 09:14:57

wrf中提取台风中心最低气压是根据最低气压以及梯度来计算的，这是世界上比较普遍流行的方法，所以建议楼主，要么可以写一个中心梯度的脚本，要么就直接用wrf输出的ATCF

Wanglethe · 发表于 2019-7-18 09:45:32

chensiqi 发表于 2019-7-18 09:12
我之前用过一个脚本直接根据中心最低气压提取的台风路径比你这个震荡的还要厉害，当时看了整体气压的变化是 ...

我这个问题应该不是另外低压的干扰，因为自己根据最小slp定的路径还比较符合，但移动网格中slp的值变化是不连续的，其实直接用ATCF画台风路径和强度变化很方便，我担心的是用wrf_user_getvar提取slp出错，会不会其他变量也出错呢？
附上自己提取最小slp定台风中心的脚本：

begin
;**************read best track data**********************************************
f0 ="/data/data1/xwang/Katrina/try_ncl/Katrina_track/best_track_Katrina(NHC).txt"
btdata = asciiread(f0,-1,"string")
; length1 = dimsizes(data1)-1 ;Skip first row of "data1" because it's just a header line
; print(length1)
lat0 = stringtofloat(str_get_field(btdata(1::), 5," ")) ;Use a space (" ") as a delimiter in str_get_field.
lon0 = -stringtofloat(str_get_field(btdata(1::), 6," ")) ;The first field is field=1 (unlike str_get_cols, in which the first column is column=0).
p0 = stringtofloat(str_get_field(btdata(1::), 8," ")) ; data(1::)----the second line to the last line(not sure)
printVarSummary(lat0)
;*************read wrfout d01 data**************************************************
a1 =addfile("/data/data1/xwang/WRF/WRF/test/em_real/wrfout2/wrfout_d01_2005-08-27_00:00:00","r")
times = wrf_user_getvar(a1,"times",-1) ; get all times in the file
nt = dimsizes(times) ; number of times in the file
print(nt)
lat1 = new(nt,float)
lon1 = new(nt,float)
do n=0,nt-1
ps = wrf_user_getvar(a1,"slp",n)
nps = dimsizes(ps) ; nps=(/*,*/),nps(0)---south_north(lat); nps(1)---west_east(lon)
; print(nps)
psmin = ps(0,0)
do i=0,nps(0)-1
do j=0,nps(1)-1
if(ps(i,j).lt.psmin)then
psmin =ps(i,j)
center_x=j
center_y=i
end if
end do
end do
lat1(n) = a1->XLAT(n,center_y,center_x)
lon1(n) = a1->XLONG(n,center_y,center_x)
end do
printVarSummary(lon1)
delete([/times,nt,ps,nps,psmin,center_x,center_y/]) ;!!!Must exist! If not, the d02 reading will go wrong.
;*************read wrfout d02 data**************************************************
a2 =addfile("/data/data1/xwang/WRF/WRF/test/em_real/wrfout2/wrfout_d02_2005-08-27_00:00:00","r")
times = wrf_user_getvar(a2,"times",-1) ; get all times in the file
nt = dimsizes(times) ; number of times in the file
print(nt)
lat2 = new(nt,float)
lon2 = new(nt,float)
do n=0,nt-1
ps = wrf_user_getvar(a2,"slp",n)
nps = dimsizes(ps) ; nps=(/*,*/),nps(0)---south_north(lat); nps(1)---west_east(lon)
; print(nps)
psmin = ps(0,0)
do i=0,nps(0)-1
do j=0,nps(1)-1
if(ps(i,j).lt.psmin)then
psmin =ps(i,j)
center_x=j
center_y=i
end if
end do
end do
lat2(n) = a2->XLAT(n,center_y,center_x)
lon2(n) = a2->XLONG(n,center_y,center_x)
end do
printVarSummary(lon2)
;***************plot Katrina track*****************************
wks = gsn_open_wks("png","Katrina_track(wrfout)") ; send graphics to PNG file
res = True
res@gsnDraw = False ; don't draw
res@gsnFrame = False ; don't advance frame
res@gsnMaximize = True
res@mpGridLineDashPattern = 1 ; lat/lon lines dashed
res@mpGridLatSpacingF = 5
res@mpGridLonSpacingF = 5
res@mpGridAndLimbDrawOrder = "Postdraw" ; Draw grid first
res@mpGridLineColor = "black"
res@mpGridMaskMode = "MaskLand"
res@mpDataBaseVersion = "MediumRes" ; better and more map outlines
res@mpDataSetName = "Earth..4"
res@mpOutlineBoundarySets = "AllBoundaries"
res@mpOutlineOn = True
res@mpGeophysicalLineThicknessF = 2.
res@mpNationalLineThicknessF = 2.
res@mpFillOn = False
res@mpPerimOn = True
;res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
res@pmTickMarkDisplayMode = "Always"
res@mpLimitMode = "LatLon" ; select subregion
res@mpMinLatF = 20
res@mpMaxLatF = 45
res@mpMinLonF = -100
res@mpMaxLonF = -70
res@tmYROn = False ; turn off right and top tickmarks
res@tmXTOn = False
res@tiMainString = "Katrina track" ; title
res@tiMainFontHeightF = 0.02
map = gsn_csm_map(wks,res) ; create map
; Set up some legend resources.
lgres = True
lgres@lgLineColors = (/"black","red","blue"/)
lgres@lgLineThicknessF = 4.
lgres@lgLabelFontHeightF = .08 ; set the legend label font thickness
lgres@vpWidthF = 0.15 ; width of legend (NDC)
lgres@vpHeightF = 0.15 ; height of legend (NDC)
lgres@lgMonoDashIndex = True
lgres@lgPerimColor = "black" ; draw the box perimeter in orange
lgres@lgPerimThicknessF = 3.0 ; thicken the box perimeter
labels = (/"best_track","wrfout_d01","wrfout_d02"/)
; ; Set up some legend resources.
; lgres = True
; lgres@lgLineColors = (/"black","red"/)
; lgres@lgLineThicknessF = 4.
; lgres@lgLabelFontHeightF = .08 ; set the legend label font thickness
; lgres@vpWidthF = 0.15 ; width of legend (NDC)
; lgres@vpHeightF = 0.15 ; height of legend (NDC)
; lgres@lgMonoDashIndex = True
; lgres@lgPerimColor = "black" ; draw the box perimeter in orange
; lgres@lgPerimThicknessF = 3.0 ; thicken the box perimeter
; labels = (/"best_track","wrfout_d01"/)
; Create the legend.
lbid = gsn_create_legend(wks,3,labels,lgres) ; create legend
; Set up resources to attach legend to map.
amres = True
amres@amParallelPosF = 0.4 ; positive move legend to the right
amres@amOrthogonalPosF = -0.32 ; positive move the legend down
annoid1 = gsn_add_annotation(map,lbid,amres) ; attach legend to plot
; Add text of every 6 hours
; txres = True
; txres@txFontHeightF = 0.015
; txres@txFontColor = "black"
; text1 = gsn_add_text(wks,map,"06",xp(0,0)+0.1,yp(0,0)+0.1,txres)
; text2 = gsn_add_text(wks,map,"12",xp(0,1)+0.15,yp(0,1),txres)
; text3 = gsn_add_text(wks,map,"18",xp(0,2)+0.15,yp(0,2),txres)
; text4 = gsn_add_text(wks,map,"00",xp(0,3), yp(0,3)+0.1,txres)
; Add trajectory lines.
pres = True ; polyline resources
pres@gsLineThicknessF = 6.0 ; line thickness
pres@gsLineColor = "black"
line0 = gsn_add_polyline(wks,map,lon0(:),lat0(:),pres) ; draw the traj
pres = True ; polyline resources
pres@gsLineColor = "red"
line1 = gsn_add_polyline(wks,map,lon1(:),lat1(:),pres) ; draw the traj
pres = True ; polyline resources
pres@gsLineColor = "blue"
line2 = gsn_add_polyline(wks,map,lon2(:),lat2(:),pres) ; draw the traj
; ; Add markers to the trajectories.
; mres0 = True ; marker resources for best track
; mres0@gsMarkerIndex = 16 ; marker style (filled circle)
; mres0@gsMarkerSizeF = 8.0 ; marker size
; mres0@gsMarkerColor = "black" ; maker color
; markers0=gsn_add_polymarker(wks,map,lon0(:),lat0(:),mres0)
; mres1 = True ; marker resources for best track
; mres1@gsMarkerIndex = 16 ; marker style (filled circle)
; mres1@gsMarkerSizeF = 8.0 ; marker size
; mres1@gsMarkerColor = "red" ; maker color
; markers1=gsn_add_polymarker(wks,map,lon1(:),lat1(:),mres1)
; mres2 = True ; marker resources for best track
; mres2@gsMarkerIndex = 16 ; marker style (filled circle)
; mres2@gsMarkerSizeF = 8.0 ; marker size
; mres2@gsMarkerColor = "blue" ; maker color
; markers2=gsn_add_polymarker(wks,map,lon2(:),lat2(:),mres2)
draw(map)
frame(wks)
end

复制代码

以及画最小slp及最大风速随时间变化的脚本：

;*************read min slp and max surface wind from wrfout_d03* data**************************************************
diri1 ="/data/data1/xwang/WRF/WRF/test/em_real/wrfout/"
fs1 =systemfunc(" ls "+diri1+"wrfout_d03* ")
inputfiles =addfiles(fs1,"r")
x =new(40,integer)
do i = 0, 39
x(i) =i+1
end do
allt =new(40,string)
minslp = new(40,float)
maxwnd = new(40,float)
do ifile = 0,4 ; not contain "wrfout_d03_2005-09-01"
a1 = inputfiles[ifile]
; times = wrf_user_getvar(a1,"times",-1) ; get all times in the file. No attributes of variables in wrfout*.nc
; nt = dimsizes(times) ; number of times in the file
; ;____wrfout*.nc中变量的time该coordinate并无赋值--> 该方法不可取！_____________________
; slp = wrf_user_getvar(a1,"slp",-1)
; printVarSummary(slp)
; times = slp&time ; get time information(type:float)
; nt = dimsizes(times) ; number of times in the file
;____________________________________________________________________________________
;____读取wrfout*.nc中的times，格式为string，不能直接作为坐标轴变量，需与numeric values(如x)配合作为坐标___________
times_in_file = wrf_user_list_times(a1) ;!!!Get time information (type:string, 2005-08-27_00:00:00)
day = str_get_cols(times_in_file,8,9) ;Strip out the day
hour = str_get_cols(times_in_file,11,12) ;Strip out the hour
times =(/hour+"/"+day/)
printVarSummary(times)
print(times)
nt = dimsizes(times) ; number of times in the file
;__Available in version 6.4.0 and later________________________________________________
; instring = "time"
; format = "%D_%H"
; times = cd_inv_string(instring, format)
;______________________________________________________________________________________
mslp = new(nt,float)
mwnd = new(nt,float)
do n=0,nt-1
slp = wrf_user_getvar(a1,"slp",n)
mslp(n) = min(slp)
u10 = wrf_user_getvar(a1,"U10",n) ; u at 10 m, mass point
v10 = wrf_user_getvar(a1,"V10",n) ; v at 10 m, mass point
wnd = sqrt(u10^2+v10^2)/0.5144 ;compute and convert to 'knot' units
mwnd(n) = max(wnd)
end do
istrt =ifile*nt
iend =(ifile+1)*nt-1
allt(istrt:iend) =times(:)
minslp(istrt:iend) =mslp(:)
maxwnd(istrt:iend) =mwnd(:)
end do
printVarSummary(allt)
print(allt)
printVarSummary(minslp)
printVarSummary(maxwnd)
;***************plot*****************************
wks = gsn_open_wks("png","Minimum_SLP") ; send graphics to PNG file
res = True
;---Titles
res@tiMainString = "Minimum SLP of Katrina"
res@tiMainFontHeightF = 0.02
res@tiYAxisFontHeightF = 0.02
res@tiXAxisString = "Time/Date (UTC)"
res@tiYAxisString = "Minimum SLP (hPa)"
;---TickMark
res@tmXBMode = "Explicit" ; Define own tick mark labels. "Automatic","Mannual","Explicit"
res@tmXBValues = x(0:39:8) ; location of explicit labels(x轴实际变量). "Explicit"
res@tmXBLabels = allt(0:39:8) ; labels are the locations(x轴显示的坐标). "Explicit"
res@tmLabelAutoStride = True ; control the step automatically
res@tmXTOn = False ; turn off the top tick marks
res@xyLineThicknesses = 2 ; increase line thickness
res@xyLineColor = "blue" ; set line color
plot = gsn_csm_xy (wks,x(:),minslp(:),res) ; create plot
delete([/wks,res/])
;__________________________________________________
wks = gsn_open_wks("png","Maximum_surfacewnd") ; send graphics to PNG file
res = True
;---Titles
res@tiMainString = "Maximum Surface Wind of Katrina"
res@tiMainFontHeightF = 0.02
res@tiYAxisFontHeightF = 0.02
res@tiXAxisString = "Time/Date (UTC)"
res@tiYAxisString = "Minimum SLP (hPa)"
res@tmXBMode = "Explicit" ; Define own tick mark labels.
res@tmXBValues = x(0:39:8) ; location of explicit labels(x轴实际变量)
res@tmXBLabels = allt(0:39:8) ; labels are the locations(x轴显示的坐标)
res@tmXTOn = False ; turn off the top tick marks
res@xyLineThicknesses = 2 ; increase line thickness
res@xyLineColor = "red" ; set line color
plot = gsn_csm_xy (wks,x(:),maxwnd(:),res) ; create plot
end

复制代码

Wanglethe · 发表于 2019-7-18 10:06:39

本帖最后由 Wanglethe 于 2019-7-18 10:08 编辑

其实我这个问题关键在于：为什么wrf_user_getvar从移动网格提取出的slp场在某些时刻出现很大的震荡误差？

E:\mywork\Katrina\1.png

Wanglethe · 发表于 2019-7-18 10:40:55

chensiqi 发表于 2019-7-18 09:14
wrf中提取台风中心最低气压是根据最低气压以及梯度来计算的，这是世界上比较普遍流行的方法，所以建议楼主 ...

谢谢您的建议，还想请问您一下，“wrf中提取台风中心最低气压是根据最低气压以及梯度来计算的”这个是从哪里看到的呢？WRF用户手册里似乎没有提到。

landachao · 发表于 2019-7-18 16:55:39

大神麻烦问一下，可以把你的县边界shp文件给我分享一下吗，特别需要麻烦了！

Wanglethe · 发表于 2019-7-18 18:51:33

landachao 发表于 2019-7-18 16:55
大神麻烦问一下，可以把你的县边界shp文件给我分享一下吗，特别需要麻烦了！

shp文件是啥？我不知道……

		自动登录	找回密码
密码			立即注册

[作图] 用NCL提取wrfout中的SLP，出现很大震荡

登录后查看更多精彩内容~