TN波通量相邻地区箭头矛盾

candypower

请教一下各位大佬，我想计算2002年3-6月北半球500百帕环流场上的波活动通量，位势高度场用的是ERA5月平均数据中2002年3-6月这四个月500百帕的平均位势高度，气候平均态用的是ERA5月平均数据中2001-2020年年3-6月这四个月的500百帕的平均位势高度，但是画图出来的结果是有些地区箭头方向是相反的，想知道这是为什么，恳请各位不吝赐教。
代码如下：lon就是0到360，lat就是90到-90，era5数据的经纬度
clear; clc; close all
lon=ncread('D:\数据\Copernicus Climate Change Service\ERA5\全球\data_1.nc', 'longitude');
lat=ncread('D:\数据\Copernicus Climate Change Service\ERA5\全球\data_1.nc', 'latitude');
z = ncread('D:\数据\Copernicus Climate Change Service\ERA5\全球\gh_rh_t_u_v_w_on_100_200_500_700_850_1000.nc', 'z');
z = reshape(squeeze(z(:, :, 4, :)), [length(lon), length(lat), 12, 20]); z = z/9.8;
u = ncread('D:\数据\Copernicus Climate Change Service\ERA5\全球\gh_rh_t_u_v_w_on_100_200_500_700_850_1000.nc', 'u');
u = reshape(squeeze(u(:, :, 4, :)), [length(lon), length(lat), 12, 20]);
v = ncread('D:\数据\Copernicus Climate Change Service\ERA5\全球\gh_rh_t_u_v_w_on_100_200_500_700_850_1000.nc', 'v');
v = reshape(squeeze(v(:, :, 4, :)), [length(lon), length(lat), 12, 20]);
%%
z_anml = squeeze(mean(z(:, :, 3:6, 2), 3, 'omitnan'));
z_mean = squeeze(mean(z(:, :, [size=21.3333px]3:6, :), [3 4], 'omitnan'));
u_mean = squeeze(mean(u(:, :, [size=21.3333px]3:6, :), [3 4], 'omitnan'));
v_mean = squeeze(mean(v(:, :, [size=21.3333px]3:6, :), [3 4], 'omitnan'));

level=500;
[waf_x, waf_y, psi_a]=tn_waf(z_anml, z_mean, u_mean, v_mean, level, lon, lat);
waf_x(sqrt(waf_x.^2+waf_y.^2)<2)=NaN;
waf_y(sqrt(waf_x.^2+waf_y.^2)<2)=NaN;
%%
figure(1)
[X, Y] = meshgrid(lon, lat);
m_proj('Equidistant Cylindrical', 'long', [0 360], 'lat', [15 80]);
hold on
psi_a1=psi_a./1e6;
m_pcolor(X, Y, psi_a1');
m_coast('color','k');
m_quiver(X(1:20:end, 1:20:end), Y(1:20:end, 1:20:end), (waf_x(1:20:end, 1:20:end)'), (waf_y(1:20:end, 1:20:end)'), 2, 'k');
setPivot(0);
m_grid('gridcolor', 'none', 'box', 'on', 'fonts', 15);
colorbar;

%%
function [waf_x,waf_y,psi_a]=tn_waf(hgt,hgt_m,uwnd_m,vwnd_m,level,lon,lat)
%% input variables
% hgt: [lon,lat] % 位势高度的原始场
% hgt_m: [lon lat] % 位势高度的气候态
% uwnd_m: [lon lat] % 风速的气候态
% vwnd_m: [lon lat] %
% level: pressure/(1000hPa) % 气压层
% lon
% lat
%% output variables
% waf_x: T-N wave activity flux zonal component
% waf_y: T-N wave activity flux meridional component
% psi: the geostrophic stream function
%% constant
[nlon, nlat]=size(hgt);

a=6.37*10^6; % earth's radius
Omega=7.292*10^-5; % earth's rotation rate
f=2*Omega.*sind(lat); % Coriolis parameter
p=level/1000;

dy=lat(2)-lat(1);
dphi=deg2rad(dy); % 差分
dx=lon(2)-lon(1);
dlamda=deg2rad(dx); % 差分

cosphi=cosd(lat);
Uspeed=sqrt(uwnd_m.^2+vwnd_m.^2);
% psi anom
f=repmat(f',nlon,1);
psi_a=((hgt-hgt_m).*9.8)./f;
%%
% dpsi/dphi
dpsi_dphi=nan(nlon,nlat);
dpsi_dphi(:,2:nlat-1)=(psi_a(:,3:nlat)-psi_a(:,1:nlat-2))./(2*dphi);
% d2psi/d2phi
d2psi_d2phi=nan(nlon,nlat);
d2psi_d2phi(:,2:nlat-1)=(psi_a(:,3:nlat)-2.*psi_a(:,2:nlat-1)+psi_a(:,1:nlat-2))./(dphi.^2);
% dpsi/dlamda
dpsi_dlamda=nan(nlon,nlat);
dpsi_dlamda(2:nlon-1,:)=(psi_a(3:nlon,:)-psi_a(1:nlon-2,:))./(2*dlamda);
% d2psi/d2lamda
d2psi_d2lamda=nan(nlon,nlat);
d2psi_d2lamda(2:nlon-1,:)=(psi_a(3:nlon,:)-2.*psi_a(2:nlon-1,:)+psi_a(1:nlon-2,:))./(dlamda.^2);
% d2psi/(dlamda dphi)
d2psi_dlamdadphi=nan(nlon,nlat);
d2psi_dlamdadphi(2:nlon-1,2:nlat-1)=(dpsi_dphi(3:nlon,2:nlat-1)-dpsi_dphi(1:nlon-2,2:nlat-1))./(2*dlamda);
%% waf_x
a2=a^2;
cosphi=repmat(cosphi',nlon,1);
cosphi2=cosphi.^2;
% waf x
wafx_term1=(uwnd_m./(a2.*cosphi2)).*((dpsi_dlamda).^2 - psi_a.*(d2psi_d2lamda));
wafx_term2=(vwnd_m./(a2.*cosphi)).*((dpsi_dlamda).*(dpsi_dphi) - psi_a.*(d2psi_dlamdadphi));
waf_x=((p.*cosphi)./(2.*Uspeed)).*(wafx_term1+wafx_term2);
% waf y
wafy_term1=(uwnd_m./(a2.*cosphi)).*((dpsi_dlamda).*(dpsi_dphi) - psi_a.*(d2psi_dlamdadphi));
wafy_term2=(vwnd_m./(a2)).*((dpsi_dphi).^2 - psi_a.*(d2psi_d2phi));
waf_y=((p.*cosphi)./(2.*Uspeed)).*(wafy_term1+wafy_term2);

waf_x(:,lat>=-15 & lat<=15)=NaN;
waf_y(:,lat>=-15 & lat<=15)=NaN;
psi_a(:,lat>=-15 & lat<=15)=NaN;
end
出图如下：