smtnatl_overview: Overview of SMT-NATL simulation

Analysis of SMT-NATL experiment for SAB poster 2020. Do not change anymore.

# 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

%%javascript
IPython.notebook.kernel.execute('nb_name = "' + IPython.notebook.notebook_name + '"')

nb_name

'smtnatl_overview.ipynb'

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
tb
sys
scipy
matplotlib
cartopy
xarray
done loading
IconData
plotting
view
calc
tb
IconData
plotting
view
quickplots
quickplots

import xarray as xr
import dask
import dask.array as da
from dask.distributed import Client, LocalCluster
import multiprocessing

def savefig(txt, dpi=150, format='pdf'):
    fbase = nb_name.split('/')[-1][:-6]
    fpath = f'{path_fig}{fbase}_{txt}.{format}'
    plt.savefig(fpath, dpi=dpi)
    print('saved figure: %s' % (fpath))
    return

Define simulation

run = 'ngSMT_tke'
savefig = False
path_fig = '../pics/'
nnf=0

gname = 'smt'
lev = 'L128'

path_data = f'/mnt/lustre01/work/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/mpim/m300602/work/icon/grids/smt/smt_tgrid.nc'
fpath_Tri = '/mnt/lustre01/work/mh0033/m300602/tmp/Tri.pkl'

path_grid     = f'/mnt/lustre01/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 = '/mnt/lustre01/work/mh0033/m300602/proj_vmix/icon/icon_ckdtree/rectgrids/smt_res0.02_180W-180E_90S-90N.npz'

# 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]

(16, 51.5)

Load the grid

%%time
IcD = pyic.IconData(
                    fname        = f'{run}_vort_f_50m_*.nc',
                    path_data    = '/mnt/lustre01/work/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 /mnt/lustre01/work/mh0033/m300602/icon/grids/smt/ckdtree/sections/. :::
::: Warning: no section found.:::
-v-: list of variables and time steps
CPU times: user 228 ms, sys: 334 ms, total: 562 ms
Wall time: 3.16 s

%%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 3.93 s, sys: 13.7 s, total: 17.6 s
Wall time: 19 s

res = np.sqrt(cell_area_p)

%%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 3.54 s, sys: 8.26 s, total: 11.8 s
Wall time: 12.5 s

lon_reg_3

[-68.5, -66.5]

%%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 609 ms, sys: 519 ms, total: 1.13 s
Wall time: 1.19 s

Plot resolution

fpath_ckdtree_res = '/mnt/lustre01/work/mh0033/m300602/proj_vmix/icon/icon_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', radius_of_influence=1.).reshape(lat_010deg.size, lon_010deg.size)
res_010deg[res_010deg==0.] = np.ma.masked

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()]
hm1 = pyic.shade(lon_010deg, lat_010deg, res_010deg/1e3,
                 ax=ax, cax=cax, clim=clim, contfs=contfs,
                 transform=ccrs_proj, rasterized=False,
                 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')

saved figure: ../pics/smtnatl_overview_resolution.pdf

Load the data

# --- 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

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]

numpy.datetime64('2010-03-09T01:00:00')

vort = vort[ind_reg]

# 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)

(34.0, 36.0)

Overview maps

%%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 588 ms, sys: 511 ms, total: 1.1 s
Wall time: 1.1 s

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)

# --- 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]

numpy.datetime64('2010-03-09T01:00:00')

# --- 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', radius_of_influence=1.).reshape(lat_002deg.size, lon_002deg.size)
ro_002deg[ro_002deg==0.] = np.ma.masked

# --- 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]

import maps_icon_smt_vort as smt

# 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)

%%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_overview_main.jpg.
Done with file ../pics/smtnatl_overview_main.jpg.
CPU times: user 32.1 s, sys: 2.01 s, total: 34.1 s
Wall time: 25.5 s

Calculate dissipation

vertex_of_cell.shape, vertex_of_cell.shape

((59799625, 3), (59799625, 3))

# edges_of_vertex = edges_of_vertex[vertex_of_cell,:]
# edge_orientation = edge_orientation[vertex_of_cell,:]

%%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 7.98 s, sys: 6.86 s, total: 14.8 s
Wall time: 14.8 s

# --- 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)
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

%%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 9.56 s, sys: 5.77 s, total: 15.3 s
Wall time: 15.3 s

A4 = 0.04 * 500**3
diss *= A4

%%time
diss_i = diss[edge_of_cell].sum(axis=1)/3.

CPU times: user 3.38 s, sys: 1.31 s, total: 4.69 s
Wall time: 4.68 s

diss_ireg = diss_i[ind_reg]
divv_ireg = div_v[ind_reg]

curl_v.shape, vlon.shape, clon.shape, diss_i.shape

((30008985,), (30008985,), (59799625,), (59799625,))

%%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 1.82 s, sys: 860 ms, total: 2.68 s
Wall time: 2.67 s

%%time
ro_selfi = ro_self[vertex_of_cell].sum(axis=1)/3.

CPU times: user 3.29 s, sys: 1.4 s, total: 4.69 s
Wall time: 4.68 s

# fpath_ckdtree = '/mnt/lustre01/work/mh0033/m300602/icon/grids/smt/new_ckdtree/rectgrids/smt_res0.02_180W-180E_90S-90N.npz'

%%time
diss_002deg = pyic.apply_ckdtree(diss_i, fpath_ckdtree, coordinates='clat clon', radius_of_influence=1.).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', radius_of_influence=1.).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', 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]

CPU times: user 14 s, sys: 7.87 s, total: 21.8 s
Wall time: 21.8 s

# 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]

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)

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')

saved figure: ../pics/smtnatl_overview_dissipation.pdf

Histogram of Ro

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)

1.0
1.0

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

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

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
saved figure: ../pics/smtnatl_overview_histogram.pdf

Remap velocities

grid_sphere_radius = 6.371229e6

%%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 13 s, sys: 5.78 s, total: 18.8 s
Wall time: 18.7 s

%%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 14.2 s, sys: 7.12 s, total: 21.3 s
Wall time: 21.3 s

edge2cell_coeff_cc.shape, ve.shape

((59799625, 3, 3), (89813639,))

edge2cell_coeff_cc.shape

(59799625, 3, 3)

%%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 22.3 s, sys: 5.53 s, total: 27.8 s
Wall time: 27.8 s

%%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 8.84 s, sys: 5.82 s, total: 14.7 s
Wall time: 14.7 s

uo_002deg = pyic.apply_ckdtree(uo, fpath_ckdtree, coordinates='clat clon', radius_of_influence=1.).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', radius_of_influence=1.).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')

Region plots

# 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)

path_scratch = '/scratch/m/m300602/tmp/'

Load temperature

# --- 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')

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'

path_data = '/mnt/lustre01/work/mh0287/users/leonidas/icon/ngSMT/results/2010-03/'

levs = np.concatenate((np.arange(0,48,2),np.arange(48,48+40,1))).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]]

(array([ 0,  2,  4,  6,  8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32,
        34, 36, 38, 40, 42, 44, 46, 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]),
 64,
 695.25)

%%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/64; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T001_sp
kk = 1/64; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T003_sp
kk = 2/64; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T005_sp
kk = 3/64; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T007_sp
kk = 4/64; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T009_sp
kk = 5/64; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T011_sp
kk = 6/64; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T013_sp
kk = 7/64; fname = ngSMT_tke_T_S_sp_001-016_20100309T010000Z.nc; var = T015_sp
kk = 8/64; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T017_sp
kk = 9/64; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T019_sp
kk = 10/64; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T021_sp
kk = 11/64; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T023_sp
kk = 12/64; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T025_sp
kk = 13/64; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T027_sp
kk = 14/64; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T029_sp
kk = 15/64; fname = ngSMT_tke_T_S_sp_017-032_20100309T010000Z.nc; var = T031_sp
kk = 16/64; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T033_sp
kk = 17/64; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T035_sp
kk = 18/64; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T037_sp
kk = 19/64; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T039_sp
kk = 20/64; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T041_sp
kk = 21/64; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T043_sp
kk = 22/64; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T045_sp
kk = 23/64; fname = ngSMT_tke_T_S_sp_033-048_20100309T010000Z.nc; var = T047_sp
kk = 24/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T049_sp
kk = 25/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T050_sp
kk = 26/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T051_sp
kk = 27/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T052_sp
kk = 28/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T053_sp
kk = 29/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T054_sp
kk = 30/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T055_sp
kk = 31/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T056_sp
kk = 32/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T057_sp
kk = 33/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T058_sp
kk = 34/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T059_sp
kk = 35/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T060_sp
kk = 36/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T061_sp
kk = 37/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T062_sp
kk = 38/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T063_sp
kk = 39/64; fname = ngSMT_tke_T_S_sp_049-064_20100309T010000Z.nc; var = T064_sp
kk = 40/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T065_sp
kk = 41/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T066_sp
kk = 42/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T067_sp
kk = 43/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T068_sp
kk = 44/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T069_sp
kk = 45/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T070_sp
kk = 46/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T071_sp
kk = 47/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T072_sp
kk = 48/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T073_sp
kk = 49/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T074_sp
kk = 50/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T075_sp
kk = 51/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T076_sp
kk = 52/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T077_sp
kk = 53/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T078_sp
kk = 54/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T079_sp
kk = 55/64; fname = ngSMT_tke_T_S_sp_065-080_20100309T010000Z.nc; var = T080_sp
kk = 56/64; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T081_sp
kk = 57/64; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T082_sp
kk = 58/64; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T083_sp
kk = 59/64; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T084_sp
kk = 60/64; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T085_sp
kk = 61/64; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T086_sp
kk = 62/64; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T087_sp
kk = 63/64; fname = ngSMT_tke_T_S_sp_081-096_20100309T010000Z.nc; var = T088_sp
CPU times: user 45 s, sys: 29.6 s, total: 1min 14s
Wall time: 1min 28s

var

'T088_sp'

temp_reg.shape, ind_reg2.shape

((64, 5615241), (5615241,))

# temp_002deg = pyic.apply_ckdtree(temp, fpath_ckdtree, coordinates='clat clon', radius_of_influence=1.).reshape(lat_002deg.size, lon_002deg.size)
# temp_002deg[temp_002deg==0.] = np.ma.masked

# i1 = (lon<lon_reg[0]).sum()
# i2 = (lon<lon_reg[1]).sum()
# j1 = (lat<lat_reg[0]).sum()
# j2 = (lat<lat_reg[1]).sum()

# temp_reg = temp_002deg[j1:j2,i1:i2]

Load other fields

fpath = '/scratch/m/m300602/tmp/pp_dissipation_new.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,:]

8 51.5

depthc_sec[:10], depthc[:10]

(masked_array(data=[ 2.5,  9.5, 15.5, 21.5, 27.5, 33.5, 39.5, 45.5, 51.5,
                    57.5],
              mask=False,
        fill_value=1e+20,
             dtype=float32),
 array([ 2.5,  6.5,  9.5, 12.5, 15.5, 18.5, 21.5, 24.5, 27.5, 30.5]))

f.variables.keys()

dict_keys(['depthc', 'diss_i', 'vort_i', 'div_v', 'uo', 'vo'])

ind_reg2.shape, uo.shape, temp_lev.shape

((5615241,), (5615241,), (5615241,))

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)

lon_reg

[-68.5, -66.5]

lon_o, lat_o = np.meshgrid(lon, lat)
lon_o = lon_o.flatten()
lat_o = lat_o.flatten()
lon_o.shape, lat_o.shape

((80000,), (80000,))

clon_reg.max(), clon_reg.min(), lon_o.max(), lon_o.min()

(-66.50003403874472, -68.49996742078073, -66.5, -68.5)

%%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 35.5 s, sys: 1.02 s, total: 36.5 s
Wall time: 36.5 s

# 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)

# 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)

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')

Load section data

%%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 34.6 s, sys: 789 ms, total: 35.3 s
Wall time: 35.3 s

%%time
kk = 0
fpath = '/scratch/m/m300602/tmp/pp_dissipation_new.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 889 ms, sys: 1.35 s, total: 2.23 s
Wall time: 2.23 s

# # 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)

Plotting

lon.shape, lat.shape, tempi.shape, vort_ii.shape

((400,), (200,), (200, 400), (200, 400))

depthc_sec

masked_array(data=[  2.5 ,   9.5 ,  15.5 ,  21.5 ,  27.5 ,  33.5 ,  39.5 ,
                    45.5 ,  51.5 ,  57.5 ,  63.5 ,  69.5 ,  75.5 ,  81.5 ,
                    87.5 ,  93.5 ,  99.5 , 105.5 , 111.5 , 117.5 , 123.5 ,
                   129.5 , 136.25, 143.25, 150.25, 154.  , 158.  , 162.  ,
                   166.  , 170.  , 174.25, 178.75, 183.25, 187.75, 192.5 ,
                   197.5 , 202.75, 208.5 , 214.75, 221.5 , 228.75, 236.5 ,
                   244.75, 254.  , 264.75, 277.  , 290.5 , 305.  , 320.25,
                   336.25, 353.  , 370.5 , 389.  , 408.5 , 429.  , 450.5 ,
                   473.  , 496.5 , 521.75, 549.75, 580.75, 615.  , 653.  ,
                   695.25],
             mask=False,
       fill_value=1e+20,
            dtype=float32)

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_sec, 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_sec, -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_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 ax in hca:
    ax.set_ylim(depthc_sec.max(),0)
hca[0].set_ylim(700,0)

(700.0, 0.0)

kmld = (temp_sec>(temp_sec[0:1,:]-0.5)).sum(axis=0)
mld = depthc_sec[kmld]

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,0.2e-6]
climD=[-9,-6]

# --- 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='Blues', logplot=True)
ax.set_title('log$_{10}$(kin. energy dissipation) [m$^2$/s$^3$]')

# --- sections
ii+=1; ax=hca[ii]; cax=hcb[ii]
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='Blues', 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')

saved figure: ../pics/smtnatl_overview_closeup.pdf

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_ii, np.ma.log10(-diss_ii), s=2)

<matplotlib.collections.PathCollection at 0x2b2d41e27550>

Hor. averaged dissipation

# 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)

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)

fpath = '/scratch/m/m300602/tmp/pp_dissipation_new.nc'
f = Dataset(fpath, 'r')
diss_i = f.variables['diss_i'][:]
f.close()

nz, nc = diss_i.shape

diss_i.shape, ind_reg.shape

((64, 5615241), (5615241,))

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

diss_hint = (diss_i*cell_area_p[np.newaxis,ind_reg]).sum(axis=1)/cell_area_p[ind_reg].sum()

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')

saved figure: ../pics/smtnatl_overview_dissipation_have.png