% This code reads surface temperature 
% data from the MPI_ECHAM5 model scenario A2,


file='tas_A1.nc';
lonvarid = 'lon';
latvarid = 'lat';
Tvarid = 'tas';

ncid = netcdf.open(file,'NC_NOWRITE')
numlonvarid=netcdf.inqVarID(ncid,lonvarid)
numlatvarid=netcdf.inqVarID(ncid,latvarid)
numTvarid=netcdf.inqVarID(ncid,Tvarid)


% Getting longitude data
lonnc = netcdf.getVar(ncid,numlonvarid);

% Getting longitude data
latnc = netcdf.getVar(ncid,numlatvarid);

% Getting the temperature data
T = netcdf.getVar(ncid,numTvarid);

% You could also get a subset of the data between lat1, lat2, lon1, lon2,
% time1, time2...we won't do this.
% pr = getnc(file,Tvarid,[time1,lat1,lon1],[time2,lat2,lon2],[1,1,1],-2);

% Now let's map the data for the first time step

%Get the time slice
T2 = squeeze(T(:,:,1));
T1 = double(T2)

%Let's assume that missing data = 9999;
T1(T1==9999) = NaN;


% We need to make a matrix of latitudes and longitudes
%The longitudes must be changed to -180 180
lon2 = lonnc;
lon2(lon2>180) = lon2(lon2>180)-360;

%Now let's make the matrices to match the size of T1
lonP = repmat(lon2,1,size(latnc,1));
latP = repmat(latnc,1,size(lonnc,1))';

%Define the characteristics of the map
% use projlist('mapprojection') to get a list of projections
axesm('MapProjection','mercator','Origin',[0 180 0],'FLatLimit',[20 60],...
    'FLonLimit',[-180 180], 'MLineLocation',20, ...
    'PLineLocation',20, 'MeridianLabel','on', 'ParallelLabel','on',...
    'LabelFormat','signed')
gridm on;
framem on;
% Plot coast lines, the coast comes with Matlab as well as lat and long
load coast
h = patchm(lat,long,'b','FaceColor','none');
%Plot the actual data
surfm(latP,lonP,T1)
%We can make the color scale different
caxis([200 310])
%Change the colormap colors
set(gcf,'Colormap',hsv)