function ice = calc_effective_viscosity_linearvisc( grid, ice)
% Calculate the effective viscosity using the artificial flow law

%% ak-grid
for i = 1: grid.nx
  for j = 1: grid.ny
    for k = 1: grid.nz
      zp = 1500 - grid.zeta( k) * ice.Hi_a( i,j);
      rc = grid.ijk2rc( i,j,k);
      ice.eta_ak( rc) = ice.a + zp;
    end
  end
end

%% bks-grid
Hi_b = map_ay_cy_2D( grid, map_ax_cx_2D( grid, ice.Hi_a));
for i = 1: grid.nx
  for j = 1: grid.ny
    for ks = 1: grid.nz-1
      zp = 1500 - grid.zeta_stag( ks) * Hi_b( i,j);
      rc = grid.ijks2rc( i,j,ks);
      ice.eta_bks( rc) = ice.a + zp;
    end
  end
end

%% Map eta from the ak-grid to the bk-grid
ice.eta_bk     = grid.M_map_ak_bk  * ice.eta_ak;
ice.deta_dz_bk = grid.M_ddz_bks_bk * ice.eta_bks;

end
