function masked_subarray, array, mask, index

Sz = size(array)

if Sz[0] eq 2 then N = 1 else N = Sz[3]

scany = where(total(mask, 1) ne 0)
scanx = where(total(mask, 2) ne 0)

Xmin = min(scanx, max = Xmax)
Ymin = min(scany, max = Ymax)

index =  [Xmin, Ymin]

frag = array[Xmin:Xmax, Ymin:Ymax, *]

fragmask = mask[Xmin:Xmax, Ymin:Ymax]

for j=0, N-1 do frag[*,*,j] = frag[*,*,j]*fragmask

return, frag

end