请使用matlab2021绘制一幅stull所著的大气边界层过程图，要求包含稳定边界层/对流边界层/残余层等，横坐标为时间，纵坐标为高度。请给出具体代码

以下是使用 Matlab 2021 绘制 Stull 大气边界层过程图的示例代码，包括稳定边界层、对流边界层和残余层的垂直分布：

% 定义常数
g = 9.81; % 重力加速度，单位：m/s^2
Cp = 1004; % 空气定压比热，单位：J/(kg K)
Lv = 2.5e6; % 水的汽化潜热，单位：J/kg
Rd = 287; % 干空气气体常数，单位：J/(kg K)
Rv = 461; % 水蒸气气体常数，单位：J/(kg K)
T0 = 273; % 零度绝对温度，单位：K
p0 = 100000; % 参考压强，单位：Pa
rho0 = 1.2; % 参考密度，单位：kg/m^3
alpha = -0.0065; % 温度随高度变化率，单位：K/m

% 定义时间和高度范围
t1 = 0;
t2 = 24*3600;
dt = 3600;
t = t1:dt:t2;
z1 = 0;
z2 = 3000;
dz = 50;
z = z1:dz:z2;

% 初始化参数数组
T = zeros(length(z), length(t));
p = zeros(length(z), length(t));
rho = zeros(length(z), length(t));
theta = zeros(length(z), length(t));
es = zeros(length(z), length(t));
e = zeros(length(z), length(t));
rh = zeros(length(z), length(t));
q = zeros(length(z), length(t));
L = zeros(length(t), 1);
H = zeros(length(t), 1);
zinv = zeros(length(t), 1);

% 计算初始状态
T(:,1) = 300; % 地面温度，单位：K
p(:,1) = p0; % 地面压强，单位：Pa
rho(:,1) = rho0; % 参考密度，单位：kg/m^3
theta(:,1) = T(:,1).*(p0./p(:,1)).^(-Rd/Cp); % 计算位势温度
es(:,1) = 611*exp(Lv/Rv*(1/T0-1./T(:,1))); % 饱和水汽压力，单位：Pa
e(:,1) = es(:,1).*0.6; % 实际水汽压力，假设相对湿度为60%
rh(:,1) = e(:,1)./es(:,1); % 计算相对湿度
q(:,1) = rh(:,1).*e(:,1)./(Rd*T(:,1)); % 计算比湿

% 循环计算边界层过程
for i = 2:length(t)
    % 计算温度和压强
    T(:,i) = T(:,i-1) + alpha*dz; % 根据温度随高度变化率计算温度
    p(:,i) = p(:,i-1)*exp(-g./(Rd*T(:,i))*dz); % 根据理想气体状态方程计算压强
    rho(:,i) = p(:,i)./(Rd*T(:,i)); % 根据理想气体状态方程计算密度
    theta(:,i) = T(:,i).*(p0./p(:,i)).^(-Rd/Cp); % 计算位势温度
    
    % 计算饱和水汽压和相对湿度
    es(:,i) = 611*exp(Lv/Rv*(1/T0-1./T(:,i))); % 饱和水汽压力，单位：Pa
    e(:,i) = es(:,i).*0.6; % 实际水汽压力，假设相对湿度为60%
    rh(:,i) = e(:,i)./es(:,i); % 计算相对湿度
    q(:,i) = rh(:,i).*e(:,i)./(Rd*T(:,i)); % 计算比湿
    
    % 计算湍流动能和潜热通量
    ustar = 0.2; % 假设摩擦速度为0.2 m/s
    wstar = ustar*0.4; % 根据Monin-Obukhov相似理论计算垂直速度标准差
    L(i) = -ustar^3*rho0/(k*g*(q(:,i)*Cp+0.61*q(:,i)*Lv/Rd).*wstar); % 计算稳定度参数
    if L(i) > 0 % 稳定边界层
        H(i) = -ustar^3/(k*g*L(i)); % 计算边界层深度
        zinv(i) = z(1) + (z2-z(1))/H(i)*tanh(H(i)/z2*z); % 计算折射指数
        if mod(i, 6) == 2 % 每隔6个小时绘制一次稳定边界层图像
            figure;
            plot(t/3600, zinv(i,:), 'k');
            xlabel('时间 (h)');
            ylabel('高度 (m)');
            title('稳定边界层垂直分布');
        end
    else % 对流边界层
        H(i) = 1000; % 假设对流层深度为1000 m
        if mod(i, 6) == 2 % 每隔6个小时绘制一次对流边界层图像
            figure;
            contourf(t/3600, z, theta', 50, 'LineStyle', 'none');
            xlabel('时间 (h)');
            ylabel('高度 (m)');
            title('对流边界层垂直分布');
            colorbar;
        end
    end
end

% 计算残余层
zres = 10000; % 残余层顶高度，单位：m
Tres = T(end,end); % 残余层顶温度，单位：K
zresvec = z2:dz:zres; % 残余层高度数组
Tresvec = Tres*ones(size(zresvec)); % 残余层温度数组

% 绘制残余层图像
figure;
plot(t/3600, zinv(end,:), 'k', t/3600, zresvec, 'k--');
xlabel('时间 (h)');
ylabel('高度 (m)');
title('残余层垂直分布');
legend('折射指数', '残余层顶');

这段代码会生成多个图形窗口，包括稳定边界层、对流边界层和残余层的垂直分布。您可以根据需要修改时间和高度范围、计算参数和绘图样式。