In [ ]:
# Uses python3
# Notebook with images at https://gist.github.com/adcroft/ebdcf46b3df13c42fa4aeda399dd5a1a
import numpy
import scipy.io.netcdf
import matplotlib.pyplot as plt
%matplotlib inline
In [ ]:
icz=scipy.io.netcdf.netcdf_file('z/Initial_state.nc', mmap=False)
ziz=-icz.variables['Interface'][:]
zcz=-icz.variables['Layer'][:]
icr=scipy.io.netcdf.netcdf_file('rho/Initial_state.nc', mmap=False)
xq=numpy.concatenate(([0],icz.variables['lonq'][:]))
xc=icz.variables['lonh'][:]
In [ ]:
# Initial conditions
plt.figure(figsize=(14,5));
plt.subplot(121)
plt.pcolormesh(xq,ziz,icz.variables['Temp'][0,:,0,:]); plt.colorbar();
plt.contour(xc,zcz,icz.variables['Salt'][0,:,0,:],numpy.arange(34.5,35.5,.05),colors='k');
plt.xlabel('x [km]'); plt.ylabel('z* [m]'); plt.title('z* coordinate I.C.s');
plt.subplot(122)
z = numpy.zeros((ziz.shape[0],xq.shape[0])); z[:,0] = icr.variables['eta'][-1,:,0,0]; z[:,-1] = icr.variables['eta'][-1,:,0,-1]
z[:,1:-1] = 0.5 * ( icr.variables['eta'][-1,:,0,:-1] + icr.variables['eta'][-1,:,0,1:] ) # interfaces at u-points
plt.pcolormesh(xq,z,icr.variables['Temp'][0,:,0,:]); plt.colorbar();
z = 0.5 * ( icr.variables['eta'][0,:-1,0,:] + icr.variables['eta'][0,1:,0,:] )
plt.contour(xc+0*z,z,icr.variables['Salt'][0,:,0,:],numpy.arange(34.5,35.5,.05),colors='k');
plt.xlabel('x [km]'); plt.ylabel('z* [m]'); plt.title(r'$\rho$ coordinate I.C.s');
plt.figure(figsize=(14,5));
plt.subplot(121)
plt.plot(icz.variables['Salt'][0,:,0,:].flatten()+1e3,icz.variables['Temp'][0,:,0,:].flatten(),'.');
plt.xlabel('$\sigma$ [kg/m$^3$]'); plt.ylabel(r'$\theta$ [conc.]'); plt.title('z* coordinate I.C.s'); plt.xlim(1034.5,1035.5);
plt.subplot(122)
plt.plot(icr.variables['Salt'][0,:,0,:].flatten()+1e3,icr.variables['Temp'][0,:,0,:].flatten(),'.'); plt.xlim(1034.5,1035.5);
plt.xlabel('$\sigma$ [kg/m$^3$]'); plt.ylabel(r'$\theta$ [conc.]'); plt.title(r'$\rho$ coordinate I.C.s');
In [ ]:
znc=scipy.io.netcdf.netcdf_file('z/prog.nc', mmap=False)
rnc=scipy.io.netcdf.netcdf_file('rho/prog.nc', mmap=False)
In [ ]:
rnc.variables['temp'].shape
In [ ]:
# Evolution/final state
plt.figure(figsize=(14,5));
plt.subplot(121)
plt.pcolormesh(xq,ziz,znc.variables['temp'][-1,:,0,:]); plt.colorbar();
plt.contour(xc,zcz,znc.variables['salt'][-1,:,0,:],numpy.arange(34.5,35.5,.05),colors='k');
plt.xlabel('x [km]'); plt.ylabel('z* [m]'); plt.title('z* coordinate');
plt.subplot(122)
z = numpy.zeros((ziz.shape[0],xq.shape[0])); z[:,0] = rnc.variables['e'][-1,:,0,0]; z[:,-1] = rnc.variables['e'][-1,:,0,-1]
z[:,1:-1] = 0.5 * ( rnc.variables['e'][-1,:,0,:-1] + rnc.variables['e'][-1,:,0,1:] ) # interfaces at u-points
plt.pcolormesh(xq,z,rnc.variables['temp'][-1,:,0,:], cmap='jet_r'); plt.colorbar();
z = 0.5 * ( rnc.variables['e'][-1,:-1,0,:] + rnc.variables['e'][-1,1:,0,:] ) # middle of layer at h-points
plt.contour(xc+0*z,z,rnc.variables['salt'][-1,:,0,:]+.001,numpy.arange(34.5,35.5,.05),colors='k');
plt.xlabel('x [km]'); plt.ylabel('z* [m]'); plt.title('z* coordinate');
plt.figure(figsize=(14,5));
plt.subplot(121)
plt.plot(icz.variables['Salt'][0,:,0,:].flatten()+1e3,icz.variables['Temp'][0,:,0,:].flatten(),'c.',label='t=0 day');
n=1; plt.plot(znc.variables['salt'][n,:,0,:].flatten()+1e3,znc.variables['temp'][n,:,0,:].flatten(),'r.',hold=True,
label='t=%i day'%(znc.variables['Time'][n]));
n=-1; plt.plot(znc.variables['salt'][n,:,0,:].flatten()+1e3,znc.variables['temp'][n,:,0,:].flatten(),'k.',hold=True,
label='t=%i day'%(znc.variables['Time'][-1]));
plt.xlabel('$\sigma$ [kg/m$^3$]'); plt.ylabel(r'$\theta$ [conc.]'); plt.xlim(1034.5,1035.5);
plt.legend(loc='upper left'); plt.title('z* coordinate');
plt.subplot(122)
plt.plot(icr.variables['Salt'][0,:,0,:].flatten()+1e3,icr.variables['Temp'][0,:,0,:].flatten(),'c.',label='t=0 day');
n=1; plt.plot(rnc.variables['salt'][n,:,0,:].flatten()+1e3,rnc.variables['temp'][n,:,0,:].flatten(),'r.',hold=True,
label='t=%i day'%(rnc.variables['Time'][n]));
n=-1; plt.plot(rnc.variables['salt'][n,:,0,:].flatten()+1e3,rnc.variables['temp'][n,:,0,:].flatten(),'k.',hold=True,
label='t=%i day'%(rnc.variables['Time'][-1]));
plt.xlabel('$\sigma$ [kg/m$^3$]'); plt.ylabel(r'$\theta$ [conc.]'); plt.xlim(1034.5,1035.5);
plt.legend(loc='upper left'); plt.title(r'$\rho$ coordinate');