脚本 geost_adjust_2.m 与函数 yprim_adj_2.m 一起,通过使用傅里叶展开扩展问题 M7.4,以检查孤立初始高度扰动形式为 h₀(x)=-h_m/[1+(x/L)²] 的地转调整。这里使用的版本采用 64 个傅里叶模式,并采用快速傅里叶变换算法(FFT)。(FFT 中有 128 个模式,但其中只有一半提供真实信息。)在这种情况下,模型可能仅运行 5 天的积分时间(与前一种情况相比需要大量的计算)。选择初始的经向尺度为 500 千米的扰动,并针对纬度为 15°、30°、45°、60°、75° 和 90° 的情况运行模型(共六次运行)。研究每种情况的动画。注意经向流完全以重力波的形式传播,远离初始扰动。经向流不仅具有传播的重力波分量,还具有气旋环流的定常部分。使用 ginput 命令,通过测量从中心西侧的负速度最大值到中心东侧的正速度最大值的距离,估计最终地转流(v 分量)的经向尺度。绘制显示经向尺度随纬度变化的曲线。(注意在此示例中,gH=400 m² s⁻²。)
其中 geost_adjust_2.m代码为:disp('barotropic Rossby adjustment problem')
disp('initial condition h(x) = -hm/(1+(x/L)^2 ')
clear all
close all
lat = input('give latitude in degrees ');
cor = 2*7.2921e-5*sin(pi*lat/180);
runtime = input('integration time in days ');
time = runtime*24*3600;
L = input('give zonal scale of initial disturbance in km (200 km minimum) ');
Lx = 16000; % Lx is the width of the domain in km
k = 2*pi/(Lx*1.e3); % lowest zonal wavenumber in m-1
csq = 20^2; % shallow water speed squared csq = gH
alph = 1.e-5; % damping rate for u component
N = 128; % number of modes for Fourier transform
x = linspace(0, Lx, N);
xm = Lx/2; % location of maximum disturbance in km
% set the initial conditions on u(x), v(x), p(x)
p0 = -400./(1+((x-xm)/L).^2);
% Fourier transform p0(x) to get P0(s)
disp('please wait while time integration is performed')
P0 = fft(p0,N);
U0 = zeros(size(P0));
V0 =zeros(size(P0));
options = []; % place holder
for s = 1: N/2;
ls = +i*k*(s-1);
als = alph*(s/N)^2;
% solution is saved once every two hours
[t,y] = ode45('yprim_adj_2',[0:7200:time],[U0(s) V0(s) P0(s)],...
options,cor,csq,ls,alph);
U(:,s) = y(:,1);
V(:,s) = y(:,2);
P(:,s) = y(:,3);
end
dt = length(t);
for j = 1:dt;
u(j,:) = real(ifft(U(j,:),N));
v(j,:) = real(ifft(V(j,:),N));
p(j,:) = real(ifft(P(j,:),N));
end
figure(1)
axis square
set(gca,'NextPlot','replacechildren')
for j = 1:dt
subplot(3,1,1), plot(x,u(j,:))
axis([0 Lx -6 6])
xlabel('x (km)'), ylabel('u velocity (m/s)')
str1 = ['time = ' num2str(t(j)/3600) ' hours'];
title(str1, 'Fontsize',12)
subplot(3,1,2), plot(x,v(j,:))
axis([0 Lx -3 3])
xlabel('x (km)'), ylabel( 'v velocity (m/s)')
subplot(3,1,3), plot(x,p(j,:))
axis([0 Lx -400 100])
xlabel('x (km)'), ylabel('height (m)')
H= gcf;
M(:,j) = getframe(H);
end
movie(H,M,1,8)
disp('to repeat movie enter command movie(H,M,2) ')
yprim_adj_2.m 代码为:function dy = yprim_adj_2(t,y,flag,cor,csq,ls,alph)
% ode function for barotropic geostrophic adjustment problem 2
dy = zeros(3,1);
dy(1) = cor*y(2) -ls*y(3) -alph*y(1);
dy(2) = -cor*y(1);
dy(3) =-csq*ls*y(1);
求解一下这个问题
最新发布