Does iterating the complex function $z\mapsto\frac{2\sqrt z}{1+z}$ always converge?
This is a very deep and interesting problem which was sort of completely solved by Gauss. The following is heavily borrowed from the paper The Arithmetic-geometric Mean of Gauss by David A. Cox which appeared in L'Enseignement Mathématique, Vol 30, 1984, pages 275-330.
Gauss considers the more general problem of agm of two complex numbers. Let us then assume that $a, b\in\mathbb {C} $ such that $ab\neq 0$ and $a\neq \pm b$ and let us define the AGM recurrence $$a_0=a,b_0=b,a_{n+1}=\frac{a_n+b_n}{2},b_{n+1}=(a_nb_n)^{1/2}\tag{1}$$ We have to fix the ambiguity involved in choosing square root now. Let us then say that a square root $b_1$ is the right choice for square root of $ab$ if $$|a_1-b_1|\leq|a_1+b_1|$$ and in case of equality $b_1/a_1$ must have positive imaginary part.
A pair of sequences $\{a_n\}, \{b_n\} $ defined by recurrences in $(1)$ is called good if $b_{n+1}$ is the right choice for $(a_nb_n) ^{1/2}$ for all but finitely many $n\geq 0$.
Cox mentions the following result in his paper
Theorem 1: If $a, b$ are complex numbers with $ab\neq 0,a\neq\pm b$ and $\{a_n\}, \{b_n\} $ are sequences defined by $(1)$ then they both converge to the same value. This common limit of both sequences is non-zero if and only if the pair of sequences is good.
This solves your problem that $z_n=b_n/a_n\to 1$ if the right branch of square root is chosen every time except for finitely many values of $n$.
But there is a lot more to come. Since we can make right or wrong choices of square root at each iteration, the limit of these sequences will depend on these choices. A complex number $\mu$ is called a value of AGM of $a, b$ and written $\mu=M(a, b) $ if there exist a good pair of sequences $\{a_n\}, \{b_n\} $ defined by $(1)$ and having a common limit $\mu$.
Thus based on allowed finite number of wrong choices of square roots there are a countable number of values of $M(a, b) $. Out of these there is a special one called the simplest value which is based on making right choice for square root in every iteration.
Gauss did some investigation to characterize all the values of $M(a, b) $ and Cox gives the corresponding result as
Theorem 2: Let $a, b$ be complex numbers with $ab\neq 0,a\neq\pm b$ and $|a|\geq |b|$. And further let $\mu, \lambda$ denote the simplest values of $M(a, b), M(a+b, a-b) $ respectively. Then all values $\mu'$ of $M(a, b) $ are given by $$\frac{1}{\mu'}=\frac {d} {\mu} +\frac{ic} {\lambda} $$ where $c, d$ are any arbitrary integers comprime to each other and $c\equiv 0\pmod {4},d\equiv 1\pmod {4}$.
The proof involves all the ideas related to modular functions, modular group, fundamental region etc and it is an interesting read. Cox also says that Gauss knew a lot of this material and gives many historical details in his paper.
The answer to your question is yes. The iterations of $f$ always converge. Furthermore the convergence is uniform towards 1 on any compact subset excluding a neighborhood of the origin and $-1$, i.e. on a subset of the following form: ($0<r<<R<+\infty$) $$K=K_{r,R}=\{z\in {\Bbb C} : |z|\geq r, |z-1|\geq r, |z|\leq R\}.$$
Complex dynamics is a useful tool in this context. If you let $H=\{z: {\rm Re\ } z >0\}$, then $K'=f(K)$ is a compact subset of $H$ and $f(K')\subset K'$. The right half-plane $H$ is what is known as a hyperbolic domain. If you look at the wiki page for Schwartz_lemma, about halfway down the page, you will find that a half-plane admits a Poincare metric $d_H$ which is contracted by any holomorphic map of $H$ into itself. The metric (cf. wiki-page for Poincaré-metric) on the right half-plane takes the explicit form: $$ d(z,w) = 2\tanh^{-1} \frac{|z-w|}{|z+w|}.$$
If $f$ is not an automorphism of $H$ (not our case) then $f$ is a strict contraction on compact subsets, in particular the set $K'$ above. More precisely, there is $\theta=\theta({r,R})<1$ so that for $z_1,z_2\in K'$ one has $$ \frac{|f(z_1)-f(z_2)|}{|f(z_1)+f(z_2)|} \leq \theta \frac{|z_1-z_2|}{|z_1+z_2|}.$$ Taking $z_2=1=f(z_2)$ and iterating the inequality, you see that the distance between $f^n(z_1)$ and $1$ goes exponentially fast to zero (in fact super-exponentially fast, since $f'(1)=0$). The claim follows by letting $r\rightarrow 0$ and $R\rightarrow +\infty$.
This function at $z$ is the same as at $1/z$, excluding the case $\mathbb R\ni z<0$, which is the only case where $\sqrt{1/z}\neq1/\sqrt z$.
$$\frac{2\sqrt{1/z}}{1+1/z}=\frac{2\cdot1/\sqrt z}{1/z+1}\cdot\frac zz=\frac{2\sqrt z}{1+z}.$$
Since $|\theta_2|\leq\pi/4$, if also $|z_2|\leq1$, then the circular sector represented by these inequalities is contained in the unit disk centred at $1$, so the sequence converges. (See the last part of the Question.) On the other hand, if $|z_2|>1$, then we can replace $z_2$ with $1/z_2$ without affecting the later terms, so again the sequence converges.
That answers the Question as asked.
Now I'll also answer a related question: Is there a well-defined arithmetic-geometric mean of two complex numbers?
Given $a_0,b_0\in\mathbb C$, define the sequences
$$a_{n+1}=\text{AM}(a_n,b_n)=\frac{a_n+b_n}{2},\quad b_{n+1}=\text{GM}(a_n,b_n)=\begin{cases}0,&a_n=0\\a_n\sqrt{\frac{b_n}{a_n}},&a_n\neq0\end{cases}$$
(again using the principal square root). This particular form of GM, rather than $\sqrt{ab}$, ensures that the result is "between" $a$ and $b$; if $a$ and $b$ are linearly independent over $\mathbb R$, then GM$(a,b)$ is a conical combination of them. (That's a linear combination with non-negative coefficients.) It behaves nicely with rotations: $\text{GM}(ka,kb)=k\,\text{GM}(a,b)$ for any $k\in\mathbb C$.
In any of the cases where $a_0b_0(a_0\!^2-b_0\!^2)=0$, it's easy to show that both sequences converge, to $0$ or to $a_0=b_0$. So let's assume that $a_0b_0(a_0\!^2-b_0\!^2)\neq0$.
Define $z_n=b_n/a_n$; then the formula for $z_{n+1}$ is exactly that in the OP, and we've already shown that $\lim_{n\to\infty}z_n=1$.
From the triangle inequality,
$$|a_{n+1}|\leq\frac{|a_n|+|b_n|}{2}\leq\frac{\max(|a_n|,|b_n|)+\max(|a_n|,|b_n|)}{2}=\max(|a_n|,|b_n|),$$
and similarly $|b_{n+1}|=\sqrt{|a_n||b_n|}\leq\max(|a_n|,|b_n|)$, so
$$\max(|a_{n+1}|,|b_{n+1}|)\leq\max(|a_n|,|b_n|)\leq\max(|a_{n-1}|,|b_{n-1}|)\leq\cdots\leq\max(|a_0|,|b_0|).$$
Thus, we see that both sequences are bounded, so
$$|a_n-b_n|=|a_n||1-z_n|\leq\max(|a_0|,|b_0|)\,|1-z_n|\to0;$$
if either one converges, they must both converge to the same value.
Indeed one of them does converge: $|a_{n+1}-a_n|=\tfrac12|a_n-b_n|\to0$, and for $k>1$,
$$|a_{n+k}-a_n|=\left|\sum_{j=0}^{k-1}(a_{n+j+1}-a_{n+j})\right|$$
$$\leq\sum_{j=0}^{k-1}|a_{n+j+1}-a_{n+j}|$$
$$=\frac12\sum_{j=0}^{k-1}|a_{n+j}-b_{n+j}|$$
$$\leq\frac12\max(|a_0|,|b_0|)\sum_{j=0}^{k-1}|1-z_{n+j}|$$
$$\leq\frac12\max(|a_0|,|b_0|)\sum_{j=0}^\infty|1-z_{n+j}|.$$
From the last part of the Question and the first part of this Answer, for $n\geq3$ we have $|1-z_n|<1$ and $|1-z_{n+j}|\leq|1-z_n|^{2^j}$, so
$$\sum_{j=0}^\infty|1-z_{n+j}|\leq\sum_{j=0}^\infty|1-z_n|^{2^j}$$
$$\leq\sum_{l=1}^\infty|1-z_n|^l$$
$$=\frac{|1-z_n|}{1-|1-z_n|}\to0.$$
Therefore, by Cauchy's criterion, $a_n$ converges. We may call the common limit
$$\lim_{n\to\infty}a_n=\lim_{n\to\infty}b_n=\text{AGM}(a_0,b_0).$$