In [1]:
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr

Details of the reference ($\beta$-plan) experiment

  • MPI with NP_XI$\times$NP_ETA=5$\times$5
  • 60x50x10 grid points
  • $\Delta x \times \Delta y \times \Delta z$ = 60km $\times$ 56km $\times$ 500m --> 3600km $\times$ 2800km $\times$ 5000m
  • $\Delta t_{bci}$ = 9600 sec
  • $\Delta t_{btp}$ = $\frac{\Delta t_{bci}}{65}$ $\approx 148$sec
  • $f_0$=10$^{-4}$ s$^{-1}$;
  • $\beta$=2$\times$10$^{-11}$ (m s)$^{-1}$
  • bottom drag = 3$\times$10$^{-4}$ m$^2$ s$^{-1}$
  • no-slip bdy cd (gamma2=-1) -- in croco.in
  • VISC2 = DIFF2 = 1000 m$^2$s$^{-1}$
  • 20 years run with 30 model state output
  • PV and (2D) KE budget as snapshots every 30 days and averaged every 5 years
  • $\frac{\Delta t_{btp}}{\Delta x} \sqrt{gH} \leqslant 0.89$   --->>>   $\Delta t_{btp} \leqslant 225 sec$
  • $\Delta t_{btp} \leqslant \frac{\Delta t_{btp}}{65} \leqslant 14627 sec$

Load data

In [2]:
saveF = False

#-- directories --
dir_in = '/home2/datawork/qjamet/CROCO_Inria/RUNS/BASIN'
dir_fig = '/home2/datahome/qjamet/Figures/BASIN'

rrun = ['1001', '1002', '1002b', '1003', '1004', '1005']
nrun = len(rrun)
nname = [r'Ref ($\beta$-plan)', r'f-plan', 'no rotation' ,'free-slip ; bottom drag=3e-5 m2/s', 'no UV_ADV', 'High-Res']

dfs = {}
for iii in range(nrun):
    print("-- load file for run: ", rrun[iii], end="\r")
    dfs[iii]  = xr.open_dataset("%s/run_%s/basin_his.nc" % (dir_in, rrun[iii]))
    
#-- vorticity budgets --
dfs_vrt = {}
for iii in range(nrun):
    print("-- load file for run: ", rrun[iii], end="\r")
    dfs_vrt[iii]  = xr.open_dataset("%s/run_%s/basin_diags_vrt_avg.nc" % (dir_in, rrun[iii]))

-- load file for run:  1005b

First look

In [7]:
ttt = -1
fig, ax = plt.subplots(figsize=(12,15), nrows=3, ncols=2)
#
irun=0
for ip in range(3):
    for jp in range(2):
        cs=ax[ip,jp].contourf(dfs[irun].x_rho, dfs[irun].y_rho, dfs[irun].zeta[ttt, :, :], \
                              20, cmap='RdBu_r', levels=np.arange(-0.5, 0.55, 0.05), extend='both')
        ax[ip,jp].set(title=nname[irun], xlabel='xi-dir [km]', ylabel='eta-dir [km]')
        fig.colorbar(ax=ax[ip, jp], mappable=cs)
        irun+=1
    
if saveF:
    figN = str("BASIN_zeta_f-plan_beta-plan")
    fig1.savefig(str("%s/%s.png" % (dir_fig, figN)), bbox_inches='tight')
    fig1.savefig(str("%s/%s.pdf" % (dir_fig, figN)), bbox_inches='tight')

Barotropic vorticity (2d, vertically integrated) budgets

From https://www.jgula.fr/Croco/diagnostics_croco.pdf

vrt_rate = vrt_xadv + vrt_yadv + vrt_prsgrd + vrt_cor + vrt_vmix + vrt_hmix + vrt_nudg

  • vrt_rate = rate of change of barotropic vorticity
  • vrt_xadv = contribution of advection + implicit dissipation along xi-axis + grid curvature terms (see CURVGRID)
  • vrt_hdiff = implicit dissipation along xi- and eta- axis (already included in vrt xadv+ vrt yadv)
  • vrt_cor = planetary vorticity advection
  • vrt_prsgrd = bottom Pressure torque
  • vrt_hmix = contribution of Horizontal diffusion (explicit)
  • vrt_vmix = contribution of Vertical mixing = vrt Wind + vrt Drag
  • vrt_nudg = contribution of Nudging, restoring, boundary conditions, etc.
  • vrt_Wind = Wind stress curl (already included in vrt_vmix)
  • vrt_Drag = Bottom drag curl (already included in vrt_vmix)
In [33]:
dfs_vrt[0]
Out[33]:
<xarray.Dataset>
Dimensions:     (time: 4, auxil: 4, eta_v: 51, xi_u: 61)
Coordinates:
  * time        (time) float64 1.555e+08 3.11e+08 4.666e+08 6.221e+08
Dimensions without coordinates: auxil, eta_v, xi_u
Data variables: (12/13)
    time_step   (time, auxil) int32 16201 0 61 1 32401 0 ... 181 3 64801 0 241 4
    scrum_time  (time) float64 7.776e+07 2.333e+08 3.888e+08 5.443e+08
    vrt_xadv    (time, eta_v, xi_u) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    vrt_yadv    (time, eta_v, xi_u) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    vrt_hdiff   (time, eta_v, xi_u) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    vrt_cor     (time, eta_v, xi_u) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    ...          ...
    vrt_hmix    (time, eta_v, xi_u) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    vrt_vmix    (time, eta_v, xi_u) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    vrt_rate    (time, eta_v, xi_u) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    vrt_nudg    (time, eta_v, xi_u) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    vrt_wind    (time, eta_v, xi_u) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    vrt_drag    (time, eta_v, xi_u) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
Attributes: (12/74)
    type:           CROCO history file
    title:          \tBasin Example
    date:           
    rst_file:       basin_rst.nc
    his_file:       basin_his.nc
    avg_file:       basin_avg.nc
    ...             ...
    levbfrc:        1
    levbfrc_expl:   Deepest level for body-force stress
    gamma2:         -1.0
    gamma2_expl:    Slipperiness parameter
    SRCS:           main.F step.F read_inp.F timers_roms.F init_scalars.F ini...
    CPP-options:    BASIN MPI SOLVE3D UV_COR UV_ADV UV_HADV_UP3 UV_VIS2 UV_VA...

First, check that PV budget is closed

In [34]:
ttt = -1
#- beta-plane -
vrt_rate = dfs_vrt[0].vrt_rate[ttt, :, :]
vrt_rhs  = dfs_vrt[0].vrt_xadv[ttt, :, :] + dfs_vrt[0].vrt_yadv[ttt, :, :] + dfs_vrt[0].vrt_prsgrd[ttt, :, :] + \
            dfs_vrt[0].vrt_cor[ttt, :, :] + dfs_vrt[0].vrt_vmix[ttt, :, :] + dfs_vrt[0].vrt_hmix[ttt, :, :]
In [35]:
jjj = 20
scaleF = 1e3
fig1 = plt.figure(figsize=(12, 10))
fig1.clf()
#
ax1 = fig1.add_subplot(2, 2, 1)
cs1 = ax1.contourf(vrt_rate, 20, cmap='RdBu_r')
plt.colorbar(cs1)
ax1.plot(np.arange(daVrt.dims['xi_u']), jjj*np.ones(daVrt.dims['xi_u']), 'k.')
ax1.set_title(r'1/ $\partial_t (\overline{PV}^{z})$')
#
ax2 = fig1.add_subplot(2, 2, 2)
cs2 = ax2.contourf(vrt_rhs, 20, cmap='RdBu_r')
plt.colorbar(cs2)
ax2.set_title(r'2/ $\sum RHS$')
#
ax3 = fig1.add_subplot(2, 2, 3)
cs3 = ax3.contourf(vrt_rate-vrt_rhs, 20, cmap='RdBu_r')
plt.colorbar(cs3)
ax3.set_title(r'Residual (1/-2/)')
#
ax4 = fig1.add_subplot(2, 2, 4)
p0 = ax4.plot(vrt_rate[jjj, :], 'k', lw=3)
p1 = ax4.plot(vrt_rhs[jjj, :], 'r.')
p2 = ax4.plot((vrt_rate[jjj, :]-vrt_rhs[jjj, :])*scaleF, 'g-')
ax4.legend((p0[0], p1[0], p2[0]), \
           (r"$\partial_t \overline{PV}^{z}$",\
            r"$\sum RHS$", \
            str(r"Residual ($\times$ %01.00e)" % (1/scaleF) )))
ax4.set_title(r"Residual -- X section")
plt.grid()

PV budegt for the $\beta$-plan case

In [43]:
irun = 0
#- beta-plane -
llev1 = np.arange(-1.0, 1.1, 0.1)*5e-10
#
fig1 = plt.figure(figsize=(15, 10))
fig1.clf()
fig1.suptitle(nname[irun], fontsize=30)
#- rate -
ax1 = fig1.add_subplot(2, 3, 1)
cs1 = ax1.contourf(dfs_vrt[irun].vrt_rate[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax1.set_title(r'Rate')
#- wind stress -
ax2 = fig1.add_subplot(2, 3, 2)
cs2 = ax2.contourf(dfs_vrt[irun].vrt_wind[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax2.set_title(r'Wind stress')
#- planetary vorticity advection -
ax3 = fig1.add_subplot(2, 3, 3)
cs3 = ax3.contourf(dfs_vrt[irun].vrt_cor[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax3.set_title(r'Planetary vorticity advection')
#- NL advection -
ax4 = fig1.add_subplot(2, 3, 4)
cs4 = ax4.contourf(dfs_vrt[irun].vrt_xadv[ttt, :, :]+dfs_vrt[irun].vrt_yadv[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax4.set_title(r'NL advection')
#- horizontal mixing -
ax5 = fig1.add_subplot(2, 3, 5)
cs5 = ax5.contourf(dfs_vrt[irun].vrt_hmix[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax5.set_title(r'Hz dissipation')
#- bottom drag -
ax6 = fig1.add_subplot(2, 3, 6)
cs6 = ax6.contourf(dfs_vrt[irun].vrt_drag[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax6.set_title(r'Bottom drag')
#
cbax1 = fig1.add_axes([0.93, 0.25, 0.01, 0.5])
cb1 = fig1.colorbar(cs1, ax=ax1, orientation='vertical', cax=cbax1)

PV budegt for the $f$-plan case

In [47]:
irun = 1
#- beta-plane -
llev1 = np.arange(-1.0, 1.1, 0.1)*5e-10
#
fig1 = plt.figure(figsize=(15, 10))
fig1.clf()
fig1.suptitle(nname[irun], fontsize=30)
#- rate -
ax1 = fig1.add_subplot(2, 3, 1)
cs1 = ax1.contourf(dfs_vrt[irun].vrt_rate[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax1.set_title(r'Rate')
#- wind stress -
ax2 = fig1.add_subplot(2, 3, 2)
cs2 = ax2.contourf(dfs_vrt[irun].vrt_wind[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax2.set_title(r'Wind stress')
#- planetary vorticity advection -
ax3 = fig1.add_subplot(2, 3, 3)
cs3 = ax3.contourf(dfs_vrt[irun].vrt_cor[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax3.set_title(r'Planetary vorticity advection')
#- NL advection -
ax4 = fig1.add_subplot(2, 3, 4)
cs4 = ax4.contourf(dfs_vrt[irun].vrt_xadv[ttt, :, :]+dfs_vrt[irun].vrt_yadv[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax4.set_title(r'NL advection')
#- horizontal mixing -
ax5 = fig1.add_subplot(2, 3, 5)
cs5 = ax5.contourf(dfs_vrt[irun].vrt_hmix[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax5.set_title(r'Hz dissipation')
#- bottom drag -
ax6 = fig1.add_subplot(2, 3, 6)
cs6 = ax6.contourf(dfs_vrt[irun].vrt_drag[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax6.set_title(r'Bottom drag')
#
cbax1 = fig1.add_axes([0.93, 0.25, 0.01, 0.5])
cb1 = fig1.colorbar(cs1, ax=ax1, orientation='vertical', cax=cbax1)

PV budegt with reduced bottom drag and free-slip bdy cd (instead of no-slip)

In [4]:
irun = 2
#- beta-plane -
llev1 = np.arange(-1.0, 1.1, 0.1)*5e-10
#
fig1 = plt.figure(figsize=(15, 10))
fig1.clf()
fig1.suptitle(nname[irun], fontsize=30)
#- rate -
ax1 = fig1.add_subplot(2, 3, 1)
cs1 = ax1.contourf(dfs_vrt[irun].vrt_rate[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax1.set_title(r'Rate')
#- wind stress -
ax2 = fig1.add_subplot(2, 3, 2)
cs2 = ax2.contourf(dfs_vrt[irun].vrt_wind[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax2.set_title(r'Wind stress')
#- planetary vorticity advection -
ax3 = fig1.add_subplot(2, 3, 3)
cs3 = ax3.contourf(dfs_vrt[irun].vrt_cor[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax3.set_title(r'Planetary vorticity advection')
#- NL advection -
ax4 = fig1.add_subplot(2, 3, 4)
cs4 = ax4.contourf(dfs_vrt[irun].vrt_xadv[ttt, :, :]+dfs_vrt[irun].vrt_yadv[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax4.set_title(r'NL advection')
#- horizontal mixing -
ax5 = fig1.add_subplot(2, 3, 5)
cs5 = ax5.contourf(dfs_vrt[irun].vrt_hmix[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax5.set_title(r'Hz dissipation')
#- bottom drag -
ax6 = fig1.add_subplot(2, 3, 6)
cs6 = ax6.contourf(dfs_vrt[irun].vrt_drag[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax6.set_title(r'Bottom drag')
#
cbax1 = fig1.add_axes([0.93, 0.25, 0.01, 0.5])
cb1 = fig1.colorbar(cs1, ax=ax1, orientation='vertical', cax=cbax1)

PV budegt with no NON-LINEAR term

In [5]:
irun = 3
#- beta-plane -
llev1 = np.arange(-1.0, 1.1, 0.1)*5e-10
#
fig1 = plt.figure(figsize=(15, 10))
fig1.clf()
fig1.suptitle(nname[irun], fontsize=30)
#- rate -
ax1 = fig1.add_subplot(2, 3, 1)
cs1 = ax1.contourf(dfs_vrt[irun].vrt_rate[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax1.set_title(r'Rate')
#- wind stress -
ax2 = fig1.add_subplot(2, 3, 2)
cs2 = ax2.contourf(dfs_vrt[irun].vrt_wind[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax2.set_title(r'Wind stress')
#- planetary vorticity advection -
ax3 = fig1.add_subplot(2, 3, 3)
cs3 = ax3.contourf(dfs_vrt[irun].vrt_cor[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax3.set_title(r'Planetary vorticity advection')
#- NL advection -
ax4 = fig1.add_subplot(2, 3, 4)
cs4 = ax4.contourf(dfs_vrt[irun].vrt_xadv[ttt, :, :]+dfs_vrt[irun].vrt_yadv[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax4.set_title(r'NL advection')
#- horizontal mixing -
ax5 = fig1.add_subplot(2, 3, 5)
cs5 = ax5.contourf(dfs_vrt[irun].vrt_hmix[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax5.set_title(r'Hz dissipation')
#- bottom drag -
ax6 = fig1.add_subplot(2, 3, 6)
cs6 = ax6.contourf(dfs_vrt[irun].vrt_drag[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax6.set_title(r'Bottom drag')
#
cbax1 = fig1.add_axes([0.93, 0.25, 0.01, 0.5])
cb1 = fig1.colorbar(cs1, ax=ax1, orientation='vertical', cax=cbax1)

PV budegt with increased resolution and decreased visc/diff

  • maximum btp time step: $\Delta t_{btp} \leqslant 225 sec$
  • with 65 btp time step in 1 bci time-step $\Delta t_{bci} \leqslant 7390 sec$
  • reduced diss and diffus to 100 m$^{2}$ s$^{-1}$
In [6]:
irun = 4
#- 
llev1 = np.arange(-1.0, 1.1, 0.1)*5e-10
#
fig1 = plt.figure(figsize=(15, 10))
fig1.clf()
fig1.suptitle(nname[irun], fontsize=30)
#- rate -
ax1 = fig1.add_subplot(2, 3, 1)
cs1 = ax1.contourf(dfs_vrt[irun].vrt_rate[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax1.set_title(r'Rate')
#- wind stress -
ax2 = fig1.add_subplot(2, 3, 2)
cs2 = ax2.contourf(dfs_vrt[irun].vrt_wind[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax2.set_title(r'Wind stress')
#- planetary vorticity advection -
ax3 = fig1.add_subplot(2, 3, 3)
cs3 = ax3.contourf(dfs_vrt[irun].vrt_cor[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax3.set_title(r'Planetary vorticity advection')
#- NL advection -
ax4 = fig1.add_subplot(2, 3, 4)
cs4 = ax4.contourf(dfs_vrt[irun].vrt_xadv[ttt, :, :]+dfs_vrt[irun].vrt_yadv[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax4.set_title(r'NL advection')
#- horizontal mixing -
ax5 = fig1.add_subplot(2, 3, 5)
cs5 = ax5.contourf(dfs_vrt[irun].vrt_hmix[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax5.set_title(r'Hz dissipation')
#- bottom drag -
ax6 = fig1.add_subplot(2, 3, 6)
cs6 = ax6.contourf(dfs_vrt[irun].vrt_drag[ttt, :, :], cmap='RdBu_r', levels=llev1, extend='both')
ax6.set_title(r'Bottom drag')
#
cbax1 = fig1.add_axes([0.93, 0.25, 0.01, 0.5])
cb1 = fig1.colorbar(cs1, ax=ax1, orientation='vertical', cax=cbax1)

Conclusion

beta-plane

  • Leading order balance between wind forcing and planetary vorticity advection in the interior ocean
  • Particular dynamics along the western boundary

f-plane

  • Leading order balance between wind forcing and bottom drag in the interior ocean
  • Particular dynamics along boundaries

no UV_ADV

  • This is the test case of Stommel (1948)