Skip to content

GitLab

Commit 86cb7ee8 authored by Veit Lüschow's avatar Veit Lüschow
Browse files

about to iniotialize density in geostrophic balance

parent d1b0be36
Hide whitespace changes
Inline Side-by-side
Showing with 73 additions and 54 deletions
veros/core/streamfunction/solve_stream.py
View file @ 86cb7ee8
......@@ -64,7 +64,7 @@ def solve_streamfunction(vs):
line_forc = allocate(vs, ('isle',))
if vs.nisle > 1:
# calculate island integrals of forcing, keep psi constant on island 1
# ulate island integrals of forcing, keep psi constant on island 1
line_forc[1:] = utilities.line_integrals(vs, fpx[..., np.newaxis],
fpy[..., np.newaxis], kind='same')[1:]
......
veros/setup/zonjet/plot.py
View file @ 86cb7ee8
......@@ -2,6 +2,9 @@ import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
from IPython import get_ipython
import matplotlib
get_ipython().magic('matplotlib inline')
#xr.open_dataset('zonjet.snapshot.nc').u.isel(x,Time=20).plot()
......@@ -48,7 +51,7 @@ class section:
pass
def fileinfo(self):
print('Info: ', self.load().data_vars)
print('Info: ', self.load())
try:
......@@ -56,9 +59,9 @@ try:
except NameError:
pass
p = section('zonjet.snapshot.nc','u')
p = section('zonjet.averages.nc','rho')
p.fileinfo()
p.select(20,'zt',9).plot()
#a.plot(vmax=0.5, vmin=-0.5)
p.select(0,'xt',13).plot.contourf(levels=11)
veros/setup/zonjet/zonjet.py
View file @ 86cb7ee8
#!/usr/bin/env python
from IPython import get_ipython
get_ipython().magic('reset -sf')
from side_tools import *
clean()
from veros import VerosSetup, veros_method
from veros.tools import cli
......@@ -24,11 +29,11 @@ class ZONJETSetup(VerosSetup):
def set_parameter(self, vs):
vs.identifier = 'zonjet'
vs.nx, vs.ny, vs.nz = 15, 30, 10
vs.dt_mom = 4800
vs.nx, vs.ny, vs.nz = 40, 40, 15
vs.dt_mom = 2400
vs.dt_tracer = 86400 / 2.
vs.runlen = 86400 * 365
vs.runlen = 86400 * 51
vs.runlen = 86400 * 365 * 3
#vs.runlen = 86400 * 110
#vs.runlen = 100
vs.coord_degree = True
......@@ -45,14 +50,17 @@ class ZONJETSetup(VerosSetup):
vs.enable_skew_diffusion = True
vs.enable_hor_friction = True
vs.A_h = (2 * vs.degtom)**3 * 2e-11
vs.enable_hor_friction_cos_scaling = True
vs.A_h = (2 * vs.degtom)**3 * 1e-13
vs.enable_hor_friction_cos_scaling = False
vs.hor_friction_cosPower = 1
vs.enable_bottom_friction = False
vs.r_bot = 1e-5
vs.enable_implicit_vert_friction = True
vs.enable_noslip_lateral = True
vs.enable_implicit_vert_friction = False
#vs.enable_tke = True
vs.c_k = 0.1
......@@ -65,7 +73,7 @@ class ZONJETSetup(VerosSetup):
#vs.enable_kappaH_profile = True
# vs.enable_tke_superbee_advection = True
vs.K_gm_0 = 1000.0
vs.K_gm_0 = 100.0
# vs.enable_eke = True
vs.eke_k_max = 1e4
vs.eke_c_k = 0.4
......@@ -73,8 +81,8 @@ class ZONJETSetup(VerosSetup):
vs.eke_cross = 2.
vs.eke_crhin = 1.0
vs.eke_lmin = 100.0
vs.enable_eke_superbee_advection = True
vs.enable_eke_isopycnal_diffusion = True
vs.enable_eke_superbee_advection = False
vs.enable_eke_isopycnal_diffusion = False
vs.enable_idemix = False
......@@ -84,14 +92,13 @@ class ZONJETSetup(VerosSetup):
@veros_method
def set_grid(self, vs):
#ddz = np.array([50., 100., 150., 200., 250., 300., 350., 400.,
# 450., 500., 550., 600., 650., 700., 750.])
vs.dxt[...] = 2.0
vs.dyt[...] = 2.0
ddz = np.array([50., 70., 100., 140., 190., 240., 290., 340.,
390., 440., 490., 540., 590., 640., 690.])
vs.dxt[...] = 1.0
vs.dyt[...] = 1.0
vs.x_origin = 0.0
vs.y_origin = -40.0
#vs.dzt[...] = ddz[::-1] / 2.5
vs.dzt[...] = 20
vs.y_origin = 0.0
vs.dzt[...] = ddz[::-1] / 2.5
@veros_method
def set_coriolis(self, vs):
......@@ -100,51 +107,59 @@ class ZONJETSetup(VerosSetup):
@veros_method
def set_topography(self, vs):
x, y = np.meshgrid(vs.xt, vs.yt, indexing='ij')
vs.kbot[...] = np.logical_or(x > 100., y < 20).astype(np.int)
#vs.kbot[...] = np.logical_or(x > 100., y < 30).astype(np.int)
vs.kbot[...] = np.int(1)
@veros_method
def set_initial_conditions(self, vs):
# initial conditions
# design initial profile for u that is also used for restoring
zscale = 50
yvariance = 5e1
zscale = 500
yvariance = 2e1
ycenter = vs.yt[0] + ( vs.yt[-1] - vs.yt[0] ) / 2
yy,zz = np.meshgrid(vs.yt,vs.zt,indexing='ij')
uz = np.exp(zz * (1/zscale) ) # exponential decay in z
uy = np.exp(- ( (yy-(ycenter) )**2 / yvariance ) ) # gaussian in y
uu = uz * uy
vs.u_ini = allocate(vs, ('xu','yt','zt'))
vs.u_ini = 2 * 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
#plt.contourf(vs.u_ini[:,:,13],10)
#plt.colorbar()
#raise SystemExit(0)
vs.u[...,0:2] = np.tile(vs.u_ini[...,np.newaxis],[1,1,1,2])
#density
vs.temp[:, :, :, 0:2] = ((1 - vs.zt[None, None, :] / vs.zw[0]) * 10 * vs.maskT)[..., None]
vs.salt[:, :, :, 0:2] = 35.0 * vs.maskT[..., None]
rho_background = density.get_rho(vs, vs.salt[...,0], vs.temp[...,0], np.abs(vs.zt)) * vs.maskT
"""
Density profile
"""
zscale=300
rho_back = ( 0.5 - np.exp( vs.zt * 1/zscale )) * 10
#rho_background = density.get_rho(vs, vs.salt[...,0], vs.temp[...,0], np.abs(vs.zt)) * vs.maskT
dzu = np.gradient(vs.u_ini,20)[2] * np.tile(vs.coriolis_t[...,np.newaxis],[1,1,vs.nz]) * (1025-1024) / 10
dzu = np.gradient(vs.u_ini,vs.zt,axis=2) * np.tile(vs.coriolis_t[...,np.newaxis],[1,1,vs.nz]) * 1024 / 10
# cumsum to get density field
#print(vs.yt[-1]-vs.yt[0] * )
sumi = np.cumsum(dzu * 2 * 111000 ,axis=1)
tmp = rho_background[:,2,:].copy()
print(vs.dyt )
sumi = np.cumsum(dzu * vs.dyt[0],axis=1) + rho_back
#tmp = rho_background[:,2,:].copy()
vs.rho_ini = allocate(vs, ('xt', 'yt', 'zt'))
vs.rho_ini = np.tile(tmp[:,np.newaxis,:],[1,vs.ny+4,1]) + sumi
vs.rho_ini = sumi
vs.rho[...,0] = vs.rho_ini[...]
vs.rho[...,1] = vs.rho_ini[...]
vs.rho[...,2] = vs.rho_ini[...]
#import matplotlib.pyplot as plt
#plt.contourf(yy,zz,vs.rho_ini[5,:,:],100)
#plt.colorbar()
#raise SystemExit(0)
vs.temp[...,0] = - 1024 / ( vs.rho_ini*1.67e-4) + 9.85
import matplotlib.pyplot as plt
#plt.plot(vs.zt,rho_back)
plt.contourf(yy,zz,vs.temp[5,:,:,0],10)
plt.colorbar()
raise SystemExit(0)
# wind stress forcing
yt_min = global_min(vs, vs.yt.min())
......@@ -174,16 +189,19 @@ class ZONJETSetup(VerosSetup):
xtmp = np.zeros(vs.xu.shape[0])
xtmp[3:-2] = vs.xu[3:-2]
t_tmp = 1 / (15 * vs.dt_mom) * ( 1-np.tanh((xtmp[-3]-xtmp)/10) )
t_tmp = 1 / (10 * vs.dt_mom) * ( 1-np.tanh((xtmp[-3]-xtmp)/10) )
t_tmp [-2:] = t_tmp[-3]
t_restoring = np.tile((t_tmp[...,np.newaxis,np.newaxis]),[1,vs.ny+4,vs.nz])
t_restoring = 1 / (vs.dt_mom * 15)
#t_restoring = 1 / (vs.dt_mom * 5)
t_restoring_tr = vs.dt_tracer * 1e-6
#import matplotlib.pyplot as plt
#plt.plot(t_tmp)
#plt.pcolormesh(t_restoring[:,:,0]);plt.colorbar()
#raise SystemExit(0)
......@@ -194,8 +212,6 @@ class ZONJETSetup(VerosSetup):
vs.du_rest[...] = vs.maskU * t_restoring * (vs.u_ini - vs.u[...,vs.taum1])
#vs.du[...,vs.taup1] = vs.du[...,vs.taup1]
#+ vs.maskU * t_restoring * (vs.u_ini - vs.u[...,vs.tau])
#vs.rho[...,vs.tau] = vs.rho[...,vs.tau]
......@@ -216,13 +232,13 @@ class ZONJETSetup(VerosSetup):
@veros_method
def set_diagnostics(self, vs):
vs.diagnostics['snapshot'].output_frequency = 86400.
vs.diagnostics['snapshot'].output_frequency = 86400. * 50
#vs.diagnostics['snapshot'].output_frequency = 1.
vs.diagnostics['snapshot'].output_variables = (
'temp', 'psi', 'u', 'v', 'w', 'kbot','rho','salt'
)
#vs.diagnostics['snapshot'].output_variables = (
# 'temp', 'psi', 'u', 'v', 'w', 'kbot','rho','salt'
#)
vs.diagnostics['averages'].output_variables = (
'salt', 'temp', 'u', 'v', 'w', 'psi', 'surface_taux', 'surface_tauy'
'salt', 'temp', 'u', 'v', 'w', 'psi', 'rho', 'surface_tauy'
)
vs.diagnostics['averages'].output_frequency = 50 * 86400.
vs.diagnostics['averages'].sampling_frequency = vs.dt_tracer * 10
......
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