tutorial_barotropic_gyre/input/gendata.py

ce0d9af5ea Jeff*0001 import numpy as np
                0002 from numpy import cos, sin, pi
                0003
                0004 Ho = 5000  # ocean depth in meters
                0005 nx = 62    # number of gridpoints in x-direction
                0006 ny = 62    # number of gridpoints in y-direction
                0007 xo = 0     # origin in x,y for ocean domain
                0008 yo = 0     # (i.e. southwestern corner of ocean domain)
                0009 dx = 20    # grid spacing in x (km)
                0010 dy = 20    # grid spacing in y (km)
                0011 xeast  = xo + (nx-2)*dx   # eastern extent of ocean domain
                0012 ynorth = yo + (ny-2)*dy   # northern extent of ocean domain
                0013
                0014 # Flat bottom at z=-Ho
                0015 h = -Ho * np.ones((ny, nx))
                0016
                0017 # Walls (surrounding domain); generate bathymetry file
                0018 h[:, [0,-1]] = 0   # set ocean depth to zero at east and west walls
                0019 h[[0,-1], :] = 0   # set ocean depth to zero at south and north walls
                0020 # save as single-precision (float32) with big-endian byte ordering
                0021 h.astype('>f4').tofile('bathy.bin')
                0022
                0023 # Ocean domain extends from (xo,yo) to (xeast,ynorth)
                0024 # (i.e. the ocean spans nx-2, ny-2 grid cells)
                0025 # out-of-box-config: xo=yo=0, dx=dy=20 km, ocean extent (0,0)-(1200,1200) km
                0026 # model domain includes a land cell surrounding the ocean domain
                0027 # The full model domain cell centers are located at:
                0028 #    XC[0,:] = -10, +10, ..., +1210 (km)
                0029 #    YC[:,0] = -10, +10, ..., +1210 (km)
                0030 # and full model domain cell corners are located at:
                0031 #    XG[0,:] = -20, 0, ..., 1200 [, 1220] (km)
                0032 #    YG[:,0] = -20, 0, ..., 1200 [, 1220] (km)
                0033 # where the last value in brackets is not included in the MITgcm grid variable
                0034 # and reflects the eastern and northern edge of the model domain respectively.
                0035 # See section 2.11.4 of the MITgcm users manual.
                0036
                0037 # Zonal wind-stress, located at u-points (see section 2.11.4)
                0038 # here we non-dimensionalize: 0 at southern and western ocean boundary
                0039 # to 1.0 at eastern and northern ocean boundary
                0040 # for the purpose of applying sinusoidal-shaped wind stress curve
                0041 tauMax = 0.1  # wind stress maximum
                0042 x = (np.arange(nx)-1) / (nx-2)   # x-coordinate, located at XG points
                0043 y = (np.arange(ny)-.5) / (ny-2)  # y-coordinate, located at YC points
                0044 Y, X = np.meshgrid(y, x, indexing='ij')
                0045
                0046 tau = -tauMax * cos(Y*pi) # generate file for -cos(y) profile between 0-1200km
                0047 tau.astype('>f4').tofile('windx_cosy.bin')
                0048 tau = tauMax * sin(Y*pi) # generate file for +sin(y) profile between 0-1200km
                0049 tau.astype('>f4').tofile('windx_siny.bin')
                0050
                0051 # Meridional wind-stress, if desired, would be located at v-points (XC, YG)