well, i did not see carefully you did a cross reference with external software. called as surface traction, a force magnitude applied to group of face element. as i know FreeCAD also doing the same feature approach, you may post a specific question at they forums.
it seems the forces are automatically converted or applied as concentrated load with equivalent magnitude depending on tributary areas of a nodes and element type (especially in quadratic). summation are required for element face sharing the same nodes. so, there’s no modification of user subroutine of dload.f as you describe, also i did not see any hints for these task in user defined loading documentations.
rather than took all above long steps, alternatively one can use distributing coupling features or rigid body keywords if additional stiffness can deal to be ignores.