Addressing NaN in averaging to fill in blanks without changing previous values

I have 21 matrices (21 different models, 241x121 of gridded temperature data that I have read in from netCDF files. Using GrADS I had gotten an average figure by adding and dividing by 21, but this gives me some blank areas in the tropics, which should not be there.

In Matlab, I have read in the 12 files, concatenated them to create a 241x121x21 3D matrix to then take the nanmean of that file. The regions that were blank in the tropics are now filled in, but regions that had previously had been filled in now have a different value.

Below is my code: I am not sure why values are changing when addressing the NaNs should only fill in the missing values that grads cannot address.

x1=ncread('filepath/name.nc','t1i'); %read in all 21 files

hist_temp=cat(3,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21); %combine 21 matrices
CompH=nanmean(hist_temp,3); %this is the average with the nanmean so blanks from grads figure should be filled.

My overall goal is to address the NaN that may be causing blanks /white space where there should be values in the grads image without changing the values. I have attached a figure of the two results.

Thank you!

[Image Analyst]

Are you sure it's the actual data that is different and not just the colormap?

[Vesp]

I wrote the results as netCDF file and plotted the two figures with the same scaling and got different results. I ended up finding a solution where I used Matlab to conver the NaN to zeros and created a matrix of the sum across the 3D matrix (matrixA). I also created another matrix with the count of non-nan(countx) and wrote those as netCDF files. The final step of dividing matrixA/countx in grads and plotted. The figures using this method have been correct or similar to the method in using grads, but the missing values are now filled in.