This notebook outlines the steps needed to change the resolution of the 2D domain.
import numpy as np
import matplotlib.pyplot as plt
import netCDF4 as nc
from salishsea_tools import nc_tools
%matplotlib inline
To change the vertical resolution, you need to work through the followinh steps.
The vertical resolution is set through a number of options in the namelist.domain.
ppsur = 999999.
ppa0 = 999999.
ppa1 = 999999.
ppkth = 25.
ppacr = 3.
ldbletanh = .FALSE. ! use double tanh true or false
ppa2 = 0 Double tanh function parameters
ppkth2 = 0
ppacr2 = 0
These parameters set the grid spacing and depth at each grid level as defined by:
$$ z(k) = ppsur + ppa0k + ppa1 ppacr \log(\cosh(\frac{k-ppkth}{ppacr})) + ppa2 ppacr2 \log\cosh(\frac{k-ppkth2} {ppacr2})) $$$$ e3(k) = ppa0k + ppa1 \tanh(\frac{k-ppkth}{ppacr}) + ppa2\tanh(\frac{k-ppkth2}{ppacr2}) $$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). Note that the second hyperbolic tan is only used if ldbletanh=.TRUE.
Note: ppsur should be set so that z(1) = 0
You should experirment with these settings in a notebook. Here is an example of a double tanh profile.
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))
ks = np.arange(1,41) # Remember Fortan starts counting at 1.
#Parameters
ppa1 = 0.5
ppacr = 2.
ppkth = 12
ppa2 = 6.7
ppacr2 =2.2
ppkth2 = 13.
#ensure first grid cell has spacing 1
ppa0 = 1 -NEMO_tanh(ks[0], ppa1, ppkth, ppacr) - NEMO_tanh(ks[0], ppa2, ppkth2, ppacr2)
#ensure first grid point is at the surface (0)
ppsur = - (ppa0 + NEMO_tanh_integral(ks[0],ppa1,ppkth,ppacr) + NEMO_tanh_integral(ks[0],ppa2,ppkth2,ppacr2))
#define grid spacing (dzs) and grid levels (zs)
dzs = ppa0 + NEMO_tanh(ks, ppa1, ppkth, ppacr) + NEMO_tanh(ks, ppa2, ppkth2, ppacr2)
zs = ppsur + ppa0*ks + NEMO_tanh_integral(ks,ppa1,ppkth,ppacr) + NEMO_tanh_integral(ks,ppa2,ppkth2,ppacr2)
#plot
fig,axs=plt.subplots(1,3,figsize=(15,5))
for ax in axs:
ax.plot(dzs,zs,'-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: ',zs[0]
print 'First grid cell thickness', dzs[0]
print 'Bottom:', zs[-1]
ppsur -83.1898878903 ppa0 8.19973820108 Surface: 1.42108547152e-14 First grid cell thickness 1.0 Bottom: 428.789503531
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
Copy this file into the workspace where you'll be interpolating the boundary conditions and initial conditions.
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,:]
xold = f.variables['nav_lon'][:]
yold = f.variables['nav_lat'][:]
print gdept_old
[ 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]
Open new mesh file
mesh = nc.Dataset('/data/nsoontie/MEOPAR/2Ddomain/vertical_resolution/dbl_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.50003511 1.50045013 2.50187922 3.50584388 4.51615047 5.54231215 6.60804415 7.77197409 9.17587852 11.14406872 14.32758141 19.71521759 28.14229202 39.53160858 52.97109604 67.5062027 82.54394531 97.79561615 113.13538361 128.51092529 143.90090942 159.29672241 174.69488525 190.09399414 205.49346924 220.89311218 236.2928009 251.69252014 267.09225464 282.49200439 297.89172363 313.29147339 328.69119263 344.09094238 359.49069214 374.89041138 390.29016113 405.68991089 421.08963013 436.48937988] [ 0.50003508 1.50045009 2.50187922 3.50584399 4.51615035 5.54231202 6.60804398 7.77197408 9.17587888 11.14406885 14.32758163 19.71521847 28.14229296 39.53160899 52.97109699 67.50620429 82.54394591 97.79561875 113.13538312 128.5109245 143.90091533 159.29672913 174.69488759 190.0939899 205.49347212 220.89310727 236.29280398 251.69252548 267.09225695 282.49199245 297.89172956 313.29146732 328.69120534 344.09094347 359.49068165 374.89041983 390.29015803 405.68989623 421.08963443 436.48937263]
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, Told[0,:,j,i])
Snew[0,:,j,i] = np.interp(gdept[:,j,i], gdept_old, Sold[0,:,j,i])
Check it makes sense
vmin=0;vmax=14
fig,axs =plt.subplots(1,2)
mesh=axs[0].pcolormesh(np.arange(Tnew.shape[-1]),gdept[:,5,:],Tnew[0,:,5,:],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[:],Told[0,:,5,:],vmin=vmin,vmax=vmax)
plt.colorbar(mesh,ax=axs[1])
axs[1].set_title('Old')
<matplotlib.text.Text at 0x7f75e079b1d0>
vmin=20;vmax=32
fig,axs =plt.subplots(1,2)
mesh=axs[0].pcolormesh(np.arange(Snew.shape[-1]),gdept[:,5,:],Snew[0,:,5,:],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[:],Sold[0,:,5,:],vmin=vmin,vmax=vmax)
plt.colorbar(mesh,ax=axs[1])
axs[1].set_title('Old')
<matplotlib.text.Text at 0x7f75e02b5d10>
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.ipynb',
nc_filepath=filename,
comment='Interpolated T+S initial condtions for alternate vertical grid spacing')
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/dbl_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.ipynb references: REQUIRED history: [2015-09-04 11:43:14] Created netCDF4 zlib=True dataset. comment: Interpolated T+S initial condtions for alternate vertical grid spacing <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/dbl_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.50003511 1.50045013 2.50187922 3.50584388 4.51615047 5.54231215 6.60804415 7.77197409 9.17587852 11.14406872 14.32758141 19.71521759 28.14229202 39.53160858 52.97109604 67.5062027 82.54394531 97.79561615 113.13538361 128.51092529 143.90090942 159.29672241 174.69488525 190.09399414 205.49346924 220.89311218 236.2928009 251.69252014 267.09225464 282.49200439 297.89172363 313.29147339 328.69119263 344.09094238 359.49069214 374.89041138 390.29016113 405.68991089 421.08963013 436.48937988] [[ 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_old[:], Told[t,:,j,i])
Snew[t,:,j,i] = np.interp(gdept_1d[:], gdept_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_old[:],Told[:,:,0,10].T,vmin=vmin,vmax=vmax)
plt.colorbar(mesh,ax=axs[1])
axs[1].set_title('Old')
<matplotlib.text.Text at 0x7f75dfaa5b90>
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_old[:],Sold[:,:,0,10].T,vmin=vmin,vmax=vmax)
plt.colorbar(mesh,ax=axs[1])
axs[1].set_title('Old')
<matplotlib.text.Text at 0x7f75df7f3e90>
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',
notebook_name='Changing resolution',
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/dbl_tanh/TS_OBC.nc')
file format: NETCDF4 Conventions: CF-1.6 title: Temperature and Salinty Boundary Conditions 2D domain - interpolated to new vertical grid institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia source: https://bitbucket.org/salishsea/2d-domain/src/tip/Changing resolution.ipynb references: REQUIRED history: [2015-09-04 11:43:15] Created netCDF4 zlib=True dataset. comment: intepolated for new vertical grid spacing
Quick check
f=nc.Dataset('/data/nsoontie/MEOPAR/2Ddomain/vertical_resolution/dbl_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.5011 33.6646 [ 0.50003511 1.50045013 2.50187922 3.50584388 4.51615047 5.54231215 6.60804415 7.77197409 9.17587852 11.14406872 14.32758141 19.71521759 28.14229202 39.53160858 52.97109604 67.5062027 82.54394531 97.79561615 113.13538361 128.51092529 143.90090942 159.29672241 174.69488525 190.09399414 205.49346924 220.89311218 236.2928009 251.69252014 267.09225464 282.49200439 297.89172363 313.29147339 328.69119263 344.09094238 359.49069214 374.89041138 390.29016113 405.68991089 421.08963013 436.48937988] <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 institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia source: https://bitbucket.org/salishsea/2d-domain/src/tip/Changing resolution.ipynb references: REQUIRED history: [2015-09-04 11:43:15] 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!