Improve speed
Currently the speed of MPR fell far behind that of its previous version that was part of the mHM code (ca. 50x for test basin). There are some possible reasons for that:
 the soil parameters (now 50% of all calculations) of the previous mHM were determined based on a lookuptable with soil classes and their properties and determined for existing classes only (read in of integer ascii grid)
 the (spatial) upscaling did not rely on weights, lowres cell boundaries needed to align with highres boundaries
 2Dupscaling was done in one go, other dimensions did not have missing values and were always at first positions in order of dimensions
which contrasts the new approach:
 we now do that fully distributed without prior classification (this also means, we read float arrays)
 we allow for nonaligned cell boundaries in high and lowresolution
 we allow for any order of dimensions and also allow missing values within all dimensions (e.g. a missing value in a soil horizon)
So the idea is to check the current code for bottlenecks and identify alternative algorithms. After that we might think about shortcuts that can be used in cases when data fulfill the requirements of mHM data (e.g. dimension order (x,y, remaining dim) and no missing values within a x,yslice). I was wondering, how @schaefed might help us here?!
The scheme is currently as follows for the example of the case where TargetVar = SourceVar1 + SourceVar2

read in the SourceVar1 and SourceVar2

store the nDimmask in a flat array, store the data a flat, packed array, store the dimensions

calculate temporary TargetVar (new data, use mask and dimension of SourceVars)

check for dimension order of TargetVar, if different than that of SourceVar, transpose temporary TargetVar

upscale each dimension individually (loop over the target dimensions)
 check if source dimension and target dimension differ (case1: no difference  continue, case2: source dim not existing  broadcast dimensions, case3: difference  upscale the dimension (see b)))
 compare dimensions, for each target dim id: calculate vector of source dim ids, vector of weights for source dim ids, number of subcells (if same pair of dimensions were used before, skip and use previous information)
 reinsert missing data in 1D flat data array (so has same shape with 1D flat mask), initialize empty missing_value_weight array in same shape
 loop over the n data chunks (n is the product of dimension sizes after the current dimension)
 reshape to 2D and transpose each chunk (new shape is (size of dimension, product of sizes after the before current dimension))
 loop now over 2nd dim, resulting vector is then a slice with the length of our current dimension
 if this vector contains nans, recalculate basic dimension weights, ids and subcells (see 5b) and missing_value_weight h) loop over target cells and apply upscaling function, yield new data