- 积分
- 4529
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2018-5-6
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 呆妹小霸王 于 2021-11-3 09:58 编辑
对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;
|
-
-
|