
if viz

    if ~exist('nodes','var')==1
        warning(['Variable ',char(39),'nodes',char(39),' does not exist. ',...
            'Skipping rest of visualization routine.'])
        return
    end


    if exist('presHeads4D','var')==1
        % do nothing
    else
        if exist('presHeads','var')==1
            presHeads4D = conv2dto4d(presHeads,nRows,nCols,nLays,nodes);
        else
            try
                [presHeads,presHeadTimes] = readOutput(pwd,'h.out',true);
                presHeads4D = conv2dto4d(presHeads,nRows,nCols,nLays,nodes);
            catch
                disp('Error processing pressure head information-- can''t visualize root extraction.')
            end
        end
    end


    if ~exist('presHeadTimes','var')==1


        presHeads4D = repmat(NaN,[size(nodes.ID),numel(presHeadTimes)]);
        for iPrintTime=1:numel(presHeadTimes)
            tmp = repmat(NaN,size(nodes.ID));
            tmp(1:nodes.n) = presHeads(nodes.ID(:),iPrintTime);
            presHeads4D(:,:,:,iPrintTime) = tmp;
        end
        clear iPrintTime
        clear nPrintTimes
        clear tmp
    end

    if ~(exist('presHeads4D','var')==1)
        disp('variable presHeads4D does not exist -- can''t plot root extraction.')
        return
    end



    hFig = subplotScreen;
    subplot(2,1,1,'fontsize',16)
    set(hFig,'numbertitle','off','name','h - root extraction',...
        'paperpositionmode','auto','inverthardcopy','off')

    [iRow,iCol,iLay] = find(nodes.X==monitor.x&...
                            nodes.Y==monitor.y&...
                            nodes.Layer==monitor.layer);

    plot(presHeadTimes,shiftdim(presHeads4D(iRow,iCol,iLay,:)),'-bo',...
        'linewidth',2,'markersize',5,'markerfacecolor','b','markeredgecolor','b')
    title(['pressure head [', swmsLengthUnits,'] for node at (X,Y,Z)=(',...
        num2str(nodes.X(iRow,iCol,iLay)),',',...
        num2str(nodes.Y(iRow,iCol,iLay)),',',...
        num2str(nodes.Z(iRow,iCol,iLay)),') [',swmsLengthUnits,']'])
    ylabel(['pressure head [', swmsLengthUnits,'] '])
    xlabel(['time [',swmsTimeUnits,']'])
    set(gca,'xlim',[min(presHeadTimes),max(presHeadTimes)],...
        'xgrid','on','linewidth',2)

    subplot(2,1,2,'fontsize',16)
    x = presHeadTimes;
    dx = x(2:end)-x(1:end-1);
    y = shiftdim(presHeads4D(iRow,iCol,iLay,:))';
    dy = y(2:end)-y(1:end-1);

    plot(mean([x(1:end-1);x(2:end)],1),dy./dx,'bo',...
        'linewidth',2,'markersize',5,'markerfacecolor','b','markeredgecolor','b')
    title(['dh/dt [', swmsLengthUnits,'/',swmsTimeUnits,'] for node at (X,Y,Z)=(',...
        num2str(nodes.X(iRow,iCol,iLay)),',',...
        num2str(nodes.Y(iRow,iCol,iLay)),',',...
        num2str(nodes.Z(iRow,iCol,iLay)),') [',swmsLengthUnits,']'])
    set(gca,'xlim',[min(presHeadTimes),max(presHeadTimes)],...
        'xgrid','on','ygrid','on','linewidth',2)
    ylabel(['dh/dt [', swmsLengthUnits,'/',swmsTimeUnits,'] '])
    xlabel(['time [',swmsTimeUnits,']'])


    clear x dx y dy



    set(hFig,'paperunits','centimeter')
    pos = get(hFig,'position');
    r = pos(3)/pos(4);
    maxWidth = 26;
    set(hFig,'paperposition',[1,1,maxWidth,maxWidth/r])



    if saveFigs
        print(gcf,'root-extraction-h','-depsc','-r300','-loose')
        print(gcf,'root-extraction-h','-dpng','-r300')    
    end



    clear iRow
    clear iCol
    clear iLay
    clear r
    clear maxWidth
    clear pos
    clear hFig

end % if viz
    