function grid = calc_gradient_operators_scaled( grid)
% Calculate gradient operators in [x_hat,y_hat,zeta]

% Calculate the basic mapping and gradient operations between the different
% staggered grids; all other operations can be expressed as combinations of
% these.
grid = calc_gradient_operators_scaled_x( grid);
grid = calc_gradient_operators_scaled_y( grid);
grid = calc_gradient_operators_scaled_z( grid);

% Special ones for the surface and basal boundary conditions
grid = calc_surface_operators( grid);
grid = calc_basal_operators(   grid);

% Operators between different grids
grid.M_ddx_hat_ak_bk  = grid.M_map_ay_cy * grid.M_ddx_hat_ax_cx;
grid.M_ddy_hat_ak_bk  = grid.M_map_ax_cx * grid.M_ddy_hat_ay_cy;
grid.M_ddx_hat_bk_ak  = grid.M_map_cy_ay * grid.M_ddx_hat_cx_ax;
grid.M_ddy_hat_bk_ak  = grid.M_map_cx_ax * grid.M_ddy_hat_cy_ay;
grid.M_ddx_hat_bk_bk  =                    grid.M_ddx_hat_ax_ax;
grid.M_ddy_hat_bk_bk  =                    grid.M_ddy_hat_ay_ay;

grid.M_map_ak_bk      = grid.M_map_ax_cx * grid.M_map_ay_cy;
grid.M_map_bk_ak      = grid.M_map_cx_ax * grid.M_map_cy_ay;
grid.M_map_ak_bks     = grid.M_map_az_cz * grid.M_map_ak_bk;
grid.M_map_bks_ak     = grid.M_map_bk_ak * grid.M_map_cz_az;

  function grid = calc_gradient_operators_scaled_x( grid)
    
    %% regular to regular (or staggered to staggered)
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 2;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        for k = 1: grid.nz
          
          i_left = i-1;
          if i_left == 0; i_left = grid.nx; end
          i_right = i+1;
          if i_right == grid.nx+1; i_right = 1; end
          
          % left-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i_left,j,k);
          vals( n) = -1 / (2 * grid.Deltax);
          
          % right-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i_right,j,k);
          vals( n) = 1 / (2 * grid.Deltax);
      
        end
      end
    end

    % Create sparse matrices from lists
    grid.M_ddx_hat_ax_ax = sparse( rows, cols, vals, nrows, ncols);

    %% staggered to regular
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 2;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals_map = zeros( nnz,1);
    vals_ddx = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        for k = 1: grid.nz
          
          i_left = i-1;
          if i_left == 0; i_left = grid.nx; end
          i_right = i;
          
          % left-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i_left,j,k);
          vals_map( n) = 1/2;
          vals_ddx( n) = -1 / grid.Deltax;
          
          % right-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i_right,j,k);
          vals_map( n) = 1/2;
          vals_ddx( n) = 1 / grid.Deltax;
      
        end
      end
    end

    % Create sparse matrices from lists
    grid.M_map_cx_ax     = sparse( rows, cols, vals_map, nrows, ncols);
    grid.M_ddx_hat_cx_ax = sparse( rows, cols, vals_ddx, nrows, ncols);

    %% regular to staggered
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 2;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals_map = zeros( nnz,1);
    vals_ddx = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        for k = 1: grid.nz
          
          i_left = i;
          i_right = i+1;
          if i_right == grid.nx+1; i_right = 1; end
          
          % left-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i_left,j,k);
          vals_map( n) = 1/2;
          vals_ddx( n) = -1 / grid.Deltax;
          
          % right-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i_right,j,k);
          vals_map( n) = 1/2;
          vals_ddx( n) = 1 / grid.Deltax;
      
        end
      end
    end

    % Create sparse matrices from lists
    grid.M_map_ax_cx     = sparse( rows, cols, vals_map, nrows, ncols);
    grid.M_ddx_hat_ax_cx = sparse( rows, cols, vals_ddx, nrows, ncols);

  end
  function grid = calc_gradient_operators_scaled_y( grid)
    
    %% regular to regular (or staggered to staggered)
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 2;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        for k = 1: grid.nz
          
          j_left = j-1;
          if j_left == 0; j_left = grid.ny; end
          j_right = j+1;
          if j_right == grid.ny+1; j_right = 1; end
          
          % left-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i,j_left,k);
          vals( n) = -1 / (2 * grid.Deltay);
          
          % right-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i,j_right,k);
          vals( n) = 1 / (2 * grid.Deltay);
      
        end
      end
    end

    % Create sparse matrices from lists
    grid.M_ddy_hat_ay_ay = sparse( rows, cols, vals, nrows, ncols);

    %% staggered to regular
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 2;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals_map = zeros( nnz,1);
    vals_ddy = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        for k = 1: grid.nz
          
          j_left = j-1;
          if j_left == 0; j_left = grid.ny; end
          j_right = j;
          
          % left-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i,j_left,k);
          vals_map( n) = 1/2;
          vals_ddy( n) = -1 / grid.Deltay;
          
          % right-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i,j_right,k);
          vals_map( n) = 1/2;
          vals_ddy( n) = 1 / grid.Deltay;
      
        end
      end
    end

    % Create sparse matrices from lists
    grid.M_map_cy_ay     = sparse( rows, cols, vals_map, nrows, ncols);
    grid.M_ddy_hat_cy_ay = sparse( rows, cols, vals_ddy, nrows, ncols);

    %% regular to staggered
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 2;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals_map = zeros( nnz,1);
    vals_ddy = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        for k = 1: grid.nz
          
          j_left = j;
          j_right = j+1;
          if j_right == grid.ny+1; j_right = 1; end
          
          % left-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i,j_left,k);
          vals_map( n) = 1/2;
          vals_ddy( n) = -1 / grid.Deltay;
          
          % right-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i,j_right,k);
          vals_map( n) = 1/2;
          vals_ddy( n) = 1 / grid.Deltay;
      
        end
      end
    end

    % Create sparse matrices from lists
    grid.M_map_ay_cy     = sparse( rows, cols, vals_map, nrows, ncols);
    grid.M_ddy_hat_ay_cy = sparse( rows, cols, vals_ddy, nrows, ncols);

  end
  function grid = calc_gradient_operators_scaled_z( grid)
    
    %% regular to regular (or staggered to staggered)
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 2;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        for k = 1: grid.nz
          
          if k==1
            k_left  = 1;
            k_right = 2;
            v       = 1 / grid.Deltazeta;
          elseif k == grid.nz
            k_left  = grid.nz - 1;
            k_right = grid.nz;
            v       = 1 / grid.Deltazeta;
          else
            k_left  = k-1;
            k_right = k+1;
            v       = 1 / (2 * grid.Deltazeta);
          end
          
          % left-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i,j,k_left);
          vals( n) = -v;
          
          % right-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i,j,k_right);
          vals( n) = v;
      
        end
      end
    end

    % Create sparse matrices from lists
    grid.M_ddzeta_az_az = sparse( rows, cols, vals, nrows, ncols);
    
    %% regular to regular (or staggered to staggered)
    % three-point one-sided scheme at the surface
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 2;
    
    nnz = nrows * nnz_per_row + 1;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        for k = 1: grid.nz
          
          if k==1
            % 1
            n = n+1;
            rows( n) = grid.ijk2rc( i,j,1);
            cols( n) = grid.ijk2rc( i,j,1);
            vals( n) = -3/2 / grid.Deltazeta;
            % 2
            n = n+1;
            rows( n) = grid.ijk2rc( i,j,1);
            cols( n) = grid.ijk2rc( i,j,2);
            vals( n) = 2 / grid.Deltazeta;
            % 3
            n = n+1;
            rows( n) = grid.ijk2rc( i,j,1);
            cols( n) = grid.ijk2rc( i,j,3);
            vals( n) = -1/2 / grid.Deltazeta;
            continue
          elseif k == grid.nz
            k_left  = grid.nz - 1;
            k_right = grid.nz;
            v       = 1 / grid.Deltazeta;
          else
            k_left  = k-1;
            k_right = k+1;
            v       = 1 / (2 * grid.Deltazeta);
          end
          
          % left-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i,j,k_left);
          vals( n) = -v;
          
          % right-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijk2rc( i,j,k_right);
          vals( n) = v;
      
        end
      end
    end

    % Create sparse matrices from lists
    grid.M_ddzeta_az_az_threepoint_onesided = sparse( rows, cols, vals, nrows, ncols);

    %% staggered to regular
    nrows       = grid.nijk;
    ncols       = grid.nijks;
    nnz_per_row = 2;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals_map = zeros( nnz,1);
    vals_ddz = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        for k = 1: grid.nz
          
          if k==1
            ks_left  = 1;
            ks_right = 2;
          elseif k == grid.nz
            ks_left  = grid.nz - 2;
            ks_right = grid.nz - 1;
          else
            ks_left  = k-1;
            ks_right = k;
          end
          
          % left-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijks2rc( i,j,ks_left);
          vals_map( n) = 1/2;
          vals_ddz( n) = -1 / grid.Deltazeta;
          
          % right-hand neighbour
          n = n+1;
          rows( n) = grid.ijk2rc( i,j,k);
          cols( n) = grid.ijks2rc( i,j,ks_right);
          vals_map( n) = 1/2;
          vals_ddz( n) = 1 / grid.Deltazeta;
      
        end
      end
    end

    % Create sparse matrices from lists
    grid.M_map_cz_az    = sparse( rows, cols, vals_map, nrows, ncols);
    grid.M_ddzeta_cz_az = sparse( rows, cols, vals_ddz, nrows, ncols);

    %% regular to staggered
    nrows       = grid.nijks;
    ncols       = grid.nijk;
    nnz_per_row = 2;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals_map = zeros( nnz,1);
    vals_ddz = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        for ks = 1: grid.nz-1
          
          k_left  = ks;
          k_right = ks+1;
          
          % left-hand neighbour
          n = n+1;
          rows( n) = grid.ijks2rc( i,j,ks);
          cols( n) = grid.ijk2rc( i,j,k_left);
          vals_map( n) = 1/2;
          vals_ddz( n) = -1 / grid.Deltazeta;
          
          % right-hand neighbour
          n = n+1;
          rows( n) = grid.ijks2rc( i,j,ks);
          cols( n) = grid.ijk2rc( i,j,k_right);
          vals_map( n) = 1/2;
          vals_ddz( n) = 1 / grid.Deltazeta;
      
        end
      end
    end

    % Create sparse matrices from lists
    grid.M_map_az_cz    = sparse( rows, cols, vals_map, nrows, ncols);
    grid.M_ddzeta_az_cz = sparse( rows, cols, vals_ddz, nrows, ncols);

  end

  function grid = calc_surface_operators( grid)
    % Special ones for the surface boundary condition
    
    %% u( 1)
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 1;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        n = n+1;
        rows( n) = grid.ijk2rc( i,j,1);
        cols( n) = grid.ijk2rc( i,j,1);
        vals( n) = 1;
      end
    end

    % Create sparse matrices from lists
    grid.M_u_1 = sparse( rows( 1:n), cols( 1:n), vals( 1:n), nrows, ncols);
    
    %% u( 2)
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 1;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        n = n+1;
        rows( n) = grid.ijk2rc( i,j,1);
        cols( n) = grid.ijk2rc( i,j,2);
        vals( n) = 1;
      end
    end

    % Create sparse matrices from lists
    grid.M_u_2 = sparse( rows( 1:n), cols( 1:n), vals( 1:n), nrows, ncols);
    
    %% u( 3)
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 1;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        n = n+1;
        rows( n) = grid.ijk2rc( i,j,1);
        cols( n) = grid.ijk2rc( i,j,3);
        vals( n) = 1;
      end
    end

    % Create sparse matrices from lists
    grid.M_u_3 = sparse( rows( 1:n), cols( 1:n), vals( 1:n), nrows, ncols);
    
    %% u( 4)
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 1;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        n = n+1;
        rows( n) = grid.ijk2rc( i,j,1);
        cols( n) = grid.ijk2rc( i,j,4);
        vals( n) = 1;
      end
    end

    % Create sparse matrices from lists
    grid.M_u_4 = sparse( rows( 1:n), cols( 1:n), vals( 1:n), nrows, ncols);
    
    %% u( 5)
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 1;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        n = n+1;
        rows( n) = grid.ijk2rc( i,j,1);
        cols( n) = grid.ijk2rc( i,j,5);
        vals( n) = 1;
      end
    end

    % Create sparse matrices from lists
    grid.M_u_5 = sparse( rows( 1:n), cols( 1:n), vals( 1:n), nrows, ncols);
    
  end
  function grid = calc_basal_operators(   grid)
    % Special ones for the basal boundary condition
    
    %% u( nz)
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 1;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        n = n+1;
        rows( n) = grid.ijk2rc( i,j,grid.nz  );
        cols( n) = grid.ijk2rc( i,j,grid.nz  );
        vals( n) = 1;
      end
    end

    % Create sparse matrices from lists
    grid.M_u_nz = sparse( rows( 1:n), cols( 1:n), vals( 1:n), nrows, ncols);
    
    %% u( nz-1)
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 1;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        n = n+1;
        rows( n) = grid.ijk2rc( i,j,grid.nz  );
        cols( n) = grid.ijk2rc( i,j,grid.nz-1);
        vals( n) = 1;
      end
    end

    % Create sparse matrices from lists
    grid.M_u_nzm1 = sparse( rows( 1:n), cols( 1:n), vals( 1:n), nrows, ncols);
    
    %% u( nz-2)
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 1;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        n = n+1;
        rows( n) = grid.ijk2rc( i,j,grid.nz  );
        cols( n) = grid.ijk2rc( i,j,grid.nz-2);
        vals( n) = 1;
      end
    end

    % Create sparse matrices from lists
    grid.M_u_nzm2 = sparse( rows( 1:n), cols( 1:n), vals( 1:n), nrows, ncols);
    
    %% u( nz-3)
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 1;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        n = n+1;
        rows( n) = grid.ijk2rc( i,j,grid.nz  );
        cols( n) = grid.ijk2rc( i,j,grid.nz-3);
        vals( n) = 1;
      end
    end

    % Create sparse matrices from lists
    grid.M_u_nzm3 = sparse( rows( 1:n), cols( 1:n), vals( 1:n), nrows, ncols);
    
    %% u( nz-4)
    nrows       = grid.nijk;
    ncols       = grid.nijk;
    nnz_per_row = 1;
    
    nnz = nrows * nnz_per_row;
    rows = zeros( nnz,1);
    cols = zeros( nnz,1);
    vals = zeros( nnz,1);
    
    n = 0;
    for i = 1: grid.nx
      for j = 1: grid.ny
        n = n+1;
        rows( n) = grid.ijk2rc( i,j,grid.nz  );
        cols( n) = grid.ijk2rc( i,j,grid.nz-2);
        vals( n) = 1;
      end
    end

    % Create sparse matrices from lists
    grid.M_u_nzm4 = sparse( rows( 1:n), cols( 1:n), vals( 1:n), nrows, ncols);
    
  end

end