- 积分
- 1602
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2016-9-5
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
请教,用wrfout的最小slp画台风路径有很大的偏差。
路径图
ncl程序也贴在下面。谢谢大家
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/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
; DATES
;date = (/1512,1600,1612,1700,1712,1800,1812,1900/)
;ndate = dimsizes(date)
;sdate = sprinti("%4.0i",date)
; Experiment name (for legend)
EXP = (/"EXP_I"/) ; (/"EXP_I","EXP_II","EXP_III"/)
nexp = dimsizes(EXP)
; To get lat/lon info.
a = addfile("/home/2009091118/wrfout_d01_2009-09-11_18:00:00.nc","r")
lat2d = a->XLAT(0,:,:)
lon2d = a->XLONG(0,:,:)
dimll = dimsizes(lat2d)
nlat = dimll(0)
mlon = dimll(1)
t = wrf_user_getvar(a,"times",-1)
dimt = dimsizes(t)
;print(t)
;print(dimt)
ndate = 28
; Sea Level Pressure
slp = wrf_user_getvar(a,"slp",0)
dims = dimsizes(slp)
; Array for track
time = new(ndate,string)
imin = new(ndate,integer)
jmin = new(ndate,integer)
smin = new(ndate,integer)
slpmin = new(ndate,float)
newlon = new(ndate,float)
newlat = new(ndate,float)
; =======
; ndate
; =======
; fs = systemfunc("ls wrfout*00")
; nfs= dimsizes(fs)
;if(nfs .ne. ndate) then
; print("Check input data:"+nfs+" .ne. "+ndate)
;end if
nn = 6
; do ifs=0,nfs-1
do fi=0,27
; f = addfile(fs(ifs)+".nc","r")
; time(ifs) = wrf_user_list_times(f)
; print(time(ifs))
; slp2d = wrf_user_getvar(f,"slp",0)
slp2d = wrf_user_getvar(a,"slp",nn+fi*6)
ctime = t(nn+fi*6)
print(ctime)
; We need to convert 2-D array to 1-D array to find the minima.
slp1d = ndtooned(slp2d)
;smin(ifs) = minind(slp1d)
smin(fi) = minind(slp1d)
slpmin(fi) = min(slp1d)
; Convert the index for 1-D array back to the indeces for 2-D array.
minij = ind_resolve(ind(slp1d.eq.min(slp2d)),dims)
;imin(ifs) = minij(0,0)
;jmin(ifs) = minij(0,1)
imin(fi) = minij(0,0)
jmin(fi) = minij(0,1)
tlon = lon2d(minij(0,0),minij(0,1))
; print(time(ifs)+" : "+min(slp2d)+" ("+imin(ifs)+","+jmin(ifs)+")")
end do
;
print(smin)
; Graphics section
wks=gsn_open_wks("png","track") ; Open PS file.
gsn_define_colormap(wks,"BlGrYeOrReVi200") ; Change color map.
res = True
res@gsnDraw = False ; Turn off draw.
res@gsnFrame = False ; Turn off frame advance.
res@gsnMaximize = True ; Maximize plot in frame.
res@tiMainString = " ";"Hurricane Isabel" ; Main title
WRF_map_c(a,res,0) ; Set up map resources
; (plot options)
plot = gsn_csm_map(wks,res) ; Create a map.
; Set up resources for polymarkers.
gsres = True
gsres@gsMarkerIndex = 16 ; filled dot
;gsres@gsMarkerSizeF = 0.005 ; default - 0.007
cols = (/5,160,40/)
; Set up resources for polylines.
res_lines = True
res_lines@gsLineThicknessF = 3. ; 3x as thick
dot = new(ndate,graphic) ; Make sure each gsn_add_polyxxx call
line = new(ndate,graphic) ; is assigned to a unique variable.
; Loop through each date and add polylines to the plot.
do i = 0,ndate-2
res_lines@gsLineColor = cols(0)
xx=(/lon2d(imin(i),jmin(i)),lon2d(imin(i+1),jmin(i+1))/)
yy=(/lat2d(imin(i),jmin(i)),lat2d(imin(i+1),jmin(i+1))/)
line(i) = gsn_add_polyline(wks,plot,xx,yy,res_lines)
end do
do j = 0 ,ndate - 1
newlon(j) = lon2d(imin(j),jmin(j))
newlat(j) = lat2d(imin(j),jmin(j))
end do
asciiwrite ("track.txt" ,newlon +" "+ newlat)
asciiwrite ("slp.txt" ,slpmin)
lon1d = ndtooned(lon2d)
lat1d = ndtooned(lat2d)
; Loop through each date and add polymarkers to the plot.
do i = 0,ndate-1
print("dot:"+lon1d(smin(i))+","+lat1d(smin(i)))
gsres@gsMarkerColor = cols(0)
dot(i)=gsn_add_polymarker(wks,plot,lon1d(smin(i)),lat1d(smin(i)),gsres)
end do
; Date (Legend)
txres = True
txres@txFontHeightF = 0.015
txres@txFontColor = cols(0)
txid1 = new(ndate,graphic)
; Loop through each date and draw a text string on the plot.
; do i = 0, ndate-1
; txres@txJust = "CenterRight"
; ix = smin(i) - 4
; print("Eye:"+ix)
; if(i.eq.1) then
; txres@txJust = "CenterLeft"
; ix = ix + 8
; end if
; txid1(i) = gsn_add_text(wks,plot,sdate(i),lon1d(ix),lat1d(ix),txres)
; end do
; Add marker and text for legend. (Or you can just use "pmLegend" instead.)
txres@txJust = "CenterLeft"
; txid2 = new(nexp,graphic)
; pmid2 = new(nexp,graphic)
; do i = 0,nexp-1
; gsres@gsMarkerColor = cols(i)
; txres@txFontColor = cols(i)
; ii = ((/129,119,109/)) ; ilat
; jj = ((/110,110,110/)) ; jlon
; ji = ii*mlon+jj ; col x row
; pmid2(i) = gsn_add_polymarker(wks,plot,lon1d(ji(i)),lat1d(ji(i)),gsres)
; txid2(i) = gsn_add_text(wks,plot,EXP(i),lon1d(ji(i)+5),lat1d(ji(i)),txres)
; end do
draw(plot)
frame(wks)
end
|
|