通过差分的方法做p坐标系涡度方程求助

无语的人

我用ncl的差分公式center_finite_diff_n把涡度方程的每一项都计算出来，但是最后算出来以后，发现基本每一项都是同号的，也就是说这些项加起来是一个相对大项，不满足涡度方程，我不知道问题出在哪里，是我对涡度方程各项的理解出错了还是程序编错了？希望有大神帮我解答一下，谢谢！


我的数据有u，v，w三个分量的数据，平面维是全球的，经纬度间隔都是1.5°
时间维是2009 July17 0 6 12 18    2009 July18 0 6 12 18

高度维是850 925 1000
也就是三个分量的维数都是[time | 8] x [levels | 3] x [latitude | 121] x [longitude | 240]


这是我的程序：
begin
f1=addfile("2009u.nc","r")
f2=addfile("2009v.nc","r")
f3=addfile("2009w.nc","r")
u=f1->u
v=f2->v
w=f3->w
time=f1->time
levels=f1->levels
lat=f1->latitude
lon=f1->longitude
lat=lat*get_pi("double")/180

;垂直方向相对涡度
d=111000*1.5            ;111000是地球经纬度1度的长度
uy=center_finite_diff_n(u,d,False,0,2)
vx=new(dimsizes(u),"double")
do i=0,dimsizes(lat)-1
        vx(:,:,i,:)=center_finite_diff_n(v(:,:,i,:),abs(d*cos(lat(i))),True,0,2)
end do
vor=vx-uy
copy_VarCoords(u,vor)

;相对涡度局地变化项
dt=3600*6
vor1=center_finite_diff_n(vor,dt,False,0,0)
copy_VarCoords(u,vor1)

;x方向相对涡度平流项
vor2x=new(dimsizes(u),"double")
do i=0,dimsizes(lat)-1
        vor2x(:,:,i,:)=center_finite_diff_n(vor(:,:,i,:),abs(d*cos(lat(i))),True,0,2)*u(:,:,i,:)
end do
copy_VarCoords(u,vor2x)

;y方向相对涡度平流项
vor2y=center_finite_diff_n(vor,d,False,0,2)*v
copy_VarCoords(u,vor2y)

;p方向相对涡度平流项
dp=75*10     ；低层100hPa大概相当于1000米
vor2p=center_finite_diff_n(vor,dp,False,0,1)*w
copy_VarCoords(u,vor2p)

;总相对涡度平流项
vor2=vor2x+vor2y+vor2p
copy_VarCoords(u,vor2)

;地转涡度平流项
f=new(dimsizes(u),"double")
do i=0,dimsizes(lat)-1
        f(:,:,i,:)=2*sin(lat(i))*7.292*10^(-5)           ；7.292*10^(-5)  是地球自转角速度
end do
vor3=center_finite_diff_n(f,d,False,0,2)*v
copy_VarCoords(u,vor3)

;散度项(相对涡度+地转涡度)
ux=new(dimsizes(u),"double")
do i=0,dimsizes(lat)-1
        ux(:,:,i,:)=center_finite_diff_n(u(:,:,i,:),abs(d*cos(lat(i))),True,0,2)
end do
vy=center_finite_diff_n(v,d,False,0,2)
vor4=(ux+vy)*(vor+f)
div2=ux(0,1,:,:)+vy(0,1,:,:)
copy_VarCoords(u,vor4)

;扭转项
wy=center_finite_diff_n(w,d,False,0,2)
wx=new(dimsizes(u),"double")
do i=0,dimsizes(lat)-1
        wx(:,:,i,:)=center_finite_diff_n(w(:,:,i,:),abs(d*cos(lat(i))),True,0,2)
end do
up=center_finite_diff_n(u,dp,False,0,1)
vp=center_finite_diff_n(v,dp,False,0,1)
vor5=wy*up-wx*vp
copy_VarCoords(u,vor5)

;a是将涡度方程中所有项相加，理论上等于0
a=vor1+vor2+vor3+vor4-vor5

end

繁星若尘

动力气象学，涡度方程的应用？

无语的人

基本的项我都基本验证过了，涡度和散度算出来都跟用ncl自带的方程uv2dv_cfd(u,v,u&latitude,u&longitude,2)与uv2vr_cfd(u,v,u&latitude,u&longitude,2)算出来的一致，基本排除差分方法出错，难道是涡度方程的哪项有陷阱不能直接差分？

无语的人

繁星若尘 发表于 2018-12-15 08:56
动力气象学，涡度方程的应用？

是的

边界浪人

你好，请问你这个涡度方程问题解决了吗，我也遇到了类似的问题

无语的人

边界浪人 发表于 2019-12-7 17:18
你好，请问你这个涡度方程问题解决了吗，我也遇到了类似的问题

应该没什么问题，应该是由于差分不能完全代表微分存在误差导致的，是正常的。

y-yang

你好，我刚学习ncl，毕业论文需要用涡度方程进行诊断，刚好看到你发了篇关于涡度方程的帖子，冒昧直接拿来用了，但运行过程中，提示错误，你方便帮我解答一下吗？
fatal:Could not coerce values for operation
fatal:["Execute.c":8637]:Execute: Error occurred at or near line 82 in file /cygdrive/d/NCL/sx0                       3/vor.ncl

82：vor5=wy*up-wx*vp

ABCD

刚好要用，谢谢啦哈哈~

deemo7

所以楼主，最后您的方程是对的吗？还是说需要用积分重新计算

无语的人

deemo7 发表于 2020-12-10 10:48
所以楼主，最后您的方程是对的吗？还是说需要用积分重新计算

应该是对的