GABLS1 LES Intercomparison Study for Stable Boundary Layers: Sensitivity with respect to Spatial Resolution
Last updated: June 2026
For case setup and physical parameters, see the Description notebook.
Resolutions compared: \(64^3\), \(128^3\), \(256^3\), \(384^3\); SGS: LASDD-SM (DP). Averaging window: 8–9 h.
Setup
The next cells load Python packages, locate the simulation outputs, and define the grid and averaging window used throughout the notebook.
[1]:
import os
import re
import glob
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
Output directories
[2]:
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
if optSGS == 1:
OutputDir1 = BaseDir / 'examples/SBL_GABLS1/runs/64x64x64_LASDD_SM_DP/output'
OutputDir2 = BaseDir / 'examples/SBL_GABLS1/runs/128x128x128_LASDD_SM_DP/output'
OutputDir3 = BaseDir / 'examples/SBL_GABLS1/runs/256x256x256_LASDD_SM_DP/output'
OutputDir4 = BaseDir / 'examples/SBL_GABLS1/runs/384x384x384_LASDD_SM_SP/output'
elif optSGS == 2:
OutputDir1 = BaseDir / 'examples/SBL_GABLS1/runs/64x64x64_LASDD_WL_DP/output'
OutputDir2 = BaseDir / 'examples/SBL_GABLS1/runs/128x128x128_LASDD_WL_DP/output'
OutputDir3 = BaseDir / 'examples/SBL_GABLS1/runs/256x256x256_LASDD_WL_DP/output'
OutputDir4 = BaseDir / 'examples/SBL_GABLS1/runs/384x384x384_LASDD_WL_DP/output'
elif optSGS == 3:
OutputDir1 = BaseDir / 'examples/SBL_GABLS1/runs/64x64x64_LAD_SM_DP/output'
OutputDir2 = BaseDir / 'examples/SBL_GABLS1/runs/128x128x128_LAD_SM_DP/output'
OutputDir3 = BaseDir / 'examples/SBL_GABLS1/runs/256x256x256_LAD_SM_DP/output'
OutputDir4 = BaseDir / 'examples/SBL_GABLS1/runs/384x384x384_LAD_SM_DP/output'
elif optSGS == 4:
OutputDir1 = BaseDir / 'examples/SBL_GABLS1/runs/64x64x64_LAD_WL_DP/output'
OutputDir2 = BaseDir / 'examples/SBL_GABLS1/runs/128x128x128_LAD_WL_DP/output'
OutputDir3 = BaseDir / 'examples/SBL_GABLS1/runs/256x256x256_LAD_WL_DP/output'
OutputDir4 = BaseDir / 'examples/SBL_GABLS1/runs/384x384x384_LAD_WL_DP/output'
Case configuration
[3]:
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 — GABLS1 quasi-steady state (hours 8–9)
T_start = 8 * 3600 # s
T_end = 9 * 3600 # s
Derived grid and averaging indices
[4]:
# 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 479 – 539
Statistics loader
[5]:
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
[6]:
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 : 540 files
128^3 : 540 files
256^3 : 540 files
384^3 : 181 files
Temporally averaged profiles
[7]:
(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 window starts after available files; plotting NaN placeholders for nz=384.
Averaging over 61 files (8.0–9.0 h)
[8]:
plt.rcParams.update({
"text.usetex": True,
"font.size": 14,
"axes.labelsize": 16,
"xtick.labelsize": 12,
"ytick.labelsize": 12
})
[9]:
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 Wind and Hodograph
The first comparison shows the mean streamwise and cross-stream wind components, the wind-speed magnitude, and the hodograph over the 8-9 h averaging window.
[10]:
fig, axs = plt.subplots(2, 2, figsize=(10, 10), constrained_layout=True)
axs = axs.ravel()
run_styles = {
'64x64x64': {'color': 'red', 'linestyle': '-'},
'128x128x128': {'color': 'blue', 'linestyle': '-'},
'256x256x256':{'color': 'green', 'linestyle': '-'},
'384x384x384':{'color': 'black', 'linestyle': '-'},
}
def plot_run_profile(ax, x, z, xlabel, run_label):
style = run_styles[run_label]
ax.plot(x, z, color=style['color'], linestyle=style['linestyle'], linewidth=2, label=run_label)
ax.set_xlabel(xlabel)
ax.set_ylabel(r"$z$ (m)")
for lbl, U, V, S, z in [('64x64x64', U_avg_1, V_avg_1, S_avg_1, z_1),
('128x128x128', U_avg_2, V_avg_2, S_avg_2, z_2),
('256x256x256', U_avg_3, V_avg_3, S_avg_3, z_3),
('384x384x384', U_avg_4, V_avg_4, S_avg_4, z_4)]:
plot_run_profile(axs[0], U, z, r"$U$ (m/s)", lbl)
plot_run_profile(axs[1], V, z, r"$V$ (m/s)", lbl)
plot_run_profile(axs[2], S, z, r"Wind Speed (m/s)", lbl)
for lbl, U, V, style in [('64x64x64', U_avg_1, V_avg_1, run_styles['64x64x64']),
('128x128x128', U_avg_2, V_avg_2, run_styles['128x128x128']),
('256x256x256', U_avg_3, V_avg_3, run_styles['256x256x256']),
('384x384x384', U_avg_4, V_avg_4, run_styles['384x384x384'])]:
axs[3].plot(U, V, color=style['color'], linestyle='-', marker='o',
linewidth=2, markersize=3, label=lbl)
axs[3].set_xlabel(r"$U$ (m/s)")
axs[3].set_ylabel(r"$V$ (m/s)")
axs[3].set_title('Hodograph')
axs[3].set_aspect('equal')
sgs_names = {1: 'LASDD-SM', 2: 'LASDD-WL', 3: 'LAD-SM', 4: 'LAD-WL'}
for ax in axs:
ax.grid()
ax.legend(frameon=False)
fig.suptitle(f"Mean Wind Profiles and Hodograph (8--9 h average): {sgs_names[optSGS]} model", fontsize=18)
plt.show()
Mean Potential Temperature and Temperature Variance
The two panels compare the horizontally averaged potential-temperature profile and the resolved temperature variance over the 8-9 h averaging window.
[11]:
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',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.grid()
ax.legend(frameon=False)
sgs_names = {1: 'LASDD-SM', 2: 'LASDD-WL', 3: 'LAD-SM', 4: 'LAD-WL'}
fig.suptitle(f"Potential Temperature Statistics (8--9 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.
[12]:
fig, axs = plt.subplots(1, 3, figsize=(15, 5), constrained_layout=True)
plot_profile(u2_avg_1, z_1, xlabel=r"Resolved $\sigma_u^2$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[0], label='64x64x64')
plot_profile(v2_avg_1, z_1, xlabel=r"Resolved $\sigma_v^2$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[1], label='64x64x64')
plot_profile(w2_avg_1, z_w_1, xlabel=r"Resolved $\sigma_w^2$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[2], label='64x64x64')
plot_profile(u2_avg_2, z_2, xlabel=r"Resolved $\sigma_u^2$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[0], label='128x128x128')
plot_profile(v2_avg_2, z_2, xlabel=r"Resolved $\sigma_v^2$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[1], label='128x128x128')
plot_profile(w2_avg_2, z_w_2, xlabel=r"Resolved $\sigma_w^2$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[2], label='128x128x128')
plot_profile(u2_avg_3, z_3, xlabel=r"Resolved $\sigma_u^2$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[0], label='256x256x256')
plot_profile(v2_avg_3, z_3, xlabel=r"Resolved $\sigma_v^2$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[1], label='256x256x256')
plot_profile(w2_avg_3, z_w_3, xlabel=r"Resolved $\sigma_w^2$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[2], label='256x256x256')
plot_profile(u2_avg_4, z_4, xlabel=r"Resolved $\sigma_u^2$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[0], label='384x384x384')
plot_profile(v2_avg_4, z_4, xlabel=r"Resolved $\sigma_v^2$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[1], label='384x384x384')
plot_profile(w2_avg_4, z_w_4, xlabel=r"Resolved $\sigma_w^2$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[2], label='384x384x384')
axs[0].grid(); axs[0].legend(frameon=False)
axs[1].grid(); axs[1].legend(frameon=False)
axs[2].grid(); axs[2].legend(frameon=False)
sgs_names = {1: 'LASDD-SM', 2: 'LASDD-WL', 3: 'LAD-SM', 4: 'LAD-WL'}
fig.suptitle(f"Profiles of Velocity Variances (8--9 h average): {sgs_names[optSGS]} model", fontsize=18)
plt.show()
Resolved Momentum Fluxes
These profiles show the resolved turbulent momentum fluxes before adding the SGS contribution.
[13]:
fig, axs = plt.subplots(1, 3, figsize=(15, 5), constrained_layout=True)
plot_profile(uv_avg_1, z_1, xlabel=r"$\langle u'v' \rangle$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[0], label='64x64x64')
plot_profile(uw_avg_1, z_w_1, xlabel=r"$\langle u'w' \rangle$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[1], label='64x64x64')
plot_profile(vw_avg_1, z_w_1, xlabel=r"$\langle v'w' \rangle$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[2], label='64x64x64')
plot_profile(uv_avg_2, z_2, xlabel=r"$\langle u'v' \rangle$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[0], label='128x128x128')
plot_profile(uw_avg_2, z_w_2, xlabel=r"$\langle u'w' \rangle$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[1], label='128x128x128')
plot_profile(vw_avg_2, z_w_2, xlabel=r"$\langle v'w' \rangle$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[2], label='128x128x128')
plot_profile(uv_avg_3, z_3, xlabel=r"$\langle u'v' \rangle$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[0], label='256x256x256')
plot_profile(uw_avg_3, z_w_3, xlabel=r"$\langle u'w' \rangle$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[1], label='256x256x256')
plot_profile(vw_avg_3, z_w_3, xlabel=r"$\langle v'w' \rangle$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[2], label='256x256x256')
plot_profile(uv_avg_4, z_4, xlabel=r"$\langle u'v' \rangle$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[0], label='384x384x384')
plot_profile(uw_avg_4, z_w_4, xlabel=r"$\langle u'w' \rangle$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[1], label='384x384x384')
plot_profile(vw_avg_4, z_w_4, xlabel=r"$\langle v'w' \rangle$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[2], label='384x384x384')
axs[0].grid(); axs[0].legend(frameon=False)
axs[1].grid(); axs[1].legend(frameon=False)
axs[2].grid(); axs[2].legend(frameon=False)
sgs_names = {1: 'LASDD-SM', 2: 'LASDD-WL', 3: 'LAD-SM', 4: 'LAD-WL'}
fig.suptitle(f"Profiles of Resolved Momentum Fluxes (8--9 h average): {sgs_names[optSGS]} model", fontsize=18)
plt.show()
SGS Stresses
The SGS stress profiles show the modeled stress contribution retained by each grid resolution.
[14]:
fig, axs = plt.subplots(1, 3, figsize=(15, 5), constrained_layout=True)
plot_profile(txy_avg_1, z_1, xlabel=r"$\tau_{xy}$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[0], label='64x64x64')
plot_profile(txz_avg_1, z_w_1, xlabel=r"$\tau_{xz}$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[1], label='64x64x64')
plot_profile(tyz_avg_1, z_w_1, xlabel=r"$\tau_{yz}$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[2], label='64x64x64')
plot_profile(txy_avg_2, z_2, xlabel=r"$\tau_{xy}$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[0], label='128x128x128')
plot_profile(txz_avg_2, z_w_2, xlabel=r"$\tau_{xz}$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[1], label='128x128x128')
plot_profile(tyz_avg_2, z_w_2, xlabel=r"$\tau_{yz}$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[2], label='128x128x128')
plot_profile(txy_avg_3, z_3, xlabel=r"$\tau_{xy}$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[0], label='256x256x256')
plot_profile(txz_avg_3, z_w_3, xlabel=r"$\tau_{xz}$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[1], label='256x256x256')
plot_profile(tyz_avg_3, z_w_3, xlabel=r"$\tau_{yz}$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[2], label='256x256x256')
plot_profile(txy_avg_4, z_4, xlabel=r"$\tau_{xy}$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[0], label='384x384x384')
plot_profile(txz_avg_4, z_w_4, xlabel=r"$\tau_{xz}$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[1], label='384x384x384')
plot_profile(tyz_avg_4, z_w_4, xlabel=r"$\tau_{yz}$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[2], label='384x384x384')
axs[0].grid(); axs[0].legend(frameon=False)
axs[1].grid(); axs[1].legend(frameon=False)
axs[2].grid(); axs[2].legend(frameon=False)
sgs_names = {1: 'LASDD-SM', 2: 'LASDD-WL', 3: 'LAD-SM', 4: 'LAD-WL'}
fig.suptitle(f"Profiles of SGS Stresses (8--9 h average): {sgs_names[optSGS]} model", fontsize=18)
plt.show()
Total Momentum Fluxes
The total vertical momentum flux combines resolved and SGS contributions.
[15]:
fig, axs = plt.subplots(1, 2, figsize=(10, 5), constrained_layout=True)
plot_profile(uw_tot_1, z_w_1, xlabel=r"$\langle u'w' \rangle + \tau_{xz}$ (m$^2$/s$^2$)", linestyle='-r', label='64x64x64', ax=axs[0])
plot_profile(uw_tot_2, z_w_2, xlabel=r"$\langle u'w' \rangle + \tau_{xz}$ (m$^2$/s$^2$)", linestyle='-b', label='128x128x128', ax=axs[0])
plot_profile(uw_tot_3, z_w_3, xlabel=r"$\langle u'w' \rangle + \tau_{xz}$ (m$^2$/s$^2$)", linestyle='-g', label='256x256x256', ax=axs[0])
plot_profile(uw_tot_4, z_w_4, xlabel=r"$\langle u'w' \rangle + \tau_{xz}$ (m$^2$/s$^2$)", linestyle='-k', label='384x384x384', ax=axs[0])
axs[0].grid()
axs[0].legend(frameon=False)
plot_profile(vw_tot_1, z_w_1, xlabel=r"$\langle v'w' \rangle + \tau_{yz}$ (m$^2$/s$^2$)", linestyle='-r', label='64x64x64', ax=axs[1])
plot_profile(vw_tot_2, z_w_2, xlabel=r"$\langle v'w' \rangle + \tau_{yz}$ (m$^2$/s$^2$)", linestyle='-b', label='128x128x128', ax=axs[1])
plot_profile(vw_tot_3, z_w_3, xlabel=r"$\langle v'w' \rangle + \tau_{yz}$ (m$^2$/s$^2$)", linestyle='-g', label='256x256x256', ax=axs[1])
plot_profile(vw_tot_4, z_w_4, xlabel=r"$\langle v'w' \rangle + \tau_{yz}$ (m$^2$/s$^2$)", linestyle='-k', label='384x384x384', ax=axs[1])
axs[1].grid()
axs[1].legend(frameon=False)
sgs_names = {1: 'LASDD-SM', 2: 'LASDD-WL', 3: 'LAD-SM', 4: 'LAD-WL'}
fig.suptitle(f"Total Momentum Fluxes: Resolved + SGS (8--9 h average): {sgs_names[optSGS]} model", fontsize=18)
plt.show()
Vertical Heat Fluxes
Resolved, SGS, and total heat-flux profiles are compared across the three grid resolutions.
[16]:
fig, axs = plt.subplots(1, 3, figsize=(15, 5), constrained_layout=True)
plot_profile(wTH_avg_1, z_w_1, xlabel=r"Resolved $\langle w'\Theta' \rangle$ (K m/s)", linestyle='-r', ax=axs[0], label='64x64x64')
plot_profile(wTH_avg_2, z_w_2, xlabel=r"Resolved $\langle w'\Theta' \rangle$ (K m/s)", linestyle='-b', ax=axs[0], label='128x128x128')
plot_profile(wTH_avg_3, z_w_3, xlabel=r"Resolved $\langle w'\Theta' \rangle$ (K m/s)", linestyle='-g', ax=axs[0], label='256x256x256')
plot_profile(wTH_avg_4, z_w_4, xlabel=r"Resolved $\langle w'\Theta' \rangle$ (K m/s)", linestyle='-k', ax=axs[0], label='384x384x384')
axs[0].grid(); axs[0].legend(frameon=False)
plot_profile(qz_avg_1, z_w_1, xlabel=r"SGS $q_z$ (K m/s)", linestyle='-r', ax=axs[1], label='64x64x64')
plot_profile(qz_avg_2, z_w_2, xlabel=r"SGS $q_z$ (K m/s)", linestyle='-b', ax=axs[1], label='128x128x128')
plot_profile(qz_avg_3, z_w_3, xlabel=r"SGS $q_z$ (K m/s)", linestyle='-g', ax=axs[1], label='256x256x256')
plot_profile(qz_avg_4, z_w_4, xlabel=r"SGS $q_z$ (K m/s)", linestyle='-k', ax=axs[1], label='384x384x384')
axs[1].grid(); axs[1].legend(frameon=False)
plot_profile(wTH_tot_1, z_w_1, xlabel=r"Total $\langle w'\Theta' \rangle$ (K m/s)", linestyle='-r', label='64x64x64', ax=axs[2])
plot_profile(wTH_tot_2, z_w_2, xlabel=r"Total $\langle w'\Theta' \rangle$ (K m/s)", linestyle='-b', label='128x128x128', ax=axs[2])
plot_profile(wTH_tot_3, z_w_3, xlabel=r"Total $\langle w'\Theta' \rangle$ (K m/s)", linestyle='-g', label='256x256x256', ax=axs[2])
plot_profile(wTH_tot_4, z_w_4, xlabel=r"Total $\langle w'\Theta' \rangle$ (K m/s)", linestyle='-k', label='384x384x384', ax=axs[2])
axs[2].grid(); axs[2].legend(frameon=False)
sgs_names = {1: 'LASDD-SM', 2: 'LASDD-WL', 3: 'LAD-SM', 4: 'LAD-WL'}
fig.suptitle(f"Profiles of Vertical Heat Flux (8--9 h average): {sgs_names[optSGS]} model", fontsize=18)
plt.show()
