Matlab script for comparing SOURCE load scenarios

clear all
close all

figure

set(gcf,'Position',[195 19 901 697])

file2 = 'http://dapds00.nci.org.au/thredds/dodsC/fx3/gbr4_bgc_GBR4_H2p0_B2p0_Cpre_Dcrt/gbr4_bgc_simple_2011-02.nc'
file1 = 'http://dapds00.nci.org.au/thredds/dodsC/fx3/gbr4_bgc_GBR4_H2p0_B2p0_Chyd_Dcrt/gbr4_bgc_simple_2011-02.nc'

filename = file1;

y_centre = ncread(filename,'latitude');
x_centre = ncread(filename,'longitude');

ylim = [-21 -18];
xlim = [146 150];

botz = ncread(filename,'botz');

vars = {'Chl_a_sum','DIN','EFI'}
varsname = {'Total Chl \it a \rm [mg/m3]','DIN [mg/m3]','TSS [mg/m3]'};
varsclim = [2 14 8];

sims = {'SOURCE Baseline','SOUCRE Pre-Ind.','Baseline - Pre-Ind.'}

tt = 1; % time index
ll = 45; % layer index

time = double(ncread(file1,'time',tt,1,1))+datenum(1990,1,1);

rows=3;
cols=3;
vspc=0.01;
hspc=0.005;
left=0.02;
right = 0.02;
bot=0.02;
top = 0.05;
wdth= (1.0 - left - right-(cols-1)*hspc)/cols;
hgt= (1.0 - bot -top -(rows-1)*vspc)/rows;

for row = 1:3
var1 = double(ncread(file1,char(vars(row)),[1 1 ll tt],[inf inf 1 1],[1 1 1 1]));
var2 = double(ncread(file2,char(vars(row)),[1 1 ll tt],[inf inf 1 1],[1 1 1 1]));
if strcmp(char(vars(row)),'EFI')
var1 = var1*1000;
var2 = var2*1000;
end

load ../misc/eastcoast.dat
figure(1)

col = 1;
h = axes('position',[left+(col-1)*(wdth+hspc) (1.0-top-row*(hgt+vspc)) wdth hgt]);

pcolor(x_centre,y_centre,var1)

axis equal % Sets the aspect ratio so that the data units are the same in
% every direction.
shading flat % Removes the black lines that are a default of pcolor.
pos = get(gca,'Position');
c=colorbar('horiz');set(gca,'clim',[0 varsclim(row)]);
set(c,'Position',[left+(col-1)*(wdth+hspc)+0.03 (1.0-top-row*(hgt+vspc))+0.05 0.1 0.02]);
set(get(c,'Title'),'String',char(varsname(row)));
set(gca,'xlim',xlim,'ylim',ylim);
set(gca,'xticklabel','','yticklabel','');
set(gca,'Position',pos);
if row==1
title(char(sims(1)));
end

text(xlim(1) + 0.4*diff(xlim),ylim(1)+0.05*diff(ylim),datestr(time,'dd-mm-yyyy'));

hold on

plot(eastcoast(:,1),eastcoast(:,2),'-k')

col = 2;
h = axes('position',[left+(col-1)*(wdth+hspc) (1.0-top-row*(hgt+vspc)) wdth hgt]);

var2(var2<1e-16) = NaN;
pcolor(x_centre,y_centre,var2)

axis equal % Sets the aspect ratio so that the data units are the same in
% every direction.
shading flat % Removes the black lines that are a default of pcolor.
pos = get(gca,'Position');
c=colorbar('horiz');set(gca,'clim',[0 varsclim(row)]);
set(c,'Position',[left+(col-1)*(wdth+hspc)+0.03 (1.0-top-row*(hgt+vspc))+0.05 0.1 0.02]);
set(get(c,'Title'),'String',char(varsname(row)));
set(gca,'xlim',xlim,'ylim',ylim);
set(gca,'xticklabel','','yticklabel','');
hold on

text(xlim(1) + 0.4*diff(xlim),ylim(1)+0.05*diff(ylim),datestr(time,'dd-mm-yyyy'));

if row ==1
title(char(sims(2)));
end
plot(eastcoast(:,1),eastcoast(:,2),'-k')

col = 3;
h = axes('position',[left+(col-1)*(wdth+hspc) (1.0-top-row*(hgt+vspc)) wdth hgt]);

pcolor(x_centre,y_centre,var1-var2)

axis equal % Sets the aspect ratio so that the data units are the same in
% every direction.
shading flat % Removes the black lines that are a default of pcolor.
pos = get(gca,'Position');
c=colorbar('horiz');set(gca,'clim',[-varsclim(row) varsclim(row)]/10);
set(c,'Position',[left+(col-1)*(wdth+hspc)+0.03 (1.0-top-row*(hgt+vspc))+0.05 0.1 0.02]);
set(get(c,'Title'),'String',char(varsname(row)));
set(gca,'xlim',xlim,'ylim',ylim);
set(gca,'xticklabel','','yticklabel','');
hold on
text(xlim(1) + 0.4*diff(xlim),ylim(1)+0.05*diff(ylim),datestr(time,'dd-mm-yyyy'));
plot(eastcoast(:,1),eastcoast(:,2),'-k')
if (row==1)
title(char(sims(3)));
end
end

print('eReefs_scenario_analysis_seagrass','-dpng','-r600')