The prepro subroutine for BLZ2D seems to have been copy/pasted from BLZ3D and contains discrepancies (notably NNODE = 8 when NNODE should be = 4)
The prepro subroutine needs to be checked and corrected.

Plane or axisymmetric state.
For the axisymmetric state, the symmetry axis is Y. The element is defined by 4 nodes.
Element type: 8
Revised by: Y. Zhu (1991)

Prepro: BLZ2DA.F
Lagamine: BLZ2DB.F

Note: The prepro subroutine for BLZ2D seems to have been copy/pasted from BLZ3D. However, it is highly uncertain that all the features available for BLZ3D are as well available for BLZ2D. The uncertain features are written below in italic.

Title (A5)
TITLE	“BLZ2D” in columns 1 to 5
Control data (2I5)
NELEM	Number of elements
INDPP	0 → Static 1 → Dynamic
INSHE	= 0 for automatic calculation of shear locking parameter
	= 1 if shear coefficient taken into account
ILOAX	= 0 for global axis computation ☛ Objectivity must be verified in the material law (with Jaumann correction) ☛ No rotation of material axes
	< 0 for computation with constant and symetrical velocity gradients pseudo local axes : use of local axes on the time step but no evolution of the local axes on the following time step ☛ Objectivity is verified ☛ No rotation of material axes
	> 0 for computation with local axes ☛ Objectivity is verified ☛ Rotation of material axes
	units: = 1 for rotations incorporated in local tangent matrix Not available = 2 apply final rotation to local tangent matrix = 3 apply initial rotation to local tangent matrix = 4 compute tangent matrix through global perturbation method
	tens (only for ILOAX>0): = 0 for local axes e1, e2, e3 initially parallel to global axes ex, ey, ez = 1 for local axes e1, e2 given (and e3=e1∧e2) = 2 for local axes e1, e2 initially in the plane (ex, ey) forming an angle θ with ex, ey (and e3=e1∧e2) = 3 same as 1 with different local axes for each element = 4 same as 2 with different local axes for each element
INSIG	0 → No initial stresses 1 → $\sigma_Y=\sigma_{Y0}+yd\sigma_Y$ 3 → See PLXLS
Specific weight - Only if INDPP = 1 (3G10.0)
WSPE(1)	= specific weight in X direction
WSPE(2)	= specific weight in Y direction
WSPE(3)	= density
Shear coefficient - Only if INSHE > 0 (G10.0)
PARSHE	Shear locking coefficient ∈ [0,1] - close to 0: avoid shear locking but higher risk of hourglass modes (use for thin elements in flexion) - close to 1: avoid hourglass modes but higher risk of shear locking (use for cubic elements in shear)
Definition of the elements (2I5/4I5)
NINTE	Number of integration points (1, 2, or 4)
LMATE	Number of the material law
NODES(4)	List of nodes
SIG0(6)	List of initial stresses (Only if INSIG ≠ 0) If INSIG = 1 (3G10.0): $\sigma_{Y0}$ - Effective stress $\sigma_Y$ at the axes origin DSIGY - Gradient of effective stress along axis OY AK0 - ratio $\sigma_X/\sigma_Y$ If INSIG = 3 (6G10.0): See PLXLS