爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 31074|回复: 26

请教,用wrfout的最小slp画台风路径出现问题

[复制链接]

新浪微博达人勋

发表于 2019-9-4 15:47:40 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

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



密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2019-9-7 13:05:53 | 显示全部楼层
没人那我就自己回答了,北面应该有干扰,划定了一个小一点的区域求最小slp就没问题了。核心代码贴在下面。
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)
maxlat = 31.0
minlat = 5.0
maxlon = 130.0
minlon = 160.0
opt = True
   loc  = wrf_user_ll_to_ij(a,(/minlon,maxlon/),(/minlat,maxlat/),opt)
   latlon = wrf_user_ij_to_ll(a,loc(0,:),loc(1,:),opt)
   loc = loc-1
  slpnew = slp2d(loc(1,0):loc(1,1),loc(0,0):loc(0,1))
  dimss = dimsizes(slpnew)
ctime = t(nn+fi*6)
print(ctime)
; We need to convert 2-D array to 1-D array to find the minima.
    slp1dnew     = ndtooned(slpnew)
    slp1d = ndtooned(slp2d)
    ;smin(ifs) = minind(slp1d)
    smin(fi) = ind(slp1d.eq.min(slp1dnew))
    slpmin(fi) = min(slp1dnew)
; Convert the index for 1-D array back to the indeces for 2-D array.
    minij     = ind_resolve(ind(slp1d.eq.min(slpnew)),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
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2019-12-23 19:25:59 | 显示全部楼层
llwang 发表于 2019-9-7 13:05
没人那我就自己回答了,北面应该有干扰,划定了一个小一点的区域求最小slp就没问题了。核心代码贴在下面。
...

您好,我用这个脚本会报错,请问该如何解决?谢谢
fatal:["NclFile.c":2100]:Subscript out of range, error in subscript #0
fatal:["Execute.c":8637]:Execute: Error occurred at or near line 657 in file $NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl

fatal:["Execute.c":8637]:Execute: Error occurred at or near line 1127 in file $NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl

fatal:["Execute.c":8637]:Execute: Error occurred at or near line 59 in file WRF_track_4.ncl
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-1-3 12:45:18 | 显示全部楼层
real-slimshady 发表于 2019-12-23 19:25
您好,我用这个脚本会报错,请问该如何解决?谢谢
fatal:["NclFile.c":2100]:Subscript out of range, e ...

可能是你定的经纬度范围太大了,大于原本模式输出的范围,你缩小一点试一下。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-2-24 11:08:23 | 显示全部楼层
llwang 发表于 2020-1-3 12:45
可能是你定的经纬度范围太大了,大于原本模式输出的范围,你缩小一点试一下。

我想请问一下经纬度设置在脚本中的哪条指令
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-3-21 10:34:42 | 显示全部楼层
不错不错,借鉴了楼主的code,感谢
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-3-23 17:05:07 | 显示全部楼层
同借鉴了楼主的脚本,太感谢了!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-7-26 17:13:59 | 显示全部楼层
zzq125 发表于 2020-2-24 11:08
我想请问一下经纬度设置在脚本中的哪条指令

请问您的问题解决了吗
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-2-24 21:22:25 | 显示全部楼层
请问楼主,nn和fi分别代表什么物理量?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-3-3 17:20:09 | 显示全部楼层
wangzhen 发表于 2021-2-24 21:22
请问楼主,nn和fi分别代表什么物理量?

我的模式数据是一小时一输出,在一个文件里面包含很多个时刻的数据。nn和fi是为了每6小时选一个数据。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表