MITgcm_contrib/tides/create_bc_grid.m

% 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 
1	dimitri	1.1	% function create_bc_grid(dir_in, nx, ny, nz, cbdry)
2			%
3			% Input: dir_in, directory that has MITgcm grid files
4			% (following files are expected,
5			% YC.data. XC.data, hFacC.data
6			% RF.data, AngleCS.data, AngleSN.data )
7			% nx, first dimension of MITgcm grid, (west-east)
8			% ny, second dimension of MITgcm grid (south-north)
9			% nz, dimension in z (or r) direction
10			% cbdry, character for bdry to generate: s, n, w, e, (south, north, west, east)
11			%
12			% Output: No output
13			% The script will generate a grid file named
14			% bdry_grid_[snwe].mat
15			%
16			%
17			% XC WANG 25jan2013
18			% Modified 01feb2013
19			%
20			%
21			function create_bc_grid(dir_in,nx,ny,nz,cbdry)
22
23			% addpath /data17/home/xiao/Mat_Tools/MITgcm_Tool
24
25			% Default location of boundary
26			ny_north = ny -1 ; ny_south=1;
27			nx_west = 1; nx_east=nx-1 ;
28
29			%dir_in=['/data17/XCW_output/Grid/data/'];
30			% Three boundaries, South, North and West
31			%ibc = 2;
32			%c_bdry3 = {'s', 'n', 'w', 'e'} ;
33			%c_bdry = c_bdry3{ibc} ;
34			%%bdry_file=['bdry_grid_',c_bdry, '.mat.test'];
35			%bdry_file=['bdry_grid_',c_bdry, '.mat'];
36
37			ufile = ['YC.data'];
38			ufile1 = [dir_in, ufile];
39			lat = readbin(ufile1, [nx ny]) ;
40
41			ufile = ['XC.data'];
42			ufile1 = [dir_in, ufile];
43			lon = readbin(ufile1, [nx ny]) ;
44
45			%ufile = ['Depth.data'];
46			%ufile1 = [dir_in, ufile];
47			%depth = readbin(ufile1, [nx ny]) ;
48
49			rf=-readbin([dir_in, 'RF.data'], nz+1) ;
50			thk=diff(rf) ;
51
52			hfac =readbin([dir_in, 'hFacC.data'], [nx ny nz]) ;
53			%obmask = squeeze(hfac(:, ny_north, :)) ;
54			thk3d = zeros(size(hfac)) ;
55			for k =1:nz
56			temp1 = squeeze(hfac(:,:, k)) ;
57			thk3d(:,:, k) = thk(k) * temp1 ;
58			end
59			depth = sum(thk3d, 3) ;
60			%depth_bc = depth2d(:, ny_north) ;
61
62			anglefile = ['AngleCS.data'] ;
63			tfile1 = [dir_in, anglefile];
64			AngleCS = readbin(tfile1, [nx ny]) ;
65
66			anglefile = ['AngleSN.data'] ;
67			tfile1 = [dir_in, anglefile];
68			AngleSN = readbin(tfile1, [nx ny]) ;
69
70			%%bdry_file=['bdry_grid_',c_bdry, '.mat.test'];
71			bdry_file=['bdry_grid_',cbdry, '.mat'];
72			switch cbdry
73			case 's'
74			lat_bc=lat(:, ny_south);
75			lon_bc=lon(:, ny_south) ;
76			depth_bc=depth(:, ny_south) ;
77			AngleCS_bc = AngleCS(:,ny_south);
78			AngleSN_bc = AngleSN(:,ny_south);
79			%bdry_file=['bdry_grid_',c_bdry, '.mat'];
80			case 'n'
81			lat_bc=lat(:, ny_north);
82			lon_bc=lon(:, ny_north) ;
83			depth_bc=depth(:, ny_north) ;
84			AngleCS_bc = AngleCS(:,ny_north);
85			AngleSN_bc = AngleSN(:,ny_north);
86			case 'w'
87			lat_bc=lat(nx_west,:);
88			lon_bc=lon(nx_west,:) ;
89			depth_bc=depth(nx_west,:) ;
90			AngleCS_bc = AngleCS(nx_west,:);
91			AngleSN_bc = AngleSN(nx_west,:);
92			case 'e'
93			lat_bc=lat(nx_east,:);
94			lon_bc=lon(nx_east,:) ;
95			depth_bc=depth(nx_east,:) ;
96			AngleCS_bc = AngleCS(nx_east,:);
97			AngleSN_bc = AngleSN(nx_east,:);
98			otherwise
99			display('No such selection ', cbdry)
100			return
101			end
102
103			%rad2deg = 180/pi;
104			%angle = atan2(AngleSN, AngleCS);
105			%angle = angle*rad2deg ;
106
107			save(bdry_file, 'lat_bc', 'lon_bc', 'depth_bc', 'AngleCS_bc', 'AngleSN_bc') ;
108
109			return