
if viz
    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 groundwater.')
            end
        end
    end


    nTimes = size(presHeads4D,4);

    hFig = subplotScreen;
    set(hFig,'numbertitle','off','name','groundwater',...
        'inverthardcopy','off','paperpositionmode','auto')

    TMP = presHeads4D(:,:,1,:);
    minGW = 0;
    maxGW = max([0,max(TMP(:))]);
    if maxGW==minGW
        minGW = -1;
        maxGW = 1;
    end

    minh = min(presHeads4D(:));
    maxh = max(presHeads4D(:));
    if minh==maxh
        minh = minh-1;
        maxh = minh+1;
    end

    clear TMP

    for iTime = 1:nTimes

        subplot(1,2,1,'fontsize',16)
        imagesc(nodes.X(1,:,1),nodes.Y(:,1,1),presHeads4D(:,:,1,iTime),...
            [minh,maxh])

        xlabel(['x [',swmsLengthUnits,']'])
        ylabel(['y [',swmsLengthUnits,']'])
        title(['pressure head [',swmsLengthUnits,'] @ time = ',num2str(presHeadTimes(iTime)),' [',swmsTimeUnits,']'])
        set(gca,'ydir','normal','xtick',nodes.X(1,:,1),...
                                'ytick',nodes.Y(:,1,1))
        axis image
        colorbar


        subplot(1,2,2,'fontsize',16)
        imagesc(nodes.X(1,:,1),nodes.Y(:,1,1),...
            (presHeads4D(:,:,1,iTime)>=0).*presHeads4D(:,:,1,iTime),...
            [minGW,maxGW])

        xlabel(['x [',swmsLengthUnits,']'])
        ylabel(['y [',swmsLengthUnits,']'])
        title(['transient saturation [',swmsLengthUnits,'] @ time = ',num2str(presHeadTimes(iTime)),' [',swmsTimeUnits,']'])
        set(gca,'ydir','normal','xtick',nodes.X(1,:,1),...
                                'ytick',nodes.Y(:,1,1))
        axis image
        colorbar

        drawnow



        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(hFig,['transient-saturation_',sprintf('%03d',iTime)],'-depsc','-r300','-loose')
            print(hFig,['transient-saturation_',sprintf('%03d',iTime)],'-dpng','-r300')    
        else
            pause(0.05)
        end

    end



    clear nTimes
    clear iTime
    clear hFig
    clear maxGW
    clear minGW
    clear r
    clear maxWidth
    clear pos
    clear minh
    clear maxh
    
end % if viz    
    