Interpolating grid with NaN values for plotting with contour

I have a table I create from a larger dataset which gives me the y, x, and z values that I want to plot. I create a grid and then interpolate the values, but I don't think the function interpolation along the rows correctly since I can see data as break-up points around Latitude 78-80 in the plot (see image), this makes me think that the interpolation is not been done correctly along the rows. Does anyone have a tip on how to smooth this data?

aou_df = df_1994.pivot_table(index='CTDPRS', columns = 'LATITUDE', values='AOU')
    aou_df = aou_df.interpolate(method='linear', limit_area='inside', axis =0 )
##Plotting AOU 1994
y =  ([   4.2,    4.7,    4.8,    4.9,    5.4,    9.1,    9.6,    9.7,
                    10.0,   10.1,
                  ...
                  3568.2, 3608.6, 3818.6, 3824.9, 3866.7, 3979.1, 4013.4, 4133.1,
                  4159.3, 4287.3],

x= ([72.13,  73.0, 73.49, 73.98,  74.5,  75.0, 75.45, 75.75, 75.94,
              76.62, 77.33, 77.78, 78.14, 78.15, 78.98, 79.98, 80.15, 80.16,
              80.33, 80.71, 81.24, 81.58, 82.47, 83.17, 84.06, 84.85, 85.89,
              87.16, 88.06, 88.79, 88.86, 88.95, 89.02,  90.0],

z = [[-12.29372749,          nan,          nan, ...,          nan,
             nan,          nan],
   [         nan,          nan, -43.41465869, ...,          nan,
             nan,          nan],
   [         nan, -54.49999783,          nan, ...,          nan,
             nan,          nan],
   ...,
   [         nan,          nan,          nan, ...,          nan,
             nan,          nan],
   [         nan,          nan,          nan, ...,          nan,
     55.87256821,          nan],
   [         nan,          nan,          nan, ...,  55.39665852,
             nan,  55.05005376]])


xi, yi = np.meshgrid(x,y,indexing='ij')

#from matplotlib.colors import LogNorm
plt.figure(figsize=(25,10))
levels = np.linspace(-135,135)
#cbar = plt.colorbar(ticks=(-85,-65,-45,-25,-5,15,35,55,75,95, 115,135))

plt.contourf(xi,yi,z, cmap = 'jet', levels=levels,vmin=-135, vmax=135)
plt.gca().invert_yaxis()
plt.gca().invert_xaxis()
cbar = plt.colorbar(ticks=(-135,-110,-85,-60,-35,0,35,60,85,110,135), extend= 'both')
cbar.set_label('AOU', fontsize=18)
cbar.ax.tick_params(labelsize=18)
plt.xlabel('LAT',fontsize=18)
plt.ylabel('Pressure (dbar)' ,fontsize=18)
plt.ylim(bottom = 1000)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.plot(x,range(len(x)),'gD', clip_on=False, markersize=10)
#plt.xlim(left = 80)

Ok, the values are in a table... I'm sorry... I updated the code... Try to run...

Try this... You can to pass method='Linear' and method='Cubic' to griddata

from scipy.interpolate import griddata

### getting the valid data
xdata, ydata = np.meshgrid(x, y)
not_nan = ~np.isnan(xdata) & ~np.isnan(ydata) & ~np.isnan(z)
i, j = np.argwhere(not_nan).T
xval = xdata[i, j]
yval = ydata[i, j]
zval = z[i, j]
xy = np.column_stack((xval, yval))

### mesh
nx_mesh = 100
ny_mesh = 100
xi, yi = np.meshgrid(np.linspace(xval.min(), xval.max(), nx_mesh), 
                     np.linspace(yval.min(), yval.max(), ny_mesh))

### interpolation
zi = griddata(xy, zval, (xi, yi), method='nearest')

### Pass xi, yi and zi to contour....

EDIT: This is the result from Joao's code using method='nearest'

This is the 'linear' result: