
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('thetas4D','var')==1
        % do nothing
    else
        if exist('thetas','var')==1
            thetas4D = conv2dto4d(thetas,nRows,nCols,nLays,nodes);
        else
            try
                [thetas,thetasTimes] = readOutput(pwd,'th.out',true);
                thetas4D = conv2dto4d(thetas,nRows,nCols,nLays,nodes);
            catch
                disp('Error processing pressure head information-- can''t visualize root extraction.')
            end
        end
    end


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



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

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

    plot(thetasTimes,shiftdim(thetas4D(iRow,iCol,iLay,:)),'-mo',...
        'linewidth',2,'markersize',5,'markerfacecolor','m','markeredgecolor','m')
    title(['water content [-] 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(['water content [-] '])
    xlabel(['time [',swmsTimeUnits,']'])
    set(gca,'xlim',[min(thetasTimes),max(thetasTimes)],...
        'xgrid','on','linewidth',2)

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

    plot(mean([x(1:end-1);x(2:end)],1),dy./dx,'mo',...
        'linewidth',2,'markersize',5,'markerfacecolor','m','markeredgecolor','m')
    title(['d\theta/dt [',swmsTimeUnits,'^{-1}] 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(thetasTimes),max(thetasTimes)],...
        'xgrid','on','ygrid','on','linewidth',2)
    ylabel(['d\theta/dt [',swmsTimeUnits,'^{-1}] '])
    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-theta','-depsc','-r300','-loose')
        print(gcf,'root-extraction-theta','-dpng','-r300')    
    end



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


end % if viz
