docs: def _make_meridional_diffusion_matrix(K, lataxis):
      """Calls :func:`_make_diffusion_matrix` with appropriate weights for
      the meridional diffusion case.
  
      :param array K:         dimensionless diffusivities at cell boundaries
                              of diffusion axis ``lataxis``
      :param axis lataxis:    latitude axis where diffusion is occuring
  
      Weights are computed as the following:
  
      .. math::
  
          \\begin{array}{ll}
              w_1 &= \\cos(\\textrm{bounds}) \\\\
                  &= \\left[ \\cos(b_0), \\cos(b_1), \\cos(b_2), \\ ... \\ , \\cos(b_{n-1}), \\cos(b_n) \\right] \\\\
              w_2 &= \\cos(\\textrm{points}) \\\\
                  &= \\left[ \\cos(p_0), \\cos(p_1), \\cos(p_2), \\ ... \\ , \\cos(p_{n-1}) \\right]
          \\end{array}
  
      when bounds and points from ``lataxis`` are written as
  
      .. math::
  
          \\begin{array}{ll}
              \\textrm{bounds}   &= [b_0, b_1, b_2, \\ ... \\ , b_{n-1}, b_{n}] \\\\
              \\textrm{points}   &= [p_0, p_1, p_2, \\ ... \\ , p_{n-1}]
          \\end{array}
  
      Giving this input to :func:`_make_diffusion_matrix` results in a matrix like:
  
      .. math::
  
          \\textrm{diffTriDiag}=
          \\left[ \\begin{array}{cccccc}
          1+\\frac{u_1 }{\\cos(p_0)} & -\\frac{u_1}{\\cos(p_0)} & 0 &  & ... & 0  \\\\
          -\\frac{u_1}{\\cos(p_1)} & 1+\\frac{u_1 + u_2}{\\cos(p_1)} & -\\frac{u_2}{\\cos(b_1)} & 0 & ... & 0 \\\\
          0 & -\\frac{u_2}{\\cos(p_2)}  & 1+\\frac{u_2 + u_3}{\\cos(p_2)} & -\\frac{u_3}{\\cos(p_2)} &... & 0  \\\\
            &  & \\ddots & \\ddots & \\ddots & \\\\
          0 & 0 & ... & -\\frac{u_{n-2}}{\\cos(p_{n-2})}  & 1+\\frac{u_{n-2} + u_{n-1}}{\\cos(p_{n-2})} & -\\frac{u_{n-1}}{\\cos(p_{n-2})} \\\\
          0 & 0 & ... & 0 & -\\frac{u_{n-1}}{\\cos(p_{n-1})}  & 1+\\frac{u_{n-1}}{\\cos(p_{n-1})} \\\\
          \\end{array} \\right]
  
      with the substitue of:
  
      .. math::
  
          u_i = \\cos(b_i) K_i
  
      """
      phi_stag = np.deg2rad(lataxis.bounds)
      phi = np.deg2rad(lataxis.points)
      weight1 = np.cos(phi_stag)
      weight2 = np.cos(phi)
      diag = _make_diffusion_matrix(K, weight1, weight2)
      return diag, weight1, weight2