function grid = calc_gradient_operators_xyz( grid, ice)
% Calculate gradient operators in [x,y,z] (i.e. geometry-dependent operators)

% d/dx: ak-bk, bk-ak, bk-bk
grid.M_ddx_ak_bk  = grid.M_ddx_hat_ak_bk + diag_sparse( ice.dzeta_dx_bk ) * grid.M_ddzeta_az_az;
grid.M_ddx_bk_ak  = grid.M_ddx_hat_bk_ak + diag_sparse( ice.dzeta_dx_ak ) * grid.M_ddzeta_az_az;
grid.M_ddx_bk_bk  = grid.M_ddx_hat_bk_bk + diag_sparse( ice.dzeta_dx_bk ) * grid.M_ddzeta_az_az;

% d/dy: ak-bk, bk-ak, bk-bk
grid.M_ddy_ak_bk  = grid.M_ddy_hat_ak_bk + diag_sparse( ice.dzeta_dy_bk ) * grid.M_ddzeta_az_az;
grid.M_ddy_bk_ak  = grid.M_ddy_hat_bk_ak + diag_sparse( ice.dzeta_dy_ak ) * grid.M_ddzeta_az_az;
grid.M_ddy_bk_bk  = grid.M_ddy_hat_bk_bk + diag_sparse( ice.dzeta_dy_bk ) * grid.M_ddzeta_az_az;

% d/dz: bk-bks, bks-bk, bk-bk
grid.M_ddz_bk_bks =                        diag_sparse( ice.dzeta_dz_bks) * grid.M_ddzeta_az_cz;
grid.M_ddz_bks_bk =                        diag_sparse( ice.dzeta_dz_bk ) * grid.M_ddzeta_cz_az;
grid.M_ddz_bk_bk  =                        diag_sparse( ice.dzeta_dz_bk ) * grid.M_ddzeta_az_az;
grid.M_ddz_bk_bk_threepoint_onesided =     diag_sparse( ice.dzeta_dz_bk ) * grid.M_ddzeta_az_az_threepoint_onesided;

end
