Interpolate 3D surfaces along Z axis to extrapolate new surfaces

I am trying to interpolate three 3D surfaces to extrapolate other 3D surfaces along the Z-axis.

In the code below I show my three surfaces. What I need is to 1) create N surfaces in between them, and 2) possibly extrapolate other surfaces below and above the first and the last surfaces.

I tried with RegularGridInterpolator, but I am able only to derive some points. Is there a way to obtain entire surfaces (set of points), instead that only one point?

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate

mat1 = np.array([     [60, 50, 40, 30],
                      [55, 45, 35, 25],
                      [44, 34, 24, 14],
                      [38, 28, 18,  8] ])

mat2 = mat1 + 20

mat3 = mat1 + 40

x = np.linspace(0,1,4)
y = np.linspace(0,100,4)
tot = np.array([mat1, mat2, mat3])

x2 = np.linspace(0,1,3)
y2 = np.linspace(0,100,4)
z2 = np.linspace(0,100,4)

fig, ax = plt.subplots(figsize=(6,6),subplot_kw=dict(projection='3d'))

# plot orgiginal surfaces
for zg in tot:
    xg, yg = np.meshgrid(x, y)
    s = ax.plot_surface(xg, yg, zg, alpha=0.8)

# grid interpolator
xg2, yg2, zg2 = np.meshgrid(x2, y2, z2)
interp = interpolate.RegularGridInterpolator((x2, y2, z2), tot)

print(interp([0.5,50,50])) # this gives me only one point z-point in x=0.5 and y=50