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
出图如下：

膘膘

首先，你看2002年3-6月的TN就有些不合理。中纬度波动位置受急流经向移动影响显著，春夏之间急流经向位置变化显著，必然导致纬向波动不明显。从HGT的结果看，程序基本没有问题。此外，欧亚大陆上为SRP或者CGT，由于是500hPa的结果，必然在青藏高原附近信号较弱。另外，北大西洋上为SRP上游或者NAO，由于没有地形，所以信号较强。

qingqi

同蹲，我也有这种情况

candypower

qingqi 发表于 2025-4-2 10:24
同蹲，我也有这种情况

你好请问你解决了吗

膘膘

地转风是用位势高度算的吗？

candypower

膘膘 发表于 2025-4-11 14:24
地转风是用位势高度算的吗？

对的，我看有帖子的评论说严格意义上就是要用位势高度算，但是也有人说其他人画的都是用uv风速算的

candypower

膘膘 发表于 2025-4-11 14:24
地转风是用位势高度算的吗？

对的，我看有帖子的评论说严格意义上就是要用位势高度算，但是也有人说其他人画的都是用uv风速算的

膘膘

candypower 发表于 2025-4-14 11:34
对的，我看有帖子的评论说严格意义上就是要用位势高度算，但是也有人说其他人画的都是用uv风速算的

你用实际风算算看

candypower

膘膘 发表于 2025-4-14 16:44
你用实际风算算看

画出来效果还是不太好：
流函数用的是http://bbs.06climate.com/forum.php?mod=viewthread&tid=20777这个帖子中的，代码改写了http://bbs.06climate.com/forum.p ... =%B2%A8%CD%A8%C1%BF这个帖子中的，如下：
lon = ncread('D:\数据\Copernicus Climate Change Service\ERA5\全球\data_1.nc', 'longitude');
lat = ncread('D:\数据\Copernicus Climate Change Service\ERA5\全球\data_1.nc', 'latitude');
lon(lon<270) = lon(lon<270) + 360;
lon = lon - 360;
[lon, order] = sort(lon);
lon_area = [-90 90];
lat_area = [20 86];
ind_lat_area = find(lat>=lat_area(1)&lat<lat_area(2));
ind_lon_area = find(lon>=lon_area(1)&lon<=lon_area(2));
%% 读取数据(将数据恢复到lon*lat，mat格式的数据为纬度*经度)
z500 = load('D:\Code\mat格式数据\气象\全球\z500.mat').z500;
z500 = z500(:, order, :, :); z500 = permute(z500, [2 1 3 4]);
u500 = load('D:\Code\mat格式数据\气象\全球\u500.mat').u500;
u500 = u500(:, order, :, :); u500 = permute(u500, [2 1 3 4]);
v500 = load('D:\Code\mat格式数据\气象\全球\v500.mat').v500;
v500 = v500(:, order, :, :); v500 = permute(v500, [2 1 3 4]);
%% 定义常数
%gas constant
Ra = 290;
%earth radius
a = 6400000;
%纬度和经度梯度
[X, Y] = meshgrid(lat, lon);
[~,dlon] = gradient(Y);
[dlat, ~] = gradient(X);
dlon = dlon*pi/180;
dlat = dlat*pi/180;
%纬度余弦
coslat = cosd(X);
sinlat = sind(X);
%科氏常数
f = sw_f(X);
f0 = sw_f(43);
lev = 500;
%% 处理数据
z_anom = mean(z500(:, :, 3:6, 2), 3, 'omitnan') - mean(z500(:, :, 3:6, :), [3 4], 'omitnan');
u_clim = mean(u500(:, :, 3:6, :), [3 4], 'omitnan');
v_clim = mean(v500(:, :, 3:6, :), [3 4], 'omitnan');
u_anom = mean(u500(:, :, 3:6, 2), 3, 'omitnan') - mean(u500(:, :, 3:6, :), [3 4], 'omitnan');
v_anom = mean(v500(:, :, 3:6, 2), 3, 'omitnan') - mean(v500(:, :, 3:6, :), [3 4], 'omitnan');
phil_anom = z_anom *9.8;  
%% QG streamfunction
% psia = phil_anom ./ f ;
[psia, upsi, vpsi]=psi_streamfunction(Y', X', u_anom', v_anom');
psia = psia';
[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 = (lev/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;
Wx = coeff ./ (a*a*coslat) .* (u_clim.*termxu./coslat + v_clim.*termxv);
Wy = coeff ./ (a*a) .* (u_clim./coslat.*termxv + v_clim.*termyv);
Wx(magU<6)=nan;
Wy(magU<6)=nan;
Wx = Wx(ind_lon_area, ind_lat_area);
Wy = Wy(ind_lon_area, ind_lat_area);
psia = psia(ind_lon_area, ind_lat_area);
psia = psia./1e6;

candypower

candypower 发表于 2025-4-14 21:05
画出来效果还是不太好：
流函数用的是http://bbs.06climate.com/forum.php?mod=viewthread&tid=20777这 ...

画图代码如下：
[X, Y] = meshgrid(lon(ind_lon_area), lat(ind_lat_area));
m_proj('Miller Cylindrical', 'long', lon_area, 'lat', [20 80]);
hold on;
m_pcolor(X, Y, psia');
hold on;
shading interp
m_coast('linewidth', 0.8, 'color', [0.7 0.7 0.7]);
hold on;
m_grid('gridcolor', 'none', 'box', 'on', 'fontsize', 15);
hold on;
% m_quiver(X(1:20:end, 1:20:end), Y(1:20:end, 1:20:end), Fx(1:20:end, 1:20:end), (Fy(1:20:end, 1:20:end))', 5);
m_vec(5, X(1:20:end, 1:20:end), Y(1:20:end, 1:20:end), (Wx(1:20:end, 1:20:end))', (Wy(1:20:end, 1:20:end))', ...
    'headlength', 1.5, 'shaftwidth', 0.2, 'edgeclip', 'on');
hold on;
cm = m_colmap('diverging');
colormap(cm);
setPivot(0);
cb = colorbar;