Use indices from MAX(A,[],dim) to get sum over only part of another array B

Suppose A and B are X x Y x Z arrays, with the 3rd dimension being vertical height. I want to determine where, vertically, A falls below a threshold. Then, I want to do a vertical sum in B ONLY below that height. Currently:

 X = 14; Y = 28; Z = 10;
 A=rand(X,Y,Z); B=rand(X,Y,Z);
 threshold=mean(A(:));
 [~,I] = max(A < threshold ,[],3); 
 I = squeeze(I); % I now has dimensions X x Y 
 for i = 1:size(I,1)
      for j = 1:size(I,2) 
          result(i,j) = nansum(B(i,j,1:I(i,j)),3);
      end
  end

The loop is extremely slow when this method is implemented for very large datasets. There must be a way to vectorize this code! Thank you in advance.

Edit: the process for finding I works well. It is only the loop that needs to be improved.

[Azzi Abdelmalek]

You said I size is X x Y . What is size(I,3)?

[Eowyn Baughman]

Sorry. I was neglecting to mention my actual data has a time dimension. I'll edit the question to be more accurate to the dummy example.