The boundaries of 2+1D fermionic topological orders
$2+1$D bosonic topological orders can be characterized by the $S,T$ matrices that encode the statistics of topological excitations. In particular, the $S,T$ matrices can be used to systematically obtain the gapped boundaries of bosonic topological orders. Such an approach, however, does not naively apply to fermionic topological orders (FTOs). In this work, we propose a systematic approach to obtain the gapped boundaries of $2+1$D abelian FTOs. The main trick is to construct a bosonic extension in which the fermionic excitation is "condensed" to form the associated FTOs. Here we choose the parent bosonic topological order to be the $\mathbb{Z}_2$ topological order, which indeed has a fermionic excitation. Such a construction allows us to find an explicit correspondence between abelian FTOs (described by odd $K$-matrix $K_F$) and the "fermion-" condensed $\mathbb{Z}_2$ topological orders (described by even $K$-matrix $K_B$). This provides a systematic algorithm to obtain the modular covariant boundary partition functions as well as the boundary topological excitations of abelian FTOs. For example, the $ν=1-\frac{1}{m}$ Laughlin's states have exactly one type of gapped boundary when $m$ is a square, whose boundary excitations form a $\mathbb{Z}_{2}\times\mathbb{Z}_{\sqrt{m}}$ fusion ring. Our approach can be easily generalized to obtain gapped and gapless boundaries of non-abelian fermionic topological orders.