found new configuration with meridional flow

veros/core/momentum.py

@@ -73,8 +73,6 @@ def momentum(vs):
     """
     add momentum restoring that is computed in set_forcing
     """
-    vs.du[:, :, :, vs.tau] += vs.du_rest
-    vs.dv[:, :, :, vs.tau] += vs.dv_rest
 
     with vs.timers['friction']:
         """
@@ -120,6 +118,9 @@ def momentum(vs):
     """
     external mode
     """
+    vs.du_mix += vs.du_rest
+    vs.dv_mix += vs.dv_rest
+    
     with vs.timers['pressure']:
         streamfunction.solve_streamfunction(vs)
 

veros/settings.py

@@ -17,6 +17,7 @@ SETTINGS = OrderedDict([
     # Logical switches for general model setup
     ('coord_degree', Setting(False, bool, 'either spherical (True) or cartesian (False) coordinates')),
     ('enable_cyclic_x', Setting(False, bool, 'enable cyclic boundary conditions')),
+    ('enable_rotate_f', Setting(False, bool, 'change coriolis gradient to zonal')),
     ('eq_of_state_type', Setting(1, int, 'equation of state: 1: linear, 3: nonlinear with comp., 5: TEOS')),
     ('enable_implicit_vert_friction', Setting(False, bool, 'enable implicit vertical friction')),
     ('enable_explicit_vert_friction', Setting(False, bool, 'enable explicit vertical friction')),

veros/setup/acc/acc.py

@@ -136,7 +136,7 @@ class ACCSetup(VerosSetup):
         vs.diagnostics['averages'].output_variables = (
             'salt', 'temp', 'u', 'v', 'w', 'psi', 'surface_taux', 'surface_tauy','rho'
         )
-        vs.diagnostics['averages'].output_frequency = 365 * 86400.
+        vs.diagnostics['averages'].output_frequency = 50 * 86400.
         vs.diagnostics['averages'].sampling_frequency = vs.dt_tracer * 10
         vs.diagnostics['overturning'].output_frequency = 365 * 86400. / 48.
         vs.diagnostics['overturning'].sampling_frequency = vs.dt_tracer * 10

veros/setup/merjet/__init__.py

+from veros.setup.merjet.merjet import merjet

veros/setup/merjet/merjet.py

+#!/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
+from veros.variables import allocate
+from veros.distributed import global_min, global_max
+import numpy as np
+
+from veros.core import density
+
+class MERJETSetup(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 = '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.runlen = 100
+
+        vs.coord_degree = True
+        vs.enable_cyclic_x = True
+        vs.enable_rotate_f = True
+
+
+        vs.congr_epsilon = 1e-12
+        vs.congr_max_iterations = 5000
+
+        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.enable_hor_friction = True
+        vs.A_h = (2 * vs.degtom)**3 * 1e-13
+        vs.A_h = 10
+        # print(vs.A_h)
+        # raise SystemExit(0)
+
+        vs.enable_bottom_friction = False
+        vs.r_bot = 1e-4
+
+        vs.enable_noslip_lateral = True
+
+
+        # vs.enable_implicit_vert_friction = True
+
+        #vs.enable_tke = True
+        vs.c_k = 0.1
+        vs.c_eps = 0.7
+        vs.alpha_tke = 30.0
+        vs.mxl_min = 1e-8
+        vs.tke_mxl_choice = 2
+        vs.kappaM_min = 1e-4
+        vs.kappaH_min = 1e-4
+        vs.enable_kappaH_profile = False
+        # vs.enable_tke_superbee_advection = True
+
+        vs.K_gm_0 = 1000.0
+       # vs.enable_eke = True
+        vs.eke_k_max = 1e4
+        vs.eke_c_k = 0.4
+        vs.eke_c_eps = 0.5
+        vs.eke_cross = 2.
+        vs.eke_crhin = 1.0
+        vs.eke_lmin = 100.0
+        vs.enable_eke_superbee_advection = False
+        vs.enable_eke_isopycnal_diffusion = False
+
+        vs.enable_idemix = False
+
+        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.x_origin = 5.0
+        vs.y_origin = 10
+        # vs.dzt[...] = ddz[::-1] / 2.5
+        vs.dzt[:] = 60
+    
+
+    @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.1
+            xtmp[0] = xtmp[2]- 2 * 0.1
+            xtmp[-1] = xtmp[-3] + 2 * 0.1
+            xtmp[-2] = xtmp[-3] + 0.1
+
+            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()
+                
+        # 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)
+        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)
+        # 
+        # import matplotlib.pyplot as plt
+        # plt.contourf(vs.kbot)
+        # plt.colorbar()
+        # raise SystemExit(0)
+        
+    @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
+        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])
+    
+        
+        # import matplotlib.pyplot as plt
+        # plt.contourf(yy,zz,vs.u_ini[0,:,:],10)
+        # plt.colorbar()
+        # raise SystemExit(0)
+
+        vs.u[...,0:2] = np.tile(vs.u_ini[...,np.newaxis],[1,1,1,2])
+        
+        """
+        Density profile
+        """
+        #background temp stratification
+        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  
+        
+        # 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,:]
+        
+        vs.rho_ini = allocate(vs, ('xt', 'yt', 'zt'))
+        vs.rho_ini = sumi.copy()
+
+        vs.temp_ini = vs.temp[:,:,:,0].copy()
+        
+                
+        # vs.temp[...,0:2] = vs.temp_ini[...,None] * vs.maskT[..., None]
+        
+        # vs.salt[:, :, :, 0:2] = 35.0 * vs.maskT[..., None]
+        
+        # rho = density.get_rho(vs, 35, vs.temp[...,1], np.abs(vs.zt))
+        
+        # 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.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.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)
+        
+
+    @veros_method
+    def set_diagnostics(self, vs):
+        vs.diagnostics['snapshot'].output_frequency = 5. * 86400.
+        #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 = 365 * 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 = MERJETSetup(*args, **kwargs)              
+    simulation.setup()
+    simulation.run()
+
+
+if __name__ == '__main__':
+    # Execute only if run as a script
+    run()
+    
+    
+    

veros/setup/merjet/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('../acc/acc.averages.nc','u')
+# p = section('good/20200401_1.averages.nc', 'u')
+    
+
+p = section('merjet.snapshot.nc', 'u')
+
+#%%
+p.select(2,'zt',30).plot.contourf(levels=15)
+
+#%%
+
+p = section('merjet.averages.nc', 'u')
+step = 3
+xpos = 5
+
+plt.figure()
+p.select(step,'xu',xpos).plot.contour(levels=15,vmin=-0.6,vmax=0.6,colors="black")
+# 
+p = section('merjet.averages.nc', 'rho')
+p.select(step,'xt',xpos).plot.contourf(levels=16,cmap='RdBu_r')
+
+"""
+Analyze transports
+"""
+p = section('merjet.snapshot.nc', 'u')
+
+trans = np.sum(p.var.isel(xu=1),axis=(1,2))
+plt.figure()
+plt.plot(p.var.Time,trans)
+plt.figure()
+trans = np.sum(p.var,axis=(1,2))
+plt.contourf(trans)
+plt.colorbar()
+# plt.show()
+
+

veros/setup/north_atlantic/north_atlantic.py

@@ -242,7 +242,7 @@ class NorthAtlanticSetup(VerosSetup):
         vs.diagnostics['averages'].output_frequency = 3600. * 24 * 10
         vs.diagnostics['averages'].sampling_frequency = vs.dt_tracer
         vs.diagnostics['averages'].output_variables = ['temp', 'salt', 'u', 'v', 'w',
-                                                       'surface_taux', 'surface_tauy', 'psi']
+                                                       'surface_taux', 'surface_tauy', 'rest_tscl']
         vs.diagnostics['cfl_monitor'].output_frequency = vs.dt_tracer * 10
 
     @veros_method

veros/setup/zonjet/plot.py

@@ -43,7 +43,7 @@ class section:
             print('Axis does not exist, choose new from: \n', self.var.coords)
             
 
-            
+
         
         return self.data
         
@@ -59,9 +59,10 @@ try:
 except NameError:
     pass
 
-p = section('zonjet.averages.nc','rho')
+# p = section('../acc/acc.averages.nc','u')
+p = section('zonjet.averages.nc', 'temp')
 
 p.fileinfo()
 
-p.select(0,'xt',13).plot.contourf(levels=11)
+p.select(6,'zt',14).plot.contourf(levels=100)
 

veros/setup/zonjet/zonjet.py

@@ -29,11 +29,11 @@ class ZONJETSetup(VerosSetup):
     def set_parameter(self, vs):
         vs.identifier = 'zonjet'
 
-        vs.nx, vs.ny, vs.nz = 40, 40, 15
-        vs.dt_mom = 2400
+        vs.nx, vs.ny, vs.nz = 60, 30, 15
+        vs.dt_mom = 600 # normally 2400
         vs.dt_tracer = 86400 / 2.
-        vs.runlen = 86400 * 365 * 3
-        #vs.runlen = 86400 * 110
+        vs.runlen = 86400 * 365 * 1
+        # vs.runlen = 86400 * 160
         #vs.runlen = 100
 
         vs.coord_degree = True
@@ -42,7 +42,7 @@ class ZONJETSetup(VerosSetup):
         vs.congr_epsilon = 1e-12
         vs.congr_max_iterations = 5000
 
-        vs.enable_neutral_diffusion = True
+        vs.enable_neutral_diffusion = False
         vs.K_iso_0 = 1000.0
         vs.K_iso_steep = 500.0
         vs.iso_dslope = 0.005
@@ -51,16 +51,17 @@ class ZONJETSetup(VerosSetup):
 
         vs.enable_hor_friction = True
         vs.A_h = (2 * vs.degtom)**3 * 1e-13
-        vs.enable_hor_friction_cos_scaling = False
-        vs.hor_friction_cosPower = 1
+        vs.A_h = 1e4
+        # print(vs.A_h)
+        # raise SystemExit(0)
 
-        vs.enable_bottom_friction = False
-        vs.r_bot = 1e-5
+        vs.enable_bottom_friction = True
+        vs.r_bot = 1e-4
 
-        vs.enable_noslip_lateral = True
+        # vs.enable_noslip_lateral = True
 
 
-        vs.enable_implicit_vert_friction = False
+        # vs.enable_implicit_vert_friction = True
 
         #vs.enable_tke = True
         vs.c_k = 0.1
@@ -68,12 +69,12 @@ class ZONJETSetup(VerosSetup):
         vs.alpha_tke = 30.0
         vs.mxl_min = 1e-8
         vs.tke_mxl_choice = 2
-        vs.kappaM_min = 2e-4
+        vs.kappaM_min = 1e-5
         vs.kappaH_min = 2e-5
-        #vs.enable_kappaH_profile = True
+        vs.enable_kappaH_profile = False
         # vs.enable_tke_superbee_advection = True
 
-        vs.K_gm_0 = 100.0
+        vs.K_gm_0 = 1000.0
        # vs.enable_eke = True
         vs.eke_k_max = 1e4
         vs.eke_c_k = 0.4
@@ -86,9 +87,11 @@ class ZONJETSetup(VerosSetup):
 
         vs.enable_idemix = False
 
-        vs.eq_of_state_type = 3
+        vs.eq_of_state_type = 1
         
-        vs.force_overwrite = 'True' # added for easier configuring the setting
+        vs.force_overwrite = True # added for easier configuring the setting
+        
+        vs.enable_tempsalt_sources = True 
 
     @veros_method
     def set_grid(self, vs):
@@ -97,13 +100,14 @@ class ZONJETSetup(VerosSetup):
         vs.dxt[...] = 1.0
         vs.dyt[...] = 1.0
         vs.x_origin = 0.0
-        vs.y_origin = 0.0
+        vs.y_origin = 10.
         vs.dzt[...] = ddz[::-1] / 2.5
+        # vs.dzt[:] = 20
+    
 
     @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')
@@ -112,54 +116,73 @@ class ZONJETSetup(VerosSetup):
         
     @veros_method
     def set_initial_conditions(self, vs):
+        
+
         # initial conditions
         
         # design initial profile for u that is also used for restoring
-        zscale = 500
+        zscale = 400
         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
+        # 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 = 0.5 * np.tile(uu[np.newaxis,:,:],[vs.xu.shape[0] , 1 , 1])
+        vs.u_ini = 1 * 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)
+        # import matplotlib.pyplot as plt