LES Intercomparison Study for CBL: Sensitivity with respect to Spatial Resolution
Last updated: May 2026
For case setup and physical parameters, see the Description notebook.
Resolutions compared: \(64^3\), \(128^3\), \(256^3\), \(384^3\); SGS: selected by the user.
Setup
The next cells load Python packages, locate the simulation outputs, and define the grid and averaging window used throughout the notebook.
[169]:
import os
import re
import glob
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
Output directories
[170]:
from pathlib import Path
# Base directory (jaxalfa/)
def find_repo_root(start=None):
path = Path(start or ('__file__' in globals() and __file__) or Path.cwd()).resolve()
for candidate in (path, *path.parents):
if (candidate / 'examples').is_dir() and (candidate / 'docs').is_dir():
return candidate
raise FileNotFoundError('Could not locate jaxalfa repository root')
BaseDir = find_repo_root()
def read_config(run_dir):
cfg = {}
exec((run_dir / 'Config.py').read_text(), cfg)
return cfg
optSGS = 1 # LASDD-SM: 1, LASDD-WL: 2, LAD-SM: 3, LAD-WL: 4
sgs_names = {1: 'LASDD-SM', 2: 'LASDD-WL', 3: 'LAD-SM', 4: 'LAD-WL'}
run_styles = {
'64x64x64': {'color': 'red', 'linestyle': '-'},
'128x128x128': {'color': 'blue', 'linestyle': '-'},
'256x256x256': {'color': 'green', 'linestyle': '-'},
'384x384x384 (SP)': {'color': 'black', 'linestyle': '-'},
}
# Resolution sensitivity: 64^3, 128^3, 256^3, 384^3 runs (double precision)
_sgs = {1: 'LASDD_SM', 2: 'LASDD_WL', 3: 'LAD_SM', 4: 'LAD_WL'}
OutputDir1 = BaseDir / f'examples/CBL_N91/runs/64x64x64_{_sgs[optSGS]}_DP/output'
OutputDir2 = BaseDir / f'examples/CBL_N91/runs/128x128x128_{_sgs[optSGS]}_DP/output'
OutputDir3 = BaseDir / f'examples/CBL_N91/runs/256x256x256_{_sgs[optSGS]}_DP/output'
OutputDir4 = BaseDir / f'examples/CBL_N91/runs/384x384x384_{_sgs[optSGS]}_SP/output'
Case configuration
[171]:
cfg_1 = read_config(OutputDir1.parent)
cfg_2 = read_config(OutputDir2.parent)
cfg_3 = read_config(OutputDir3.parent)
cfg_4 = read_config(OutputDir4.parent)
nz_1 = int(cfg_1['nz'])
nz_2 = int(cfg_2['nz'])
nz_3 = int(cfg_3['nz'])
nz_4 = int(cfg_4['nz'])
l_z = float(cfg_1['l_z'])
z_damping = float(cfg_1.get('z_damping', np.nan))
RelaxTime = float(cfg_1.get('RelaxTime', np.nan))
OutputInterval_sec = float(cfg_1.get('OutputInterval_sec', 60.0))
# Averaging window — CBL quasi-steady state (hours 3.0-3.5)
T_start = 3.0 * 3600 # s
T_end = 3.5 * 3600 # s
Derived grid and averaging indices
[172]:
# Half levels — u, v, TH
z_1 = np.array([(k + 0.5) * l_z / (nz_1 - 1) for k in range(nz_1)])
z_2 = np.array([(k + 0.5) * l_z / (nz_2 - 1) for k in range(nz_2)])
z_3 = np.array([(k + 0.5) * l_z / (nz_3 - 1) for k in range(nz_3)])
z_4 = np.array([(k + 0.5) * l_z / (nz_4 - 1) for k in range(nz_4)])
# Full levels — w, uw, vw, wTH, qz
z_w_1 = np.array([k * l_z / (nz_1 - 1) for k in range(nz_1)])
z_w_2 = np.array([k * l_z / (nz_2 - 1) for k in range(nz_2)])
z_w_3 = np.array([k * l_z / (nz_3 - 1) for k in range(nz_3)])
z_w_4 = np.array([k * l_z / (nz_4 - 1) for k in range(nz_4)])
# File indices for the averaging window (same for all runs: same OutputInterval_sec)
T_start_index = int(T_start / OutputInterval_sec) - 1
T_end_index = int(T_end / OutputInterval_sec) - 1
print(f'Averaging window: file indices {T_start_index} – {T_end_index}')
Averaging window: file indices 179 – 209
Statistics loader
[173]:
def LoadStatsAverage(stat_files, T_start_index, T_end_index, nz_expected):
if len(stat_files) == 0:
print(f'No statistics files available; plotting NaN placeholders for nz={nz_expected}.')
nan = np.full(nz_expected, np.nan)
return tuple(nan.copy() for _ in range(15))
U = []; V = []; TH = []
u2 = []; v2 = []; w2 = []; TH2 = []
uv = []; uw = []; vw = []
txy = []; txz = []; tyz = []
wTH = []; qz = []
for f in stat_files:
with np.load(f) as d:
U.append(d['U']); V.append(d['V']); TH.append(d['TH'])
u2.append(d['u2']); v2.append(d['v2']); w2.append(d['w2'])
TH2.append(d['TH2'])
uv.append(d['uv']); uw.append(d['uw']); vw.append(d['vw'])
txy.append(d['txy']); txz.append(d['txz']); tyz.append(d['tyz'])
wTH.append(d['wTH']); qz.append(d['qz'])
U = np.array(U); V = np.array(V); TH = np.array(TH)
u2 = np.array(u2); v2 = np.array(v2); w2 = np.array(w2); TH2 = np.array(TH2)
uv = np.array(uv); uw = np.array(uw); vw = np.array(vw)
txy = np.array(txy); txz = np.array(txz); tyz = np.array(tyz)
wTH = np.array(wTH); qz = np.array(qz)
sl = slice(T_start_index, min(T_end_index + 1, len(stat_files)))
if sl.start >= len(stat_files):
print(f'Averaging window starts after available files; plotting NaN placeholders for nz={nz_expected}.')
nan = np.full(nz_expected, np.nan)
return tuple(nan.copy() for _ in range(15))
return (
np.mean(U[sl], axis=0), np.mean(V[sl], axis=0), np.mean(TH[sl], axis=0),
np.mean(u2[sl], axis=0), np.mean(v2[sl], axis=0), np.mean(w2[sl], axis=0),
np.mean(TH2[sl], axis=0),
np.mean(uv[sl], axis=0), np.mean(uw[sl], axis=0), np.mean(vw[sl], axis=0),
np.mean(txy[sl], axis=0), np.mean(txz[sl], axis=0), np.mean(tyz[sl], axis=0),
np.mean(wTH[sl], axis=0), np.mean(qz[sl], axis=0)
)
Available statistics files
[174]:
def get_stat_files(output_dir):
files = sorted(
glob.glob(str(output_dir / 'ALFA_Statistics_Iteration_*.npz')),
key=lambda x: int(re.search(r'Iteration_(\d+)', x).group(1))
)
return files
StatFiles1 = get_stat_files(OutputDir1)
StatFiles2 = get_stat_files(OutputDir2)
StatFiles3 = get_stat_files(OutputDir3)
StatFiles4 = get_stat_files(OutputDir4)
print(f'64^3 : {len(StatFiles1)} files')
print(f'128^3 : {len(StatFiles2)} files')
print(f'256^3 : {len(StatFiles3)} files')
print(f'384^3 : {len(StatFiles4)} files')
64^3 : 210 files
128^3 : 210 files
256^3 : 210 files
384^3 : 210 files
Temporally averaged profiles
[175]:
(U_avg_1, V_avg_1, TH_avg_1,
u2_avg_1, v2_avg_1, w2_avg_1, TH2_avg_1,
uv_avg_1, uw_avg_1, vw_avg_1,
txy_avg_1, txz_avg_1, tyz_avg_1,
wTH_avg_1, qz_avg_1) = LoadStatsAverage(StatFiles1, T_start_index, T_end_index, nz_1)
(U_avg_2, V_avg_2, TH_avg_2,
u2_avg_2, v2_avg_2, w2_avg_2, TH2_avg_2,
uv_avg_2, uw_avg_2, vw_avg_2,
txy_avg_2, txz_avg_2, tyz_avg_2,
wTH_avg_2, qz_avg_2) = LoadStatsAverage(StatFiles2, T_start_index, T_end_index, nz_2)
(U_avg_3, V_avg_3, TH_avg_3,
u2_avg_3, v2_avg_3, w2_avg_3, TH2_avg_3,
uv_avg_3, uw_avg_3, vw_avg_3,
txy_avg_3, txz_avg_3, tyz_avg_3,
wTH_avg_3, qz_avg_3) = LoadStatsAverage(StatFiles3, T_start_index, T_end_index, nz_3)
(U_avg_4, V_avg_4, TH_avg_4,
u2_avg_4, v2_avg_4, w2_avg_4, TH2_avg_4,
uv_avg_4, uw_avg_4, vw_avg_4,
txy_avg_4, txz_avg_4, tyz_avg_4,
wTH_avg_4, qz_avg_4) = LoadStatsAverage(StatFiles4, T_start_index, T_end_index, nz_4)
S_avg_1 = np.sqrt(U_avg_1**2 + V_avg_1**2)
uw_tot_1 = uw_avg_1 + txz_avg_1
vw_tot_1 = vw_avg_1 + tyz_avg_1
wTH_tot_1 = wTH_avg_1 + qz_avg_1
S_avg_2 = np.sqrt(U_avg_2**2 + V_avg_2**2)
uw_tot_2 = uw_avg_2 + txz_avg_2
vw_tot_2 = vw_avg_2 + tyz_avg_2
wTH_tot_2 = wTH_avg_2 + qz_avg_2
S_avg_3 = np.sqrt(U_avg_3**2 + V_avg_3**2)
uw_tot_3 = uw_avg_3 + txz_avg_3
vw_tot_3 = vw_avg_3 + tyz_avg_3
wTH_tot_3 = wTH_avg_3 + qz_avg_3
S_avg_4 = np.sqrt(U_avg_4**2 + V_avg_4**2)
uw_tot_4 = uw_avg_4 + txz_avg_4
vw_tot_4 = vw_avg_4 + tyz_avg_4
wTH_tot_4 = wTH_avg_4 + qz_avg_4
print(f'Averaging over {T_end_index - T_start_index + 1} files '
f'({T_start/3600:.1f}–{T_end/3600:.1f} h)')
Averaging over 31 files (3.0–3.5 h)
[176]:
plt.rcParams.update({
"text.usetex": True,
"font.size": 14,
"axes.labelsize": 16,
"xtick.labelsize": 12,
"ytick.labelsize": 12
})
[177]:
def plot_profile(x, z, xlabel, ylabel=r"$z$ (m)", linestyle='-k', label=None, ax=None):
if ax is None:
fig, ax = plt.subplots(figsize=(5, 6), constrained_layout=True)
ax.plot(x, z, linestyle, linewidth=2, label=label)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.grid(False)
Mean Potential Temperature and Temperature Variance
The two panels compare the horizontally averaged potential-temperature profile and the resolved temperature variance over the 3.0-3.5 h averaging window.
[178]:
fig, axs = plt.subplots(1, 2, figsize=(10, 5), constrained_layout=True)
for lbl, TH, TH2, z in [('64x64x64', TH_avg_1, TH2_avg_1, z_1),
('128x128x128', TH_avg_2, TH2_avg_2, z_2),
('256x256x256', TH_avg_3, TH2_avg_3, z_3),
('384x384x384 (SP)', TH_avg_4, TH2_avg_4, z_4)]:
style = run_styles[lbl]
axs[0].plot(TH, z, color=style['color'], linestyle=style['linestyle'], linewidth=2, label=lbl)
axs[1].plot(TH2, z, color=style['color'], linestyle=style['linestyle'], linewidth=2, label=lbl)
axs[0].set_xlabel(r"$\langle \theta \rangle$ (K)")
axs[0].set_ylabel(r"$z$ (m)")
axs[0].set_title("Mean Potential Temperature")
axs[1].set_xlabel(r"$\langle \theta^{\prime 2} \rangle$ (K$^2$)")
axs[1].set_ylabel(r"$z$ (m)")
axs[1].set_title("Temperature Variance")
for ax in axs:
ax.axhline(y=z_damping, color='k', linestyle=':', linewidth=1.5)
ax.grid()
ax.legend(frameon=False)
fig.suptitle(f"Potential Temperature Statistics (3.0--3.5 h average): {sgs_names[optSGS]} model", fontsize=18)
plt.show()
Resolved Velocity Variances
The resolved variance profiles indicate how the resolved turbulent kinetic energy is distributed among the three velocity components.
[179]:
fig, axs = plt.subplots(1, 3, figsize=(15, 5), constrained_layout=True)
for lbl, u2, v2, w2, z, z_w in [
('64x64x64', u2_avg_1, v2_avg_1, w2_avg_1, z_1, z_w_1),
('128x128x128', u2_avg_2, v2_avg_2, w2_avg_2, z_2, z_w_2),
('256x256x256', u2_avg_3, v2_avg_3, w2_avg_3, z_3, z_w_3),
('384x384x384 (SP)', u2_avg_4, v2_avg_4, w2_avg_4, z_4, z_w_4),
]:
style = run_styles[lbl]
axs[0].plot(u2, z, color=style['color'], linestyle=style['linestyle'], linewidth=2, label=lbl)
axs[1].plot(v2, z, color=style['color'], linestyle=style['linestyle'], linewidth=2, label=lbl)
axs[2].plot(w2, z_w, color=style['color'], linestyle=style['linestyle'], linewidth=2, label=lbl)
axs[0].set_xlabel(r"Resolved $\sigma_u^2$ (m$^2$/s$^2$)"); axs[0].set_ylabel(r"$z$ (m)")
axs[1].set_xlabel(r"Resolved $\sigma_v^2$ (m$^2$/s$^2$)"); axs[1].set_ylabel(r"$z$ (m)")
axs[2].set_xlabel(r"Resolved $\sigma_w^2$ (m$^2$/s$^2$)"); axs[2].set_ylabel(r"$z$ (m)")
for ax in axs:
ax.axhline(y=z_damping, color='k', linestyle=':', linewidth=1.5)
ax.grid()
ax.legend(frameon=False)
fig.suptitle(f"Resolved Velocity Variances (3.0--3.5 h average): {sgs_names[optSGS]} model", fontsize=18)
plt.show()
Vertical Heat Fluxes
Resolved, SGS, and total heat-flux profiles are compared across the four grid resolutions.
[180]:
fig, axs = plt.subplots(1, 3, figsize=(15, 5), constrained_layout=True)
for lbl, wTH, qz, wTH_tot, z_w in [
('64x64x64', wTH_avg_1, qz_avg_1, wTH_tot_1, z_w_1),
('128x128x128', wTH_avg_2, qz_avg_2, wTH_tot_2, z_w_2),
('256x256x256', wTH_avg_3, qz_avg_3, wTH_tot_3, z_w_3),
('384x384x384 (SP)', wTH_avg_4, qz_avg_4, wTH_tot_4, z_w_4),
]:
style = run_styles[lbl]
axs[0].plot(wTH, z_w, color=style['color'], linestyle=style['linestyle'], linewidth=2, label=lbl)
axs[1].plot(qz, z_w, color=style['color'], linestyle=style['linestyle'], linewidth=2, label=lbl)
axs[2].plot(wTH_tot, z_w, color=style['color'], linestyle=style['linestyle'], linewidth=2, label=lbl)
axs[0].set_xlabel(r"Resolved $\langle w'\Theta' \rangle$ (K m/s)"); axs[0].set_ylabel(r"$z$ (m)"); axs[0].set_title("Resolved")
axs[1].set_xlabel(r"SGS $\langle w'\Theta' \rangle$ (K m/s)"); axs[1].set_ylabel(r"$z$ (m)"); axs[1].set_title("SGS")
axs[2].set_xlabel(r"Total $\langle w'\Theta' \rangle$ (K m/s)"); axs[2].set_ylabel(r"$z$ (m)"); axs[2].set_title("Total (Resolved + SGS)")
for ax in axs:
ax.axhline(y=z_damping, color='k', linestyle=':', linewidth=1.5)
ax.grid()
ax.legend(frameon=False)
fig.suptitle(f"Vertical Heat Fluxes (3.0--3.5 h average): {sgs_names[optSGS]} model", fontsize=18)
plt.show()
