%code to read hycom .a files and create plots of speed within the layers that contain
%the depth 100mab (not interpolated to 100mab because these deep layers are usually 
%thick slabs with the 100mab depth being in the lowest layer).

%directory with HYCOM .a files
hdir='/hycom/ftp/datasets/GOMl0.04/expt_50.1/data/';

%%%% dimensions -- these can be read from a .b file
IDM=541;
JDM=385;
KDM=36;

IJDM=IDM*JDM;
npad=4096-mod(IJDM,4096);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gridfile='//hycom/ftp/datasets/GOMl0.04/expt_50.1/topo/regional.grid.a';
topofile='/hycom/ftp/datasets/GOMl0.04/expt_50.1/topo/depth_GOMl0.04_03i.a';


grid_fid=fopen(gridfile,'r');
depth_fid=fopen(topofile,'r');

[plon,count]=fread(grid_fid,IJDM,'float32','ieee-be');
fseek(grid_fid,4*(npad+IJDM),-1);
[plat,count]=fread(grid_fid,IJDM,'float32','ieee-be');

plon=reshape(plon,IDM,JDM);
plon=plon(:,1);
plat=reshape(plat,IDM,JDM);
plat=plat(1,:);

fseek(depth_fid,6*4*(npad+IJDM),-1);
[depth,count]=fread(depth_fid,IJDM,'float32','ieee-be');
y=find(depth>1e10);
depth(y)=nan;
depth=reshape(depth,IDM,JDM)';

fclose(grid_fid);
fclose(depth_fid);

%variables for vectorizing 100mab computations
depth3d=permute(repmat(depth,[1 1 KDM]),[3 1 2]);
[ii jj]=meshgrid([1:IDM],[1:JDM]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% field indices in the archv file%%%%%%%%%%
%% Start from zero
%% These are hard-coded here but can be read from the .b file

SSH_INDEX=1;
UBARO_INDEX=7;
VBARO_INDEX=8;
UVEL_INDEX=9;
VVEL_INDEX=10;
LT_INDEX=11;
TEMP_INDEX=12;
SALT_INDEX=13;

nfields=8; %number of fields per layer

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LOOP THROUGH Years, DAYS, HOURS
%%%%%%
day1=datenum(1993,1,1,0,0,0);
day2=datenum(2013,12,31,0,0,0);

for day=day1:day2;

  hr=0;

  [year,month,mday,hrtmp,mi,sc]=datevec(day);
  yday=day-datenum(year,1,1,0,0,0)+1;
  sday=num2str(yday,'%3.3i');
  shr=num2str(hr,'%2.2i');
  syear=num2str(year,'%4.4i');
  archfile= [hdir syear '/501_archv.' syear '_' sday '_' shr '.a'];
  disp(archfile)
  
  if exist(archfile)

   archv_fid=fopen(archfile,'r');
 
   fseek(archv_fid,(UBARO_INDEX)*4*(npad+IJDM),-1);
   [a,count]=fread(archv_fid,IJDM,'float32','ieee-be');
   y=find(a>1e10);
   a(y)=nan;
   a=reshape(a,IDM,JDM)';
   ubaro(:,:)=a;

   fseek(archv_fid,(VBARO_INDEX)*4*(npad+IJDM),-1);
   [a,count]=fread(archv_fid,IJDM,'float32','ieee-be');
   y=find(a>1e10);
   a(y)=nan;
   a=reshape(a,IDM,JDM)';
   vbaro(:,:)=a;

   fseek(archv_fid,(SSH_INDEX)*4*(npad+IJDM),-1);
   [a,count]=fread(archv_fid,IJDM,'float32','ieee-be');
   y=find(a>1e10);
   a(y)=nan;
   a=reshape(a,IDM,JDM)';
   ssh(:,:)=a/9.806; %SSH stored internally scaled by g

   for k=0:KDM-1
     fseek(archv_fid,(k*nfields+UVEL_INDEX)*4*(npad+IJDM),-1);
     [a,count]=fread(archv_fid,IJDM,'float32','ieee-be');
     y=find(a>1e10);
     a(y)=nan;
     a=reshape(a,IDM,JDM)';
     uvel(k+1,:,:)=a+ubaro;
   end

   for k=0:KDM-1
     fseek(archv_fid,(k*nfields+VVEL_INDEX)*4*(npad+IJDM),-1);
     [a,count]=fread(archv_fid,IJDM,'float32','ieee-be');
     y=find(a>1e10);
     a(y)=nan;
     a=reshape(a,IDM,JDM)';
     vvel((k+1),:,:)=a+vbaro;
   end

   for k=0:KDM-1
     fseek(archv_fid,(k*nfields+LT_INDEX)*4*(npad+IJDM),-1);
     [a,count]=fread(archv_fid,IJDM,'float32','ieee-be');
     y=find(a>1e10);
     a(y)=nan;
     a=reshape(a,IDM,JDM)';
     th((k+1),:,:)=a/9806.;
   end

   fclose(archv_fid);
   %co-locate u/v
   %for Hij, surrounding upts are ui,j and ui+1,j.  vpts are vi,j-1 and vi,j
   uvel(isnan(uvel))=0;
   vvel(isnan(vvel))=0;
   th2=th;
   th2(isnan(th))=0;
   thu=th2;
   thu(:,:,2:end)=0.5*(th2(:,:,1:end-1)+th2(:,:,2:end));
   thv=th2;
   thv(:,2:end,:)=0.5*(th2(:,1:end-1,:)+th2(:,2:end,:));
   u=uvel*nan;
   u(:,:,1:end-1)=0.5.*(uvel(:,:,1:end-1).*thu(:,:,1:end-1)+uvel(:,:,2:end).*thu(:,:,2:end))./th(:,:,1:end-1);
   v=vvel*nan;
   v(:,1:end-1,:)=0.5.*(vvel(:,1:end-1,:).*thv(:,1:end-1,:)+vvel(:,2:end,:).*thv(:,2:end,:))./th(:,1:end-1,:);

   
   layer_int=cumsum(th,1); %layer interfaces

   kb=squeeze(sum(depth3d-layer_int>100,1))+1;
   inds=sub2ind(size(u),kb,jj,ii);
   ubot=u(inds);
   vbot=v(inds);
   ubot(depth<500)=nan;
   vbot(depth<500)=nan;
   spdbot=sqrt(ubot.*ubot+vbot.*vbot);
  
   clf;
   pcolor(plon,plat,spdbot); shading interp; colorbar; title(datestr(day));
   disp('press a key to continue...')
   pause;

  end %if exist filename

end %day
