A few observations and conjectures. First clearly $f$ is order $1$ and type $1$ too since $f(r)<e^r, r >0$ but $f(r)>>e^{(1-\epsilon)r}, r>0$, the last following from the fact that $$\lim_{n \to \infty} c^n\exp(\ln^2(n+1))=0$$ for all $0<c<1$ since $$-n |\log c|+\log^2(n+1) \to -\infty$$
In particular, $f$ has infinitely many zeroes $z_k$ which when ordered by increasing modulus form a sequence of genus $1$, so $\sum 1/|z_k|=\infty, \sum 1/|z_k|^{1+\epsilon}<\infty$ and $f(z)=e^{az+b}\Pi(1-\frac{z}{z_k}), -1 < a <1, b \in \mathbb R$ where the product is taken in the Cauchy sense as is not absolutely convergent.
Also, $f$ has a diagram $D$ (a convex compact set) and an indicator function $h(\theta)=\limsup_{r \to \infty}\frac{\log |f(re^{i\theta})|}{r}, -\pi \le \theta \le \pi$ which is its support function ($h(\theta)$ is the oriented projection distance of $D$ on the ray $\theta$, meaning the oriented distance from the origin to the farthest projection point of $D$ on the ray).
Note that $|h(\theta)| \le 1, h(0)=1, h(\pi)>-1$ by general theory since only $e^z$ can have $h(\pi)=-1$ when the function is type $1$ and $h(0)=1$ like here. We will show that $h(\pm \pi/2)=0$ which implies that the indicator is a segment $[b,1], 0 \le b <1$ so $$h(\theta)=\frac{1-b}{2}|\cos \theta|+\frac{1+b}{2}\cos \theta, h(\pi)=-b$$ Hence $f$ will decrease exponentially to infinity on the negative axis if $b>0$ which I kind of suspect but cannot prove.
Similarly, I suspect but cannot prove that $f$ is bounded on the imaginary axis which immediately gives the above results about $h$ and also implies that most zeroes of $f$ (in a definite sense of most) lie close to the imaginary line, so, in particular, the product representation above becomes absolutely convergent when we group the conjugate zeroes together (clearly $z_k$ root implies $\bar z_k$ root), but again this is conjecture only and is not implied by the fact that $\limsup_{r \to \infty}\frac{\log |f(\pm ir)|}{r}=0$ which I can prove as below.
Consider the (modified as the usual one is given as a function in $1/w$) Borel transform $$F(w)= \sum_{n= 0}^{\infty} \frac{w^n}{\exp(\ln^2(n+1))}, |w|<1$$ Polya's Theorem tells us that the convex hull of the singularities of $F$ inverted $w \to 1/w$ is the reflection of $D$ in the real axis (also called the conjugate diagram for obvious reasons and indeed it would be the convex hull of the singularities for the regular Borel transform).
Note that $F$ is singular at $1$ by Pringsheim theorem (coefficients positive and radius convergence $1$) but it must have other singularities too since only exponentials have point diagrams.
The Theorem of Le Roy and Lindelof gives analytic continuation for $G(w)=\sum g(n)w^n$ where $g(z)$ is an analytic function defined on some right half plane which grows slowly enough as $|z| \to \infty$. Here we can take $g(z)=\exp (-\log^2 (z+1)), \Re z \ge 0$, which actually decays as $\Re (-\log^2(z+1))=-\log^2|z+1|+\arg^2(z+1) \to -\infty, |z| \to \infty$ in the half plane, so it follows that $F(w)$ can be analytically continued except possibly outside the half line $[1, \infty]$. Note that $\infty$ may or may not be a singularity of $F$ as a multivalued function (in other words taking all continuations on paths from the original defining power series) and I suspect that it is not but cannot prove it. In particular if $[1,1/b], 0 \le b <1$ is the convex hull of singularities of $F$ in the above sense, the conjugate diagram of $f$ is $[b,1]$, but the reflection in the real axis of it is itself so the diagram of $f$ is $[b,1]$ which immediately implies the form of $h$ above by elementary geometry.
In conclusion we can say that $\limsup_{r \to \infty}\frac{\log |f(\pm ir)|}{r}=0$ and $\limsup_{r \to \infty}\frac{|\log f(-r)|}{r}=-b, 0 \le b <1$, but sadly not more as of now. If the things I suspect, so the boundness of $f$ on the imaginary axis and the analyticity of the modified Borel transform at infinity, are true then we could pretty much answer most of the OP questions as it would follow that $f$ goes (negative) exponentially to $0$ at $-\infty$ and most of its zeroes lie close to the imaginary axis, only negligible numbers being away from it.
Edit later:
As an aside, I would mention that the difficulties in directly proving anything about the behavior of $f$ on negative or purely imaginary numbers lies in the behavior of the Taylor series for $e^{-x}, x>0$ and $\cos x, \sin x, x\in \mathbb R$ whose partial sums actually increase very fast in absolute value up to a point and then they decrease. So for example if we take instead of $\sum z^n/n!$ the function $\sum z^n/(2n)!$ such a function cannot be bounded except on at most a half line (this one is actually highly unbounded on any half line but we can fiddle with signs and make a similar one bounded on a half line), and if we take something like $\sum z^n/(3n)!$ then its minimum modulus actually grows exponentially on a subsequence on the negative reals or on any half line like say the imaginary ones, so such an $f$ will be highly unbounded on any half line. Now in our case, the coefficients are just a little less than $1/n!$ exponentially speaking, since they are (eventually) greater than any $c^n/n!, 0<c<1$, so it is highly plausible that the behavior of $e^z$ repeats itself on both the imaginary line and on the negative reals, though as we saw $f$ cannot decay faster than $e^{-br}$ for $-r \to -\infty$ and some $0 \le b<1$, so already there is an imbalance between the growth at infinity which is close to $e^r$ and the (possible) decay at $-\infty$ and partial summation arguments are not fine enough (or at least I can not make them so) to prove even boundness there
