<jupyter_start><jupyter_text>* Calculate mouth Kelvin number.<jupyter_code>import numpy as np
import netCDF4 as nc
from salishsea_tools import (nc_tools,viz_tools,tidetools,geo_tools)
import matplotlib.pyplot as plt
import FroudeNumber as FN
import warnings
warnings.filterwarnings("ignore")
import seaborn as sns
sns.set(style="whitegrid")
sns.set_context("notebook", font_scale=1.5, rc={"lines.linewidth": 2.5})
%matplotlib inline
grid6 = nc.Dataset('/ocean/jieliu/research/meopar/river-treatment/bathy_meter_SalishSea6.nc')
X = grid6.variables['nav_lon'][:, :]
Y = grid6.variables['nav_lat'][:, :]
bathy = grid6.variables['Bathymetry'][:, :]
## tmask
mesh = nc.Dataset('/data/jieliu/MEOPAR/river-treatment/oct8_101e061e05/mesh_mask.nc')
tmask = mesh.variables['tmask'][0,:,380:510,240:397]
np_mask = np.abs(1-tmask) 
e3t = mesh.variables['e3t'][0,:,380:510,240:397]
## May results
all_T = nc.Dataset('/data/jieliu/MEOPAR/SurfaceCurrent/May2015combineall/may2015all_T.nc','r')
all_U = nc.Dataset('/data/jieliu/MEOPAR/SurfaceCurrent/May2015combineall/may2015all_U.nc','r')
## Oct results
octall_T = nc.Dataset('/data/jieliu/MEOPAR/SurfaceCurrent/Oct2014combineall/oct2014all_T.nc','r')
octall_U = nc.Dataset('/data/jieliu/MEOPAR/SurfaceCurrent/Oct2014combineall/oct2014all_T.nc','r')
## Jan results
janall_T = nc.Dataset('/data/jieliu/MEOPAR/SurfaceCurrent/Jan2016combineall/jan2016all_T.nc','r')
janall_U = nc.Dataset('/data/jieliu/MEOPAR/SurfaceCurrent/Jan2016combineall/jan2016all_U.nc','r')<jupyter_output><empty_output><jupyter_text>Functions<jupyter_code>def calculate_density(t, s):
    """Caluclates the density given temperature in deg C (t)
    and salinity in psu (s).

    returns the density as an array (rho)
    """
    rho = (
        999.842594 + 6.793952e-2 * t
        - 9.095290e-3 * t*t + 1.001685e-4 * t*t*t
        - 1.120083e-6 * t*t*t*t + 6.536332e-9 * t*t*t*t*t
        + 8.24493e-1 * s - 4.0899e-3 * t*s
        + 7.6438e-5 * t*t*s - 8.2467e-7 * t*t*t*s
        + 5.3875e-9 * t*t*t*t*s - 5.72466e-3 * s**1.5
        + 1.0227e-4 * t*s**1.5 - 1.6546e-6 * t*t*s**1.5
        + 4.8314e-4 * s*s)
    return rho
def Calculate_ReducedGravity(rho,dep,j,i):
    """
    calculate reduced gravity, g' = g(rho0 - integral(rho dz)/integral(dz))/rho0
    
    rho: density calculated based on T and S
    
    dep: level of depth selected
    
    j: along-strait index 
    
    i: across-strait index
    
    return: g_prime
    """
    g = 9.80665
    rho0 = 1023
    rho_dz = np.zeros([SSH.shape[0],dep])
    dz = np.zeros([SSH.shape[0],dep])
    for t in range(SSH.shape[0]):
        for z in range(dep):
            adj_ssh = 1+SSH[t,j,i]/np.sum(e3t[:,j,i]*tmask[:,j,i],axis = 0)
            rho_dz[t,z] = rho[t,z] * e3t[z, j,i]* adj_ssh * tmask[z,j,i]
            dz[t,z]= e3t[z, j,i]* adj_ssh * tmask[z,j,i]
    integral_rho_dz = np.nansum(rho_dz, axis = 1)
    integral_dz = np.nansum(dz,axis = 1)
    g_prime = np.sqrt(g*(np.nanmean(rho0-integral_rho_dz/integral_dz))/rho0)
    return g_prime
def Calculate_RossbyDeformationRadius(g_prime, h):
    """
    """
    f = 1e-4
    RDR = np.sqrt(g_prime*h)/f
    return RDR<jupyter_output><empty_output><jupyter_text>calculation<jupyter_code>## In May
T = all_T.variables['votemper'][:,:,40,72] ## select mouth grid point
S = all_T.variables['vosaline'][:,:,40,72]
DEPTH = all_T.variables['deptht']
SSH = all_T.variables['sossheig']
rho = calculate_density(T,S)
g_prime_may = Calculate_ReducedGravity(rho,8,40,72)
RDR_may = Calculate_RossbyDeformationRadius(g_prime_may,8)
print('Rossby deformation radius is ', RDR_may/1e3, ' km in May')
## In oct
T = octall_T.variables['votemper'][:,:,40,72]
S = octall_T.variables['vosaline'][:,:,40,72]
SSH = octall_T.variables['sossheig']
rho = calculate_density(T,S)
g_prime_oct = Calculate_ReducedGravity(rho,8,40,72)
RDR_oct = Calculate_RossbyDeformationRadius(g_prime_oct,8)
print('Rossby deformation radius is ', RDR_oct/1e3, ' km in Oct')
# In Jan
T = janall_T.variables['votemper'][:,:,40,72]
S = janall_T.variables['vosaline'][:,:,40,72]
SSH = janall_T.variables['sossheig']
rho = calculate_density(T,S)
g_prime_jan = Calculate_ReducedGravity(rho,8,40,72)
RDR_jan = Calculate_RossbyDeformationRadius(g_prime_jan,8)
print('Rossby deformation radius is ', RDR_jan/1e3, ' km in Jan')<jupyter_output>Rossby deformation radius is  13.820482326  km in Jan
{'fraction_non_alphanumeric': 0.16235401877719258, 'fraction_numerical': 0.06709411495305702, 'mean_word_length': 4.480552070263488, 'pattern_counts': {'":': 1, '<': 12, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '461988', 'license': 'Apache-2.0', 'n_tokens_mistral': 2018, 'n_tokens_neox': 1841, 'n_words': 287, 'path': 'jie/surface_current/CalculateKelvinNumber.ipynb', 'repo_name': 'SalishSeaCast/analysis-jie', 'size': '4367'}