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