clc
clear all
close all

cn51_err_read_path = '/discover/nobackup/jkolassa/CN51_config_debug/carbon_error_files/';
cn51_err_old_read_path = '/discover/nobackup/jkolassa/CN51_config_debug/carbon_error_files_curr/';

%-----------------------------------
% read location data
%-----------------------------------

cn51_data_file = '/discover/nobackup/jkolassa/CN51_config_debug/GEOSldas_CN51_local_280/output/SMAP_EASEv2_M36/cat/ens0000/Y1980/M01/GEOSldas_CN51_local_280.tavg24_1d_lnd_Nt.monthly.198001.nc4';

cn51_lat = ncread(cn51_data_file,'lat');
cn51_lon = ncread(cn51_data_file,'lon');

%----------------------
% read CN51 column data
%----------------------

% carbon imbalance
cn51_errcb_col_file = [cn51_err_read_path 'carbon_error_col_year1'];
fid = fopen(cn51_errcb_col_file,'r');
data = fscanf(fid,'%i %f');
fclose(fid);

cn51_errcb_col = reshape(data,[2,length(data)/2]);
cn51_errcb_col = cn51_errcb_col';


cn51_errcb_old_col_file = [cn51_err_old_read_path 'carbon_error_col_year1'];
fid = fopen(cn51_errcb_old_col_file,'r');
data = fscanf(fid,'%i %f');
fclose(fid);

cn51_errcb_old_col = reshape(data,[2,length(data)/2]);
cn51_errcb_old_col = cn51_errcb_old_col';

%-------------------
% transform matrices
%-------------------

for i = 1:max(cn51_errcb_col(:,1))

    col_index = find(cn51_errcb_col(:,1)==i);
    cn51_errcb_col_mat(i,:) = cn51_errcb_col(col_index,2);
    cn51_max_errcb(i) = max(abs(cn51_errcb_col_mat(i,:)));
    cn51_mean_errcb(i) = mean(cn51_errcb_col_mat(i,:));

    col_index_2 = find(cn51_errcb_old_col(:,1)==i);
    cn51_errcb_old_col_mat(i,:) = cn51_errcb_old_col(col_index_2,2);
    cn51_max_errcb_old(i) = max(abs(cn51_errcb_old_col_mat(i,:)));
    cn51_mean_errcb_old(i) = mean(cn51_errcb_old_col_mat(i,:));
end

max_plot_mat = [cn51_max_errcb' cn51_max_errcb_old'];

figure
bar(max_plot_mat)
legend('w/ fix','w/ no fix')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15]) 
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/balance/max_carbon_balance'],'png')

figure
plot(cn51_max_errcb,'*b')
hold on
plot(cn51_max_errcb_old,'*r')
legend('w/ fix','w/ no fix')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/balance/max_carbon_balance_scatter'],'png')

figure
plot(cn51_mean_errcb,'*b')
hold on
plot(cn51_mean_errcb_old,'*r')
legend('w/ fix','w/ no fix')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/balance/mean_carbon_balance_scatter'],'png')
