Skip to content

GitLab

Commit 028a045f authored by Veit Lüschow's avatar Veit Lüschow
Browse files

latest

parent bc9ca6d1
Hide whitespace changes
Inline Side-by-side
Showing with 85 additions and 132 deletions
side_tools.py
View file @ 028a045f
......@@ -17,6 +17,6 @@ def clean():
print(test)
for item in test:
if item.endswith(".nc") or item.endswith('.h5'):
if item.endswith(".nc"):
os.remove( os.path.join( directory, item ) )
print('deleted ', item)
\ No newline at end of file
veros/setup/merjet/merjet.py
View file @ 028a045f
#!/usr/bin/env python
#!/usr/bin/env veros
from IPython import get_ipython
get_ipython().magic('reset -sf')
......@@ -28,12 +28,15 @@ class MERJETSetup(VerosSetup):
@veros_method
def set_parameter(self, vs):
vs.identifier = 'merjet'
vs.nx, vs.ny, vs.nz = 50, 35, 60
vs.dt_mom = 150 # normally 2400
vs.dt_tracer = 86400 / 32
vs.runlen = 86400 * 365 * 1
vs.runlen = 86400 * 100
# vs.restart_input_filename='merjet_9637.restart.h5'
vs.nx, vs.ny, vs.nz = 150, 30, 30
vs.dt_mom = 75 # normally 2400
vs.dt_tracer = 86400 / 16
vs.runlen = 86400 * 365 * 2
# vs.runlen = 86400 * 300
# vs.restart_frequency = 86400 * 100
# vs.runlen = 100
vs.coord_degree = True
......@@ -57,7 +60,7 @@ class MERJETSetup(VerosSetup):
# print(vs.A_h)
# raise SystemExit(0)
vs.enable_bottom_friction = False
vs.enable_bottom_friction = True
vs.r_bot = 1e-4
vs.enable_noslip_lateral = True
......@@ -91,20 +94,25 @@ class MERJETSetup(VerosSetup):
vs.eq_of_state_type = 1
vs.force_overwrite = True # added for easier configuring the setting
vs.enable_tempsalt_sources = True
@veros_method
def set_grid(self, vs):
ddz = np.array([50., 70., 100., 140., 190., 240., 290., 340.,
390., 440., 490., 540., 590., 640., 690.])
vs.dxt[...] = 0.1
vs.dyt[...] = 0.1
vs.dzt[...] = np.array([240.,200.,160,140,120,100,90,80,70,60,50,40,40,40,40, \
40,40,40,40,50,60,70,80,90,100,120,140,160,200,240])
vs.dxt[...] = 0.25
vs.dyt[...] = 0.25
vs.x_origin = 5.0
vs.y_origin = 10
# vs.dzt[...] = ddz[::-1] / 2.5
vs.dzt[:] = 60
# import matplotlib.pyplot as plt
# print(ddz.shape)
# plt.plot(vs.dzt)
# plt.plot(np.cumsum(vs.dzt))
# plt.show()
# raise SystemExit(0)
@veros_method
......@@ -112,39 +120,32 @@ class MERJETSetup(VerosSetup):
vs.coriolis_t[:, :] = 2 * vs.omega * np.sin(vs.yt[None, :] / 180. * vs.pi)
if vs.enable_rotate_f:
xtmp = vs.xt.copy()
xtmp[1] = xtmp[2]-0.1
xtmp[0] = xtmp[2]- 2 * 0.1
xtmp[-1] = xtmp[-3] + 2 * 0.1
xtmp[-2] = xtmp[-3] + 0.1
xtmp[1] = xtmp[2]-0.25
xtmp[0] = xtmp[2]- 2 * 0.25
xtmp[-1] = xtmp[-3] + 2 * 0.25
xtmp[-2] = xtmp[-3] + 0.25
vs.coriolis_t = 2 * vs.omega * np.sin(np.tile(xtmp[:,None],(1,vs.ny+4)) / 180. * vs.pi)
# print((np.sin(vs.xt[:,None] /180*vs.pi)).shape)
# print(vs.coriolis_t.shape)
# import matplotlib.pyplot as plt
# plt.contourf(vs.coriolis_t)
# # plt.plot(xtmp)
# print(vs.xt)
# print(vs.yt)
# plt.colorbar()
# plt.plot(vs.zt)
# raise SystemExit(0)
@veros_method
def set_topography(self, vs):
x, y = np.meshgrid(vs.xt, vs.yt, indexing='ij')
vs.kbot[:,:-4] = np.int(2)
vs.kbot[:,-4:] = np.int(32)
vs.kbot[:,-5] = np.int(24)
vs.kbot[:,-6] = np.int(16)
#%% steep slope topography, uniform
vs.kbot[:,-4:] = np.int(14)
vs.kbot[:,-5] = np.int(8)
vs.kbot[:,-6] = np.int(6)
vs.kbot[:,-7:-6] = np.int(4)
vs.kbot[:,-8:-7] = np.int(3)
vs.kbot[:,-4:] = np.int(0)
vs.kbot[:,-8:-7] = np.int(8)
# vs.kbot[-20:,7:8] = np.int(8)
# print(vs.yt)
# vs.kbot[20:80,4:] = np.int(2)
#
# vs.kbot[:,:] = np.int(1)
# import matplotlib.pyplot as plt
# plt.contourf(vs.kbot)
# plt.colorbar()
......@@ -153,22 +154,20 @@ class MERJETSetup(VerosSetup):
@veros_method
def set_initial_conditions(self, vs):
# initial conditions
# design initial profile for u that is also used for restoring
zcenter = 1900
zvariance = 6e5
yvariance = 0.07
ycenter = 11.7 # normally 12.7
"""
Flow initial and restoring profile
"""
zcenter = 1500
zvariance = 3e5
yvariance = 0.7
ycenter = 15.5 # normally 16
yy,zz = np.meshgrid(vs.yt,vs.zt,indexing='ij')
uz = np.exp(- ( (zz+(zcenter) )**2 / zvariance ) ) # exponential decay in z
# uz = (vs.zt[0]-zz)/vs.zt[0]
uy = np.exp(- ( (yy-(ycenter) )**2 / yvariance ) ) # gaussian in y
uu = uz * uy
vs.u_ini = allocate(vs, ('xu','yt','zt'))
vs.u_ini = -1 * np.tile(uu[np.newaxis,:,:],[vs.xu.shape[0] , 1 , 1])
vs.u_ini = -0.5 * np.tile(uu[np.newaxis,:,:],[vs.xu.shape[0] , 1 , 1])
# import matplotlib.pyplot as plt
......@@ -185,97 +184,41 @@ class MERJETSetup(VerosSetup):
vs.temp[:, :, :, 0:2] = ((1 - vs.zt[None, None, :] / vs.zw[0]) * 15 * vs.maskT)[..., None]
vs.salt[:, :, :, 0:2] = 35.0 * vs.maskT[..., None]
temp = vs.temp[...,0].copy()
salt = vs.salt[...,0].copy()
# rho_back = ( 1 - np.exp( vs.zt * 1/zscale )) * 2
rho_back = density.get_rho(vs, salt, temp, np.abs(vs.zt)) * vs.maskT
vs.temp_ini = vs.temp[:,:,:,0].copy()
# dzu = np.gradient(vs.u_ini,vs.zt,axis=2) * np.tile(vs.coriolis_t[...,np.newaxis],[1,1,vs.nz]) * 1024 / 10
dzu = np.gradient(vs.u_ini,vs.zt,axis=2) * vs.coriolis_t[0,np.int(ycenter)] * 1024 / 10
# cumsum to get density field
#
#print(vs.dyt )
sumi = allocate(vs,('xt','yt','zt'))
sumi[:,2:-2,:] = np.cumsum(dzu[:,2:-2,:]*10000,axis=1)+ rho_back[:,2:-2,:]
"""
Restoring timescales
"""
vs.rho_ini = allocate(vs, ('xt', 'yt', 'zt'))
vs.rho_ini = sumi.copy()
xtmp = np.zeros(vs.xt.shape[0])
xtmp[3:-2] = vs.xt[3:-2]
vs.temp_ini = vs.temp[:,:,:,0].copy()
# t_tmp = 2-np.tanh((xtmp[-3]-xtmp)/1) + np.tanh((-xtmp)/2)
t_tmp = 1-np.tanh((xtmp[-3]-xtmp[:])/2)
t_tmp [-2:] = t_tmp[-3]
# vs.temp[...,0:2] = vs.temp_ini[...,None] * vs.maskT[..., None]
ft = np.flipud(t_tmp)
# vs.salt[:, :, :, 0:2] = 35.0 * vs.maskT[..., None]
vs.t_restoring = allocate(vs,('xt','yt','zt'))
vs.t_restoring_tr = allocate(vs,('xt','yt','zt'))
vs.t_restoring = 1 / (1e2 * vs.dt_mom) * np.tile((((ft))[...,np.newaxis,np.newaxis]),[1,vs.ny+4,vs.nz])
# rho = density.get_rho(vs, 35, vs.temp[...,1], np.abs(vs.zt))
vs.t_restoring_tr[:,:,-1] = 2e-5
vs.t_restoring_tr[:,:,-2] = 1e-6
# import matplotlib.pyplot as plt
# # # plt.plot(vs.zt,rho_back)
# # # plt.figure()
# plt.contourf(yy,zz,vs.rho_ini[10,:,:]-rho_back[10,:,:])
# plt.plot(ft)
# plt.contourf(vs.t_restoring[:,:,1])
# plt.colorbar()
# # # plt.contourf(vs.rho_ini[0,:,:])
# raise SystemExit(0)
# wind stress forcing
yt_min = global_min(vs, vs.yt.min())
yu_min = global_min(vs, vs.yu.min())
yt_max = global_max(vs, vs.yt.max())
yu_max = global_max(vs, vs.yu.max())
taux = allocate(vs, ('yt',))
vs.surface_taux[:, :] = 0
# surface heatflux forcing
vs._t_star = allocate(vs, ('yt',), fill=15)
vs._t_star = 0
vs._t_rest = 0
@veros_method
def set_forcing(self, vs):
#restoring timescale only
xtmp = np.zeros(vs.xu.shape[0])
xtmp[3:-2] = vs.xu[3:-2]
t_tmp = 2-np.tanh((xtmp[-3]-xtmp)/1) + np.tanh((-xtmp)/2)
t_tmp [-2:] = t_tmp[-3]
# t_tmp = np.tanh((-xtmp)/5)
# t_tmp [-3:] = t_tmp[-4]
# t_tmp = np.exp(- ( (xtmp-35 )**2 / 8 ) ) # gaussian in y
t_restoring = 1 / (1e2 * vs.dt_mom) * np.tile((t_tmp[...,np.newaxis,np.newaxis]),[1,vs.ny+4,vs.nz])
t_restoring_tr = 1e-5 * np.tile((t_tmp[...,np.newaxis,np.newaxis]),[1,vs.ny+4,vs.nz])
t_restoring_tr[:,:,vs.nz-1] = 2e-5
# t_restoring_tr[:,2,:] = 1e-5
vs.du_rest[...] = vs.maskU * vs.t_restoring * (vs.u_ini[...] - vs.u[...,vs.tau])
# vs.dv_rest = 0.
# t_restoring = 1 / (1000 * vs.dt_mom)
# t_restoring_tr[3:4,:,:] = 1e-10
vs.du_rest[...] = vs.maskU * t_restoring * (vs.u_ini[...] - vs.u[...,vs.tau])
# print('hi')
# tscl = allocate(vs,('xt','yt','zt'))
# tscl [3:4,:,:] = 1e-7
# vs.u_sources = 0.
# vs.u_sources = vs.maskU * t_restoring * (vs.u_ini - vs.u[...,vs.taum1])
# vs.temp[... , vs.tau] = vs.temp[... , vs.tau]
# + t_restoring_tr * (vs.temp_ini - vs.temp[... , vs.taum1])
vs.temp_source[...] = vs.maskT * t_restoring_tr * (vs.temp_ini - vs.temp[... , vs.tau])
vs.salt_source[...] = vs.maskT * t_restoring_tr * (35 - vs.salt[... , vs.tau])
# import matplotlib.pyplot as plt
# plt.plot(t_tmp)
# plt.contourf(t_restoring_tr[:,:,1])
# plt.colorbar()
# raise SystemExit(0)
vs.temp_source[...] = vs.maskT * vs.t_restoring_tr * (vs.temp_ini - vs.temp[... , vs.tau])
vs.salt_source[...] = vs.maskT * vs.t_restoring_tr * (35 - vs.salt[... , vs.tau])
@veros_method
......@@ -292,7 +235,7 @@ class MERJETSetup(VerosSetup):
vs.diagnostics['averages'].sampling_frequency = vs.dt_tracer * 10.
# vs.diagnostics['overturning'].output_frequency = 365 * 86400.
# vs.diagnostics['overturning'].sampling_frequency = vs.dt_tracer * 1000
# vs.diagnostics['tracer_monitor'].output_frequency = 365 * 86400.
vs.diagnostics['tracer_monitor'].output_frequency = 10 * 86400.
# vs.diagnostics['energy'].output_frequency = 365 * 86400.
# vs.diagnostics['energy'].sampling_frequency = vs.dt_tracer * 10
#vs.diagnostics['snapshot'].output_variables = (
......
veros/setup/merjet/plot.py
View file @ 028a045f
......@@ -58,24 +58,34 @@ except NameError:
pass
# p = section('../acc/acc.averages.nc','u')
# p = section('good/20200401_1.averages.nc', 'u')
# p = section('/Users/veitlueschow/veros_exps/merjet-good/20200401_4.averages.nc', 'u')
# #%%
# v = xr.open_dataset('merjet.snapshot.nc').v.isel(zt=30)
# u = xr.open_dataset('merjet.snapshot.nc').u.isel(zt=30)
# zeta = np.gradient(u,axis=1) - np.gradient(v,axis=2)
# #%%
# plt.contourf(zeta[0,:,:],np.linspace(-0.3,0.3,10),cmap = 'RdBu_r',extend='both')
# plt.colorbar()
p = section('merjet.snapshot.nc', 'u')
#%%
p.select(2,'zt',30).plot.contourf(levels=15)
p = section('merjet.snapshot.nc', 'u')
p.select(133,'zt',15).plot.contourf(levels=15)
#%%
p = section('merjet.averages.nc', 'u')
step = 3
p = section('merjet.snapshot.nc', 'u')
step = 1
xpos = 5
plt.figure()
p.select(step,'xu',xpos).plot.contour(levels=15,vmin=-0.6,vmax=0.6,colors="black")
p.select(step,'xu',xpos).plot.contour(levels=15,vmin=-0.2,vmax=0.2,colors="black")
#
p = section('merjet.averages.nc', 'rho')
p = section('merjet.snapshot.nc', 'rho')
p.select(step,'xt',xpos).plot.contourf(levels=16,cmap='RdBu_r')
"""
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment