- 积分
- 1095
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-4-3
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 SG哭晕在厕所 于 2017-10-21 08:08 编辑
最近要画一张某区域内沙氏指数的水平分布图,
参照坛子里某个gs脚本写了个NCL脚本,用的是wrfout数据,参考吕新刚和周志强的一篇文章更正了一下公式,但是出图后数值明显不对,大部分区域都在-10以下了,自己也对着文献思考了很久,大约知道问题可能出在循环里,但是还是不知道怎么解决,只能上论坛来求助大神们了,不知道设置个回帖奖励能不能防沉。。。。
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
a = addfile("I:\RRTMG\wrfout_d02_2016-11-03_18_00_00","r")
type = "png"
wks = gsn_open_wks(type,"E:/si")
; times= wrf_user_list_times(a); get times in the file
; ntimes = dimsizes(times) ; number of times in the file
tc = wrf_user_getvar(a,"tc",-1) ; T in C
u= wrf_user_getvar(a,"ua",-1) ; u averaged to mass points
v= wrf_user_getvar(a,"va",-1) ; v averaged to mass points
p= wrf_user_getvar(a, "pressure",-1) ;
rh = wrf_user_getvar(a,"rh",-1) ; relative humidity;
z = wrf_user_getvar(a, "z",-1) ; grid point height
pressure_levels = (/ 925.,850.,700.,500.,200./) ; pressure levels to plot
nlevels = dimsizes(pressure_levels) ; number of 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)
u_plane= wrf_user_intrp3d( u,p,"h",pressure_levels,0.,False)
v_plane= wrf_user_intrp3d( v,p,"h",pressure_levels,0.,False)
p_plane= wrf_user_intrp3d( p,p,"h",pressure_levels,0.,False)
p1=z_plane(0,1,:,:)
t1=tc_plane(0,1,:,:)
rh1=rh_plane(0,1,:,:)
td1=t1-(14.55+0.114*t1)*(1-0.01*rh1) + ((2.5+0.007*t1)*(1-0.01*rh1))^3+ (15.9+0.37*t1)*(1-0.01*rh1)^14
tk1=tc_plane(0,1,:,:)+273.16
est1=6.112*10^((7.5*t1)/(t1+237.3))
etd1=6.112*10^((7.5*td1)/(td1+237.3))
e1=etd1
u1=(etd1/est1)*100
tl1=2840/(3.5*log(tk1)-log(e1)-4.805)+55
r1=622*(e1/(p1-e1))
p2=z_plane(0,3,:,:)
t2=tc_plane(0,3,:,:)
rh2=rh_plane(0,3,:,:)
td2=t2-(14.55+0.114*t2)*(1-0.01*rh2) +((2.5+0.007*t2)*(1-0.01*rh2))^3+ (15.9+0.37*t2)*(1-0.01*rh2)^14
tk2=tc_plane(0,3,:,:)+273.16
est2=6.112*10^((7.5*t2)/(t2+237.3))
etd2=6.112*10^((7.5*td2)/(td2+237.3))
e2=etd2
u2=(etd2/est2)*100
tl2=2840/(3.5*log(tk2)-log(e2)-4.805)+55
r2=622*(e2/(p2-e2))
T3=tl1-273.15 ;绝对温度转换为摄氏温度设为气块初始温度
a1=tk1*(1000/p1)^(0.2854*(1-0.00028*r1))*exp((3376/tl1-2.54)*r1*(1+0.00081*r1)) ;计算850hpa的相当位温(K)
akb=new((/177,162/),float)
tk=new((/177,162/),float)
est=new((/177,162/),float)
etd=new((/177,162/),float)
e=new((/177,162/),float)
us=new((/177,162/),float)
tl=new((/177,162/),float)
r=new((/177,162/),float)
a2=new((/177,162/),float)
akb=t2
printVarSummary(a1)
do i=0,176
do j=0,161
n=1
do while (akb(i,j).gt.0.01)
T3(i,j)=T3(i,j)+0.01*n
tk(i,j)=T3(i,j)+273.15
;print(T3(i,j)+" "+tk(i,j))
est(i,j)=6.112*10^((7.5*T3(i,j))/(T3(i,j)+237.3))
; print(T3(i,j)+" "+est(i,j))
etd(i,j)=6.112*10^((7.5*td2(i,j))/(td2(i,j)+237.3))
e(i,j)=etd(i,j)
us(i,j)=(etd(i,j)/est(i,j))*100
; print(T3(i,j)+" "+us(i,j))
tl(i,j)=2840/(3.5*log(tk(i,j))-log(e(i,j))-4.805)+55
; print(T3(i,j)+" "+tl(i,j))
r(i,j)=622*(e(i,j)/(p2(i,j)-e(i,j)))
a2(i,j)=tk1(i,j)*(1000/p1(i,j))^(0.2854*(1-0.00028*r(i,j)))*exp((3376/tl1(i,j)-2.54)*r(i,j)*(1+0.00081*r(i,j)))
; print(T3(i,j)+" "+a2(i,j))
akb(i,j)=a2(i,j)-a1(i,j)
end do
end do
end do
t=T3
si=t2-t
;print(si)
f = addfile("I:\RRTMG\wrfout_d02_2016-11-03_11_00_00","r")
lat = f->XLAT(0,:,50)
lat@unit = "degree_north"
lon = f->XLONG(0,50,:)
lon@unit = "degree_east"
si!0 = "lat"
si!1 = "lon"
si&lat = lat
si&lon = lon
res = True ; Variable to hold plot options
res@gsnDraw = False
res@gsnFrame = False
res@gsnAddCyclic = False
res@mpOutlineOn = True
res@cnFillOn = True ; Turn on contour fill.
res@cnLinesOn = True
res@cnLineLabelsOn = False
;
res@cnLevelSelectionMode = "ManualLevels" ; manually select levels
; res@cnLevelSpacingF = 1000 ; contour spacing
;
res@cnMinLevelValF = 95000 ; min level
;
res@cnMaxLevelValF = 103000 ; max leve
;
res@lbOrientation = "Vertical" ; Move labelbar to side.
res@gsnMajorLatSpacing = 3
res@gsnMajorLonSpacing = 3
res@mpGeophysicalLineThicknessF = 3.
res@mpProvincialLineThicknessF = 3.
res@pmTickMarkDisplayMode = "Always"
res@mpDataBaseVersion = "MediumRes"
res@mpDataSetName = "Earth..4"
res@mpOutlineSpecifiers = (/"China:states"/)
res@mpLandFillColor = 0
res@lbLabelBarOn = True
res@mpLimitMode = "Latlon"
res@mpMinLatF = 34
res@mpMaxLatF = 43
res@mpMinLonF = 112
res@mpMaxLonF = 122
;
res@cnFillDrawOrder = "PostDraw"
res@mpOceanFillColor = 0
res@mpInlandWaterFillColor = 0
plot = gsn_csm_contour_map(wks, si(:,:), res)
draw(plot)
frame(wks)
end
|
|