NCL和grads计算涡度和散度的函数差异很大

xo2001

最近想用NCL画图，又怕函数用错，于是就用NCL，grads和自己编程三种方法对比，发现结果差异很大，比如散度，NCL用的是div=uv2dv_cfd(u,v,lat,lon,2)这个函数，自己编程用的也是中央差分的方法，grads用的是hdivg（u，v）的函数，但是结果，自己编程的和NCL在lat=-60到60度之间基本一致，高纬是有些差异。但是grads计算得到的差异就太大了。第一个 是grads画的，第二个是NCL，第三个是自己编程算的，有没有人知道为什么？求助~计算的步骤主要是
grads：
sdfopen  /work/cal/temp3/1999/uwnd.1999.nc
sdfopen /work/cal/temp3/1999/vwnd.1999.nc
d hdivg(uwnd,vwnd.2)

NCL：
ufile=addfile("/work/cal/temp3/1999/uwnd.1999.nc","r")
vfile=addfile("/work/cal/temp3/1999/vwnd.1999.nc","r")
u=ufile->uwnd(178,6,:,:)*0.01+202.66
v=vfile->vwnd(178,6,:,:)*0.01+202.66
lat=ufile->lat
lon=ufile->lon
div1=uv2dv_cfd(u,v,lat,lon,2)
plot=gsn_csm_contour_map(wks,div,res)

编程：
do i=2,nx-1
do j=2,ny-1
dx(i,j)=6400000*5/180*3.14159*cos((90-2.5*(j-1))/180*3.14159)    !2.5度一个格距，中央查分用到前后两个点5度距离
dy(i,j)=6400000*5/180*3.14159                                                           !2.5度一个格距，中央查分用上下两个点5度距离
div(i,j)=  (u(i+1,j)-u(i-1,j))/dx(i,j)+(v(i,j-1)-v(i,j+1))/dy(i,j)                       !计算散度
enddo
enddo

qixiangren34

是个有心人啊，继续找找原因。

river

楼主图太小，看不清楚。贴几张大图吧，还有把你几种方法的过程，用的程序也都放上来吧。大家可以来个大讨论

nanqiyuan

嗯。。对。有心人啊

xo2001

qixiangren34 发表于 2013-6-20 18:25

                               
登录/注册后可看大图

是个有心人啊，继续找找原因。

编程的和NCL的都应该是中央差分法，看了风场的实况，感觉NCL和编程的更接近实况。而涡度我用三个方法计算，NCL和编程也是基本一致，grads虽然有差异，但是没散度这么离谱。

xo2001

river 发表于 2013-6-20 20:06

                               
登录/注册后可看大图

楼主图太小，看不清楚。贴几张大图吧，还有把你几种方法的过程，用的程序也都放上来吧。大家可以来个大讨论

贴好了，我觉得可能是我的grads做错了，或者函数做了一些平滑什么的，因为编程应该是没错的，NCL高纬度可能算法不一样，但是低纬度基本一致。

xo2001

nanqiyuan 发表于 2013-6-20 21:04

                               
登录/注册后可看大图

嗯。。对。有心人啊

最近想用NCL直接计算就省得编程了，所以就想把几个常用的函数了解了解~谁知道刚做个散度涡度就有问题...

尽头的尽头

恳请楼主把错误的原因告诉大家把，方便其他人学习和避免出现这样的错误，非常感谢

xo2001

尽头的尽头 发表于 2013-6-20 22:41

                               
登录/注册后可看大图

恳请楼主把错误的原因告诉大家把，方便其他人学习和避免出现这样的错误，非常感谢

主要是不知道什么错误，所以把我的步骤贴出来给大家看看，不知道是我弄错了还是两个软件就是存在差异。版主有没有遇到过这种情况？上面回复那个弄好了是我把图和步骤贴粗来了...不是弄明白了...

kongfeng0824

NCL好一些。
是ncl内置的算法吗？