分享一个画湿位涡的脚本

1649518749

ncl画湿位涡官网没有直接给出，官网只有画位涡的，论坛上也找不到。
脚本是本人根据公式编的，如果有错误，欢迎大家指出，好做修改。
;   Example script to produce Vorticity plots from WRF ARW model data
;   Novemner 2008

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
  a = addfile("/home/fc/datafile/wrfoutdata/wrfout_d01_2016-08-17_18:00:00","r")

; We generate plots, but what kind do we prefer?
;  type = "x11"
; type = "pdf"
  type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"mpv")
;gsn_define_colormap(wks,"rainbow")
; Set some basic resources
  res = True
  res@MainTitle                   = "REAL-TIME WRF"

  pltres = True
  mpres = True
;>==============set  for   map====================================<
        m=asciiread ("/home/fc/datafile/barmes.txt", (/25,29,3/) , "float")      
        lat=fspan(10,70,25)
        lon=fspan(60,130,29)
        lat@units="degrees_north"
        lon@units="degrees_east"  
        m!0="lat"
        m!1="lon"
        m&lat=lat
        m&lon=lon

         mpres@mpLimitMode       = "LatLon"
         mpres@mpMinLatF         = 15         
         mpres@mpMaxLatF         = 55
         mpres@mpMinLonF         = 72
         mpres@mpMaxLonF         = 136

        mpres@mpDataBaseVersion="MediumRes"    ;MediumRes  ;;; Ncarg4_1            
        mpres@mpDataSetName="Earth..4"                        
        mpres@mpOutlineOn            = True
        mpres@mpOutlineSpecifiers=(/"China:states","Taiwan","Disputed area between India and China","India:Arunachal Pradesh"/)      
        mpres@mpOutlineBoundarySets ="NoBoundaries"   
        mpres@mpAreaMaskingOn = True                          
        mpres@mpOceanFillColor = "white"                           
        mpres@mpLandFillColor         = "white"
        mpres@mpInlandWaterFillColor = "white"  
        mpres@mpNationalLineColor     = "Black"  ;guojie yan se
        mpres@mpUSStateLineColor          = "Black" ;shengjie yanse
        mpres@mpGeophysicalLineColor  = "Black"  ; kao shui diyu yan se
        mpres@mpNationalLineThicknessF    = 1.5
        mpres@mpGeophysicalLineThicknessF = 1.5
        mpres@mpUSStateLineThicknessF     = 1.5
;>============================================================<


; What times and how many time steps are in the data set?
  times  = wrf_user_getvar(a,"times",-1) ; get times in the file
  ntimes = dimsizes(times)         ; number of times in the file

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  do it = 0,ntimes-1,2             ; TIME LOOP

    print("Working on time: " + times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots


  ; Get the data
    avo   = wrf_user_getvar(a,"avo",it)               ;absolute vorticity
    eth   = wrf_user_getvar(a,"eth",it)                 ;Equivalent Potential Temperature；官网上只有相当位温，不知是不是假相当位温。
    p     = wrf_user_getvar(a,"pressure",it)
     ua     = wrf_user_getvar(a,"ua",it)
   va   = wrf_user_getvar(a,"va",it)

  ; Interpolate to pressure
    eth925_plane = wrf_user_intrp3d(eth,p,"h",925.,0,False)
    ua925_plane = wrf_user_intrp3d(ua,p,"h",925.,0,False)
   va925_plane = wrf_user_intrp3d(va,p,"h",925.,0,False)

    eth850_plane = wrf_user_intrp3d(eth,p,"h",850.,0,False)
    ua850_plane = wrf_user_intrp3d(ua,p,"h",850.,0,False)
   va850_plane = wrf_user_intrp3d(va,p,"h",850.,0,False)

    R=6370949.0
    pi=3.14159
    g=9.8
    deth=eth925_plane-eth850_plane
    dp=(925-850)*100
    du=ua925_plane-ua850_plane
    dv=va925_plane-va850_plane
    dy=2*R*pi/180.
do lati=15., 55.
    dx=2.0*R*cos(lati*pi/180.0)*pi/180.     
    f=2*7.292*sin(lati*3.14159/180)*0.00001
    mpv1=-g*((dv/dx-du/dy)+f)*deth/dp
    mpv2=g*((dv/dp)*(deth/dx)-(du/dp)*(deth/dy))
    mpv=mpv1+mpv2
    mpv=mpv*10^6
end do
    ; Plotting options
      opts = res                        
      opts@cnFillOn = True  
      opts@lbTitleString = "MPV"
       opts@FieldTitle="850hPa MPV"
      opts@lbTitleString = "m^2*k/s*kg"     
     ; opts@gsnSpreadColorEnd = -3  ; End third from the last color in color map
      opts@ContourParameters = (/ -4., 2., 1./)
      opts@cnFillColors = (/"navyblue","blue4","blue", "steelblue2","cadetblue1","red","red3","red4"/)
      contour = wrf_contour(a,wks,mpv,opts)
      delete(opts)

    ; MAKE PLOTS                                       
    ;  plot = wrf_map_overlays(a,wks,(/contour_a/),pltres,mpres)
      plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  end do        ; END OF TIME LOOP
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
end

ayzqs

谢谢分享，试了一下，有些不解，模拟区域外还有数据

1649518749

ayzqs 发表于 2016-8-26 17:31
谢谢分享，试了一下，有些不解，模拟区域外还有数据

啊？谢谢关注，麻烦你上传张图片看看呗。看看哪的问题。

andrewsoong

楼主的分享精神值得赞扬，原创精神亦是！

1649518749

andrewsoong 发表于 2016-8-26 17:40
楼主的分享精神值得赞扬，原创精神亦是！

我都不知道对不对。

andrewsoong

1649518749 发表于 2016-8-26 17:43
我都不知道对不对。

的确是区域外有数据，我觉得可能每个人的lat不一样导致的吧

1649518749

andrewsoong 发表于 2016-8-26 17:44
的确是区域外有数据，我觉得可能每个人的lat不一样导致的吧

你们要把范围改了呀、opts@ContourParameters = (/ -4., 2., 1./)这是我调过之后的

1649518749

andrewsoong 发表于 2016-8-26 17:44
的确是区域外有数据，我觉得可能每个人的lat不一样导致的吧

那你们画出来的图长啥样啊

andrewsoong

1649518749 发表于 2016-8-26 18:02
那你们画出来的图长啥样啊

完全超出了模拟区域，建议楼主检查脚本

1649518749

andrewsoong 发表于 2016-8-26 18:03
完全超出了模拟区域，建议楼主检查脚本

我没明白你们所说的完全超出模拟区域是什么情况，能上张图吗，我看看，再来修改