There is some contrarian thinking reemerging recently with explicit-shape arrays in Fortran which is what you have at the moment. Unfortunately though there is little support for defensive/safe programming against it in the language.
Until that might appear (in a distant future), my personal suggestion will be to follow the route set forth with “modern” Fortran starting with the Fortran 90 standard revision.
Which is to use the combination of explicit interfaces readily achieved via module procedures and also assumed-shape arrays.
Also, to think along the lines of library development as a collection of codes toward a particular “physics” (or math, economics, etc.) problem in what might be a multi"physics" solution/application i.e., many effects get considered together. Then in one or more such collections, it is possible the codes like to work with particular lower bounds that are different from unity. In that case, look at alternate options to achieve such lower bounds such as
- named constants (parameter attribute) at the
MODULElevel, - module entities that have setter methods (and getter if need be),
- parameterized derived types where the length/kind type parameter facilitates lower bound handling via a OO “class” type of approach.
Followed by the use of some run-time assertion approach (may be a simple assert subroutine) to ensure the shapes of arrays in the procedures are as expected: see this, it may already be on the anvil with Fortran stdlib