Chernoff-Hoeffding bounds for the number of nonzeros in a submatrix
Okay, here is a full answer.

We will use the fact, that any bipartite graph of maximum degree $d$ can be broken into (at most) $d$ matchings.

In our case, this means that we can split $A$ into (at most) $\sqrt{k}$ disjoint sets of elements $S_i\subseteq\{(i,j)\in[n]\times[n] \mid A_{i,j}=1\}$ of size at most $n$ such that any every element in each set has a unique row and column.
(This uses that a $1$ in $A$ can be seen as an edge in the $n\times n$ bipartite graph that is the matrix.)

Now let $X$ be the total number of $1$s in your $n/\sqrt k\times n/\sqrt k$ matrix, $A^\ell$, and let $X_i$ be the number of elements coming from $S_i$.
Then $X\ge t\implies \exists_i\, X_i\ge t/\sqrt k$ and so
$$\Pr[X\ge t]\le \sqrt k \Pr[X_1 \ge t/\sqrt k]$$
by the union bound.
Now, because elements in $S_1$ don't share any rows/columns, it is easy to analyze using your original approach.

In particular, reorder the rows and columns using the matching, so that the first element in $S_1$ is $(1,1)$ the next is $(2,2)$ and so on.
Let $r_i$ be the random variable that is 1 if row $i$ is chosen and $c_i$ likewise. Then we want to upper bound $\Pr[\sum_i r_ic_i \ge t]$.
Notice $E[r_ic_ir_jc_j] = E[r_ir_j]E[c_ic_j] \le E[r_i]^2E[c_i]^2$ by Maclaurin’s Inequality, so
$$\exp(t\sum_i r_ic_i) = \sum_k t^k/k!\,E\left(\sum_i r_ic_i\right)^k \le \sum_k t^k/k!\,E\left(\sum_i b_i\right)^k = \exp(t\sum_i b_i)$$
where $b_i$ is a normal binomial variable with $p=(n/\sqrt{k}/n)^2=1/k$.

Finally, we then get

$$\Pr[X\ge t]\le \sqrt k \Pr[X_1 \ge t/\sqrt k] \le \sqrt k \exp\left(-2\left(\frac{t-p n \sqrt{k}}{n\sqrt k}\right)^2 n\right)$$

or more simply

$$\Pr[X\ge \lambda\sqrt{nk}+n/\sqrt k]\le \sqrt k \exp(-2\lambda^2)$$

or if we use the Hoeffding bound for small values of $p$ (large $k$):

$$\Pr[X\ge \lambda\sqrt{nk}+n/\sqrt k]\le \sqrt k \exp\left(\frac{-\lambda^2}{2/k+\lambda/\sqrt n}\right)$$

Note that if we had simply assumed all elements were independent, we would have gotten the bound $\Pr[X\ge \lambda\sqrt{n\sqrt{k}}+n/\sqrt k]\le \exp(-2\lambda^2)$. Equal to ours except for a factor $k^{1/4}$ in the exponent.