function [A_surf, b_surf] = calc_stiffness_matrix_BPA_BC_surf_fivepoint_onesided( grid, ice)
% Calculate stiffness matrix for the surface boundary condition
%
% At the ice surface (k=1), the zero-stress boundary condition implies that:
%
% [1]   2 dh/dx (2 du/dx + dv/dy) + dh/dy (du/dy + dv/dx) - du/dz = 0
%
% Discretising du/dz using a five-point, one-sided differencing scheme yields:
%
% [2]   du/dz = dzeta/dz (-25/12 u( 1) + 48/12 u( 2) - 36/12 u( 3) + 16/12 u( 4) - 3/12 u( 5)) / dzeta
%
% Substituting [2] into [1] yields:
%
% [3]   2 dh/dx (2 du/dx + dv/dy) + dh/dy (du/dy + dv/dx) - dzeta/dz (-11/6 u( 1) + 18/6 u( 2) - 9/6 u( 3) + 2/6 u( 4)) / dzeta = 0
% 

% Get dh/dx, dh/dy in vector form
dh_dx_bk = zeros( grid.nijk, 1);
dh_dy_bk = zeros( grid.nijk, 1);

for i = 1: grid.nx
for j = 1: grid.ny
for k = 1: grid.nz
  rc = grid.ijk2rc( i,j,k);
  dh_dx_bk( rc) = ice.dh_dx_b( i,j);
  dh_dy_bk( rc) = ice.dh_dy_b( i,j);
end
end
end

%% Assemble stiffness matrix for Eq. 1

% 2 dh/dx ( 2 du/dx + dv/dy ) + dh/dy (   du/dy + dv/dx )
Au1 = 2 * diag_sparse( dh_dx_bk) * (2 * grid.M_ddx_bk_bk * grid.M_map_ijku2ijk + grid.M_ddy_bk_bk * grid.M_map_ijkv2ijk) + ...
          diag_sparse( dh_dy_bk) * (    grid.M_ddy_bk_bk * grid.M_map_ijku2ijk + grid.M_ddx_bk_bk * grid.M_map_ijkv2ijk);
        
% -du/dz
Au2 = diag_sparse( ice.dzeta_dz_bk) * (...
  -25/12 * grid.M_u_1 + ...
   48/12 * grid.M_u_2 + ...
  -36/12 * grid.M_u_3 + ...
   16/12 * grid.M_u_4 + ...
   -3/12 * grid.M_u_5) / grid.Deltazeta * grid.M_map_ijku2ijk;
    
% Combine
Au = grid.M_map_ijk2ijku * (Au1 + Au2);

%% Assemble stiffness matrix for Eq. 2

% 2 dh/dy ( 2 dv/dy + du/dx ) + dh/dx (   dv/dx + du/dy )
Av1 = 2 * diag_sparse( dh_dy_bk) * (2 * grid.M_ddy_bk_bk * grid.M_map_ijkv2ijk + grid.M_ddx_bk_bk * grid.M_map_ijku2ijk) + ...
          diag_sparse( dh_dx_bk) * (    grid.M_ddx_bk_bk * grid.M_map_ijkv2ijk + grid.M_ddy_bk_bk * grid.M_map_ijku2ijk);
        
% -dv/dz
Av2 = diag_sparse( ice.dzeta_dz_bk) * (...
  -25/12 * grid.M_u_1 + ...
   48/12 * grid.M_u_2 + ...
  -36/12 * grid.M_u_3 + ...
   16/12 * grid.M_u_4 + ...
   -3/12 * grid.M_u_5) / grid.Deltazeta * grid.M_map_ijkv2ijk;
    
% Combine
Av = grid.M_map_ijk2ijkv * (Av1 + Av2);

%% Combine Eqs. 1 and 2
A_surf = Au + Av;

% Fill load vector
b_surf = zeros( grid.nijkuv, 1);

end