function coasts(detail,colr,nancolr,edgecolr,pacflag);

% coasts(detail,colr,nancolr,edgecolr,pacflag);

% coasts - overlay the high-resolution HSHHS coasts in
%   /phod008/lumpkin/bathystuff/gshhs.
% colr=color of coastline.  
% detail=0 to 10; 0-9 plots progressively more of the 
%   polygons at increasing resolution, 10 plots everything
%   visible within the current axes.
% pacflag=1 to repeat coasts for longitudes through +360+30

load hrcoasts;

if ~exist('colr')
    colr=0;
end;
if (colr==1 & length(colr)==1)
    col=[140 120 80]/170;  % beigeish-brown
	%col=[.5 .7 .5];		% green-grey
elseif (colr==0) col=[.6 .6 .6];		% grey
elseif (~length(colr)) col=[1 1 1];
else col=colr;
end;

if ~exist('detail') detail=1; end;
% calculate the number of polygons(continents, islands, lakes)
%   to include in the plot
if (detail==0)
    numpolys=30;
elseif (detail>0 & detail<=9)   % steady increase to detail 9
    numpolys=500*detail;
elseif (detail==10)
    numpolys=length(qp);
end;

if ~exist('nancolr')
    nancolr=[1 1 1]; %white
end;
if ~length(nancolr)
    nancolr=[1 1 1]; %white
end;

if ~exist('edgecolr')
    if (detail<2 & sum(col)<3) 
        edgecolr='none';
    else edgecolr='k';
    end;
end;

if ~exist('pacflag');
    lim=get(gca,'xlim');
    if (max(lim)>200) pacflag=1;
    else pacflag=0;
    end;
end;

isheld=ishold;
xlim=get(gca,'xlim');
ylim=get(gca,'ylim');

skipind=max(1,floor(1e4/numpolys));

hold on;
for i=2:numpolys;
    a=qp(i-1)+1:qp(i)-1;
    x=coast_x(a);
    y=coast_y(a);
    
    % subsample the polygons according to detail level
    if (length(x)>20)
        x=x(1:skipind:length(x));
        y=y(1:skipind:length(y));
    end;
    
    if (min(x)>180 & ~pacflag)
        x=x-360;
    elseif (max(x)<0 & pacflag)
        x=x+360;
    end;
    
    % don't plot if detail level is >7 and polygon
    %   falls outside current limits of plot
    
    if (detail>7) 
        if (max(x)<min(xlim) | min(x)>max(xlim) | ...
            max(y)<min(ylim) | min(y)>max(ylim))
                x=[]; y=[];
        else
            q=find(x>min(xlim)-.1 & x<max(xlim)+.1 & ...
                y>min(ylim)-.1 & y<max(ylim)+.1);
            x=x(q); y=y(q);
        end;
    end;
    if (length(x))
        
            
        antflag=0;
    
        if (i==2 & ~pacflag)
            % repeat a bit of Siberia west of Alaska
            q=find(x>175);
            tx=x(q)-360;
            ty=y(q);
            if (sum(col)<3)    % all-white: plot edges only
                patch(tx,ty,col,'edgecolor',edgecolr);
            else
                q=find(abs(diff(ty))>5);
                tx(q)=NaN;ty(q)=NaN;
                plot(tx,ty,'color',edgecolr);
            end;
            
        elseif (i==8); % Antarctica
            % flip, so coastline runs west-to-east
            x=flipud(x); y=flipud(y);
            if ~pacflag
                q1=find(x>180); 
                q2=find(x<=180);
                tx=[x(q2)];
                ty=[y(q2)];
                if (sum(col)<3)
                    tx=[tx;181;181;-181;-181];
                    ty=[ty;last(y(q2));-91;-91;y(q1(1))];
                else
                    tx=[tx;180;NaN]; ty=[ty;last(ty);NaN];
                end;
                tx=[tx;x(q1)-360];
                ty=[ty;y(q1)];
                x=tx; y=ty;
            else
                tx=[x;x+360];
                ty=[y;y];
                if (sum(col)<3)
                    tx=[tx;721;721;-1;-1];
                    ty=[ty;last(ty);-91;-91;y(1)];
                else
                    tx=[tx;720;NaN]; ty=[ty;last(ty);NaN];
                end;
                x=tx; y=ty;
            end;
            antflag=1;
        end;
        polycol=col;
        if (inland(i))
            polycol=nancolr;
        else
             polycol=col;
        end;
        if (sum(col)<3) % patch normally (with edgecolor)
                        % if white hasn't been specified as
                        % the fill color
            patch(x,y,polycol,'edgecolor',edgecolr);
            if (min(x)<360 & pacflag & ~antflag)
                patch(x+360,y,polycol,'edgecolor',edgecolr);
            end;
        else    % all-white: plot edges only
            plot(x,y,'color',edgecolr);
        end;
    
    end;    % loop setting x=[] if (a) detail>5, 
            % (b) coast poly outside current xlim,ylim
end;


if ~isheld
    hold off;
end;
