In [ ]:
import math
import numpy
import matplotlib.pyplot as plt
import scipy.io
%matplotlib inline
In [ ]:
zNODC = numpy.concatenate((numpy.arange(0,100,5),numpy.arange(100,500,25),numpy.arange(500,2000,50),numpy.arange(2000,6600,100)))
dzNODC = zNODC[1:] - zNODC[:-1]
zCM21 = numpy.array([0, 10.1, 20.2, 30.3, 40.5, 50.7, 61, 71.5, 82.2, 93.1, 104.4, 116.1,
128.5, 141.8, 156.4, 172.7, 191.3, 213, 238.9, 270.3, 308.8, 356.2,
414.6, 485.9, 571.9, 673.8, 791.9, 925.9, 1074.6, 1236.2, 1408.8, 1590.3,
1778.9, 1973, 2171.3, 2372.7, 2576.4, 2781.8, 2988.5, 3196.1, 3404.4,
3613.1, 3822.2, 4031.5, 4241, 4450.7, 4660.5, 4870.3, 5080.2, 5290.1,
5499.1])
dzCM21 = zCM21[1:] - zCM21[:-1]
dzOM4 = numpy.array([2, 2, 2, 2, 2.01, 2.01, 2.02, 2.03, 2.05, 2.08, 2.11, 2.15, 2.2, 2.27,
2.34, 2.44, 2.55, 2.69, 2.85, 3.04, 3.27, 3.54, 3.85, 4.22, 4.66, 5.18,
5.79, 6.52, 7.37, 8.37, 9.55, 10.94, 12.57, 14.48, 16.72, 19.33, 22.36,
25.87, 29.91, 34.53, 39.79, 45.72, 52.37, 59.76, 67.89, 76.74, 86.29,
96.47, 107.2, 118.35, 129.81, 141.42, 153.01, 164.41, 175.47, 186.01,
195.9, 205.01, 213.27, 220.6, 226.99, 232.43, 236.96, 240.63, 243.52,
245.72, 247.33, 248.45, 249.18, 249.62, 249.86, 249.96, 249.99, 250, 250])
dzDunne = numpy.array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2.1, 2.23, 2.4, 2.61, 2.85, 3.13, 3.46, 3.84, 4.28, 4.8, 5.39, 6.07,
6.86, 7.77, 8.83, 10.04, 11.44, 13.06, 14.92, 17.06, 19.53, 22.37, 25.62,
29.36, 33.64, 38.52, 44.09, 50.43, 57.62, 65.76, 74.93, 85.24, 96.78,
109.64, 123.9, 139.62, 156.85, 175.59, 195.82, 217.45, 240.32, 264.18,
288.63, 313.13, 336.87, 358.61, 376.17, 376.17, 1000, 1000])
def dz2(N,Htotal,Hmin,p,f2=0.,prec=0.01):
dz = Hmin * numpy.ones( N )
k = numpy.linspace(0, N, N).astype(numpy.float64)/N
dzp = k**p * ( 1 - f2*k )
dzp = (Htotal - dz.sum()) * ( dzp / dzp.sum() )
dzp = prec * numpy.round( dzp / prec)
dzp = (Htotal - dz.sum()) * ( dzp / dzp.sum() )
dzp = prec * numpy.round( dzp / prec)
dzp[-1] += Htotal - (dz+dzp).sum()
dzp = prec * numpy.round( dzp / prec)
return dz + dzp
dzOM4new=dz2(75,6500,2,6)
print(dzOM4new)
print(dzOM4new.sum())
In [ ]:
def zFromDz(dz):
"""
Sums dz to return zInterface and zCenter.
"""
zInt = numpy.zeros(dz.shape[0]+1)
zInt[1:] = -numpy.cumsum(dz)
zCenter = 0.5*( zInt[:-1] + zInt[1:] )
return zInt, zCenter
def plot1Dz(dz, scale=1, fmt='.', hold=True, label=None):
"""Plots dz as bars in z-space"""
zI, zC = zFromDz(dz)
plt.errorbar(dz/scale,zC,yerr=dz/2,fmt=fmt,hold=hold,label=label)
plt.gca().invert_yaxis()
def plotRatio(dz, scale=1, fmt='.', hold=True, label=None):
"""Plots dz as bars in z-space"""
zI, zC = zFromDz(dz)
plt.plot(dz[1:]/dz[:-1],zI[1:-1],hold=hold,label=label)
plt.gca().invert_yaxis()
def plotDzs():
plot1Dz(dzCM21, label='CM2.1', hold=False)
#plot1Dz(dzNODC, label='NODC')
plot1Dz(dzOM4, label='OM4')
plot1Dz(dzDunne, label='OM4 JPD')
plot1Dz(dzOM4new, label='OM4 new')
plt.xlabel('$\Delta$z [m]'); plt.ylabel('z [m]')
plt.legend(loc='upper right')
def plotRatios():
plotRatio(dzCM21, label='CM2.1', hold=False)
#plotRatio(dzNODC, label='NODC')
plotRatio(dzOM4, label='OM4')
plotRatio(dzDunne, label='OM4 JPD')
plotRatio(dzOM4new, label='OM4 new')
plt.xlabel('Ratio $\Delta$z$_{k+1}$/$\Delta$z$_k$'); plt.ylabel('z [m]')
plt.legend(loc='upper right')
plt.figure(figsize=(16,7))
plt.subplot(133);plotDzs();plt.ylim(-6500,0); plt.xlim(0,1010);
plt.subplot(132);plotDzs();plt.ylim(-1000,0); plt.xlim(0,150);
plt.subplot(131);plotDzs();plt.ylim(-300,0); plt.xlim(0,40);
plt.figure(figsize=(16,7))
plt.subplot(133);plotRatios();plt.ylim(-6500,0); plt.xlim(.97,1.4);
plt.subplot(132);plotRatios();plt.ylim(-1000,0); plt.xlim(.97,1.4);
plt.subplot(131);plotRatios();plt.ylim(-300,0); plt.xlim(.97,1.4);
In [ ]:
fname = 'vgrid_75_2m_575m.nc'
nc = scipy.io.netcdf_file(fname,'w')
nc.createDimension('nz',dzOM4new.shape[0])
nc.filename = fname
v = nc.createVariable('dz','double',('nz',))
v.units = 'm'
v.long_name = 'z* coordinate level thickness'
v[:] = dzOM4new[:]
nc.close()
In [ ]:
zt = numpy.array([2.5, 10, 20, 30, 50, 75, 100, 125, 150, 200, 250, 300, 400, 500, 600,
700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1750, 2000, 2500,
3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500 ])
zw = numpy.array([0, 5, 15, 25, 40, 62.5, 87.5, 112.5, 137.5, 175, 225, 275, 350, 450,
550, 650, 750, 850, 950, 1050, 1150, 1250, 1350, 1450, 1625, 1875, 2250,
2750, 3250, 3750, 4250, 4750, 5250, 5750, 6250, 6750 ])
In [ ]:
plt.figure(); plt.errorbar(zw[1:]-zw[:-1],-zt,yerr=0.5*(zw[1:]-zw[:-1]),fmt='.')
plt.figure(); plt.plot( 0.5*(zw[1:]+zw[:-1]) - zt, -zt); plt.ylim(-3000,0)
In [ ]:
fname = 'analysis_vgrid_lev35.v1.nc'
nc = scipy.io.netcdf_file(fname,'w')
nc.createDimension('zt',zt.shape[0])
nc.createDimension('zw',zw.shape[0])
nc.filename = fname
v1 = nc.createVariable('zt','double',('zt',))
v2 = nc.createVariable('zw','double',('zw',))
v1.units = 'm'
v2.units = 'm'
v1.long_name = 'Diagnostic depth coordinate level position'
v1.comment = 'Used for diagnostics only, based on WOA09 standard levels'
v2.long_name = 'Diagnostic depth coordinate interface position'
v2.comment = 'Used for diagnostics only, based on WOA09 standard levels'
v1[:] = zt[:]
v2[:] = zw[:]
nc.close()