Modify quadrature used for ice shelf viscosity (#468)

Previously, when ice_viscosity_compute == MODEL, ice shelf viscosity was calculated using 1 quadrature point per cell; however, this quadrature point was not centered in the cell.
-Added a routine bilinear_shape_fn_grid_1qp, which is used to instead calculate viscosity for a single cell-centered quadrature point
-Added an option to define parameter ice_viscosity_compute = MODEL_QUADRATURE, where viscosity is calculated at the same four quadrature points used during the SSA solution (the typical approach in finite element codes).

Note that when using one quadrature point (ice_viscosity_compute == MODEL), array ice_visc(:,:) is the cell-centered ice viscosity. When using four quadrature points (ice_viscosity_compute == MODEL_QUADRATURE), ice_visc(:,:) is the cell-centered ice viscosity divided by the effective stress, Ee, which varies between each quadrature point and is saved in a separate array CS%Ee(:,:,:). In the SSA, ice viscosity is calculated for a cell with indices i,j as ice_visc(i,j) * Ee, where Ee=1 if using one quadrature point for viscosity and Ee=CS%Ee(i,j,k) if using four quad points (where k the current quad point). Also note the post_data call for ice_visc when using four quad points, where Ice_visc is outputted for visualization from a cell as ice_visc(i,j)*Ee_av(i,j), where Ee_av(i,j) is the average Ee in the cell.