% Define paths and variables
datapath = './';
varnames = {'AODVIS', 'colccn.3'};
nvar = length(varnames);
filenames = {'v2_ndg_cdnc_ssat_diag_simple_pi.eam.h0.2011-01.nc', ...
             'v2_ndg_cdnc_ssat_diag_simple_pd.eam.h0.2011-01.nc'};
nfile = length(filenames);
nPlotColumn = 5;

% Create figure
fig = figure('Position', [100, 100, 1200, 800]);

% Define fill value for avoiding division by zero
FillValue = -999;

% Get lat/lon info and check files for consistent dimensions
for ifile = 1:nfile
    ncfile = fullfile(datapath, filenames{ifile});
    if ifile == 1
        area = ncread(ncfile, 'area');
        ncol = length(area);
        lon = ncread(ncfile, 'lon');
        lat = ncread(ncfile, 'lat');
    else
        area_check = ncread(ncfile, 'area');
        ncol1 = length(area_check);
        if ncol1 ~= ncol
            error('Error: two files have different numbers of grid columns. Abort.');
        end
    end
end

% Process each variable
for ivar = 1:nvar
    array4plotting = zeros(nPlotColumn, ncol);
    
    % Read data from files
    for ifile = 1:nfile
        ncfile = fullfile(datapath, filenames{ifile});
        array4plotting(ifile,:) = ncread(ncfile, varnames{ivar});
    end
    
    % Calculate differences and relative differences
    % Absolute difference
    array4plotting(3,:) = array4plotting(2,:) - array4plotting(1,:);
    
    % Relative difference: (test - ctrl)/ctrl
    tmp = array4plotting(1,:);
    denom = tmp;
    denom(tmp == 0) = FillValue;
    array4plotting(4,:) = (array4plotting(2,:) - array4plotting(1,:))./denom;
    array4plotting(4, denom == FillValue) = FillValue;
    
    % Relative difference: (test - ctrl)*2/(test + ctrl)
    tmp = 0.5 * (array4plotting(1,:) + array4plotting(2,:));
    denom = tmp;
    denom(tmp == 0) = FillValue;
    array4plotting(5,:) = (array4plotting(2,:) - array4plotting(1,:))./denom;
    array4plotting(5, denom == FillValue) = FillValue;
    
    % Create plots for this variable
    for ic = 1:nPlotColumn
        subplot_idx = (ivar-1)*nPlotColumn + ic;
        subplot(nvar, nPlotColumn, subplot_idx);
        
        % Create trimesh/trisurf for unstructured grid
        % We need to create a Delaunay triangulation for unstructured grid
        tri = delaunayTriangulation(lon, lat);
        
        switch ic
            case 1
                title_str = [varnames{ivar} ' - Ctrl'];
                cmap = parula;
            case 2
                title_str = [varnames{ivar} ' - Test'];
                cmap = parula;
            case 3
                title_str = 'Difference: test - ctrl';
                cmap = redblue(64); % Custom red-blue colormap (see function below)
                
                % Symmetric color limits
                clim_val = max(abs([min(array4plotting(ic,:)), max(array4plotting(ic,:))]));
                if ~isnan(clim_val) && clim_val > 0
                    caxis([-clim_val, clim_val]);
                end
            case 4
                title_str = 'Rel diff: (test - ctrl)/ctrl';
                cmap = redblue(64);
                
                % Symmetric color limits
                valid_data = array4plotting(ic,:) ~= FillValue;
                if any(valid_data)
                    clim_val = max(abs([min(array4plotting(ic,valid_data)), max(array4plotting(ic,valid_data))]));
                    if ~isnan(clim_val) && clim_val > 0
                        caxis([-clim_val, clim_val]);
                    end
                end
            case 5
                title_str = 'Rel diff: (test - ctrl)*2/(test + ctrl)';
                cmap = redblue(64);
                
                % Symmetric color limits
                valid_data = array4plotting(ic,:) ~= FillValue;
                if any(valid_data)
                    clim_val = max(abs([min(array4plotting(ic,valid_data)), max(array4plotting(ic,valid_data))]));
                    if ~isnan(clim_val) && clim_val > 0
                        caxis([-clim_val, clim_val]);
                    end
                end
        end
        
        % Plot the data on the unstructured grid
        trisurf(tri, lon, lat, zeros(size(lon)), array4plotting(ic,:), 'EdgeColor', 'none');
        view(2);  % 2D view
        colormap(cmap);
        colorbar;
        
        % Map appearance
        axis equal tight;
        title(title_str);
        xlabel('Longitude');
        ylabel('Latitude');
        
        % Set world map boundaries
        xlim([-180, 180]);
        ylim([-90, 90]);
        
        % Add grid lines
        grid on;
        
        % Add coastlines using mapping toolbox if available
        if license('test', 'Mapping_Toolbox')
            hold on;
            coastlines = load('coastlines');
            plot(coastlines.coastlon, coastlines.coastlat, 'k-', 'LineWidth', 0.5);
            hold off;
        end
    end
end

% Adjust layout
tight_layout(fig);

% Save figure
sa% Define paths and variables
datapath = './';
varnames = {'AODVIS', 'colccn.3'};
nvar = length(varnames);
filenames = {'v2_ndg_cdnc_ssat_diag_simple_pi.eam.h0.2011-01.nc', ...
             'v2_ndg_cdnc_ssat_diag_simple_pd.eam.h0.2011-01.nc'};
nfile = length(filenames);
nPlotColumn = 5;

% Create figure
fig = figure('Position', [100, 100, 1200, 800]);

% Define fill value for avoiding division by zero
FillValue = -999;

% Get lat/lon info and check files for consistent dimensions
for ifile = 1:nfile
    ncfile = fullfile(datapath, filenames{ifile});
    if ifile == 1
        area = ncread(ncfile, 'area');
        ncol = length(area);
        lon = ncread(ncfile, 'lon');
        lat = ncread(ncfile, 'lat');
    else
        area_check = ncread(ncfile, 'area');
        ncol1 = length(area_check);
        if ncol1 ~= ncol
            error('Error: two files have different numbers of grid columns. Abort.');
        end
    end
end

% Process each variable
for ivar = 1:nvar
    array4plotting = zeros(nPlotColumn, ncol);
    
    % Read data from files
    for ifile = 1:nfile
        ncfile = fullfile(datapath, filenames{ifile});
        array4plotting(ifile,:) = ncread(ncfile, varnames{ivar});
    end
    
    % Calculate differences and relative differences
    % Absolute difference
    array4plotting(3,:) = array4plotting(2,:) - array4plotting(1,:);
    
    % Relative difference: (test - ctrl)/ctrl
    tmp = array4plotting(1,:);
    denom = tmp;
    denom(tmp == 0) = FillValue;
    array4plotting(4,:) = (array4plotting(2,:) - array4plotting(1,:))./denom;
    array4plotting(4, denom == FillValue) = FillValue;
    
    % Relative difference: (test - ctrl)*2/(test + ctrl)
    tmp = 0.5 * (array4plotting(1,:) + array4plotting(2,:));
    denom = tmp;
    denom(tmp == 0) = FillValue;
    array4plotting(5,:) = (array4plotting(2,:) - array4plotting(1,:))./denom;
    array4plotting(5, denom == FillValue) = FillValue;
    
    % Create plots for this variable
    for ic = 1:nPlotColumn
        subplot_idx = (ivar-1)*nPlotColumn + ic;
        subplot(nvar, nPlotColumn, subplot_idx);
        
        % Create trimesh/trisurf for unstructured grid
        % We need to create a Delaunay triangulation for unstructured grid
        tri = delaunayTriangulation(lon, lat);
        
        switch ic
            case 1
                title_str = [varnames{ivar} ' - Ctrl'];
                cmap = parula;
            case 2
                title_str = [varnames{ivar} ' - Test'];
                cmap = parula;
            case 3
                title_str = 'Difference: test - ctrl';
                cmap = redblue(64); % Custom red-blue colormap (see function below)
                
                % Symmetric color limits
                clim_val = max(abs([min(array4plotting(ic,:)), max(array4plotting(ic,:))]));
                if ~isnan(clim_val) && clim_val > 0
                    caxis([-clim_val, clim_val]);
                end
            case 4
                title_str = 'Rel diff: (test - ctrl)/ctrl';
                cmap = redblue(64);
                
                % Symmetric color limits
                valid_data = array4plotting(ic,:) ~= FillValue;
                if any(valid_data)
                    clim_val = max(abs([min(array4plotting(ic,valid_data)), max(array4plotting(ic,valid_data))]));
                    if ~isnan(clim_val) && clim_val > 0
                        caxis([-clim_val, clim_val]);
                    end
                end
            case 5
                title_str = 'Rel diff: (test - ctrl)*2/(test + ctrl)';
                cmap = redblue(64);
                
                % Symmetric color limits
                valid_data = array4plotting(ic,:) ~= FillValue;
                if any(valid_data)
                    clim_val = max(abs([min(array4plotting(ic,valid_data)), max(array4plotting(ic,valid_data))]));
                    if ~isnan(clim_val) && clim_val > 0
                        caxis([-clim_val, clim_val]);
                    end
                end
        end
        
        % Plot the data on the unstructured grid
        trisurf(tri, lon, lat, zeros(size(lon)), array4plotting(ic,:), 'EdgeColor', 'none');
        view(2);  % 2D view
        colormap(cmap);
        colorbar;
        
        % Map appearance
        axis equal tight;
        title(title_str);
        xlabel('Longitude');
        ylabel('Latitude');
        
        % Set world map boundaries
        xlim([-180, 180]);
        ylim([-90, 90]);
        
        % Add grid lines
        grid on;
        
        % Add coastlines using mapping toolbox if available
        if license('test', 'Mapping_Toolbox')
            hold on;
            coastlines = load('coastlines');
            plot(coastlines.coastlon, coastlines.coastlat, 'k-', 'LineWidth', 0.5);
            hold off;
        end
    end
end

% Adjust layout
tight_layout(fig);

% Save figure
saveas(fig, './fig_compare_two_runs_MATLAB.pdf');

function cmap = redblue(m)
    % Custom red-blue colormap similar to NCL's "BlueDarkRed18"
    if nargin < 1
        m = 64;
    end
    
    % Create a colormap from blue to white to red
    n = ceil(m/2);
    r = [linspace(0, 1, n); ones(n, 1)']';
    g = [linspace(0, 1, n); linspace(1, 0, n)]';
    b = [ones(n, 1)'; linspace(1, 0, n)]';
    
    cmap = [r g b];
    
    % Ensure the colormap has exactly m rows
    if size(cmap, 1) > m
        cmap = cmap(1:m, :);
    end
end

function tight_layout(fig)
    % Simple function to improve subplot spacing
    p = get(fig, 'Position');
    set(fig, 'Position', [p(1) p(2) p(3) p(4)]);
    set(fig, 'PaperPositionMode', 'auto');
    
    % Adjust subplot spacing
    hspace = 0.1;  % Horizontal space between subplots
    vspace = 0.15; % Vertical space between subplots
    
    axs = findall(fig, 'type', 'axes');
    for i = 1:length(axs)
        pos = get(axs(i), 'Position');
        set(axs(i), 'Position', [pos(1) pos(2) pos(3)*(1-hspace) pos(4)*(1-vspace)]);
    end
en% Define paths and variables
datapath = './';
varnames = {'AODVIS', 'colccn.3'};
nvar = length(varnames);
filenames = {'v2_ndg_cdnc_ssat_diag_simple_pi.eam.h0.2011-01.nc', ...
             'v2_ndg_cdnc_ssat_diag_simple_pd.eam.h0.2011-01.nc'};
nfile = length(filenames);
nPlotColumn = 5;

% Create figure
fig = figure('Position', [100, 100, 1200, 800]);

% Define fill value for avoiding division by zero
FillValue = -999;

% Get lat/lon info and check files for consistent dimensions
for ifile = 1:nfile
    ncfile = fullfile(datapath, filenames{ifile});
    if ifile == 1
        area = ncread(ncfile, 'area');
        ncol = length(area);
        lon = ncread(ncfile, 'lon');
        lat = ncread(ncfile, 'lat');
    else
        area_check = ncread(ncfile, 'area');
        ncol1 = length(area_check);
        if ncol1 ~= ncol
            error('Error: two files have different numbers of grid columns. Abort.');
        end
    end
end

% Process each variable
for ivar = 1:nvar
    array4plotting = zeros(nPlotColumn, ncol);
    
    % Read data from files
    for ifile = 1:nfile
        ncfile = fullfile(datapath, filenames{ifile});
        array4plotting(ifile,:) = ncread(ncfile, varnames{ivar});
    end
    
    % Calculate differences and relative differences
    % Absolute difference
    array4plotting(3,:) = array4plotting(2,:) - array4plotting(1,:);
    
    % Relative difference: (test - ctrl)/ctrl
    tmp = array4plotting(1,:);
    denom = tmp;
    denom(tmp == 0) = FillValue;
    array4plotting(4,:) = (array4plotting(2,:) - array4plotting(1,:))./denom;
    array4plotting(4, denom == FillValue) = FillValue;
    
    % Relative difference: (test - ctrl)*2/(test + ctrl)
    tmp = 0.5 * (array4plotting(1,:) + array4plotting(2,:));
    denom = tmp;
    denom(tmp == 0) = FillValue;
    array4plotting(5,:) = (array4plotting(2,:) - array4plotting(1,:))./denom;
    array4plotting(5, denom == FillValue) = FillValue;
    
    % Create plots for this variable
    for ic = 1:nPlotColumn
        subplot_idx = (ivar-1)*nPlotColumn + ic;
        subplot(nvar, nPlotColumn, subplot_idx);
        
        % Create trimesh/trisurf for unstructured grid
        % We need to create a Delaunay triangulation for unstructured grid
        tri = delaunayTriangulation(lon, lat);
        
        switch ic
            case 1
                title_str = [varnames{ivar} ' - Ctrl'];
                cmap = parula;
            case 2
                title_str = [varnames{ivar} ' - Test'];
                cmap = parula;
            case 3
                title_str = 'Difference: test - ctrl';
                cmap = redblue(64); % Custom red-blue colormap (see function below)
                
                % Symmetric color limits
                clim_val = max(abs([min(array4plotting(ic,:)), max(array4plotting(ic,:))]));
                if ~isnan(clim_val) && clim_val > 0
                    caxis([-clim_val, clim_val]);
                end
            case 4
                title_str = 'Rel diff: (test - ctrl)/ctrl';
                cmap = redblue(64);
                
                % Symmetric color limits
                valid_data = array4plotting(ic,:) ~= FillValue;
                if any(valid_data)
                    clim_val = max(abs([min(array4plotting(ic,valid_data)), max(array4plotting(ic,valid_data))]));
                    if ~isnan(clim_val) && clim_val > 0
                        caxis([-clim_val, clim_val]);
                    end
                end
            case 5
                title_str = 'Rel diff: (test - ctrl)*2/(test + ctrl)';
                cmap = redblue(64);
                
                % Symmetric color limits
                valid_data = array4plotting(ic,:) ~= FillValue;
                if any(valid_data)
                    clim_val = max(abs([min(array4plotting(ic,valid_data)), max(array4plotting(ic,valid_data))]));
                    if ~isnan(clim_val) && clim_val > 0
                        caxis([-clim_val, clim_val]);
                    end
                end
        end
        
        % Plot the data on the unstructured grid
        trisurf(tri, lon, lat, zeros(size(lon)), array4plotting(ic,:), 'EdgeColor', 'none');
        view(2);  % 2D view
        colormap(cmap);
        colorbar;
        
        % Map appearance
        axis equal tight;
        title(title_str);
        xlabel('Longitude');
        ylabel('Latitude');
        
        % Set world map boundaries
        xlim([-180, 180]);
        ylim([-90, 90]);
        
        % Add grid lines
        grid on;
        
        % Add coastlines using mapping toolbox if available
        if license('test', 'Mapping_Toolbox')
            hold on;
            coastlines = load('coastlines');
            plot(coastlines.coastlon, coastlines.coastlat, 'k-', 'LineWidth', 0.5);
            hold off;
        end
    end
end

% Adjust layout
tight_layout(fig);

% Save figure
saveas(fig, './fig_compare_two_runs_MATLAB.pdf');

function cmap = redblue(m)
    % Custom red-blue colormap similar to NCL's "BlueDarkRed18"
    if nargin < 1
        m = 64;
    end
    
    % Create a colormap from blue to white to red
    n = ceil(m/2);
    r = [linspace(0, 1, n); ones(n, 1)']';
    g = [linspace(0, 1, n); linspace(1, 0, n)]';
    b = [ones(n, 1)'; linspace(1, 0, n)]';
    
    cmap = [r g b];
    
    % Ensure the colormap has exactly m rows
    if size(cmap, 1) > m
        cmap = cmap(1:m, :);
    end
end

function tight_layout(fig)
    % Simple function to improve subplot spacing
    p = get(fig, 'Position');
    set(fig, 'Position', [p(1) p(2) p(3) p(4)]);
    set(fig, 'PaperPositionMode', 'auto');
    
    % Adjust subplot spacing
    hspace = 0.1;  % Horizontal space between subplots
    vspace = 0.15; % Vertical space between subplots
    
    axs = findall(fig, 'type', 'axes');
    for i = 1:length(axs)
        pos = get(axs(i), 'Position');
        set(axs(i), 'Position', [pos(1) pos(2) pos(3)*(1-hspace) pos(4)*(1-vspace)]);
    end
enddveas(fig, './fig_compare_two_runs_MATLAB.pdf');

function cmap = redblue(m)
    % Custom red-blue colormap similar to NCL's "BlueDarkRed18"
    if nargin < 1
        m = 64;
    end
    
    % Create a colormap from blue to white to red
    n = ceil(m/2);
    r = [linspace(0, 1, n); ones(n, 1)']';
    g = [linspace(0, 1, n); linspace(1, 0, n)]';
    b = [ones(n, 1)'; linspace(1, 0, n)]';
    
    cmap = [r g b];
    
    % Ensure the colormap has exactly m rows
    if size(cmap, 1) > m
        cmap = cmap(1:m, :);
    end
end

function tight_layout(fig)
    % Simple function to improve subplot spacing
    p = get(fig, 'Position');
    set(fig, 'Position', [p(1) p(2) p(3) p(4)]);
    set(fig, 'PaperPositionMode', 'auto');
    
    % Adjust subplot spacing
    hspace = 0.1;  % Horizontal space between subplots
    vspace = 0.15; % Vertical space between subplots
    
    axs = findall(fig, 'type', 'axes');
    for i = 1:length(axs)
        pos = get(axs(i), 'Position');
        set(axs(i), 'Position', [pos(1) pos(2) pos(3)*(1-hspace) pos(4)*(1-vspace)]);
    end
end