tutorial_reentrant_channel/input/gendata.50km.py

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