function STATM = speedy_stats(varargin)    %varargin=(var2,0,360,-90,90)';% glb

var=varargin{1};  % var2 (var=reshape(a,96,48,2) where 1 model eof1 and 2 is ncep eof1)

lon1=varargin{2}; % 0
lon2=varargin{3}; % 360
lat1=varargin{4}; % -90
lat2=varargin{5}; % 90

lat=[-87.159   -83.479 -79.777 -76.070 -72.362 -68.652 -64.942 -61.232   -57.521   -53.810   -50.099   -46.389   -42.678   -38.967 -35.256   -31.545   -27.833   -24.122   -20.411   -16.700   -12.989    -9.278  -5.567    -1.856     1.856     5.567     9.278    12.989    16.700    20.411  24.122    27.833    31.545    35.256    38.967    42.678    46.389    50.099 53.810    57.521    61.232    64.942    68.652    72.362    76.070    79.777  83.479    87.159];
lon=0:3.75:(360-3.75);

latinds=find((lat >= lat1)&(lat <= lat2));
loninds=find((lon >= lon1)&(lon <= lon2));
nlat=size(latinds,2);
nlon=size(loninds,2);

weightfun1=1.0e+05*[0.0847  0.1948  0.3055  0.4147  0.5222  0.6275  0.7301  0.8297  0.9259  1.0181  1.1059  1.1892  1.2677  1.3407  1.4081  1.4698  1.5251  1.5739  1.6162  1.6517  1.6804  1.7019  1.7164  1.7238  1.7238  1.7164  1.7019  1.6804  1.6517  1.6162  1.5739  1.5251  1.4698  1.4081  1.3407  1.2677  1.1892  1.1059  1.0181  0.9259  0.8297  0.7301  0.6275  0.5222  0.4147  0.3055  0.1948  0.0847];

for k=1:96;
  area(k,:) = weightfun1;
end;
areatot=sum(sum(area(loninds,latinds)));
N=nlon*nlat;
  
%ncep
var_aave(2)=sum(sum(var(loninds,latinds,1).*area(loninds,latinds)))/areatot;
vec(2,:)=reshape(var(loninds,latinds,1),nlon*nlat,1);
areavec=reshape(area(loninds,latinds),nlon*nlat,1);

st(2) = sqrt(sum(  areavec'.*(vec(2,:)-var_aave(1) ).^2)  /areatot);
rms(2) = sqrt(sum(  areavec'.*( ( vec(2,:)-var_aave(2) )-( vec(2,:)-var_aave(2) )).^2)  /areatot);
me(2) = var_aave(2);
co(2) = sum(  areavec'.*( ( vec(2,:)-var_aave(2) ).*( vec(2,:)-var_aave(2) )))/st(2)/st(2)/areatot;

%ncep
for k=1
  var_aave(1)=sum(sum(var(loninds,latinds,0).*area(loninds,latinds)))/areatot;
  vec(1,:)=reshape(var(loninds,latinds,0),nlon*nlat,1);
  
  st(1) = sqrt(sum(  areavec'.*(vec(1,:)-var_aave(1) ).^2)  /areatot);
  me(1) = var_aave(1);
  rms1k) = sqrt(sum(  areavec'.*( ( vec(1,:)-var_aave(1) )-( vec(2,:)-var_aave(2) )).^2) /areatot);
  co(1) = sum(  areavec'.*( ( vec(1,:)-var_aave(1) ).*( vec(2,:)-var_aave(2) )))/st(1)/st(2)/areatot;  
%  rmsenorm(1) = sqrt(sum(  areavec'.*( ( vec(1,:)-var_aave(1) )-( vec(2,:)-var_aave(2) ) ).^2  ) /areatot)/sqrt(st(1).^2+st(2).^2);
end;

STATM(1,:) = me;
STATM(2,:) = st;
STATM(3,:) = rms;
STATM(4,:) = co;
%STATM(5,:) = rmsenorm;