cd ../
startup
cd HYDRO
skip1=0;
if skip1
skip=0;
if skip
load /home/abatz/PROCESSING/PROCESSNARR/tmax4kmprism
tmax=permute(tmax,[1 2 4 3]);
load /home/abatz/PROCESSING/PROCESSNARR/tmin4kmprism
tmin=permute(tmin,[1 2 4 3]);
load /home/abatz/PROCESSING/PROCESSNARR/ppt4kmprism
ppt=permute(ppt,[1 2 4 3]);ppt=ppt/100;
cd ../
[x,y]=latlonPRISM;
y=flipdim(y,1);
cd HYDRO
load us_25m.dem
elev=flipdim(us_25m,1);clear us_25m;
load prismmonthlyradwinds
dwsrfprism=flipdim(dwsrfprism,1);
vprism=flipdim(vprism,1);
s1=size(tmax);
PET=PETSTUFF(dwsrfprism,vprism,tmax,tmin,y,x,elev);
save PETUS -v7.3 PET
disp('Computing Monthly Water Balance')
end
load /home/abatz/PROCESSING/PROCESSNARR/ppt4kmprism
ppt=permute(ppt,[1 2 4 3]);
load PETUS
WB = ppt - PET;
clear ppt PET
f=find((WB(:,:,1))>-9999);
% time scales of interest 1-12,15,18,24,30,36,48,60,72
sz=size(WB);
WB=shiftdim(WB,2);
WB=WB(:,:,f);
WB=shiftdim(WB,2);
ds=[1:12 15 18 24 30 36 48 60 72];
for i=1:20
for k=1:5
s1=1+100000*(k-1):100000*k;
if max(s1)>size(WB,1) s1=1+100000*(k-1):size(WB,1);end
spei(s1,:,:)=calculateSPI(WB(s1,:,:),ds(i));
end
save(['SPEI_',num2str(ds(i))],'-v7.3','spei');
clear spei
end

disp('Loop Complete');
%%
%PLOTS
end
for i=1:20
ds=[1:12 15 18 24 30 36 48 60 72];
load floc2
DATA=NaN*single(ones(12,116,621,1405));
load(['SPEI_',num2str(ds(i))])
DATA(:,:,f)=shiftdim(spei,1);clear spei;
DATA=shiftdim(DATA,2);
save(['/storage/DATA/DROUGHT/SPEI_',num2str(ds(i))],'-v7.3','DATA');
end