This study presents a nonlinear viscoelastic material model incorporating a progressive damage framework with an iterative algorithm for glacial/freshwater polycrystalline ice subject to compressive impact load induced by ice-structure interaction. The damage model accounts for microcracking, dynamic recrystallisation, pressure melting, and high-shear elastic failure with a pressure- and rate-dependent convex failure locus. The constitutive laws are written in Fortran and implemented as a vectorised user material (VUMAT) in Abaqus with three different numerical methods, Lagrangian FEM, coupled FEM-SPH, and ALE-FEM. The constitutive model together with the implemented numerical methods are validated against two different types of laboratory-scale physical tests, indentation of cone-shaped ice and triaxial creep. The proposed model implemented with the coupled FEM-SPH method enables simulation of the cyclic transition from solid-like intact ice to the fluid-like pulverised/granular substance, progressively extruded with viscous rheology.