NCL画假相当位温图

天玖Y

本人用fnl资料计算冬季某一时间段内，200hPa-1000hPa的假相当位温，参照了论坛里很多的脚本，包括grads和ncl的，但是计算出来的值一直不对；同时也参考了一些冰面是、水面上计算的方法，改了很多版本，就是找不到问题出在哪里，请大家帮忙看看

以下是我的脚本：
begin

fi = addfile("./new01.grib2", "r")
time  = fi->initial_time0_hours
;
YYYYMMDDHH = tostring(cd_calendar(time, -3))
dims = dimsizes(YYYYMMDDHH)
cst=new((dims),string)
cst1=new((dims),string)
do i = 0,dims-1
utc = str_get_cols(YYYYMMDDHH(i), 0, 7)
hr = str_get_cols(YYYYMMDDHH(i), 8, 9)
cst(i) = systemfunc("date +'%Y%m%d%H' -d " + "'" + utc + " " + hr + " " + 8 + " hour'")
cst1(i) = str_get_cols(cst(i), 2, 9)
end do


;读取数据
rh=fi->RH_P0_L100_GLL0(:,12:30,{20:40},{113})
;print(rh(0,0,0))   ;rh=17.2%
lev=fi->lv_ISBL0
levv=lev/100
;print(levv)
dims_lev = dimsizes(levv)
;print(levv)
t1=fi->TMP_P0_L100_GLL0(:,12:30,{20:40},{113})
tt=t1-273.16
copy_VarCoords(t1,tt)

eqt=new((/20,19,21/),float)
do k = 0,dims-1
        do n = 0,dims_lev-13
                do i= 0,20
                                   ;es = 6.112*exp((17.67*tt(k,n,i))/(tt(k,n,i)+273.16-29.65))
                                es = where( tt(k,n,i).lt.0, 6.1078*exp(21.87* tt(k,n,i)/( t1(k,n,i)-7.66)),\
                                   6.1078*exp(17.2693882*tt(k,n,i)/(t1(k,n,i)-35.86)))
                                qs = 0.62197*es/(levv(n)-0.378*es)  
                                q  = qs*rh(k,n,i)/100.
                                ;print(q)
                                ;e=prs*q/(0.62197+0.378*q)+1e-10       ;1e-10=1×10^(-10)=1/10000000000
                                ee = (levv(n)*q)/(0.62197+q)   
                                ;ee = (levv(n)*q)/(0.62197+q)+10^(-10)
                                tlcl=55.0+2840.0/(3.5*log(t1(k,n,i))-log(ee)-4.805)
                                theta=t1(k,n,i)*((1000./levv(n))^(0.2854*(1.0-0.28*q)))
                                eqt(k,n,i)=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))
                end do
        end do
end do

end



图：

随缘

且不看别的，首先你让ncl做这样的三层嵌套循环，计算效率就慢了好多

许市

随缘 发表于 2019-2-21 07:49
且不看别的，首先你让ncl做这样的三层嵌套循环，计算效率就慢了好多

的确循环效率太低，但她这个为甚出图错误？我用了其他方法出图为正确的了

许市

; 假相当位温

begin
fi = addfile("./new.grib2", "r")
time  = fi->initial_time0_hours
;
YYYYMMDDHH = tostring(cd_calendar(time, -3))
dims = dimsizes(YYYYMMDDHH)
cst=new((dims),string)
cst1=new((dims),string)
do i = 0,dims-1
    utc = str_get_cols(YYYYMMDDHH(i), 0, 7)
    hr = str_get_cols(YYYYMMDDHH(i), 8, 9)
    cst(i) = systemfunc("date +'%Y%m%d%H' -d " + "'" + utc + " " + hr + " " + 8 + " hour'")
    cst1(i) = str_get_cols(cst(i), 2, 9)
end do
;读取数据
rh=fi->RH_P0_L100_GLL0(:,10:30,:,:)
;print(rh(0,0,0))   ;rh=17.2%
lev=fi->lv_ISBL0(10:30)
levv=lev/100
;print(levv)
dims_lev = dimsizes(levv)
;print(levv)
tk=fi->TMP_P0_L100_GLL0(:,10:30,:,:)
tc=tk-273.16
copy_VarCoords(tk,tc)
printVarSummary(tc)
levv4d=conform(tc,levv,1)
;print (levv4d(0,20,0:10,0:10))
copy_VarCoords(tk,levv4d)
printVarSummary(levv4d)
eqt=new((/dims,dims_lev,181,360/),float)
copy_VarCoords(tk,eqt)
;es = 6.112*exp((17.67*tt(k,n,:,:))/(tt(k,n,:,:)+273.16-29.65))
;es = where( tt(k,n,:,:).lt.0, 6.1078*exp(21.87* tt(k,n,:,:)/( t1(k,n,:,:)-7.66)),\
; 6.1078*exp(17.2693882*tt(k,n,:,:)/(t1(k,n,:,:)-35.86)))
;qs = 0.62197*es/(levv(n)-0.378*es)    ;prs 为高度，设置高度循环
;q  = qs*rh(k,n,:,:)/100.
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, levv4d-0.278*exp(17.26*(tk-273.16)/(tk-35.86)),\
levv4d-0.278*exp(21.87*(tk-273.16)/(tk-7.66)))
q=a1/b1/100
ee = (levv4d*q)/(0.62197+q)   
tlcl=55.0+2840.0/(3.5*log(tk)-log(ee)-4.805)
theta=tk*((1000./levv4d)^(0.2854*(1.0-0.28*q)))
eqt=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))

;eqt!0="time"
;eqt!1="levels"
;eqt!2="lat"
;eqt!3="lon"
;copy_VarCoords(tk,eqt)



res      = True
res@gsnDraw        = False          ; do not draw
res@gsnFrame       = False          ; do not advance the frame
res@trYReverse   = True     ; Reverse the Y values.
res@cnFillOn          =False               ; turns on color fill
res@cnFillPalette     = "BlueDarkRed18"             ; set color map
res@cnLineLabelsOn    = True               ; no contour labels
;res@cnLevelSelectionMode = "ManualLevels"  ; manual levels
;res@cnMinLevelValF       = -1.2
;res@cnMaxLevelValF       =  1.2
;res@cnLevelSpacingF      =  0.1


; x-bottom axis changes
;res@tmXBMode       = "Explicit" ; Define own tick mark labels.
;res@tmXBValues     = eqt&lat      ; location of explicit labels
;res@tmXBLabels     = cst1



;res@tmXBLabelAngleF      = 45.            ; change label angle
;res@tmXBMinorOn          = False          ; No minor tick marks.
;res@tmXBLabelJust        = "CenterCenter" ; label justification
;res@tmXBLabelFontHeightF = .015           ; Font size




;******************画图****************************
    wks = gsn_open_wks("png", "tj")

    res@gsnLeftString = "tjj"

    plot  = gsn_csm_pres_hgt(wks,eqt(0,:,{20:40},{113}),res)

getvalues plot                        ; get the X/Y axis min/max for use in the loglin plot
    "trXMinF"    : trxmin
    "trXMaxF"    : trxmax
    "trYMinF"    : trymin
    "trYMaxF"    : trymax
  end getvalues

  loglin = create "logling" logLinPlotClass wks     ; draw a loglin plot, with expanded X/Y axis
    "trXMinF"        : trxmin-1
    "trXMaxF"        : trxmax+1
    "trYMinF"        : trymin
    "trYMaxF"        : trymax
    "trYReverse"     : True
    "vpXF"           : .15            ; set the X-axis NDC starting point
    "vpYF"           : .8             ; set the Y-axis NDC starting point
    "vpWidthF"       : .7             ; set the width of the plot in NDC units
    "vpHeightF"      : .45            ; set the height of the plot in NDC units
   
  end create
  overlay(loglin,plot)              ; overlay plot with the loglin plot
  draw(loglin)                      ; draw the plot
  frame(wks)                  
end
http://bbs.06climate.com/forum.php?mod=attachment&aid=ODExNTd8YjZkNmNkYTk4ZDkyMDJjODliM2YzNGZkNjU1YjdjMTh8MTczMjQzMDY1NA%3D%3D&request=yes&_f=.png

伽蓝鸟

NCL有函数计算相当位温吧，那个跟假相当位温在数值上非常非常接近了

许市

伽蓝鸟 发表于 2019-2-27 10:12
NCL有函数计算相当位温吧，那个跟假相当位温在数值上非常非常接近了

哪个函数?

伽蓝鸟

许市 发表于 2019-3-16 10:09
哪个函数?

你上官网查查看？我自己不用NCL。。之前用MATLAB写了个

大仓鼠

许市 发表于 2019-2-21 11:51
; 假相当位温

begin

请问脚本第三行time  = fi->initial_time0_hours表示什么，运行的时候显示我的grib2数据里没initial_time0_hours这一变量，我打开数据说明倒是每个变量的信息最后一行都有initial_time这一项，但是没法读取呀？
还有一个问题就是这是话哪一个经度上的剖面图呢,在哪里调整呢？

yuyumage

是不是pot_vort_isobaric？？？？

天玖Y

大仓鼠 发表于 2019-4-14 23:43
请问脚本第三行time  = fi->initial_time0_hours表示什么，运行的时候显示我的grib2数据里没initial_time ...

在最后画图的plot里面设置经度就行