This is challenging because it's a large dataset, and because the shapes of your top and bot arrays are a little awkward. Before I show the interp2 solution, I thought it might be helpful to show a loop version which is slower and not matrix oriented, but maybe helps think about what math needs to be done. This assumes that both top and bot's x,y correspond exactly to the xq and yq values, which they seem to do in this dataset.
thistop = top(top(:,1) == xq(i,j) & top(:,2) == yq(i,j),3);
thisbot = bot(bot(:,1) == xq(i,j) & bot(:,2) == yq(i,j),3);
ind = zq(i,j,:) > thistop | zq(i,j,:) < thisbot;
cq(i,j,squeeze(ind))=nan;
For the interp2 version, we really want a matrix of z values, with x's and y's that look like what meshgrid produces. I'm not sure of a simple way to get that in a generic way based on your matrices, but it's apparent that (for top at least) the x,y values are reshaped versions of such a matrix. (This also appears to be the case for bot, except that the first row is a duplicate).
xi_top=reshape(top(:,1),80,80)';
yi_top=reshape(top(:,2),80,80)';
zi_top=reshape(top(:,3),80,80)';
this makes a matrix for x,y,z where for any row/column the x,y co-ordinates are xi(row,col),yi(row,col) and the z co-ordinate of the top is zi(row,col).
Now we're ready to interpolate:
top_i=interp2(xi_top,yi_top,zi_top,xq,yq,'nearest');
The result top_i is the top that corresponds to each xq,yq. We'll do the same for bot:
xi_bot=reshape(bot(2:end,1),80,80)';
yi_bot=reshape(bot(2:end,2),80,80)';
zi_bot=reshape(bot(2:end,3),80,80)';
bot_i=interp2(xi_bot,yi_bot,zi_bot,xq,yq,'nearest');
Now we have something really useful: bot_i and top_i are co-ordinates of the top and bottom in a shape that exactly matches zq and cq. So marking NaN is as easy as:
cq(zq<bot_i | zq>top_i)=nan;
