function [A_base, b_base] = calc_stiffness_matrix_BPA_BC_base_ghostpoint( grid, ice)
% Calculate stiffness matrix for the basal boundary condition

% 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

% d/dx ( 2 eta ( 2 du/dx + dv/dy ))
Au1 = grid.M_ddx_ak_bk * 2 * diag_sparse( ice.eta_ak) * (2 * grid.M_ddx_bk_ak * grid.M_map_ijku2ijk + grid.M_ddy_bk_ak * grid.M_map_ijkv2ijk);
% d/dy (   eta (   du/dy + dv/dx ))
Au2 = grid.M_ddy_ak_bk *     diag_sparse( ice.eta_ak) * (    grid.M_ddy_bk_ak * grid.M_map_ijku2ijk + grid.M_ddx_bk_ak * grid.M_map_ijkv2ijk);

% 2 db/dx ( 2 du/dx + dv/dy ) + db/dy (   du/dy + dv/dx ) + beta / eta u
BCu = 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) + ...
          diag_sparse( beta_over_eta_bk) * grid.M_map_ijku2ijk;

% deta/dz 2 db/dx ( 2 du/dx + dv/dy ) + db/dy (   du/dy + dv/dx ) + beta / eta u
Au3 = diag_sparse( ice.deta_dz_bk) * BCu;
        
% (dzeta/dz)^2 2 eta / Deltazeta^2 [u( nz-1) - u( nz) + Deltazeta/(dzeta/dz) ( 2 db/dx (..) + db/dy (..) + beta / eta u)]
Au4 = diag_sparse( 2 * ice.dzeta_dz_bk.^2 .* ice.eta_bk) / grid.Deltazeta^2 * ((grid.M_u_nzm1 - grid.M_u_nz) * grid.M_map_ijku2ijk + ...
      diag_sparse( grid.Deltazeta ./ ice.dzeta_dz_bk) * BCu);
    
% Combine
Au = grid.M_map_ijk2ijku * (Au1 + Au2 + Au3 + Au4);

%% Assemble stiffness matrix for Eq. 2

% d/dy ( 2 eta ( 2 dv/dy + du/dx ))
Av1 = grid.M_ddy_ak_bk * 2 * diag_sparse( ice.eta_ak) * (2 * grid.M_ddy_bk_ak * grid.M_map_ijkv2ijk + grid.M_ddx_bk_ak * grid.M_map_ijku2ijk);
% d/dx (   eta (   dv/dx + du/dy ))
Av2 = grid.M_ddx_ak_bk *     diag_sparse( ice.eta_ak) * (    grid.M_ddx_bk_ak * grid.M_map_ijkv2ijk + grid.M_ddy_bk_ak * grid.M_map_ijku2ijk);

% 2 db/dy ( 2 dv/dy + du/dx ) + db/dx (   dv/dx + du/dy ) + beta / eta v
BCv = 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) + ...
          diag_sparse( beta_over_eta_bk) * grid.M_map_ijkv2ijk;

% deta/dz 2 db/dy ( 2 dv/dy + du/dx ) + db/dx (   dv/dx + du/dy ) + beta / eta v
Av3 = diag_sparse( ice.deta_dz_bk) * BCv;
        
% (dzeta/dz)^2 2 eta / Deltazeta^2 [v( nz-1) - v( nz) + Deltazeta/(dzeta/dz) ( 2 dh/dy (..) + dh/dx (..) + beta / eta v)]
Av4 = diag_sparse( 2 * ice.dzeta_dz_bk.^2 .* ice.eta_bk) / grid.Deltazeta^2 * ((grid.M_u_nzm1 - grid.M_u_nz) * grid.M_map_ijkv2ijk + ...
      diag_sparse( grid.Deltazeta ./ ice.dzeta_dz_bk) * BCv);
    
% Combine
Av = grid.M_map_ijk2ijkv * (Av1 + Av2 + Av3 + Av4);

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

% Fill load vector
b_base = zeros( grid.nijkuv, 1);
for i = 1: grid.nx
for j = 1: grid.ny
for k = 1: grid.nz
  rcu = grid.ijkuv2rc( i,j,k,1);
  rcv = grid.ijkuv2rc( i,j,k,2);
  b_base( rcu) = ice.rho * ice.g * ice.db_dx_b( i,j);
  b_base( rcv) = ice.rho * ice.g * ice.db_dy_b( i,j);
end
end
end

end