
% Different options when computing SOMs

% LAND + OCEAN
%fln='output_som4_wt_y_d_1_18-Apr-2018_16:17:50.mat';
%fln='output_som4_wt_y_d_2_18-Apr-2018_15:58:50.mat';
%fln='output_som6_wt_y_d_1_18-Apr-2018_15:29:20.mat';
%fln='output_som6_wt_y_d_2_18-Apr-2018_15:56:38.mat';
%fln='output_som8_wt_y_d_1_18-Apr-2018_17:03:04.mat';
%fln='output_som8_wt_y_d_2_18-Apr-2018_17:11:25.mat';
%fln='output_som9_wt_y_d_1_18-Apr-2018_16:54:38.mat';
%fln='output_som9_wt_y_d_2_18-Apr-2018_16:51:02.mat';
%fln='output_som12_wt_y_d_1_18-Apr-2018_17:22:36.mat';
%fln='output_som12_wt_y_d_2_18-Apr-2018_17:18:55.mat';
%fln='output_som12_wt_y_d_1_nodetrend_19-Apr-2018_14:25:02.mat';

% OCEAN ONLY
%fln='output_som4_d_1_01-May-2018_13:47:27.mat';
%fln='output_som6_d_1_02-May-2018_20:06:26.mat';
%fln='output_som8_d_1_02-May-2018_15:03:05.mat';
%fln='output_som9_d_1_01-May-2018_13:39:50.mat';
%fln='output_som12_d_1_02-May-2018_15:05:12.mat';
%fln='output_som15_d_1_02-May-2018_15:11:26.mat';

% For Hun
%fln='output_som10_25-Jan-2019_12:10:47.mat';
%fln='output_som12_25-Jan-2019_12:25:17.mat';
%fln='output_som8_25-Jan-2019_12:28:40.mat';
%fln='output_som6_yrs8501850_09-May-2019_13:38:12.mat';
%fln='output_som10_yrs8501850_09-May-2019_13:40:18.mat';
%fln='output_som12_yrs8501850_09-May-2019_13:53:03.mat';
%fln='output_som8_yrs8501850_09-May-2019_13:59:18.mat';

% SA & NA period of analysis
%fln='output_som6_yrs10001925_13-Sep-2019_13:34:26.mat';
%fln='output_som8_yrs10001925_13-Sep-2019_13:28:57.mat';
%fln='output_som12_yrs10001925_13-Sep-2019_13:33:49.mat';
%fln='output_som16_yrs10001925_13-Sep-2019_14:14:34.mat';
%fln='output_phyda_djf_som6_yrs10001925_18-Dec-2019_14:28:44.mat';
%fln='output_phyda_djf_som8_yrs10001925_19-Dec-2019_11:10:50.mat';
%fln='output_phyda_djf_som12_yrs10001925_19-Dec-2019_11:28:09.mat';
fln='output_phyda_djf_som16_yrs10001925_19-Dec-2019_11:29:22.mat';


fln0=fln;
%load(['./output_som/' fln])
load([fln])

%---------------------------
%      PLOT PATTERNS
%---------------------------

inputf='cesm_latlon.nc';

% define the output variable
outputf=['soms_sst_pdsi_' fln0(1:end-13) '.nc'];
%outputf=['soms_sst_pdsi_' fln0(1:end-13) '_matlabtest.nc'];
%outputf='test_output.nc';inputf='ccsm4_latlon.nc';
% grab lat/lon from model, put it into the output file
eval(['!ncks -A -v lat,lon ',inputf,' ',outputf]);

% PDSI
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:49:45.mat';
efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:40:23.mat';
pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';load([pth,efl],'Xa_m','xmeta','state_tp','reconYrs')
%analyrs0=850:1850;
%[~,~,ir0] = intersect(analyrs0,r_o:r_f);

k=find(state_tp=='q');
Xvar=Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
Xpdsi=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),reconYrs);
Xpdsi(xmeta{k}.lndidx,:)=Xvar;
Xpdsi=reshape(Xpdsi,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs);

% Detrend data for composites 
D0=reshape(Xpdsi(:,:,ir0),length(lat)*length(lon),length(ir0));
dt_pdsi=detrend(D0')';% Detrend the data, if desired

% Add time dimension
nccreate(outputf,'time','Dimensions',{'time',length(analyrs0)});
ncwrite(outputf,'time',analyrs0);
ncwriteatt(outputf,'time','units','Annual mean starting from April to the next calendar year March');
ncwriteatt(outputf,'time','long_name','Years CE');

nccreate(outputf,'bmus','Dimensions',{'time',length(analyrs0)});
ncwrite(outputf,'bmus',bmus);
ncwriteatt(outputf,'bmus','long_name','Best matching units');
ncwriteatt(outputf,'bmus','units','SOM node number');



for i=1:somRows*somCols

   % Get composites of PDSI over the SOMs
   cind=find(bmus==i);
   cmp_pdsi=permute(reshape(mean(dt_pdsi(:,cind),2),length(lat),length(lon)),[2 1]);
   cmp_pdsi(cmp_pdsi==NaN)=1e20;
   cnm=['cmp_pdsi_' num2str(i)];

   tmp_ptrn=permute(reshape(ptrns(i,:),length(lat),length(lon)),[2 1]);
   tmp_ptrn(tmp_ptrn==NaN)=1e20;
   pnm=['tmp_ptrn_' num2str(i)];
   
   nccreate(outputf,cnm,'Dimensions',{'lon',length(lon),'lat',length(lat)},'FillValue',1e20);
   ncwrite(outputf,cnm,cmp_pdsi);
   ncwriteatt(outputf,cnm,'missing_value',1e20);

   nccreate(outputf,pnm,'Dimensions',{'lon',length(lon),'lat',length(lat)},'FillValue',1e20);
   ncwrite(outputf,pnm,tmp_ptrn);
   ncwriteatt(outputf,pnm,'missing_value',1e20);

end

% Add note about what kind of SOM map it is
eval(['!ncatted -O --glb_att_add note="SOM rows = ' num2str(somRows) ', SOM columns = ' num2str(somCols) '" ' outputf]);
 

 

