%
%   This is a program to read NCEP-DOE II reanalysis monthly mean 
% data using matlab's own functions, for any question plz contact
% yangzhang@nju.edu.cn
%
%----------------------------------------------------------------
%
%--Below are parameters you have to input or vary to finish Q1
%--Input #1: the full path of the data file, 
%            for temperature, use data file: air.mon.mean.nc
%            for zonal wind, use data file: uwnd.mon.mean.nc
path_data='C:\YOURPATH\uwnd.mon.mean.nc';
%
%--Input #2: the year you want to average 
%            This data convers from 1979/1 -- 2009/12
%            e.g. averaged from 1981 to 2000, input: year=[1981 2000]
%            DO NOT input any value smaller than 1980 (as winter including 
%            Dec. of the previous year) and make sure that the first input 
%            value is NOT LARGER than the latter one. 
year=[1980 1981];
%
%--Input #3: the season you want to average
%            Winter (DJF), input: season=1
%            Spring (MAM), input: season=2
%            Summer (JJA), input: season=3
%            Fall   (SON), input: season=4
%            All year,     input: season=0
%
season=0;
%
%--Input #4: the vertical level you want to plot, check air_readme.txt and
%            uwnd_readme.txt for the corresponding vertical level
%            e.g. if plot at 850hpa, input: level=3
level=6;
%
%--Input #5: full path to access the landmask data
%
path_land='C:\YOURPATH\land.nc';
%
%
%--YOU ARE DONE WITH THE INPUT!!! Just wait for the pics!
%
%--------------------------------------------------------------
%--STEP1: read nc files
%
ncid=netcdf.open(path_data,'NOWRITE');
[ndims,nvars,ngatts,unlimdimid] = netcdf.inq(ncid);
% number of dimensions, number of variables, number of global attributes, 
% and the identity of the unlimited dimension, if any.
%
%---read info on lon,lat,lev
dimname=cell(1,5);
for i=1:ndims
 [dimname{i}, dimlen(i)] = netcdf.inqDim(ncid,i-1);
 %returns the name and length of a dimension
 %length of lon,lat,lev,time 
end
lev=netcdf.getVar(ncid,0);
lat=netcdf.getVar(ncid,1);
lon=netcdf.getVar(ncid,2);
nbd=netcdf.getVar(ncid,3); % a dimension only appears in reanalysis II 
%
%----Begins to read the variables
%--Read parameters
[varname,xtype,dimids,natts] = netcdf.inqVar(ncid,5);
%returns the name, datatype, dimensions IDs, 
%and the number of attributes of the variable identified by varid.
%--ABOVE COMMAND IS OPTIONAL
%
scale=netcdf.getatt(ncid,5,'scale_factor');
offset=netcdf.getatt(ncid,5,'add_offset');
%--NOTE: var_real=var_read*scale+offset
%
if (season==0) 
    m_int=(year(1)-1979)*12; d_mon=(year(2)-year(1)+1)*12;
    tmp=double(...
        netcdf.getVar(...
        ncid,5,[0 0 level-1 m_int],[dimlen(1) dimlen(2) 1 d_mon])...
        ).*scale+offset;
else
    d_mon=(year(2)-year(1)+1)*3;
    i_mon=0;
    tmp=zeros(dimlen(1),dimlen(2),1,d_mon);
    for y_ind=year(1):year(2)
        tmp(:,:,:,i_mon+(1:3))=double(...
            netcdf.getVar(...
            ncid,5,[0 0 level-1 (y_ind-1979)*12+(season-1)*3-1],...
            [dimlen(1) dimlen(2) 1 3])...
            ).*scale+offset; 
        i_mon=i_mon+3;
    end   
end
var=squeeze(mean(tmp,4));
clear tmp;
%
%-----DONE WITH THE DATA READING, START THE PLOT
%
t_name=cell(5,1);
t_name{1}='Annual mean';  
t_name{2}='DJF mean'; t_name{3}='MAM mean'; 
t_name{4}='JJA mean'; t_name{5}='SON mean'; 
%
figure(2)
subplot(3,3,1:6)
[C,h]=contourf(lon,lat,var','LineStyle','none');hold on;
p=get(h,'Children'); 
c=get(p,'Cdata'); 
set(p([c{:}]<0),'LineStyle','--');
set(p([c{:}]==0),'LineStyle','-','Linewidth',3);
clabel(C,'fontsize',12,'color','k');
ylim([-80,80]);
set(gca,'Fontweight','bold','Fontsize',14); 
set(gca,'ytick',[-80 -60 -45 -30 -15 0 15 30 45 60 80]);
set(gca,'yticklabel',...
         {'SP','60S','','30S','','EQ','','30N','','60N','NP'});
title([t_name{season+1},' at ',num2str(lev(level)),' hpa from ',num2str(year(1)),' to ',num2str(year(2))]);
%
%----add the land sea distribution
ncid_land=netcdf.open(path_land,'NOWRITE');
land=netcdf.getVar(ncid_land,3);
v=[-1 1];
contour(lon,lat,land',v,'-k','Linewidth',2);
%
%---specify the range of colormap
v_max=max(max(var)); v_min=min(min(var));
caxis([v_min v_max]);
%