
% 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';

% SA & NA time period
%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';

% For Hun
%fln='output_som12_25-Jul-2018_05:12:34.mat';

%fln='output_som12_23-Jan-2019_20:14:26.mat';
%fln='output_som10_24-Jan-2019_16:15:05.mat';
%fln='output_som8_23-Jan-2019_03:21:03.mat';


% SA and NA analysis
%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';
%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';

fln0=fln;
%load(['./som-output/' fln])
load([fln])

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

load('./colormaps/distinct_colors_20.mat')
%load('./colormaps/WhiteYellowOrangeRed.mat')
C=load('./colormaps/NCV_blue_red.mat');cmap_br=C.cmap;
C=load('./colormaps/MPL_BrBG.mat');cmap_bg=C.cmap;


make_figs='y';
plot_sp='n';% make spatial plots?


%analyrs4=800:1925;
analyrs4=1000:1925;
[~,~,ir4] = intersect(analyrs4,r_o:r_f);

% % PULL OUT SPATIAL MEAN OVER AMERICAN SOUTHWEST
% %load('/d1/nsteiger/climate-data/cesm-lme/Nrth_A_Swthwst_mask.mat')
% varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/010'];
% sftlf=permute(ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'landmask'),[2 1]);
% %lat=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lat');
% %lon=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lon');
% % NASW = box of 235:255 E (125:105 W); 31:42 N (just the US, not Mexico)
% % Find nearest lat/lon equal to or outside bounds
% lonext=[235 255];latext=[31 42];
% ln_a=find((abs(lon-lonext(1))-min(abs(lon-lonext(1))))<1e-5);
% ln_b=find((abs(lon-lonext(2))-min(abs(lon-lonext(2))))<1e-5);
% lt_a=find((abs(lat-latext(1))-min(abs(lat-latext(1))))<1e-5);
% lt_b=find((abs(lat-latext(2))-min(abs(lat-latext(2))))<1e-5);
% nasw=zeros(size(sftlf));
% nasw(lt_a:lt_b,ln_a:ln_b)=sftlf(lt_a:lt_b,ln_a:ln_b);
% 
% nasw(nasw==0)=NaN;
% Xv_msk=bsxfun(@times,Xpdsi,nasw);
% A=cosd(repmat(xmeta{k}.lat,[1 length(xmeta{k}.lon)]));
% Xa_di=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX
% 
% % How to compute the megadroughts? 1 => 11 year running mean, 2 => 2S2E
% %ds4=mdrought(Xa_di(ir4),1,'y'); % 'y' indicates that the mean is removed
% ds4=mdrought(Xa_di(ir4),1,'y'); % 'y' indicates that the mean is removed
% % Most severe droughts have severity < -10 and all occur before 1600 (last one ends in 1590)
% si=find(ds4(:,3)<-10);
% %si=find(ds4(:,2)+analyrs1(1)<analyrs1(end));
% ds0=ds4(si,:);
% 

%-------------------------------
% LOAD DROUGHT VARIABLE
%load('svmd_ens_90pagree.mat')
%load('svmd_ens_50pagree.mat')
%ds0=svmd;



% if plot_sp=='y'
% 
% % Detrend data for composites 
% D0=reshape(Xpdsi(:,:,ir0),length(lat)*length(lon),length(ir0));
% dt_pdsi=detrend(D0')';% Detrend the data, if desired
% 
% for i=1:somRows*somCols
% 
%    % Get composites of PDSI over the SOMs
%    cind=find(bmus==i);
%    varX=reshape(mean(dt_pdsi(:,cind),2),length(lat),length(lon));
%    varZ=cat(2,varX,varX(:,1));
% 
%    s=load('coast');mlat=s.lat;mlon=s.long;
% 
%    figure
%    hold on
%    %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
%    h1=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
%    %h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 max(lat)]);
%    title(['Pattern ' num2str(i)],'fontsize',18,'Interpreter','none')
%    pcolorm(lat,lon,reshape(ptrns(i,:),length(lat),length(lon)));
%    h2=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
%    pcolorm(lat,[lon; 360],varZ);
%    %pcolorm(lat,[lon; 360],varZ); % add additional value to avoid line in pcolor plots
%    plotm(mlat,mlon,'color',[0.5 0.5 0.5])
%    colormap(h1,cmap_br) 
%    colormap(h2,cmap_bg) 
%    set(gca,'fontsize',18)
%    %set(get(colorbar('location','southoutside'),'xlabel'),'string',sn{j},'fontsize',14)
%    %set(get(colorbar('location','southoutside'),'xlabel'),'string','std units','fontsize',18)
%    caxis(h1,[-1.5 1.5])
%    caxis(h2,[-1.25 1.25])
%    %cbh = findobj( 0, 'tag', 'Colorbar' ); % and
%    %set(cbh,'YTick',[-1.5:0.25:1.5])
%    %delete( cbh )
%    tightmap
%    hold off
% 
% 
% end
% 
% end

%=====================================
% LOAD ENSEMBLE MEGADROUGHT ESTIMATES

% Load both south american and north america drought indices
load('chile_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
xa_di_sa=Xa_di;
xa_di_ens_sa=Xa_di_ens;
ds_sa=svmd;
mdagr_sa=ens_drght;
sa_yrs=analyrs;

load('nasw_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
xa_di_na=Xa_di;
xa_di_ens_na=Xa_di_ens;
ds_na=svmd;
mdagr_na=ens_drght;
na_yrs=analyrs;

% Droughts in the real data mean
ds0_n=mdrought(xa_di_na(analyrs),1,'y');
ds0_s=mdrought(xa_di_sa(analyrs),1,'y');

% Droughts too frequent, throw out based on relative threshold
%dcut=67;
dcut=50;
np67=prctile(ds0_n(:,3),dcut);
rmv=find(ds0_n(:,3)>np67);
ds0_n(rmv,:)=[];

sp67=prctile(ds0_s(:,3),dcut);
rmv=find(ds0_s(:,3)>sp67);
ds0_s(rmv,:)=[];


% MAKE DROUGHT VECTORS
% North America
dn0=ds0_n(:,1:2)';dn0=dn0(:);

% Fill out all years of droughts in vector
k=1;dn1=[];
while k<length(dn0)
   dn1=cat(2,dn1,dn0(k):dn0(k+1));
   k=k+2;
end

% South America
ds0=ds0_s(:,1:2)';ds0=ds0(:);

% Fill out all years of droughts in vector
k=1;ds1=[];
while k<length(ds0)
   ds1=cat(2,ds1,ds0(k):ds0(k+1));
   k=k+2;
end


% Figure out indices of neither, NA, SA, both
inthr=setdiff(1:length(analyrs),union(dn1,ds1)); % neither
[~,ina,isa]=setxor(dn1,ds1); % indices of NA only and SA only
ibth=intersect(dn1,ds1); % both

%------------------------------------------------------------------------------------------
% Get BMUs for Mdroughts during 800-1600, for non-drought yrs 800-1600, all years 1601-1925
%------------------------------------------------------------------------------------------

% reference frequency
if length(bmus)~=length(analyrs);error('bmu years not consistent with droughts');end
rff=histcounts(bmus,0.5:1:(somRows*somCols+0.5))/length(bmus);

% Frequency during either, 'both', neither
dfc_b=histcounts(bmus(ibth),0.5:1:(somRows*somCols+0.5))/length(ibth);
dfc_n=histcounts(bmus(ina),0.5:1:(somRows*somCols+0.5))/length(ina);
dfc_s=histcounts(bmus(isa),0.5:1:(somRows*somCols+0.5))/length(isa);
dfc_t=histcounts(bmus(inthr),0.5:1:(somRows*somCols+0.5))/length(inthr);

% Frequency change
disp('Frequency change: both')
(dfc_b-rff)./rff
disp('Frequency change: nh')
(dfc_n-rff)./rff
disp('Frequency change: sh')
(dfc_s-rff)./rff
disp('Frequency change: neither')
(dfc_t-rff)./rff

return

% rsz=500;
% a2fc_w=zeros(rsz,somRows*somCols);
% a2fc_d=zeros(rsz,somRows*somCols);
% for i=1:rsz
%    % Sample the reference period to get ranges for droughts and pluvials
%    wrs=randsample(fltm,length(wf),'true');
%    drs=randsample(fltm,length(df),'true');
%    a2fc_w(i,:)=histcounts(wrs,0.5:1:(somRows*somCols+0.5))/length(wrs);
%    a2fc_d(i,:)=histcounts(drs,0.5:1:(somRows*somCols+0.5))/length(drs);
% end


% BMU frequency during megadroughts (note droughts are determined analysis period relative)
% UNCOMMENT BELOW FOR SINGLE LOCATION MEGADROUGHTS
% df=[];
% for i=1:length(ds0)
%   df0=ds0(i,1):ds0(i,2);
%   df=cat(2,df,df0);
% end
% dfc=histcounts(bmus(df),0.5:1:(somRows*somCols+0.5))/length(df);
% 
% % Frequency during non-droughts
% [~,ay1,~] = intersect(analyrs0,analyrs1);
% wf=setdiff(ay1(:),df(:));
% wfc=histcounts(bmus(wf),0.5:1:(somRows*somCols+0.5))/length(wf);
% 
% % Frequency during 800-1600
% efc=histcounts(bmus(ay1),0.5:1:(somRows*somCols+0.5))/length(ay1);
% 
% % Frequency during 1600-1925 period
% [~,ay2,~] = intersect(analyrs0,analyrs2);
% fltm=bmus(ay2);
% rff=histcounts(fltm,0.5:1:(somRows*somCols+0.5))/length(fltm);
% 
% rsz=500;
% a2fc_w=zeros(rsz,somRows*somCols);
% a2fc_d=zeros(rsz,somRows*somCols);
% for i=1:rsz
%    % Sample the reference period to get ranges for droughts and pluvials
%    wrs=randsample(fltm,length(wf),'true');
%    drs=randsample(fltm,length(df),'true');
%    a2fc_w(i,:)=histcounts(wrs,0.5:1:(somRows*somCols+0.5))/length(wrs);
%    a2fc_d(i,:)=histcounts(drs,0.5:1:(somRows*somCols+0.5))/length(drs);
% end

% Get percentiles of the above uncertainty estimates
% THESE ARE WRONG...
%a2fc_wp=prctile(a2fc_w,[5,50,95],1);
%a2fc_dp=prctile(a2fc_d,[5,50,95],1);

 
% Frequency for each drought
%ds_pfrq=zeros(length(ds0),somRows*somCols);
%for i=1:length(ds0)
%  ds_pfrq(i,:)=histcounts(bmus(ds0(i,1):ds0(i,2)),0.5:1:(somRows*somCols+0.5))/length(ds0(i,1):ds0(i,2));
%end

disp(fln0)
disp('Drought percent change of the frequency of the pattern')
pd=(dfc-rff)./rff;
disp(pd)



%----------------------------
% PLOT BMUS FOR EACH PATTERN
%----------------------------

%C=load('./colormaps/MPL_BrBG.mat');
%brn=C.cmap(1,:);grn=C.cmap(end,:);

figure
hold on
% Mean freq during mega D
%errorbar(1:somRows*somCols,rff,2*std(a2fc_w),'o','markeredgecolor',[0.5 0.5 0.5],...
%'color',grn,'linewidth',2)
errorbar(1:somRows*somCols,mean(a2fc_d),2*std(a2fc_d),'o','markeredgecolor',[0.5 0.5 0.5],...
'color',[0.5 0.5 0.5],'linewidth',2)
% Mean freq 1601-1925
%plot(1:somRows*somCols,efc,'p','markersize',10,'color',clrs.blue,'markerfacecolor',clrs.blue)
plot(1:somRows*somCols,dfc,'s','markersize',7,'color',clrs.red,'markerfacecolor',clrs.red)
plot(1:somRows*somCols,wfc,'^','markersize',7,'color',clrs.black,'markerfacecolor',clrs.black)
%plot(1:somRows*somCols,efc,'p','markersize',10,'color',clrs.blue,'markerfacecolor',clrs.blue)
%legend('1601-1925','800-1600','Megadrought 800-1600','Non-drought 800-1600','location','northoutside')
legend('1601-1925','Megadrought 800-1600','Non-drought 800-1600','location','northoutside')
xlabel('Node number','fontsize',18)
ylabel('Frequency','fontsize',18)
xticks(1:somRows*somCols)
xlim([0.5 somRows*somCols+0.5])
%ax = gca;ax.YAxisLocation='right';
set(gca,'fontsize',16)
fig = gcf;fig.PaperUnits = 'inches';
%fig.PaperPosition = [0 2 3.5 5]; % [left bottom width height]
fig.PaperPosition = [0 2 4.5 5]; % [left bottom width height]
box on

if  make_figs=='y'
   print('-dpng',['/home/nsteiger/Dropbox/megadrought-nasw-2019/som' num2str(somRows*somCols) '_ptrnfrq_' timestp2 '.png'],'-r300')
end















% ORIGINAL VERSION OF THE CODE
% 
% if plot_p=='y'
%    % Pull out proxy locations
%    %plat=nan(length(proxy),1);
%    %plon=nan(length(proxy),1);
%    j=1;k=1;l=1;
%    for i=1:length(proxy)
%       if ismember(proxy{i}.ptype,{'tree','Tree Rings'})
% 	 tlat(j)=proxy{i}.location(1);
% 	 tlon(j)=proxy{i}.location(2);
% 	 j=j+1;
%       elseif ismember(proxy{i}.ptype,{'coral','sclerosponge','Corals and Sclerosponges','Bivalve'})
% 	 clat(k)=proxy{i}.location(1);
% 	 clon(k)=proxy{i}.location(2);
% 	 k=k+1;
%       else % Other proxy type 
% 	 olat(l)=proxy{i}.location(1);
% 	 olon(l)=proxy{i}.location(2);
% 	 l=l+1;
%       end 
%    end
% end
% 
% if plot_sp=='y'
% 
% for i=1:somRows*somCols
% 
% 
%    s=load('coast');mlat=s.lat;mlon=s.long;
% 
%    figure
%    hold on
%    %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
%    h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
%    %h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 max(lat)]);
%    title(['Pattern ' num2str(i)],'fontsize',18,'Interpreter','none')
%    pcolorm(lat,lon,reshape(ptrns(i,:),length(lat),length(lon)));
%    %pcolorm(lat,[lon; 360],varZ); % add additional value to avoid line in pcolor plots
%    plotm(mlat,mlon,'color',[0.5 0.5 0.5])
%    if plot_p=='y';
%       h1=plotm(tlat,tlon,'s','color',clrs.green);
%       h2=plotm(clat,clon,'*','color',clrs.purple);
%       h3=plotm(olat,olon,'^','color',clrs.black);
%       %legend([h1,h2,h3],{'Tree-Rings','Corals','Others'},'location','southeast')
%    end 
%    colormap(cmap) % colormap from 'WhiteYellowOrangeRed'
%    set(gca,'fontsize',18)
%    %set(get(colorbar('location','southoutside'),'xlabel'),'string',sn{j},'fontsize',14)
%    set(get(colorbar('location','southoutside'),'xlabel'),'string','std units','fontsize',18)
%    caxis([-1.5 1.5])
%    %cbh = findobj( 0, 'tag', 'Colorbar' ); % and
%    %set(cbh,'YTick',[-1.5:0.25:1.5])
%    %delete( cbh )
%    hold off
% 
%    if  make_figs=='y'
%       print('-dpng',['./figs/som' num2str(somRows*somCols) '_ptrn_' num2str(i) '_' timestp2 '.png'],'-r300')
%    end
% 
% end
% 
% end
% 
% 
% % figure
% % for i=1:somRows*somCols
% %   subplot(somRows,somCols,i)
% %   h=pcolor(reshape(ptrns(i,:),length(lat),length(lon)));
% %   set(h,'edgecolor','none')
% % end
% % %ht=suptitle(['' fn(1:enddd-13)]);
% % %set(ht,'interpreter','none')
%  
%  
% %----------------------------------------------
% %  COMPARE SOM FREQUENCIES BY PERIOD
% %----------------------------------------------
% 
% 
% % Get BMUs for Mdroughts during 800-1600, for non-drought yrs 800-1600, all years 1601-1925
% 
% % BMU frequency during megadroughts (note droughts are determined analysis period relative)
% df=[];
% for i=1:length(ds0)
%   df0=ds0(i,1):ds0(i,2);
%   df=cat(2,df,df0);
% end
% dfc=histcounts(bmus(df),0.5:1:(somRows*somCols+0.5))/length(df);
% 
% % Frequency during non-droughts
% [~,ay1,~] = intersect(analyrs0,analyrs1);
% wf=setdiff(ay1(:),df(:));
% wfc=histcounts(bmus(wf),0.5:1:(somRows*somCols+0.5))/length(wf);
% 
% % Frequency during 1600-1925 period
% [~,ay2,~] = intersect(analyrs0,analyrs2);
% a2fc=histcounts(bmus(ay2),0.5:1:(somRows*somCols+0.5))/length(ay2);
%  
% % Frequency for each drought
% ds_pfrq=zeros(length(ds0),somRows*somCols);
% for i=1:length(ds0)
%   ds_pfrq(i,:)=histcounts(bmus(ds0(i,1):ds0(i,2)),0.5:1:(somRows*somCols+0.5))/length(ds0(i,1):ds0(i,2));
% end
% 
% disp(fln0)
% disp('Drought percent change of the frequency of the pattern')
% pd=(dfc-a2fc)./a2fc;
% disp(pd)
% 
% %----------------------------
% % PLOT BMUS FOR EACH PATTERN
% %----------------------------
% 
% C=load('./colormaps/MPL_BrBG.mat');
% 
% figure
% hold on
% % Mean freq during mega D
% b=bar(1:somRows*somCols,a2fc);
% b.FaceColor=[0.5 0.5 0.5];
% % Mean freq 1601-1925
% plot(1:somRows*somCols,dfc,'o','linewidth',2,'color',C.cmap(1,:),'markerfacecolor',C.cmap(1,:))
% plot(1:somRows*somCols,wfc,'^','linewidth',2,'color',C.cmap(end,:))
% legend('1601-1925','Megadrought 800-1600','Non-drought 800-1600','location','south')
% %title('Ptrn. Freq. MDrought')
% xlabel('Pattern number','fontsize',18)
% ylabel('Frequency','fontsize',18)
% xticks(1:somRows*somCols)
% set(gca,'fontsize',16)
% box on
% if  make_figs=='y'
%    print('-dpng',['./figs/som' num2str(somRows*somCols) '_ptrnfrq_' timestp2 '.png'],'-r300')
% end
% 
% 
% 
% %figur
% %plot(analyrs0,bmus,'or')
% %ylabel('Pattern')
% %yticks(1:6)
% %xlabel('Year')
% %box on
% %if  make_figs=='y'
% %   print('-dpng',['./figs/ptrn_frq_timeseries' timestp2 '.png',])
% %end
% 
% % Model the transition frequencies as Markov chain
% 
% % Compute transition matrix
% %[T1] = transM(bmus,1);
% 
% % Markov chain object
% %mc = dtmc(T1);
% 
% % Directed graph with color = probabilities
% %figure
% %graphplot(mc,'ColorEdges',true);
% 
% 
% 
% %----------------------------------------------
% %  COMPOSITES OF PDSI OVER THE SOMS
% %----------------------------------------------
% 
% if plot_sp=='y'
% 
% % Detrend data for composites 
% D0=reshape(Xpdsi(:,:,ir0),length(lat)*length(lon),length(ir0));
% dt_pdsi=detrend(D0')';% Detrend the data, if desired
% 
% 
% for i=1:somRows*somCols
% 
%    cind=find(bmus==i);
%    varX=reshape(mean(dt_pdsi(:,cind),2),length(lat),length(lon));
% 
%    figure
%    hold on
%    s=load('coast');mlat=s.lat;mlon=s.long;
%    %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
%    %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[120 33]);
%    h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
%    varZ=cat(2,varX,varX(:,1));
%    pcolorm(lat,[lon; 360],varZ);
%    title(['Pattern ' num2str(i)],'fontsize',18,'Interpreter','none')
%    %cptcmap('GMT_red2green')
%    %cptcmap('GMT_drywet')
%    %C=importdata('MPL_BrBG.txt');
%    C=load('./colormaps/MPL_BrBG.mat');
%    colormap(C.cmap)
%    caxis([-1.25 1.25])
%    set(gca,'fontsize',16)
%    set(get(colorbar('location','southoutside'),'xlabel'),'string','PDSI','fontsize',16)
%    cbh = findobj( 0, 'tag', 'Colorbar' );
%    %set(cbh,'YTick',[-1:0.2:1])
%    plotm(mlat,mlon,'color',[0.5 0.5 0.5])
%    tightmap
% 
%    if make_figs=='y';print('-dpng',['./figs/cmpst_pdsi_som' num2str(somRows*somCols) '_ptrn_' num2str(i) '_' timestp2 '.png'],'-r300'); end
% end
% 
% end
% 







