
%-----------------------------------------------
% SAVE SEASONAL CYCLE OF PRECIP FOR NETCDF FIGS
%-----------------------------------------------

% Easter island lats/lons: -27.116667, -109.366667
eic_lat=-27.116667;
eic_lon=360-109.366667; % 250.6333

dt=(datenum('01-Jan-1979'):datenum('31-Dec-2017'))';

% lf_t2m='/d1/nsteiger/climate-data/era5/era5_t2m_scpac_daily_1979_2017.nc';
% %lf_t2m='/d1/nsteiger/climate-data/era5/era5_t2m_swpac_daily_1979_2017.nc';
lf_t2m='/d1/nsteiger/climate-data/era5/era5_t2m_sh_daily_1979_2017.nc';
% t2m=ncread(lf_t2m,'t2m'); % larger area
% t2m=flip(permute(t2m,[2 1 3]),1);
% 
slat=double(ncread(lf_t2m,'latitude'));
%slat=flip(slat);
slon=double(ncread(lf_t2m,'longitude'));
% 
mslp=ncread('/d1/nsteiger/climate-data/era5/era5_msl_sh_daily_1979_2017.nc','msl'); % larger area
% %mslp=ncread('/d1/nsteiger/climate-data/era5/era5_msl_swpac_daily_1979_2017.nc','msl');
% mslp=ncread('/d1/nsteiger/climate-data/era5/era5_msl_scpac_daily_1979_2017.nc','msl'); % larger area
%mslp=flip(permute(mslp,[2 1 3]),1);
mslp=permute(mslp,[2 1 3]); % don't flip right-side up
% 
%pcp=ncread('/d1/nsteiger/climate-data/era5/era5_tp_sh_daily_1979_2017.nc','tp'); % larger area
% %pcp=ncread('/d1/nsteiger/climate-data/era5/era5_tp_swpac_daily_1979_2017.nc','tp');
%pcp=ncread('/d1/nsteiger/climate-data/era5/era5_tp_scpac_daily_1979_2017.nc','tp'); % larger area
%pcp=flip(permute(pcp,[2 1 3]),1);
%pcp=permute(pcp,[2 1 3]); % don't flip right-side up
%pcp(pcp<0)=0; % specify 0 for missing values

% Monthly mean of spatial fields
%[pcp_mn,mdt_mon]=mean_ts(pcp,dt,'nansum','month');
[pcp_mn,mdt_mon]=mean_ts(mslp,dt,'mean','month');

% mean of all months
xmn=mean(pcp_mn,3);

% Reshape data
k=1;l=12;
pcp_mx=zeros(size(pcp_mn,1),size(pcp_mn,2),12,length(mdt_mon)/12);
for n=1:length(mdt_mon)/12
    pcp_mx(:,:,:,n)=pcp_mn(:,:,k:l);
    k=k+12; l=l+12;
end

% Monthly mean
mn_pcp=permute(squeeze(mean(pcp_mx,4)),[2 1 3]);
%mn_pcp(mn_pcp==NaN)=1e20; % any fill values...

% Mean of monthly means
pcp_nomn=bsxfun(@minus,mn_pcp,permute(mean(pcp_mn,3),[2 1]));


% define input and output variable
inputf='era5_sh_latlon.nc';
%outputf=['era5_sh_pcp_climatology.nc'];
outputf=['era5_sh_mslp_climatology.nc'];

% grab lat/lon from model, put it into the output file
eval(['!ncks -A -v latitude,longitude ',inputf,' ',outputf]);


% Add time dimension
nccreate(outputf,'time','Dimensions',{'time',12});
ncwrite(outputf,'time',1:12);
ncwriteatt(outputf,'time','units','Monthly mean');
ncwriteatt(outputf,'time','long_name','Month');



for i=1:12

   mnpcp0=squeeze(mn_pcp(:,:,i));
   %cnm=['pcp_mon_' num2str(i)];
   cnm=['mslp_mon_' num2str(i)];
   
   %nccreate(outputf,cnm,'Dimensions',{'longitude',length(slon),'latitude',length(slat)},'FillValue',1e20);
   nccreate(outputf,cnm,'Dimensions',{'longitude',length(slon),'latitude',length(slat)});
   ncwrite(outputf,cnm,mnpcp0);
%   ncwriteatt(outputf,cnm,'missing_value',1e20);

   mnpcp1=squeeze(pcp_nomn(:,:,i));
   %cnm1=['pcp_mon_anom_' num2str(i)];
   cnm1=['mslp_mon_anom_' num2str(i)];
   
   nccreate(outputf,cnm1,'Dimensions',{'longitude',length(slon),'latitude',length(slat)});
   ncwrite(outputf,cnm1,mnpcp1);

end

% Add mean of all years
cnm1=['mslp_allmean'];

nccreate(outputf,cnm1,'Dimensions',{'longitude',length(slon),'latitude',length(slat)});
ncwrite(outputf,cnm1,permute(xmn,[2 1]));



% 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]);
 









