%code for reading HYCOM file

% directory where the hycom .a output files exist
hdir='/hycom/ftp/datasets/GOMl0.04/expt_50.1/data/';

%paths to the model grid and topography
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';

%%%% dimensions  
IDM=541;
JDM=385;
KDM=36;

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% field indices in the archv file%%%%%%%%%%
%% These can be determined from the .b files

SSH_INDEX=1;
UBARO_INDEX=7; %barotropic velocity
VBARO_INDEX=8;
UVEL_INDEX=9;  %baroclinic velocity (need to add barotropic)
VVEL_INDEX=10;
LT_INDEX=11;   %layer thickness
TEMP_INDEX=12; %potential temperature
SALT_INDEX=13; %salinity

nfields=8; %number of fields per layer

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%read model grid and bathymetry

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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; %only read 00UTC files

  [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) %some files may be missing from the archive

   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 is stored 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; %add barotropic to get total velocity
   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; %add barotropic to get total velocity
   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.; %layer thicknesses are stored scaled by g rho0
   end

   for k=0:KDM-1
     fseek(archv_fid,(k*nfields+TEMP_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)';
     temp((k+1),:,:)=a; %potential temperature
   end

   for k=0:KDM-1
     fseek(archv_fid,(k*nfields+SALT_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)';
     salt((k+1),:,:)=a; %salinity
   end

   %Note, this experiment also saved momentum viscosity and diffusion terms, as well as montgomery
   %potential and various surface and mixed layer fields.  See .b file for variables.

   fclose(archv_fid);

   %Variables are saved on the C-grid.  Here's some fast code for co-locating
   %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;

   %first, interpolate layer thicknesses from grid cell centers to u/v points
   %then average transports u*h and v*h
   %finally, divide by the grid cell thickness

   th2=th;
   th2(isnan(th))=0;
   thu=th2;
   thu(:,:,2:end)=0.5*(th2(:,:,1:end-1)+th2(:,:,2:end));
   thv=th2;
   thv(:,1:end-1,:)=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(:,2:end,:)=0.5.*(vvel(:,1:end-1,:).*thv(:,1:end-1,:)+vvel(:,2:end,:).*thv(:,2:end,:))./th(:,2:end,:);

   

  end %if exist filename
end %day
