function ice = calc_zeta_gradients( grid, ice)
% Calculate gradients of zeta

% Calculate horizontal gradients of Hi and Hs
dHi_dx_a = ddx_ax_ax_2D( grid, ice.Hi_a);
dHi_dy_a = ddy_ay_ay_2D( grid, ice.Hi_a);
dHs_dx_a = ddx_ax_ax_2D( grid, ice.Hs_a);
dHs_dy_a = ddy_ay_ay_2D( grid, ice.Hs_a);

Hi_cx    = map_ax_cx_2D( grid, ice.Hi_a);
Hi_cy    = map_ay_cy_2D( grid, ice.Hi_a);
Hi_b     = map_ax_cx_2D( grid, Hi_cy);
Hs_cx    = map_ax_cx_2D( grid, ice.Hs_a);
Hs_cy    = map_ay_cy_2D( grid, ice.Hs_a);

dHi_dx_b = ddx_ax_cx_2D( grid, Hi_cy);
dHi_dy_b = ddy_ay_cy_2D( grid, Hi_cx);
dHs_dx_b = ddx_ax_cx_2D( grid, Hs_cy);
dHs_dy_b = ddy_ay_cy_2D( grid, Hs_cx);

%% On the ak-grid

for i = 1: grid.nx
  for j = 1: grid.ny
    
    Hi = ice.Hi_a( i,j);
    
    dHi_dx = dHi_dx_a( i,j);
    dHs_dx = dHs_dx_a( i,j);
    dHi_dy = dHi_dy_a( i,j);
    dHs_dy = dHs_dy_a( i,j);
    
    for k = 1: grid.nz
      rc = grid.ijk2rc( i,j,k);
      ice.dzeta_dx_ak( rc) = ( 1 / Hi) * (dHs_dx - grid.zeta( k) * dHi_dx);
      ice.dzeta_dy_ak( rc) = ( 1 / Hi) * (dHs_dy - grid.zeta( k) * dHi_dy);
      ice.dzeta_dz_ak( rc) = (-1 / Hi);
    end
    
  end
end

%% On the bk-grid

for i = 1: grid.nx
  for j = 1: grid.ny
    
    Hi = Hi_b( i,j);
    
    dHi_dx = dHi_dx_b( i,j);
    dHs_dx = dHs_dx_b( i,j);
    dHi_dy = dHi_dy_b( i,j);
    dHs_dy = dHs_dy_b( i,j);
    
    for k = 1: grid.nz
      rc = grid.ijk2rc( i,j,k);
      ice.dzeta_dx_bk( rc) = ( 1 / Hi) * (dHs_dx - grid.zeta( k) * dHi_dx);
      ice.dzeta_dy_bk( rc) = ( 1 / Hi) * (dHs_dy - grid.zeta( k) * dHi_dy);
      ice.dzeta_dz_bk( rc) = (-1 / Hi);
    end
    
  end
end

%% On the bks-grid

for i = 1: grid.nx
  for j = 1: grid.ny
    
    Hi = Hi_b( i,j);
    
    dHi_dx = dHi_dx_b( i,j);
    dHs_dx = dHs_dx_b( i,j);
    dHi_dy = dHi_dy_b( i,j);
    dHs_dy = dHs_dy_b( i,j);
    
    for ks = 1: grid.nz-1
      rc = grid.ijks2rc( i,j,ks);
      ice.dzeta_dx_bks( rc) = ( 1 / Hi) * (dHs_dx - grid.zeta_stag( ks) * dHi_dx);
      ice.dzeta_dy_bks( rc) = ( 1 / Hi) * (dHs_dy - grid.zeta_stag( ks) * dHi_dy);
      ice.dzeta_dz_bks( rc) = (-1 / Hi);
    end
    
  end
end

end