clear;close all;clc;

% This script is to validate other variables

addpath('/Users/xudo627/donghui/CODE/dengwirda-inpoly-355c57c/');

xv=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','xv');
yv=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','yv');
xc=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','xc');
yc=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','yc');


files1 = dir('/Volumes/LaCie/PROJECTS/Pantanal/outputs_calibration_CLMCRUNCEPv7/Pantanal_*.elm.*.nc');
files0 = dir('/Volumes/LaCie/PROJECTS/Pantanal/outputs_default_CLMCRUNCEPv7/Pantanal_*.elm.*.nc');
evaluate_zwt = 1;
evaluate_sm  = 0;
evaluate_et  = 0;

if evaluate_et
    
    [lon_gsim,lat_gsim,Sh,yr,mo,da,mu_gsim,sd,cv] = get_GSIM_discharge('BR_0002835',1);
    
    for i = 1 : length(files1)
        filename = fullfile(files1(i).folder,files1(i).name);
        filename0 = fullfile(files0(i).folder,files0(i).name);
        disp(['i = ' num2str(i) ': ' filename]);
        if i == 1
            et  = ncread(filename,'EFLX_LH_TOT');
            ro  = ncread(filename,'QRUNOFF');
            tws = ncread(filename,'TWS');
            et0 = ncread(filename0,'EFLX_LH_TOT');
            et0 = et0(1:191);
            ro0 = ncread(filename0,'QRUNOFF');
            ro0 = ro0(1:191);
            tws0 = ncread(filename0,'TWS');
            tws0 = tws0(1:191);
        else
            et = cat(2,et,ncread(filename,'EFLX_LH_TOT'));
            tmp = ncread(filename0,'EFLX_LH_TOT');
            et0 = cat(2,et0,tmp(1:191));
            ro = cat(2,ro,ncread(filename,'QRUNOFF'));
            tmp = ncread(filename0,'QRUNOFF');
            ro0 = cat(2,ro0,tmp(1:191));
            tws = cat(2,tws,ncread(filename,'TWS'));
            tmp = ncread(filename0,'TWS');
            tws0 = cat(2,tws0,tmp(1:191));
        end
    end
    %tws(tws < 0) = 0;
    
    figure; set(gcf,'Position',[10 10 1200 1000]);
    subplot(2,2,1);
    patch(xv,yv,nanmean(et,2),'LineStyle','none');  colorbar; caxis([0 110]); colormap(gca,jet(100));
    subplot(2,2,2);
    patch(xv,yv,nanmean(et0,2),'LineStyle','none'); colorbar; caxis([0 110]); colormap(gca,jet(100));
    subplot(2,2,[3 4]);
    plot(1:384,nanmean(et,1),'b-','LineWidth',2); hold on;
    plot(1:384,nanmean(et0,1),'r-','LineWidth',2);
    
    figure;
    plot(nanmean(reshape(nanmean(et,1)',[12 32]),1),'b-','LineWidth',2); hold on;
    plot(nanmean(reshape(nanmean(et0,1)',[12 32]),1),'r--','LineWidth',2);
    
    figure;
    plot(nanmean(reshape(nanmean(ro,1)',[12 32]),1),'b-','LineWidth',2); hold on;
    plot(nanmean(reshape(nanmean(ro0,1)',[12 32]),1),'r--','LineWidth',2);
    
    baseline  = nanmean(tws(:,(2004-1979)*12+1 : (2009-1979+1)*12),2);
    baseline0 = nanmean(tws0(:,(2004-1979)*12+1 : (2009-1979+1)*12),2);
    
    lon = ncread('~/Downloads/twsa_0.5x0.5.nc','lon');
    lat = ncread('~/Downloads/twsa_0.5x0.5.nc','lat');
    twsa= ncread('~/Downloads/twsa_0.5x0.5.nc','twsa');
    
    [lon,lat] = meshgrid(lon,lat);
    lon = lon';
    lat = lat';
    
    in = inpoly2([lon(:) lat(:)],[Sh.X' Sh.Y']);
    
    twsa_pantanal = NaN(nansum(in(:)),size(twsa,3));
    for i = 1 : size(twsa,3)
        tmp = twsa(:,:,i);
        twsa_pantanal(:,i) = tmp(in);
    end
    
    figure;
    subplot(2,1,1);
    plot(nanmean(tws,1)-nanmean(baseline),'b-','LineWidth',2); hold on;
    plot(nanmean(tws0,1)-nanmean(baseline0),'r--','LineWidth',2);
    subplot(2,1,2);
    plot(1979:2010,nanmean(reshape(nanmean(tws-baseline,1)',[12 32]),1),'b-','LineWidth',2); hold on;
    plot(1979:2010,nanmean(reshape(nanmean(tws0-baseline0,1)',[12 32]),1),'r--','LineWidth',2);
    plot(2002:2014,nanmean(reshape(nanmean(twsa_pantanal,1)',[12 13]),1),'k--','LineWidth',2);
   
    figure;
    subplot(1,2,1);
    patch(xv,yv,nanmean(tws,2),'LineStyle','none'); colorbar;
    subplot(1,2,2);
    patch(xv,yv,nanmean(tws0,2),'LineStyle','none'); colorbar;
end

if evaluate_sm
    
    [lon_gsim,lat_gsim,Sh,yr,mo,da,mu_gsim,sd,cv] = get_GSIM_discharge('BR_0002835',1);
    
    for i = 1 : length(files1)
        filename = fullfile(files1(i).folder,files1(i).name);
        filename0 = fullfile(files0(i).folder,files0(i).name);
        disp(['i = ' num2str(i) ': ' filename]);
        if i == 1
            sm  = ncread(filename,'SOILWATER_10CM')./100;
            sm0 = ncread(filename0,'SOILWATER_10CM')./100;
            sm0 = sm0(1:191);
        else
            sm = cat(2,sm,ncread(filename,'SOILWATER_10CM')./100);
            tmp = ncread(filename0,'SOILWATER_10CM')./100;
            sm0 = cat(2,sm0,tmp(1:191));
        end
    end
    SMsurf = ncread('/Volumes/LaCie/DATA/SMsurf_1980-2020_GLEAM_v3.5a_MO.nc','SMsurf');
    lon = ncread('/Volumes/LaCie/DATA/SMsurf_1980-2020_GLEAM_v3.5a_MO.nc','lon');
    lat = ncread('/Volumes/LaCie/DATA/SMsurf_1980-2020_GLEAM_v3.5a_MO.nc','lat');
    [lon, lat] = meshgrid(lon,lat);
%     lon = lon'; lat = lat';
    in = inpoly2([lon(:) lat(:)],[Sh.X' Sh.Y']);
    
    SMsurf   = SMsurf(:,:,1:372);
    sm_gleam = NaN(nansum(in),size(SMsurf,3));
    for i = 1 : size(SMsurf,3)
        tmp = SMsurf(:,:,i);
        sm_gleam(:,i) = tmp(in);
    end
    xc2 = lon(in); yc2 = lat(in);
    xv2 = NaN(4,length(xc2)); yv2 = NaN(4,length(yc2));
    xv2(1,:) = xc2 - 0.25/2;xv2(2,:) = xc2 + 0.25/2;xv2(3,:) = xc2 + 0.25/2;xv2(4,:) = xc2 - 0.25/2;
    yv2(1,:) = yc2 - 0.25/2;yv2(2,:) = yc2 - 0.25/2;yv2(3,:) = yc2 + 0.25/2;yv2(4,:) = yc2 + 0.25/2;
    
    figure; set(gcf,'Position',[10 10 1200 1000]);
    subplot(2,2,1);
    patch(xv,yv,nanmean(sm,2),'LineStyle','none'); colorbar; %caxis([0 60]); colormap(gca,jet(10));
    subplot(2,2,2);
    patch(xv2,yv2,nanmean(sm_gleam,2),'LineStyle','none'); colorbar;%caxis([0 60]);  colormap(gca,jet(10));
    subplot(2,2,3);
    patch(xv,yv,nanmean(sm0,2),'LineStyle','none'); colorbar;
    subplot(2,2,4);
    plot(1:384,nanmean(sm,1),'b-','LineWidth',2); hold on;
    plot(1:384,nanmean(sm0,1),'r-','LineWidth',2);
    plot(13:384,nanmean(sm_gleam,1),'k-','LineWidth',2);
    
end
    
if evaluate_zwt
    for i = 1 : length(files1)
        filename = fullfile(files1(i).folder,files1(i).name);
        filename0= fullfile(files0(i).folder,files0(i).name);
        disp(['i = ' num2str(i) ': ' filename]);
        if i == 1
            zwt = ncread(filename,'ZWT');
            zwt = zwt(1:191);
            
            zwt0 = ncread(filename0,'ZWT');
            zwt0 = zwt0(1:191);
        else
            tmp = ncread(filename,'ZWT');
            zwt = cat(2,zwt,tmp(1:191));
            
            tmp = ncread(filename0,'ZWT');
            tmp = tmp(1:191);
            zwt0 = cat(2,zwt0,tmp(1:191));
            
        end
    end

    % Load Fan wtd data
    load('/Users/xudo627/projects/ELM_Runoff_Sensitivity/global/WTD_Fan_el_al_2013_05.mat');
    lon = ncread('/Users/xudo627/projects/ELM_Runoff_Sensitivity/global/inputdata/domain_lnd_uq_1_01_c210331.nc','xc');
    lat = ncread('/Users/xudo627/projects/ELM_Runoff_Sensitivity/global/inputdata/domain_lnd_uq_1_01_c210331.nc','yc');
    lon = lon(1:70302);
    lat = lat(1:70302);

    zwt_fan = NaN(length(xc),1);
    for i = 1 : length(xc)
        ind = find(lon == xc(i) & lat == yc(i));
        if isempty(ind)
            disp(['Cannot find cell for i = ' num2str(i)]);
        else
            zwt_fan(i) = -wtd05(ind);
        end
    end
    save('zwt_fan.mat','zwt_fan');

    figure; set(gcf,'Position',[10 10 1200 400]);
    subplot(2,2,1);
    patch(xv,yv,nanmean(zwt,2),'LineStyle','none'); colorbar; caxis([0 60]); colormap(gca,jet(10));
    subplot(2,2,2);
    patch(xv,yv,nanmean(zwt0,2),'LineStyle','none'); colorbar; caxis([0 60]); colormap(gca,jet(10));
    subplot(2,2,3);
    patch(xv,yv,zwt_fan,'LineStyle','none'); colorbar;caxis([0 60]);  colormap(gca,jet(10));
    subplot(2,2,4);
    patch(xv,yv,nanmean(zwt,2)-zwt_fan,'LineStyle','none'); colorbar;caxis([-30 30]);  colormap(gca,blue2red(11));
end