%% standardized Extra-Tropical Index sETI Analysis Script

% --- Descripción ---

% Este script principal realiza un análisis del Índice Extratropical Estandarizado (SETI),

% --- Input ---

% - Datos brutos de temperatura (formato: [Ejemplo: formato de datos]).
% - Datos brutos de presión atmosférica (formato: [Ejemplo: formato de datos]).

% --- Output ---

% - Conjuntos de datos detallados del SETI, que incluyen:
% - SETI normalizado en formato tabular.
% - Gráficos de tendencias y patrones climáticos.
% - Informes estadísticos para la interpretación de resultados.

% --- Modificaciones ---

% [xxxx-xx-xx]: 
% - se cambio nanmean() por mean(, 'omitnan') 
% - se cambio las fechas de los datos MJO de HH 12 a HH 00
% - se quitaron los ticks de la barra de colores (length 0), dejando solo los labels
% - se cambio el orden de la linea segmentada negra del panel inferior para
% que no se reajuste automaticamente
% - agrego la siguiente linea: clear temp, en la seccion para sincronizar
% las temperaturas

% [2023-12-26]: 

% [2023-06-06]: 

% Autor:
% [Tu Nombre]
% [Tu Afiliación]
% [Fecha Inicial]

%%

clear all; clc

set(0,'defaultfigurecolor',[1 1 1])

addpath([pwd, '/../output'])

% time limits settings
auxset.requestday       = datenum(date);% current date
auxset.timeinterval     = [auxset.requestday-21 auxset.requestday+16];% x axis limits
auxset.timeticks        = datenum( auxset.timeinterval(1):2:auxset.timeinterval(2) )';% xticks 

%% gfs analysis 21 days before

% We verify that the GFS operational file is complete
file = 'forecast_16_days_after_oper.nc';
time = ncread(file, 'time');% forecast last 65 times (6 hourly timestep)
if length(time)<65
    error('The file has less than 65 records, try the next execution!, this error message was made on purpose...')
end

file = 'gfs_data.nc';% main file GFS analysis
time = ncread(file, 'time');% 6 hourly timestep

reftime = datenum(2019, 8, 13, 0, 0, 0);% reference date gfs_data.nc file
out2    =  datenum(reftime + hours(time));

% spatial coordinates
lon = ncread(file, 'lon');
lon = double(lon);
lat = ncread(file, 'lat');
lat = double(lat);

% west center
lat1      = find(lat == -52);
lon1      = find(lon == 14);
startLocW = [lon1-16 lat1-8 1 1];% the resolution is 0.5 degrees, if this is not correct change the positions to extract! 
countW    = [2*16+1 2*8+1 1 Inf]; 
strideW   = [1 1 1 1];% the file is [lon lat lev time] 

gh             = ncread(file, 'gh', startLocW, countW, strideW);% read geopotential height
gh             = double(squeeze(gh));
auxvar         = NaN(size(gh));
gh_west        = gh;
gh_backup_west = gh;

rm   = 40;% days to average
cpos = find(out2 == datenum(out2(end))-16);% current day position (00 UTC)

for t = 1 + cpos - 61*4:cpos% 40 days running mean 61 days before
    for i = 1:size(gh, 1)
        for j = 1:size(gh, 2)
            auxmean         = mean(gh(i, j, t-rm*4:t-1), 3, 'omitnan');
            auxstd          = std(gh(i, j, t-rm*4:t-1), [], 3, 'omitnan');
            auxvar(i, j, t) = (gh(i, j, t) - auxmean) ./ auxstd;
        end
    end
end
ghwest = mean(mean(auxvar, 1, 'omitnan'), 2, 'omitnan');
ghwest = squeeze(ghwest);

% east center
lat2      = find(lat == -58);
lon2      = find(lon == 66);
startLocE = [lon2-16 lat2-8 1 1]; 
countE    = [2*16+1 2*8+1 1 Inf]; 
strideE   = [1 1 1 1]; 

gh             = ncread(file, 'gh', startLocE, countE, strideE);
gh             = double(squeeze(gh));
auxvar            = NaN(size(gh));
gh_east        = gh;
gh_backup_east = gh;

for t = 1 + cpos - 61*4:cpos% 40 days running mean 61 days before
    for i = 1:size(gh, 1)
        for j = 1:size(gh, 2)
            auxmean         = mean(gh(i, j, t-rm*4:t-1), 3, 'omitnan');
            auxstd          = std(gh(i, j, t-rm*4:t-1), [], 3, 'omitnan');
            auxvar(i, j, t) = (gh(i, j, t) - auxmean) ./ auxstd;
        end
    end
end
gheast = mean(mean(auxvar, 1, 'omitnan'), 2, 'omitnan');
gheast = squeeze(gheast);

% calc of ETI
eti = gheast - ghwest;% east-west

seti = NaN(size(eti));
for t = cpos-21*4:cpos% 40 dias running mean 21 days before
    auxmean = mean(eti(t-rm*4:t-1), 1, 'omitnan');
    auxstd  = std(eti(t-rm*4:t-1), [], 1, 'omitnan');
    seti(t) = (eti(t) - auxmean) ./ auxstd;
end

anl_interval = cpos-21*4:cpos;% interval 21 days before

series.eti_anl       = eti(anl_interval);% backup serie!
series.seti_anl      = seti(anl_interval);% backup serie!
series.time_anl      = out2(anl_interval);% backup serie!
series.time_forecast = out2(end-64:end);% backup serie!

clear dataseti

auy               = series.seti_anl;
auxtime           = series.time_anl;
Time              = datetime(datevec(auxtime));
Time.Format       = 'uuuu-MM-dd HH:mm:ss';
dataseti{1}.value = timetable(Time, auy);

clear i j aux a file lat lon lat1 lat2 lon1 lon2 t time hrs anl_interval gheast ghwest eti seti out

%% gfs forecast 16 days after (gfs operational)

file = 'forecast_16_days_after_oper.nc';

if isfile(['../output/', file])
    disp(['GFS operationnal exists.'])
    
    % west center
    gh_w = ncread(file, 'gh', startLocW, countW, strideW);
    gh_w = double(squeeze(gh_w));
    if size(gh_w, 3)==65
        
        auxvar                       = NaN(size(gh_west));
        gh_west(:, :, end-64:end) = gh_w(:, :, end-64:end);
        
        for t=1+size(gh_west(:, :, 1:end-64), 3)-40*4:size(gh_west, 3)
            for i=1:size(gh_west, 1)
                for j=1:size(gh_west, 2)
                    auxmean         = mean(gh_west(i, j, t-rm*4:t-1), 3, 'omitnan');
                    auxstd          = std(gh_west(i, j, t-rm*4:t-1), [], 3, 'omitnan');
                    auxvar(i, j, t) = (gh_west(i, j, t) - auxmean) ./ auxstd;
                end
            end
        end
        ghwest = mean(mean(auxvar, 1, 'omitnan'), 2, 'omitnan');
        ghwest = squeeze(ghwest);
        
        % east center
        gh_e = ncread(file, 'gh', startLocE, countE, strideE);
        gh_e = double(squeeze(gh_e));
        
        auxvar = NaN(size(gh_east));
        gh_east(:, :, end-64:end) = gh_e(:, :, end-64:end);
        
        for t=1+size(gh_east(:, :, 1:end-64), 3)-40*4:size(gh_east, 3)
            for i=1:size(gh_east, 1)
                for j=1:size(gh_east, 2)
                    auxmean         = mean(gh_east(i,j,t-rm*4:t-1),3, 'omitnan');
                    auxstd          = std(gh_east(i,j,t-rm*4:t-1),[],3, 'omitnan');                    
                    auxvar(i, j, t) = (gh_east(i,j,t) - auxmean) ./ auxstd;
                end
            end
        end
        gheast = mean(mean(auxvar, 1, 'omitnan'), 2, 'omitnan');
        gheast = squeeze(gheast);
        
        eti = gheast - ghwest;
        
        seti = NaN(size(eti));
        for t=cpos:size(gh, 3)
            auxmean = mean(eti(t-rm*4:t-1),1, 'omitnan');
            auxstd  = std(eti(t-rm*4:t-1),[],1, 'omitnan');
            seti(t) = (eti(t) - auxmean) ./ auxstd;
        end
        
        series.eti_oper = eti(end-64:end);
        series.seti_oper = seti(end-64:end);
        
        auy     = series.seti_oper;
        auxtime = series.time_forecast;
    else
        disp(['GFS operationnal has less than 65 records... is not included.'])
        auy     = NaN;
        auxtime = series.time_anl(end);% current date
    end
else
    disp(['GFS operationnal does not exist.'])
    auy     = NaN;
    auxtime = series.time_anl(end);% current date
end

Time              = datetime(datevec(auxtime));
Time.Format       = 'uuuu-MM-dd HH:mm:ss';
dataseti{2}.value = timetable(Time, auy);

clear i j t seti eti file gheast ghwest

%% gfs forecast 16 days after (ensemble average)
   
file = ['forecast_16_days_after_eavg.nc'];

if isfile(['../output/', file])
    disp(['Ensemble average exists.'])
    
    % west center
    gh_w = ncread(file, 'gh', startLocW, countW, strideW);
    gh_w = double(squeeze(gh_w));
    if size(gh_w, 3) == 65
        
        auxvar = NaN(size(gh_west));
        gh_west(:, :, end-64:end) = gh_w(:, :, end-64:end);
        
        for t = 1+size(gh_west(:, :, 1:end-64),3)-40*4:size(gh_west, 3)
            for i=1:size(gh_west, 1)
                for j=1:size(gh_west, 2)
                    auxmean         = mean(gh_west(i, j, t-rm*4:t-1), 3, 'omitnan');
                    auxstd          = std(gh_west(i, j, t-rm*4:t-1), [], 3, 'omitnan');
                    auxvar(i, j, t) = (gh_west(i, j, t) - auxmean) ./ auxstd;
                end
            end
        end
        ghwest = mean(mean(auxvar, 1, 'omitnan'), 2, 'omitnan');
        ghwest = squeeze(ghwest);
        
        % east center
        gh_e = ncread(file, 'gh', startLocE, countE, strideE);
        gh_e = double(squeeze(gh_e));
        
        auxvar = NaN(size(gh_east));
        gh_east(:, :, end-64:end) = gh_e(:, :, end-64:end);
        
        for t = 1+size(gh_east(:, :, 1:end-64), 3)-40*4:size(gh_east, 3)
            for i=1:size(gh_east, 1)
                for j=1:size(gh_east, 2)
                    auxmean         = mean(gh_east(i, j, t-rm*4:t-1), 3, 'omitnan');
                    auxstd          = std(gh_east(i, j, t-rm*4:t-1), [], 3, 'omitnan');
                    auxvar(i, j, t) = (gh_east(i, j, t) - auxmean) ./ auxstd;
                end
            end
        end
        gheast = mean(mean(auxvar, 1, 'omitnan'), 2, 'omitnan');
        gheast = squeeze(gheast);
        
        eti = gheast - ghwest;
        
        seti=NaN(size(eti));
        for t=cpos:size(gh, 3)
            auxmean = mean(eti(t-rm*4:t-1), 1, 'omitnan');
            auxstd  = std(eti(t-rm*4:t-1), [], 1, 'omitnan');
            seti(t) = (eti(t) - auxmean) ./ auxstd;
        end
        
        series.eti_eavg = eti(end-64:end);
        series.seti_eavg = seti(end-64:end);
        
        auy     = series.seti_eavg;
        auxtime = series.time_forecast;
    else
        disp(['Ensemble average has less than 65 records... is not included.'])
        auy     = NaN;
        auxtime = series.time_anl(end);% current date
    end
else
    disp(['Ensemble average does not exist.'])
    auy     = NaN;
    auxtime = series.time_anl(end);% current date
end

Time              = datetime(datevec(auxtime));
Time.Format       = 'uuuu-MM-dd HH:mm:ss';
dataseti{3}.value = timetable(Time, auy);

clear i j t seti eti file gheast ghwest

%% loop ensemble gfs forecast 16 days after

for em=1:20    
    file = ['forecast_16_days_after_e', num2str(em, '%02d'), '.nc'];
    
    eval(['series.eti_e', num2str(em, '%02d'), ' = [];']);
    eval(['series.seti_e', num2str(em, '%02d'), ' = [];']);
    if isfile(['../output/', file])
        disp(['Ensemble member ', num2str(em, '%02d'), ' exists.'])
        
        % west center
        gh_w = ncread(file, 'gh', startLocW, countW, strideW);
        gh_w = double(squeeze(gh_w));
        
        if size(gh_w, 3) == 65
            auxvar = NaN(size(gh_west));
            gh_west(:, :, end-64:end) = gh_w(:, :, end-64:end);
            
            for t=1+size(gh_west(:, :, 1:end-64), 3)-40*4:size(gh_west, 3)
                for i=1:size(gh_west, 1)
                    for j=1:size(gh_west, 2)
                        auxmean         = mean(gh_west(i, j, t-rm*4:t-1), 3, 'omitnan');
                        auxstd          = std(gh_west(i, j, t-rm*4:t-1), [], 3, 'omitnan');
                        auxvar(i, j, t) = (gh_west(i, j, t) - auxmean) ./ auxstd;
                    end
                end
            end
            ghwest = mean(mean(auxvar, 1, 'omitnan'), 2, 'omitnan');
            ghwest = squeeze(ghwest);
            
            % east center
            gh_e = ncread(file, 'gh', startLocE, countE, strideE);
            gh_e = double(squeeze(gh_e));
            
            auxvar = NaN(size(gh_east));
            gh_east(:, :, end-64:end) = gh_e(:, :, end-64:end);
            
            for t=1+size(gh_east(:, :, 1:end-64), 3)-40*4:size(gh_east, 3)
                for i=1:size(gh_east, 1)
                    for j=1:size(gh_east, 2)
                        auxmean         = mean(gh_east(i, j, t-rm*4:t-1), 3, 'omitnan');
                        auxstd          = std(gh_east(i, j, t-rm*4:t-1), [], 3, 'omitnan');
                        auxvar(i, j, t) = (gh_east(i, j, t) - auxmean) ./ auxstd;
                    end
                end
            end
            gheast = mean(mean(auxvar, 1, 'omitnan'), 2, 'omitnan');
            gheast = squeeze(gheast);
            
            eti = gheast - ghwest;
            
            seti = NaN(size(eti));
            for t=cpos:size(gh, 3)
                auxmean = mean(eti(t-rm*4:t-1), 1, 'omitnan');
                auxstd  = std(eti(t-rm*4:t-1), [], 1, 'omitnan');
                seti(t) = (eti(t) - auxmean) ./ auxstd;
            end
                        
            eval(['series.eti_e', num2str(em, '%02d'), ' = eti(end-64:end);']);
            eval(['series.seti_e', num2str(em, '%02d'), ' = seti(end-64:end);']); 
            
            auy     = seti(end-64:end);
            auxtime = series.time_forecast;            
        else
            disp(['Ensemble member ', num2str(em, '%02d'),' has less than 65 records... is not included.'])            
            auy     = NaN;
            auxtime = series.time_anl(end);% current date
        end
    else
        disp(['Ensemble member ', num2str(em, '%02d'), ' does not exist.'])
        auy     = NaN;
        auxtime = series.time_anl(end);% current date
    end
    Time                   = datetime(datevec(auxtime));
    Time.Format            = 'uuuu-MM-dd HH:mm:ss';
    dataseti{3 + em}.value = timetable(Time, auy);
end

clear i j t seti eti file gheast ghwest em

%% synchronize seti series

seti = synchronize(dataseti{1}.value,...
    dataseti{2}.value,...
    dataseti{3}.value,...
    dataseti{4}.value,...
    dataseti{5}.value,...
    dataseti{6}.value,...
    dataseti{7}.value,...
    dataseti{8}.value,...
    dataseti{9}.value,...
    dataseti{10}.value,...
    dataseti{11}.value,...
    dataseti{12}.value,...
    dataseti{13}.value,...
    dataseti{14}.value,...
    dataseti{15}.value,...
    dataseti{16}.value,...
    dataseti{17}.value,...
    dataseti{18}.value,...
    dataseti{19}.value,...
    dataseti{20}.value,...
    dataseti{21}.value,...
    dataseti{22}.value,...
    dataseti{23}.value);

TT2 = seti;
TT3 = timetable2table(TT2); 
writetable(TT3,'../output/SETI_all.csv');
writetable(TT3,'../upload/SETI_all.csv');

%% Temperature forecast 

% info
name = {'santiago','curico','chillan','temuco','valdivia','futaleufu'};

rlat = [ -33.4450,...
    -34.9664,...
    -36.5872,...
    -38.7700,...
    -39.6506,...
    -43.1889 ];

rlon  = [ -70.6828,...
    -71.2167,...
    -72.0400,...
    -72.6369,...
    -73.0808,...
    -71.8525 ];

rlon = rlon + 360;

% OPERATIONAL
clear tmp
clear aux
clear minimox positionx minimoy positiony mas_cercanox mas_cercanoy
for i = 1:length(rlon)
    file    = ['temp_forecast_16_days_after_oper.nc'];
    if isfile(['../output/', file])
        auxt     = ncread(file, 'time');
        auxtime = datenum( datenum(date) + hours(auxt) );
        
        lon = ncread(file, 'lon');
        lat = ncread(file, 'lat');
        
        x = lon;
        y = lat;
        
        [minimox(i), positionx(i)] = min(abs(x - rlon(i)));
        [minimoy(i), positiony(i)] = min(abs(y - rlat(i)));
        mas_cercanox(i) = x(positionx(i));
        mas_cercanoy(i) = y(positiony(i));
        
        aux{i} = ncread(file, '2t', [positionx(i) positiony(i) 1 1], [1 1 1 inf]);
        aux{i} = squeeze(aux{i}) - 273.15;
    else
        auxtime = auxset.timeinterval(end);
        aux{i}  = NaN;
    end
end

% curico serie
% curico{1}.value = aux{2};
% curico{1}.time  = auxtime;
auy             = aux{2};
Time            = datetime(datevec(auxtime));
Time.Format     = 'uuuu-MM-dd HH:mm:ss';
curi{1} = timetable(Time, auy);

% central chile mean
auy             = mean([aux{1}, aux{2}, aux{3}], 2, 'omitnan');
Time            = datetime(datevec(auxtime));
Time.Format     = 'uuuu-MM-dd HH:mm:ss';
tmp{1}.value{1} = timetable(Time, auy);

% south central chile mean
auy             = mean([aux{4}, aux{5}, aux{6}], 2, 'omitnan');
Time            = datetime(datevec(auxtime));
Time.Format     = 'uuuu-MM-dd HH:mm:ss';
tmp{1}.value{2} = timetable(Time, auy);

% EAVG
%clear tmp
clear aux
clear minimox positionx minimoy positiony mas_cercanox mas_cercanoy
for i = 1:length(rlon)
    file    = ['temp_forecast_16_days_after_eavg.nc'];
    if isfile(['../output/', file])
        auxt     = ncread(file, 'time');
        auxtime = datenum( datenum(date) + hours(auxt) );
        
        lon = ncread(file, 'lon');
        lat = ncread(file, 'lat');
        
        x = lon;
        y = lat;
        
        [minimox(i), positionx(i)] = min(abs(x - rlon(i)));
        [minimoy(i), positiony(i)] = min(abs(y - rlat(i)));
        mas_cercanox(i) = x(positionx(i));
        mas_cercanoy(i) = y(positiony(i));
        
        aux{i} = ncread(file, '2t', [positionx(i) positiony(i) 1 1], [1 1 1 inf]);
        aux{i} = squeeze(aux{i}) - 273.15;
    else
        auxtime = auxset.timeinterval(end);
        aux{i}  = NaN;
    end
end

% curico serie
% curico{1}.value = aux{2};
% curico{1}.time  = auxtime;
auy             = aux{2};
Time            = datetime(datevec(auxtime));
Time.Format     = 'uuuu-MM-dd HH:mm:ss';
curi{2} = timetable(Time, auy);

% central chile mean
auy             = mean([aux{1}, aux{2}, aux{3}], 2, 'omitnan');
Time            = datetime(datevec(auxtime));
Time.Format     = 'uuuu-MM-dd HH:mm:ss';
tmp{2}.value{1} = timetable(Time, auy);

% south central chile mean
auy             = mean([aux{4}, aux{5}, aux{6}], 2, 'omitnan');
Time            = datetime(datevec(auxtime));
Time.Format     = 'uuuu-MM-dd HH:mm:ss';
tmp{2}.value{2} = timetable(Time, auy);

% ENSEMBLE
for em=1:20
    clear aux
    clear minimox positionx minimoy positiony mas_cercanox mas_cercanoy
    for i = 1:length(rlon)
        file    = ['temp_forecast_16_days_after_e', num2str(em, '%02d'), '.nc'];
        if isfile(['../output/', file])
            auxt     = ncread(file, 'time');
            auxtime = datenum( datenum(date) + hours(auxt) );
            
            lon = ncread(file, 'lon');
            lat = ncread(file, 'lat');
            
            x = lon;
            y = lat;
            [minimox(i), positionx(i)] = min(abs(x - rlon(i)));
            [minimoy(i), positiony(i)] = min(abs(y - rlat(i)));
            mas_cercanox(i) = x(positionx(i));
            mas_cercanoy(i) = y(positiony(i));
            
            aux{i} = ncread(file, '2t', [positionx(i) positiony(i) 1 1], [1 1 1 inf]);
            aux{i} = squeeze(aux{i}) - 273.15;
        else
            auxtime = auxset.timeinterval(end);
            aux{i}  = NaN;
        end
    end
    
    % curico serie
%     curico{em + 1}.value = aux{2};
%     curico{em + 1}.time  = auxtime;
    
    auy             = aux{2};
    Time            = datetime(datevec(auxtime));
    Time.Format     = 'uuuu-MM-dd HH:mm:ss';
    curi{em + 2} = timetable(Time, auy);
    
    % central chile mean
    auy                  = mean([aux{1}, aux{2}, aux{3}], 2, 'omitnan');
    Time                 = datetime(datevec(auxtime));
    Time.Format          = 'uuuu-MM-dd HH:mm:ss';
    tmp{em + 2}.value{1} = timetable(Time, auy);
    
    % south central chile mean
    auy                  = mean([aux{4}, aux{5}, aux{6}], 2, 'omitnan');
    Time                 = datetime(datevec(auxtime));
    Time.Format          = 'uuuu-MM-dd HH:mm:ss';
    tmp{em + 2}.value{2} = timetable(Time, auy);
end

%% synchronize temperature series

clear temp

temp.up = synchronize(tmp{1}.value{1},...
    tmp{2}.value{1},...
    tmp{3}.value{1},...
    tmp{4}.value{1},...
    tmp{5}.value{1},...
    tmp{6}.value{1},...
    tmp{7}.value{1},...
    tmp{8}.value{1},...
    tmp{9}.value{1},...
    tmp{10}.value{1},...
    tmp{11}.value{1},...
    tmp{12}.value{1},...
    tmp{13}.value{1},...
    tmp{14}.value{1},...
    tmp{15}.value{1},...
    tmp{16}.value{1},...
    tmp{17}.value{1},...
    tmp{18}.value{1},...
    tmp{19}.value{1},...
    tmp{20}.value{1},...
    tmp{21}.value{1},...
    tmp{22}.value{1});

temp.bottom = synchronize(tmp{1}.value{2},...
    tmp{2}.value{2},...
    tmp{3}.value{2},...
    tmp{4}.value{2},...
    tmp{5}.value{2},...
    tmp{6}.value{2},...
    tmp{7}.value{2},...
    tmp{8}.value{2},...
    tmp{9}.value{2},...
    tmp{10}.value{2},...
    tmp{11}.value{2},...
    tmp{12}.value{2},...
    tmp{13}.value{2},...
    tmp{14}.value{2},...
    tmp{15}.value{2},...
    tmp{16}.value{2},...
    tmp{17}.value{2},...
    tmp{18}.value{2},...
    tmp{19}.value{2},...
    tmp{20}.value{2},...
    tmp{21}.value{2},...
    tmp{22}.value{2});

TT2 = temp.up;
TT3 = timetable2table(TT2);
writetable(TT3,'../output/tmp_up.csv');
writetable(TT3,'../upload/tmp_up.csv');

TT2 = temp.bottom;
TT3 = timetable2table(TT2);
writetable(TT3,'../output/tmp_bottom.csv');
writetable(TT3,'../upload/tmp_bottom.csv');

%% figure settings (top panel)

h = figure('Position', [100, 100, 960, 620]);

ax(1) = subplot(2, 1, 1);

aux = timetable2table(seti);

% seti analysis
x   = datenum( table2array( aux(:, 1) ) );
y   = table2array( aux(:, 2) );
plot(ax(1), x, y, 'linewidth', 1.5, 'Color', [0 0.3 0.9])
hold on

% seti ensemble
for i = 5:size(aux, 2)
    x   = datenum( table2array( aux(:, 1) ) );
    y   = table2array( aux(:, i) );
    plot(ax(1), x, y, 'linewidth', 1.5, 'Color', [0.9 0.6 0.5])
end

% seti operational
x   = datenum( table2array( aux(:, 1) ) );
y   = table2array( aux(:, 3) );
pl2 = plot(ax(1), x, y, 'linewidth', 1.5, 'Color', [0 0.5 0.2]);


% seti average
x   = datenum( table2array( aux(:, 1) ) );
y   = table2array( aux(:, 4) );
pl1 = plot(ax(1), x, y, 'linewidth', 1.5, 'Color', [1 0 0]);

% settings

% title
%t = title(['From GFS analysis and forecast (21+16 days) ',datestr(out2(end-64),'dd-mmm-yyyy HH:MM:SS')],'FontWeight','b','FontSize',13);
%t1 = get(t, 'position');
%set(t, 'position', [t1(1) t1(2)+2.3 t1(3)])

% yaxis
%ylabel({'Standardized';'Extra Tropical Index [SETI]'},'FontWeight','b','FontSize',12)
ylabel({'SETI [\sigma]'},'FontWeight','b','FontSize',12)
ylim([-5 5])

NumTicks = 11;
L = get(gca,'YLim');
set(gca,'YTick',linspace(L(1),L(2),NumTicks))

% xaxis
% xlabel([],'FontWeight','b','FontSize',12)

xlim(auxset.timeinterval)
set(ax(1), 'XTick', auxset.timeticks)
datetick(ax(1),'x','dd','keeplimits','keepticks')

% general 
set(ax(1), 'TickDir', 'out');
grid(ax(1), 'on')

% plot
plot([1 1]*out2(end-64), ylim, '--', 'Color', [0 0 0.3], 'linewidth', 1.5)

% horizontal patch settings
limx = xlim;
f1 = patch([limx(1) limx(2) limx(2) limx(1)],[-1 -1 1 1],[0.5608    0.6980    0.5373],'LineStyle','none','FaceAlpha',0.4);
uistack(f1,'bottom')

% legend
if (~isempty(pl1) & ~isempty(pl2))
    legend([pl1 pl2],'GEFS avg','GFS','location','southeast')
elseif (~isempty(pl1) & isempty(pl2))
    legend([pl1],'GEFS avg','location','southeast')
elseif (isempty(pl1) & ~isempty(pl2))
    legend([pl2],'GFS','location','southeast')
else
end

dim = [.14 .61 .3 .3];

texto = {['From GFS analysis and forecast (21+16 days) '];...
    ['                    ',datestr(out2(end-64),'dd-mmm-yyyy HH:MM:SS')]};
ann = annotation('textbox',dim,'String',texto,...
    'backgroundcolor','w','FitBoxToText','on');
ann.FontWeight = 'bold';
ann.FontSize = 10;

%% figure settings (bottom panel)

% read MJO data
fname = [pwd,'/../data/series/amp.nc'];% amplitude
time1 = ncread(fname, 'T');
time1 = datenum(datetime(time1, 'convertfrom', 'juliandate'));
time1 = datevec(time1);
time1 = time1(:, 1:3);
time1 = datenum(time1);
amp1 = ncread(fname, 'amplitude');

fname = [pwd,'/../data/series/pha.nc'];% phase
time2 = ncread(fname, 'T');
time2 = datenum(datetime(time2, 'convertfrom', 'juliandate'));
time2 = datevec(time2);
time2 = time2(:, 1:3);
time2 = datenum(time2);
pha2 = ncread(fname, 'phase');

auxtime = [datenum(date)-100: datenum(date)]';
% auxtime = datenum(auxtime + duration(12, 0, 0));

amp = NaN(size(auxtime));
pha = NaN(size(auxtime));

for i = 1:length(auxtime)
    idx1 = auxtime(i) == time1;
    idx2 = auxtime(i) == time2;
    
    if ~isempty(time1(idx1))
        amp(i) = amp1(idx1);
    end
    if ~isempty(time2(idx2))
        pha(i) = pha2(idx2);
    end
end

auxtime = auxtime(~isnan(pha));
amp     = amp(~isnan(pha));
pha     = pha(~isnan(pha));

auxtime = auxtime(~isnan(amp));
pha     = pha(~isnan(amp));
amp     = amp(~isnan(amp));

nphase = 8;

auxcmap = [0.1961    0.5333    0.7412;
           0.4000    0.7608    0.6471;
           0.6706    0.8667    0.6431;
           0.9020    0.9608    0.5961;
           0.9961    0.8784    0.5451;
           0.9569    0.4275    0.2627;
           0.7098    0.2314    0.1294;
           0.8118    0.5020    0.8745];

ax(2) = subplot(2, 1, 2);

% yaxis right 
yyaxis right
ax(2).YAxis(2).Color = 'k';
f1 = patch([limx(1) limx(2) limx(2) limx(1)], [0 0 1 1], [0.5608    0.6980    0.5373], 'LineStyle', 'none', 'FaceAlpha', 0.4);
hold on
if exist('auxtime','var')
    if (~isempty(auxtime))
        for t=1:size(auxtime,1)
            if pha(t)>8
            else
                scatter(auxtime(t), amp(t), 'ok', 'SizeData', 80, 'MarkerFaceColor', [auxcmap(pha(t),:)], 'MarkerEdgeColor', 'k')
                hold on
            end
        end
    end
end
ylabel(['MJO Amplitude'], 'FontWeight', 'b', 'FontSize', 12)
ylim([0 6])
NumTicks = 7;
L = get(gca,'YLim');
set(gca,'YTick',linspace(L(1),L(2),NumTicks))

% yaxis left
yyaxis left
ax(2).YAxis(1).Color = 'k';

clear p

% codigo estacion general freire
cod = 340031;

% read Tx
data = readtable([pwd,'/../data/series/T_max'], 'HeaderLines', 0);
data = table2array(data);

cods = NaN(length(data(1,:)),1);
for i = 2:length(data(1,:))
    cods(i) = str2double(data(1,i));    
end

idx = find(cods == cod);

if (~isempty(idx))
    clear time
    clear value
    cont = 1;
    for i = 16:length(data(:,1))
        time(cont,1) = datenum(data{i,1});
        value(cont,1) = str2double(data(i,idx));
        cont = cont + 1;
    end    
    value(value==-9999)=NaN;
    p(1) = plot(time,value,'-','linewidth',1.5,'Color',[1 0.3 0.3],'DisplayName','T max');
    hold on
end

% read Tm
data = readtable([pwd,'/../data/series/T_mean'], 'HeaderLines', 0);
data = table2array(data);

cods = NaN(length(data(1,:)),1);
for i = 2:length(data(1,:))
    cods(i) = str2double(data(1,i));    
end

idx = find(cods == cod);

if (~isempty(idx))
    clear time
    clear value
    cont = 1;
    for i = 16:length(data(:,1))
        time(cont,1) = datenum((data{i,1}));
        value(cont,1) = str2double(data(i,idx));
        cont = cont + 1;
    end    
    value(value==-9999)=NaN;
    p(2) = plot(time,value,'-','linewidth',1.5,'Color',[0.2 0.5 0.5],'DisplayName','T mean');
    hold on
end

% read Tn
data = readtable([pwd,'/../data/series/T_min'], 'HeaderLines', 0);
data = table2array(data);

cods = NaN(length(data(1,:)),1);
for i = 2:length(data(1,:))
    cods(i) = str2double(data(1,i));    
end

idx = find(cods == cod);

if (~isempty(idx))
    clear time
    clear value
    cont = 1;
    for i = 16:length(data(:,1))
        time(cont,1) = datenum(data{i,1});
        value(cont,1) = str2double(data(i,idx));
        cont = cont + 1;
    end    
    value(value==-9999)=NaN;
    p(3) = plot(time,value,'-','linewidth',1.5,'Color',[0.1 0 0.8],'DisplayName','T min');
    hold on
end


% plot forecast curico temperature
for i = 3:length(curi)
    plot(datenum(curi{i}.Time), curi{i}.auy, '-', 'linewidth', 1.5, 'Color', [0.9 0.6 0.5]);
end
i = 1;
plot(datenum(curi{i}.Time), curi{i}.auy, '-', 'linewidth', 1.5, 'Color', [0 0.5 0.2]);
i = 2;
plot(datenum(curi{i}.Time), curi{i}.auy, '-', 'linewidth', 1.5, 'Color', [1 0 0]);

curim = synchronize(curi{1},...
    curi{2},...
    curi{3},...
    curi{4},...
    curi{5},...
    curi{6},...
    curi{7},...
    curi{8},...
    curi{9},...
    curi{10},...
    curi{11},...
    curi{12},...
    curi{13},...
    curi{14},...
    curi{15},...
    curi{16},...
    curi{17},...
    curi{18},...
    curi{19},...
    curi{20},...
    curi{21},...
    curi{22});

TT2 = curim;
TT3 = timetable2table(TT2);

% 
aux = TT3;
x   = datenum( table2array( aux(:, 1) ) );
y   = table2array( aux(:, 3:end) );

% boxplot(y', 'positions', x)

% pl2 = plot(ax(1), x, y, 'linewidth', 1.5, 'Color', [0 0.5 0.2]);



altura = data{5,idx};
latitud = abs(str2double(data{6,idx}));
longitud = abs(str2double(data{7,idx}));
name_estacion = data{4,idx};


% settings

% title
% t=title(['Curico DGA 34.9811°S, 71.2361°W, 195m'], 'FontWeight', 'b', 'FontSize', 10.3);
t = title([name_estacion,' ',num2str(latitud), '°S, ', num2str(longitud), '°W, ', altura, 'm'], 'FontWeight', 'b', 'FontSize', 10.3);
set(t, 'horizontalAlignment', 'right')
set(t, 'units', 'normalized')
t1 = get(t, 'position');
set(t, 'position', [1 t1(2) t1(3)])

% yaxis
ylabel(['Temperature [°C]'],'FontWeight','b','FontSize',12)
% ylim([-10 40])

% xaxis
xlim(auxset.timeinterval)
set(ax(2), 'XTick', auxset.timeticks)
datetick('x', 'dd', 'keeplimits', 'keepticks')

% general
set(ax(2), 'TickDir', 'out');
grid(ax(2), 'on')

% colorbar
cb = colorbar('Southoutside', 'FontWeight', 'bold', 'Linewidth', 1, 'XTickLabel', 0:8);
ylabel(cb, 'MJO Phase', 'FontWeight', 'bold', 'FontSize', 12);

% Establecer la longitud de los ticks en 0
cb.TickLength = 0;

colormap(auxcmap)
caxis([0 8])

% plot
YLIM = ylim;
plot([1 1]*out2(end-64), YLIM, '--', 'Color', [0 0 0.3], 'linewidth', 1.5)

% legend
if (~isempty(p))
    legend([p(1) p(2) p(3)],'T max','T mean','T min', 'location', 'northwest')
end

% 'DisplayName',
% 'HandleVisibility', 'off'

% save
saveas(h,[pwd,'/../figures/SETI forecast From GFS ',datestr(out2(end-64),'yyyy-mm-dd HH:MM:SS'),'.png'])
close(h)

return
