- 积分
- 3230
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2020-3-19
- 最后登录
- 1970-1-1
![未绑定新浪微博用户 新浪微博达人勋](source/plugin/sina_login/img/gray.png)
|
20金钱
本帖最后由 HANSEN 于 2024-4-17 15:44 编辑
各位UU,最近在研究用WRF输出的高分辨数据进行涡度方程的诊断,一般来说因为摩擦或者次网格误差因素,方程的左右两边通常不等,会有残差项,但是查询一些文献,都指出残差项通常是小项,可以忽略不计。但是我在诊断的过程中发现,残差项的量级几乎与其他几项量级相当了,这让我百思不得其解。
我这里是随便拿了一个暴雨个例做测试,模拟2022年7月15日陕西一次暴雨事件,模式用FNL驱动,只设置一层网格,分辨率12km,垂直方向设置了53层,如图1给出了模拟的暴雨分布,WRF的模拟是成功的。模式时间步长为60s,同时我也输出了最高分辨率60s的wrfout数据,因为考虑到 偏R/偏t(R为相对涡度)希望t尽量趋近于0,我把wrfout的数据首先插值到P坐标系,然后插值成经纬度网格(0.1°),选取了暴雨中心区域计算涡度各项,如图2得到了各项的垂直剖面图,绿色虚线为残差项,紫色局地变化项,红色平流项,橙色散度项,另外还有一些小项不做说明。
问题来了:为什么这个残差项会这么大呢?
正常来说,我分辨率设置的已经非常高了,涡度方程的右边各项应该尽可能等于左边叭,我不能接受残差项的跟其他各项几乎相当!!
是不是我的脚本写错了,跪求各位大佬帮我看看!或者有没有做过这种分析的大佬传授一下经验
脚本如下:
begin
path = "~/WRF4.2-6/wrf/run/wrfout_d01_2022-07-14_18:00:00"
a = addfile(path,"r")
times = wrf_user_getvar(a,"times",-1) ; get times in the file
lat2d = a->XLAT(0,:,:)
lon2d = a->XLONG(0,:,:)
minlat = min(lat2d)
maxlat = max(lat2d)
minlon = min(lon2d)
maxlon = max(lon2d)
lon = fspan(minlon, maxlon, toint((maxlon-minlon)/0.1)+1)
lon!0 = "lon"
lon&lon = lon
lon@units = "degrees_east"
lon@long_name = "Longitude"
lat = fspan(minlat, maxlat, toint((maxlat-minlat)/0.1)+1)
lat!0 = "lat"
lat&lat = lat
lat@units = "degrees_north"
lat@long_name = "Latitude"
pp = new(13, "float")
pp = ispan(800, 500, 25)
pp@units = "hPa"
duration = 3
k= 0
do hh = 599,601
ua = wrf_user_getvar(a,"ua",hh)
va = wrf_user_getvar(a,"va",hh)
wa = wrf_user_getvar(a,"wa",hh)
prs = wrf_user_getvar(a,"pressure",hh)
u1 = wrf_user_interp_level(ua,prs,pp,False)
u2 = rcm2rgrid_Wrap(lat2d, lon2d, u1, lat, lon, 1)
delete(prs)
pp1 = pp * 100.
prs = conform_dims(dimsizes(u2), pp1, 0)
printVarSummary(prs)
print(prs(5,:,:))
v1 = wrf_user_interp_level(va,prs,pp,False)
v2 = rcm2rgrid_Wrap(lat2d, lon2d, v1, lat, lon, 1)
w1 = wrf_user_interp_level(wa,prs,pp,False)
w2 = rcm2rgrid_Wrap(lat2d, lon2d, w1, lat, lon, 1)
if (hh.eq.599) then
dims = dimsizes(u2)
VOR = new((/duration,dims(0),dims(1),dims(2)/),float)
UU = new((/duration,dims(0),dims(1),dims(2)/),float)
VV = new((/duration,dims(0),dims(1),dims(2)/),float)
WW = new((/duration,dims(0),dims(1),dims(2)/),float)
end if
UU(k,:,:,:) = u2
VV(k,:,:,:) = v2
WW(k,:,:,:) = w2
k = k + 1
end do
delete(prs)
DIV = uv2dv_cfd(UU, VV, lat, lon, 2)
VOR = uv2vr_cfd(UU, VV, lat, lon, 2)
dt = 60
dRdt = center_finite_diff_n(VOR, dt, False, 0, 0) * 1000000000.
pp1 = pp * 100.
prs = conform_dims(dimsizes(u2), pp1, 0)
nlat = dimsizes(lat)
nlon = dimsizes(lon)
nlevel = dimsizes(pp1)
dx = new((/nlevel,nlat,nlon/), float)
dy = dx
f = dx
dlon = (lon(1)-lon(0)) * 0.0174533
do nx = 0,nlat-1
dy(:,nx,:) = 6378388.*lat(nx)*0.0174533
f(:,nx,:) = coriolis_param(lat(nx))
do ny = 0,nlon-1
dx(:,nx,ny) = 6378388.*cos(0.0174533*lat(nx))*dlon * ny
end do
end do
delete([/ua,va,wa,u1,u2,v1,v2,w1,w2/])
do j = 1,1
div = DIV(j,:,:,:)
vor = VOR(j,:,:,:)
u = UU(j,:,:,:)
v = VV(j,:,:,:)
w = WW(j,:,:,:)
vor!0 = "level"
vor!1 = "lat"
vor!2 = "lon"
vor&level = pp
vor&lat = lat
vor&lon = lon
;--------- Term 1: -V·▽R ---------
dRdx = center_finite_diff_n(vor, dx, False, 0, 2)
dRdy = center_finite_diff_n(vor, dy, False, 0, 1)
term1 = -(u*dRdx + v*dRdy) * 1000000000.
;--------- Term 2: -w*dRdp ---------
term2 = - w * center_finite_diff_n(vor, prs, False, 0, 0) * 1000000000.
;--------- Term 3: -(R+f)*D ---------
term3 = - (vor+f) * div * 1000000000.
;--------- Term 4: -v*dfdy ---------
dfdy = center_finite_diff_n(f, dy, False, 0, 1)
term4 = -v*dfdy * 1000000000.
;--------- Term 5: dwdy*dudp-dwdx*dvdp ---------
dwdx = center_finite_diff_n(w, dx, False, 0, 2)
dwdy = center_finite_diff_n(w, dy, False, 0, 1)
dudp = center_finite_diff_n(u, prs, False, 0, 0)
dvdp = center_finite_diff_n(v, prs, False, 0, 0)
term5 = (dwdy*dudp-dwdx*dvdp) * 1000000000.
copy_VarCoords(vor, term1)
copy_VarCoords(vor, term2)
copy_VarCoords(vor, term3)
copy_VarCoords(vor, term4)
copy_VarCoords(vor, term5)
copy_VarCoords(vor, dRdt(0,:,:,:))
maxlat1 = 36.
minlat1 = 35.
maxlon1 = 110.
minlon1 = 109.
dRdt_avg = dim_avg_n_Wrap(dRdt(j,:,{minlat1:maxlat1},{minlon1:maxlon1}), (/1,2/))
term1_avg = dim_avg_n_Wrap(term1(:,{minlat1:maxlat1},{minlon1:maxlon1}), (/1,2/))
term2_avg = dim_avg_n_Wrap(term2(:,{minlat1:maxlat1},{minlon1:maxlon1}), (/1,2/))
term3_avg = dim_avg_n_Wrap(term3(:,{minlat1:maxlat1},{minlon1:maxlon1}), (/1,2/))
term4_avg = dim_avg_n_Wrap(term4(:,{minlat1:maxlat1},{minlon1:maxlon1}), (/1,2/))
term5_avg = dim_avg_n_Wrap(term5(:,{minlat1:maxlat1},{minlon1:maxlon1}), (/1,2/))
resi = dRdt_avg - term1_avg - term2_avg - term3_avg - term4_avg - term5_avg
total = term1_avg + term2_avg + term3_avg + term4_avg + term5_avg
copy_VarCoords(vor, resi)
wks = gsn_open_wks("png", "vertical_vor_eqt_test")
gsn_define_colormap(wks, "Cat12")
lev = toint(pp1 / 100.)
lev1 =lev(::-1)
res = True
res@gsnFrame = False
res@gsnDraw = False
res@gsnMaximize = True
res@tiMainString = ""
res@tiMainFontHeightF = 0.03
res@tiXAxisString = "10~S~-9~N~ "
res@tiYAxisString = "Pressure Level (hPa)"
res@tiXAxisFontHeightF = 0.015
res@tiYAxisFontHeightF = 0.015
res@trXMinF = -15
res@trXMaxF = 15
res@trYMaxF = 800
res@trYMinF = 500
res@xyLineThicknesses = 5
res@xyDashPatterns = 0
res@tmYLMode = "Explicit"
res@tmYLValues = lev(0::2)
res@tmYLLabels = tostring(lev1(0::2))
res@tmYLLabelFontHeightF = 0.015
res@tmYLMajorLengthF = 0.015
res@tmYLMajorOutwardLengthF = 0.015
res@tmXBLabelFontHeightF = 0.015
res@gsnXRefLine = 0
res@gsnXRefLineDashPattern = 1
res@tmXTOn = False
res@tmYROn = False
res@xyLineColors = 11
plot1 = gsn_csm_xy(wks,dRdt_avg,lev(::-1),res)
res@xyLineColors = 13
plot2 = gsn_csm_xy(wks,term1_avg,lev(::-1),res)
overlay(plot1, plot2)
res@xyLineColors = 12
plot3 = gsn_csm_xy(wks,term2_avg,lev(::-1),res)
overlay(plot1, plot3)
res@xyLineColors = 3
plot4 = gsn_csm_xy(wks,term3_avg,lev(::-1),res)
overlay(plot1, plot4)
res@xyLineColors = 5
plot6 = gsn_csm_xy(wks,term4_avg,lev(::-1),res)
overlay(plot1, plot6)
res@xyLineColors = 9
plot7 = gsn_csm_xy(wks,term5_avg,lev(::-1),res)
overlay(plot1, plot7)
res@xyDashPatterns = 1
res@xyLineColors = 7
plot5 = gsn_csm_xy(wks,resi,lev(::-1),res)
overlay(plot1, plot5)
draw(plot1)
frame(wks)
end do
end
|
-
图1 暴雨WRF模拟
-
图2 诊断的各项垂直分布
-
图3 涡度方程
|