tutorial_reentrant_channel/input/gendata.5km.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 eddying resolution: dx=dy=5 km (Cartesian)
                0007 # to exactly match the bathymetry and forcing in the coarse-res setup
                0008
                0009 # grid depths generated Using the hyperbolic tangent method of
                0010 # Stewart et al. (2017) DOI: 10.1016/j.ocemod.2017.03.012
                0011 # to design an optimal grid.
                0012 # https://github.com/kialstewart/vertical_grid_for_ocean_models
                0013
                0014 dr = np.array([5.48716549,   6.19462098,   6.99291201,   7.89353689, \
                0015                8.90937723,  10.05483267,  11.34595414,  12.80056778, \
                0016               14.43837763,  16.28102917,  18.35210877,  20.67704362, \
                0017               23.28285446,  26.1976981 ,  29.45012046,  33.06792588, \
                0018               37.07656002,  41.496912  ,  46.34247864,  51.61592052, \
                0019               57.30518684,  63.37960847,  69.78661289,  76.44996107, \
                0020               83.27047568,  90.13003112,  96.89898027, 103.44631852, \
                0021              109.65099217, 115.4122275 , 120.65692923, 125.34295968, \
                0022              129.45821977, 133.01641219, 136.05088105, 138.60793752, \
                0023              140.74074276, 142.50436556, 143.95220912, 145.133724  , \
                0024              146.09317287, 146.86917206, 147.49475454, 147.99774783, \
                0025              148.40131516, 148.72455653, 148.98310489, 149.18968055, \
                0026              149.35458582])
                0027 nx = 200
                0028 ny = 400
                0029 nr = len(dr)
                0030 rF = (np.insert(np.cumsum(dr), 0, 0)) # z-coordinates of vertical cell faces
                0031 z = np.diff(rF)/2 + rF[:-1]           # z-coordinates of vertical cell centers
                0032 H = -np.sum(dr)                       # max depth of vertical grid
                0033
                0034 # bathymetry -- flat bottom of depth H (m) with idealized mid-depth ridge
                0035 bump_max = 2000.   # peak height of ridge above flat bottom depth
                0036 bathy = H * np.ones([ny, nx])
                0037 bump=np.zeros([ny, nx])
                0038 # sinusoidal bump running N-S through middle of domain
                0039 # this is hard-coded for nx=200, ny=400 resolution
                0040 r1 = bump_max * sin(np.linspace(0,pi,80))
                0041 r2 = np.reshape((np.linspace(0,1,52)), [52,1])  # create linear ramp for center notch
                0042 bump[:,55:135] = r1
                0043 # linearly lower bump height toward center notch
                0044 bump[134:186,:] = bump[134:186,:] * np.flip(r2[:])
                0045 bump[204:256,:] = bump[204:256,:] * r2[:]
                0046 bump[185:205,:] = 0.0  # notch; in these latitude bands, contours of f/H are unblocked
                0047 bathy = bathy + bump;
                0048 bathy[0:10,:] = 0.                    # wall at southern boundary, matching coarse-res
                0049 bathy.astype('>f4').tofile('bathy.5km.bin')
                0050
                0051 # zonal wind stress file  (matched to coarse-res profile)
                0052 taux_max = 0.2
                0053 taux = np.zeros([ny, nx])             # at (XG,YC) points
                0054 taux[5:395,:] = np.reshape(taux_max * sin(np.linspace(0,pi,390)), [390,1])
                0055 taux.astype('>f4').tofile('zonal_wind.5km.bin')
                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 = np.zeros([nr, ny, nx])
                0061 rbcs_mask[1:,-10:,:] = 1.0
                0062 rbcs_mask[1:,-20:-10,:] = 0.25
                0063 rbcs_mask.astype('>f4').tofile('T_relax_mask.5km.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.05,39.95,ny), [ny,1]) # at (XC,YC) points
                0070 sst_relax.astype('>f4').tofile('SST_relax.5km.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.5km.bin')