This notebook generates new boundary conditions and initial conditions for the triple tanh vertical resolution.
import numpy as np
import matplotlib.pyplot as plt
import netCDF4 as nc
from salishsea_tools import nc_tools
%matplotlib inline
The triple tanh profile has been hard coded into domzgr.F90. The profile is
$$ z(k) = ppsur + ppa0k + ppa1 ppacr \log(\cosh(\frac{k-ppkth}{ppacr})) + ppa2 ppacr2 \log\cosh(\frac{k-ppkth2} {ppacr2})) + ppa3 ppacr3 \log\cosh(\frac{k-ppkth3} {ppacr3}))$$$$ e3(k) = ppa0 + ppa1 \tanh(\frac{k-ppkth}{ppacr}) + ppa2\tanh(\frac{k-ppkth2}{ppacr2}) + ppa3\tanh(\frac{k-ppkth3}{ppacr3}) $$where $k$ is the grid level, $z(k)$ is the depth of the grid level, and $e3(k)$ is the grid spacing at that grid level (NEMO manual).
def NEMO_tanh(ks, a, kth, acr):
"""A generic tanh used in NEMO"""
return a* np.tanh((ks-kth)/acr)
def NEMO_tanh_integral(ks, a, kth, acr):
"""A generic integrated tanh used in NEMO grid"""
return a*acr*np.log(np.cosh((ks-kth)/acr))
Quick visualization
ks = np.arange(1,41)
#translated from above
ppa1 = 3/2.; ppkth = 16.7; ppacr = 6.;
ppa2 = 13/2.; ppkth2 = 22.; ppacr2 = 5.5;
ppa3 = 16/2.; ppkth3 = 33.877; ppacr3 = 5.;
#setting e3(1) = 1
ppa0 = 1 - (NEMO_tanh(ks[0], ppa1, ppkth, ppacr) + NEMO_tanh(ks[0], ppa2, ppkth2, ppacr2) +
NEMO_tanh(ks[0], ppa3, ppkth3, ppacr3) )
#setting z(1) = 0
ppsur = -(NEMO_tanh_integral(ks[0], ppa1, ppkth, ppacr) + NEMO_tanh_integral(ks[0], ppa2, ppkth2, ppacr2) +
(NEMO_tanh_integral(ks[0], ppa3, ppkth3, ppacr3) + ppa0*ks[0]))
print 'ppa0', ppa0
print 'ppsur', ppsur
dz = ppa0 + NEMO_tanh(ks, ppa1, ppkth, ppacr) + NEMO_tanh(ks, ppa2, ppkth2, ppacr2) + NEMO_tanh(ks, ppa3, ppkth3, ppacr3)
z = ppsur + ppa0*ks + (NEMO_tanh_integral(ks, ppa1, ppkth, ppacr) + NEMO_tanh_integral(ks, ppa2, ppkth2, ppacr2) +
NEMO_tanh_integral(ks, ppa3, ppkth3, ppacr3))
ppa0 16.9777762082 ppsur -381.364772219
#plot
fig,axs=plt.subplots(1,3,figsize=(15,5))
for ax in axs:
ax.plot(dz,z,'-o')
ax.grid()
ax.set_xlabel('Grid spacing (m)')
ax.set_ylabel('Depth (m)')
axs[0].set_ylim([450,0])
axs[1].set_ylim([80,0])
axs[2].set_ylim([30,0])
print 'ppsur', ppsur
print 'ppa0', ppa0
print 'Surface: ',z[0]
print 'First grid cell thickness', dz[0]
print 'Bottom:', z[-1]
ppsur -381.364772219 ppa0 16.9777762082 Surface: 0.0 First grid cell thickness 1.0 Bottom: 443.304590924
Now that that you have decided on a different vertical grid structure, add the above parameters to your namelist.domain and set up a "fake" simulation. The purpose of the fake simulation is to generate a new mesh_mask.nc file which contains the grids for your new domain. Once you have the new grid, you can interpolate the boundary and initial conditions.
Run NEMO with nn_msh = 1 in namelist.domain. You only need to run a short simulation (one time step). This will generate a file called mesh_mask.nc
This file was copied to /data/nsoontie/MEOPAR/2Ddomain/vertical_resolution/trpltanh
The new mesh gives us the z-levels needed for interpolation. Choose your initial condition file to interpolate.
Open old dataset
#open old file
fold = nc.Dataset('/data/nsoontie/MEOPAR/2Ddomain/initial_conditions/TS_uniform_36.nc')
Sold = fold.variables['vosaline'][:]
Told = fold.variables['votemper'][:]
#masking
Sold = np.ma.masked_values(Sold,0)
Told = np.ma.masked_values(Told,0)
#old grid
f=nc.Dataset('/data/nsoontie/MEOPAR/2Ddomain/grid/mesh_mask.nc')
gdept_old = f.variables['gdept_0'][0,:,:,:]
gdept_1d_old = f.variables['gdept_1d'][0,:]
xold = f.variables['nav_lon'][:]
yold = f.variables['nav_lat'][:]
print gdept_1d_old[:]
print gdept_1d_old.shape
[ 0.50000027 1.50000314 2.5000115 3.50003055 4.50007042 5.50015083 6.50031022 7.50062342 8.50123623 9.50243254 10.5047653 11.50931127 12.51816684 13.53541212 14.56898216 15.63428737 16.76117342 18.00713456 19.48178514 21.38997868 24.10025665 28.22991514 34.68575798 44.51772486 58.48433368 76.58558445 98.06295924 121.8665184 147.08945807 173.11448217 199.57304923 226.26030574 253.06663733 279.93454976 306.83419736 333.75016973 360.6745318 387.60320347 414.53408835 441.46610968] (40,)
Open new mesh file
mesh = nc.Dataset('/data/nsoontie/MEOPAR/2Ddomain/vertical_resolution/trpl_tanh/mesh_mask.nc')
gdept = mesh.variables['gdept_0'][0,:,:,:]
e3t = mesh.variables['e3t_0'][0,:,:,:]
gdept_1d = mesh.variables['gdept_1d'][0,:]
# multi dimensional because of patrial cells
print gdept[:,0,0]
print gdept_1d
[ 0.50100231 1.51016855 2.53205109 3.57176328 4.6364336 5.73595858 6.88400364 8.09929371 9.40723228 10.84182644 12.44784164 14.28297424 16.41976166 18.94689369 21.9696846 25.6096077 30.0026474 35.29554749 41.63795471 49.16830826 57.99394608 68.17064667 79.69096375 92.4890213 106.46230316 121.5032959 137.53115845 154.51651001 172.49569702 191.5723114 211.90315247 233.66589355 257.01101685 282.01119995 308.62915039 336.71963501 366.06262207 396.40994263 427.52453613 459.20376587] [ 0.50100231 1.51016859 2.53205114 3.57176335 4.63643351 5.73595859 6.88400344 8.0992936 9.40723195 10.84182681 12.44784199 14.28297417 16.41976131 18.94689327 21.96968528 25.60960816 30.00264702 35.29554825 41.63795305 49.16830922 57.99394791 68.17065039 79.69096548 92.48901994 106.4623063 121.50329221 137.53115246 154.51651255 172.49569794 191.57230836 211.90315488 233.66589683 257.01101119 282.01118741 308.62915991 336.71962219 366.0626321 396.4099436 427.52452947 459.20375832]
Plot to check
plt.plot(e3t[:,0,0], gdept_1d)
plt.grid()
plt.ylim([450,0])
plt.xlabel('grid spacing (m)')
plt.ylabel('depth (m)')
<matplotlib.text.Text at 0x7fcbdbf0a950>
Interpolate
Tnew = np.copy(Told)
Snew = np.copy(Sold)
for j in np.arange(Tnew.shape[2]):
for i in np.arange(Tnew.shape[3]):
Tnew[0,:,j,i] = np.interp(gdept[:,j,i], gdept_old[:,j,i], Told[0,:,j,i])
Snew[0,:,j,i] = np.interp(gdept[:,j,i], gdept_old[:,j,i], Sold[0,:,j,i])
Check it makes sense
vmin=0;vmax=14
fig,axs =plt.subplots(1,2)
j=5
mesh=axs[0].pcolormesh(np.arange(Tnew.shape[-1]),gdept[:,j,:],Tnew[0,:,j,:],vmin=vmin,vmax=vmax)
plt.colorbar(mesh,ax=axs[0])
axs[0].set_title('New')
mesh=axs[1].pcolormesh(np.arange(Told.shape[-1]),gdept_old[:,j,:],Told[0,:,j,:],vmin=vmin,vmax=vmax)
plt.colorbar(mesh,ax=axs[1])
axs[1].set_title('Old')
<matplotlib.text.Text at 0x7fcbdbec9c50>
vmin=20;vmax=32
fig,axs =plt.subplots(1,2)
j=5
mesh=axs[0].pcolormesh(np.arange(Snew.shape[-1]),gdept[:,j,:],Snew[0,:,j,:],vmin=vmin,vmax=vmax)
plt.colorbar(mesh,ax=axs[0])
axs[0].set_title('New')
mesh=axs[1].pcolormesh(np.arange(Sold.shape[-1]),gdept_old[:,j,:],Sold[0,:,j,:],vmin=vmin,vmax=vmax)
plt.colorbar(mesh,ax=axs[1])
axs[1].set_title('Old')
<matplotlib.text.Text at 0x7fcbdbc7f650>
Save in a new file
def save_interopolated_netcdf_initial(filename):
"""Save the new interpolate T+S initial conditions in a netcdf file"""
new_TS = nc.Dataset(filename, 'w')
nc_tools.init_dataset_attrs(
new_TS,
title='2D T+S initialization - interpolated',
notebook_name='Changing Resolution-trpltanh.ipynb',
nc_filepath=filename,
comment='Interpolated T+S initial condtions for alternate vertical grid spacing - trpl tanh')
new_TS.createDimension('y', 10)
new_TS.createDimension('x', 1100)
new_TS.createDimension('deptht',size = len(gdept_old))
new_TS.createDimension('time_counter', None)
nc_tools.show_dimensions(new_TS)
# variables
y = new_TS.createVariable('y', 'float32', ('y','x'))
y.long_name = 'spanwise'
y.units = 'km'
y[:] =yold
x = new_TS.createVariable('x', 'float32', ('y','x'))
x.long_name = 'streamwise'
x.units = 'km'
x[:] = xold
deptht = new_TS.createVariable('deptht', 'float32', ('deptht'))
deptht.long_name = 'Depth'
deptht.units = 'm'
deptht.positive = 'down'
deptht[:] = gdept_1d
time_counter = new_TS.createVariable('time_counter', 'float32', ('time_counter'))
time_counter.units = 'seconds since 2003-08-09 0:00:00'
time_counter.long_name = 'Time axis'
vosaline = new_TS.createVariable('vosaline', 'float32',
('time_counter','deptht','y','x'))
vosaline.units = 'none'
vosaline.long_name = 'Practical Salinity'
vosaline.coordinates = 'x y deptht time_counter'
vosaline.grid = '2D bathy'
vosaline[0] =Snew
votemper = new_TS.createVariable('votemper', 'float32',
('time_counter','deptht','y','x'))
votemper.units = 'degC'
votemper.long_name = 'Temperature'
votemper.coordinates = 'x y deptht time_counter'
votemper[0] = Tnew
new_TS.close()
save_interopolated_netcdf_initial('/data/nsoontie/MEOPAR/2Ddomain/vertical_resolution/trpl_tanh/TS_uniform_36.nc')
file format: NETCDF4 Conventions: CF-1.6 title: 2D T+S initialization - interpolated institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia source: https://bitbucket.org/salishsea/2d-domain/src/tip/Changing Resolution-trpltanh.ipynb references: REQUIRED history: [2015-10-09 11:12:53] Created netCDF4 zlib=True dataset. comment: Interpolated T+S initial condtions for alternate vertical grid spacing - trpl tanh <type 'netCDF4._netCDF4.Dimension'>: name = 'y', size = 10 <type 'netCDF4._netCDF4.Dimension'>: name = 'x', size = 1100 <type 'netCDF4._netCDF4.Dimension'>: name = 'deptht', size = 40 <type 'netCDF4._netCDF4.Dimension'> (unlimited): name = 'time_counter', size = 0
Quick check that variables were saved.
f=nc.Dataset('/data/nsoontie/MEOPAR/2Ddomain/vertical_resolution/trpl_tanh/TS_uniform_36.nc')
depth=f.variables['deptht'][:]
x=f.variables['x'][:]
S = f.variables['vosaline'][:]
print depth
print x
print S.shape, S.min(), S.max()
[ 0.50100231 1.51016855 2.53205109 3.57176328 4.6364336 5.73595858 6.88400364 8.09929371 9.40723228 10.84182644 12.44784164 14.28297424 16.41976166 18.94689369 21.9696846 25.6096077 30.0026474 35.29554749 41.63795471 49.16830826 57.99394608 68.17064667 79.69096375 92.4890213 106.46230316 121.5032959 137.53115845 154.51651001 172.49569702 191.5723114 211.90315247 233.66589355 257.01101685 282.01119995 308.62915039 336.71963501 366.06262207 396.40994263 427.52453613 459.20376587] [[ 0.00000000e+00 5.00000000e-01 1.00000000e+00 ..., 5.48500000e+02 5.49000000e+02 5.49500000e+02] [ 0.00000000e+00 5.00000000e-01 1.00000000e+00 ..., 5.48500000e+02 5.49000000e+02 5.49500000e+02] [ 0.00000000e+00 5.00000000e-01 1.00000000e+00 ..., 5.48500000e+02 5.49000000e+02 5.49500000e+02] ..., [ 0.00000000e+00 5.00000000e-01 1.00000000e+00 ..., 5.48500000e+02 5.49000000e+02 5.49500000e+02] [ 0.00000000e+00 5.00000000e-01 1.00000000e+00 ..., 5.48500000e+02 5.49000000e+02 5.49500000e+02] [ 0.00000000e+00 5.00000000e-01 1.00000000e+00 ..., 5.48500000e+02 5.49000000e+02 5.49500000e+02]] (1, 40, 10, 1100) 0.0 31.3702
The forcing files for temperature/salinity all need to be interpolated onto the new vertical grid as well.
Open old TS boundary condition
#open old file
fold = nc.Dataset('/data/nsoontie/MEOPAR/2Ddomain/boundary_conditions/TS_OBC.nc')
Sold = fold.variables['vosaline'][:]
Told = fold.variables['votemper'][:]
print fold
times = fold.variables['time_counter'][:]
#Will use gdept_old for depths
<type 'netCDF4._netCDF4.Dataset'> root group (NETCDF4 data model, file format UNDEFINED): Conventions: CF-1.6 title: Temperature and Salinty Boundary Conditions 2D domain institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia source: https://bitbucket.org/salishsea/2d-domain/src/tip/Generate T+S Forcing - NEMO3.6.ipynb references: REQUIRED history: [2015-08-05 16:23:30] Created netCDF4 zlib=True dataset. comment: based on average values across mouth of JdF and 3D weekly climatology dimensions(sizes): xb(80), yb(1), time_counter(52), deptht(40) variables(dimensions): float32 deptht(deptht), float32 time_counter(time_counter), float32 votemper(time_counter,deptht,yb,xb), float32 vosaline(time_counter,deptht,yb,xb), int32 nbidta(yb,xb), int32 nbjdta(yb,xb), int32 nbrdta(yb,xb) groups:
Interpolate onto new grid - gdept_1d
Tnew = np.copy(Told)
Snew = np.copy(Sold)
for t in np.arange(Tnew.shape[0]):
for j in np.arange(Tnew.shape[2]):
for i in np.arange(Tnew.shape[3]):
Tnew[t,:,j,i] = np.interp(gdept_1d[:], gdept_1d_old[:], Told[t,:,j,i])
Snew[t,:,j,i] = np.interp(gdept_1d[:], gdept_1d_old[:], Sold[t,:,j,i])
Quick check
vmin=0;vmax=14
fig,axs =plt.subplots(1,2)
mesh=axs[0].pcolormesh(np.arange(Tnew.shape[0]),gdept_1d[:],Tnew[:,:,0,10].T,vmin=vmin,vmax=vmax)
plt.colorbar(mesh,ax=axs[0])
axs[0].set_title('New')
mesh=axs[1].pcolormesh(np.arange(Told.shape[0]),gdept_1d_old[:],Told[:,:,0,10].T,vmin=vmin,vmax=vmax)
plt.colorbar(mesh,ax=axs[1])
axs[1].set_title('Old')
<matplotlib.text.Text at 0x7fcbdb0d19d0>
vmin=28;vmax=34
fig,axs =plt.subplots(1,2)
mesh=axs[0].pcolormesh(np.arange(Snew.shape[0]),gdept_1d[:],Snew[:,:,0,10].T,vmin=vmin,vmax=vmax)
plt.colorbar(mesh,ax=axs[0])
axs[0].set_title('New')
mesh=axs[1].pcolormesh(np.arange(Sold.shape[0]),gdept_1d_old[:],Sold[:,:,0,10].T,vmin=vmin,vmax=vmax)
plt.colorbar(mesh,ax=axs[1])
axs[1].set_title('Old')
<matplotlib.text.Text at 0x7fcbdae17710>
Save netcdf
def save_interopolated_netcdf_boundary(filename):
Ny=8#remeber masked edges
nrim=10
ndepth=gdept_1d.shape[0]
nemo = nc.Dataset(filename, 'w', zlib=True)
#start and end points
length_rim =nrim
lengthi=Ny*length_rim
#time and depth
depth_levels =ndepth
# dataset attributes
nc_tools.init_dataset_attrs(
nemo,
title='Temperature and Salinty Boundary Conditions 2D domain - interpolated to new vertical grid - trpl tanh',
notebook_name='Changing resolution-trpltanh',
nc_filepath=filename,
comment='intepolated for new vertical grid spacing')
# dimensions
nemo.createDimension('xb', lengthi)
nemo.createDimension('yb', 1)
nemo.createDimension('time_counter', None)
nemo.createDimension('deptht', depth_levels)
# variables
# deptht
deptht = nemo.createVariable('deptht', 'float32', ('deptht',))
deptht.long_name = 'Vertical T Levels'
deptht.units = 'm'
deptht.positive = 'down'
deptht[:]=gdept_1d[:]
# time_counter
time_counter = nemo.createVariable('time_counter', 'float32', ('time_counter'))
time_counter.long_name = 'Time axis'
time_counter.axis = 'T'
time_counter.units = 'weeks since beginning of year'
time_counter[:]=times[:]
# votemper
votemper = nemo.createVariable('votemper', 'float32',
('time_counter','deptht','yb','xb'))
votemper.units = 'degC'
votemper.long_name = 'Temperature'
votemper[:]=Tnew[:]
# vosaline
vosaline = nemo.createVariable('vosaline', 'float32',
('time_counter','deptht','yb','xb'))
vosaline.units = 1
vosaline.long_name = 'Practical Salinity'
vosaline[:]=Snew
# nbidta, ndjdta, ndrdta
nbidta = nemo.createVariable('nbidta', 'int32' , ('yb','xb'))
nbidta.long_name = 'i grid position'
nbidta.units = 1
nbjdta = nemo.createVariable('nbjdta', 'int32' , ('yb','xb'))
nbjdta.long_name = 'j grid position'
nbjdta.units = 1
nbrdta = nemo.createVariable('nbrdta', 'int32' , ('yb','xb'))
nbrdta.long_name = 'position from boundary'
nbrdta.units = 1
for ir in range(length_rim):
nbidta[0,ir*Ny:(ir+1)*Ny] = ir
nbjdta[0,ir*Ny:(ir+1)*Ny] = range(Ny)
nbrdta[0,ir*Ny:(ir+1)*Ny] = ir
nemo.close()
save_interopolated_netcdf_boundary('/data/nsoontie/MEOPAR/2Ddomain/vertical_resolution/trpl_tanh/TS_OBC.nc')
file format: NETCDF4 Conventions: CF-1.6 title: Temperature and Salinty Boundary Conditions 2D domain - interpolated to new vertical grid - trpl tanh institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia source: https://bitbucket.org/salishsea/2d-domain/src/tip/Changing resolution-trpltanh.ipynb references: REQUIRED history: [2015-10-09 11:15:58] Created netCDF4 zlib=True dataset. comment: intepolated for new vertical grid spacing
Quick check
f=nc.Dataset('/data/nsoontie/MEOPAR/2Ddomain/vertical_resolution/trpl_tanh/TS_OBC.nc')
depth=f.variables['deptht'][:]
sal=f.variables['vosaline'][:]
print sal.shape, sal.min(), sal.max()
print depth
print f
(52, 40, 1, 80) 29.5012 33.6646 [ 0.50100231 1.51016855 2.53205109 3.57176328 4.6364336 5.73595858 6.88400364 8.09929371 9.40723228 10.84182644 12.44784164 14.28297424 16.41976166 18.94689369 21.9696846 25.6096077 30.0026474 35.29554749 41.63795471 49.16830826 57.99394608 68.17064667 79.69096375 92.4890213 106.46230316 121.5032959 137.53115845 154.51651001 172.49569702 191.5723114 211.90315247 233.66589355 257.01101685 282.01119995 308.62915039 336.71963501 366.06262207 396.40994263 427.52453613 459.20376587] <type 'netCDF4._netCDF4.Dataset'> root group (NETCDF4 data model, file format UNDEFINED): Conventions: CF-1.6 title: Temperature and Salinty Boundary Conditions 2D domain - interpolated to new vertical grid - trpl tanh institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia source: https://bitbucket.org/salishsea/2d-domain/src/tip/Changing resolution-trpltanh.ipynb references: REQUIRED history: [2015-10-09 11:15:58] Created netCDF4 zlib=True dataset. comment: intepolated for new vertical grid spacing dimensions(sizes): xb(80), yb(1), time_counter(52), deptht(40) variables(dimensions): float32 deptht(deptht), float32 time_counter(time_counter), float32 votemper(time_counter,deptht,yb,xb), float32 vosaline(time_counter,deptht,yb,xb), int32 nbidta(yb,xb), int32 nbjdta(yb,xb), int32 nbrdta(yb,xb) groups:
All done! Now try to run with the new vertical resolution!