%matplotlib inline
from matplotlib import pylab
import matplotlib.pyplot as plt
import netCDF4 as NC
import numpy as np

fT = NC.Dataset('../../Results/Oct18OneDay/SalishSea_1h_20021018_20021018_grid_T.nc','r')
depth = fT.variables['deptht'][:]

fW = NC.Dataset('../../Results/Oct18OneDay/SalishSea_1h_20021018_20021018_grid_W.nc','r')
wvel = fW.variables['vovecrtz'][:]
mw = wvel == 0
wvel_masked = np.ma.array(wvel,mask=mw)

level = 28
print 0.5*(depth[28]+depth[29])
plt.figure(figsize=(10,10))
plt.subplot(2,2,1)
plt.pcolormesh(wvel_masked[15,level,300:370,220:320])
plt.colorbar()
plt.subplot(2,2,2)
plt.pcolormesh(wvel_masked[2,level,300:370,220:320])
plt.colorbar()
plt.subplot(2,2,3)
plt.pcolormesh(wvel_masked[4,level,300:370,220:320])
plt.colorbar()
plt.subplot(2,2,4)
plt.pcolormesh(wvel_masked[21,level,300:370,220:320])
plt.colorbar()

alpha = np.zeros(39)
for level in np.arange(1,39):
    print level,depth[level],np.max(wvel_masked[:,level]), np.min(wvel_masked[:,level]),(depth[level+1]-depth[level-1])/2./np.max(wvel_masked[:,level]),(depth[level+1]-depth[level-1])/2./np.min(wvel_masked[:,level])
    alpha[level] = np.maximum(-(depth[level]-depth[level-1])/np.min(wvel_masked[:,level]),(depth[level]-depth[level-1])/np.max(wvel_masked[:,level]))

plt.plot(np.log(alpha[1:39]))
plt.title("Log dz/w")
plt.xlabel("Vertical Level")
print np.min(alpha[1:39])