function map_plot_sub(lats,lons,vals,titleStr,zString,varargin)

%--------------------------------------------------------------------------
%
% This function plots the specified data on a global map. Input data can
% either be provided individually as vectors/strings or as cell arrays of
% vectors and strings to make multiple plots at a time.
%
%
% INPUT:
%       lats                -        vector of latitudes of the data points
%
%       lons                -        vector of longitudes of the data points
%
%       vals                -        vector of data values
%
%       titleStr            -        string containing title of the plot
%
%       zString             -        string specifying the label for the
%                                    colorbar, i.e. name and unit of z-data
%
%       cLimits             -        optional input argument; two-element 
%                                    vector specifying the limits of the 
%                                    color axis (or cell array of vectors
%                                    in case of multiple plots)
%
%       fName               -        optional input argument; string or
%                                    cell array of string specifying the
%                                    name(s) under which the figure(s) are
%                                    to be saved
%
% HISTORY:
%       09/05/2011      file created                        Jana Kolassa
%--------------------------------------------------------------------------

pixelsize = 4;

%load('NaturalColormap');

riverPlot = 'off';
soilMoisture = 'no';

if strcmp(riverPlot,'on')==1
    
    rdFile = '/Volumes/data1/Data/HydroSheds/rivNet_global';
    rdData = dlmread(rdFile);
    
    thresh = 0.01*(max(rdData)-min(rdData)) + min(rdData);
    
    load('cellNr2Coord_Mat');
    rdLats = cellNrMatrix(:,2);
    rdLons = cellNrMatrix(:,3);
    
    rdIndex = NaN*ones(length(lats),1);
    
    for i = 1:length(lats)
        
        tmpIndex = (rdLats == lats(i)) + (rdLons==lons(i));
        rdIndex(i) = find(tmpIndex==2);
        
    end
    
    rdData = rdData(rdIndex);
    rdData(rdData<thresh) = NaN;
    rdData = isnan(rdData) == 0;
end



% load map data

load('lats');
load('lons');

if iscell(vals) == 1
    
    
    for i = 1:length(vals)
        
        if isempty(varargin) == 0
            
            if length(varargin) == 2
                
                if ischar(varargin{1}{1}) == 0
                    caxis_max = varargin{1}{i}(2);
                    caxis_min = varargin{1}{i}(1);
                    fName = varargin{2}{i};
                else
                    caxis_max = varargin{2}{i}(2);
                    caxis_min = varargin{2}{i}(1);
                    fName = varargin{1}{i};
                end
                
            elseif length(varargin) == 1
                
                if ischar(varargin{1}{1}) == 1
                    fName = varargin{1}{i};
                    caxis_max = max(max(vals{i}));
                    caxis_min = min(min(vals{i}));
                else
                    caxis_max = varargin{1}{i}(2);
                    caxis_min = varargin{1}{i}(1);
                    fName = titleStr{i};
                end
            end
            
            
            
            
        else
            caxis_max = max(max(vals{i}));
            caxis_min = min(min(vals{i}));
            fName = titleStr{i};
        end
        
        if iscell(lats) == 1
            latitudes = lats{i};
            longitudes = lons{i};
        else
            latitudes = lats;
            longitudes = lons;
        end
        
        % reshape matrices into vectors
        
        dim1 = size(vals{i},1);
        dim2 = size(vals{i},2);
        
        if (dim1 > 1) && (dim2 > 1)
            
            vals{i} = reshape(vals{i},dim1*dim2,1);
            latitudes = reshape(latitudes,dim1*dim2,1);
            longitudes = reshape(longitudes,dim1*dim2,1);
            
        end
        
        % plot data
        
       % figure;
        
        Lat_min = min(latitudes);
        Lat_max = max(latitudes);
        Lon_min = min(longitudes);
        Lon_max = max(longitudes);
        
        longitudes(find(longitudes > 180))=longitudes(find(longitudes > 180))-360;
        
        if Lat_min == -99
            Lat_min = min(latitudes);
            Lat_max = max(latitudes);
            Lon_min = min(longitudes);
            Lon_max = max(longitudes);
            Lat_min = Lat_min-1;
            Lat_max = Lat_max+1;
            Lon_min = Lon_min-1;
            Lon_max = Lon_max+1;
        end
        
        val_min = min(vals{i});
        val_max = max(vals{i});
        
        % trace coastlines
       % fig1=figure
        plot(long(:,1),lat(:,1), 'k-', 'MarkerSize', 5);
        hold on;
        
        colormap('jet')
        cmap = colormap;
        if strcmp(soilMoisture,'yes')==1
            cmap = mycmap;
            set(gcf,'Colormap',mycmap);
        end
        clen = length(cmap);
        range_c = caxis_max - caxis_min;
        thecolor = min(max(round((vals{i}(:)-caxis_min)*(clen-1)/range_c),0),clen);
        
        
        
        
        for il=clen:-1:1
            list_color = find(thecolor  == il);
            plot(longitudes(list_color(:)),latitudes(list_color(:)),'.','MarkerSize',  pixelsize, 'color', cmap(il,:));
        end
        
        a=gca();
        set(a,'DataAspectRatioMode','auto');
        
        caxis([caxis_min caxis_max]);
        
        h = colorbar;
        cbfreeze(h);
        cblabel(zString{i},'FontSize',12,'Rotation',270);
        %set(get(h,'ylabel'),'String', zString{i}, 'Rotation', 270,'VerticalAlignment','Bottom','fontsize',12);
        xlabel('Longitude [deg]','fontsize',12);
        ylabel('Latitude [deg]','fontsize',12);
        axis equal
        axis([Lon_min Lon_max Lat_min Lat_max])
        set(gca,'FontSize',12);
        set(gca,'LineWidth',1.5);
        freezeColors(gca)
        title(titleStr{i},'fontsize',12);
        hold on
        
        if strcmp(riverPlot,'on')==1
            
            plot(longitudes(rdData),latitudes(rdData),'.k','MarkerSize',3)
        end
        hold off;

        
        % save figures to figure folder
        
        
%         figurePath = [pwd '/figures/'];
%         figureName = [figurePath fName];
%         
%         
%         saveas(gcf, [figureName '.jpg'],'jpg');
        
     end
    
else
    
        if isempty(varargin) == 0
            
            if length(varargin) == 2
                
                if ischar(varargin{1}) == 0
                    caxis_max = varargin{1}(2);
                    caxis_min = varargin{1}(1);
                    fName = varargin{2};
                else
                    caxis_max = varargin{2}(2);
                    caxis_min = varargin{2}(1);
                    fName = varargin{1};
                end
            elseif length(varargin) == 1
                if ischar(varargin{1}) == 1
                    fName = varargin{1};
                    caxis_max = max(max(vals));
                    caxis_min = min(min(vals));
                else
                    caxis_max = varargin{1}(2);
                    caxis_min = varargin{1}(1);
                    fName = titleStr;
                end
            end
            
            
            
            
        else
            caxis_max = max(max(vals));
            caxis_min = min(min(vals));
            fName = titleStr;
        end
    

    
    
    
    % reshape matrices into vectors
    
    dim1 = size(vals,1);
    dim2 = size(vals,2);
    
    if (dim1 > 1) && (dim2 > 1)
        
        vals = reshape(vals,dim1*dim2,1);
        lats = reshape(lats,dim1*dim2,1);
        lons = reshape(lons,dim1*dim2,1);
        
    end
    
    % plot data
    
    
    
    Lat_min = min(lats);
    Lat_max = max(lats);
    Lon_min = min(lons);
    Lon_max = max(lons);
    
    lons(find(lons > 180))=lons(find(lons > 180))-360;
    
    if Lat_min == -99
        Lat_min = min(lats);
        Lat_max = max(lats);
        Lon_min = min(lons);
        Lon_max = max(lons);
        Lat_min = Lat_min-1;
        Lat_max = Lat_max+1;
        Lon_min = Lon_min-1;
        Lon_max = Lon_max+1;
    end
    
    val_min = min(vals);
    val_max = max(vals);
    
    % trace coastlines
    %fig1=figure
    plot(long(:,1),lat(:,1), 'k-', 'MarkerSize', 5);
    hold on;
    
    colormap('jet')
    %colormap(flipud(colormap))
    cmap = colormap;
    if strcmp(soilMoisture,'yes')==1
        cmap = mycmap;
        set(gcf,'Colormap',mycmap)
    end
    clen = length(cmap);
    range_c = caxis_max - caxis_min;
    thecolor = min(max(round((vals(:)-caxis_min)*(clen-1)/range_c),0),clen);
    
    
    
    
    for il=clen:-1:1
        list_color = find(thecolor  == il);
        plot(lons(list_color(:)),lats(list_color(:)),'.','MarkerSize',  pixelsize, 'color', cmap(il,:));
    end
    
    a=gca();
    set(a,'DataAspectRatioMode','auto');
    
    caxis([caxis_min caxis_max]);
    
    h = colorbar;
    cbfreeze(h);
    cblabel(zString,'FontSize',12,'Rotation',270);
    %set(get(h,'ylabel'),'String', zString, 'Rotation', 270,'VerticalAlignment','Bottom','fontsize',12);
    xlabel('Longitude [deg]','fontsize',12);
    ylabel('Latitude [deg]','fontsize',12);
    axis equal
    axis([Lon_min Lon_max Lat_min Lat_max])
    set(gca,'FontSize',12);
    set(gca,'LineWidth',1.5);
    freezeColors(gca)
    title(titleStr,'fontsize',12);
    
    hold on
    
    if strcmp(riverPlot,'on')==1
        plot(lons(rdData),lats(rdData),'.k','MarkerSize',3)
    end
    
    hold off;
    
    
    % save figures to figure folder
    
    
%     figurePath = [pwd '/figures/'];
%     figureName = [figurePath fName];
%     
%     
%     saveas(gcf, [figureName '.jpg'],'jpg');
    
end