Comparison for new and old MPAS test case to make sure the model gives the same results in both cases. Nodebook calculates a model skill score for temperature at 10 random locations.¶
In [1]:
from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
import math
import random
from matplotlib import transforms
%cd /net/scratch3/lconlon/MPAS-Model
#import data from control run
my_nc_file='temp_control.nc'
fh = Dataset(my_nc_file, mode='r')
T_control = fh.variables['TEMP'][:]
depth=fh.variables['depth_t'][:]
lon=fh.variables['t_lon'][:]
lat=fh.variables['t_lat'][:]
fh.close()
#import new data
my_nc_file_new='PTemp.Jan_p3.filled.mpas100levs.160127.nc'
fh = Dataset(my_nc_file_new, mode='r')
T_new = fh.variables['TEMP'][:]
fh.close()
# Choose 10 random points for comparison between the 2 runs and get depth, lat, lon for those points
rand_val_lon=np.zeros(10,dtype=np.int8)
rand_val_lat=np.zeros(10,dtype=np.int8)
for x in range(10):
rand_val_lon[x] = random.randint(1,361)
rand_val_lat[x] = random.randint(1,180)
lon_sub=lon[rand_val_lon]
lat_sub=lat[rand_val_lat]
T_new_sub=T_new[:,rand_val_lat,rand_val_lon]
T_control_sub=T_control[:,rand_val_lat,rand_val_lon]
In [2]:
#Calculate model skill score between the new and old run for the 10 random locations
# Here, it compares temp, but could switch to another variable
WS=np.zeros(10,dtype=np.int8)
for x in range(10):
x1=T_control_sub[:,x]
x2=T_new_sub[:,x]
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[x]=np.abs(1-fraction)
# if model skill is less that 0.9, the test fails.
if WS[x] < 0.90:
print("fail")
if WS[x] > 0.90:
print("pass")
base = plt.gca().transData
rot = transforms.Affine2D().rotate_deg(270)
# plot temperature for new and old runs at the random points to visualize differences
plt.plot(abs(depth), x1,'r',transform= rot + base)
plt.plot(abs(depth), x2,'b',transform= rot + base)
plt.legend(['control','new'])
plt.title(['Model skill=',WS])
plt.ylabel('depth')
plt.xlabel('temperature')
# print skill scores for all random locations
print ('model skill for temperature profile at 10 locations is:')
WS
Out[2]:
In [ ]: