Incorrect colorbar range in subplot

I have 4 countour plots that I want to put on one subplot, each with their own colorbar and range. The problem I'm having is that, when I do that, some of the subplots display incorrect colorbar. The colors don't go the full range, so instead of having red for max value it's only green. In the plot it's correct though, and when I do the plot separately (so, not as a part of subplot) it works well.

N=5;
figure

subplot(2,2,1)    %Surface pressure and wind
levels=(960:5:1030);
contourf(longitude,latitude,msl'./100,levels)
contourcmap(960:5:1031,'jet','colorbar','on','TitleString', 'hPa');
hold on
geoshow(coastlat,coastlon,'Color','k', 'LineWidth',1.5)
quiver(longitude(1:N:end),latitude(1:N:end),...
u10(1:N:end,1:N:end,1)',v10(1:N:end,1:N:end,1)','k')
title('MSLP')
set(gca,'XLim',[-20,35],...
        'YLim',[30,75])
    
subplot(2,2,2)     %Temperature and wind at 850hPa
contourf(longitude,latitude,t(:,:,3)'-272.15,-5:5:40)
contourcmap(-5:5:41,'jet','colorbar','on','TitleString', '°C');
hold on
geoshow(coastlat,coastlon,'Color','k', 'LineWidth',1.5)
quiver(longitude(1:N:end),latitude(1:N:end),...
u(1:N:end,1:N:end,3)',v(1:N:end,1:N:end,3)','k')
title('Temperature (850hPa)')
set(gca,'XLim',[-20,35],...
        'YLim',[30,75])
    
subplot(2,2,3)     %Rel. humidity and wind at 700hPa
contourf(longitude,latitude,r(:,:,2)',0:10:100)
contourcmap(0:10:101,'jet','colorbar','on','TitleString', '%');
hold on
geoshow(coastlat,coastlon,'Color','k', 'LineWidth',1.5)
quiver(longitude(1:N:end),latitude(1:N:end),...
u(1:N:end,1:N:end,2)',v(1:N:end,1:N:end,2)','k')
title('Relative humidity (700hPa)')
set(gca,'XLim',[-20,35],...
        'YLim',[30,75])
    
subplot(2,2,4)     %geopotential and wind at 700hPa
contourf(longitude,latitude,z(:,:,1)',52000:500:60000)
contourcmap(52000:500:60001,'jet','colorbar','on','TitleString', 'J/kg');
hold on
geoshow(coastlat,coastlon,'Color','k', 'LineWidth',1.5)
quiver(longitude(1:N:end),latitude(1:N:end),...
u(1:N:end,1:N:end,1)',v(1:N:end,1:N:end,1)','k')
title('Geopotential (500hPa)')
set(gca,'XLim',[-20,35],...
        'YLim',[30,75])

What I get:

What Temperature and humidity scale should look like:

I get that somehow the subplot is affecting the display of the colorbars, but I can't figure out how to make them do what I want.

Edit:I've attached the data I'm using for this.

Answer 1

[Iñigo Escanciano]

It may have to do with CDataMapping being "direct" or the Clim property of the axes being changed for some reason. It would be great having the data to see the values you have are within the limits of these variables and what the plots are really plotting.

[Marin Vojkovic]

Not sure how to upload data on here, but the values are within the ranges. The ranges are chosen to be a bit bigger than the data. The idea is to compare several panels like this, so I want to have a common scale.

Edit: Figured it out. The data is in attachment.