KPP/Large Eddy simulation Comparison

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

/lustre/scratch3/turquoise/lconlon/coastal_test_cases/RUN_ifort.grizzly_hdf5_mpirun_cvmix_compare_ocean
/lustre/scratch3/turquoise/lconlon

#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

fail
model skill for temperature profile is:

0.14229683510646618

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

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

Text(0.5,1,u"['fail', 'Model skill=', 0.14229683510646618]")

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

 Entrainment depth percent difference=
[23.58078603]

# 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

fail
model skill for entrainment depth is:

0.4118265192298509