fix(plotutil): number of plottable layers should be from util3d.shape[0]

[mwtoews]

Previous behavior assumed that .plot() on a util3d object had as many layers as model.modelgrid.nlay. However, this isn't always correct, e.g. with BCF's vcont dataset, which is nlay - 1:

flopy/flopy/modflow/mfbcf.py

Lines 209 to 217 in aca1446

	if model.nlay > 1:
	self.vcont = Util3d(
	model,
	(nlay - 1, nrow, ncol),
	np.float32,
	vcont,
	"Vertical Conductance",
	locat=self.unit_number[0],
	)

this would cause ml.bcf6.vcont.plot() (or ml.plot()) to raise "IndexError: list index out of range"

---

Also improved test_model_dot_plot to complete in ~2ish seconds instead of ~22 seconds (compare 20 plots in a smaller model to previous 226 plots from a more complex model). The return object from .plot() is also checked, as this is either a list or plt.Axes object, depending on the number of plots.

[codecov]

Codecov Report

Merging #1077 (4917b3b) into develop (aca1446) will decrease coverage by 0.413%.
The diff coverage is 100.000%.

@@              Coverage Diff              @@
##           develop     #1077       +/-   ##
=============================================
- Coverage   71.764%   71.351%   -0.414%     
=============================================
  Files          225       225               
  Lines        51893     51887        -6     
=============================================
- Hits         37241     37022      -219     
- Misses       14652     14865      +213     

Impacted Files	Coverage Δ
flopy/plot/plotutil.py	46.175% <100.000%> (-1.767%)	⬇️
flopy/mt3d/mtsft.py	6.569% <0.000%> (-72.263%)	⬇️
flopy/utils/formattedfile.py	88.311% <0.000%> (-2.598%)	⬇️
flopy/utils/mfreadnam.py	75.912% <0.000%> (-1.460%)	⬇️
flopy/modflow/mflmt.py	91.578% <0.000%> (-1.053%)	⬇️
flopy/mbase.py	70.858% <0.000%> (+0.120%)	⬆️
flopy/utils/util_array.py	83.724% <0.000%> (+0.227%)	⬆️
flopy/mf6/modflow/mfgwtfmi.py	61.904% <0.000%> (+15.476%)	⬆️