- 积分
- 2850
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-10-26
- 最后登录
- 1970-1-1
|
发表于 2016-6-9 14:27:14
|
显示全部楼层
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
f = addfile("../wrfout_d02_2009-09-16_06:00:00.nc","r")
times = wrf_user_list_times(f) ; get times in the filewww.mnmuc.org)
ntimes = dimsizes(times) ; number of times in the file
tc = wrf_user_getvar(f,"tc",-1)
p = wrf_user_getvar(f, "pressure",-1) ;
rh = wrf_user_getvar(f,"rh",-1)
z = wrf_user_getvar(f, "z",-1)
pressure_levels = (/ 1000.,950.,900.,850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,350.,300.,250.,200.,150.,100./)
nlevels = dimsizes(pressure_levels)
tc_plane = wrf_user_intrp3d(tc,p,"h",pressure_levels,0.,False)
z_plane = wrf_user_intrp3d(z,p,"h",pressure_levels,0.,False)
rh_plane = wrf_user_intrp3d(rh,p,"h",pressure_levels,0.,False)
p_plane = wrf_user_intrp3d(p,p,"h",pressure_levels,0.,False)
;t = f->tc
;rh = f->rh
;pre_7 = f->lv_ISBL3 ;(100,1000) 氓<8d><95>盲陆<8d>hPa
;#################################################################################################
;猫庐隆莽庐<97>氓<81><87>莽<9b>赂氓陆<93>盲陆<8d>忙赂?
;==========================
;printVarSummary(rh)
;printVarSummary(t)
tk=tc_plane+273.166
printVarSummary(pressure_levels)
p = conform(tc_plane, pressure_levels, 1)
a1 = where(tk .gt. 263.0, 0.622*6.11*exp(17.26*(tk-273.16)/(tk-35.86)), \
0.622*6.11*exp(21.87*(tk-273.16)/(tk-7.66))) ;
b1= where(tk .gt. 263.0, P-0.278*exp(17.26*(tk-273.16)/(tk-35.86)),\
P-0.278*exp(21.87*(tk-273.16)/(tk-7.66)))
;prs = conform(t,pre_7,0)
;printVarSummary(p)
;es = 6.112*exp(17.67*(t-273.15)/(t-29.65))
;q = rh*(0.62197*es/(prs-es))/100.
q = a1/b1
q1 = q*rh_plane
e = p*q1/100./(0.62197+q1/100.)+1e-10
tlcl = 55.0+2840.0/(3.5* log(tk)-log(e)-4.805)
theta = tk*((1000./p)^(0.2854*(1.0-0.28*q1/100.)))
eqt = theta*exp(((3376./tlcl)-2.54)*q1/100.0*(1.0+0.81*q/100.0))
copy_VarCoords(t,eqt)
printVarSummary(eqt)
; system("rm -f eqttest.nc") ; remove any pre-existing file
; ncdf = addfile("eqttest.nc" ,"c") ; open output netCDF file
; ncdf->eqt = eqt
;##########################################################
;PLOT
;=========================
wks = gsn_open_wks("x11","h_lat_7-21-00-test")
do it =0, ntimes-1
res = True
res@gsnDraw = False ; Do not draw plot
res@gsnFrame = False ; Do not advance frome
res@gsnMaximize = True
res@cnLevelSelectionMode = "ManualLevels" ; manual contouring
res@cnMinLevelValF = 326 ; set min contour level
res@cnMaxLevelValF = 408 ; set max contour level
res@cnLevelSpacingF = 2 ; set contour spacing
res@cnInfoLabelOn = False
res@cnLineThicknessF = 2
res@cnLineLabelDensityF = 0.7
;res@cnSmoothingOn = True
;res@cnSmoothingDistanceF = 0.005
;res@cnSmoothingTensionF = -2.5
plot3 = gsn_csm_pres_hgt(wks,eqt(it,:,{30.5:32.5},{120}),res)
draw(plot3)
frame(wks)
end do
end
请教一下,这个脚本画假相当位温剖面图,总是报错,报错的信息为fatal:Dimension sizes of left hand side and right hand side of assignment do not match
fatal:["Execute.c":8567]:Execute: Error occurred at or near line 31 in file test_theta.ncl,
31行为标红的那行:p = conform(tc_plane, pressure_levels, 1)
请问该如何解决呢,麻烦了啊 |
|