Showing equation(3.68) for SalishSeaLake as we vary different paramters (eg h, q, L)
import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from salishsea_tools import (nc_tools, viz_tools, gsw_calls, visualisations)
import numpy.ma as ma
%matplotlib inline
mesh_mask = nc.Dataset('/home/vdo/MEOPAR/NEMO-forcing/grid/lake_meshmask201702.nc')
lakebathy = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/lake_bathy201702.nc')
tau = 0.2
rho = 1022
rhoprime = 1024
h = 30
hprime = 200 - h
epsilon = (rhoprime - rho) / rhoprime
H = h + hprime
g = 9.8
c1 = np.sqrt( g * H)
c2 = np.sqrt( epsilon * g * h * hprime / (h + hprime))
ustar = np.sqrt(tau / rho)
f = 1e-4
R1 = c1 / f
R2 = c2 / f
delta = (np.square(ustar) * t / c1 *
(np.exp(x / R1) + np.sqrt( epsilon * np.square(H) / (h * hprime)) * np.exp(x / R2) ))
deltaprime = -(np.square(ustar) * t * hprime / (c2 * (h + hprime) ) *
(np.exp(x / R2) + np.sqrt(epsilon*hprime/h)*np.exp(x/R2)))
x = np.arange(-10000,0)
delta = (np.square(ustar) / c1 *
(np.exp(x / R1) + np.sqrt( epsilon * np.square(H) / (h * hprime)) * np.exp(x / R2) )) * 86400
deltaprime = -(np.square(ustar) * hprime / (c2 * (h + hprime) ) *
(np.exp(x / R2) + np.sqrt(epsilon*hprime/h)*np.exp(x/R1))) *86400
plt.plot(x,delta,'r-', x, deltaprime, 'g-')
plt.grid('on')
R1
442718.87242357316
R2
6986.314747418699
# January with vtau = 0.2, using max and min densities
Q0 = Jan_den.max()
C = 10**(-3)
Qa = 1
g = 9.8
deltaq = Jan_den.max() - Jan_den.min()
h1 = 76.58558654785156
L = 200000
k = 1
gprime = deltaq * g / Q0
ustar = np.sqrt(np.square(14.2) * Qa * C / Q0)
W = gprime * np.square(h1) / (np.square(ustar) * L)
deltae = k*h1/W
deltae
16.850409226740055
# June with vtau = 0.2, using max and min densities
Q0 = Jun_den.max()
C = 10**(-3)
Qa = 1
g = 9.8
deltaq = Jun_den.max() - Jun_den.min()
h1 = 76.58558654785156
L = 200000
k = 1
gprime = deltaq * g / Q0
ustar = np.sqrt(np.square(14.2) * Qa * C / Q0)
W = gprime * np.square(h1) / (np.square(ustar) * L)
deltae = k*h1/W
deltae
7.437888536719055
# Jan with vtau = 0.4, using max and min densities
Q0 = Jan_den.max()
C = 10**(-3)
Qa = 1
g = 9.8
deltaq = Jan_den.max() - Jan_den.min()
h1 = 76.58558654785156
L = 200000
k = 1
gprime = deltaq * g / Q0
ustar = np.sqrt(np.square(20) * Qa * C / Q0)
W = gprime * np.square(h1) / (np.square(ustar) * L)
deltae = k*h1/W
deltae
33.426719354770988
# June with vtau = 0.4, using max and min densities
Q0 = Jun_den.max()
C = 10**(-3)
Qa = 1
g = 9.8
deltaq = Jun_den.max() - Jun_den.min()
h1 = 76.58558654785156
L = 200000
k = 1
gprime = deltaq * g / Q0
ustar = np.sqrt(np.square(20) * Qa * C / Q0)
W = gprime * np.square(h1) / (np.square(ustar) * L)
deltae = k*h1/W
deltae
14.754787813368489
Q0 = Jun_den.max()
C = 10**(-3)
Qa = 1
g = 9.8
deltaq = np.linspace(5,20,100)
h1 = 40
L = 200000
k = 1
gprime = deltaq * g / Q0
ustar = np.sqrt(np.square(20) * Qa * C / Q0)
W = gprime * np.square(h1) / (np.square(ustar) * L)
deltae = k*h1/W
plt.plot(deltaq, deltae)
plt.grid('on')
plt.title('Varying delta Q from 5 to 20, h1 = 40, L = 200km, June')
<matplotlib.text.Text at 0x7fe185a45080>
Q0 = Jan_den.max()
C = 10**(-3)
Qa = 1
g = 9.8
deltaq = np.linspace(5,20,100)
h1 = 40
L = 200000
k = 1
gprime = deltaq * g / Q0
ustar = np.sqrt(np.square(20) * Qa * C / Q0)
W = gprime * np.square(h1) / (np.square(ustar) * L)
deltae = k*h1/W
plt.plot(deltaq, deltae)
plt.grid('on')
plt.title('Varying delta Q from 5 to 20, h1 = 40, L = 200km, January')
<matplotlib.text.Text at 0x7fe185aa09e8>
Q0 = Jun_den.max()
C = 10**(-3)
Qa = 1
g = 9.8
deltaq = 12.5
h1 = np.linspace(20,60,100)
L = 200000
k = 1
gprime = deltaq * g / Q0
ustar = np.sqrt(np.square(20) * Qa * C / Q0)
W = gprime * np.square(h1) / (np.square(ustar) * L)
deltae = k*h1/W
plt.plot(h1, deltae)
plt.grid('on')
plt.title('Varying h1 from 20 to 60, delta Q = 12.5, L = 200km, June')
<matplotlib.text.Text at 0x7fe1883b1cf8>
Q0 = Jan_den.max()
C = 10**(-3)
Qa = 1
g = 9.8
deltaq = 12.5
h1 = np.linspace(20,60,100)
L = 200000
k = 1
gprime = deltaq * g / Q0
ustar = np.sqrt(np.square(20) * Qa * C / Q0)
W = gprime * np.square(h1) / (np.square(ustar) * L)
deltae = k*h1/W
plt.plot(h1, deltae)
plt.grid('on')
plt.title('Varying h1 from 20 to 60, delta Q = 12.5, L = 200km, January')
<matplotlib.text.Text at 0x7fe18839ae80>
Q0 = Jun_den.max()
C = 10**(-3)
Qa = 1
g = 9.8
deltaq = 12.5
h1 = 40
L = np.linspace(150000,250000,1000)
k = 1
gprime = deltaq * g / Q0
ustar = np.sqrt(np.square(20) * Qa * C / Q0)
W = gprime * np.square(h1) / (np.square(ustar) * L)
deltae = k*h1/W
plt.plot(L, deltae)
plt.grid('on')
plt.title('Varying L from 150K to 250K, delta Q = 12.5, h1 = 40, June')
<matplotlib.text.Text at 0x7fe1858e2d68>
Q0 = Jan_den.max()
C = 10**(-3)
Qa = 1
g = 9.8
deltaq = 12.5
h1 = 40
L = np.linspace(150000,250000,1000)
k = 1
gprime = deltaq * g / Q0
ustar = np.sqrt(np.square(20) * Qa * C / Q0)
W = gprime * np.square(h1) / (np.square(ustar) * L)
deltae = k*h1/W
plt.plot(L, deltae)
plt.grid('on')
plt.title('Varying L from 150K to 250K, delta Q = 12.5, h1 = 40, January')
<matplotlib.text.Text at 0x7fe185880c18>
June = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/initial_strat/summer2016_201702.nc')
June
<class 'netCDF4._netCDF4.Dataset'> root group (NETCDF4 data model, file format HDF5): coordinates: time_counter dimensions(sizes): deptht(40), y(898), x(398) variables(dimensions): float64 vosaline(deptht,y,x), float64 votemper(deptht,y,x), float32 deptht(deptht), float32 nav_lat(y,x), float32 nav_lon(y,x), int64 time_counter() groups:
mesh_mask = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/mesh_mask201702.nc')
bathy = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
plt.contourf(np.ma.masked_array(June.variables['votemper'][0,...],
mask = bathy.variables['Bathymetry'][:].mask))
plt.colorbar()
<matplotlib.colorbar.Colorbar at 0x7f7ecc558e80>
Jan = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/initial_strat/winter2017_201702.nc')
plt.contourf(np.ma.masked_array(Jan.variables['votemper'][0,...],
mask = bathy.variables['Bathymetry'][:].mask))
plt.colorbar()
<matplotlib.colorbar.Colorbar at 0x7f7ecc37ccc0>
def calc_rho(Sal, TempC, P):
sqrSal = np.sqrt(Sal)
R1 = ((((6.536332e-9 * TempC - 1.120083e-6) * TempC + 1.001685e-4)
* TempC - 9.095290e-3) * TempC + 6.793952e-2) * TempC - 28.263737
R2 = (((5.3875e-9 * TempC - 8.2467e-7) * TempC + 7.6438e-5)
* TempC - 4.0899e-3) * TempC + 8.24493e-1
R3 = (-1.6546e-6 * TempC + 1.0227e-4) * TempC - 5.72466e-3
SIG = (4.8314e-4 * Sal + R3 * sqrSal + R2) * Sal + R1
V350P = 1.0 / 1028.1063
SVA = -SIG * V350P / (1028.1063 + SIG)
rho = 28.106331 - SVA / (V350P * (V350P + SVA)) + 1000
return rho
Jan_den = np.ma.masked_array(calc_rho(Jan.variables['vosaline'][:],
Jan.variables['votemper'][:],
Jan.variables['deptht'][:][:, np.newaxis, np.newaxis]),
mask = 1-mesh_mask.variables['tmask'][0,...])
plt.contourf(np.arange(0,398), Jan.variables['deptht'][:],Jan_den[:,500,:])
plt.ylim(100,0)
plt.xlim(200,350)
plt.grid('on')
Jun_den = np.ma.masked_array(calc_rho(June.variables['vosaline'][:],
June.variables['votemper'][:],
June.variables['deptht'][:][:, np.newaxis, np.newaxis]),
mask = 1-mesh_mask.variables['tmask'][0,...])
plt.contourf(np.arange(0,398), June.variables['deptht'][:], Jun_den[:,500,:])
plt.colorbar()
plt.ylim(100,1)
plt.xlim(200,350)
plt.grid('on')
Jun_den.max()
1024.3840748161324
Jun_den.min()
1017.159972304221
Jan_den.max()
1024.3937942281859
Jan_den.min()
1021.2050248316444
Jun_den[:26, ...].mean()
1021.5478781651364
Jun_den[26:,...].mean()
1024.2532018354989
Jan_den[:26,...].mean()
1022.7565199808918
Jan_den[26:,...].mean()
1024.3013007326404
Equations show we are limited in the displacement range. Our tests stay in the range ish.