[MITgcm-support] RBCS extending to bottom & shelf dry points in lake setup
Fateme Sharifi
fateme.sharifi1391 at gmail.com
Mon Sep 1 01:49:00 EDT 2025
Dear MITgcm users,
I am setting up a small lake experiment with *Gaussian bathymetry, a dry
shelf band, and a shallow “moat” region (within 400 m of shoreline, depth ≤
5 m) where I apply RBCS restoring toward 4 °C*.
My input files are generated with a MATLAB script (attached) that writes:
-
Depth.bin (outside/shelf = 0 m),
-
hFacC.bin with proper partial cells,
-
rbcs_mask.bin limited to shallow moat layers only (z_mid ≤ moat_depth,
hFacC > 0),
-
rbcs_Tr1_fld.bin set to 4 °C inside the mask.
At model start (t=1), the RBCS mask behaves correctly (restricted to the
top ~5 m near the boundary). However, after just 2–3 timesteps, the RBCS
influence appears to extend vertically to the bottom (see attached figure
with Y–Z sections at t=1 and t=3). This was unexpected because the written
rbcs_mask.bin only has ones in shallow moat layers.
I also noticed some sensitivity to shelf dry points: despite setting Depth=0
and checking that hFacC=0 in the shelf region, I worry that something in
RBCS or partial cells may be activating those points during runtime.
*My questions:*
1.
Has anyone seen RBCS restoring “leak” into deeper layers, even if the
input mask is strictly shallow?
2.
Is there an additional runtime parameter (e.g. in data.rbcs or data)
needed to prevent the restoring from being applied throughout the column?
3.
What is the recommended way to guarantee that RBCS remains limited to
the intended top layers throughout the run (not just in the input mask)?
4.
Any best practices for handling dry shelf points with RBCS in lakes or
enclosed basins?
I have attached both my MATLAB setup script (
setup_gaussian_lake_rbcs_revised_v2.m) and the figure (
vertical_section_temperature.jpg).
Best regards,
Fatemeh Sadat Sharifi
Ph.D. Candidate at
The Leibniz Institute of Freshwater
Ecology and Inland Fisheries (IGB)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20250901/2c8b6f30/attachment-0001.html>
-------------- next part --------------
% setup_gaussian_lake_rbcs_revised_v2.m
% Fully revised MATLAB script for MITgcm lake setup
% - Keeps user parameters unchanged
% - Gaussian bathymetry, dry shelf (last 500 m) and dry outside
% - Moat: band within 400 m of the shoreline inside lake and only in shallow layers (<=5 m)
% - Initial T = 4 degC ONLY in the moat (for MITgcm binaries we do NOT write NaNs)
% - BUT the script also writes separate NaN-masked MATLAB files (and PNGs) for
% "lake-setup" use/inspection where outside-lake values are set to NaN.
% - RBC mask/target written for iterations: base (no suffix), .0000001000, and .0000000000..max_iter
% - Writes T.init, S.init (no NaNs in binaries!), Depth.bin, hFacC.bin, delR.bin, rbcs_mask.bin, rbcs_Tr1_fld.bin
% - Includes a robust shelf diagnostic (plots and text) proving the shelf is dry
% - Appends a short MOC helper function to compute MOC from GM_PsiY, dxG and deepFacF
clear; close all; clc;
%% ---------------- User Parameters (kept unchanged) ----------------
R = 3440; dx = 20; Nz = 30;
dz_top = 0.5; dz_bot = 2.0; Hmax = 57;
sigma = R/2.5; shelf_width = 500; h_min_edge = 5;
T_surf = 0.001; T_bot = 2.0; T_target = 4.0; % T_target used only in moat layers
moat_width = 400; moat_depth = 5.0; % moat band width and max depth for RBC
S_bg = 1e-4; hFacMin = 0.0; min_thick = 1e-6;
prec = 'real*8'; ieee = 'ieee-be';
max_iter = 10; % RBC iterative write (0..max_iter)
%% ---------------- Horizontal Grid ----------------
x = -R:dx:R; y = -R:dx:R; Nx = numel(x); Ny = numel(y);
[XC,YC] = meshgrid(x,y); Rgrid = sqrt(XC.^2 + YC.^2);
delX = dx * ones(Ny,Nx); delY = dx * ones(Ny,Nx);
%% ---------------- Vertical Grid ----------------
rawDelZ = linspace(dz_top,dz_bot,Nz)';
delZ = rawDelZ * (Hmax / sum(rawDelZ)); % stretches to match Hmax exactly
z_bot = cumsum(delZ);
z_top = [0; z_bot(1:end-1)];
z_mid = 0.5 * (z_top + z_bot);
%% ---------------- Bathymetry (Gaussian, dry shelf & dry outside) ----------------
% Full Gaussian bathymetry
bathy = Hmax .* exp(-(Rgrid.^2) ./ (2 * sigma^2));
% Make the outer 'shelf_width' completely dry (Depth = 0)
bathy(Rgrid >= (R - shelf_width)) = 0; % dry shelf band
bathy(Rgrid > R) = 0; % dry beyond domain radius (outside lake)
% Remove any very thin puddles near inner/outer edges
bathy(bathy <= h_min_edge) = 0;
% lake mask (logical) where bathy > 0
lake_mask = bathy > 0;
%% ---------------- hFacC (partial cells) ----------------
hFacC = zeros(Ny,Nx,Nz);
for k = 1:Nz
fully = bathy >= z_bot(k);
empty = bathy <= (z_top(k) + min_thick);
part = ~(fully | empty);
layer = zeros(Ny,Nx);
layer(fully) = 1;
layer(part) = (bathy(part) - z_top(k)) ./ delZ(k);
% enforce physical bounds and hFacMin behavior
layer(layer < hFacMin) = 0;
layer(layer > 1) = 1;
layer(empty) = 0;
hFacC(:,:,k) = layer;
end
%% ---------------- Moat (2D band inside shoreline) ----------------
% Effective shoreline radius of the wet lake (should equal ~R - shelf_width)
if any(lake_mask,'all')
R_eff = max(Rgrid(lake_mask));
else
error('No wet cells found: check R, dx, shelf_width and bathymetry settings');
end
moat_inner = max(R_eff - moat_width, 0); % inner edge of moat
moat2d = (Rgrid >= moat_inner) & (Rgrid <= R_eff) & lake_mask; % strictly inside lake
%% ---------------- RBC 3D mask (STRICT: only shallow layers & wet) ----------------
% Key fix: ensure rbcs never goes to bottom. rbcs limited to layers whose mid-point <= moat_depth
rbcs3d = false(Ny,Nx,Nz);
for k = 1:Nz
shallow_layer = (z_mid(k) <= moat_depth); % choose layers whose midpoint is <= moat_depth
rbcs3d(:,:,k) = moat2d & shallow_layer & (hFacC(:,:,k) > 0);
end
% Ensure RBC never includes points where bathy==0 (shelf or outside)
rbcs3d = rbcs3d & repmat(lake_mask, [1 1 Nz]);
% Extra safety: ensure rbcs3d top/bottom indices are within hFacC>0 (already used) and do not include
% levels where hFacC is extremely small (below min_thick via hFacMin)
%% ---------------- Initial T/S (SAFE for binaries) ----------------
% MITgcm requires numeric values (no NaN) in binary input files.
T_write = repmat(T_surf, [Ny,Nx,Nz]);
S_write = repmat(S_bg, [Ny,Nx,Nz]);
for k = 1:Nz
wet = (hFacC(:,:,k) > 0);
% linear vertical increase inside lake to bottom (only where wet)
T_write(:,:,k) = T_write(:,:,k) + (T_bot - T_surf) * (z_mid(k)/Hmax) .* double(wet);
end
% Set the moat target temperature ONLY where rbcs3d is true (for writing rbcs_Tr1_fld)
% But keep T_write (initial) having moat temperature in moat layers so T.init reflects the
% intended initial moat = 4 degC in those cells
T_write(rbcs3d) = T_target;
% S_write stays constant (small S_bg) everywhere; dry cells remain S_bg
% Create plotting / "lake-setup" copies that *do* contain NaN outside lake for your other tools
T_NaN = T_write; S_NaN = S_write; Depth_NaN = bathy;
T_NaN(repmat(~lake_mask, [1 1 Nz])) = NaN; % for external inspection (NOT for MITgcm)
S_NaN(repmat(~lake_mask, [1 1 Nz])) = NaN;
Depth_NaN(~lake_mask) = NaN;
% Save NaN-masked .mat files for your `lake-setup` tools that expect NaNs
save('T_init_withNaN.mat','T_NaN','-v7.3');
save('S_init_withNaN.mat','S_NaN','-v7.3');
save('Depth_withNaN.mat','Depth_NaN','-v7.3');
%% ---------------- RBC mask/target arrays for writing ----------------
rbcs_T_write = repmat(T_surf, [Ny,Nx,Nz]); % default (not used by RBC where mask=0)
rbcs_T_write(rbcs3d) = T_target; % target only in moat shallow layers
rbcs_mask_write = zeros(Ny,Nx,Nz);
rbcs_mask_write(rbcs3d) = 1; % mask exactly matches moat layers
% For plotting, create NaN-masked versions
rbcs_mask_plot = rbcs_mask_write;
rbcs_mask_plot(repmat(~lake_mask, [1 1 Nz])) = NaN; % show NaN outside lake for display only
%% ---------------- Debugging / Sanity Checks for Shelf (Depth=0 -> hFacC=0) ----------------
% Identify shelf band indices (radial between R - shelf_width and R)
shelf_idx = (Rgrid >= (R - shelf_width)) & (Rgrid <= R);
% Check bathy in shelf region (should be exactly zero)
num_nonzero_bathy_on_shelf = nnz(bathy(shelf_idx) > 0);
% Check any hFacC > 0 inside shelf (should be zero)
hFacC_on_shelf = false(Nz,1);
for k = 1:Nz
tmp = hFacC(:,:,k);
if any(tmp(shelf_idx) > 0)
hFacC_on_shelf(k) = true;
else
hFacC_on_shelf(k) = false;
end
end
any_hFacC_on_shelf = any(hFacC_on_shelf);
% Compose diagnostic messages
diag_msg = {};
if num_nonzero_bathy_on_shelf > 0 || any_hFacC_on_shelf
diag_msg{end+1} = sprintf('WARNING: Detected %d non-zero bathy cells inside shelf band.', num_nonzero_bathy_on_shelf);
if any_hFacC_on_shelf
diag_msg{end+1} = 'WARNING: Some shelf grid cells have hFacC > 0 (unexpected).';
end
else
diag_msg{end+1} = 'Shelf diagnostic: OK — bathy==0 and hFacC==0 everywhere in the configured shelf band.';
end
% Save diagnostic text
fid = fopen('shelf_diagnostic.txt','w');
for ii = 1:numel(diag_msg)
fprintf(fid, '%s\n', diag_msg{ii});
end
fclose(fid);
% print to console also
for ii = 1:numel(diag_msg)
fprintf('%s\n', diag_msg{ii});
end
% Also produce a small diagnostic PNG explicitly showing shelf mask & any problematic cells
fig_shelf = figure('Color','w','Name','Shelf diagnostic','NumberTitle','off','Units','normalized','Position',[0.1 0.1 0.4 0.4]);
subplot(1,2,1);
bathy_plot = bathy; bathy_plot(~lake_mask) = NaN; pcolor(XC,YC,bathy_plot); shading flat; axis equal tight; colorbar; title('Bathymetry (NaN outside lake)');
subplot(1,2,2);
imagesc(XC(1,:),YC(:,1), double(shelf_idx)); axis equal tight; set(gca,'YDir','normal'); colorbar; title('Shelf band (1 = shelf region)');
saveas(fig_shelf,'shelf_diagnostic.png');
%% ---------------- Write MITgcm binaries (SAFE: no NaNs) ----------------
write_binary('delX.bin', delX, prec, ieee);
write_binary('delY.bin', delY, prec, ieee);
write_binary('XC.bin', XC, prec, ieee);
write_binary('YC.bin', YC, prec, ieee);
write_binary('Depth.bin', bathy, prec, ieee); % outside lake = 0, shelf = 0
write_binary('hFacC.bin', hFacC, prec, ieee);
% IMPORTANT: MITgcm expects 1-D delR (Nz elements)
write_binary('delR.bin', delZ, prec, ieee);
% Initial conditions (NO NaNs in files!)
write_binary('T.init', T_write, prec, ieee);
write_binary('S.init', S_write, prec, ieee);
% RBC base files
write_binary('rbcs_mask.bin', rbcs_mask_write, prec, ieee);
write_binary('rbcs_Tr1_fld.bin',rbcs_T_write, prec, ieee);
% Exactly one file with suffix '0000001000'
iter_str = sprintf('%010d', 1000);
write_binary(['rbcs_Tr1_fld.bin.' iter_str], rbcs_T_write, prec, ieee);
write_binary(['rbcs_mask.bin.' iter_str], rbcs_mask_write, prec, ieee);
% Then files for each iteration 0..max_iter
for i = 0:max_iter
iter_str = sprintf('%010d', i);
write_binary(['rbcs_Tr1_fld.bin.' iter_str], rbcs_T_write, prec, ieee);
write_binary(['rbcs_mask.bin.' iter_str], rbcs_mask_write, prec, ieee);
end
%% ---------------- Quick 2D verification (and save) ----------------
fig1 = figure('Color','w','Name','Quick 2D checks','NumberTitle','off','Units','normalized','Position',[0.05 0.05 0.7 0.7]);
subplot(2,2,1); bathy_plot = bathy; bathy_plot(~lake_mask) = NaN; pcolor(XC,YC,bathy_plot); shading flat; axis equal tight; colorbar; title('Bathymetry (NaN outside lake)');
subplot(2,2,2); pcolor(XC,YC,double(moat2d)); shading flat; axis equal tight; colorbar; title('Moat Mask (2D)');
subplot(2,2,3); pcolor(XC,YC,squeeze(hFacC(:,:,1))); shading flat; axis equal tight; colorbar; title('Top Layer hFacC');
subplot(2,2,4); pcolor(XC,YC,rbcs_mask_write(:,:,1)); shading flat; axis equal tight; colorbar; title('RBC Mask Top Layer');
saveas(fig1,'quick_2D_checks.png');
% Also save the NaN-masked T/S for visual check
save('T_S_plot_versions.mat','T_NaN','S_NaN','-v7.3');
%% ---------------- Vertical X-Z and Y-Z slices (with RBC contour) ----------------
jmid = round(Ny/2); imid = round(Nx/2);
fig2 = figure('Color','w','Name','X-Z Vertical with RBC','NumberTitle','off','Units','normalized','Position',[0.05 0.1 0.4 0.6]);
imagesc(x, z_mid, squeeze(T_NaN(:, jmid, :))'); set(gca, 'YDir', 'normal'); hold on;
contour(x, z_mid, squeeze(rbcs_mask_write(:, jmid, :))', [0.5 0.5], 'r', 'LineWidth', 1.5);
xlabel('X (m)'); ylabel('Z (m)'); colorbar; title('X-Z Vertical with RBC (outside lake masked)');
yline(moat_depth, 'r--', 'Moat depth (m)', 'LabelHorizontalAlignment', 'left');
saveas(fig2,'XZ_RBC_vertical.png');
fig3 = figure('Color','w','Name','Y-Z Vertical with RBC','NumberTitle','off','Units','normalized','Position',[0.45 0.1 0.4 0.6]);
imagesc(y, z_mid, squeeze(T_NaN(imid, :, :))'); set(gca, 'YDir', 'normal'); hold on;
contour(y, z_mid, squeeze(rbcs_mask_write(imid, :, :))', [0.5 0.5], 'r', 'LineWidth', 1.5);
xlabel('Y (m)'); ylabel('Z (m)'); colorbar; title('Y-Z Vertical with RBC (outside lake masked)');
yline(moat_depth, 'r--', 'Moat depth (m)', 'LabelHorizontalAlignment', 'left');
saveas(fig3,'YZ_RBC_vertical.png');
%% ---------------- Radial cross-section over centerline (RBC layers + moat points) ----------------
radial_mask = jmid; % along +X centerline (index along y)
T_cut = squeeze(T_NaN(:, radial_mask, :))'; % Nz x Nx
rbcs_cut= squeeze(rbcs_mask_write(:, radial_mask, :))';
fig4 = figure('Color','w','Name','Radial cross-section (moat)','NumberTitle','off','Units','normalized','Position',[0.1 0.1 0.8 0.5]);
imagesc(x, z_mid, T_cut); set(gca, 'YDir', 'normal'); hold on;
contour(x, z_mid, rbcs_cut, [0.5 0.5], 'w', 'LineWidth', 2);
% mark moat 2D intersection on surface
moat_line = moat2d(:, radial_mask);
moat_x = x(moat_line);
plot(moat_x, 0.5 * ones(size(moat_x)), 'k.', 'MarkerSize', 6);
xlabel('X (m)'); ylabel('Z (m)'); colorbar; title('Radial X-Z: RBC layers (white) and moat points');
saveas(fig4,'radial_cross_section_moat.png');
%% ---------------- Combined 3D Bathymetry + RBC Temperature overlay ----------------
idx_lin = find(rbcs3d);
[i_r, j_c, k_lev] = ind2sub(size(rbcs3d), idx_lin);
lin2d = sub2ind([Ny Nx], i_r, j_c);
x_rbcs = XC(lin2d);
y_rbcs = YC(lin2d);
z_rbcs = z_mid(k_lev);
T_rbcs = T_write(idx_lin);
fig5 = figure('Color','w','Name','Combined 3D Bathymetry + RBC','NumberTitle','off','Units','normalized','Position',[0.05 0.05 0.8 0.8]); hold on;
surf(XC, YC, bathy, 'EdgeColor', 'none', 'FaceAlpha', 0.6); shading interp;
scatter3(x_rbcs, y_rbcs, z_rbcs, 40, T_rbcs, 'filled');
colorbar; caxis([T_surf T_target]); xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)');
title('Combined 3D Bathymetry + RBC Temperature (outside lake masked)');
view(45,30); axis equal; grid on;
% shoreline and moat inner edge
theta = linspace(0, 2*pi, 500);
x_shore = R_eff * cos(theta); y_shore = R_eff * sin(theta);
plot3(x_shore, y_shore, 0.5 * ones(size(theta)), 'k--', 'LineWidth', 1.2); % shoreline
x_mi = moat_inner * cos(theta); y_mi = moat_inner * sin(theta);
plot3(x_mi, y_mi, 0.5 * ones(size(theta)), 'r--', 'LineWidth', 1.2); % moat inner edge
saveas(fig5,'combined3D_bathy_RBC.png');
%% ---------------- Automatic 3D RBC-check (layers colored by temperature) ----------------
hFac_vals = hFacC(sub2ind([Ny Nx Nz], i_r, j_c, k_lev));
marker_sizes = max(6, 60 * reshape(hFac_vals, [], 1)); % size ~ local hFac thickness
fig6 = figure('Color','w','Name','Automatic RBC 3D check','NumberTitle','off','Units','normalized','Position',[0.05 0.05 0.8 0.8]); hold on;
surf(XC, YC, bathy, 'EdgeColor', 'none', 'FaceAlpha', 0.25); shading interp;
scatter3(x_rbcs, y_rbcs, z_rbcs, marker_sizes, T_rbcs, 'filled');
xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)');
title('RBC 3D overlay: colored by temperature (moat top layers only)');
colorbar; caxis([T_surf T_target]); view(45,30); axis equal; grid on;
saveas(fig6,'automatic_RBC_3D_check.png');
%% ---------------- Presentation 3D ----------------
fig7 = figure('Color','w','Name','Presentation 3D','NumberTitle','off','Units','normalized','Position',[0.05 0.05 0.8 0.8]); hold on;
surf(XC, YC, bathy, 'EdgeColor', 'none', 'FaceAlpha', 0.6); shading interp;
plot3(x_shore, y_shore, 0.5 * ones(size(theta)), 'k--', 'LineWidth', 1.2);
plot3(x_mi, y_mi, 0.5 * ones(size(theta)), 'r--', 'LineWidth', 1.2);
scatter3(x_rbcs, y_rbcs, z_rbcs, 40, T_rbcs, 'filled');
xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)');
title('3D Bathymetry + Moat (RBC) overlay — shore (black) / moat-inner (red)');
colorbar; caxis([T_surf T_target]); view(45,30); axis equal; grid on;
saveas(fig7,'presentation_3D_moat_RBC.png');
%% ---------------- Write summary text: which vertical levels are in RBC ----------------
rbcs_levels = find(squeeze(any(any(rbcs3d,1),2))); % indices of vertical levels included
fileID = fopen('rbcs_levels_summary.txt','w');
fprintf(fileID,'RBC Levels Summary\n');
fprintf(fileID,'===================\n');
fprintf(fileID,'Moat parameters: moat_width = %.1f m, moat_depth = %.2f m\n\n', moat_width, moat_depth);
fprintf(fileID,'Found %d vertical levels that include RBC (top->down):\n\n', numel(rbcs_levels));
fprintf(fileID,' k (idx) z_top (m) z_mid (m) z_bot (m) dZ (m)\n');
fprintf(fileID,'------------------------------------------------------\n');
for ii = 1:numel(rbcs_levels)
k = rbcs_levels(ii);
fprintf(fileID,'%5d %9.3f %9.3f %9.3f %8.3f\n', k, z_top(k), z_mid(k), z_bot(k), delZ(k));
end
fprintf(fileID,'\nTotal number of RBC 3D cells: %d\n', numel(idx_lin));
fclose(fileID);
fprintf('Wrote rbcs_levels_summary.txt\n');
%% ---------------- Short console summary ----------------
fprintf('\nShort console summary:\n');
fprintf('----------------------\n');
fprintf('Moat 2D cells (count) : %d\n', nnz(moat2d));
fprintf('RBC 3D cells (count) : %d\n', numel(idx_lin));
fprintf('Vertical levels with RBC : %s\n', mat2str(rbcs_levels'));
fprintf('RBC midpoint depths (m) : %s\n', mat2str(z_mid(rbcs_levels')));
fprintf('Binary files written (examples): T.init, S.init, hFacC.bin, delR.bin, rbcs_mask.bin, rbcs_Tr1_fld.bin\n');
fprintf('Verification PNGs written : quick_2D_checks.png, XZ_RBC_vertical.png, YZ_RBC_vertical.png, radial_cross_section_moat.png, combined3D_bathy_RBC.png, automatic_RBC_3D_check.png, presentation_3D_moat_RBC.png\n\n');
%% ---------------- Helper Function ----------------
function write_binary(fname, data, prec, ieee)
% Writes data to a binary file with the requested precision & endianness
fid = fopen(fname, 'w', ieee);
if fid == -1
error('Cannot open %s for writing', fname);
end
fwrite(fid, data, prec);
fclose(fid);
fprintf('Wrote %s\n', fname);
end
%% ---------------- MOC helper (MATLAB) ----------------
% USAGE: [MOC_bolus_layerwise, MOC_bolus_fromGM] = compute_MOC_from_GM_PsiY('GM_PsiY_file.mat', dxG, deepFacF)
% expects GM_PsiY(i,j,k) (same dims as model diagnostics), dxG(i,j) and deepFacF(k)
function [MOC_layer, MOC_fromGM, check] = compute_MOC_from_GM_PsiY(GM_PsiY, dxG, deepFacF)
% If GM_PsiY is a filename, load it
if ischar(GM_PsiY) || isstring(GM_PsiY)
data = load(GM_PsiY);
% try common variable names
if isfield(data,'GM_PsiY')
GM = data.GM_PsiY;
elseif isfield(data,'GM_Psi')
GM = data.GM_Psi;
else
error('Could not find GM_PsiY variable inside file');
end
else
GM = GM_PsiY; % already in workspace
end
% dims: (i,j,k). dxG should be (i,j). deepFacF should be (k,1)
[ni,nj,nk] = size(GM);
if ~isequal(size(dxG), [ni nj])
error('dxG must have dimensions [ni nj] equal to first two dims of GM');
end
if numel(deepFacF) ~= nk
error('deepFacF length must equal vertical levels nk');
end
% Layer-wise MOC (double sum over longitude i and sum from bottom to level k):
% MOC_bolus(j,k) = sum_i GM_PsiY(i,j,k) * dxG(i,j) * deepFacF(k)
MOC_fromGM = zeros(nj, nk);
for k = 1:nk
% Multiply GM at level k by dxG and deepFacF(k)
tmp = squeeze(GM(:,:,k)) .* dxG * deepFacF(k); % elementwise
MOC_fromGM(:,k) = sum(tmp,1); % sum over i -> 1 x nj
end
% Alternative approach: integrate -d/dz GM_PsiY vertically (layerwise)
% numerical vertical derivative: dGM_dz(k) ~ (GM(:,:,k) - GM(:,:,k+1)) / dz_k
% but we don't have dz here; instead we can compare cumulative sums
% Layerwise approach (sum GM over i then multiply by deepFacF to get volume transport)
MOC_layer = zeros(nj,nk);
for k = 1:nk
% integrate from bottom to k of GM(:,:,k') * dxG * deepFacF(k')
accum = zeros(1,nj);
for kk = 1:k
accum = accum + sum(squeeze(GM(:,:,kk)) .* dxG,1) * deepFacF(kk);
end
MOC_layer(:,k) = accum';
end
% Final check: compare two routes (they are related but may differ by sign/definition). Provide simple stats
check.max_abs_diff = max(abs(MOC_layer(:) - MOC_fromGM(:)));
check.mean_abs_diff = mean(abs(MOC_layer(:) - MOC_fromGM(:)));
end
-------------- next part --------------
A non-text attachment was scrubbed...
Name: vertical_section_temperature.jpg
Type: image/jpeg
Size: 126037 bytes
Desc: not available
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20250901/2c8b6f30/attachment-0001.jpg>
More information about the MITgcm-support
mailing list