Warning, /verification/tutorial_reentrant_channel/input/gendata.50km.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 62748458 on 2021-08-11 00:15:22 UTC
62748458e3 Jeff*0001 % generate support data files for tutorial Southern Ocean Reentrant Channel
0002 % hard-coded for coarse (non-eddying) resolution: dx=dy=50 km (Cartesian)
0003
0004 % grid depths generated Using the hyperbolic tangent method of
0005 % Stewart et al. (2017) DOI: 10.1016/j.ocemod.2017.03.012
0006 % to design an optimal grid.
0007 % https://github.com/kialstewart/vertical_grid_for_ocean_models
0008
0009 dr = [ 5.48716549, 6.19462098, 6.99291201, 7.89353689, ...
0010 8.90937723, 10.05483267, 11.34595414, 12.80056778, ...
0011 14.43837763, 16.28102917, 18.35210877, 20.67704362, ...
0012 23.28285446, 26.1976981 , 29.45012046, 33.06792588, ...
0013 37.07656002, 41.496912 , 46.34247864, 51.61592052, ...
0014 57.30518684, 63.37960847, 69.78661289, 76.44996107, ...
0015 83.27047568, 90.13003112, 96.89898027, 103.44631852, ...
0016 109.65099217, 115.4122275 , 120.65692923, 125.34295968, ...
0017 129.45821977, 133.01641219, 136.05088105, 138.60793752, ...
0018 140.74074276, 142.50436556, 143.95220912, 145.133724 , ...
0019 146.09317287, 146.86917206, 147.49475454, 147.99774783, ...
0020 148.40131516, 148.72455653, 148.98310489, 149.18968055, ...
0021 149.35458582];
0022 nx = 20;
0023 ny = 40;
0024 nr = length(dr);
0025 rF = -[0 cumsum(dr)]; % z-coordinates of vertical cell faces
0026 z = rF(1:end-1) + diff(rF)/2; % z-coordinates of vertical cell centers
0027 H = -sum(dr); % max depth of vertical grid
0028
0029 % bathymetry -- flat bottom of depth H (m) with idealized mid-depth ridge
0030 bump_max = 2000; % peak height of ridge above flat bottom depth
0031 bathy = H .* ones(nx, ny);
0032 bump = zeros(nx, ny);
0033 % sinusoidal bump running N-S through middle of domain
0034 % this is hard-coded for nx=20, ny=40 resolution
0035 r1 = bump_max*repmat(sin(0:pi/8:pi),[ny 1])';
0036 r2 = 0:.2:1; % create linear ramp for center notch
0037 bump(6:14,:) = r1;
0038 % linearly lower bump height toward center notch
0039 bump(6:14,14:19) = bump(6:14,14:19) .* fliplr(repmat(r2, [9 1]));
0040 bump(6:14,21:26) = bump(6:14,21:26) .* repmat(r2, [9 1]);
0041 bump(6:14,19:21) = 0; % notch; in these lat bands, contours of f/H are unblocked
0042 bathy = bathy + bump;
0043 bathy(:,1) = 0; % wall at southern boundary
0044 fid=fopen('bathy.50km.bin', 'w', 'b');
0045 fwrite(fid, bathy, 'float32');
0046 fclose(fid);
0047
0048 % zonal wind stress file
0049 taux_max = 0.2;
0050 taux = repmat(taux_max * sin(0:pi/(ny-1):pi), [nx 1]); % at (XG,YC) points
0051 fid = fopen('zonal_wind.50km.bin', 'w', 'b');
0052 fwrite(fid, taux, 'float32');
0053 fclose(fid);
0054
0055 % 3-D mask for RBCS temperature relaxation
0056 % note we implement "sponge layer" with full restoring at N boundary row (y=40),
0057 % weaker restoring at row just south of N boundary (y=39)
0058 % mask set to zero in surface layer (where core model SST restoring applied)
0059 rbcs_mask = zeros(nx, ny, nr);
0060 rbcs_mask(:,end,2:end) = 1.0;
0061 rbcs_mask(:,end-1,2:end) = 0.25;
0062 fid=fopen('T_relax_mask.50km.bin', 'w', 'b');
0063 fwrite(fid, rbcs_mask, 'float32');
0064 fclose(fid);
0065
0066 % 2-D SST field for relaxation, linear ramp between Tmin and Tmax
0067 Tmax = 10.0;
0068 Tmin = -2.0;
0069 sst_relax = repmat(Tmin + (Tmax-Tmin)*(0.5:39.5)/40, [nx 1]); % at (XC,YC) points
0070 fid=fopen('SST_relax.50km.bin', 'w', 'b');
0071 fwrite(fid, sst_relax, 'float32');
0072 fclose(fid);
0073
0074 % 3-D Temperature field for initial conditions and RBCS northern wall profile
0075 h = 500; % e-folding scale for temperature decrease with depth
0076 T_surf = sst_relax; % use 2-D SST relaxation field for surface values
0077 zscale = permute(repmat((exp(z/h) - exp(H/h)) / (1-exp(H/h)), [nx 1 ny]), [ 1 3 2]);
0078 T_3D = repmat(T_surf - Tmin, [1 1 nr]) .* zscale + Tmin;
0079 fid=fopen('temperature.50km.bin', 'w', 'b');
0080 fwrite(fid, T_3D, 'float32');
0081 fclose(fid);
