dwbc2 produces reasonable 1deg dwbc

batch.sh

+#!/bin/bash
+#SBATCH -J veros              # Specify job name
+#SBATCH -p prepost             # Use partition prepost
+#SBATCH -N 1                   # Specify number of nodes
+#SBATCH -n 1                  # Specify max. number of tasks to be invoked
+#SBATCH --mem-per-cpu=5300     # Set memory required per allocated CPU
+#SBATCH -t 08:00:00            # Set a limit on the total run time
+#SBATCH -A mh0256              # Charge resources on this project account
+#SBATCH -o my_job.o%j          # File name for standard and error output
+
+# Execute a serial program, e.g.
+source activate veros
+python merjet.py

veros/core/momentum.py

@@ -70,6 +70,14 @@ def momentum(vs):
     vs.du[:, :, :, vs.tau] += vs.du_adv
     vs.dv[:, :, :, vs.tau] += vs.dv_adv
     
+    if vs.dwbc2:
+        print('linear at boundary')
+        # vs.du[:,:1,:,vs.tau] = vs.du[:,:1,:,vs.tau] - vs.du_adv[:,:1,:]
+        vs.du[:,:3,:,vs.tau] -= vs.du_adv[:,:3,:]
+        vs.du[:,-4:,:,vs.tau] -= vs.du_adv[:,-4:,:]
+        vs.dv[:,:3,:,vs.tau] -= vs.dv_adv[:,:3,:]
+        vs.dv[:,-4:,:,vs.tau] -= vs.dv_adv[:,-4:,:]
+    
     """
     add momentum restoring that is computed in set_forcing
     """

veros/core/utilities.py

@@ -17,6 +17,49 @@ def enforce_boundaries(vs, arr, local=False):
         return
 
     exchange_overlap(vs, arr, ['xt', 'yt'])
+    
+@veros_method
+def inflow_outflow(vs, u, v, temp, local=False):
+    from ..distributed import exchange_cyclic_boundaries, exchange_overlap
+    from ..decorators import CONTEXT
+
+    print('enforce inflow outflow')
+    # raise SystemExit(0)
+
+
+    # v[:,:3,:] = vs.v_ini[:,4:7,:].copy()
+    # v[:,-4:,:] = vs.v_ini[:,4:8,:].copy()
+
+    # temp[:,:3,:] = vs.temp_ini[:,4:7,:].copy()
+    temp[:,2,:] = vs.temp_ini[:,5,:].copy()
+
+    temp[:,-3,:] = vs.temp_ini[:,-5,:].copy()
+    # print(vs.v_ini[:,4:7,:].shape,'INI')
+    # print(v[:,-3:,:].shape,'UPPER')
+    # raise SystemExit(0)
+
+    # u[:,:3,:] = 0.
+    # u[:,-4:,:] = 0.
+
+    # arr[:,-1,:] = arr[:,-2,:].copy()
+    # exchange_overlap(vs, arr, ['xt', 'yt'])
+
+@veros_method
+def obc_boundaries(vs,arr):
+    arr[:,0,:] = arr[:,2,:].copy()
+    arr[:,1,:] = arr[:,2,:].copy()
+    arr[:,-1,:] = arr[:,-4,:].copy()
+    arr[:,-2,:] = arr[:,-4,:].copy()
+    arr[:,-3,:] = arr[:,-4,:].copy()
+
+@veros_method
+def set_obc_temp_salt(vs):
+    print('do phase speed stuff')
+
+    vs.temp[:,2,:,vs.taup1] = vs.temp_ini[:,5,:].copy()
+    vs.temp[:,-4,:,vs.taup1] = vs.temp_ini[:,-5,:].copy()
+    print('restore temp and salt')
+
 
 
 @veros_method(inline=True)

veros/settings.py

@@ -37,6 +37,7 @@ SETTINGS = OrderedDict([
     ('enable_store_bottom_friction_tke', Setting(False, bool, 'transfer dissipated energy by bottom/rayleig fric. to TKE, else transfer to internal waves')),
     ('enable_store_cabbeling_heat', Setting(False, bool, 'transfer non-linear mixing terms to potential enthalpy, else transfer to TKE and EKE')),
     ('enable_noslip_lateral', Setting(False, bool, 'enable lateral no-slip boundary conditions in harmonic- and biharmonic friction.')),
+    ('dwbc', Setting(False, bool, 'choose for inflow outflow configuration')),
 
     # External mode
     ('congr_epsilon', Setting(1e-12, float, 'convergence criteria for Poisson solver')),

veros/setup/dwbc1/__init__.py

+from veros.setup.dwbc1.dwbc1 import dwbc1

veros/setup/dwbc1/dwbc.py

+#!/usr/bin/env veros
+#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
+from veros.variables import allocate
+from veros.distributed import global_min, global_max
+from veros.core import utilities
+
+import numpy as np
+
+from veros.core import density
+
+class DWBC1Setup(VerosSetup):
+    """A model using spherical coordinates with a partially closed domain representing the Atlantic and ACC.
+
+    Wind forcing over the channel part and buoyancy relaxation drive a large-scale meridional overturning circulation.
+
+    This setup demonstrates:
+     - setting up an idealized geometry
+     - updating surface forcings
+     - basic usage of diagnostics
+
+    `Adapted from pyOM2 <https://wiki.cen.uni-hamburg.de/ifm/TO/pyOM2/ACC%202>`_.
+    """
+    @veros_method
+    def set_parameter(self, vs):
+        vs.identifier = 'dwbc1'
+        vs.dwbc = True
+        
+        vs.enable_biharmonic_mixing = True
+        vs.K_hbi = 1e12
+        vs.enable_biharmonic_friction = True
+        vs.A_hbi = 1e10
+        
+        
+        # vs.restart_input_filename='merjet_9637.restart.h5'
+
+        vs.nx, vs.ny, vs.nz = 30, 35, 30
+        vs.dt_mom = 300 # normally 2400
+        vs.dt_tracer = 86400 / 16
+        vs.runlen = 86400 * 365 * 2
+        vs.runlen = 86400 * 10
+        # vs.restart_frequency = 86400 * 100
+        # vs.runlen = 100
+
+        vs.coord_degree = True
+
+
+        vs.congr_epsilon = 1e-12
+        vs.congr_max_iterations = 5000
+
+        vs.enable_neutral_diffusion = True
+        vs.K_iso_0 = 1000.0
+        vs.K_iso_steep = 500.0
+        vs.iso_dslope = 0.005
+        vs.iso_slopec = 0.01
+        vs.enable_skew_diffusion = True
+
+        vs.enable_hor_friction = True
+        vs.A_h = (2 * vs.degtom)**3 * 1e-13
+        vs.A_h = 1000
+        vs.K_h = 1000
+        
+        vs.kappaM = 1e-4
+        vs.kappaH = 1e-4
+        # raise SystemExit(0)
+
+        # vs.enable_bottom_friction = True
+        vs.r_bot = 1e-3
+
+        vs.enable_noslip_lateral = True
+
+
+        # vs.enable_implicit_vert_friction = True
+
+
+        vs.K_gm_0 = 1000.0
+ 
+        vs.eq_of_state_type = 1
+        
+        
+        vs.enable_tempsalt_sources = True 
+
+    @veros_method
+    def set_grid(self, vs):
+        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[...] = 1
+        vs.dyt[...] = 1
+        vs.x_origin = 5.0
+        vs.y_origin = 10
+        
+        # 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
+    def set_coriolis(self, vs):        
+        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.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)
+
+    @veros_method
+    def set_topography(self, vs):
+        x, y = np.meshgrid(vs.xt, vs.yt, indexing='ij')        
+        vs.kbot[:,:] = np.int(1)
+        vs.kbot[:5,:] = np.int(0)
+    
+        # import matplotlib.pyplot as plt
+        # plt.contourf(vs.kbot)
+        # plt.colorbar()
+        # raise SystemExit(0)
+        
+    @veros_method
+    def set_initial_conditions(self, vs):
+        
+        """
+        Flow initial and restoring profile
+        """
+        zcenter = 1500
+        zvariance = 3e5
+        xvariance = 2.5
+        xcenter =  9 # normally 16
+        xx,zz = np.meshgrid(vs.xt,vs.zt,indexing='ij')
+        uz = np.exp(- ( (zz+(zcenter) )**2 / zvariance ) ) # exponential decay in z
+        ux = np.exp(- ( (xx-(xcenter) )**2 / xvariance ) ) # gaussian in x
+        
+        uu = uz * ux
+        
+        print(uu.shape)
+        print(vs.yu.shape[0])
+        vs.v_ini = allocate(vs, ('xt','yu','zt'))
+        print(vs.v_ini.shape, 'vs_v_ini')
+        print(zz.shape, 'zz')
+
+        vs.v_ini = -0.2 * np.tile(uu[:,np.newaxis,:],[1, vs.yu.shape[0] , 1])
+        
+        # vs.v_ini[:,:3,:] = 0.
+        # vs.v_ini[:,-4:,:] = 0.
+    
+        
+        # import matplotlib.pyplot as plt
+        # plt.contourf(xx,zz,vs.v_ini[:,1,:])
+        # plt.colorbar()
+        # raise SystemExit(0)
+
+        vs.v[...,0:2] = np.tile(vs.v_ini[...,np.newaxis],[1,1,1,2])
+        
+        """
+        Density profile in thermal wind balance with flow condition
+        """
+        #background temp stratification
+        vs.temp[:, :, :, 0:2] = ((1 - vs.zt[None, None, :] / vs.zw[0]) * 22 * vs.maskT)[..., None]
+        vs.salt[:, :, :, 0:2] = 35.0 * vs.maskT[..., None]
+        
+        # import matplotlib.pyplot as plt
+        # plt.plot(vs.zt,vs.temp[2,2,:,0])
+        # plt.show()
+        # raise SystemExit(0)
+        
+        
+        temp = vs.temp[...,0].copy()
+        salt = vs.salt[...,0].copy()
+        rho_back = density.get_rho(vs, salt, temp, np.abs(vs.zt)) * vs.maskT  
+        
+        # dzu = np.gradient(vs.u_ini,vs.zt,axis=2) * np.tile(vs.coriolis_t[...,np.newaxis],[1,1,vs.nz]) * 1024 / 10
+        dzv = -1 * np.gradient(vs.v_ini,vs.zt,axis=2) * vs.coriolis_t[0,np.int(xcenter)] * 1024 / 10
+        # cumsum to get density field
+        # 
+
+        #print(vs.dyt )
+        sumi = allocate(vs,('xt','yt','zt'))
+        sumi[2:-2,:,:] = np.cumsum(dzv[2:-2,:,:]*100000,axis=0)+ rho_back[2:-2,:,:]
+        
+        vs.rho_ini = allocate(vs, ('xt', 'yt', 'zt'))
+        vs.rho_ini = sumi.copy()
+
+        vs.temp_ini = allocate(vs, ('xt', 'yt', 'zt'))
+        vs.temp_ini[2:-2,:,:] = (- (vs.rho_ini[2:-2,:,:] / 1024.) / (1.67e-4) + 9.85 )
+        
+                
+        vs.temp[...,0:2] = vs.temp_ini[...,None] * vs.maskT[..., None]
+        
+        """
+        Restoring timescales
+        """
+        
+        ytmp = np.zeros(vs.yt.shape[0])
+        ytmp[2:-2] = vs.yt[2:-2]
+
+        # t_tmp =  2-np.tanh((xtmp[-3]-xtmp)/1) + np.tanh((-xtmp)/2)
+        t_tmp = 1-np.tanh((ytmp[-4]-ytmp[:]))
+        t_tmp [-2:] = t_tmp[-3]
+        
+        ft = np.flipud(t_tmp)
+        
+        vs.t_restoring = allocate(vs,('xt','yt','zt'))
+        vs.t_restoring_tr = allocate(vs,('xt','yt','zt'))
+
+        # vs.t_restoring_tr = 2e-5 * np.tile((ft+t_tmp)[np.newaxis,:,np.newaxis],[vs.nx+4,1,vs.nz])
+        
+
+        # vs.t_restoring[-3:,:,:] = 1 / (1e1 * vs.dt_mom)
+        # vs.t_restoring[-5:-3,:,:] = 1 / (2e1 * vs.dt_mom)
+       
+        vs.t_restoring_tr[:,:,-1] = 5e-4
+        vs.t_restoring_tr[:,:,-2] = 2e-4
+        
+        # vs.t_restoring_tr[-3:,:,:] = 2e-5
+        # vs.t_restoring_tr[-5:-3,:,:] = 1e-5
+
+        # import matplotlib.pyplot as plt
+        # plt.plot(t_tmp+ft)
+        # xx,yy = np.meshgrid(vs.xt,vs.yt,indexing='ij')
+
+        # plt.contourf(xx,yy,vs.t_restoring_tr[:,:,5])
+        # plt.contourf(vs.temp_ini[:,-3,:],10)
+        # plt.contour(vs.v_ini[:,10,:],5,colors='black')
+
+        # plt.colorbar()
+        # raise SystemExit(0)
+        
+        """
+        Boundaries for OBC
+        """
+        if vs.dwbc:
+           utilities.obc_boundaries(vs,vs.temp_ini)
+           utilities.obc_boundaries(vs,vs.temp[:,:,:,0:2])
+           utilities.obc_boundaries(vs,vs.v_ini)
+           utilities.obc_boundaries(vs,vs.v[:,:,:,0:2])
+
+    @veros_method
+    def set_forcing(self, vs):
+        # vs.dv_rest[...] = vs.maskV * vs.t_restoring * (vs.v_ini[...] - vs.v[...,vs.tau])
+
+        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
+    def set_diagnostics(self, vs):
+        vs.diagnostics['snapshot'].output_frequency = 1
+        #vs.diagnostics['snapshot'].output_frequency = 1.        
+        vs.diagnostics['snapshot'].output_variables = (
+            'temp', 'psi', 'u', 'v', 'w', 'coriolis_t','rho','salt'
+        )    
+        vs.diagnostics['averages'].output_variables = (
+            'salt', 'temp', 'u', 'v', 'w', 'psi', 'rho', 'du_rest'
+        )
+        vs.diagnostics['averages'].output_frequency = 20. * 86400.
+        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 = 10 * 86400.
+        # vs.diagnostics['energy'].output_frequency = 365 * 86400.
+        # vs.diagnostics['energy'].sampling_frequency = vs.dt_tracer * 10
+        #vs.diagnostics['snapshot'].output_variables = (
+        #        'kbot')
+
+    def after_timestep(self, vs):
+        pass
+
+
+@cli
+def run(*args, **kwargs):
+    simulation = DWBC1Setup(*args, **kwargs)              
+    simulation.setup()
+    simulation.run()
+
+
+if __name__ == '__main__':
+    # Execute only if run as a script
+    run()
+    
+    
+    

veros/setup/dwbc1/plot.py

+#%%
+import xarray as xr
+import matplotlib.pyplot as plt
+import numpy as np
+
+from IPython import get_ipython
+get_ipython().magic('reset -sf')
+
+import matplotlib
+get_ipython().magic('matplotlib inline')
+
+
+class section:
+    def __init__(self, filename, varname):
+        self.filename = filename
+        self.varname = varname
+        try:
+            self.var = self.load()[varname]
+        except KeyError:
+            print('Variable does not exist, choose new from:')
+            print(self.load().data_vars)
+            
+        self.data = 'something'
+        
+    def load(self):
+        return xr.open_dataset(self.filename)
+    
+    
+    def select(self, timestep, axis, position):
+        try: 
+            self.var[axis]
+            if axis == 'xu':
+                self.data = self.var.isel(xu=position).isel(Time=timestep)
+            elif axis == 'yt':
+                self.data = self.var.isel(yt=position).isel(Time=timestep)
+            elif axis == 'zt':
+                self.data = self.var.isel(zt=position).isel(Time=timestep)
+            elif axis == 'xt':
+                self.data = self.var.isel(xt=position).isel(Time=timestep)
+        except KeyError:
+            print('Axis does not exist, choose new from: \n', self.var.coords)
+            
+
+
+        
+        return self.data
+        
+    def figure(self):
+        pass
+    
+    def fileinfo(self):
+        print('Info: ', self.load())
+
+
+try:
+    del p
+except NameError:
+    pass
+
+p = section('dwbc1.snapshot.nc','temp')
+
+
+plt.figure()
+p.select(28,'zt',15).plot.contourf()
+# 
+
+# plt.colorbar()
+# plt.show()
+
+

veros/setup/dwbc2/__init__.py

+from veros.setup.dwbc1.dwbc1 import dwbc1

veros/setup/dwbc2/dwbc2.py

+#!/usr/bin/env veros
+#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
+from veros.variables import allocate
+from veros.distributed import global_min, global_max
+from veros.core import utilities
+
+import numpy as np
+
+from veros.core import density
+
+class DWBC2Setup(VerosSetup):
+    """A model using spherical coordinates with a partially closed domain representing the Atlantic and ACC.
+
+    Wind forcing over the channel part and buoyancy relaxation drive a large-scale meridional overturning circulation.
+
+    This setup demonstrates:
+     - setting up an idealized geometry
+     - updating surface forcings
+     - basic usage of diagnostics
+
+    `Adapted from pyOM2 <https://wiki.cen.uni-hamburg.de/ifm/TO/pyOM2/ACC%202>`_.
+    """
+    @veros_method
+    def set_parameter(self, vs):
+        
+        """
+        GENERAL
+        """
+        vs.identifier = 'dwbc2'
+        vs.dwbc2 = True
+        vs.eq_of_state_type = 1
+        vs.enable_tempsalt_sources = True 
+        
+        # vs.restart_input_filename='merjet_9637.restart.h5'
+        
+        
+        """
+        GRID, TIMESTEP,...
+        """
+        vs.nx, vs.ny, vs.nz = 20, 35, 30
+        vs.dt_mom = 1200 # normally 2400
+        vs.dt_tracer = 86400 / 4
+        vs.runlen = 86400 * 365 * 1
+        # vs.runlen = 86400 * 100
+        # vs.restart_frequency = 86400 * 100
+        # vs.runlen = 100
+
+        vs.coord_degree = True
+        vs.congr_epsilon = 1e-12
+        vs.congr_max_iterations = 5000
+        
+        
+        """
+        Friction, Diffusion, ...
+        """
+        
+        'BOUNDARIES'
+        vs.enable_noslip_lateral = True
+        
+        ' BIHARMONIC'
+        vs.enable_biharmonic_mixing = True
+        vs.K_hbi = 1e10
+        vs.enable_biharmonic_friction = True
+        vs.A_hbi = 1e10
+        
+        'HARMONIC'
+        vs.enable_hor_friction = True
+        vs.A_h = 1000
+        vs.K_h = 1000
+        
+        'VERTICAL'
+        vs.enable_implicit_vert_friction = True
+        vs.kappaM = 1e-4
+        vs.kappaH = 1e-4
+    
+        'GM /REDI'
+        vs.enable_neutral_diffusion = False
+        vs.K_iso_0 = 1000.0
+        vs.K_iso_steep = 500.0
+        vs.iso_dslope = 0.005
+        vs.iso_slopec = 0.01
+        vs.enable_skew_diffusion = False
+        vs.K_gm_0 = 1000.0
+
+        'BOTTOM'
+        vs.enable_bottom_friction = True
+        vs.r_bot = 1e-3
+
+
+
+
+    @veros_method
+    def set_grid(self, vs):
+        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[...] = 1
+        vs.dyt[...] = 1
+        vs.x_origin = 5.0
+        vs.y_origin = 10
+        
+        # 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
+    def set_coriolis(self, vs):        
+        vs.coriolis_t[:, :] = 2 * vs.omega * np.sin(vs.yt[None, :] / 180. * vs.pi)
+
+    @veros_method
+    def set_topography(self, vs):
+        x, y = np.meshgrid(vs.xt, vs.yt, indexing='ij')        
+        vs.kbot[:,:] = np.int(1)
+        vs.kbot[:3,:] = np.int(0)
+    
+        # import matplotlib.pyplot as plt
+        # plt.contourf(vs.kbot)
+        # plt.colorbar()
+        # raise SystemExit(0)
+        
+    @veros_method
+    def set_initial_conditions(self, vs):
+        
+        """
+        Flow initial and restoring profile
+        """
+        zcenter = 1500
+        zvariance = 3e5
+        xvariance = 2.
+        xcenter =  7.5 # normally 16
+        xx,zz = np.meshgrid(vs.xt,vs.zt,indexing='ij')
+        uz = np.exp(- ( (zz+(zcenter) )**2 / zvariance ) ) # exponential decay in z
+        ux = np.exp(- ( (xx-(xcenter) )**2 / xvariance ) ) # gaussian in x
+        uu = uz * ux
+        
+        vs.v_ini = allocate(vs, ('xt','yu','zt'))
+        vs.v_ini = -0.2 * np.tile(uu[:,np.newaxis,:],[1, vs.yu.shape[0] , 1]) * vs.maskV
+
+        # vs.v[...,0:2] = np.tile(vs.v_ini[...,np.newaxis],[1,1,1,2])
+        
+        # vs.v[:,:2,:,0:2] = 0.
+        # vs.v[:,-3:,:,0:2] = 0.
+        
+        
+        """
+        Density profile in thermal wind balance with flow condition
+        """
+        #background temp stratification
+        vs.temp[:, :, :, 0:2] = ((1 - vs.zt[None, None, :] / vs.zw[0]) * 22 * vs.maskT)[..., None]
+        vs.salt[:, :, :, 0:2] = 35.0 * vs.maskT[..., None]
+        
+        vs.temp_surface = allocate(vs,('xt','yt','zt'))
+        vs.temp_surface = vs.temp[:,:,:,0].copy()
+        
+        # import matplotlib.pyplot as plt
+        # plt.pcolormesh(xx,zz,vs.v[:,5,:,0])
+        # plt.colorbar()
+        # plt.show()
+        # raise SystemExit(0)
+        
+        
+        temp = vs.temp[...,0].copy()
+        salt = vs.salt[...,0].copy()
+        rho_back = density.get_rho(vs, salt, temp, np.abs(vs.zt)) * vs.maskT  
+        
+        # dzu = np.gradient(vs.u_ini,vs.zt,axis=2) * np.tile(vs.coriolis_t[...,np.newaxis],[1,1,vs.nz]) * 1024 / 10
+        dzv = -1 * np.gradient(vs.v_ini,vs.zt,axis=2) * vs.coriolis_t[0,np.int(xcenter)] * 1024 / 10 * vs.maskV
+        # cumsum to get density field
+        # 
+
+        #print(vs.dyt )
+        sumi = allocate(vs,('xt','yt','zt'))
+        sumi[:,:,:] = np.cumsum(dzv[:,:,:]*vs.dxt[0],axis=0)+ rho_back[:,:,:]
+        
+        vs.rho_ini = allocate(vs, ('xt', 'yt', 'zt'))
+        vs.rho_ini = sumi.copy()
+
+        vs.temp_ini = allocate(vs, ('xt', 'yt', 'zt'))
+        vs.temp_ini[:,:,:] = (- (vs.rho_ini[:,:,:] / 1024.) / (1.67e-4) + 9.85 )
+        
+                
+        vs.temp[...,0:2] = vs.temp_ini[...,None] * vs.maskT[..., None]
+        
+        tmp = np.gradient(vs.temp_ini[...],axis=0)
+        """
+        Boundaries for OBC
+        """
+        if vs.dwbc2:
+           utilities.obc_boundaries(vs,vs.temp_ini)
+           utilities.obc_boundaries(vs,vs.temp[:,:,:,0:2])
+           utilities.obc_boundaries(vs,vs.v_ini)
+           utilities.obc_boundaries(vs,vs.v[:,:,:,0:2])
+           utilities.obc_boundaries(vs,tmp)
+           utilities.obc_boundaries(vs,dzv)
+
+        # import matplotlib.pyplot as plt
+        # plt.plot(t_tmp+ft)
+        # xx,yy = np.meshgrid(vs.xt,vs.yt,indexing='ij')
+
+        # plt.contourf(xx,yy,vs.temp_ini[:,:,15],20)
+        # plt.contourf(xx,yy,tmp[:,:,10])
+        # plt.colorbar()
+        # plt.contour(xx,zz,vs.v_ini[:,10,:],5,colors='black')
+
+        # raise SystemExit(0)
+        
+
+           
+        """
+        Restoring timescales
+        """
+        
+        ytmp = np.zeros(vs.yt.shape[0])
+        ytmp[2:-2] = vs.yt[2:-2]
+
+        # t_tmp =  2-np.tanh((xtmp[-3]-xtmp)/1) + np.tanh((-xtmp)/2)
+        t_tmp = 1-np.tanh((ytmp[-4]-ytmp[:]))
+        t_tmp [-2:] = t_tmp[-3]
+        
+        ft = np.flipud(t_tmp)
+        
+        vs.t_restoring = allocate(vs,('xt','yt','zt'))
+        vs.t_restoring_tr = allocate(vs,('xt','yt','zt'))
+
+        # vs.t_restoring_tr = 2e-5 * np.tile((ft+t_tmp)[np.newaxis,:,np.newaxis],[vs.nx+4,1,vs.nz])
+        
+
+        # vs.t_restoring[-3:,:,:] = 1 / (1e1 * vs.dt_mom)
+        # vs.t_restoring[-5:-3,:,:] = 1 / (2e1 * vs.dt_mom)
+       
+        vs.t_restoring_tr[:,:,-1] = 2e-5
+        vs.t_restoring_tr[:,:,-2] = 1e-5
+
+    @veros_method
+    def set_forcing(self, vs):
+        pass
+        # vs.dv_rest[...] = vs.maskV * vs.t_restoring * (vs.v_ini[...] - vs.v[...,vs.tau])
+
+        # vs.temp_source[...] = vs.maskT * vs.t_restoring_tr * (vs.temp_surface - vs.temp[... , vs.tau])
+        # vs.salt_source[...] = vs.maskT * vs.t_restoring_tr * (35 - vs.salt[... , vs.tau])
+