Let me summarize and condense what @IosefPinelis has said.
We are to estimate $$f(y,L) = \int_L^\infty e^{-t-y/t} dt = \sqrt{y} g(L/\sqrt{y},\sqrt{y}),$$
where $g(r,x)=\int_r^\infty e^{-x(s+1/s)} ds.$ By a substitution of variables $t = r+1/r$,
$$g(r,x) = \frac{1}{2} I(r,x) + \begin{cases} 0 &\text{if $r\geq 1$,}\\ J(r,x)& \text{if $0<r<1$,}\end{cases}$$
where $$I(r,x) = \int_{r+1/r}^\infty e^{-x t} (F(t)+1) dt,\;\;\;\;\;\;\;
J(r,x) = \int_2^{r+1/r} e^{-x t} F(t) dt$$ and $F(t) = \frac{t}{\sqrt{t^2-4}}$. (We are using the fact that $F(t)$ is an odd function.) Letting $t = u+2$, we see that
$$I(r,x) = e^{-2 x} \int_{u_r}^\infty e^{-x u} (G(u)+1) du,
\;\;\;\;J(r,x)= e^{-2 x}\int_0^{u_r} e^{-x u} G(u) du,$$
where $G(u) = F(u+2) = \frac{1}{\sqrt{u}} \frac{1+u/2}{\sqrt{1+u/4}}$ for $u\geq 0$ and $u_r = r+1/r-2$.
We can give different upper bounds on $G(u)$. @IosefPinelis used the bound $H(u)\leq H(u_r) + (u-u_r) H'(u_r)$ for $H(u) = \sqrt{u} G(u)$ and $u\geq u_r$ (by concavity of $H(u)$) and, in particular, $H(u)\leq 1 + \frac{3 u}{8}$ for $u\geq 0$; he showed how to obtain bounds on $g(r,x)$ in consequence. He also mentioned once a bound that is useful, if crude, viz., given by $G(u)\leq 1/\sqrt{u} + 1$ for $u\geq 0$. (It can overestimate $G(u)$ by up to almost 50%, and $1+G(u)$ by up to almost 30%.) Let us work with this crude bound in detail.
We obtain
$$\begin{aligned}I(r,x)&\leq e^{-2 x}\int_{u_r}^\infty e^{-x u} \left(2 + \frac{1}{\sqrt{u}}\right) du = \frac{e^{-x (u_r+2)}}{x} + \frac{e^{-2 x}}{\sqrt{x}} \int_{x u_r}^\infty e^{-v} \frac{dv}{\sqrt{v}} \\ &= \frac{e^{-x (r+1/r)}}{x} + \frac{2 e^{-2 x}}{\sqrt{x}} \int_{\sqrt{xu_r}}^\infty e^{-t^2} dt
= \left(\frac{1}{x} + \frac{2}{\sqrt{x}} M(\sqrt{x u_r}) \right)e^{-x (r+1/r)},\end{aligned}$$
where $M(x)$ is Mills' ratio $M(x) = e^{x^2} \int_x^\infty e^{-t^2} dt$.
The special case $u_r=0$ in that, according to an answer by @openletter.mousetail.nl, it leads to a fancy special function called a leaky aquifer function (too fancy for the NIST handbook, apparently). Of course we can just use the bound above.
Just as above,
$$\begin{aligned}J(r,x) &= e^{-2 x} \int_0^{u_r} e^{-x u} \left(1 + \frac{1}{\sqrt{u}}\right) du = e^{-2 x} \frac{1-e^{-x u_r}}{x} + \frac{e^{-2 x}}{\sqrt{x}} \int_0^{x u_r} e^{-v} \frac{dv}{\sqrt{v}} \\ &= \frac{e^{-2 x}-e^{-x (r+1/r)}}{x} + \frac{2 e^{-2 x}}{\sqrt{x}} \int_0^{\sqrt{xu_r}} e^{-t^2} dt
\leq \frac{e^{-2 x}}{x} + \sqrt{\frac{\pi}{x}} e^{-2 x} -
\left(\frac{1}{x} + \frac{2}{\sqrt{x}} M(\sqrt{x u_r}) \right) e^{-x (r+1/r)}.\end{aligned}$$
We conclude that
$$g(r,x) = \begin{cases}\left(\frac{1}{2 x} + \frac{M(\sqrt{x u_r})}{\sqrt{x}}\right) e^{-x (r+1/r)} &\text{if $r\geq 1$,}\\ \left(\frac{1}{x} + \sqrt{\frac{\pi}{x}}\right) e^{-2 x} - \left(\frac{1}{2 x} + \frac{M(\sqrt{x u_r})}{\sqrt{x}}\right) e^{-x (r+1/r)}& \text{if $0<r<1$,}
\end{cases}$$
For $x\geq 0$, $\frac{1}{x+\sqrt{x^2+2}} < M(x)\leq \frac{1}{x+\sqrt{x^2+(4/\pi)}}$ and $\frac{\sqrt{\pi}}{2 \sqrt{\pi} x + 2} \leq M(x)<\frac{1}{x+1}$ (NIST Handbook, section 7.8).
Of course we also have the upper bound $$g(r,x)\leq g(0,x) = \int_{0}^{\infty} e^{-x (s+1/s)} ds = 2 K_1(2 x),$$ where $K_1$ is the modified Bessel function of the second kind. The leading-order terms of $2 K_1(2 x)$ are $\sqrt{\frac{\pi}{x}} e^{-2 x} \left(1 + \frac{3}{8 x} + \dotsc\right)$. As @igorKhavkine points out, this is very good for $r$ close to $0$ (or $r\in (0,1)$ and $x$ large).
I'll add some remarksLet's talk about an application, in part so that we can see what happens for some specific range of values of the parameters.
Assume we have a bitRiemann-Hypothesis check up to height $H$, i.e., assume we know that $\Re s = 1/2$ for every non-trivial zero $s$ of $\zeta(s)$ with $|\Im s|\leq H$. (This is currently known rigorously for $H \leq 3\cdot 10^{12}$ (Platt-Trudgian, 2021).) Assume as well that we have a zero-free region of classical form: for every non-trivial zero $s$ of $\zeta(s)$, $\Re s \leq 1 - \frac{1}{C \Im s}$. (This is known with $C = 5.573\dotsc$ (Mossinghoff-Trudgian, 2015); there's a very recent improvement with $C = 5.558\dotsc$ (Mossinghoff-Trudgian-Yang).)
The problem is to estimate sums of the form
$$\sum_{n\leq X} \Lambda(n) \log \frac{X}{n},$$
where $\Lambda(n)$ is the von Mangoldt function. (Here $\log \frac{X}{n}$ is a very natural smoothing that often pops up of its own accord.) By standard Mellin-transform magic and the assumptions above, this problem reduces to estimating $$I = \int_H^\infty X^{1 - \frac{1}{C \log \tau}} \frac{dt}{t^2}.$$
(I am not forgetting a factor of $\log t$ here; there are good explicit density results that enable us to do this (Kadiri 2013, in the style of Bohr-Landau) and even more (Kadiri-Lumley-Ng 2018).)
By the substitution $\tau = e^t$,
$$I = \int_{\log H}^\infty X^{1 - \frac{1}{C t}} \frac{dt}{e^t}
= X \int_{\log H}^\infty e^{- t- \frac{(\log X)/C}{t}} dt.$$
In other words, $I = X f(y,L)$ with $y = \frac{\log X}{C}$, $L = \log H$, and so $$I = x g\left(\frac{\log H}{x},x\right)\cdot X$$
with $x = \sqrt{\frac{\log X}{C}}$.
We can assume $X\geq 10^{19}$ or so (as estimates for sums with smaller $X$ can be obtained by brute force). Hence, we can assume
$$x\geq \sqrt{\frac{\log 10^{19}}{C}} = 2.8018\dotsc,$$ $$r \leq \frac{\log H}{x} \leq 10.254.$$
