- 积分
- 4422
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2020-5-18
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2022-3-27 19:56:09
|
显示全部楼层
t-logp图吗?
这是脚本,我写的很乱的,你凑合看
;**************************************************
; 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 ... es_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------------;
diri = "D:/1aaa研究生/work3_dblw_LLP/data/20210531-0603_ERA5/"
fili = "download_pre_2021_06_02_06.grib"
; force a 'time' dimension
setfileoption("grb","SingleElementDimensions","Initial_time")
f = addfile(diri+fili,"r")
p = f->lv_ISBL1 ; ( lv_ISBL0 )
p@units = "hPa" ; skewT, mixhum_ptrh use mb (hPa)
time = f->initial_time0 ; yyyymmddhh.hh_frac
; RUC grid point locations
lat2d = f->g0_lat_2 ; ( ygrid_0, xgrid_0 )
lon2d = f->g0_lon_3
print("lat2d: min="+min(lat2d)+" ; max="+max(lat2d))
print("lon2d: min="+min(lon2d)+" ; max="+max(lon2d))
; --- Specify one or more locations
; lat = (/ 25 , 55 /)
; lon = (/110 ,140 /)
; npts = dimsizes(lat)
; lat=42
; lon=121 ;dry
lat=45
lon=122.5 ;rainy
; lat=50
; lon=125
; (47,132)
; lat=47
; lon=132
;*************************
; create plot(s)
;*************************
skewtOpts = True
skewtOpts@DrawColAreaFill = True ; default is False
dataOpts = True
dataOpts@PrintZ = True
; do n=0,npts-1 ; loop over each grid pt
; find grid point nearest the user specified location
; nm = getind_latlon2d (lat2d,lon2d, lat, lon) ;查找最接近用户指定的纬度/经度坐标对的索引(下标)
; nn = nm(0,0)
; mm = nm(0,1)
; print("location=("+lat(n)+","+lon(n)+") grid=("+lat2d(nn,mm)+","+lon2d(nn,mm)+")")
; (40,125)
tk = f->T_GDS0_ISBL(0,:,200,500)
z = f->Z_GDS0_ISBL(0,:,200,500) ;位势高度,需要/9.8得到m
rh = f->R_GDS0_ISBL(0,:,200,500)
u = f->U_GDS0_ISBL(0,:,200,500)
v = f->V_GDS0_ISBL(0,:,200,500)
; (47,132)
; tk = f->T_GDS0_ISBL(0,:,172,528)
; z = f->Z_GDS0_ISBL(0,:,172,528)
; rh = f->R_GDS0_ISBL(0,:,172,528)
; u = f->U_GDS0_ISBL(0,:,172,528)
; v = f->V_GDS0_ISBL(0,:,172,528)
; change units and calculate needed variables
z =z/9.8
z@units="m"
;print(z)
tc = tk - 273.15
tc@units= "degC"
q = mixhum_ptrh (p, tk, rh, 2)
; q =q*1000
q@units = "kg/kg"
tdc = dewtemp_trh(tk,rh) - 273.15
tdc@units = "degC"
wspd = sqrt(u^2 + v^2)
wdir = wind_direction(u,v,0)
itime= toint(time)
skewtOpts@tiMainString = time+" ("+lat+","+lon+")"
; each location will have a different file name
wks = gsn_open_wks ("png", "06_02_06_skewt_("+lat2d(200)+","+lon2d(500)+")")
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 |
|