tutorial_reentrant_channel/input/gendata.5km.m

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