请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 19267|回复: 18

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

[复制链接]

新浪微博达人勋

发表于 2019-7-17 16:53:33 | 显示全部楼层 |阅读模式

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

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

x
最近用NCL处理模拟出wrfout数据,画出的台风路径(用最小海平面气压slp定中心)尚可,画最大表面风速的时间序列也正常,但画台风最小中心气压slp随时间变化的序列图时,发现整体虽然有先减小后增大的趋势,但某些时刻slp有很大的震荡。于是画了每个时刻的SLP平面图,发现震荡时刻对应的平面图上,虽然还能显现出台风环流(这是为什么定中心还可以),但slp的值都偏大很多,且整体比较混乱。
2.png 1.png Katrina_track(wrfout)2.png
又模拟了两次,每次出来的wrfout画slp都出现这种振荡错误,而且用ncview查看psfc这个变量随时间变化很连续,应该可以排除模拟的问题。因为我WRF模拟采用的是vortex following,出问题的也都是最里层移动网格的slp场,外面固定网格提取出的slp都没有震荡,现在考虑:是否NCL提取slp的命令slp = wrf_user_getvar(a,"slp",-1),当对应移动网格的设置时计算出现问题?请问各位大神,还有没有另外计算slp的方法呢?


来自群组: 南京大学风云英华
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2019-7-17 17:07:42 | 显示全部楼层
NCL提取slp的命令wrf_user_getvar定义在$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl中,里面显示的计算方法是:
  1. if( variable .eq. "slp" ) then
  2.        if(isfilevar(nc_file,"T")) then
  3.          ;; first compute theta - function wrf_tk needs theta and pressure (Pa) on input
  4.          ;; THEN compute sea level pressure, from qv, p (Pa), tk, z
  5.          if ( time .eq. -1 ) then
  6.            if(ISFILE) then
  7.              T      = nc_file->T
  8.              P      = nc_file->P
  9.              PB     = nc_file->PB
  10.              QVAPOR = nc_file->QVAPOR
  11.              PH     = nc_file->PH
  12.              PHB    = nc_file->PHB
  13.            else
  14.              T      = file_handle[:]->T
  15.              P      = file_handle[:]->P
  16.              PB     = file_handle[:]->PB
  17.              QVAPOR = file_handle[:]->QVAPOR
  18.              PH     = file_handle[:]->PH
  19.              PHB    = file_handle[:]->PHB
  20.            end if
  21.          else
  22.            if(ISFILE) then
  23.              T      = nc_file->T(time_in,:,:,:)
  24.              P      = nc_file->P(time_in,:,:,:)
  25.              PB     = nc_file->PB(time_in,:,:,:)
  26.              QVAPOR = nc_file->QVAPOR(time_in,:,:,:)
  27.              PH     = nc_file->PH(time_in,:,:,:)
  28.              PHB    = nc_file->PHB(time_in,:,:,:)
  29.            else
  30.              T      = file_handle[:]->T(time_in,:,:,:)
  31.              P      = file_handle[:]->P(time_in,:,:,:)
  32.              PB     = file_handle[:]->PB(time_in,:,:,:)
  33.              QVAPOR = file_handle[:]->QVAPOR(time_in,:,:,:)
  34.              PH     = file_handle[:]->PH(time_in,:,:,:)
  35.              PHB    = file_handle[:]->PHB(time_in,:,:,:)
  36.            end if
  37.          end if
  38.          T = T + 300.
  39.          P = P + PB
  40.          QVAPOR = QVAPOR > 0.000
  41.          PH    = ( PH + PHB ) / 9.81
  42.          z = wrf_user_unstagger(PH,PH@stagger)

  43.          tk = wrf_tk( P , T )    ; calculate TK
  44.   [b]       slp   = wrf_slp( z, tk, P, QVAPOR )  ; calculate slp[color=Red][/color][/b]
  45.          delete_VarAtts(T,(/"description","units"/))
  46.          copy_VarAtts(T,slp)
复制代码

所以其实ncl计算slp是通过提取wrfout中的z, tk, P, QVAPOR计算得出,tk是通过 P , T计算得出。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-7-17 17:11:24 | 显示全部楼层
附上我画时间序列代码中提取slp的部分,其实十分简单,画slp平面图也是这几步。
  1.     mslp = new(nt,float)   
  2.     mwnd = new(nt,float)

  3.     do n=0,nt-1

  4.         slp = wrf_user_getvar(a1,"slp",n)
  5.         mslp(n) = min(slp)

  6.         u10 = wrf_user_getvar(a1,"U10",n)    ; u at 10 m, mass point
  7.         v10 = wrf_user_getvar(a1,"V10",n)    ; v at 10 m, mass point
  8.         wnd = sqrt(u10^2+v10^2)/0.5144 ;compute and convert to 'knot' units
  9.         mwnd(n) = max(wnd)

  10.     end do
复制代码
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-7-18 09:12:50 | 显示全部楼层
我之前用过一个脚本直接根据中心最低气压提取的台风路径比你这个震荡的还要厉害,当时看了整体气压的变化是连续的,但是在模拟台风后期,d01出现了又一个小低压,你是否方便把所有脚本都放上来,或许是,台风中心最低气压转为经纬度时出了错误
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-7-18 09:14:57 | 显示全部楼层
wrf中提取台风中心最低气压是根据最低气压以及梯度来计算的,这是世界上比较普遍流行的方法,所以建议楼主,要么可以写一个中心梯度的脚本,要么就直接用wrf输出的ATCF
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

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

我这个问题应该不是另外低压的干扰,因为自己根据最小slp定的路径还比较符合,但移动网格中slp的值变化是不连续的,其实直接用ATCF画台风路径和强度变化很方便,我担心的是用wrf_user_getvar提取slp出错,会不会其他变量也出错呢?
附上自己提取最小slp定台风中心的脚本:
  1. begin
  2. ;**************read best track data**********************************************
  3.   f0 ="/data/data1/xwang/Katrina/try_ncl/Katrina_track/best_track_Katrina(NHC).txt"
  4.   btdata  = asciiread(f0,-1,"string")

  5.   ; length1 = dimsizes(data1)-1 ;Skip first row of "data1" because it's just a header line
  6.   ; print(length1)
  7.   lat0 = stringtofloat(str_get_field(btdata(1::), 5," "))  ;Use a space (" ") as a delimiter in str_get_field.
  8.   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).
  9.   p0 = stringtofloat(str_get_field(btdata(1::), 8," "))    ; data(1::)----the second line to the last line(not sure)
  10.   printVarSummary(lat0)

  11. ;*************read wrfout d01 data**************************************************
  12. a1 =addfile("/data/data1/xwang/WRF/WRF/test/em_real/wrfout2/wrfout_d01_2005-08-27_00:00:00","r")


  13. times = wrf_user_getvar(a1,"times",-1)  ; get all times in the file
  14. nt = dimsizes(times)         ; number of times in the file
  15. print(nt)


  16. lat1 = new(nt,float)   
  17. lon1 = new(nt,float)


  18. do n=0,nt-1

  19.     ps = wrf_user_getvar(a1,"slp",n)
  20.     nps = dimsizes(ps) ; nps=(/*,*/),nps(0)---south_north(lat); nps(1)---west_east(lon)
  21.     ; print(nps)
  22.     psmin = ps(0,0)


  23.     do i=0,nps(0)-1
  24.       do j=0,nps(1)-1

  25.         if(ps(i,j).lt.psmin)then
  26.           psmin =ps(i,j)
  27.           center_x=j
  28.           center_y=i
  29.         end if

  30.       end do
  31.     end do

  32.     lat1(n)  = a1->XLAT(n,center_y,center_x)
  33.     lon1(n) = a1->XLONG(n,center_y,center_x)

  34. end do
  35. printVarSummary(lon1)
  36. delete([/times,nt,ps,nps,psmin,center_x,center_y/]) ;!!!Must exist! If not, the d02 reading will go wrong.


  37. ;*************read wrfout d02 data**************************************************
  38. a2 =addfile("/data/data1/xwang/WRF/WRF/test/em_real/wrfout2/wrfout_d02_2005-08-27_00:00:00","r")


  39. times = wrf_user_getvar(a2,"times",-1)  ; get all times in the file
  40. nt = dimsizes(times)         ; number of times in the file
  41. print(nt)


  42. lat2 = new(nt,float)   
  43. lon2 = new(nt,float)


  44. do n=0,nt-1

  45.     ps = wrf_user_getvar(a2,"slp",n)
  46.     nps = dimsizes(ps) ; nps=(/*,*/),nps(0)---south_north(lat); nps(1)---west_east(lon)
  47.     ; print(nps)
  48.     psmin = ps(0,0)


  49.     do i=0,nps(0)-1
  50.       do j=0,nps(1)-1

  51.         if(ps(i,j).lt.psmin)then
  52.           psmin =ps(i,j)
  53.           center_x=j
  54.           center_y=i
  55.         end if

  56.       end do
  57.     end do

  58.     lat2(n)  = a2->XLAT(n,center_y,center_x)
  59.     lon2(n) = a2->XLONG(n,center_y,center_x)

  60. end do
  61. printVarSummary(lon2)

  62. ;***************plot Katrina track*****************************
  63.    wks      = gsn_open_wks("png","Katrina_track(wrfout)")       ; send graphics to PNG file

  64.    res             = True
  65.    res@gsnDraw     = False                         ; don't draw
  66.    res@gsnFrame    = False                         ; don't advance frame
  67.    res@gsnMaximize = True

  68.    res@mpGridLineDashPattern  = 1                  ; lat/lon lines dashed
  69.    res@mpGridLatSpacingF      = 5
  70.    res@mpGridLonSpacingF      = 5
  71.    res@mpGridAndLimbDrawOrder = "Postdraw"          ; Draw grid first
  72.    res@mpGridLineColor        = "black"
  73.    res@mpGridMaskMode         = "MaskLand"

  74.    res@mpDataBaseVersion      = "MediumRes"       ; better and more map outlines
  75.    res@mpDataSetName          = "Earth..4"
  76.    res@mpOutlineBoundarySets  = "AllBoundaries"
  77.    res@mpOutlineOn            = True
  78.    res@mpGeophysicalLineThicknessF = 2.
  79.    res@mpNationalLineThicknessF    = 2.

  80.    res@mpFillOn              = False
  81.    res@mpPerimOn             = True
  82.    ;res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
  83.    res@pmTickMarkDisplayMode = "Always"

  84.    res@mpLimitMode = "LatLon"        ; select subregion
  85.    res@mpMinLatF   = 20
  86.    res@mpMaxLatF   = 45               
  87.    res@mpMinLonF   = -100
  88.    res@mpMaxLonF   = -70

  89.    res@tmYROn      = False     ; turn off right and top tickmarks
  90.    res@tmXTOn      = False

  91.    res@tiMainString      = "Katrina track"  ; title
  92.    res@tiMainFontHeightF = 0.02

  93.    map = gsn_csm_map(wks,res)                     ; create map

  94. ; Set up some legend resources.
  95.    lgres                    = True
  96.    lgres@lgLineColors       = (/"black","red","blue"/)
  97.    lgres@lgLineThicknessF   = 4.
  98.    lgres@lgLabelFontHeightF = .08            ; set the legend label font thickness
  99.    lgres@vpWidthF           = 0.15           ; width of legend (NDC)
  100.    lgres@vpHeightF          = 0.15            ; height of legend (NDC)
  101.    lgres@lgMonoDashIndex    = True   
  102.    lgres@lgPerimColor       = "black"       ; draw the box perimeter in orange
  103.    lgres@lgPerimThicknessF = 3.0            ; thicken the box perimeter
  104.    labels = (/"best_track","wrfout_d01","wrfout_d02"/)

  105.    ; ; Set up some legend resources.
  106.    ; lgres                    = True
  107.    ; lgres@lgLineColors       = (/"black","red"/)
  108.    ; lgres@lgLineThicknessF   = 4.
  109.    ; lgres@lgLabelFontHeightF = .08            ; set the legend label font thickness
  110.    ; lgres@vpWidthF           = 0.15           ; width of legend (NDC)
  111.    ; lgres@vpHeightF          = 0.15            ; height of legend (NDC)
  112.    ; lgres@lgMonoDashIndex    = True   
  113.    ; lgres@lgPerimColor       = "black"       ; draw the box perimeter in orange
  114.    ; lgres@lgPerimThicknessF = 3.0            ; thicken the box perimeter
  115.    ; labels = (/"best_track","wrfout_d01"/)

  116. ; Create the legend.
  117.    lbid   = gsn_create_legend(wks,3,labels,lgres)         ; create legend

  118. ; Set up resources to attach legend to map.
  119.    amres = True
  120.    amres@amParallelPosF   =  0.4           ; positive move legend to the right
  121.    amres@amOrthogonalPosF = -0.32                 ; positive move the legend down
  122.    annoid1 = gsn_add_annotation(map,lbid,amres)   ; attach legend to plot

  123. ; Add text of every 6 hours  
  124. ;   txres               = True
  125. ;   txres@txFontHeightF = 0.015        
  126. ;   txres@txFontColor   = "black"
  127. ;   text1 = gsn_add_text(wks,map,"06",xp(0,0)+0.1,yp(0,0)+0.1,txres)
  128. ;   text2 = gsn_add_text(wks,map,"12",xp(0,1)+0.15,yp(0,1),txres)
  129. ;   text3 = gsn_add_text(wks,map,"18",xp(0,2)+0.15,yp(0,2),txres)
  130. ;   text4 = gsn_add_text(wks,map,"00",xp(0,3),    yp(0,3)+0.1,txres)

  131. ; Add trajectory lines.
  132.    pres                  = True               ; polyline resources
  133.    pres@gsLineThicknessF = 6.0                ; line thickness
  134.    pres@gsLineColor      = "black"
  135.    line0 = gsn_add_polyline(wks,map,lon0(:),lat0(:),pres)      ; draw the traj

  136.    pres                  = True               ; polyline resources
  137.    pres@gsLineColor      = "red"
  138.    line1 = gsn_add_polyline(wks,map,lon1(:),lat1(:),pres)      ; draw the traj

  139.    pres                  = True               ; polyline resources
  140.    pres@gsLineColor      = "blue"
  141.    line2 = gsn_add_polyline(wks,map,lon2(:),lat2(:),pres)      ; draw the traj


  142. ; ; Add markers to the trajectories.
  143. ;    mres0                = True         ; marker resources for best track
  144. ;    mres0@gsMarkerIndex  = 16           ; marker style (filled circle)
  145. ;    mres0@gsMarkerSizeF  = 8.0          ; marker size
  146. ;    mres0@gsMarkerColor  = "black"      ; maker color
  147. ;    markers0=gsn_add_polymarker(wks,map,lon0(:),lat0(:),mres0)

  148. ;    mres1                = True         ; marker resources for best track
  149. ;    mres1@gsMarkerIndex  = 16           ; marker style (filled circle)
  150. ;    mres1@gsMarkerSizeF  = 8.0          ; marker size
  151. ;    mres1@gsMarkerColor  = "red"      ; maker color
  152. ;    markers1=gsn_add_polymarker(wks,map,lon1(:),lat1(:),mres1)

  153. ;    mres2                = True         ; marker resources for best track
  154. ;    mres2@gsMarkerIndex  = 16           ; marker style (filled circle)
  155. ;    mres2@gsMarkerSizeF  = 8.0          ; marker size
  156. ;    mres2@gsMarkerColor  = "blue"      ; maker color
  157. ;    markers2=gsn_add_polymarker(wks,map,lon2(:),lat2(:),mres2)

  158.    draw(map)                                          
  159.    frame(wks)                                         
  160.   
  161. end
复制代码


以及画最小slp及最大风速随时间变化的脚本:
  1. ;*************read min slp and max surface wind from wrfout_d03* data**************************************************
  2. diri1     ="/data/data1/xwang/WRF/WRF/test/em_real/wrfout/"
  3. fs1       =systemfunc(" ls "+diri1+"wrfout_d03* ")
  4. inputfiles =addfiles(fs1,"r")

  5. x =new(40,integer)
  6. do i = 0, 39
  7.   x(i) =i+1
  8. end do

  9. allt =new(40,string)
  10. minslp = new(40,float)
  11. maxwnd = new(40,float)


  12. do ifile = 0,4  ; not contain "wrfout_d03_2005-09-01"
  13.     a1  = inputfiles[ifile]
  14.     ; times = wrf_user_getvar(a1,"times",-1)  ; get all times in the file. No attributes of variables in wrfout*.nc
  15.     ; nt = dimsizes(times)         ; number of times in the file
  16.   
  17.    ; ;____wrfout*.nc中变量的time该coordinate并无赋值--> 该方法不可取!_____________________
  18.    ;  slp = wrf_user_getvar(a1,"slp",-1)
  19.    ;  printVarSummary(slp)
  20.    ;  times = slp&time ; get time information(type:float)
  21.    ;  nt = dimsizes(times)         ; number of times in the file
  22.    ;____________________________________________________________________________________
  23.    
  24.    ;____读取wrfout*.nc中的times,格式为string,不能直接作为坐标轴变量,需与numeric values(如x)配合作为坐标___________
  25.     times_in_file  = wrf_user_list_times(a1) ;!!!Get time information (type:string, 2005-08-27_00:00:00)
  26.     day = str_get_cols(times_in_file,8,9) ;Strip out the day
  27.     hour = str_get_cols(times_in_file,11,12) ;Strip out the hour
  28.     times =(/hour+"/"+day/)
  29.     printVarSummary(times)
  30.     print(times)
  31.     nt = dimsizes(times) ; number of times in the file
  32.    
  33.     ;__Available in version 6.4.0 and later________________________________________________
  34.     ; instring = "time"
  35.     ; format   = "%D_%H"
  36.     ; times = cd_inv_string(instring, format)
  37.     ;______________________________________________________________________________________

  38.     mslp = new(nt,float)   
  39.     mwnd = new(nt,float)

  40.     do n=0,nt-1

  41.         slp = wrf_user_getvar(a1,"slp",n)
  42.         mslp(n) = min(slp)

  43.         u10 = wrf_user_getvar(a1,"U10",n)    ; u at 10 m, mass point
  44.         v10 = wrf_user_getvar(a1,"V10",n)    ; v at 10 m, mass point
  45.         wnd = sqrt(u10^2+v10^2)/0.5144 ;compute and convert to 'knot' units
  46.         mwnd(n) = max(wnd)

  47.     end do
  48.    
  49.     istrt =ifile*nt
  50.     iend =(ifile+1)*nt-1

  51.     allt(istrt:iend) =times(:)
  52.     minslp(istrt:iend) =mslp(:)
  53.     maxwnd(istrt:iend) =mwnd(:)

  54. end do
  55. printVarSummary(allt)
  56. print(allt)
  57. printVarSummary(minslp)
  58. printVarSummary(maxwnd)

  59. ;***************plot*****************************
  60.    wks      = gsn_open_wks("png","Minimum_SLP")       ; send graphics to PNG file

  61.    res                   = True
  62.    ;---Titles
  63.    res@tiMainString           = "Minimum SLP of Katrina"
  64.    res@tiMainFontHeightF      = 0.02
  65.    res@tiYAxisFontHeightF     = 0.02
  66.    res@tiXAxisString          = "Time/Date (UTC)"
  67.    res@tiYAxisString          = "Minimum SLP (hPa)"

  68.    ;---TickMark
  69.    res@tmXBMode       = "Explicit"    ; Define own tick mark labels. "Automatic","Mannual","Explicit"
  70.    res@tmXBValues     = x(0:39:8)        ; location of explicit labels(x轴实际变量). "Explicit"
  71.    res@tmXBLabels     = allt(0:39:8)     ; labels are the locations(x轴显示的坐标). "Explicit"
  72.    res@tmLabelAutoStride = True       ; control the step automatically
  73.    res@tmXTOn         = False         ; turn off the top tick marks
  74.    res@xyLineThicknesses = 2          ; increase line thickness
  75.    res@xyLineColor    =  "blue"       ; set line color

  76.    plot  = gsn_csm_xy (wks,x(:),minslp(:),res) ; create plot         
  77.    delete([/wks,res/])
  78. ;__________________________________________________
  79.    wks      = gsn_open_wks("png","Maximum_surfacewnd")       ; send graphics to PNG file

  80.    res                   = True
  81.    ;---Titles
  82.    res@tiMainString           = "Maximum Surface Wind of Katrina"
  83.    res@tiMainFontHeightF      = 0.02
  84.    res@tiYAxisFontHeightF     = 0.02
  85.    res@tiXAxisString          = "Time/Date (UTC)"
  86.    res@tiYAxisString          = "Minimum SLP (hPa)"

  87.    res@tmXBMode       = "Explicit"    ; Define own tick mark labels.
  88.    res@tmXBValues     = x(0:39:8)        ; location of explicit labels(x轴实际变量)
  89.    res@tmXBLabels     = allt(0:39:8)     ; labels are the locations(x轴显示的坐标)
  90.    res@tmXTOn         = False         ; turn off the top tick marks
  91.    res@xyLineThicknesses = 2          ; increase line thickness
  92.    res@xyLineColor    =  "red"       ; set line color

  93.    plot  = gsn_csm_xy (wks,x(:),maxwnd(:),res) ; create plot
  94.               
  95.   
  96. end
复制代码
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-7-18 10:06:39 | 显示全部楼层
本帖最后由 Wanglethe 于 2019-7-18 10:08 编辑

其实我这个问题关键在于:为什么wrf_user_getvar从移动网格提取出的slp场在某些时刻出现很大的震荡误差?
1.png
E:\mywork\Katrina\1.png
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

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

谢谢您的建议,还想请问您一下,“wrf中提取台风中心最低气压是根据最低气压以及梯度来计算的”这个是从哪里看到的呢?WRF用户手册里似乎没有提到。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-7-18 16:55:39 | 显示全部楼层
大神麻烦问一下,可以把你的县边界shp文件给我分享一下吗,特别需要麻烦了!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-7-18 18:51:33 | 显示全部楼层
landachao 发表于 2019-7-18 16:55
大神麻烦问一下,可以把你的县边界shp文件给我分享一下吗,特别需要麻烦了!

shp文件是啥?我不知道……
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

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

本版积分规则

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

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

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