请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 13167|回复: 15

[作图] NCL画沙氏指数出现问题

[复制链接]

新浪微博达人勋

发表于 2017-10-19 18:55:53 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

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


密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-10-19 18:58:03 | 显示全部楼层
毫无办法中,自顶一下
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-10-20 11:39:08 | 显示全部楼层

回帖奖励 +3 金钱

wrfout的数据有问题吗
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-10-20 12:31:39 | 显示全部楼层
songwei 发表于 2017-10-20 11:39
wrfout的数据有问题吗

没问题,因为已经画过很多物理量了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-10-21 07:13:14 | 显示全部楼层
。。。果然还是要沉
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-10-21 07:24:20 | 显示全部楼层

回帖奖励 +10 金钱

你把图贴出来看一下吧
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-10-21 07:32:17 | 显示全部楼层
沙氏指数,第一次听到这个概念,学习了!!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-10-21 08:07:00 | 显示全部楼层
talkd 发表于 2017-10-21 07:24
你把图贴出来看一下吧

si

si
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-10-22 13:47:45 | 显示全部楼层

回帖奖励 +7 金钱

占位…………
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2017-10-26 12:36:35 | 显示全部楼层
很好的东西学习
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表