% function  create_bc_grid(dir_in, nx, ny, nz, cbdry)
% 
% Input:  dir_in, directory that has MITgcm grid files 
%         (following files are expected, 
%          YC.data. XC.data, hFacC.data
%          RF.data, AngleCS.data, AngleSN.data ) 
%         nx, first dimension of MITgcm grid, (west-east)  
%         ny, second  dimension of MITgcm grid (south-north) 
%         nz, dimension in z (or r)  direction 
%         cbdry, character for bdry to generate: s, n, w, e, (south, north, west, east)
%
% Output: No output 
%         The script will generate a grid file named
%         bdry_grid_[snwe].mat
%
% 
% XC WANG 25jan2013
% Modified 01feb2013 
% 
% 
function  create_bc_grid(dir_in,nx,ny,nz,cbdry)

% addpath /data17/home/xiao/Mat_Tools/MITgcm_Tool   

% Default location of boundary  
ny_north = ny -1 ; ny_south=1; 
nx_west = 1; nx_east=nx-1 ; 

%dir_in=['/data17/XCW_output/Grid/data/']; 
% Three boundaries, South, North and West   
%ibc = 2; 
%c_bdry3 = {'s', 'n', 'w', 'e'} ;  
%c_bdry = c_bdry3{ibc} ; 
%%bdry_file=['bdry_grid_',c_bdry, '.mat.test']; 
%bdry_file=['bdry_grid_',c_bdry, '.mat']; 

ufile = ['YC.data'];  
ufile1 = [dir_in, ufile]; 
lat = readbin(ufile1, [nx ny]) ;  

ufile = ['XC.data'];  
ufile1 = [dir_in, ufile]; 
lon = readbin(ufile1, [nx ny]) ; 

%ufile = ['Depth.data'];  
%ufile1 = [dir_in, ufile]; 
%depth = readbin(ufile1, [nx ny]) ;     

rf=-readbin([dir_in, 'RF.data'], nz+1) ;
thk=diff(rf) ;

hfac =readbin([dir_in, 'hFacC.data'], [nx ny nz]) ;
%obmask = squeeze(hfac(:, ny_north, :)) ;  
thk3d = zeros(size(hfac)) ;
for k =1:nz  
 temp1 = squeeze(hfac(:,:, k))  ; 
 thk3d(:,:, k) = thk(k) * temp1 ; 
end  
depth = sum(thk3d, 3) ;  
%depth_bc = depth2d(:, ny_north) ; 

anglefile = ['AngleCS.data'] ;  
tfile1 = [dir_in, anglefile]; 
AngleCS = readbin(tfile1, [nx ny]) ;    

anglefile = ['AngleSN.data'] ;  
tfile1 = [dir_in, anglefile]; 
AngleSN = readbin(tfile1, [nx ny]) ;      

%%bdry_file=['bdry_grid_',c_bdry, '.mat.test']; 
bdry_file=['bdry_grid_',cbdry, '.mat'];  
switch cbdry 
 case 's' 
  lat_bc=lat(:, ny_south); 
  lon_bc=lon(:, ny_south) ;
  depth_bc=depth(:, ny_south) ; 
  AngleCS_bc = AngleCS(:,ny_south); 
  AngleSN_bc = AngleSN(:,ny_south);   
%bdry_file=['bdry_grid_',c_bdry, '.mat'];  
 case 'n' 
  lat_bc=lat(:, ny_north); 
  lon_bc=lon(:, ny_north) ;
  depth_bc=depth(:, ny_north) ; 
  AngleCS_bc = AngleCS(:,ny_north); 
  AngleSN_bc = AngleSN(:,ny_north);    
 case 'w' 
  lat_bc=lat(nx_west,:); 
  lon_bc=lon(nx_west,:) ;
  depth_bc=depth(nx_west,:) ; 
  AngleCS_bc = AngleCS(nx_west,:); 
  AngleSN_bc = AngleSN(nx_west,:);     
 case 'e' 
  lat_bc=lat(nx_east,:); 
  lon_bc=lon(nx_east,:) ;
  depth_bc=depth(nx_east,:) ; 
  AngleCS_bc = AngleCS(nx_east,:); 
  AngleSN_bc = AngleSN(nx_east,:);     
 otherwise 
  display('No such selection ', cbdry)   
 return 
end 

%rad2deg = 180/pi; 
%angle = atan2(AngleSN, AngleCS); 
%angle = angle*rad2deg ; 

save(bdry_file, 'lat_bc', 'lon_bc', 'depth_bc', 'AngleCS_bc', 'AngleSN_bc') ;  

return