Ok, you can, after some calculations see that
the absolute value of your sum is less than or equal to
$$\sum_{n=-\infty}^\infty e^{\pi(-tn^2 + 2n|z|)}$$
where $t = im(\tau) >0.$
The part of the sum where $n<0$ will be less than some finite constant
depending only on $\tau,$ so we only need to worry about $n>0.$
(This you need to prove.)
Ok, so, divide that sum into two parts:
Let $N$ be the least integer greater than $4|z|/t$
(where $t = im(\tau) >0.$)
If $n\geq N$ we have $nt/4 > |z|.$
This will imply that $-tn^2 + 2n|z| \leq -n^2t/2.$
Use this on the second sum, $N...\infty$ and use some other estimations on the first part.
Edit: Continuation:
First, you should see that
$$\sum_{n=-\infty}^0 e^{\pi(-tn^2 + 2n|z|)} \leq
\sum_{n=-\infty}^0 e^{-\pi t n^2}
$$
and this clearly converges to something that only depends on $\tau.$
(The exponent will diverge to $-\infty$ rapidly, so the sum must converge.)
Now, if $n\geq N$ then
$$\sum_{n=N}^\infty e^{\pi(-tn^2 + 2n|z|)}\leq \sum_{n=N}^\infty e^{\pi(-n^2 t/2)}
\leq \sum_{n=0}^\infty e^{\pi(-n^2 t/2)}
$$
which converges to a finite number independent of $z$.
Now, you only need to do something about the last part:
$$\sum_{n=0}^{N-1} e^{\pi(-tn^2 + 2n|z|)}$$
(Hint: use $N-1< 4|z|/t$ which follows from the conditions above).
Edit: Continuation:
Ok, choose $r \in \mathbb{R}$ so that $\theta(r,\tau) \neq 0.$
By quasi-periodicity, we have for integers $m$ that
$$\theta(r+m\tau,\tau) = e^{-i\pi m^2\tau - 2im \pi r}\theta(r,\tau)$$
so by taking absolute values on each side, we see that
$$|\theta(r+m\tau,\tau)| = e^{t\pi m^2}|\theta(r,\tau)|$$
where $t = im(\tau)>0.$
Hence,
$$\lim_{m \rightarrow \infty} \frac{|\theta(r+m\tau,\tau)|}{e^{t\pi m^2}}
= |\theta(r,\tau)| > 0.$$
Thus, the order is at least 2.
They arise from an effort to find a product expansion of the elliptic functions by using the Weierstrass theorem to factor the zeros and poles.
I hope the following outline gives an intuitive motivation for the creation of the theta functions.
Naively we could write
$$\text{sn}(z)=z\pi \prod \frac{1-\frac{z}{2Kn+2K^{\prime}mi} }{ 1-\frac{z}{2Kn+2K^{\prime}mi +K^{\prime}i}}.$$
Define $\tau=\frac{K^{\prime}i}{K}$ then
$$\text{sn}(2Kz)=2Kz\prod \left( \frac{1-\frac{z}{n+m\tau} }{1-\frac{z}{n+m\tau+\frac{1}{2}\tau}}\right).$$
However neither product
$\prod\limits_{nm} 1-\frac{z}{n+m\tau}$ or $\prod\limits_{nm} 1-\frac{z}{n+m\tau+\frac{1}{2}\tau}$ converges absolutely. There are two ways to make convergent products out of these expresssions. One is to introduce convergence factors following Weierstrauss, and this leads to the $\sigma$ functions. Another way is to take the product first by $n$ and then by $m$. This results in the important theta functions.
We define
$$\Theta_1(z)=2Kz \prod\limits_m \prod\limits_n \left( 1- \frac{z}{n+m\tau} \right).$$
After some manipulation and with $q=e^{\tau \pi i}$,
$$\Theta_1(z)=\frac{2K}{\pi}\sin(\pi z)\frac{\prod\limits_{m=1}^{\infty}1-2q^{2m}\cos(2z\pi)+q^{4m}}{\prod\limits_{m=1}^{\infty}\left[ 1-q^{2m} \right]^2}$$
is well defined with absolutely convergent numerator and denomenator.
Similarly we define
$$\Theta_0(z)= \prod\limits_m \prod\limits_n \left( 1- \frac{z}{n+(m+\frac{1}{2})\tau} \right)$$ and
obtain,
$$\Theta_0(z)=\frac{\prod\limits_{m=1}^{\infty}1-2q^{2m-1}\cos(2z\pi)+q^{4m-2}}{\prod\limits_{m=1}^{\infty}\left[ 1-q^{2m-1} \right]^2}.$$
As a result we have a product expression for the function $\text{sn}(2Kz)$.
$$\text{sn}(2Kz)=\frac{\Theta_1(z)}{\Theta_0(z)}.$$
Similarly we may define the functions,
$$\Theta_2(z)=\cos(\pi z)\frac{\prod\limits_{m=1}^{\infty}1+2q^{2m}\cos(2z\pi)+q^{4m}}{\prod\limits_{m=1}^{\infty}\left[ 1+q^{2m} \right]^2}$$
and
$$\Theta_3(z)=\frac{\prod\limits_{m=1}^{\infty}1+2q^{2m-1}\cos(2z\pi)+q^{4m-2}}{\prod\limits_{m=1}^{\infty}\left[ 1+q^{2m-1} \right]^2}.$$
And we have
$$\text{cn}(2Kz)=\frac{\Theta_2(z)}{\Theta_0(z)}$$
$$\text{dn}(2Kz)=\frac{\Theta_3(z)}{\Theta_0(z)}.$$
Best Answer
UPD: the previous version contained a square which shouldn't be there.
Actually, your function is even more simply expressed in terms of $\vartheta_4$-function. Also, I prefer this notation in which $$f(y)=\vartheta_4(0,e^{-y})=\vartheta_4\Bigl(0\Bigr|\Bigl.\frac{iy}{\pi}\Bigr).$$ I.e. I use the convention $\vartheta_k(z,q)=\vartheta_k(z|\tau)$.
Then, to obtain the asymptotics as $y\rightarrow 0^+$, we need two things:
Jacobi's imaginary transformation, after which the transformed nome and half-period behave as $q'\rightarrow0$, $\tau'\rightarrow i\infty$ (instead of $q\rightarrow1$, $\tau\rightarrow0$): $$\vartheta_4\Bigl(0\Bigr|\Bigl.\frac{iy}{\pi}\Bigr)=\sqrt{\frac{\pi}{y}}\vartheta_2\Bigl(0\Bigr|\Bigl.\frac{i\pi}{y}\Bigr).$$
Series representations for theta functions (e.g. the formula (8) by the first link), which implies that $$\vartheta_2(0,q')\sim 2(q')^{\frac14}$$ as $q'\rightarrow 0$. Note that you can also obtain an arbitrary number of terms in the asymptotic expansion if you want.
Taking into account the two things above, we obtain that the leading asymptotic term is given by $$f(y\rightarrow0)\sim 2\sqrt{\frac{\pi}{y}} \exp\left\{-\frac{\pi^2}{4y}\right\}.$$