notebook.community

Edit and run

MOM6 diagnostics for KPP single column test case

Results from this notebook:

  1. Basic diagnostics of KPP boundary layer and prognostic fields.

Caveats regarding this notebook:

  1. This notebook is written for MOM6-examples/ocean_only/single_column_z.

Hopes for the use of this notebook:

  1. To provide a starting point to document single column model tests;
  2. To illustrate a self-contained iPython notebook of use for MOM6 analysis.

This iPython notebook was originally developed at NOAA/GFDL, and it is provided freely to the MOM6 community. GFDL scientists developing MOM6 make extensive use of Python for diagnostics. We solicit modifications/fixes that are useful to the MOM6 community.





In [2]:



    


import matplotlib.pyplot as plt                  # For plotting
import netCDF4                                   # For reading data
pylab.rcParams['figure.figsize'] = (17.0, 10.0)  # Large figures

path = '../single_column_z/'
deltaz = 'CM4'
fname_deltaz = '_CM4'
expt_name = 'single_column_z'

dpi=200



    











In [3]:



    


visc = netCDF4.Dataset(path+'visc.nc')
for v in visc.variables: print v,

# recall data layout is (t,z,y,x)

time = visc.variables['Time'][:]

# KPP boundary layer depth as function of time (metre)
h    = visc.variables['KPP_OBLdepth'][:,0,0]

# KPP tracer diffusivity as function of time and depth (m2/sec)
Kt   = visc.variables['KPP_Kheat'][:,:,0,0]

# KPP velocity diffusivity as function of time and depth (m2/sec)
Ku   = visc.variables['KPP_Kv'][:,:,0,0]

# surface (and penetrating) buoyancy flux, as used by [CVmix] KPP (m2/s3)
KPP_buoyFlux = visc.variables['KPP_buoyFlux'][:,:,0,0]

# Temperature tendency due to non-local transport of heat, as calculated by KPP (K/s)
KPP_dTdt = visc.variables['KPP_dTdt'][:,:,0,0]

# Saln tendency due to non-local transport of heat, as calculated by KPP (ppt/s)
KPP_dSdt = visc.variables['KPP_dSdt'][:,:,0,0]

# friction velocity used by KPP (m/s)
KPP_uStar = visc.variables['KPP_uStar'][:,0,0]

# Squared buoyancy frequency (1/s^2)
KPP_N2 = visc.variables['KPP_N2'][:,:,0,0]



    


















    






xh yh zi Time nv zl xq yq Kd_interface KPP_uStar KPP_buoyFlux KPP_QminusSW KPP_netSalt KPP_OBLdepth KPP_sigma KPP_BulkRi KPP_Ws KPP_N KPP_N2 KPP_Vt2 KPP_BulkUz2 KPP_BulkDrho KPP_Kd_in KPP_Kheat KPP_Ksalt KPP_Kv KPP_NLtransport_heat KPP_NLtransport_salt KPP_dTdt KPP_dSdt KPP_Tsurf KPP_Ssurf KPP_Usurf KPP_Vsurf MLD_003 average_T1 average_T2 average_DT Time_bounds

















In [4]:



    


prog  = netCDF4.Dataset(path+'prog.nc')
for v in prog.variables: print v,
   
# depth of cell interface     
zi = prog.variables['zi'][:]

# depth of cell center 
zl = prog.variables['zl'][:]

# zonal velocity as function of time and depth
u  = prog.variables['u'][:,:,0,0]

# zonal velocity as function of time and depth
v  = prog.variables['v'][:,:,0,0]

# temperature as function of time and depth 
temp = prog.variables['temp'][:,:,0,0]


# saln as function of time and depth 
saln = prog.variables['salt'][:,:,0,0]



    


















    






xq yh zl Time xh yq zi u v h e temp salt wd

















In [5]:



    


# thickness of cells 
dz    = zi[1:]-zi[:-1]
field = zi[1:]+dz[:]
#depth = np.ma.masked_array(data, mask=[data>400.])
depth = field 

#print dz
#print dz[1901:2000]
print zi

figure(1)
fig, ax = plt.subplots(figsize=(16, 10), dpi=dpi)
S = plt.plot(dz,-depth,'*',linewidth=3)
plt.ylim((-400,0))
plt.xlim((0, 50))
yticks = [-400, -350, -300, -250, -200, -150, -100, -50, 0]
ax.set_yticks(yticks)
#xticks = [0, .1, .2, .3, .4, .5]
#xticks = [0, .2, .4, .6, .8, 1.0, 1.2]
#xticks = [0, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260]
xticks = [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50]
ax.set_xticks(xticks)
#plt.xlim((0, 40)
ax.set_xticks(xticks)
plt.grid()
plt.xlabel('Delta z (m)',fontsize=18)
plt.ylabel('depth (m)',fontsize=18)
plt.title(r'$\Delta z$ for configuration '+deltaz,fontsize=20)
ax.set_yticklabels(["$%.1f$" % y for y in yticks], fontsize=18)
ax.set_xticklabels(["$%.1f$" % x for x in xticks], fontsize=18);



    


















    






[  0.00000000e+00   2.00000000e+00   4.00000000e+00   6.00000000e+00
   8.00000000e+00   1.00100000e+01   1.20200000e+01   1.40400000e+01
   1.60600000e+01   1.81100000e+01   2.01800000e+01   2.22700000e+01
   2.44000000e+01   2.65800000e+01   2.88200000e+01   3.11200000e+01
   3.35200000e+01   3.60200000e+01   3.86400000e+01   4.14200000e+01
   4.43700000e+01   4.75400000e+01   5.09600000e+01   5.46700000e+01
   5.87400000e+01   6.32200000e+01   6.81900000e+01   7.37400000e+01
   7.99700000e+01   8.70100000e+01   9.50000000e+01   1.04110000e+02
   1.14540000e+02   1.26520000e+02   1.40320000e+02   1.56260000e+02
   1.74680000e+02   1.96000000e+02   2.20660000e+02   2.49170000e+02
   2.82080000e+02   3.20000000e+02   3.63560000e+02   4.13430000e+02
   4.70310000e+02   5.34860000e+02   6.07770000e+02   6.89640000e+02
   7.81040000e+02   8.82430000e+02   9.94160000e+02   1.11645000e+03
   1.24938000e+03   1.39285000e+03   1.54663000e+03   1.71033000e+03
   1.88340000e+03   2.06518000e+03   2.25494000e+03   2.45183000e+03
   2.65500000e+03   2.86358000e+03   3.07672000e+03   3.29359000e+03
   3.51346000e+03   3.73565000e+03   3.95959000e+03   4.18477000e+03
   4.41081000e+03   4.63740000e+03   4.86430000e+03   5.09137000e+03
   5.31851000e+03   5.54567000e+03   5.77283000e+03   6.00000000e+03]










    








<matplotlib.figure.Figure at 0x44b19d0>










    
























In [6]:



    


# KPP related fields 

plt.subplot(2,2,1)
data  = KPP_buoyFlux
field = np.ma.masked_array(data, mask=[data==0.])
plt.contourf(time, -zi, field.T, 16, levels=numpy.arange(-2e-7, 2e-7, 1e-8))
CS = plt.plot(time, -h,'r-',linewidth=1)
plt.colorbar()
plt.xlim((0, 370))
plt.ylim((-30,0))
plt.ylabel(r'z(m) with $\Delta z$ ='+deltaz,fontsize=14)
plt.title(r'Upper ocn buoy flux $(m^2/s^3)$ for '+expt_name,fontsize=18)
#plt.xlabel('Time (days)',fontsize=14)

plt.subplot(2,2,2)
data  = log10(sqrt(KPP_N2))
field = np.ma.masked_array(data, mask=[data==0.])
plt.contourf(time, -zi, field.T, 16, levels=numpy.arange(-5,-1,.25))
#plt.contourf(time, -zi, field.T, 16)
CS = plt.plot(time, -h,'r-',linewidth=1)
plt.colorbar()
plt.xlim((0, 370))
plt.ylim((-300,0))
plt.ylabel(r'z(m) with $\Delta z$ ='+deltaz,fontsize=14)
plt.title(r'$Log_{10}(N)$ $(1/s)$ for '+expt_name,fontsize=18)
plt.xlabel('Time (days)',fontsize=14)


plt.subplot(2,2,3)
data  = KPP_buoyFlux[:,0]
CS = plt.plot(time,data,'*',linewidth=3)
plt.xlim((0, 370))
plt.xlim((0, 370))
plt.grid()
plt.title(r'Surface buoyancy flux for '+expt_name,fontsize=18)
plt.xlabel('Time (days)',fontsize=14)
plt.ylabel(r'Buoyancy flux $(m^2/s^3)$',fontsize=14)


plt.subplot(2,2,4)
CS = plt.plot(time,KPP_uStar,'*',linewidth=3)
#plt.ylim((-400,0))
plt.xlim((0, 370))
plt.grid()
plt.xlabel('Time (days)',fontsize=14)
plt.ylabel(r'$u^{*}$ (m/s))',fontsize=14)
plt.title(r'$u^{*}$ used in KPP for '+expt_name,fontsize=18)
#ax.set_yticklabels(["$%.1f$" % y for y in yticks], fontsize=18)
#ax.set_xticklabels(["$%.1f$" % x for x in xticks], fontsize=18);



    


















    






-c:16: RuntimeWarning: invalid value encountered in sqrt
-c:16: RuntimeWarning: divide by zero encountered in log10










    
























In [7]:



    


# KPP related fields 

plt.subplot(2,2,1)
data   = Kt
field  = np.ma.masked_array(data, mask=[data==0.])
plt.contourf(time, -zi, field.T, levels=numpy.arange(0,.5,.05))
plt.colorbar()
plt.ylim((-300,0))
plt.xlim((0, 370))
plt.ylabel(r'z(m) with $\Delta z$ ='+deltaz,fontsize=14)
plt.title(r'$\kappa_{\Theta} (m^2/s)$ for '+expt_name,fontsize=18)

plt.subplot(2,2,2)
data   = Ku
field  = np.ma.masked_array(data, mask=[data==0.])
plt.contourf(time, -zi, field.T, levels=numpy.arange(0,.5,.05))
plt.colorbar()
plt.xlim((0, 370))
plt.ylim((-300,0))
plt.ylabel(r'z(m) with $\Delta z$ ='+deltaz,fontsize=14)
plt.title(r'$\kappa_u (m^2/s)$ for '+expt_name,fontsize=18)

plt.subplot(2,2,3)
data  = KPP_dTdt
field = np.ma.masked_array(data, mask=[data==0.])
plt.contourf(time, -zl, field.T, 16)
CS = plt.plot(time, -h,'r-',linewidth=1)
plt.colorbar()
plt.xlim((0, 370))
plt.ylim((-300,0))
plt.ylabel(r'z(m) with $\Delta z$ ='+deltaz,fontsize=14)
plt.title(r'Non-local $\Theta$ tendency $(K/s)$ for '+expt_name,fontsize=18)
plt.xlabel('Time (days)',fontsize=14)

plt.subplot(2,2,4)
data  = KPP_dSdt
field = np.ma.masked_array(data, mask=[data==0.])
plt.contourf(time, -zl, field.T, 16)
CS = plt.plot(time, -h,'r-',linewidth=1)
plt.colorbar()
plt.xlim((0, 370))
plt.ylim((-300,0))
plt.ylabel(r'z(m) with $\Delta z$ ='+deltaz,fontsize=14)
plt.title(r'Non-local Saln tendency $(ppt/s)$ for '+expt_name,fontsize=18)
plt.xlabel('Time (days)',fontsize=14)



    


















    
Out[7]:








<matplotlib.text.Text at 0x6bec810>










    
























In [8]:



    


# Boundary layer depth and SST 


figure(4)
fig, ax = plt.subplots(figsize=(16, 10), dpi=dpi)
CS = plt.plot(time, -h,'-',linewidth=3)
plt.grid()
plt.xlabel('Time (days)',fontsize=18)
plt.ylabel(r'z(m) with $\Delta z$ ='+deltaz,fontsize=18)
plt.title(r'KPP boundary layer depth for '+expt_name,fontsize=18)
plt.xlim((0, 400))
xticks = [0, 50, 100, 150, 200, 250, 300, 350, 400]
ax.set_xticklabels(["$%.1f$" % x for x in xticks], fontsize=18);
plt.ylim((-300, 0))
yticks = [-300, -250, -200, -150, -100, -50, 0]
ax.set_yticklabels(["$%.1f$" % y for y in yticks], fontsize=18)

figure(5)
fig, ax = plt.subplots(figsize=(16, 10), dpi=dpi)
CS = plt.plot(time, temp[:,0],'-',linewidth=3)
plt.grid()
plt.xlabel('Time (days)',fontsize=18)
plt.ylabel('SST(C)',fontsize=18)
plt.title(r'SST for '+expt_name,fontsize=18)
plt.xlim((0, 400))
xticks = [0, 50, 100, 150, 200, 250, 300, 350, 400]
ax.set_xticklabels(["$%.1f$" % x for x in xticks], fontsize=18);
plt.ylim((18, 28))
yticks = [18, 20, 22, 24, 26, 28]
ax.set_yticklabels(["$%.1f$" % y for y in yticks], fontsize=18)



    


















    
Out[8]:








[<matplotlib.text.Text at 0x53f12d0>,
 <matplotlib.text.Text at 0x5f2ff90>,
 <matplotlib.text.Text at 0x5cc2750>,
 <matplotlib.text.Text at 0x53f0590>,
 <matplotlib.text.Text at 0x5710e50>,
 <matplotlib.text.Text at 0x5712a90>]










    








<matplotlib.figure.Figure at 0x6688050>










    

















    
























In [9]:



    


# u,v,temp

plt.subplot(2,2,1)
data  = u
field = np.ma.masked_array(data, mask=[data==0.])
plt.contourf(time, -zl, field.T, levels=numpy.arange(-.10,.11,.01))
plt.colorbar()
plt.ylim((-100,0))
plt.ylabel('z (m)',fontsize=14)
plt.title('u (m/s) for expt = '+expt_name,fontsize=16)

plt.subplot(2,2,2)
data  = v
field = np.ma.masked_array(data, mask=[data==0.])
plt.contourf(time, -zl, field.T, levels=numpy.arange(-.10,.11,.01))
plt.colorbar()
plt.ylim((-100,0))
plt.ylabel('z (m)',fontsize=14)
plt.title('v (m/s) for expt = '+expt_name,fontsize=16)

plt.subplot(2,2,3)
data  = temp
field = np.ma.masked_array(data, mask=[data==0.])
CS = plt.contourf(time, -zl, field.T,levels=numpy.arange(18,30,1))
plt.colorbar()
#C = plt.contour(time, -zl, temp.T, 8,linewidth=.5, colors='black')
#plt.clabel(C, inline=1, fontsize=10)
plt.ylim((-100,0))
plt.ylabel('z (m)',fontsize=14)
plt.title(r'$\Theta$ ($\degree$C) for expt = '+expt_name,fontsize=16)
plt.xlabel('Time (days)',fontsize=14);


plt.subplot(2,2,4)
data  = saln[:,:]
field = np.ma.masked_array(data, mask=[data==0.])
CS = plt.contourf(time, -zl, field.T,levels=numpy.arange(36.3,37,.01))
#CS = plt.contourf(time, -zl, field.T)
plt.colorbar()
plt.ylim((-100,0))
plt.ylabel('z (m)',fontsize=14)
plt.title(r'$S$ ($ppt$) for expt = '+expt_name,fontsize=16);
plt.xlabel('Time (days)',fontsize=14);



    


















    
























In [9]:



    






    











In [9]:



    






    











In [1]:



    






    











In [1]:



    






    











In [1]:



    






    











In [1]:



    






    











In [ ]:

Content source: aidanheerdegen/MOM6-examples

Similar notebooks:

notebook.community | gallery | about