Building the Cloud Parcel Model on Matlab - Part 2: Create the core ODE system (A. basic - completed corrected version of previous) || note: with a single uniform particle size as the input

Here is the fully updated MATLAB script with both adiabatic cooling and condensation-driven latent heat release now included in the temperature tendency.

note: with a single uniform particle size as the input. check the next post for using a PSD as the input instead of a single particle size.

---

✅ Part 1: Driver Code

% Define parameters
params_base = struct( ...
    'P', 85000, ...       % Pressure [Pa]
    'Nd', 100e6, ...      % Droplet number concentration [#/m^3]
    'kappa', 0.3, ...     % Hygroscopicity
    'rho_air', 1.2, ...   % Air density [kg/m^3]
    'wz', 1.0 ...         % Updraft velocity [m/s]
);

% Initial conditions: T0 [K], w0 [kg/kg], ri0 [m] for each droplet
T0 = 273.15 + 0.5;        % Slightly above freezing
w0 = 0.010;               % Initial water vapor mixing ratio
ri0 = 0.05e-6 * ones(50,1);  % 50 small droplets, radius = 0.05 micron
y0 = [T0; w0; ri0];       % Concatenate all initial conditions

% Call ODE solver
[t, y] = ode15s(@(t,y) parcel_ode_kappa(t, y, params_base), [0 1800], y0);

---

✅ Part 2: Core ODE System

function dydt = parcel_ode_kappa(t, y, p)
    T = y(1);
    w = y(2);
    ri = y(3:end);

    S = compute_supersat(T, w, p);
    dri_dt = droplet_growth_rate(ri, S, T, p);
    dT_dt = latent_heating_rate(ri, dri_dt, p, T);
    dw_dt = vapor_consumption_rate(ri, dri_dt, p);

    dydt = [dT_dt; dw_dt; dri_dt];
end

---

✅ Part 3: Supporting Functions

function S = compute_supersat(T, w, p)
    es = 611 * exp((17.27 * (T - 273.15)) / (T - 35.85));
    ws = 0.622 * es / (p.P - es);
    S = (w - ws) / ws;
end

function dri_dt = droplet_growth_rate(ri, S, T, p)
    G = compute_G(T, p);
    ak = kappa_kohler(ri, p);
    Seff = (1 + ak) .* S;
    dri_dt = G .* Seff ./ ri;
end

function G = compute_G(T, p)
    L = 2.5e6;           % latent heat [J/kg]
    Rv = 461;            % gas constant for water vapor [J/kg/K]
    Dv = 2.5e-5;         % diffusivity of water vapor [m^2/s]
    ka = 0.024;          % thermal conductivity [W/m/K]
    rho_w = 1000;        % density of water [kg/m^3]
    es = 611 * exp((17.27 * (T - 273.15)) / (T - 35.85));
    G = (rho_w ./ (Rv * T)) .* (Dv + ka * L^2 / (Rv * T^2));
end

function dT_dt = latent_heating_rate(ri, dri_dt, p, T)
    L = 2.5e6;           % Latent heat [J/kg]
    cp = 1005;           % Specific heat of air [J/kg/K]
    rho_w = 1000;        % Density of water [kg/m^3]
    g = 9.81;            % Gravity [m/s^2]

    % Total condensation rate (kg water / m^3 air / s)
    dmdt = sum(4 * pi * ri.^2 .* p.Nd .* dri_dt) * rho_w;

    % Heating due to condensation
    heating = (L / (p.rho_air * cp)) * dmdt;

    % Cooling due to adiabatic expansion
    if isfield(p, 'wz')
        cooling = -(g / cp) * p.wz;  % in K/s
    else
        cooling = 0;  % no vertical motion if wz not specified
    end

    dT_dt = heating + cooling;
end

function dw_dt = vapor_consumption_rate(ri, dri_dt, p)
    dw_dt = -sum(4 * pi * ri.^2 .* p.Nd .* dri_dt) / p.rho_air;
end

function ak = kappa_kohler(ri, p)
    ak = p.kappa ./ ri;
end

---

🔄 Summary of Changes

• params_base now includes wz, the vertical velocity (e.g., 1 m/s).

• latent_heating_rate now adds condensation heating and subtracts adiabatic cooling using -g/cp * wz.