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