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 |