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') mitgcm.functions.show_variables('state.all.nc') 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 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) 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') 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)') depth_1027 = mitgcm.functions.calc_surface(m.density['RHO'][:],1027,m.grid['Z'][:]) 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) 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') # 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) #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) 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)) 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') 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 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)) 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))