风应力旋度

wusc

大佬们好：
小弟刚刚接触MATLAB，最近要画一个风应力旋度进而求解SverdrupTransoprt。
参照了一些程序


得到的图很奇怪，希望得到指正，不胜感谢

wusc

clear;
    close all;
    load('lon.mat');
    load('lat.mat');
    load('tauu.mat');
    load('tauv.mat');
    R=6.3781e6; % earth's radius
    [lon , lat]   = meshgrid(lat,lon);
    [dx, ~]=gradient(lon);
    dx=dx.*(pi./180).*R.*cos(lat*pi./180);
   
    [~, dy]=gradient(lat);
    dy=dy.*(pi./180).*R;
  
    [~, du]=gradient(tauu(:,:,1,1));
    [dv, ~]=gradient(tauv(:,:,1,1));
    curlz=dv./dx-du./dy+tauu(:,:,1,1).*tan(lat.*pi./180)./R;
    curlz(abs(lat)==90)=NaN;
    contour(curlz);
这是我的代码

wusc

这是我绘制的图

callmezoro

为什么不直接尝试一下MATLAB里的curl函数呢，把经纬度meshgrid之后，直接用[curlz,~]=curl(lon,lat,taux,tauy);就可以求出旋度了，当然经纬度要转换成距离，至于后面求解SverdrupTransoprt我就不太了解了

callmezoro

我也看到了你这个方法的帖子，但是我的感觉是只要确保你的个点数据对应的不是经纬度，而是实际距离，那么地球曲率在其中产生的影响就不是那么大，即使在极地海区也是可以接受的，所以不一定非要采用球坐标系的计算方法

wusc

callmezoro 发表于 2018-4-23 15:54
为什么不直接尝试一下MATLAB里的curl函数呢，把经纬度meshgrid之后，直接用[curlz,~]=curl(lon,lat,taux,ta ...

谢谢，直接用图形和网格都不对不知道为什么

sxq

楼主解决了吗