function ice = calc_masking_matrices( grid, ice)
% Calculate masking matrices

nrows = grid.nijkuv;
ncols = grid.nijkuv;

nnz_max = grid.nx*grid.ny*grid.nz;

rows_free = zeros( nnz_max,1);
cols_free = zeros( nnz_max,1);
vals_free = zeros( nnz_max,1);
nnz_free  = 0;

rows_surf = zeros( nnz_max,1);
cols_surf = zeros( nnz_max,1);
vals_surf = zeros( nnz_max,1);
nnz_surf  = 0;

rows_base = zeros( nnz_max,1);
cols_base = zeros( nnz_max,1);
vals_base = zeros( nnz_max,1);
nnz_base  = 0;

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);
      
      if k == 1
        % Surface
        
        nnz_surf = nnz_surf + 1;
        rows_surf( nnz_surf) = rcu;
        cols_surf( nnz_surf) = rcu;
        vals_surf( nnz_surf) = 1;
        
        nnz_surf = nnz_surf + 1;
        rows_surf( nnz_surf) = rcv;
        cols_surf( nnz_surf) = rcv;
        vals_surf( nnz_surf) = 1;
        
      elseif k == grid.nz
        % Base
        
        nnz_base = nnz_base + 1;
        rows_base( nnz_base) = rcu;
        cols_base( nnz_base) = rcu;
        vals_base( nnz_base) = 1;
        
        nnz_base = nnz_base + 1;
        rows_base( nnz_base) = rcv;
        cols_base( nnz_base) = rcv;
        vals_base( nnz_base) = 1;
        
      else
        % Free column
        
        nnz_free = nnz_free + 1;
        rows_free( nnz_free) = rcu;
        cols_free( nnz_free) = rcu;
        vals_free( nnz_free) = 1;
        
        nnz_free = nnz_free + 1;
        rows_free( nnz_free) = rcv;
        cols_free( nnz_free) = rcv;
        vals_free( nnz_free) = 1;
        
      end
      
    end
  end
end

ice.M_free = sparse( rows_free( 1:nnz_free), cols_free( 1:nnz_free), vals_free( 1:nnz_free), nrows, ncols);
ice.M_surf = sparse( rows_surf( 1:nnz_surf), cols_surf( 1:nnz_surf), vals_surf( 1:nnz_surf), nrows, ncols);
ice.M_base = sparse( rows_base( 1:nnz_base), cols_base( 1:nnz_base), vals_base( 1:nnz_base), nrows, ncols);

end