求助如何在垂直方向上正确使用cdiff中央差分函数

summersleave

看了很多帖子cdiff都是在水平方向上的应用示例，垂直方向上没找到，因为需要求假相当位温垂直递减率，cdiff(thetase,z)/cdiff(p,z)

下了 1000 925 850 700 500这几层的nc 温度和相对湿度数据 计算假相当位温eqt后
gs 写成
'define prs=lev' ;*获得某层高度的气压
'set gxout contour'
'set clab forced'
'set ccolor 1'
'set clopts -1 8 0.25'
*'set cmax -1'
'set cthick 8'
*'set cint 2'
*'set clab masked'
'd cdiff(eqt,prs)'

收到报错如图
请大神们看看是什么问题
感谢

栤唥de嶶笑

eqt是怎么计算的也没贴上来啊？eqt是否各层都算了？另外，Warning提示有1个值<=0，检查下是否除数出现负值情况？

栤唥de嶶笑

可以试一下加个语句，'set lev 1000 500'试试。还有一个问题，你ctl文件里面，是否只有1000 925 850 700 500这几层？如果不是的话，要把中间的层次的eqt也算出来，例如975 900 750等等这些层次，要跟ctl保持一致

伽蓝鸟

ept是算相当位温的内置函数，你能找到这个内置函数，干嘛不去看看其他内置函数呢
madvw(W,EXPR)
This function calculates

   -d(W*EXPR)/dp
using an upwind scheme.

summersleave

'reinit'

'sdfopen I:\ncdata\se2008.nc'
'enable print  E:\2021\paper\ecdata\thetse700.gmf'
'set vpage 0 11 0 8.5'
'set parea 1.0 10.0 1.0 8.0'
*'set t 21 24'
*'set lat 10 45'
*'set lon 90 135'
'set lev 1000 500'


'set clopts -1 -1 0.14' ;*等值线标值属性 颜色 粗细 大小
'set xlopts 1 4 0.14'   ;*坐标刻度和标值的属性 颜色 粗细 大小
'set ylopts 1 4 0.14'
'set csmooth on'   ;*光滑开关

'define tc=t-273.16' ;*某高度层的摄氏温度C
'define tk=t' ;*某高度层的开氏温度K
'define rh=r' ;*某高度层的相对湿度Relative humidity%
'define prs=lev' ;*获得某层高度的气压


******************************************************************************************
*求饱和水汽压Tetens经验公式
*水面es,tk开氏温度,tc摄氏温度
if(tk>273.16)
'define es=6.1078*exp(17.2693882*tc/(tk-35.86))'
endif
*冰面es,tk开氏温度,tc摄氏温度
if(tk<=273.16)
'define es=6.1078*exp(21.8745584*tc/(tk-35.86))'
endif
*饱和比湿
'define qs=0.622*es/(prs-0.378*es)'
*用相对湿度等求比湿
'q=rh*qs/100'
*水汽压
'e=prs*q/(0.622+q)'
*凝结高度的绝对温度,tk起始面上绝对温度K，
'define tlcl=55.0+2840.0/(3.5*log(tk)-log(e)-4.805)'
*求假相当位温Bolton公式,se为开氏温度K
'define theta=tk*pow((1000/prs),(0.2854*(1.0-0.28*q)))'
'define se=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))'
'define dp=100*(prs(z-1)-prs(z+1))'
'define dse=se(z-1)-se(z+1)'


******************************************************************************************

i=21
while(i<=24)
'set t ' i
'set grads off'
'set grid off'
'set xlopts 1 5 0.25'
'set ylopts 1 5 0.25'
'set xlint 5'
'set ylint 5'
'set map 1 1 3'
'set mpdset cnworld'

'set font 5'
'set csmooth on'
'set cterp on'
*'set parea 1 10.0 1 7.8'
'set lat 15 35'
'set lon 100 125'

'set lev 700'
'set gxout contour'
'set clab forced'
'set ccolor 1'
'set clopts -1 8 0.25'
*'set cmax -1'
'set cthick 8'
*'set cint 2'
*'set clab masked'
'd cdiff(se,z)/cdiff(prs/z)'


'q time'
xx=subwrd(result,3)
'draw title  'xx'(UTC)'
*'run C:\OpenGrADS\Contents\Resources\Scripts\cbarn.gs '
*1.0 1.0 9.8 4.0

'print'
'c'
i=i+1
endwhile
'disable print'
;

summersleave

栤唥de嶶笑 发表于 2021-6-15 11:51
eqt是怎么计算的也没贴上来啊？eqt是否各层都算了？另外，Warning提示有1个值

我把gs贴上来了  这里面是哪里涉及到了log呢 提示底数负数

summersleave

栤唥de嶶笑 发表于 2021-6-15 11:54
可以试一下加个语句，'set lev 1000 500'试试。还有一个问题，你ctl文件里面，是否只有1000 925 850 700 50 ...

我担心是垂直方向层数不均等问题 重新下了1000-500百帕每隔50百帕的数据来算

summersleave

伽蓝鸟 发表于 2021-6-15 12:25
ept是算相当位温的内置函数，你能找到这个内置函数，干嘛不去看看其他内置函数呢
madvw(W,EXPR)
This fun ...

假相当位温是我参考版里的帖子求的
请问这个函数是什么呢 我没有找到

伽蓝鸟

summersleave 发表于 2021-6-15 12:42
假相当位温是我参考版里的帖子求的
请问这个函数是什么呢 我没有找到

GrADS有很多内置函数的，平流那一系列的函数可以在这里找到/opengrads2.1/Contents/Resources/Documentation/opengrads/doc/udxt/libbjt/
没事的时候也可以戳开其他几个文件夹，了解一下有哪些好用的内置函数

summersleave

伽蓝鸟 发表于 2021-6-15 13:03
GrADS有很多内置函数的，平流那一系列的函数可以在这里找到/opengrads2.1/Contents/Resources/Documentat ...

谢谢 我学习一下