function ice = calc_effective_viscosity_Glen( grid, ice)
% Calculate the effective viscosity using Glen's flow law

epsilon_sq_0 = 1e-20;
eta_min      = 1e5;

A = ice.A_Glen;
n = ice.n_Glen;

%% ak-grid

% Calculate the square of the effective strain rate epsilon
epsilon_sq = ice.du_dx_ak.^2 + ...
             ice.dv_dy_ak.^2 + ...
             ice.du_dx_ak .* ice.dv_dy_ak + ...
             0.25 * (ice.du_dy_ak + ice.dv_dx_ak).^2 + ...
             0.25 * (ice.du_dz_ak.^2 + ice.dv_dz_ak.^2) + ...
             epsilon_sq_0;

% Calculate the effective viscosity eta
ice.eta_ak = max( eta_min, 0.5 * A^(-1/n) * (epsilon_sq).^((1-n)/(2*n)));

%% bks-grid

% Calculate the square of the effective strain rate epsilon
epsilon_sq = ice.du_dx_bks.^2 + ...
             ice.dv_dy_bks.^2 + ...
             ice.du_dx_bks .* ice.dv_dy_bks + ...
             0.25 * (ice.du_dy_bks + ice.dv_dx_bks).^2 + ...
             0.25 * (ice.du_dz_bks.^2 + ice.dv_dz_bks.^2) + ...
             epsilon_sq_0;

% Calculate the effective viscosity eta
ice.eta_bks = max( eta_min, 0.5 * A^(-1/n) * (epsilon_sq).^((1-n)/(2*n)));

%% 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
