TN波通量的计算问题

呆妹小霸王

对TN波通量二维的计算，但自己出来的图，如绿色圈起来的地方，和别人示例给的y方向上有点不对
希望有人可以指点下我的程序哪里有问题


clc;clear;
load('K:\NCEP1\daily\monthly_uvhgt_from_daily.mat')  % 数据为1979~2018

phil = hgt *9.8;  %位势

%% 定义常数
%gas constant
Ra = 290;
%earth radius
a = 6400000;
%纬度和经度梯度
[LAT, LON] = meshgrid(lats, lons);
[~,dlon] = gradient(LON);
[dlat, ~] = gradient(LAT);

dlon = dlon*pi/180;
dlat = dlat*pi/180;

%纬度余弦
coslat = cosd(LAT);
sinlat = sind(LAT);

%科氏常数
f = sw_f(LAT);
f0 = sw_f(43);

lev = 250;  %9

%% 250hPa数据提取
%气候态月平均(选取90~10年作为气候态)

uu_90_10 = uu(:,:,:,133:384);
vv_90_10 = vv(:,:,:, 133:384);
phil_90_10 = phil(:,:,:, 133:384);

u_monthly = zeros(144,73,12);
v_monthly = zeros(144,73,12);
phil_monthly = zeros(144,73,12);
for it = 1:12
    u_monthly(:, :, it) = squeeze(nanmean(uu_90_10(:,:,9,it:12:end),4));
    v_monthly(:, :, it) = squeeze(nanmean(vv_90_10(:,:,9,it:12:end),4));
    phil_monthly(:, :, it) = squeeze(nanmean(phil_90_10(:,:,9,it:12:end),4));
end

%减去气候态月平均
for it = 1:40
    ua(:,:, 1+(it-1)*12:12+(it-1)*12) = squeeze(uu(:,:,9,1+(it-1)*12:12+(it-1)*12)) - u_monthly;
    va(:,:, 1+(it-1)*12:12+(it-1)*12) = squeeze(vv(:,:,9,1+(it-1)*12:12+(it-1)*12)) - v_monthly;
    phila(:,:, 1+(it-1)*12:12+(it-1)*12) = squeeze(phil(:,:,9,1+(it-1)*12:12+(it-1)*12)) - phil_monthly;
end

%% 提取11月的气候态
u_clim = u_monthly(:,:,11);
v_clim = v_monthly(:,:,11);
phil_clim = phil_monthly(:,:,11);

%% 提取14年11月的距平值
u_a = ua(:,:,431);
v_a = va(:,:,431);
phil_a = phila(:,:,431);

%% QG streamfunction
psia = phil_a ./ f ;

%%
%一阶导
[dpsia_dy, dpsia_dx] = gradient(psia);   %gradient函数第一项是纬度求导，第二项是经度
dpsia_dlat = dpsia_dy ./ dlat;
dpsia_dlon = dpsia_dx ./ dlon;

%二阶导
[ddpsia_dxy, ddpsia_dxx] = gradient(dpsia_dlon);
ddpsia_dlonlat = ddpsia_dxy ./ dlat;
ddpsia_dlonlon = ddpsia_dxx ./ dlon;

[ddpsia_dyy, ~] = gradient(dpsia_dlat);
ddpsia_dlatlat = ddpsia_dyy ./ dlat ;

%% 各项
%常数项
magU = sqrt(u_clim.^2 + v_clim.^2);
coeff = (250/1000) .* coslat ./ (2*magU);

termxu = dpsia_dlon .* dpsia_dlon - psia .* ddpsia_dlonlon;
termxv = dpsia_dlon .* dpsia_dlat - psia .* ddpsia_dlonlat;
termyv = dpsia_dlat .* dpsia_dlat - psia .* ddpsia_dlatlat;

%% xy分量

Wx = coeff ./ (a*a*coslat) .* (u_clim.*termxu./coslat + v_clim.*termxv);

Wy = coeff ./ (a*a) .* (u_clim./coslat.*termxv + v_clim.*termyv);

psia(:, 33:41)=nan;
Wx(:, 33:41) = nan;
Wy(:, 33:41) = nan;

Wx(find(magU<5))=nan;
Wy(find(magU<5))=nan;

figure;
D=3;
m_proj('miller','lon',[0 360],'lat',[-90 90]);
h1=m_quiver(LON(1:D:end,1:D:end),LAT(1:D:end,1:D:end),Wx(1:D:end,1:D:end),Wy(1:D:end,1:D:end),0,'k');hold on   
m_contourf(LON,LAT,psia,10,'linestyle','none')
cmocean('balance');
colorbar
m_grid('fontname','Times','fontsize',14); hold on
hU=get(h1,'UData');  %adjust
hV=get(h1,'VData');
qs=1;     %·放大倍数
set(h1,'UData',qs*hU,'VData',qs*hV);
m_coast;

wangshsh06

你好。想问你一下你的问题解决了吗？

膘膘

流函数计算换成实际的流函数

wangshsh06

你好，流函数和实际流函数有什么区别呢，谢谢了

呆妹小霸王

膘膘 发表于 2021-11-6 15:29
流函数计算换成实际的流函数

实际的流函数是什么？ 我是跟着那个网站的例子画的，位势φ除以f作为流函数

膘膘

呆妹小霸王 发表于 2021-11-6 17:44
实际的流函数是什么？ 我是跟着那个网站的例子画的，位势φ除以f作为流函数

用uv去算流函数，你算的是地转流函数

呆妹小霸王

膘膘 发表于 2021-11-7 13:53
用uv去算流函数，你算的是地转流函数

哦哦这样啊。这是大家都这么用吗，我看的文献比较少

膘膘

呆妹小霸王 发表于 2021-11-8 09:49
哦哦这样啊。这是大家都这么用吗，我看的文献比较少

是的，基本都是

呆妹小霸王

膘膘 发表于 2021-11-8 12:59
是的，基本都是

Okkk 谢谢你

何佳田力

你可以对照一下，由于设置的线的间隔，图像并不能完全对应