KPP/Large Eddy simulation Comparison
The physics: The K-Profile Parameterization (KPP) represents small scale (below the grid scale) vertical turbulent fluxes of heat, salt, and momentum in the ocean.
The code: KPP does not contain any prognostic equations, but utilizes diagnostic equations and scaling relationships to parameterize turbulent fluxes. KPP is known to not exhibit convergence with vertical resolution, but is robust to time step variation. In most cases there are no known analytic solutions for comparison.
The verification: Utilize a large eddy simulation as a baseline to verify the temperature tendency and entrainment depth (figure) returned by KPP, critical quantities for the ocean simulation.
Unit Isolation:CVMix module and MPAS-Ocean vertical mixing interface
In [1]:
from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
import math
%cd /lustre/scratch3/turquoise/lconlon/coastal_test_cases/RUN_ifort.grizzly_hdf5_mpirun_cvmix_compare_ocean/
#import data from Palm
my_nc_file = 'DATA_3D_NETCDF'
fh = Dataset(my_nc_file, mode='r')
u1 = fh.variables['u'][:]
v1 = fh.variables['v'][:]
w1= fh.variables['w'][:]
time1 = fh.variables['time'][:]
t1 = fh.variables['pt'][:]-273.15
s1= fh.variables['sa'][:]
zu1 = fh.variables['zu_3d'][:]
rho1 = fh.variables['rho_ocean'][:]
fh.close()
%cd /lustre/scratch3/turquoise/lconlon/
#import data from cvmix
u2=np.loadtxt('u_velocity.dat')
v2=np.loadtxt('v_velocity.dat')
t2=np.loadtxt('temperature.dat')
d2=np.loadtxt('diff2.dat')
s2=np.loadtxt('salinity.dat')
time2=np.loadtxt('time.dat')
rho2=np.loadtxt('rho.dat')
In [2]:
#Calculate model skill score between the two models
#
# Here, it compares temp, but could switch to another variable
t_palm=np.mean(np.mean(t1,2),2) #horizontal average
palm=t_palm[-1,:] #get palm output after 1 day
cvmix=t2[-1,:] #get cvmix output after 1 day
x1=palm
x2=np.flipud(cvmix)
mse = np.mean((x1 - x2)**2)
mean_standard=np.mean(x1)
denom=np.mean(abs(x2-mean_standard)+abs(x1-mean_standard))**2
fraction=mse/denom
WS=np.abs(1-fraction)
if WS < 0.90:
print("fail")
if WS > 0.90:
print("pass")
print ('model skill for temperature profile is:')
WS
Out[2]:
Here, the model fails because the cvmix run wasn't set up correctly. The notebook then goes on to calculate some stats for more information¶
In [3]:
plt.plot(abs(zu1), x1) #plot palm profile
plt.plot(abs(zu1), x2) #plot CVMIX profile
plt.legend(['Palm','CVMIX'])
plt.title(['Model skill=',WS])
plt.xlabel('depth')
plt.ylabel('temperature')
if WS < 0.95:
pf="fail"
if WS > 0.95:
pf="pass"
plt.title([pf, 'Model skill=', WS])
Out[3]:
In [221]:
#Calculate entrainment depth between the two models (depth of minimum vertical heat flux)
#First, do palm
t_palm=np.mean(np.mean(t1,2),2) #horizontal average
w_palm=np.mean(np.mean(w1,2),2) #horizontal average
palm=(np.diff(np.diff(t_palm[-1,:]*w_palm[-1,:])/np.diff(zu1))) #get palm output after 1 day
m = np.max(palm)
test=[i for i, j in enumerate(palm) if j == m]
palm_depth=np.abs(zu1[test])
# next, do cvmix
t_cvmix=t2[-1,:]
cvmix=(np.diff(np.diff(t2[-1,:]*d2[-1,0:-1])/np.diff(zu1))) #get palm output after 1 day
m = np.max(cvmix)
test=[i for i, j in enumerate(cvmix) if j == m]
cvmix_depth=np.abs(zu1[test])
#calc % difference
percent_diff=(np.abs(palm_depth-cvmix_depth)/((palm_depth+cvmix_depth)/2))*100
print(' Entrainment depth percent difference=')
print(percent_diff)
In [226]:
# calculate gravitational potential energy model skill
# PE=rho_g_h
rho_palm=np.mean(np.mean(rho1,2),2) #horizontal average
palm=np.flipud(rho_palm[-1,:]) #get palm output after 1 day
cvmix=rho2[-1,:] #get cvmix output after 1 day
x1=(palm)
x2=cvmix
mse = np.mean((x1 - x2)**2)
mean_standard=np.mean(x1)
denom=np.mean(abs(x2-mean_standard)+abs(x1-mean_standard))**2
fraction=mse/denom
WS=np.abs(1-fraction)
if WS < 0.90:
print("fail")
if WS > 0.90:
print("pass")
print ('model skill for entrainment depth is:')
WS
Out[226]: