2020-04-17

利用FNL资料画Tlog-p图，参考链接http://www.ncl.ucar.edu/Applications/skewt.shtml
代码如下：

;**************************************************
; skewt_6.ncl
;
; Concepts illustrated:
;   - Reading RUC (Rapid Update Cycle)  GRIB data
;   - Using getind_latlon2d to determine grid locations
;   - Drawing Skew-T plots at nearest grid locations  
;                    to user specified locations
;**************************************************
;
; These files are loaded by default in NCL V6.2.0 and newer
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;
; This file still has to be loaded manually
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
;***********************************************
; RUC Data downloaded from:
;     http://nomads.ncdc.noaa.gov/data.php?name=access#hires_weather_datasets
; RUC-ANL: http
;     Specifically: http://nomads.ncdc.noaa.gov/data/rucanl/201205/20120501/
;***********************************************
; The GRIB file's contents can be examined via:
;     ncl_filedump -itime ruc2anl_130_20120501_0000_000.grb2  | less
;***********************************************
 begin

; --- Read RUC GRIB file------------;
f = addfile("new.grib2", "r")
;printVarSummary(f)
time  = f->initial_time0_hours
YYYYMMDDHH = tostring(cd_calendar(time, -3))
dims = dimsizes(YYYYMMDDHH)
cst=new((dims),string)
cst1=new((dims),string)
do i = 0,dims-1
    utc = str_get_cols(YYYYMMDDHH(i), 0, 7)
   hr = str_get_cols(YYYYMMDDHH(i), 8, 9)
    cst(i) = systemfunc("date +'%Y%m%d%H' -d " + "'" + utc + " " + hr + " " + 8 + " hour'")
    cst1(i) = str_get_cols(cst(i), 2, 9)
end do
;  p     = f->lv_ISBL0               ;  ( lv_ISBL0 )
  p     = f->lv_ISBL5({10000:100000})
  print(p)
  ;time  = f->initial_time0_encoded  ; yyyymmddhh.hh_frac

                                    ; RUC grid point locations
  lat2d = f->lat_0              ; ( ygrid_0, xgrid_0 ) 
  lon2d = f->lon_0
;  xCopy = new ( (/dimsizes(lat)-1,dimsizes(lon)-1/), typeof(x))
;  do i=0,dimsizes(lat)-1
;;    do j=0,dimsizes(lon)-1
;;     lat2d(k,:,:) =lat(i)
;     lat2d(k,:,:) = lat()
;  end do
  print("lat2d: min="+min(lat2d)+"  ;  max="+max(lat2d))
  print("lon2d: min="+min(lon2d)+"  ;  max="+max(lon2d))
  printVarSummary(lat2d)
  p     = p*0.01                    ; change units
  p@units = "hPa"                   ; skewT, mixhum_ptrh use mb (hPa)
;--- Specify one or more locations
  printVarSummary(p)

nn=29.28
mm=113.5

;*************************
; create plot(s)
;*************************
do i = 0, dims-1
    skewtOpts                 = True
    skewtOpts@DrawColAreaFill = True    ; default is False
   
    dataOpts   = True
    dataOpts@PrintZ = True                          
    tk   = f->TMP_P0_L100_GLL0(i,{10000:100000},{nn},{mm})
    z    = f->HGT_P0_L100_GLL0(i,{10000:100000},{nn},{mm})
    rh   = f->RH_P0_L100_GLL0(i,{10000:100000},{nn},{mm})
    u    = f->UGRD_P0_L100_GLL0(i,{10000:100000},{nn},{mm})
    v    = f->VGRD_P0_L100_GLL0(i,{10000:100000},{nn},{mm})
    printVarSummary(tk)
    printVarSummary(rh)
    printVarSummary(u)
    printVarSummary(z)

         ; change units and calculate needed variables
    tc   = tk - 273.15
    tc@units= "degC"            

    q    = mixhum_ptrh (p, tk, rh, 2)
    q@units = "kg/kg"

    tdc  = dewtemp_trh(tk,rh) - 273.15 
    tdc@units = "degC"            

    wspd = sqrt((u*2.5)^2 + (v*2.5)^2)
    wdir = wind_direction(u,v,0)

     ;itime= toint(time)
    skewtOpts@tiMainString    = cst(i)+"("+nn+","+mm+")"
    skewtOpts@DrawFahrenheit  = False   ; default is True
          ; each location will have a different file name      
    wks  = gsn_open_wks ("png", "FNL_skewt_"+cst(i))
    skewt_bkgd = skewT_BackGround (wks, skewtOpts)
    skewt_data = skewT_PlotData   (wks, skewt_bkgd, p,tc,tdc,z \
                                       , wspd,wdir, dataOpts)
    draw (skewt_bkgd)
    draw (skewt_data)
    frame(wks)
end do
end  
;shox 肖沃特稳定指数 Showalter stability index 斜坡对流 sloping convection

FNL_skewt_2020032120.png

（21日夜间湘北狂风骤雨，多地出现冰雹天气，其中临湘还出现了超级单体）
值得注意的是FNL资料一直在更新，2020更新了lv_ISBL5垂直分布，共有34层，这样Tk，z都是34层，而rh和 u、v是31层，比lv_ISBL0（31层）多了3层，在画图时要注意挑选层次{10000:100000}

最后编辑于：2020.04.17 21:24:48