clear;clc;close all;

load('data/Forcing_Comparison.mat');
load('data/Pantanal_shp.mat');
xv0 = ncread('data/domain.lnd.360x720_gswp3.0v1.c170606.nc','xv');
yv0 = ncread('data/domain.lnd.360x720_gswp3.0v1.c170606.nc','yv');
xc  = ncread('data/domain.lnd.360x720_gswp3.0v1.c170606.nc','xc');
yc  = ncread('data/domain.lnd.360x720_gswp3.0v1.c170606.nc','yc');
xc2  = ncread('data/domain.lnd.360x720_cruncep.130305.nc','xc');
yc2  = ncread('data/domain.lnd.360x720_cruncep.130305.nc','yc');

figure;
subplot(1,2,1);
plot(xc,xc2,'bo');
subplot(1,2,2);
plot(yc,yc2,'bo');


in = inpolygon(xc,yc,S.X+360,S.Y);
xv = NaN(4,sum(in(:)));
yv = NaN(4,sum(in(:)));
for i = 1 : 4
    tmp1 = xv0(:,:,i);
    tmp2 = yv0(:,:,i);
    xv(i,:) = tmp1(in);
    yv(i,:) = tmp2(in);
end
figure; set(gcf,'Position',[10 10 1200 1000]);
subplot(4,1,1);
a = nanmean(reshape(nanmean(Prec_GSWPv1,1)',[12 12]),2);
b = nanmean(reshape(nanmean(Prec_CRUNv7,1)',[12 12]),2);
plot(nanmean(reshape(nanmean(Prec_GSWPv1,1)',[12 12]),2),'b-','LineWidth',2);  hold on;
plot(nanmean(reshape(nanmean(Prec_CRUNv7,1)',[12 12]),2),'r--','LineWidth',2); grid on;
ylabel('Precipitation [mm/hour]','FontSize',16,'FontWeight','bold');
leg = legend('GSWP3v1','CRUNCEPv7');
leg.FontSize = 16;
leg.FontWeight = 'bold';

rd1 = (nanmean(a) - nanmean(b))/nanmean(a);

subplot(4,1,2);
a = nanmean(reshape(nanmean(Solr_GSWPv1,1)',[12 12]),2);
b = nanmean(reshape(nanmean(Solr_CRUNv7,1)',[12 12]),2);
plot(nanmean(reshape(nanmean(Solr_GSWPv1,1)',[12 12]),2),'b-','LineWidth',2);  hold on;
plot(nanmean(reshape(nanmean(Solr_CRUNv7,1)',[12 12]),2),'r--','LineWidth',2); grid on;
ylabel('Short wave incoming [W/m^{2}]','FontSize',16,'FontWeight','bold');

rd2 = (nanmean(a) - nanmean(b))/nanmean(a);

subplot(4,1,3);
plot(nanmean(reshape(nanmean(TBOT_GSWPv1,1)',[12 12]),2),'b-','LineWidth',2);  hold on;
plot(nanmean(reshape(nanmean(TBOT_CRUNv7,1)',[12 12]),2),'r--','LineWidth',2); grid on;
ylabel('Temperature [K]','FontSize',16,'FontWeight','bold');

subplot(4,1,4);
a = nanmean(reshape(nanmean(QBOT_GSWPv1,1)',[12 12]),2);
b = nanmean(reshape(nanmean(QBOT_CRUNv7,1)',[12 12]),2);
plot(a,'b-','LineWidth',2);  hold on;
plot(b,'r--','LineWidth',2); grid on;
ylabel('Humidity [-]','FontSize',16,'FontWeight','bold');

rd3 = (nanmean(a) - nanmean(b))/nanmean(a);

figure;set(gcf,'Position',[10 10 1200 300]);
subplot(1,3,1);
patch(xv,yv,nanmean(Prec_GSWPv1,2).*86400,'LineStyle','none'); colorbar;axis equal; caxis([0 5]);
subplot(1,3,2);
patch(xv,yv,nanmean(Prec_CRUNv7,2).*86400,'LineStyle','none'); colorbar;axis equal; caxis([0 5]);
subplot(1,3,3);
patch(xv,yv,nanmean(Prec_GSWPv1,2).*86400 - nanmean(Prec_CRUNv7,2).*86400,'LineStyle','none'); colorbar;axis equal;
caxis([-0.6 0.6]);colormap(gca,blue2red(11));

figure; set(gcf,'Position',[10 10 1200 300]);
subplot(1,3,1);
patch(xv,yv,nanmean(Solr_GSWPv1,2),'LineStyle','none'); colorbar; axis equal;
caxis([204 226]);
subplot(1,3,2);
patch(xv,yv,nanmean(Solr_CRUNv7,2),'LineStyle','none'); colorbar; axis equal;
caxis([204 226]);
subplot(1,3,3);
patch(xv,yv,nanmean(Solr_GSWPv1,2) - nanmean(Solr_CRUNv7,2),'LineStyle','none'); colorbar; axis equal;
caxis([-20 20]); colormap(gca,blue2red(11));

figure; set(gcf,'Position',[10 10 1200 300]);
subplot(1,3,1);
patch(xv,yv,nanmean(TBOT_GSWPv1,2),'LineStyle','none'); colorbar; axis equal;
%caxis([204 226]);
subplot(1,3,2);
patch(xv,yv,nanmean(TBOT_CRUNv7,2),'LineStyle','none'); colorbar; axis equal;
%caxis([204 226]);
subplot(1,3,3);
patch(xv,yv,nanmean(TBOT_GSWPv1,2) - nanmean(TBOT_CRUNv7,2),'LineStyle','none'); colorbar; axis equal;
%caxis([-20 20]); colormap(gca,blue2red(11));

figure; set(gcf,'Position',[10 10 1200 300]);
subplot(1,3,1);
patch(xv,yv,nanmean(QBOT_GSWPv1,2),'LineStyle','none'); colorbar; axis equal;
caxis([0.01 0.02]);
subplot(1,3,2);
patch(xv,yv,nanmean(QBOT_CRUNv7,2),'LineStyle','none'); colorbar; axis equal;
caxis([0.01 0.02]);
subplot(1,3,3);
patch(xv,yv,nanmean(QBOT_GSWPv1,2) - nanmean(QBOT_CRUNv7,2),'LineStyle','none'); colorbar; axis equal;
caxis([-0.003 0.003]); colormap(gca,blue2red(11));