function [A_base, b_base] = calc_stiffness_matrix_BPA_BC_base_twopoint_onesided( grid, ice)
% Calculate stiffness matrix for the basal boundary condition
%
% At the ice base (k=nz), the zero-stress boundary condition implies that:
%
% [1]   2 dh/dx (2 du/dx + dv/dy) + dh/dy (du/dy + dv/dx) - du/dz  + beta/eta u = 0
%
% Discretising du/dz using a two-point, one-sided differencing scheme yields:
%
% [2]   du/dz = dzeta/dz (u( nz) - u( nz-1)) / dzeta
%
% Substituting [2] into [1] yields:
%
% [3]   2 dh/dx (2 du/dx + dv/dy) + dh/dy (du/dy + dv/dx) - dzeta/dz (u( nz) - u( nz-1)) / dzeta + beta/eta u( nz)= 0
% 

% Get dh/dx, dh/dy, and beta/eta in vector form
db_dx_bk         = zeros( grid.nijk, 1);
db_dy_bk         = zeros( grid.nijk, 1);
beta_over_eta_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);
  db_dx_bk(         rc) = ice.db_dx_b( i,j);
  db_dy_bk(         rc) = ice.db_dy_b( i,j);
  beta_over_eta_bk( rc) = ice.basal_friction_coefficient_b( i,j) / ice.eta_bk( rc);
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( db_dx_bk) * (2 * grid.M_ddx_bk_bk * grid.M_map_ijku2ijk + grid.M_ddy_bk_bk * grid.M_map_ijkv2ijk) + ...
          diag_sparse( db_dy_bk) * (    grid.M_ddy_bk_bk * grid.M_map_ijku2ijk + grid.M_ddx_bk_bk * grid.M_map_ijkv2ijk);
        
% -du/dz
Au2 = -1.0 * diag_sparse( ice.dzeta_dz_bk) * (grid.M_u_nz - grid.M_u_nzm1) / grid.Deltazeta * grid.M_map_ijku2ijk;

% beta/eta
Au3 = diag_sparse( beta_over_eta_bk) * grid.M_map_ijku2ijk;
    
% Combine
Au = grid.M_map_ijk2ijku * (Au1 + Au2 + Au3);

%% Assemble stiffness matrix for Eq. 2

% 2 dh/dy ( 2 dv/dy + du/dx ) + dh/dx (   dv/dx + du/dy )
Av1 = 2 * diag_sparse( db_dy_bk) * (2 * grid.M_ddy_bk_bk * grid.M_map_ijkv2ijk + grid.M_ddx_bk_bk * grid.M_map_ijku2ijk) + ...
          diag_sparse( db_dx_bk) * (    grid.M_ddx_bk_bk * grid.M_map_ijkv2ijk + grid.M_ddy_bk_bk * grid.M_map_ijku2ijk);
        
% -dv/dz
Av2 = -1.0 * diag_sparse( ice.dzeta_dz_bk) * (grid.M_u_nz - grid.M_u_nzm1) / grid.Deltazeta * grid.M_map_ijkv2ijk;

% beta/eta
Av3 = diag_sparse( beta_over_eta_bk) * grid.M_map_ijkv2ijk;
    
% Combine
Av = grid.M_map_ijk2ijkv * (Av1 + Av2 + Av3);

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

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

end