200hPa辐散风，速度势及散度脚本

vicky510475733

画的是2012年7月23日06时200hPa辐散风、速度势，以及散度
先选择适当时次，再利用公式进行计算相应分量，不过在计算散度时，数值太小我放大了百万倍，且由于经纬度是自选的，dv = uv2dv_cfd(u,v,lat,lon,2)中的lat和lon不能是数字或者数组，所以先定义了一个数组，存放经纬度信息
lat=（ /5,7.5,10,12.5,15,17.5,20,22.5,25,27.5,30,32.5,35/)
lon = (/5,7.5,10,12.5,15,17.5,20,22.5,25,27.5,30,32.5,35/)
最后计算散度dv = uv2dv_cfd(u,v,lat,lon,2)*1000000
其次选择时次方面我用了一个很笨的方法，就是print日期信息再根据顺序数核对公历日期，如果有更简便的方法请告诉我
;  date = a->time
;  date1 = cd_calendar(date, 0)
;  date2 = cd_calendar(date1, 0)


以下是脚本：
********************************
散度
********************************
begin
a = addfile("C:\Users\Vicky\Desktop\u2012.nc", "r")
b = addfile("C:\Users\Vicky\Desktop\v2012.nc", "r")
u = a->uwnd(817,9,{5:35},{100:130})
v = b->vwnd(817,9,{5:35},{100:130})  
; printVarSummary(u)
wks = gsn_open_wks("pdf","divergence")                  ; send graphics to PNG file
lat = (/5,7.5,10,12.5,15,17.5,20,22.5,25,27.5,30,32.5,35/)
lon = (/5,7.5,10,12.5,15,17.5,20,22.5,25,27.5,30,32.5,35/)
dv = uv2dv_cfd(u,v,lat,lon,2)*1000000
dv@long_name = "divergence*1000000"
dv!0 = "lat"
dv!1 = "lon"
lat@units = "degrees_E"
lon@units = "degrees_N"
printVarSummary(dv)
  res = True
  res@gsnAddCyclic            = False
  res@mpMaxLatF               = 35
  res@mpMinLatF               = 5
  res@mpMaxLonF               = 130
  res@mpMinLonF               = 100
  res@mpOutlineOn             = True   

  res@tiMainString            = "Divergence" ; title
  res@gsnRightString      = "div*1000000/s"
  res@gsnLeftString      = "200hPa"
  plot = gsn_csm_contour_map(wks,dv,res)         ; create a default plot
  end


********************************
辐散风以及速度势
********************************
begin
a = addfile("C:\Users\Vicky\Desktop\u2012.nc", "r")
b = addfile("C:\Users\Vicky\Desktop\v2012.nc", "r")
u = a->uwnd(817,9,{5:35},{100:130})
v = b->vwnd(817,9,{5:35},{100:130})  
; ***************************************************
;先将日历转换为公历，选择所需时次并核对
; ***************************************************
; date = a->time
;  date1 = cd_calendar(date, 0)
; date2 = cd_calendar(date1, 0)
    ; print(date)
    ; print(date2)
    ; printVarSummary(u)
; ***************************************************
; step one 200hPa 原始风场图
; ***************************************************
; wks = gsn_open_wks("pdf","wind")                 ; send graphics to PNG file

  ; res                         = True               ; plot mods desired
  ; res@mpMaxLatF               = 35
  ; res@mpMinLatF               = 5
  ; res@mpMaxLonF               = 130
  ; res@mpMinLonF               = 100
  ; res@mpOutlineOn             = True               ; turn on map outline

  ; res@gsnAddCyclic            =False

  ; res@tiMainString            = "Wind Barbs" ; title
  ; res@gsnRightString      = "m/s"
  ; res@gsnLeftString      = "200hPa"
  ; plot=gsn_csm_streamline_map(wks,u,v,res) ;gsn_csm_vector_map(wks,u,v,res)             ; create plot
; ***************************************************
; step two 200hPa 辐散风及CHI图
; ***************************************************
; 先计算辐散风，再将辐散风u,v分量分开
;*************************************************
  div = uv2dvG_Wrap(u,v)           ; u,v ==> divergence  

  ud    = new ( dimsizes(u), typeof(u), u@_FillValue )
  vd    = new ( dimsizes(v), typeof(v), v@_FillValue )

  dv2uvg(div,ud,vd) ; div  ==> divergent  wind components
  copy_VarCoords(u, ud )
  copy_VarCoords(u, vd )
  ud@long_name  = "Zonal Divergent Wind"
  ud@units      = u@units
  vd@long_name  = "Meridional Divergent Wind"
  vd@units      = v@units
;*************************************************
; calculate 速度势CHI
;*************************************************
  chi = ilapsG_Wrap ( div , 0)
  chi = (/chi/1e6/)
  chi@long_name = "velocity potential"
  chi@units     = "m/s"
;*************************************************
; plot results
;*************************************************   
  wks  = gsn_open_wks("pdf","201272306dvwind&CHI")                ; send graphics to PNG file

  res                      = True
  res@gsnAddCyclic            = False
  res@mpMaxLatF               = 35
  res@mpMinLatF               = 5
  res@mpMaxLonF               = 130
  res@mpMinLonF               = 100
  res@mpOutlineOn             = True   

  res@vcRefMagnitudeF         = 3.                 ; make vectors larger
  res@vcRefLengthF            = 0.050              ; ref vector length
  res@vcGlyphStyle            = "CurlyVector"      ; turn on curly vectors
  res@vcMinDistanceF          = 0.012              ; thin the vectors
  ; res@vcRefAnnoOrthogonalPosF = -1.0               ; move ref vector up

  res@gsnScalarContour     = True               ; vectors over contours
  res@gsnDraw              = False              ; don't draw
  res@gsnFrame             = False              ; don't advance frame

  res@cnLineLabelPerimOn         = True           ; line label box on
  res@cnLineLabelFontHeightF     = .013           ; label font height
  res@cnInfoLabelOn              = False          ; no contour info label

  res@tiMainString         = "Velocity Potential" ; title
  res@gsnRightString      = "VP scaled by 1e6"  ; center string
  res@gsnLeftString        = "Divergent Wind"     ; left string
  res@gsnRightStringFontHeightF   = 0.02
  res@gsnLeftStringFontHeightF   = 0.02


  res@cnHighLabelsOn             = True           ; high labels on
  res@cnHighLabelBackgroundColor = "white"        ; background white
  res@cnHighLabelPerimOn         = True           ; box on

  res@cnLowLabelsOn             = True            ; low label on
  res@cnLowLabelBackgroundColor = "white"         ; background white
  res@cnLowLabelPerimOn         = True            ; box on

  plot=gsn_csm_vector_scalar_map(wks,ud(:,:),vd(:,:),chi(:,:),res)
  plot = ZeroNegDashLineContour(plot)             ; neg lines dashed
                          ; zero line double thick
  draw(plot)
  frame(wks)
end

独一吖

楼主你这数据是均匀格点的，用的计算函数是高斯格点的。

sunflower

好像楼主分享

白衣风月倾天下

感谢分享 马一下

fengyy789

感谢楼主分享

luyunlong

想请问一下200hpa风的数据在哪下载的

vicky510475733

独一吖 发表于 2018-12-10 18:26
楼主你这数据是均匀格点的，用的计算函数是高斯格点的。

有什么不一样么？？

不会游泳的fish

楼主，请问一下ncl计算辐散风的函数看要求好像说要全球的格点数据，如果只有区域的求会有问题么？

1297712690

不会游泳的fish 发表于 2020-12-23 21:32
楼主，请问一下ncl计算辐散风的函数看要求好像说要全球的格点数据，如果只有区域的求会有问题么？

区域应该不能球谐函数展开吧？要考虑边界条件

1297712690

不会游泳的fish 发表于 2020-12-23 21:32
楼主，请问一下ncl计算辐散风的函数看要求好像说要全球的格点数据，如果只有区域的求会有问题么？

楼主这个感觉更像是区域的