bsxfun
is for just this sort of dimension expansion. The trick is that the dimensions need to match exactly, except for the dimensions you want to expand. This means your 3d matrix needs to be instead a 4d time x 1 x lat x lon matrix:
npptot_reshaped = reshape(npptot, [size(npptot,1) 1 size(npptot,2) size(npptot,3)]);
fraction = bsxfun(@rdivide, npp, npptot_reshaped);