This notebook shows some of the things that the mitgcm module is able to do.
import numpy as np
import os
import copy
import matplotlib
import matplotlib.pyplot as plt
import netCDF4
import mitgcm
%matplotlib inline
plt.rcParams['figure.figsize'] = (12, 12) # set default figure size to 12x12 inches
path_to_output = 'data/'
m = mitgcm.MITgcm_Simulation(path_to_output,'grid.all.nc')
Print out a list of the variables contained in the netCDF file.
mitgcm.functions.show_variables('state.all.nc')
[u'Xp1', u'Y', u'Z', u'X', u'Yp1', u'Zl', u'T', u'iter', u'U', u'V', u'Temp', u'S', u'Eta', u'W']
f = netCDF4.Dataset('state.all.nc')
total_time_levels = len(f.variables['T'][:])
m['time'] = copy.deepcopy(f.variables['T'][:])
f.close()
print 'time levels in data file = ', total_time_levels
time levels in data file = 2
m.zonal_velocity = mitgcm.Upoint_field(m,'state.all.nc','U',time_level=0)
m.meridional_velocity = mitgcm.Vpoint_field(m,'state.all.nc','V',time_level=0)
m.vertical_velocity = mitgcm.Wpoint_field(m,'state.all.nc','W',time_level=0)
m.temperature = mitgcm.Tracerpoint_field(m,'state.all.nc','Temp',time_level=0)
m.free_surface = mitgcm.Tracerpoint_field(m,'state.all.nc','Eta',time_level=0)
m.salinity = mitgcm.Tracerpoint_field(m,'state.all.nc','S',time_level=0)
m.density = mitgcm.Density(m,Talpha=0.0002, Sbeta=0, RhoNil=1025, cp=4000,temp_field='Temp', density_field='RHO')
m.pressure = mitgcm.Pressure(m,density_field='RHO',ETAN_field='Eta')
m.bernoulli = mitgcm.Bernoulli(m,density_field='RHO', UVEL_field='U', VVEL_field='V')
im = plt.pcolormesh(m.grid['X'][:],m.grid['Y'][:],m.density['RHO'][0,:-1,:-1]-1000,cmap='rainbow_r')
plt.colorbar(im)
<matplotlib.colorbar.Colorbar instance at 0x1117ad0e0>
The simulated glider path is shown as a black line on the figure below. It's just east of the western boundary current, about halfway up.
x = np.linspace(3e5,5e5,160)
y = np.linspace(1e6,1e6,160)
z = np.linspace(-1,-500,20)
z = np.append(z,z[::-1],)
z = np.tile(z,4)
temperature_cast = mitgcm.streamlines.extract_along_path3D(x,y,z,m.grid['X'][:],m.grid['Y'][:],m.grid['Z'][:],
m.temperature['Temp'])
im = plt.pcolormesh(m.grid['X'][:],m.grid['Y'][:],m.density['RHO'][0,:-1,:-1]-1000,cmap='rainbow_r')
plt.colorbar(im)
plt.plot(x,y,'k')
[<matplotlib.lines.Line2D at 0x112409210>]
import scipy.interpolate
X,Z = np.meshgrid(x/1e3,z)
TEMPERATURE = scipy.interpolate.griddata((x/1e3,z),temperature_cast,(X,Z))
plt.figure()
CS = plt.contour(X,Z,TEMPERATURE,np.arange(0,30,2),cmap='gnuplot',vmin=10,vmax=30)
#plt.clabel(CS,fontsize=15, use_clabeltext=True)
im = plt.scatter(x/1e3,z,c = temperature_cast,linewidths=0,s=50,cmap='gnuplot',vmin=10,vmax=30)
plt.colorbar(im)
plt.ylabel('Depth (m)',fontsize=20)
plt.xlabel('x coordinate (km)',fontsize=20)
plt.title('Temperature data from a simulated glider track (contours plotted every 2 degrees C)')
<matplotlib.text.Text at 0x112299150>
The surface is masked if the specified value isn't found at that (x,y) location.
depth_1027 = mitgcm.functions.calc_surface(m.density['RHO'][:],1027,m.grid['Z'][:])
/Users/doddridge/anaconda/lib/python2.7/site-packages/mitgcm-0.1-py2.7.egg/mitgcm/functions.py:30: RuntimeWarning: input field is not strictly monotonic in search direction. Strange results may occur. warnings.warn("input field is not strictly monotonic in search direction. Strange results may occur.", RuntimeWarning) /Users/doddridge/anaconda/lib/python2.7/site-packages/mitgcm-0.1-py2.7.egg/mitgcm/functions.py:79: RuntimeWarning: invalid value encountered in divide input_array[indices_min[:,:] + direction,indsy,indsx]))*
u_1027 = mitgcm.functions.extract_on_surface(m.zonal_velocity['U'][:,:,:-1],depth_1027,m.grid['Z'][:])
v_1027 = mitgcm.functions.extract_on_surface(m.meridional_velocity['V'][:,:-1,:],depth_1027,m.grid['Z'][:])
BP_1027 = mitgcm.functions.extract_on_surface(m.bernoulli['BP'][:],depth_1027,m.grid['Z'][:])
fig = plt.figure()
im = plt.pcolormesh(m.grid['X'][:],m.grid['Y'][:],depth_1027,vmin=-500,vmax=-0,cmap='YlGnBu_r')
plt.contour(m.grid['X'][:],m.grid['Y'][:],BP_1027,20)
plt.colorbar(im)
<matplotlib.colorbar.Colorbar instance at 0x1140cc878>
Let's calculate the northward volum flux in the lyaer between $\sigma$=26 and $\sigma$=27.
The function returns zero if either of the specified surfaces are masked at the same location.
depth_1026 = mitgcm.functions.calc_surface(m.density['RHO'][:],1026,m.grid['Z'][:])
v_volume_transport = mitgcm.functions.layer_integrate(depth_1027,depth_1026,m.grid['Z'][:],
integrand=m.meridional_velocity['V'][:,:-1,:])
v_volume_transport = v_volume_transport*m.grid['dxF'][:]/1e6
plt.pcolormesh(m.grid['X'][:],m.grid['Y'][:],v_volume_transport,cmap='seismic',vmin=-0.4,vmax=0.4)
plt.colorbar()
plt.title('Meridional volume transport in Sv per grid cell')
<matplotlib.text.Text at 0x114158a10>
# Get Bernoulli value at
x0=6e5
y0=5e5
BP_value = mitgcm.streamlines.bilinear_interp(x0,y0,BP_1027,m.grid['X'][:],
m.grid['Y'][:],len(m.grid['X'][:]),len(m.grid['Y'][:]))
BP_contour = plt.contour(m.grid['X'][:],m.grid['Y'][:],BP_1027,[BP_value])
BP_contour_levels = BP_contour.levels
BP_contour_vertices = BP_contour.collections[0].get_paths()[0].vertices
x_BP_contour = BP_contour_vertices[:,0]
y_BP_contour = BP_contour_vertices[:,1]
x_stream,y_stream,t_stream = mitgcm.streamlines.stream2(u_1027,v_1027,
x0,y0,m.grid,
t_int=8e7,delta_t=3600)
fig = plt.figure()
im = plt.pcolormesh(m.grid['X'][:],m.grid['Y'][:],depth_1027,cmap='YlGnBu_r')
plt.plot(x_stream,y_stream,'k',linewidth=2)
plt.plot(x0,y0,'r*',markersize=10)
plt.contour(m.grid['X'][:],m.grid['Y'][:],BP_1027,[BP_value],colors='y',linewidth=3)
plt.colorbar(im)
<matplotlib.colorbar.Colorbar instance at 0x114414368>
These two lines don't match, but that's expected, the bernoulli contour gives the geostrophic streamline, not the actual streamline.
#3D streamline from same starting point
z0 = mitgcm.streamlines.bilinear_interp(x0,y0,depth_1027,
m.grid['X'][:],m.grid['Y'][:],len(m.grid['X'][:]),len(m.grid['Y'][:]))
stream3d_x,stream3d_y,stream3d_z,stream3d_t = mitgcm.streamlines.stream3(m.zonal_velocity['U'][:],
m.meridional_velocity['V'][:],m.vertical_velocity['W'][:],
x0,y0,z0,m.grid,
t_int=8e7,delta_t=3600)
fig = plt.figure()
im = plt.pcolormesh(m.grid['X'][:],m.grid['Y'][:],depth_1027,cmap='YlGnBu_r')
plt.plot(x_stream,y_stream,'k',linewidth=2)
plt.plot(stream3d_x,stream3d_y,'r',linewidth=2)
plt.contour(m.grid['X'][:],m.grid['Y'][:],BP_1027,[BP_value],colors='y')
plt.colorbar(im)
<matplotlib.colorbar.Colorbar instance at 0x1161828c0>
The 3D (red) and 2D (black) streamlines don't match up perfectly because the 3D streamline was not constrained to stay on the density surface.
We can use some of the inbuilt functions to check if the density of the fluid parcel stays constant along the 3D streamline.
stream_density = mitgcm.streamlines.extract_along_path3D(stream3d_x[:],stream3d_y[:],stream3d_z[:],
m.grid['X'][:],m.grid['Y'][:],m.grid['Z'][:],
m.density['RHO'])
#other variables can be extracted too - remember to use the correct x and y vectors
cori = mitgcm.streamlines.extract_along_path2D(stream3d_x[:],stream3d_y[:],
m.grid['Xp1'][:],m.grid['Yp1'][:],m.grid['fCoriG'])
im = plt.scatter(stream3d_x[:],stream3d_y[:],c=(stream_density-1000),lw=0)
plt.colorbar(im)
plt.xlim((0,2e6))
plt.ylim((0,2e6))
(0, 2000000.0)
The density is changing along the streamline, showing that the velocity field contains a diapycnal component.
The warnings are ok to ignore - the values get cleaned up by an np.nan_to_num() call before being returned and should be supressed in the future.
m.density.take_d_dx(m,input_field='RHO',output_field='dRHO_dx')
m.density.take_d_dy(m,input_field='RHO',output_field='dRHO_dy')
m.density.take_d_dz(m,input_field='RHO',output_field='dRHO_dz')
m.density.take_d_dx(m,input_field='dRHO_dx',output_field='d2RHO_dx2')
m.density.take_d_dy(m,input_field='dRHO_dy',output_field='d2RHO_dy2')
m.density.take_d_dz(m,input_field='dRHO_dz',output_field='d2RHO_dz2')
m.temperature.take_d_dx(m,input_field='Temp',output_field='dTemp_dx')
m.temperature.take_d_dy(m,input_field='Temp',output_field='dTemp_dy')
m.temperature.take_d_dz(m,input_field='Temp',output_field='dTemp_dz')
m.temperature.take_d_dx(m,input_field='dTemp_dx',output_field='d2Temp_dx2')
m.temperature.take_d_dx(m,input_field='dTemp_dy',output_field='d2Temp_dy2')
m.temperature.take_d_dx(m,input_field='dTemp_dz',output_field='d2Temp_dz2')
m.bernoulli.take_d_dx(m,input_field='BP',output_field='dBP_dx')
m.bernoulli.take_d_dy(m,input_field='BP',output_field='dBP_dy')
/Users/doddridge/anaconda/lib/python2.7/site-packages/mitgcm-0.1-py2.7.egg/mitgcm/core.py:637: RuntimeWarning: invalid value encountered in divide wet_mask_TH[:,:,i+1]*dxC[:,i+1])) /Users/doddridge/anaconda/lib/python2.7/site-packages/mitgcm-0.1-py2.7.egg/mitgcm/core.py:674: RuntimeWarning: invalid value encountered in divide wet_mask_TH[:,j+1,:]*dyC[j+1,:]))
$\frac{\partial \rho}{\partial t} + \mathbf{u} \cdot \nabla \rho = A \nabla^{2} \rho$
We need to extract the various terms along the streamline.
u = mitgcm.streamlines.extract_along_path3D(stream3d_x[:],stream3d_y[:],stream3d_z[:],
m.grid['Xp1'][:],m.grid['Y'][:],m.grid['Z'][:],
m.zonal_velocity['U'][:])
v = mitgcm.streamlines.extract_along_path3D(stream3d_x[:],stream3d_y[:],stream3d_z[:],
m.grid['X'][:],m.grid['Yp1'][:],m.grid['Z'][:],
m.meridional_velocity['V'][:])
w = mitgcm.streamlines.extract_along_path3D(stream3d_x[:],stream3d_y[:],stream3d_z[:],
m.grid['X'][:],m.grid['Y'][:],m.grid['Zl'][:],
m.vertical_velocity['W'][:])
d_rho_dx = mitgcm.streamlines.extract_along_path3D(stream3d_x[:],stream3d_y[:],stream3d_z[:],
m.grid['X'][:],m.grid['Y'][:],m.grid['Z'][:],
m.density['dRHO_dx'][:])
d_rho_dy = mitgcm.streamlines.extract_along_path3D(stream3d_x[:],stream3d_y[:],stream3d_z[:],
m.grid['X'][:],m.grid['Y'][:],m.grid['Z'][:],
m.density['dRHO_dy'][:])
d_rho_dz = mitgcm.streamlines.extract_along_path3D(stream3d_x[:],stream3d_y[:],stream3d_z[:],
m.grid['X'][:],m.grid['Y'][:],m.grid['Z'][:],
m.density['dRHO_dz'][:])
d2_rho_dx2 = mitgcm.streamlines.extract_along_path3D(stream3d_x[:],stream3d_y[:],stream3d_z[:],
m.grid['X'][:],m.grid['Y'][:],m.grid['Z'][:],
m.density['d2RHO_dx2'][:])
d2_rho_dy2 = mitgcm.streamlines.extract_along_path3D(stream3d_x[:],stream3d_y[:],stream3d_z[:],
m.grid['X'][:],m.grid['Y'][:],m.grid['Z'][:],
m.density['d2RHO_dy2'][:])
d2_rho_dz2 = mitgcm.streamlines.extract_along_path3D(stream3d_x[:],stream3d_y[:],stream3d_z[:],
m.grid['X'][:],m.grid['Y'][:],m.grid['Z'][:],
m.density['d2RHO_dz2'][:])
cori = mitgcm.streamlines.extract_along_path2D(stream3d_x[:],stream3d_y[:],
m.grid['Xp1'][:],m.grid['Yp1'][:],m.grid['fCoriG'][:])
d_BP_dx = mitgcm.streamlines.extract_along_path3D(stream3d_x[:],stream3d_y[:],stream3d_z[:],
m.grid['X'][:],m.grid['Y'][:],m.grid['Z'][:],
m.bernoulli['dBP_dx'][:])
d_BP_dy = mitgcm.streamlines.extract_along_path3D(stream3d_x[:],stream3d_y[:],stream3d_z[:],
m.grid['X'][:],m.grid['Y'][:],m.grid['Z'][:],
m.bernoulli['dBP_dy'][:])
# calculate the geostrophic velocity along the streamline.
vel_on_contour = np.sqrt(d_BP_dx**2 + d_BP_dy**2)/(cori)
dissipation = 1e-5*d2_rho_dz2 + 100*(d2_rho_dx2 + d2_rho_dy2)
advective_term = u*d_rho_dx + v*d_rho_dy + w*d_rho_dz
This shows that advection is generally decreasing desnsity in the western boundary current
im = plt.scatter(stream3d_x[:],stream3d_y[:],c=(-advective_term),lw=0,vmin=-6e-8,vmax=6e-8,cmap='RdBu_r')
plt.plot(stream3d_x[:],stream3d_y[:],'k')
plt.colorbar(im)
plt.xlim((0,2e6))
plt.ylim((0,2e6))
(0, 2000000.0)
This term is generally increasing density in the western boundary current.
im = plt.scatter(stream3d_x[:],stream3d_y[:],c=(dissipation),lw=0,vmin=-3e-8,vmax=3e-8,cmap='RdBu_r')
plt.plot(stream3d_x[:],stream3d_y[:],'k')
plt.colorbar(im)
plt.xlim((0,2e6))
plt.ylim((0,2e6))
(0, 2000000.0)