用ERA5资料计算假相当位温出错

onenuister

ncl初学者，资料是ERA5 1000-100hPa，在家园上找了相关的脚本，画出来的图在500-550hPa和300hPa附近总是很密集且平行，代码和图片如下，希望有大佬能指点指点我啊（也有想过用pot_temp_equiv_tlcl函数，但没法运行）
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/wrf/WRFUserARW.ncl"

begin
f = addfile("E:/Trend_0.25/2020.01_02_20.nc", "r")

lon = f->longitude ;longitude=281
lat = f->latitude  ;latitude=241
prs = f->level({1000:200})     ;lev=20 100,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,925,950,975,1000
time = f->time     ;time=232
time@units = "hours since 1900-01-01 00:00:00.0"
utc_date = cd_calendar(time, 0)
year = tointeger(utc_date(:,0))
month = tointeger(utc_date(:,1))
day = tointeger(utc_date(:,2))
hour = tointeger(utc_date(:,3))
nmonth = 1
nday_start = 21
nday_end = 25
nhour_start = 0
nhour_end = 18
do iday = nday_start,nday_end
    do ihour = nhour_start,nhour_end,6
        ii = ind(month.eq.nmonth.and.day.eq.iday.and.hour.eq.ihour)
        ;print(ii)
        rh = short2flt(f->r(ii,{1000:200},{20:40},{112}))          ;相对湿度
        tk = short2flt(f->t(ii,{1000:200},{20:40},{112}))          ;温度
        V = short2flt(f->v(ii,{1000:200},{20:40},{112}))           ;径向速度
        omega = short2flt(f->w(ii,{1000:200},{20:40},{112}))       ;垂直速度

        rh@_Fillvalue = -32767
        rh@units = "%"
        tk@_Fillvalue = -32767
        ;tk@units = "K"
        V@_Fillvalue = -32767
        ;V@units = "m s**-1"
        omega@_Fillvalue = -32767
        omega@long_name = " "
        ;omega@units = "Pa s**-1"

        tc = tk - 273.16 ;转换为摄氏度
        copy_VarMeta(tk, tc)
        tc@units = "degC"
        es = new(dimsizes(tk),typeof(tk),tk@_Fillvalue)
        es := where(tk.gt.273.16,6.1078*exp(17.2693882*tc/(tk-35.86)),6.1078*exp(21.8745584*tc/(tk-7.66))) ;单位hPa  ;运用饱和水汽压Tetens经验公式
        copy_VarCoords(tk,es)


        prsconform = conform(es,prs,0) ;将气压一维数组扩展为跟es一样的二维数组
        qs = 0.622*es/(prsconform - 0.378*es) ;计算饱和比湿，单位g/g
        e = es * rh/100  ;计算水汽压，根据rh的单位原因除以100
        copy_VarCoords(rh,e)
        q = 0.622*e/(prsconform - 0.378*e) ;计算比湿
        copy_VarCoords(rh,q)
        tlcl = 55.0+2840.0/(3.5*log(tk)-log(e)-4.805) ;凝结高度的绝对温度
        theta = tk * ((1000/prsconform)^(0.2854*(1.0-0.28*q)))
        copy_VarCoords(rh,theta)
        se = theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))
        ;se = pot_temp_equiv_tlcl(prsconform,tk,tlcl,q*1000,(/0,1,0,1/)
        copy_VarCoords(rh,se)
        ; se@long_name = "potential pseudo-equivalent temperature"
        se@long_name = " "
        ; se@units = "K"
        ;print(theta(0,:))
        wks   = gsn_open_wks ("png", "thse_"+tostring(nmonth)+"_"+tostring(iday)+"_"+tostring(ihour))         
        re = True
        re@gsnDraw = False
        re@gsnFrame = False
        re@gsnLeftString = " "
        re@gsnRightString = " "
        re@cnLinesOn = False
        re@cnFillOn = True
        re@cnInfoLabelOn = False
        re@cnLevelSelectionMode = "ExplicitLevels"
        re@cnLevelSpacingF      = 0.2
        re@tiYAxisString        = "Pressure(hPa)"
        re@tiMainString         = tostring(nmonth)+"_"+tostring(iday)+"_"+tostring(ihour)
        re@cnMaxLevelValF       = 0
        re@cnLevels            = (/-1.8,-1.6,-1.4,-1.2,-1.0,-0.8,-0.6,-0.4,-0.2,0/)
        re@cnFillPalette        = "GMT_gray"
        plot_omega = gsn_csm_pres_hgt(wks, omega, re)

        res = True               
        res@gsnDraw  = False
        res@gsnFrame = False
        res@cnLevelSelectionMode = "ManualLevels"       ;根据上述最值设置等值线数值
        res@cnLevelSpacingF      = 4                  
        ;res@cnMinLevelValF       = 340               
        ;res@cnMaxLevelValF       =  420                 
        res@cnLinesOn            = True
        res@cnSmoothingOn        = True
        res@cnLineLabelsOn       = True                 ;关闭等值线标签
        res@cnFillOn             = False                 ; 打开填色
        res@cnInfoLabelOn        = False
        res@cnLineThicknessF     = 2
        res@cnLineColor          = "red"

        ;res@cnFillPalette        = "BlWhRe"             ; 绘图颜色设置
        prsconform = conform(es,prs*100,0) ;将气压转换为Pa
        omega = omega_to_w(omega,prsconform,tk) ;将垂直速度转换为m/s
        omega = omega*100
        res@vcMapDirection       = False
        res@vcPositionMode       = "ArrowTail"
        res@vcGlyphStyle         = "CurlyVector"
        res@vcLineArrowThicknessF = 3.0
        res@vcLevelSpacingF      = 4
        res@vcRefAnnoOn          = True ;绘制参考箭头
        res@vcRefLengthF         = 0.045 ;参考箭头长度
        ;res@vcRefAnnoBackgroundColor = -1 ;参考箭头背景色为透明
        res@vcRefAnnoString2     = "m/s" ;绘制参考箭头之下的字符串
        ;res@vcRefAnnoOrthogonalPosF = -0.175 ;垂直移动参考箭头
        ;res@vcRefAnnoParallelPosF = 0.1 ;水平移动参考箭头
        res@vcRefAnnoPerimOn     = False ;绘制参考箭头边框
        res@vcRefMagnitudeF      = 40    ;参考箭头所表示物理量大小

        plot_se  = gsn_csm_pres_hgt_vector(wks, se, V, omega, res)
        overlay(plot_omega, plot_se)
        draw(plot_omega)
        frame(wks)

        delete(ii)
        delete(rh)
        delete(tc)
        delete(V)
        delete(omega)
        delete(tk)
        delete(es)
        delete(prsconform)
        delete(qs)
        delete(e)
        delete(q)
        delete(tlcl)
        delete(theta)
        delete(se)
    end do
end do
end

onenuister

破案了破案了，prs = int2flt(f->level({1000:200}) )

Miss_星期八

感谢楼主困扰很久的问题

zaozedi

请问楼主用的是这个资料吗

onenuister

zaozedi 发表于 2021-8-3 11:28
请问楼主用的是这个资料吗

我用的就是这个

Уēｓ、尜玉

请问楼主知道tltc_rh_bolton和pot_temp_equiv_tlcl函数为何调用过程种无法使用么

onenuister

Уēｓ、尜玉 发表于 2021-10-24 23:13
请问楼主知道tltc_rh_bolton和pot_temp_equiv_tlcl函数为何调用过程种无法使用么

可能是ncl版本太低？我的sublime没有把这个函数高亮

小浪

onenuister 发表于 2021-4-8 14:38
破案了破案了，prs = int2flt(f->level({1000:200}) )

theta = tk * ((1000/prsconform)^(0.2854*(1.0-0.28*q)))，因为1000/pre是两个整型除以整型，可以1000.0

1271756664

楼主您好，请问阴影表示的是什么物理量，初学者，还望指教

onenuister

1271756664 发表于 2021-11-11 15:58
楼主您好，请问阴影表示的是什么物理量，初学者，还望指教

垂直速度，hPa/s