请使用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

% 定义高度范围
z1 = 0;
z2 = 3000;
dz = 50;
z = z1:dz:z2;

% 计算温度和压强
T = zeros(size(z));
p = zeros(size(z));
T(1) = 300; % 地面温度，单位：K
p(1) = p0; % 地面压强，单位：Pa
for i = 2:length(z)
    T(i) = T(i-1) + alpha*dz; % 根据温度随高度变化率计算温度
    p(i) = p(i-1)*exp(-g/(Rd*T(i))*dz); % 根据理想气体状态方程计算压强
end

% 计算密度和位势温度
rho = p./(Rd*T); % 根据理想气体状态方程计算密度
theta = T.*(p0./p).^(-Rd/Cp); % 计算位势温度

% 计算饱和水汽压和相对湿度
es = 611*exp(Lv/Rv*(1/T0-1./T)); % 饱和水汽压力，单位：Pa
e = es.*0.6; % 实际水汽压力，假设相对湿度为60%
rh = e./es; % 计算相对湿度

% 绘制大气边界层过程图
figure;
subplot(2,2,1);
plot(T,z);
xlabel('温度 (K)');
ylabel('高度 (m)');
title('温度垂直分布');

subplot(2,2,2);
plot(p/100,z);
xlabel('压强 (hPa)');
ylabel('高度 (m)');
title('压强垂直分布');

subplot(2,2,3);
plot(rho,z);
xlabel('密度 (kg/m^3)');
ylabel('高度 (m)');
title('密度垂直分布');

subplot(2,2,4);
plot(theta,z);
xlabel('位势温度 (K)');
ylabel('高度 (m)');
title('位势温度垂直分布');

这段代码会生成一个四个子图的图形窗口，分别展示温度、压强、密度和位势温度随高度变化的垂直分布。您可以根据需要修改高度范围、计算参数和绘图样式。