% initialize a few things
nx=8640; ny=4320;                          % model grid dimensions
dt=1/24;                                   % model output time step (days)
dte=datenum(2011,3,6):dt:datenum(2013,4,22,6,0,0); % model output times
pin='~dmenemen/llc_2160/regions/latlon/'; % model output files location
suf=['_' int2str(nx) 'x' int2str(ny)];
for fld={'XC','YC','XG','YG','hFacC'}      % load grid locations and mask
  fnm=[pin 'grid/' fld{1} suf];
  eval([fld{1} '=single(readbin(fnm,[nx ny]));'])
end
xc=XC(:,500);                              % longitude of grid centers
xg=XG(:,500);                              % longitude of grid corners
xc(find(xc<xc(1)))=xc(find(xc<xc(1)))+360;
xg(find(xg<xg(1)))=xg(find(xg<xg(1)))+360;
yc=YC(6000,:);                             % latitude of grid centers
yg=YG(6000,:);                             % latitude of grid corners
suf=['_' int2str(nx) 'x' int2str(ny) 'x1.'];

% Cameroon region indices
ix1=find(xc>-10.1 & xc<13.1);
iy1=find(yc>-9.1 & yc<8.1);
ix2=find(xc>-10.0 & xc<13.0);
iy2=find(yc>-9.0 & yc<8.0);

% plot hourly velocity in Cameroon region
for t=datenum(2012,1,1):dt:datenum(2013,1,1)
  fnu=[pin 'U/U' suf datestr(t,30)];
  tmp=readbin(fnu,[nx ny]);
  U=tmp(ix1,iy1);                        % zonal velocity (m/s)
  fnv=[pin 'V/V' suf datestr(t,30)];
  tmp=readbin(fnv,[nx ny]);
  V=tmp(ix1,iy1);                        % meridional velocity (m/s)
  u=interp2(yc(iy1),xg(ix1),U,YC(ix2,iy2),XC(ix2,iy2));
  v=interp2(yg(iy1),xc(ix1),V,YC(ix2,iy2),XC(ix2,iy2));
  u(481:end,163:end)=0;
  v(481:end,163:end)=0;
  clf
  subplot(211), mypcolor(u'); caxis([-1 1]/2); thincolorbar
  title(['uVel on ' datestr(t)])
  subplot(212), mypcolor(v'); caxis([-1 1]/2); thincolorbar
  title(['vVel on ' datestr(t)])
  pause(.1)
end

% experiment with quiver plots
%>>>> change t for model snapshot
t=datenum(2012,1,1);
fnu=[pin 'U/U' suf datestr(t,30)];
tmp=readbin(fnu,[nx ny]);
U=tmp(ix1,iy1);                        % zonal velocity (m/s)
fnv=[pin 'V/V' suf datestr(t,30)];
tmp=readbin(fnv,[nx ny]);
V=tmp(ix1,iy1);                        % meridional velocity (m/s)
u=interp2(yc(iy1),xg(ix1),U,YC(ix2,iy2),XC(ix2,iy2));
v=interp2(yg(iy1),xc(ix1),V,YC(ix2,iy2),XC(ix2,iy2));
u(481:end,163:end)=nan;
v(481:end,163:end)=nan;
u(find(~u))=nan;
v(find(~v))=nan;

% test 1
clf
myquiver(XC(ix2(5:10:end),iy2(5:10:end)),YC(ix2(5:10:end),iy2(5:10:end)), ...
         u(5:10:end,5:10:end),v(5:10:end,5:10:end),4,'k')  
plotland('k',12,double(XC(ix2,1)),double(YC(1,iy2))')
axis([-10 13 -9 8])

% test 2
clf
mag=sqrt(u.^2+v.^2);
mypcolor(XC(ix2,1),YC(1,iy2),mag');
%>>>> change cm for colormap
cm=cmap;
cm=jet;
cm=hot;
cm=gray;
cm(1,:)=1;
colormap(cm)
%>>>> change cx for inensity
cx=.3;
caxis([0 cx])
h=thincolorbar;
hold on
%>>>> change sc for length of arrows
sc=2;
myquiver(XC(ix2(5:10:end),iy2(5:10:end)),YC(ix2(5:10:end),iy2(5:10:end)), ...
         u(5:10:end,5:10:end),v(5:10:end,5:10:end),sc,'k')  
xt=get(h,'ytick');
set(h,'YLim',[cx/length(cm) cx],'ytick',[cx/length(cm) xt(2:end)],'yticklabel',xt)

% compute monthly mean 2012 fields
Umonthly=zeros(length(ix1),length(iy1),12);
Vmonthly=zeros(length(ix1),length(iy1),12);
for m=1:12
  n=0;
  for t=datenum(2012,m,1):dt:datenum(2012,m+1,1)
    disp(datestr(t))
    n=n+1;
    fnu=[pin 'U/U' suf datestr(t,30)];
    tmp=readbin(fnu,[nx ny]);
    Umonthly(:,:,m)=Umonthly(:,:,m)+tmp(ix1,iy1);
    fnv=[pin 'V/V' suf datestr(t,30)];
    tmp=readbin(fnv,[nx ny]);
    Vmonthly(:,:,m)=Vmonthly(:,:,m)+tmp(ix1,iy1);
  end
  Umonthly(:,:,m)=Umonthly(:,:,m)/n;
  Vmonthly(:,:,m)=Vmonthly(:,:,m)/n;
end
lon=xc(ix1); lat=yc(iy1);
fnm=[pin 'matlab/CameroonMonthly2012'];
eval(['save ' fnm ' Umonthly Vmonthly lon lat'])

%%%%%%%%%%%%%%
% plot monthly-mean fields
clear all
cd /data1/llc/llc_2160/regions/latlon/matlab
load CameroonMonthly2012
u=sum(abs(Umonthly)+abs(Vmonthly),3);
in=find(~u);
LON=lon*ones(1,length(lat));
LAT=ones(length(lon),1)*lat;
cm=cmap;
cm(1,:)=.8*[1 1 1];
colormap(cm)
cx=.5;
sc=2.5;
for m=1:12
  u=Umonthly(:,:,m);
  v=Vmonthly(:,:,m);
  mag=sqrt(u.^2+v.^2);
  mag(in)=-cx/length(cm);
  clf
  mypcolor(lon,lat,mag');
  hold on
  caxis([-cx/length(cm) cx]);
  u(in)=nan;
  v(in)=nan;
  myquiver(LON(5:10:end,5:10:end),LAT(5:10:end,5:10:end), ...
           u(5:10:end,5:10:end),v(5:10:end,5:10:end),sc,'k')
  h=thincolorbar; set(h,'YLim',[0 cx])
  title(datestr(datenum(2012,m,15),'mmm yyyy'));
  pause(1)
end