计算风应力旋度

wywen

在家园里看到有说用[curlz,~]=curl(lon,lat,taux,tauy)求出旋度，经纬度要转换成距离。但是我用距离求不出来，直接用lon,lat才能求出值，和另一种方法求出的一样，麻烦大家帮忙看看我的代码哪里有错，谢谢！！！

wywen

clc,clear,close all;
lon=ncread('fengchang40_year.nc','lon');
lat=ncread('fengchang40_year.nc','lat');
time=ncread('fengchang40_year.nc','time');
taux40_year=ncread('fengchang40_year.nc','taux40_year');
tauy40_year=ncread('fengchang40_year.nc','tauy40_year');
%方法1
curlmean_year=NaN(40,1);
for i=1:40
    taux1=taux40_year(:,:,i);
    tauy1=tauy40_year(:,:,i);
    [mm,nn]=curl(taux1,tauy1);
    curlmean_year(i)=nanmean(mm(:));
end
% 方法1（距离）
%网格距离
[lon,lat]=meshgrid(lon,lat);
d=2*pi*6371000/360;%一个纬度网格的距离，单位m，(地球半径为6371 km),大约为111km
d_lon=lon*d.*cosd(lat);
d_lat=lat*d;
curlz_year=NaN(40,1);
for i=1:40
    taux1=taux40_year(:,:,i);
    tauy1=tauy40_year(:,:,i);
    [curlz,~]=curl(d_lon',d_lat',taux1,tauy1);
    curlz_year(i)=nanmean(curlz(:));
    disp(i)
end
% 方法2
[dtauy_dx, ~] = gradient(tauy40_year);%求偏导
[~, dtaux_dy] = gradient(taux40_year);
wind_stress_curl =dtauy_dx-dtaux_dy;

curl_year=NaN(40,1);
for i=1:40
    wind_stress_curl11=wind_stress_curl(:,:,i);

   curl_year(i)=nanmean(wind_stress_curl11(:));

end

wywen

wywen 发表于 2024-9-26 18:38
clc,clear,close all;
lon=ncread('fengchang40_year.nc','lon');
lat=ncread('fengchang40_year.nc','la ...

求指导{:eb511:}

wywen

wywen 发表于 2024-9-26 18:38
clc,clear,close all;
lon=ncread('fengchang40_year.nc','lon');
lat=ncread('fengchang40_year.nc','la ...

经纬度范围比较小

MengXC

计算风应力旋度的程序家园里面是有的，看看这个：http://bbs.06climate.com/forum.p ... hlight=%D0%FD%B6%C8

wywen

MengXC 发表于 2024-9-27 09:29
计算风应力旋度的程序家园里面是有的，看看这个：http://bbs.06climate.com/forum.php?mod=viewthread&tid= ...

好的，谢谢