Gauss's Law does not apply if there is a point charge on the Gaussian surface. That is the case even if we try to adopt a convention to treat the charge as being 'inside' or 'outside' the surface - either way the law will not apply.
To see this we can calculate the total flux $\Phi$ such a charge produces through a sphere $S$. This flux is not directly well-defined since the electric field at the charge itself is infinite, however $\Phi$ does have a well-defined value if we calculate it as a limit. For any $\delta \in (0, \pi/2)$ we can calculate the total flux $\Phi_{\delta}$ through an approximation $S_{\delta}$ to $S$ which is defined by the spherical coordinates (using the Physics convention [SPH]) : $\theta \in [\delta, \pi - \delta]$ and $\varphi \in [0, 2\pi]$, since the electric field does not become infinite at any point of $S_{\delta}$. Then we can define $\Phi = \lim_{\delta \rightarrow 0^+} \Phi_{\delta}$.
There is a second reason for using this limiting process namely to be able to apply the spherical formula for surface integration [MSE] which requires the surface not to intersect the $z$-axis so that a division by $0$ doesn't occur in that formula when $\sin \theta = 0$.
Let $S$ have radius $R$ and set up $xyz$-axes so that the point charge $q$ lies on the positive $z$-axis at a distance $R$ from $O$. We will calculate the total flux $\Phi$ through $S_{\delta}$ due to $q$.
Let $\underline{\mathbf{\hat{r}}} = \sin \theta \cos \varphi \: \underline{\mathbf{i}} + \sin \theta \sin \varphi \: \underline{\mathbf{j}} + \cos \theta \: \underline{\mathbf{k}}$ be the radial unit vector, so the position vector $\overrightarrow{QP}$ of the general point $P$ on $S_{\delta}$ wrt $q$ is given by $\overrightarrow{QP} = -R \, \underline{\mathbf{k}} + R \, \underline{\mathbf{\hat{r}}}$. Then at $P$ we have :
$$
\underline{\mathbf{E}} = \frac{q}{4 \pi \epsilon_0} \cdot \frac{-R \, \underline{\mathbf{k}} + R \, \underline{\mathbf{\hat{r}}}}{| {-R} \, \underline{\mathbf{k}} + R \, \underline{\mathbf{\hat{r}}} \: |^3}
$$
where
$$
|\underline{\mathbf{\hat{r}}} - \underline{\mathbf{k}}|^2 = (1 - \cos \theta)^2 + \sin^2 \theta = 2 - 2 \cos \theta.
$$
Since a unit outward normal to $S_{\delta}$ at $P$ is $\underline{\mathbf{n}} = \underline{\mathbf{\hat{r}}}$ and since $\underline{\mathbf{k}} \cdot \underline{\mathbf{\hat{r}}} = \cos \theta$ we have :
$$
\underline{\mathbf{E}} \cdot \underline{\mathbf{n}} = \frac{q}{4 \pi \epsilon_0 R^2} \cdot \frac{1}{2 \sqrt{2}} \cdot \frac{1}{(1 - \cos \theta)^{1/2}}
$$
Using the formula for surface integration in spherical coordinates described in [MSE], with $r(\theta, \varphi) = R, \; \forall \; (\theta, \varphi) \in D_{\delta} = [\delta, \pi - \delta] \times [0, 2 \pi]$ and $\frac{\partial r}{\partial \theta} = \frac{\partial r}{\partial \varphi} = 0$ throughout $S_{\delta}$ we have :
\begin{eqnarray*}
\Phi_{\delta} = \int_{S_{\delta}} \underline{\mathbf{E}} \cdot \underline{\mathbf{n}} & = & \int_{D_{\delta}} \; \frac{q}{8 \pi \epsilon_0 R^2 \sqrt{2}} \cdot \frac{1}{(1 - \cos \theta)^{1/2}} \cdot R^2 \sin \theta \\
& = & \int_{\delta}^{\pi - \delta} \int_0^{2 \pi} \frac{q}{8 \pi \epsilon_0 \sqrt{2}} \cdot \frac{\sin \theta}{(1 - \cos \theta)^{1/2}} \; d\varphi \; d\theta \\
& = & \int_{\delta}^{\pi - \delta} \frac{q}{4 \sqrt{2} \epsilon_0} \cdot \frac{\sin \theta}{(1 - \cos \theta)^{1/2}} \; d\theta = \frac{q}{4 \sqrt{2} \epsilon_0} \cdot \left[ 2(1 - \cos \theta)^{1/2} \right]_{\delta}^{\pi - \delta} \\
\mbox{and hence} \hspace{2em} \Phi & = & \lim_{\delta \rightarrow 0^+} \Phi_{\delta} = \frac{q}{2 \sqrt{2} \epsilon_0} ( \sqrt{2} - 0 ) = \frac{q}{2 \epsilon_0}
\end{eqnarray*}
We can calculate the above integral for $q$ at any distance $R_2 \geq 0$ from $O$ :
\begin{eqnarray*}
\underline{\mathbf{E}} & = & \frac{q}{4 \pi \epsilon_0} \cdot \frac{-R_2 \, \underline{\mathbf{k}} + R \, \underline{\mathbf{\hat{r}}}}{| {-R_2} \, \underline{\mathbf{k}} + R \, \underline{\mathbf{\hat{r}}} \: |^3} = \frac{q}{4 \pi \epsilon_0 R^2} \cdot \frac{-\frac{R_2}{R}\underline{\mathbf{k}} + \underline{\mathbf{\hat{r}}}}{|-\frac{R_2}{R}\underline{\mathbf{k}} + \underline{\mathbf{\hat{r}}}|^3} \\
\mbox{where} \hspace{2em} \left|-\frac{R_2}{R}\underline{\mathbf{k}} + \underline{\mathbf{\hat{r}}}\right|^2 & = & \left(\frac{R_2}{R} - \cos \theta\right)^2 + \sin^2 \theta = 1 + \frac{R_2^2}{R^2} - 2\frac{R_2}{R} \cos \theta.
\end{eqnarray*}
Then we obtain :
\begin{eqnarray*}
\Phi_{\delta} = \int_{S_{\delta}} \underline{\mathbf{E}} \cdot \underline{\mathbf{n}} & = & \frac{q}{4 \pi \epsilon_0} \int_{D_{\delta}} \frac{-(R_2/R)\sin \theta \cos \theta + \sin \theta}{\left( 1 + (R_2/R)^2 -2 (R_2/R) \cos \theta \right)^{3/2}} \\
& = & \frac{q}{2 \epsilon_0} \cdot \left[ \left(\cos \theta - \frac{R}{R_2}\right) \left(1 + \frac{R_2^2}{R^2} -2\frac{R_2}{R}\cos \theta\right)^{-1/2} + \frac{R}{R_2}\left(1 + \frac{R_2^2}{R^2} -2\frac{R_2}{R} \cos \theta \right)^{1/2} \right]_{\delta}^{\pi - \delta}, \\
& & \mbox{using integration by parts for the first part of the integral.}
\end{eqnarray*}
Then we obtain :
$$
\Phi = \lim_{\delta \rightarrow 0^+} \Phi_{\delta} = \left\{ \begin{array}{ll} q/\epsilon_0, & \mbox{if } 0 \leq R_2 < R \\[1ex] q/2\epsilon_0, & \mbox{if } R_2 = R \\[1ex] 0, & \mbox{if } R_2 > R \end{array} \right.
$$
in keeping with Gauss's Law in the case $R_2 \neq R$. The electric flux $\Phi$ is a discontinuous function of the distance $R_2$ of $q$ from the centre of the sphere. If we had multiple point charges on the surface of $S$, and/or a continuous surface charge distribution $\sigma$ across $S$ then because of additivity of the electric field the above formula for $\Phi$ would still apply but with $q$ replaced by $Q$, the total surface charge. If $S$ was exposed only to a volume charge distribution $\rho$ then the total surface charge would be $Q = 0$.
A similar process of taking a limit can be used to calculate the average electric potential $V_{\mathrm{avg}}$ over a sphere $S$ (eg see [DJG, p117-118]) due to a point charge $q$ on its surface. The potential will be infinite at the charge itself so $V_{\mathrm{avg}}$ is not well-defined, but as we did with the flux $\Phi$ above we can use a limit to give it a well-defined value :
\begin{eqnarray*}
V_{\mathrm{avg}} & = & \lim_{\delta \rightarrow 0^+} \, \frac{1}{A_{\delta}} \, \int_{S_{\delta}} V \; dS, \hspace{2em} \mbox{where $A_{\delta}$ is the area of $S_{\delta}$} \\
& = & \lim_{\delta \rightarrow 0^+} \, \frac{1}{A_{\delta}} \, \int_{\delta}^{\pi - \delta} \int_0^{2 \pi} \; \frac{q}{4 \pi \epsilon_0} \cdot \frac{R \sin \theta}{\left( 2 - 2 \cos \theta \right)^{1/2}} \; d\varphi \; d\theta = \frac{q}{4 \pi \epsilon_0 R}
\end{eqnarray*}
and for $q$ not on $S$ at distance $R_2 \geq 0$ from $O$ :
\begin{eqnarray*}
V_{\mathrm{avg}} & = & \lim_{\delta \rightarrow 0^+} \frac{1}{A_{\delta}} \int_{\delta}^{\pi - \delta} \int_0^{2 \pi} \; \frac{q}{4 \pi \epsilon_0} \cdot \frac{R^2 \sin \theta}{\left( R^2 + R_2^2 - 2R R_2 \cos \theta \right)^{1/2}} \; d\varphi \; d\theta = \frac{q}{8 \pi \epsilon_0 R R_2} \left[ (R + R_2) - |R - R_2| \right],
\end{eqnarray*}
where in the latter integral the only reason we have used the limiting process is so the spherical surface integral formula [MSE] is valid, ie. to avoid a division by $0$ in that formula. Then we have :
$$
V_{\mathrm{avg}} = \left\{ \begin{array}{ll} \displaystyle \frac{q}{4 \pi \epsilon_0 R}, & \mbox{if $0 \leq R_2 < R$, (ie. as if $q$ at centre) } \\[3ex] \displaystyle \frac{q}{4 \pi \epsilon_0 R}, & \mbox{if $R_2 = R$, (ie. as if $q$ at centre)} \\[3ex] \displaystyle \frac{q}{4 \pi \epsilon_0 R_2}, & \mbox{if $R_2 > R$, (ie. potential at centre due to $q$)} \end{array} \right.
$$
In contrast to the flux $\Phi$, $V_{\mathrm{avg}}$ is a continuous function of $R_2$.
Another situation where the limiting process applies is in calculating the potential at a point $P$ due to a uniformly charged spherical shell of total charge $Q$ where by performing a similar integration as those above we obtain :
$$
V = \left\{ \begin{array}{ll} \displaystyle \frac{Q}{4 \pi \epsilon_0 R}, & \mbox{if } r < R \\[3ex] \displaystyle \frac{Q}{4 \pi \epsilon_0 R}, & \mbox{if } r = R \\[3ex] \displaystyle \frac{Q}{4 \pi \epsilon_0 r}, & \mbox{if } r > R \end{array} \right.
$$
where $R$ is the radius of the shell and $P$ is at distance $r$ from its centre. In each of the 3 cases the limit process avoids the division by $0$ in the spherical surface integral formula [MSE]. On the shell's surface it secondly avoids the problem of an infinite potential due to a non-zero charge density at $P$, which also produces a division by $0$. From Gauss's Law and the calculation above, the total flux $\Phi$ of this charge distribution through a Gaussian surface of radius $r$ and concentric with the shell is given by :
$$
\Phi = \left\{ \begin{array}{ll} \displaystyle 0, & \mbox{if } r < R \\[1ex] \displaystyle Q/2\epsilon_0, & \mbox{if } r = R \\[1ex] \displaystyle Q/\epsilon_0, & \mbox{if } r > R \end{array} \right.
$$
For the case $r = R$ the flux can't be calculated from Gauss's Law, but it can be calculated from the Generalized Gauss Theorem discussed in the answer [PSE] to the present question. Because of the symmetry, the electric field of the uniformly charged spherical shell must be radial and a function of $r$ only. From the latter formula its magnitude $E$ is given by :
$$
E = \left\{ \begin{array}{ll} 0, & \mbox{if } r < R \\[2ex] \displaystyle \frac{Q}{8 \pi \epsilon_0 R^2}, & \mbox{if } r = R \\[2ex] \displaystyle \frac{Q}{4 \pi \epsilon_0 r^2}, & \mbox{if } r > R \end{array} \right.
$$
However the value of $E$ for $r = R$ is purely a value calculated within the mathematical model which arises from the physical model of electrostatics and this physical model reaches the limits of its applicability when considering electric field right on a surface of zero thickness - it would be unlikely to be possible to measure such a field value in an experiment or to detect it in any way.
In practice the limiting process used above is not normally written out, and the limits $0$ and $\pi$ are just applied to the integration and the correct answer is arrived at. A square root term in the denominator which goes to zero becomes a square root term in the numerator once integrated and then going to zero is not a problem. But a check should be made mentally that a problem is not arising due to a division by $0$.
References
[SPH] Spherical coordinate system, https://en.wikipedia.org/wiki/Spherical_coordinate_system
[MSE] Maths StackExchange Answer : Surface integrals in spherical coordinates, https://math.stackexchange.com/a/4867164/417024
[DJG] David J. Griffiths (2013), Introduction To Electrodynamics, 4th Edition, Pearson Education Inc.
[PSE] Physics StackExchange Answer : Using Gauss's law when point charges lie exactly on the Gaussian surface, https://physics.stackexchange.com/a/544481/111652
