cdiff的用法

siljechen

在线等，关于cdiff的用法不熟悉。
看了下，cdiff(expr,dim），expr表示进行差分运算的量，dim表示进行差分运算的维数方向，为x,y,z中的任一个字符，
那么我现在算一个对流稳定度的判据， ，假相当位温thetase已经算好了，那这个是不是写成cdiff(thetse,z)?
，我一开始想会不会像平时一样写成cdiff(thetase,lev),grads不认，只有z, 但这是不是差值在等压面上啊？我前面有没有必要加负号。论坛里面很多都是在x,y方向进行中央差分，垂直方向的没搜到！请教各位大神。

river

去看一下官方文档，它一般都会有几个例子的，也有相应的一些解释，还是挺有用的。

cdiff

cdiff(expr,dim)

Performs a centered difference operation on expr in the direction specified by dim. The difference is done in the grid space, and no adjustment is performed for unequally spaced grids. The result value at each grid point is the value at the grid point plus one minus the value at the grid point minus one. The dim argument specifies the dimension over which the difference is to be taken, and is a single character: X, Y, Z, or T.

Result values at the grid boundaries are set to missing.

Usage Notes

Examples

The cdiff function may be used to duplicate the calculation done by the hcurl function:
define dv = cdiff(v,x)
define dx = cdiff(lon,x)*3.1416/180
define du = cdiff(u*cos(lat*3.1416/180),y)
define dy = cdiff(lat,y)*3.1416/180
display (dv/dx-du/dy)/(6.37e6*cos(lat*3.1416/180))
The above example assumes an X-Y varying dimension environment. Note that the intrinsic variables lat and lon give results in degrees and must be converted to radians in the calaculation. Also note the radius of the earth is assumed to be 6.37e6 meters thus the U and V winds are assumed to have units of m/sec.

Temperature advection can be calculated using the cdiff function as follows:
define dtx = cdiff(t,x)
define dty = cdiff(t,y)
define dx = cdiff(lon,x)*3.1416/180
define dy = cdiff(lat,y)*3.1416/180
display -1*( (u*dtx)/(cos(lat*3.1416/180)*dx) + v*dty/dy )/6.37e6
where the variable t is temperature, u is the U component of the wind, and v is the V component of the wind.

siljechen

谢谢river,我要是画出来继续贴

siljechen

继续说下我对cdiff的理解，和师妹讨论后，认为 要用中央差分函数cdiff来写，应该是写成cdiff(thetase,z)/cdiff(p,z) ,由于fnl资料里面我没找到现成的气压场p的物理量，我不知道直接用lev是不是就可以？？（帮看看这点），后来就想要不算 ，其中z可以通过fnl现成的位势高度来算，我想，所以define zz=HGTprs/9.8 ，其中hgtprs是fnl的位势高度，26层数据。这样在gs里面就写成cdiff(thetase,z)/cdiff(zz,z) ,然后出图。
是这样的。看了下文献的量级，也是0.0几，而且那时候分析应该是个稳定层结，即〉0，似乎与实况是吻合的。不过出的图有几个奇怪的地方，第一：只能算到200多，上面没有了。但我看看数据，就算是rh也有21层，能到100hPa的呀。然后再经过中央差分，也不会只有200多，不知道怎么回事。第二：一开始，我想先define a=cdiff(thetase,z), define b=cdiff(zz,z) ,d a/b, 老说我公式有错，怎么改都不行，直接d cdiff(thetase,z)/cdiff(zz,z) 才能画出来，为甚么呀为甚么？？第三：横坐标的设定，set xlab我是可以控制时间输出，但我想有两行，第一行时间，第二行日期，这个怎么设定？求解

以上贴出图和gs,请各位帮看看，首先我对cdiff的理解和用位势换算的理解对不对，第二请帮我看看我的三个问题。
*假相当位温及其垂直递减率计算，使用的是fnl格点资料计算

'reinit'
'open g:\snowstorm\fnl_20130217_00_00_c.ctl'

'set grads off'
'set grid off'
'set parea 1 10 1.5 8'
'set lon 70 135'
'set lat 15 65'
'set z 1 21'
'set t 1 12'
'set clab forced'
*'set mpdset hires cnworld'
*'set mpdset cnworld'
*'set mpdset cnriver'
*'set map 15 1 3'
'set ylopts 1 4 0.12'
'set xlopts 1 4 0.12'
*'set xyrev on'
*'set xlabs 14h,17Feb|20h,17Feb|02h,18Feb|08h,18Feb|14h,18Feb|20h,18Feb|02h,19Feb|08h,19Feb|14h,19Feb|20h,19Feb|'
*'run g:\snowstorm\gs\color.gs'
*'cbarn 0.8 1'
*假相当位温
'define tc=tmpprs-273.15'
'define rh=rhprs'
'define prs=lev'  
'define es=(6.112*exp((17.67*tc)/(tc+243.5)))'
'define qs=0.62197*es/(prs-0.378*es)'
'define q=rh*qs/100'
'define e=prs*q/(0.62197+q)+1e-10'
'define tlcl=55.0+2840.0/(3.5*log(tc+273.16)-log(e)-4.805)'
'define theta=(tc+273.16)*pow((1000/prs),(0.2854*(1.0-0.28*q)))'
'define eqt=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))'

*假相当位温的垂直递减率
'define zz=HGTprs/9.8'

'set lon 120'
'set t 2 11'
'set lat 31'
'set lev 1000 200'
'set zlog on'
'set clab forced'
'set gxout shaded'
'set csmooth on'
'set cmin 0'
'd cdiff(eqt,z)/cdiff(zz,z)'
'set gxout contour'
'd cdiff(eqt,z)/cdiff(zz,z)'

'set gxout contour'
'set cint 3'
'set cthick 6'
'set csmooth on'
*'d eqt'
*垂直速度，绝对涡度
*'d ABSVprs*1e5'
*'d vvelprs/(-10)'
'printim g:\snowstorm\gs\thetase-sid.png white'
;

siljechen

出来的图是这样的。

麦田_smile

river 发表于 2013-7-10 23:01
去看一下官方文档，它一般都会有几个例子的，也有相应的一些解释，还是挺有用的。

cdiff

您好 关于cdiff我有一点疑问 在您给的第一个例子中为什么define du = cdiff(u*cos(lat*3.1416/180),y)
而define dv = cdiff(v,x)
也就是说 为什么du要多乘以cos(lat*3.1416/180) 而dv不用呢

river

麦田_smile 发表于 2014-5-7 22:04
您好 关于cdiff我有一点疑问 在您给的第一个例子中为什么define du = cdiff(u*cos(lat*3.1416/180),y)
...

好久没学习了，我也头大了，您您老人家还是看看书上的公式吧。

278803532

第一问题，因为中央差分cdiff是有边界的，如果你想画图100hPa，那么你的层次必须超过100hPa才行，要不100hPa就没有数据的。

278803532

还有，这个我差分我建议你自己手动一层层算，这样会精确很多，

麦田_smile

278803532 发表于 2014-5-8 10:10
还有，这个我差分我建议你自己手动一层层算，这样会精确很多，

请教下 关于4楼给的例子 为什么前后不一致？有6楼那样的疑惑