smtnatl_dissipation.ipynb: Dissipation in ICON SMT-NATL#
[1]:
# Jupyter Notebook with widget matplotlib plots
#%matplotlib notebook
# Jupyter Lab with widget matplotlib plots
#%matplotlib widget
# with Jupyter and Jupyter Lab but without widget matplotlib plots
%matplotlib inline
%load_ext autoreload
%autoreload 2
[2]:
# %%javascript
# IPython.notebook.kernel.execute('nb_name = "' + IPython.notebook.notebook_name + '"')
[3]:
nb_name = 'smtnatl_dissipation'
[4]:
import sys
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from netCDF4 import Dataset, num2date
from ipdb import set_trace as mybreak
import pyicon as pyic
import cartopy.crs as ccrs
import glob
import pickle
import maps_icon_smt_temp as smt
import datetime
from matplotlib.patches import Rectangle
from icon_smt_levels import dzw, dzt, depthc, depthi
-----calc
sys glob os
numpy
netcdf
Done modules calc.
-----calc_xr
sys glob os
numpy
netcdf
xarray
Done modules calc.
-----tb
sys
json
numpy
scipy
netcdf datetime
matplotlib
mybreak
pnadas
xarray
done xarray
-----IconData
-----plotting
-----accessor
-----simulation
-----view
-----calc
-----calc_xr
-----tb
-----IconData
-----plotting
-----accessor
-----simulation
-----view
-----quickplots
-----quickplots
[5]:
import xarray as xr
import dask
import dask.array as da
from dask.distributed import Client, LocalCluster
import multiprocessing
[6]:
def savefig(txt, dpi=150, format='pdf'):
if False:
fbase = nb_name
fpath = f'{path_fig}{fbase}_{txt}.{format}'
plt.savefig(fpath, dpi=dpi)
print('saved figure: %s' % (fpath))
return
Define simulation#
[7]:
run = 'ngSMT_tke'
savefig = False
path_fig = '../pics/'
nnf=0
gname = 'smt'
lev = 'L128'
# path_data = f'/work/mh0287/users/leonidas/icon/ngSMT/results/????-??/'
path_data = '/work/mh0287/from_Mistral/mh0287/users/leonidas/icon/ngSMT/results/????-??/'
#fpath_tgrid = '/mnt/lustre01/work/mh0287/users/leonidas/icon/submeso/grid/cell_grid-OceanOnly_SubmesoNA_2500m_srtm30-icon.nc'
fpath_tgrid = '/home/m/m300602/work/icon/grids/smt/smt_tgrid.nc'
fpath_Tri = '/work/mh0033/m300602/tmp/Tri.pkl'
path_grid = f'/work/mh0033/m300602/icon/grids/{gname}/'
path_ckdtree = f'{path_grid}ckdtree/'
# fpath_ckdtree = f'{path_grid}ckdtree/rectgrids/{gname}_res0.30_180W-180E_90S-90N.npz'
fpath_ckdtree = '/work/mh0033/m300602/icon/grids/smt/ckdtree/rectgrids/smt_res0.02_180W-180E_90S-90N.npz'
[8]:
# time0 = np.datetime64('2010-02-23T19:00:00')
time0 = np.datetime64('2010-03-09T01:00:00')
lon_reg_1 = [-88, 0]
lat_reg_1 = [ 5, 68]
lon_reg_2 = [-75, -55]
# lat_reg_2 = [30, 40]
lat_reg_2 = [33, 43]
# lon_reg_3 = [-72, -68]
# lat_reg_3 = [34, 36]
# lon_reg_3 = [-69.5, -65.5]
# lat_reg_3 = [36, 38]
lon_reg_3 = [-68.5, -66.5]
lat_reg_3 = [36.5, 37.5]
lon_reg_4 = [-66., -64.]
lat_reg_4 = [41, 42]
depth0 = 50.
iz = (depthc<depth0).sum()
iz, depthc[iz]
[8]:
(16, 51.5)
Load the grid#
[9]:
%%time
IcD = pyic.IconData(
fname = f'{run}_vort_f_50m_*.nc',
path_data = '/work/mh0287/from_Mistral/mh0287/users/leonidas/icon/ngSMT/results/2010-03/',
path_grid = path_grid,
gname = gname,
lev = lev,
do_triangulation = False,
omit_last_file = True,
calc_coeff = False,
load_vertical_grid = False,
load_triangular_grid = False,
load_rectangular_grid = False,
load_variable_info = False,
verbose = True,
time_mode ='float2date',
)
-v-: set paths and fnames
-v-: set global variables
-v-: find ckdtrees etc.
::: Warning: Could not find any section-npz-file in /work/mh0033/m300602/icon/grids/smt/ckdtree/sections/. :::
::: Warning: no section found.:::
-v-: list of variables and time steps
CPU times: user 64.2 ms, sys: 416 ms, total: 480 ms
Wall time: 4.27 s
[10]:
%%time
f = Dataset(fpath_tgrid, 'r')
clon = f.variables['clon'][:] * 180./np.pi
clat = f.variables['clat'][:] * 180./np.pi
vlon = f.variables['vlon'][:] * 180./np.pi
vlat = f.variables['vlat'][:] * 180./np.pi
vertex_of_cell = f.variables['vertex_of_cell'][:].transpose()-1
edge_of_cell = f.variables['edge_of_cell'][:].transpose()-1
edges_of_vertex = f.variables['edges_of_vertex'][:].transpose()-1
dual_edge_length = f.variables['dual_edge_length'][:]
edge_orientation = f.variables['edge_orientation'][:].transpose()
edge_length = f.variables['edge_length'][:]
cell_area_p = f.variables['cell_area_p'][:]
dual_area = f.variables['dual_area'][:]
edge_length = f.variables['edge_length'][:]
adjacent_cell_of_edge = f.variables['adjacent_cell_of_edge'][:].transpose()-1
orientation_of_normal = f.variables['orientation_of_normal'][:].transpose()
edge_vertices = f.variables['edge_vertices'][:].transpose()-1
f.close()
CPU times: user 1.94 s, sys: 7.77 s, total: 9.71 s
Wall time: 13.2 s
[11]:
res = np.sqrt(cell_area_p)
[12]:
%%time
dtype = 'float32'
f = Dataset(fpath_tgrid, 'r')
elon = f.variables['elon'][:] * 180./np.pi
elat = f.variables['elat'][:] * 180./np.pi
cell_cart_vec = np.ma.zeros((clon.size,3), dtype=dtype)
cell_cart_vec[:,0] = f.variables['cell_circumcenter_cartesian_x'][:]
cell_cart_vec[:,1] = f.variables['cell_circumcenter_cartesian_y'][:]
cell_cart_vec[:,2] = f.variables['cell_circumcenter_cartesian_z'][:]
# vert_cart_vec = np.ma.zeros((vlon.size,3), dtype=dtype)
# vert_cart_vec[:,0] = f.variables['cartesian_x_vertices'][:]
# vert_cart_vec[:,1] = f.variables['cartesian_y_vertices'][:]
# vert_cart_vec[:,2] = f.variables['cartesian_z_vertices'][:]
edge_cart_vec = np.ma.zeros((elon.size,3), dtype=dtype)
edge_cart_vec[:,0] = f.variables['edge_middle_cartesian_x'][:]
edge_cart_vec[:,1] = f.variables['edge_middle_cartesian_y'][:]
edge_cart_vec[:,2] = f.variables['edge_middle_cartesian_z'][:]
# edge_prim_norm = np.ma.zeros((elon.size,3), dtype=dtype)
# edge_prim_norm[:,0] = f.variables['edge_primal_normal_cartesian_x'][:]
# edge_prim_norm[:,1] = f.variables['edge_primal_normal_cartesian_y'][:]
# edge_prim_norm[:,2] = f.variables['edge_primal_normal_cartesian_z'][:]
f.close()
CPU times: user 1.74 s, sys: 4.49 s, total: 6.23 s
Wall time: 7.51 s
[13]:
%%time
ds_tg = xr.open_mfdataset(fpath_tgrid)
res_xr = np.sqrt(ds_tg.cell_area_p).compute()
CPU times: user 184 ms, sys: 696 ms, total: 880 ms
Wall time: 1.22 s
[14]:
%%time
lon_reg = lon_reg_3
lat_reg = lat_reg_3
clon_reg, clat_reg, vertex_of_cell_reg, edge_of_cell_reg, ind_reg = pyic.crop_tripolar_grid(lon_reg, lat_reg,
clon, clat, vertex_of_cell, edge_of_cell)
Tri = matplotlib.tri.Triangulation(vlon, vlat, triangles=vertex_of_cell_reg)
CPU times: user 356 ms, sys: 288 ms, total: 644 ms
Wall time: 1.06 s
Plot resolution#
[15]:
fpath_ckdtree_res = '/work/mh0033/m300602/icon/grids/smt/ckdtree/rectgrids/smt_res0.10_180W-180E_90S-90N.npz'
ddnpz = np.load(fpath_ckdtree_res)
lon_010deg = ddnpz['lon']
lat_010deg = ddnpz['lat']
res_010deg = pyic.apply_ckdtree(res, fpath_ckdtree_res, coordinates='clat clon').reshape(lat_010deg.size, lon_010deg.size)
res_010deg[res_010deg==0.] = np.ma.masked
[16]:
res_xr
[16]:
<xarray.DataArray 'cell_area_p' (cell: 59799625)>
array([ 567.72810763, 567.67196579, 567.67694549, ..., 1221.08765225,
1220.44448307, 1220.3412435 ])
Coordinates:
clon (cell) float64 -0.9462 -0.9462 -0.9463 ... -1.253 -1.253 -1.253
clat (cell) float64 0.6369 0.637 0.6369 0.6369 ... 0.1647 0.1649 0.165
Dimensions without coordinates: cell
Attributes:
long_name: area of grid cell
units: m2
grid_type: unstructured
number_of_grid_in_reference: 1[17]:
res_xr.rename(cell='ncells').pyic.plot()
/home/m/m300602/miniconda3/envs/pyicon_py39_cartopy19/lib/python3.9/site-packages/cartopy/crs.py:228: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry.
if len(multi_line_string) > 1:
/home/m/m300602/miniconda3/envs/pyicon_py39_cartopy19/lib/python3.9/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry.
line_strings = list(multi_line_string)
/home/m/m300602/miniconda3/envs/pyicon_py39_cartopy19/lib/python3.9/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
line_strings = list(multi_line_string)
/home/m/m300602/miniconda3/envs/pyicon_py39_cartopy19/lib/python3.9/site-packages/cartopy/crs.py:280: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
for line in multi_line_string:
/home/m/m300602/miniconda3/envs/pyicon_py39_cartopy19/lib/python3.9/site-packages/cartopy/crs.py:347: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry.
if len(p_mline) > 0:
/home/m/m300602/miniconda3/envs/pyicon_py39_cartopy19/lib/python3.9/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
line_strings.extend(multi_line_string)
/home/m/m300602/miniconda3/envs/pyicon_py39_cartopy19/lib/python3.9/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry.
line_strings.extend(multi_line_string)
[18]:
ccrs_proj = ccrs.PlateCarree()
hca, hcb = pyic.arrange_axes(1,1, plot_cb=True, fig_size_fac=1.5, asp=0.5,
axlab_kw=None, projection=ccrs_proj,
)
ii=-1
ii+=1; ax=hca[ii]; cax=hcb[ii]
contfs = np.arange(0.,11e3+500.,500.)/1e3
clim = [contfs.min(), contfs.max()]
nn=10
hm1 = pyic.shade(lon_010deg[::nn], lat_010deg[::nn], res_010deg[::nn,::nn]/1e3,
ax=ax, cax=cax, clim=clim, #contfs=contfs,
# transform=ccrs_proj, rasterized=True,
# cmap=plt.cm.RdYlBu,
cbticks=contfs[::2],
)
ht = ax.set_title('ICON SMT: resolution [km]')
ax.add_patch(
Rectangle(xy=[lon_reg_2[0], lat_reg_2[0]],
width=lon_reg_2[1]-lon_reg_2[0], height=lat_reg_2[1]-lat_reg_2[0],
edgecolor='k',
facecolor='none',
linewidth=2,
transform=ccrs_proj) )
# ax.text(lon_reg_2[0], lat_reg_2[0]-5., 'region 1', ha='left', va='top', fontsize=10,
# bbox=dict(ec='none', fc='w', alpha=0.8, boxstyle='square,pad=0.3'),
# )
ax.annotate('region 1', xy=(lon_reg_2[1], lat_reg_2[0]), xycoords='data',
xytext=(-10, 20), textcoords='data',
fontsize=10,
arrowprops=dict(arrowstyle="->", facecolor='black'),
ha='center', va='center',
bbox=dict(ec='none', fc='w', alpha=0.8, boxstyle='square,pad=0.')
)
pyic.plot_settings(ax, template='global')
# savefig('resolution', format='pdf')
Vorticity / local Ro#
[18]:
# --- load times and flist
search_str = f'{run}_vort_f_50m_*.nc'
flist = np.array(glob.glob(path_data+search_str))
flist.sort()
timesv, flist_tsv, itsv = pyic.get_timesteps(flist, time_mode='float2date')
# timesv, flist_tsv, itsv = IcD.times, IcD.flist_ts, IcD.its
[19]:
itv = np.argmin((timesv-time0).astype(float)**2)
fpath = flist_tsv[itv]
f = Dataset(fpath, 'r')
var = 'vort_f_cells_50m'
vort = f.variables[var][itsv[itv],:]
f.close()
timesv[itv]
[19]:
numpy.datetime64('2010-03-09T01:00:00')
[20]:
vort = vort[ind_reg]
[21]:
# # ccrs_proj = ccrs.PlateCarree()
# ccrs_proj = None
# hca, hcb = pyic.arrange_axes(1, 1, plot_cb=True, asp=0.5, projection=ccrs_proj, fig_size_fac=3)
# ii=-1
# ii+=1; ax=hca[ii]; cax=hcb[ii]
# hm = pyic.shade(Tri, vort, ax=ax, cax=cax, clim=2)
# ax.set_xlim(lon_reg)
# ax.set_ylim(lat_reg)
Overview of vorticity#
[22]:
%%time
lon_reg = lon_reg_3
lat_reg = lat_reg_3
clon_reg, clat_reg, vertex_of_cell_reg, edge_of_cell_reg, ind_reg = pyic.crop_tripolar_grid(lon_reg, lat_reg,
clon, clat, vertex_of_cell, edge_of_cell)
Tri = matplotlib.tri.Triangulation(vlon, vlat, triangles=vertex_of_cell_reg)
CPU times: user 353 ms, sys: 320 ms, total: 672 ms
Wall time: 674 ms
[23]:
clon_reg2, clat_reg2, vertex_of_cell_reg2, edge_of_cell_reg2, ind_reg2 = pyic.crop_tripolar_grid(lon_reg_2, lat_reg_2,
clon, clat, vertex_of_cell, edge_of_cell)
[21]:
# --- load data
f = Dataset(flist_tsv[itv], 'r')
ro = f.variables['vort_f_cells_50m'][itsv[itv],:]
f.close()
# --- cut tripolar data
ro_reg = ro[ind_reg]
timesv[itv]
[21]:
numpy.datetime64('2010-03-09T01:00:00')
[22]:
# --- interpolate to rectgrid
ddnpz = np.load(fpath_ckdtree)
lon_002deg = ddnpz['lon']
lat_002deg = ddnpz['lat']
ro_002deg = pyic.apply_ckdtree(ro, fpath_ckdtree, coordinates='clat clon').reshape(lat_002deg.size, lon_002deg.size)
ro_002deg[ro_002deg==0.] = np.ma.masked
[23]:
# --- cut domain 1
i1r1 = (lon_002deg<lon_reg_1[0]).sum()
i2r1 = (lon_002deg<lon_reg_1[1]).sum()
j1r1 = (lat_002deg<lat_reg_1[0]).sum()
j2r1 = (lat_002deg<lat_reg_1[1]).sum()
lon_002deg_r1 = lon_002deg[i1r1:i2r1]
lat_002deg_r1 = lat_002deg[j1r1:j2r1]
ro_002deg_r1= ro_002deg[j1r1:j2r1,i1r1:i2r1]
# --- cut domain 2
i1r2 = (lon_002deg<lon_reg_2[0]).sum()
i2r2 = (lon_002deg<lon_reg_2[1]).sum()
j1r2 = (lat_002deg<lat_reg_2[0]).sum()
j2r2 = (lat_002deg<lat_reg_2[1]).sum()
lon_002deg_r2 = lon_002deg[i1r2:i2r2]
lat_002deg_r2 = lat_002deg[j1r2:j2r2]
ro_002deg_r2= ro_002deg[j1r2:j2r2,i1r2:i2r2]
# --- cut domain 3
i1r3 = (lon_002deg<lon_reg_3[0]).sum()
i2r3 = (lon_002deg<lon_reg_3[1]).sum()
j1r3 = (lat_002deg<lat_reg_3[0]).sum()
j2r3 = (lat_002deg<lat_reg_3[1]).sum()
lon_002deg_r3 = lon_002deg[i1r3:i2r3]
lat_002deg_r3 = lat_002deg[j1r3:j2r3]
ro_002deg_r3= ro_002deg[j1r3:j2r3,i1r3:i2r3]
# --- cut domain 3
i1r4 = (lon_002deg<lon_reg_4[0]).sum()
i2r4 = (lon_002deg<lon_reg_4[1]).sum()
j1r4 = (lat_002deg<lat_reg_4[0]).sum()
j2r4 = (lat_002deg<lat_reg_4[1]).sum()
lon_002deg_r4 = lon_002deg[i1r4:i2r4]
lat_002deg_r4 = lat_002deg[j1r4:j2r4]
ro_002deg_r4= ro_002deg[j1r4:j2r4,i1r4:i2r4]
[24]:
import maps_icon_smt_vort as smt
[25]:
# ccrs_proj = ccrs.PlateCarree()
# hca, hcb = pyic.arrange_axes(1,1,plot_cb=True, asp=0.5, fig_size_fac=2, projection=ccrs_proj)
# ii=-1
# ii+=1; ax=hca[ii]; cax=hcb[ii]
# hm3 = pyic.shade(Tri, ro_reg, ax=ax, cax=cax, clim=clim,
# transform=ccrs_proj, rasterized=False)
# ax.set_xlim(lon_reg_3)
# ax.set_ylim(lat_reg_3)
[35]:
%%time
clim = 1.5
ccrs_proj = ccrs.PlateCarree()
dpi = 250
fig, hca, cax = smt.make_axes(lon_reg_1, lat_reg_1, lon_reg_2, lat_reg_2, lon_reg_3, lat_reg_3, dpi)
ii=-1
hcb = [cax,0,0]
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm1 = pyic.shade(lon_002deg_r1, lat_002deg_r1, ro_002deg_r1, ax=ax, cax=cax, clim=clim,
transform=ccrs_proj, rasterized=False)
cax.set_ylabel('rel. vort. / plan. vort.')
ht = ax.set_title(timesv[itv])
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(lon_002deg_r2, lat_002deg_r2, ro_002deg_r2, ax=ax, cax=cax, clim=clim,
transform=ccrs_proj, rasterized=False)
ii+=1; ax=hca[ii]; cax=hcb[ii]
# hm3 = pyic.shade(lon_002deg_r3, lat_002deg_r3, ro_002deg_r3, ax=ax, cax=cax, clim=clim,
# transform=ccrs_proj, rasterized=False)
hm3 = pyic.shade(Tri, ro_reg, ax=ax, cax=cax, clim=clim,
transform=ccrs_proj, rasterized=False)
ax.set_xlim(lon_reg_3)
ax.set_ylim(lat_reg_3)
# savefig('main', dpi=dpi)
base_name = nb_name.split('/')[-1][:-6]
fpath_fig = f'{path_fig}{base_name}_main.jpg'
print(f'Writing file {fpath_fig}.')
plt.savefig(fpath_fig, dpi=dpi)
print(f'Done with file {fpath_fig}.')
Writing file ../pics/smtnatl_dissi_main.jpg.
Done with file ../pics/smtnatl_dissi_main.jpg.
CPU times: user 19.1 s, sys: 1.33 s, total: 20.4 s
Wall time: 20.5 s
(replace with pre-calculated dissipation?) Calculate dissipation#
The next cell is needed in the following
[40]:
# --- load times and flist
# search_str = f'{run}_T_S_sp_001-016_*.nc'
search_str = f'{run}_vn_sp_001-016_*.nc'
# search_str = f'{run}_vn_sp_081-096_*.nc'
flist = np.array(glob.glob(path_data+search_str))
flist.sort()
timesd, flist_tsd, itsd = pyic.get_timesteps(flist, time_mode='float2date')
itd = np.argmin((timesd-time0).astype(float)**2)
[31]:
vertex_of_cell.shape, vertex_of_cell.shape
[31]:
((59799625, 3), (59799625, 3))
[32]:
# edges_of_vertex = edges_of_vertex[vertex_of_cell,:]
# edge_orientation = edge_orientation[vertex_of_cell,:]
[33]:
%%time
rot_coeff = ( dual_edge_length[edges_of_vertex]*edge_orientation ) / dual_area[:,np.newaxis]
div_coeff = ( edge_length[edge_of_cell]*orientation_of_normal ) / cell_area_p[:,np.newaxis]
grad_coeff = (1./dual_edge_length)
rottr_coeff = (1./edge_length)
CPU times: user 4.06 s, sys: 4.01 s, total: 8.07 s
Wall time: 8.09 s
[87]:
fpath = flist_tsd[itd]
f = Dataset(fpath, 'r')
# var = 'vn001_sp'
# var = 'vn090_sp'
var = 'vn016_sp'
ve = f.variables[var][itsd[itd],:]
f.close()
print(timesd[itd])
2010-03-09T01:00:00
[36]:
%%time
curl_v = (ve[edges_of_vertex] * rot_coeff).sum(axis=1)
div_v = (ve[edge_of_cell] * div_coeff).sum(axis=1)
grad_div_v = (div_v[adjacent_cell_of_edge[:,1]]-div_v[adjacent_cell_of_edge[:,0]])*grad_coeff
curlt_curl_v = (curl_v[edge_vertices[:,1]]-curl_v[edge_vertices[:,0]])*rottr_coeff
diss = -grad_div_v**2 - curlt_curl_v**2
CPU times: user 7.21 s, sys: 3.3 s, total: 10.5 s
Wall time: 10.5 s
[37]:
A4 = 0.04 * 500**3
diss *= A4
[38]:
%%time
diss_i = diss[edge_of_cell].sum(axis=1)/3.
CPU times: user 1.34 s, sys: 694 ms, total: 2.03 s
Wall time: 2.04 s
[39]:
diss_ireg = diss_i[ind_reg]
divv_ireg = div_v[ind_reg]
[40]:
curl_v.shape, vlon.shape, clon.shape, diss_i.shape
[40]:
((30008985,), (30008985,), (59799625,), (59799625,))
[41]:
%%time
earth_angular_velocity = 7.29212e-05
fv = 2.* earth_angular_velocity * np.sin(vlat*np.pi/180.)
ro_self = curl_v/fv
CPU times: user 665 ms, sys: 477 ms, total: 1.14 s
Wall time: 1.14 s
[42]:
%%time
ro_selfi = ro_self[vertex_of_cell].sum(axis=1)/3.
CPU times: user 1.29 s, sys: 715 ms, total: 2 s
Wall time: 2 s
[43]:
# fpath_ckdtree = '/mnt/lustre01/work/mh0033/m300602/icon/grids/smt/new_ckdtree/rectgrids/smt_res0.02_180W-180E_90S-90N.npz'
[44]:
%%time
diss_002deg = pyic.apply_ckdtree(diss_i, fpath_ckdtree, coordinates='clat clon').reshape(lat_002deg.size, lon_002deg.size)
diss_002deg[diss_002deg==0.] = np.ma.masked
diss_002deg_r2= diss_002deg[j1r2:j2r2,i1r2:i2r2]
divv_002deg = pyic.apply_ckdtree(div_v, fpath_ckdtree, coordinates='clat clon').reshape(lat_002deg.size, lon_002deg.size)
divv_002deg[divv_002deg==0.] = np.ma.masked
divv_002deg_r2= divv_002deg[j1r2:j2r2,i1r2:i2r2]
rotv_002deg = pyic.apply_ckdtree(ro_selfi, fpath_ckdtree, coordinates='clat clon').reshape(lat_002deg.size, lon_002deg.size)
rotv_002deg[rotv_002deg==0.] = np.ma.masked
rotv_002deg_r2= rotv_002deg[j1r2:j2r2,i1r2:i2r2]
CPU times: user 7.54 s, sys: 5.54 s, total: 13.1 s
Wall time: 13.2 s
[45]:
# rotv_002deg = pyic.apply_ckdtree(ro_self, fpath_ckdtree, coordinates='vlat vlon', radius_of_influence=1.).reshape(lat_002deg.size, lon_002deg.size)
# rotv_002deg[rotv_002deg==0.] = np.ma.masked
# rotv_002deg_r2= rotv_002deg[j1r2:j2r2,i1r2:i2r2]
[46]:
help(pyic.arrange_axes)
Help on function arrange_axes in module pyicon.pyicon_plotting:
arrange_axes(nx, ny, sharex=True, sharey=False, xlabel='', ylabel='', do_axes_labels=True, axlab_kw={}, plot_cb=True, projection=None, asp=1.0, sasp=0.0, wax='auto', hax=4.0, dfigl=0.0, dfigr=0.0, dfigt=0.0, dfigb=0.0, daxl=1.8, daxr=0.8, daxt=0.8, daxb=1.2, dcbl=-0.5, dcbr=1.4, dcbt=0.0, dcbb=0.5, wcb=0.5, hcb='auto', fig_size_fac=1.0, f_wax=1.0, f_hax=1.0, f_wcb=1.0, f_hcb=1.0, f_dfigl=1.0, f_dfigr=1.0, f_dfigt=1.0, f_dfigb=1.0, f_daxl=1.0, f_daxr=1.0, f_daxt=1.0, f_daxb=1.0, f_dcbl=1.0, f_dcbr=1.0, f_dcbt=1.0, f_dcbb=1.0, fs_label=10.0, fs_title=12.0, fs_ticks=10.0, f_fs=1, reverse_order=False)
Plot vorticity / dissipation North Atlantic#
[47]:
ccrs_proj = ccrs.PlateCarree()
# ccrs_proj = None
hca, hcb = pyic.arrange_axes(1, 2, plot_cb=True, asp=0.5, projection=ccrs_proj, fig_size_fac=2,
sharex=False,
dfigl=-0.34, dfigr=-0.66)
ii=-1
# clim = 2e-7
# clim = 'sym'
ii+=1; ax=hca[ii]; cax=hcb[ii]
# hm = pyic.shade(Tri, diss_ireg, ax=ax, cax=cax, clim=clim)
hm2 = pyic.shade(lon_002deg_r2, lat_002deg_r2, ro_002deg_r2, ax=ax, cax=cax, clim=2,
transform=ccrs_proj,)
# hm2 = pyic.shade(lon_002deg_r2, lat_002deg_r2, rotv_002deg_r2, ax=ax, cax=cax, clim=2,
# transform=ccrs_proj, rasterized=False)
ax.set_title('local Rossby number')
ax.set_title(f'{time0}', loc='right')
# ii+=1; ax=hca[ii]; cax=hcb[ii]
# hm2 = pyic.shade(lon_002deg_r2, lat_002deg_r2, divv_002deg_r2, ax=ax, cax=cax, clim=2e-4,
# transform=ccrs_proj, rasterized=False)
# ax.set_title('div. of hor. velocity [1/s]')
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(lon_002deg_r2, lat_002deg_r2, -diss_002deg_r2, ax=ax, cax=cax, clim=[-9,-5], cmap='Blues',
transform=ccrs_proj, logplot=True)
ax.set_title('log$_{10}$(kin. energy dissipation) [m$^2$/s$^3$]')
for ax in hca:
# ax.set_xlim(lon_reg_2)
# ax.set_ylim(lat_reg_2)
pyic.plot_settings(ax, xlim=lon_reg_2, ylim=lat_reg_2)
ax.text(0.02,0.9, 'z = 50m', transform=ax.transAxes)
ax.add_patch(
Rectangle(xy=[lon_reg_3[0], lat_reg_3[0]],
width=lon_reg_3[1]-lon_reg_3[0], height=lat_reg_3[1]-lat_reg_3[0],
edgecolor='k',
facecolor='none',
linewidth=2,
transform=ccrs_proj) )
ax.text(lon_reg_3[0], lat_reg_3[0]-0.1, 'region 2', ha='left', va='top', fontsize=10,
bbox=dict(ec='none', fc='w', alpha=0.8, boxstyle='square,pad=0.'),
)
ax.add_patch(
Rectangle(xy=[lon_reg_4[0], lat_reg_4[0]],
width=lon_reg_4[1]-lon_reg_4[0], height=lat_reg_4[1]-lat_reg_4[0],
edgecolor='k',
facecolor='none',
linewidth=2,
transform=ccrs_proj) )
ax.text(lon_reg_4[0], lat_reg_4[0]-0.1, 'region 3', ha='left', va='top', fontsize=10,
bbox=dict(ec='none', fc='w', alpha=0.8, boxstyle='square,pad=0.'),
)
ax.annotate('Gulf Stream front', xy=(-68.3, 39), xycoords='data',
xytext=(-68.54, 41), textcoords='data',
fontsize=10,
arrowprops=dict(arrowstyle="->", facecolor='black'),
ha='center', va='center',
bbox=dict(ec='none', fc='w', alpha=0.8, boxstyle='square,pad=0.')
)
ax.text(-61, 35.5, 'enhanced\nsubmesoscale dynamics', ha='left', va='top', fontsize=10,
bbox=dict(ec='none', fc='w', alpha=0.8, boxstyle='square,pad=0.3'),
)
ax.text(-61, 42.5, 'reduced\nsubmesoscale dynamics', ha='left', va='top', fontsize=10,
bbox=dict(ec='none', fc='w', alpha=0.8, boxstyle='square,pad=0.3'),
)
# savefig('dissipation', format='pdf')
[48]:
res_xr.dims
[48]:
('cell',)
[49]:
fpath_ckdtree
[49]:
'/work/mh0033/m300602/icon/grids/smt/ckdtree/rectgrids/smt_res0.02_180W-180E_90S-90N.npz'
[50]:
resi = pyic.interp_to_rectgrid_xr(res_xr.rename(cell='ncells'), '/work/mh0033/m300602/icon/grids/smt/ckdtree/rectgrids/smt_res0.02_180W-180E_90S-90N.nc',
lon_reg=lon_reg_1, lat_reg=lat_reg_1)
[51]:
ccrs_proj = ccrs.PlateCarree()
# ccrs_proj = None
hca, hcb = pyic.arrange_axes(1, 1, plot_cb=True, asp=0.325, projection=ccrs_proj, fig_size_fac=1.2,
sharex=False,
dfigl=-0.34, dfigr=-0.66)
ii=-1
# clim = 2e-7
# clim = 'sym'
ii+=1; ax=hca[ii]; cax=hcb[ii]
# hm = pyic.shade(Tri, diss_ireg, ax=ax, cax=cax, clim=clim)
hm2 = pyic.shade(lon_002deg_r2, lat_002deg_r2, ro_002deg_r2, ax=ax, cax=cax, clim=2,
transform=ccrs_proj,)
# hm2 = pyic.shade(lon_002deg_r2, lat_002deg_r2, rotv_002deg_r2, ax=ax, cax=cax, clim=2,
# transform=ccrs_proj, rasterized=False)
ax.set_title('local Rossby number')
ax.set_title(f'{time0}', loc='right')
for ax in hca:
# ax.set_xlim(lon_reg_2)
# ax.set_ylim(lat_reg_2)
pyic.plot_settings(ax, xlim=lon_reg_2, ylim=lat_reg_2)
ax.text(0.02,0.9, 'z = 50m', transform=ax.transAxes)
ax.add_patch(
Rectangle(xy=[lon_reg_3[0], lat_reg_3[0]],
width=lon_reg_3[1]-lon_reg_3[0], height=lat_reg_3[1]-lat_reg_3[0],
edgecolor='k',
facecolor='none',
linewidth=2,
transform=ccrs_proj) )
ax.text(lon_reg_3[0], lat_reg_3[1]+0.1, 'zoom area', ha='left', va='bottom', fontsize=10,
bbox=dict(ec='none', fc='w', alpha=0.8, boxstyle='square,pad=0.'),
)
ax.annotate('Gulf Stream front', xy=(-68.3, 39), xycoords='data',
xytext=(-68.54, 41), textcoords='data',
fontsize=10,
arrowprops=dict(arrowstyle="->", facecolor='black'),
ha='center', va='center',
bbox=dict(ec='none', fc='w', alpha=0.8, boxstyle='square,pad=0.')
)
ax.set_ylim(35.5,42)
# savefig('dissipation', format='pdf')
Cl = ax.contour(resi.lon, resi.lat, resi, np.arange(500,1000,50), colors='0.3')
ax.clabel(Cl, fmt='%.fm', fontsize=8)
asp = ax.get_position().height*plt.gcf().get_figheight() / ( ax.get_position().width*plt.gcf().get_figwidth() )
asp
[51]:
0.32499999999999996
Histogram of Ro#
[27]:
# hca, hcb = pyic.arrange_axes(1, 1, plot_cb=False, asp=0.5, fig_size_fac=2)
# ii=-1
# ii+=1; ax=hca[ii]; cax=hcb[ii]
# hist, bins = np.histogram(ro_002deg_r3, bins=np.linspace(-4,4,51))
# hdens = hist/hist.sum()/(bins[1]-bins[0])
# print((hdens*(bins[1]-bins[0])).sum())
# ax.plot(0.5*(bins[1:]+bins[:-1]), hdens, marker=None, label='r3')
# hist, bins = np.histogram(ro_002deg_r4, bins=np.linspace(-4,4,51))
# hdens = hist/hist.sum()/(bins[1]-bins[0])
# print((hdens*(bins[1]-bins[0])).sum())
# ax.plot(0.5*(bins[1:]+bins[:-1]), hdens, marker=None, label='r4')
# ax.legend()
# ax.set_xlim(-2,3)
# ax.grid(True)
[28]:
lon_reg = lon_reg_2
lat_reg = lat_reg_2
ind_mask = ( (clon>lon_reg[0])
& (clon<=lon_reg[1])
& (clat>lat_reg[0])
& (clat<=lat_reg[1]) )
ro_r2 = ro[ind_mask]
ro_r2[ro_r2==0] = np.ma.masked
lon_reg = lon_reg_3
lat_reg = lat_reg_3
ind_mask = ( (clon>lon_reg[0])
& (clon<=lon_reg[1])
& (clat>lat_reg[0])
& (clat<=lat_reg[1]) )
ro_r3 = ro[ind_mask]
ro_r3[ro_r3==0] = np.ma.masked
lon_reg = lon_reg_4
lat_reg = lat_reg_4
ind_mask = ( (clon>lon_reg[0])
& (clon<=lon_reg[1])
& (clat>lat_reg[0])
& (clat<=lat_reg[1]) )
ro_r4 = ro[ind_mask]
ro_r4[ro_r4==0] = np.ma.masked
[29]:
def plot_histogram(ax, roi, label='asdf', c='tab:blue'):
roi2 = roi.copy()
roi2 = roi2[roi2.mask==False]
hist, bins = np.histogram(roi2, bins=np.linspace(-4,4,101))
hdens = hist/hist.sum()/(bins[1]-bins[0])
ax.plot(0.5*(bins[1:]+bins[:-1]), hdens, marker=None, label=label, color=c)
print((hdens*(bins[1]-bins[0])).sum())
return
[30]:
hca, hcb = pyic.arrange_axes(1, 1, plot_cb=False, asp=0.5, fig_size_fac=2)
ii=-1
ii+=1; ax=hca[ii]; cax=hcb[ii]
plot_histogram(ax, ro_r2, label='region 1', c='0.7')
plot_histogram(ax, ro_r3, label='region 2', c='tab:blue')
plot_histogram(ax, ro_r4, label='region 3', c='tab:orange')
ax.set_xlim(-2,3)
ax.grid(True)
ax.set_title('histogramm of local Rossby number at 50m depth')
ax.set_xlabel('Rossby number')
ax.set_ylabel('points in interval\n / all points / interval size')
ax.legend()
# ax.axvline(-1, color='k')
# ax.axvline(1, color='k')
ax.axvline(0, color='k')
plt.axvspan(-2, -1, color='b', alpha=0.1, lw=0)
plt.axvspan(1, 3, color='b', alpha=0.1, lw=0)
ax.text(-1.5, 0.3, 'symmetrically\nunstable', ha='center', va='center')
ax.text(-1.5, 0.7, 'ageostrophic\nregime', ha='center', va='center')
ax.text(1.5, 0.7, 'ageostrophic\nregime', ha='center', va='center')
# ax.annotate('local max', xy=(2, 1), xytext=(3, 1.5),
# arrowprops=dict(facecolor='black', shrink=0.05),
# )
# savefig('histogram', format='pdf')
1.0
1.0
1.0
[30]:
Text(1.5, 0.7, 'ageostrophic\nregime')
(del?) Remap velocities#
[56]:
grid_sphere_radius = 6.371229e6
[57]:
%%time
dist_vector = edge_cart_vec[edge_of_cell,:] - cell_cart_vec[:,np.newaxis,:]
norm = np.sqrt(pyic.scalar_product(dist_vector,dist_vector,dim=2))
prime_edge_length = edge_length/grid_sphere_radius
fixed_vol_norm = ( 0.5 * norm
* (prime_edge_length[edge_of_cell]))
fixed_vol_norm = fixed_vol_norm.sum(axis=1)
CPU times: user 8.92 s, sys: 6.49 s, total: 15.4 s
Wall time: 15.5 s
[58]:
%%time
dist_vector = ( edge_cart_vec[edge_of_cell,:]
- cell_cart_vec[:,np.newaxis,:] )
edge2cell_coeff_cc = ( dist_vector
* (edge_length[edge_of_cell,np.newaxis] / grid_sphere_radius)
* orientation_of_normal[:,:,np.newaxis] )
CPU times: user 10.4 s, sys: 7.05 s, total: 17.5 s
Wall time: 17.5 s
[59]:
edge2cell_coeff_cc.shape, ve.shape
[59]:
((59799625, 3, 3), (89813639,))
[60]:
edge2cell_coeff_cc.shape
[60]:
(59799625, 3, 3)
[61]:
%%time
p_vn_c = ( edge2cell_coeff_cc[:,:,:]
* ve[edge_of_cell,np.newaxis]
# * prism_thick_e[iz,IcD.edge_of_cell,np.newaxis]
* dzw[iz]
).sum(axis=1)
p_vn_c *= 1./(fixed_vol_norm[:,np.newaxis]*dzw[iz])
CPU times: user 15.8 s, sys: 5.99 s, total: 21.8 s
Wall time: 21.9 s
[62]:
%%time
sinLon = np.sin(clon*np.pi/180.)
cosLon = np.cos(clon*np.pi/180.)
sinLat = np.sin(clat*np.pi/180.)
cosLat = np.cos(clat*np.pi/180.)
u1 = p_vn_c[:,0]
u2 = p_vn_c[:,1]
u3 = p_vn_c[:,2]
uo = u2*cosLon - u1*sinLon
vo = -(u1*cosLon + u2*sinLon)*sinLat + u3*cosLat
CPU times: user 4.82 s, sys: 3.18 s, total: 8 s
Wall time: 8.02 s
[ ]:
uo_002deg = pyic.apply_ckdtree(uo, fpath_ckdtree, coordinates='clat clon').reshape(lat_002deg.size, lon_002deg.size)
uo_002deg[uo_002deg==0.] = np.ma.masked
uo_002deg_r2= uo_002deg[j1r2:j2r2,i1r2:i2r2]
vo_002deg = pyic.apply_ckdtree(vo, fpath_ckdtree, coordinates='clat clon').reshape(lat_002deg.size, lon_002deg.size)
vo_002deg[vo_002deg==0.] = np.ma.masked
vo_002deg_r2= vo_002deg[j1r2:j2r2,i1r2:i2r2]
[ ]:
ke_002deg_r2 = 0.5*(uo_002deg_r2**2+vo_002deg_r2**2)
[ ]:
ccrs_proj = ccrs.PlateCarree()
# ccrs_proj = None
hca, hcb = pyic.arrange_axes(1, 3, plot_cb=True, asp=0.5, projection=ccrs_proj, fig_size_fac=2)
ii=-1
clim = 1.5
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(lon_002deg_r2, lat_002deg_r2, uo_002deg_r2, ax=ax, cax=cax, clim=clim,
transform=ccrs_proj, rasterized=False)
ax.set_title('zon. ve. [m/s]')
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(lon_002deg_r2, lat_002deg_r2, vo_002deg_r2, ax=ax, cax=cax, clim=clim,
transform=ccrs_proj, rasterized=False)
ax.set_title('mer. ve. [m/s]')
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(lon_002deg_r2, lat_002deg_r2, ke_002deg_r2, ax=ax, cax=cax, clim=[0,1],
transform=ccrs_proj, rasterized=False)
ax.set_title('kin. energy [m$^2$/s$^2$]')
for ax in hca:
pyic.plot_settings(ax, xlim=lon_reg_2, ylim=lat_reg_2)
# savefig('dissipation', format='png')
Dissipation#
Crop grid#
[31]:
# lon_reg_2 = [-75, -55]
# lat_reg_2 = [33, 43]
lon_reg = lon_reg_3
lat_reg = lat_reg_3
lon, lat = lon_002deg_r3, lat_002deg_r3
clon_reg, clat_reg, vertex_of_cell_reg, edge_of_cell_reg, ind_reg = pyic.crop_tripolar_grid(lon_reg, lat_reg,
clon, clat, vertex_of_cell, edge_of_cell)
[32]:
# path_scratch = '/scratch/m/m300602/tmp/'
Load temperature#
[33]:
# --- load times and flist
search_str = f'{run}_T_S_sp_001-016_*.nc'
# search_str = f'{run}_vort_f_50m_*.nc'
flist = np.array(glob.glob(path_data+search_str))
flist.sort()
times, flist_ts, its = pyic.get_timesteps(flist, time_mode='float2date')
[34]:
varfile = 'T_S'
layers = ( ['sp_001-016']*16 + ['sp_017-032']*16 + ['sp_033-048']*16
+ ['sp_049-064']*16 + ['sp_065-080']*16 + ['sp_081-096']*16 + ['sp_097-112']*16)
# nz = len(layers)
tstr = '20100309T010000Z'
[35]:
path_data = '/work/mh0287/from_Mistral/mh0287/users/leonidas/icon/ngSMT/results/2010-03/'
[36]:
# levs = np.concatenate((np.arange(0,48,2),np.arange(48,48+40,1))).astype(int)
levs = np.arange(88).astype(int)
# levs.size
# np.arange(0,nz/2,2),np.arange(nz/2,nz/2+48,1)
nz = levs.size
levs, nz, depthc[levs[-1]]
[36]:
(array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67,
68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84,
85, 86, 87]),
88,
695.25)
This takes a long time (1min 20s)…
[41]:
%%time
it = np.argmin((times-time0).astype(float)**2)
temp_reg = np.ma.zeros((nz,ind_reg2.size))
for kk, lev in enumerate(levs):
fname = f'{run}_{varfile}_{layers[lev]}_{tstr}.nc'
fpath = f'{path_data}{fname}'
# var = 'vn001_sp'
# var = 'vn090_sp'
# var = 'vn016_sp'
var = f'T{lev+1:03d}_sp'
print(f'kk = {kk}/{nz}; fname = {fname}; var = {var}')
fi = Dataset(fpath, 'r')
temp = fi.variables[var][itsd[itd],:]
# print(fi.variables[var].shape)
fi.close()
temp_reg[kk,:] = temp[ind_reg2]
# f = Dataset(fpath, 'r')
# var = 'T016_sp'
# # var = 'vort_f_cells_50m'
# temp = f.variables[var][its[it],:]
# temp_reg = temp[ind_reg2]
# f.close()
kk = 0/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T001_sp
kk = 1/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T002_sp
kk = 2/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T003_sp
kk = 3/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T004_sp
kk = 4/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T005_sp
kk = 5/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T006_sp
kk = 6/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T007_sp
kk = 7/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T008_sp
kk = 8/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T009_sp
kk = 9/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T010_sp
kk = 10/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T011_sp
kk = 11/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T012_sp
kk = 12/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T013_sp
kk = 13/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T014_sp
kk = 14/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T015_sp
kk = 15/88; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T016_sp
kk = 16/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T017_sp
kk = 17/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T018_sp
kk = 18/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T019_sp
kk = 19/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T020_sp
kk = 20/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T021_sp
kk = 21/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T022_sp
kk = 22/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T023_sp
kk = 23/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T024_sp
kk = 24/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T025_sp
kk = 25/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T026_sp
kk = 26/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T027_sp
kk = 27/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T028_sp
kk = 28/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T029_sp
kk = 29/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T030_sp
kk = 30/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T031_sp
kk = 31/88; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T032_sp
kk = 32/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T033_sp
kk = 33/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T034_sp
kk = 34/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T035_sp
kk = 35/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T036_sp
kk = 36/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T037_sp
kk = 37/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T038_sp
kk = 38/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T039_sp
kk = 39/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T040_sp
kk = 40/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T041_sp
kk = 41/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T042_sp
kk = 42/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T043_sp
kk = 43/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T044_sp
kk = 44/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T045_sp
kk = 45/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T046_sp
kk = 46/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T047_sp
kk = 47/88; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T048_sp
kk = 48/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T049_sp
kk = 49/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T050_sp
kk = 50/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T051_sp
kk = 51/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T052_sp
kk = 52/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T053_sp
kk = 53/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T054_sp
kk = 54/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T055_sp
kk = 55/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T056_sp
kk = 56/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T057_sp
kk = 57/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T058_sp
kk = 58/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T059_sp
kk = 59/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T060_sp
kk = 60/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T061_sp
kk = 61/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T062_sp
kk = 62/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T063_sp
kk = 63/88; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T064_sp
kk = 64/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T065_sp
kk = 65/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T066_sp
kk = 66/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T067_sp
kk = 67/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T068_sp
kk = 68/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T069_sp
kk = 69/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T070_sp
kk = 70/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T071_sp
kk = 71/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T072_sp
kk = 72/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T073_sp
kk = 73/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T074_sp
kk = 74/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T075_sp
kk = 75/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T076_sp
kk = 76/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T077_sp
kk = 77/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T078_sp
kk = 78/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T079_sp
kk = 79/88; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T080_sp
kk = 80/88; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T081_sp
kk = 81/88; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T082_sp
kk = 82/88; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T083_sp
kk = 83/88; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T084_sp
kk = 84/88; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T085_sp
kk = 85/88; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T086_sp
kk = 86/88; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T087_sp
kk = 87/88; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T088_sp
CPU times: user 45.7 s, sys: 22.6 s, total: 1min 8s
Wall time: 1min 20s
Load data from pp script#
[42]:
fpath = '/work/mh0033/m300602/proj_smtnatl/pp/pp_dissipation_all.nc'
f = Dataset(fpath, 'r')
depthc_sec = f.variables['depthc'][:]
iz = (depthc_sec<depth0).sum()
print(iz, depthc_sec[iz])
kk = iz
diss_i = f.variables['diss_i'][kk,:]
vort_i = f.variables['vort_i'][kk,:]
div_v = f.variables['div_v'][kk,:]
uo = f.variables['uo'][kk,:]
vo = f.variables['vo'][kk,:]
f.close()
temp_lev = temp_reg[kk,:]
16 51.5
[43]:
ny = 200
nx = ny*2
lon = np.linspace(lon_reg[0], lon_reg[1], nx)
lat = np.linspace(lat_reg[0], lat_reg[1], ny)
[46]:
lon_reg
[46]:
[-68.5, -66.5]
[48]:
lon_o, lat_o = np.meshgrid(lon, lat)
lon_o = lon_o.flatten()
lat_o = lat_o.flatten()
lon_o.shape, lat_o.shape
[48]:
((80000,), (80000,))
[49]:
clon_reg.max(), clon_reg.min(), lon_o.max(), lon_o.min()
[49]:
(-66.50003403874472, -68.49996742078073, -66.5, -68.5)
[50]:
%%time
dckdtree_c, ickdtree_c = pyic.calc_ckdtree(lon_i=clon_reg2, lat_i=clat_reg2,
lon_o=lon_o, lat_o=lat_o,
n_nearest_neighbours=1,
)
CPU times: user 2.3 s, sys: 359 ms, total: 2.65 s
Wall time: 2.8 s
[51]:
# region points
diss_ii = diss_i[ickdtree_c]
diss_ii = diss_ii.reshape(lat.size, lon.size)
vort_ii = vort_i[ickdtree_c]
vort_ii = vort_ii.reshape(lat.size, lon.size)
div_vi = div_v[ickdtree_c]
div_vi = div_vi.reshape(lat.size, lon.size)
uoi = uo[ickdtree_c]
uoi = uoi.reshape(lat.size, lon.size)
voi = vo[ickdtree_c]
voi = voi.reshape(lat.size, lon.size)
tempi = temp_lev[ickdtree_c]
tempi = tempi.reshape(lat.size, lon.size)
kei = 0.5*(uoi**2+voi**2)
[52]:
# p1 = [-70.15, 34.68]
# p2 = [-68.55, 35.75]
p1 = [-68.21, 37.2]
p2 = [-66.58, 36.53]
lon_sec, lat_sec, dist_sec = pyic.derive_section_points(p1, p2, npoints=201)
[54]:
ccrs_proj = ccrs.PlateCarree()
hca, hcb = pyic.arrange_axes(2, 2, plot_cb=True, asp=0.5, projection=ccrs_proj, fig_size_fac=1.2, sharey=True)
ii=-1
clim = 1e-6
# ii+=1; ax=hca[ii]; cax=hcb[ii]
# hm2 = pyic.shade(lon, lat, uoi, ax=ax, cax=cax, clim=1.5,
# transform=ccrs_proj, rasterized=False)
# ax.set_title('zon. ve. [m/s]')
# ii+=1; ax=hca[ii]; cax=hcb[ii]
# hm2 = pyic.shade(lon, lat, ke, ax=ax, cax=cax, clim=1.5,
# transform=ccrs_proj, rasterized=False)
# ax.set_title('mer. ve. [m/s]')
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(lon, lat, tempi, ax=ax, cax=cax, clim=[17,18],
transform=ccrs_proj, rasterized=False)
ax.set_title('temperature [$^o$C]')
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(lon, lat, vort_ii, ax=ax, cax=cax, clim=2,
transform=ccrs_proj, rasterized=False)
ax.set_title('local Rossby number')
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(lon, lat, kei, ax=ax, cax=cax, clim=[-2,0],
transform=ccrs_proj, rasterized=False, cmap='Reds', logplot=True)
ax.set_title('kin. energy [m$^2$/s$^2$]')
iix, iiy = 10, 10
ax.quiver(lon[::iix], lat[::iiy], uoi[::iiy,::iix], voi[::iiy,::iix])
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm = pyic.shade(lon, lat, -diss_ii, ax=ax, cax=cax, clim=[-9,-6],
transform=ccrs_proj, rasterized=False, cmap='Blues', logplot=True)
ax.set_title(f'kin. energy dissipation at {depthc[0]}m [m$^2$/s$^3$]')
# ii+=1; ax=hca[ii]; cax=hcb[ii]
# hm2 = pyic.shade(lon, lat, div_vi, ax=ax, cax=cax, clim=2e-4,
# transform=ccrs_proj, rasterized=False)
# ax.set_title('hor. velocity divergence [1/s]')
for ax in hca:
pyic.plot_settings(ax, xlim=lon_reg_3, ylim=lat_reg_3)
ax.plot(lon_sec, lat_sec, transform=ccrs_proj, color='k')
[55]:
%%time
dckd_sec_c, ickd_sec_c = pyic.calc_ckdtree(lon_i=clon_reg2, lat_i=clat_reg2,
lon_o=lon_sec, lat_o=lat_sec,
n_nearest_neighbours=1,
)
CPU times: user 2.24 s, sys: 442 ms, total: 2.68 s
Wall time: 2.68 s
[56]:
%%time
kk = 0
fpath = '/work/mh0033/m300602/proj_smtnatl/pp/pp_dissipation_all.nc'
f = Dataset(fpath, 'r')
nsec = ickd_sec_c.size
diss_sec = np.ma.zeros((nz,nsec))
vort_sec = np.ma.zeros((nz,nsec))
temp_sec = np.ma.zeros((nz,nsec))
uo_sec = np.ma.zeros((nz,nsec))
vo_sec = np.ma.zeros((nz,nsec))
for kk in range(nz):
# print(f'kk = {kk}/{nz}')
diss_sec[kk,:] = f.variables['diss_i'][kk,:][ickd_sec_c]
vort_sec[kk,:] = f.variables['vort_i'][kk,:][ickd_sec_c]
# divv_sec[kk,:] = f.variables['div_v'][kk,:][ickd_sec_c]
uo_sec[kk,:] = f.variables['uo'][kk,:][ickd_sec_c]
vo_sec[kk,:] = f.variables['vo'][kk,:][ickd_sec_c]
temp_sec[kk,:] = temp_reg[kk,:][ickd_sec_c]
f.close()
ke_sec = 0.5*(uo_sec**2+vo_sec**2)
CPU times: user 713 ms, sys: 2.95 s, total: 3.66 s
Wall time: 24.6 s
[92]:
# # section points
# diss_ii = diss_i[ickdtree_c]
# diss_ii = diss_ii.reshape(lat.size, lon.size)
# vort_ii = vort_i[ickdtree_c]
# vort_ii = vort_ii.reshape(lat.size, lon.size)
# div_vi = div_v[ickdtree_c]
# div_vi = div_vi.reshape(lat.size, lon.size)
# uoi = uo[ickdtree_c]
# uoi = uoi.reshape(lat.size, lon.size)
# voi = vo[ickdtree_c]
# voi = voi.reshape(lat.size, lon.size)
# tempi = temp_lev[ickdtree_c]
# tempi = tempi.reshape(lat.size, lon.size)
# kei = 0.5*(uoi**2+voi**2)
Plots#
Section plots#
[57]:
hca, hcb = pyic.arrange_axes(1, 4, plot_cb=True, asp=0.5, projection=None, fig_size_fac=1.2, sharex=True, sharey=False)
ii=-1
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(dist_sec, depthc[levs], temp_sec, clevs=[11,13,15,16.6,16.7,16.8,16.9,17.1,17.2,17.3,17.4,17.5], conts=np.arange(10.6,16.6,1), ax=ax, cax=cax, clim=[10,17.5])
ax.set_title('temperature [$^o$C]')
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(dist_sec, depthc[levs], vort_sec, ax=ax, cax=cax, clim=2)
ax.set_title('local Rossby number')
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(dist_sec, depthc[levs], -diss_sec, ax=ax, cax=cax, clim=[0,1e-6], cmap='Blues')
ax.set_title(f'kin. energy dissipation at {depthc[0]}m [m$^2$/s$^3$]')
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm = pyic.shade(dist_sec, depthc[levs], ke_sec, ax=ax, cax=cax, clim=[0,2], cmap='Reds')
ax.set_title(f'kin. energy dissipation at {depthc[0]}m [m$^2$/s$^3$]')
for ax in hca:
# ax.set_ylim(depthc_sec.max(),0)
ax.set_ylim(700,0)
hca[0].set_ylim(700,0)
[57]:
(700.0, 0.0)
Horizontal and section plots of two SM eddies#
[58]:
kmld = (temp_sec>(temp_sec[0:1,:]-0.5)).sum(axis=0)
mld = depthc_sec[kmld]
[59]:
hca, hcb = pyic.arrange_axes(3, 2, plot_cb=[1,1,1,1,1,1], asp=0.5,
projection=[ccrs_proj, ccrs_proj, ccrs_proj, None, None, None],
fig_size_fac=1.2, sharex=False, sharey=True, f_daxl=0.7, dfigl=0.5, dfigb=1.5, daxb=0.3)
# hca, hcb = pyic.arrange_axes(2, 3, plot_cb=[0,1,0,1,0,1], asp=0.5,
# projection=[ccrs_proj, None, ccrs_proj, None, ccrs_proj, None],
# fig_size_fac=1.2, sharex=True, sharey=False, f_daxl=0.7, dfigl=0.5)
ii=-1
# hca = np.array(hca)[[0,2,4,1,3,5]]
# hcb = np.array(hcb)[[0,2,4,1,3,5]]
climT=[16.8,17.8]
climR=2
# climD=[0,4e-8]
climD=[-9,-7.0]
# --- hor. plots
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(lon, lat, tempi, ax=ax, cax=cax, clim=climT,
transform=ccrs_proj)
ax.set_title('temperature [$^o$C]')
iix, iiy = 10, 10
ax.quiver(lon[::iix], lat[::iiy], uoi[::iiy,::iix], voi[::iiy,::iix])
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(lon, lat, vort_ii, ax=ax, cax=cax, clim=climR,
transform=ccrs_proj)
ax.set_title('local Rossby number')
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm = pyic.shade(lon, lat, -diss_ii, ax=ax, cax=cax, clim=climD,
transform=ccrs_proj, cmap='RdBu_r', logplot=True)
ax.set_title('log$_{10}$(kin. energy dissipation) [m$^2$/s$^3$]')
# --- sections
ii+=1; ax=hca[ii]; cax=hcb[ii]
depthc_sec = depthc[levs]
hm2 = pyic.shade(dist_sec/1e3, depthc_sec, temp_sec, ax=ax, cax=cax, clim=climT)
Cl = ax.contour(dist_sec/1e3, depthc_sec, temp_sec, [11.8,12.8,13.8,14.8,15.8,16.8], colors='0.7')
ax.clabel(Cl, fmt='%.1f', fontsize=8)
# ax.set_title('temperature [$^o$C]')
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(dist_sec/1e3, depthc_sec, vort_sec, ax=ax, cax=cax, clim=climR)
# ax.set_title('local Rossby number')
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(dist_sec/1e3, depthc_sec, -diss_sec, ax=ax, cax=cax, clim=climD, cmap='RdBu_r', logplot=True)
# ax.set_title('log$_{10}$(kin. energy dissipation) [m$^2$/s$^3$]')
# ii+=1; ax=hca[ii]; cax=hcb[ii]
# hm = pyic.shade(dist_sec, depthc_sec, ke_sec, ax=ax, cax=cax, clim=[0,2], cmap='Reds')
# ax.set_title(f'kin. energy dissipation at {depthc[0]}m [m$^2$/s$^3$]')
for nn in [0,1,2]:
ax=hca[nn]
pyic.plot_settings(ax, xlim=lon_reg_3, ylim=lat_reg_3)
ax.plot(lon_sec, lat_sec, transform=ccrs_proj, color='k')
ax.annotate('ageostrophic\nspiral eddy', xy=(-67.77, 37.04), xycoords='data',
xytext=(-67.11, 37.3), textcoords='data',
fontsize=10,
arrowprops=dict(arrowstyle="simple", facecolor='k', edgecolor='w', linewidth=1),
ha='center', va='center',
bbox=dict(ec='none', fc='w', alpha=0.9, boxstyle='square,pad=0.1')
)
for nn in [3,4,5]:
ax=hca[nn]
ax.set_ylim(depthc_sec.max(),0)
# ax.set_ylabel('depth [m]')
ax.set_xlabel('distance [km]')
ax.plot(dist_sec/1e3, mld, color='k')
ax.text(0.5, 0.1, 'pycnocline', transform=ax.transAxes, ha='center', va='center', fontsize=10,
bbox=dict(ec='none', fc='w', alpha=0.9, boxstyle='square,pad=0.3'),
)
ax.text(0.5, 0.8, 'mixed layer', transform=ax.transAxes, ha='center', va='center', fontsize=10,
bbox=dict(ec='none', fc='w', alpha=0.9, boxstyle='square,pad=0.3'),
)
hca[3].set_ylabel('depth [m]')
ax=hca[0]
ax.text(0.02,0.9, f'z = {depthc_sec[iz]}m', transform=ax.transAxes, bbox=dict(fc='w', ec='none'))
# hca[5].set_xlabel('longitude')
# savefig('closeup', format='pdf')
[59]:
Text(0.02, 0.9, 'z = 51.5m')
Vorticity vs. Rossby number#
[16]:
%%time
fpath = '/work/mh0033/m300602/proj_smtnatl/pp/pp_dissipation_all.nc'
f = Dataset(fpath, 'r')
diss_i = f.variables['diss_i'][:]
vort_i = f.variables['vort_i'][:]
f.close()
CPU times: user 363 ms, sys: 3.52 s, total: 3.88 s
Wall time: 4.79 s
[26]:
cell_area_reg2 = ds_tg.cell_area[ind_reg2].compute()
[27]:
diss_i.shape, cell_area_reg2.shape
[27]:
((112, 5615241), (5615241,))
[28]:
# hca, hcb = pyic.arrange_axes(1,1, plot_cb=False, asp=0.5,
# fig_size_fac=1.2, sharex=False, sharey=True)
# ii=-1
# ii+=1; ax=hca[ii]; cax=hcb[ii]
# ax.scatter(vort_i, np.ma.log10(-diss_i), s=2)
[29]:
# hca, hcb = pyic.arrange_axes(1,1, plot_cb=False, asp=0.5,
# fig_size_fac=1.2, sharex=False, sharey=True)
# ii=-1
# ii+=1; ax=hca[ii]; cax=hcb[ii]
# ax.scatter(np.abs(vort_sec), (-diss_sec), s=2)
Histrogram of dissipation vs. Rossby number#
[30]:
depthc_xr = xr.DataArray(depthc, dims=['depth'])
[37]:
%%time
vol = cell_area_reg2 * depthc_xr
vol = vol.transpose().compute().data
CPU times: user 563 ms, sys: 804 ms, total: 1.37 s
Wall time: 1.37 s
[67]:
%%time
bins = 101
hist, bin_edges = np.histogram(vort_i, bins=bins, weights=vol)
hdens = hist/hist.sum()/(bin_edges[1]-bin_edges[0])
# nps, bin_edges = np.histogram(vort_i, bins=bins)
CPU times: user 8.15 s, sys: 1.36 s, total: 9.5 s
Wall time: 9.53 s
[74]:
hca, hcb = pyic.arrange_axes(1, 1, plot_cb=False, asp=0.5, fig_size_fac=2)
ii=-1
ii+=1; ax=hca[ii]; cax=hcb[ii]
ax.plot(0.5*(bin_edges[1:]+bin_edges[:-1]), hdens, label='reg2', c='0.7')
# plot_histogram(ax, ro_r3, label='region 2', c='tab:blue')
# plot_histogram(ax, ro_r4, label='region 3', c='tab:orange')
ax.set_xlim(-2,3)
ax.grid(True)
ax.set_title('histogramm of local Rossby')
ax.set_xlabel('Rossby number')
ax.set_ylabel('points in interval\n / all points / interval size')
ax.legend()
# # ax.axvline(-1, color='k')
# # ax.axvline(1, color='k')
# ax.axvline(0, color='k')
# plt.axvspan(-2, -1, color='b', alpha=0.1, lw=0)
# plt.axvspan(1, 3, color='b', alpha=0.1, lw=0)
# ax.text(-1.5, 0.3, 'symmetrically\nunstable', ha='center', va='center')
# ax.text(-1.5, 0.7, 'ageostrophic\nregime', ha='center', va='center')
# ax.text(1.5, 0.7, 'ageostrophic\nregime', ha='center', va='center')
# # ax.annotate('local max', xy=(2, 1), xytext=(3, 1.5),
# # arrowprops=dict(facecolor='black', shrink=0.05),
# # )
# # savefig('histogram', format='pdf')
[74]:
<matplotlib.legend.Legend at 0x7ffcc07353a0>
[75]:
%%time
bins = 40
# hist, bin_edges = np.histogram(vort_sec, bins=bins, weights=-diss_sec)
# nps, bin_edges = np.histogram(vort_sec, bins=bins)
hist, bin_edges = np.histogram(vort_i, bins=bins, weights=-diss_i*vol)
# nps, bin_edges = np.histogram(vort_i, bins=bins)
CPU times: user 8.94 s, sys: 1.91 s, total: 10.8 s
Wall time: 10.9 s
[76]:
# def plot_histogram(ax, roi, label='asdf', c='tab:blue'):
# roi2 = roi.copy()
# roi2 = roi2[roi2.mask==False]
# hist, bins = np.histogram(roi2, bins=np.linspace(-4,4,101))
# hdens = hist/hist.sum()/(bins[1]-bins[0])
# ax.plot(0.5*(bins[1:]+bins[:-1]), hdens, marker=None, label=label, color=c)
# print((hdens*(bins[1]-bins[0])).sum())
# return
[77]:
hca, hcb = pyic.arrange_axes(1,1, plot_cb=False, asp=0.3,
fig_size_fac=1., sharex=False, sharey=True)
ii=-1
ii+=1; ax=hca[ii]; cax=hcb[ii]
ax.plot(0.5*(bin_edges[1:]+bin_edges[:-1]), np.log10(hist))
ax.set_xlabel('local Rossby number')
ax.set_ylabel('mean dissipation [m$^2$/s$^3$]')
ax.set_title('mean kin. energy dissipation within Rossby number interval')
ax.grid(True)
ax.set_xlim(-10,10)
ax.axvline(0, color='k')
plt.axvspan(-2, -1, color='b', alpha=0.1, lw=0)
plt.axvspan(1, 3, color='b', alpha=0.1, lw=0)
ax.text(0.1, 0.3, 'symmetrically\nunstable', ha='center', va='center', transform=ax.transAxes)
ax.text(0.1, 0.85, 'ageostrophic\nregime', ha='center', va='center', transform=ax.transAxes)
ax.text(0.7, 0.85, 'ageostrophic\nregime', ha='center', va='center', transform=ax.transAxes)
[77]:
Text(0.7, 0.85, 'ageostrophic\nregime')
PDE of dissipation#
[41]:
# clon.shape
[41]:
(59799625,)
[131]:
# # lon_reg_2 = [-75, -55]
# # lat_reg_2 = [33, 43]
# lon_reg = lon_reg_2
# lat_reg = lat_reg_2
# clon_reg, clat_reg, vertex_of_cell_reg, edge_of_cell_reg, ind_reg = pyic.crop_tripolar_grid(lon_reg, lat_reg,
# clon, clat, vertex_of_cell, edge_of_cell)
# Tri = matplotlib.tri.Triangulation(vlon, vlat, triangles=vertex_of_cell_reg)
[132]:
# clon_reg.shape, ind_reg.shape
[132]:
((5615241,), (5615241,))
[133]:
# ind_mask = ( (clon_reg>lon_reg_3[0])
# & (clon_reg<=lon_reg_3[1])
# & (clat_reg>lat_reg_3[0])
# & (clat_reg<=lat_reg_3[1]) )
# Tri.set_mask(ind_mask==False)
[134]:
# fpath = '/work/mh0033/m300602/proj_smtnatl/pp/pp_dissipation_all.nc'
# # f = Dataset(fpath, 'r')
# # diss_i = f.variables['diss_i'][:]
# # f.close()
[135]:
# nz, nc = diss_i.shape
[136]:
# diss_i.shape, ind_reg.shape
[136]:
((112, 5615241), (5615241,))
[42]:
izs = pyic.indfind([2.5, 50., 500., 3000.], depthc)
# izs = [0]
depthc[izs]
[42]:
array([2.50000e+00, 4.85000e+01, 4.96500e+02, 3.06375e+03])
[43]:
import scipy
[50]:
%%time
bins = 120
bin_edges = np.zeros((len(izs), bins+1))
hist = np.zeros((len(izs), bins))
norm_dist = np.zeros((len(izs), bins))
# mean = np.zeros((len(izs), bins))
# var = np.zeros((len(izs), bins))
# skew = np.zeros((len(izs), bins))
# kurt = np.zeros((len(izs), bins))
for kk in range(len(izs)):
data = -diss_i[izs[kk],:]#*vol[izs[kk],:]
data = data[data!=0.]
hist[kk,:], bin_edges[kk,:] = np.histogram(np.log10(data), bins=bins, density=True)
x = 0.5*(bin_edges[kk,1:]+bin_edges[kk,:-1])
mu = np.log10(data).mean()
sigma = np.log10(data).std()
norm_dist[kk,:] = 1/(sigma * np.sqrt(2 * np.pi)) * np.exp(-(x - mu)**2 / (2 * sigma**2))
# norm_dist[kk,:] = np.random.normal(np.log10(data).mean(), np.log10(data).std(), bins)
# mean[kk,:] = scipy.stats.norm.stats()
CPU times: user 766 ms, sys: 90.2 ms, total: 856 ms
Wall time: 862 ms
[52]:
hca, hcb = pyic.arrange_axes(1,1, plot_cb=False, asp=0.8,
fig_size_fac=2., sharex=False, sharey=True)
ii=-1
cols = plt.rcParams['axes.prop_cycle'].by_key()['color']
ii+=1; ax=hca[ii]; cax=hcb[ii]
for kk in range(len(izs)):
ax.plot(0.5*(bin_edges[kk,1:]+bin_edges[kk,:-1]), hist[kk,:], label=f'{depthc[izs[kk]]} m', marker='x', linestyle='', color=cols[kk])
ax.plot(0.5*(bin_edges[kk,1:]+bin_edges[kk,:-1]), norm_dist[kk,:], color=cols[kk])
ax.set_xlabel('log$_{10}$(mean dissipation / (m$^2$ s$^{-3}$))')
ax.set_ylabel('PDF')
ax.set_title('kin. energy dissipation')
ax.grid(True)
ax.set_xlim(-12,-4)
ax.legend()
[52]:
<matplotlib.legend.Legend at 0x7ffcc1a9d100>
Hor. averaged dissipation#
[208]:
ccrs_proj = ccrs.PlateCarree()
hca, hcb = pyic.arrange_axes(2, 2, plot_cb=True, asp=0.5, projection=ccrs_proj, fig_size_fac=1.2)
ii=-1
clim = 1e-6
for kk in [0, 10, 40, 63]:
# for kk in [0]:
print(kk)
ii+=1; ax=hca[ii]; cax=hcb[ii]
hm2 = pyic.shade(Tri, diss_i[kk,:], ax=ax, cax=cax, clim=clim,
transform=ccrs_proj, rasterized=False)
ax.set_title(f'kin. energy dissipation at {depthc[kk]}m [m$^2$/s$^3$]')
for ax in hca:
pyic.plot_settings(ax, xlim=lon_reg_3, ylim=lat_reg_3)
# savefig('dissipation_diff_depth', format='png')
print('done')
0
10
40
63
done
[209]:
diss_hint = (diss_i*cell_area_p[np.newaxis,ind_reg]).sum(axis=1)/cell_area_p[ind_reg].sum()
[211]:
hca, hcb = pyic.arrange_axes(1, 1, plot_cb=False, asp=1.2, fig_size_fac=2, ylabel='depth [m]')
ii=-1
ii+=1; ax=hca[ii]; cax=hcb[ii]
ax.plot(diss_hint, depthc[:nz])
ax.set_title('hor. ave. dissipation [m$^2$/s$^3$]')
ax.set_ylim(depthc[:nz].max(),0)
ax.set_xlim(diss_hint.min(), 0)
# savefig('dissipation_have', format='png')
ax.grid(True)
[ ]: