Joint pdf of N generally correlated (absolute values of) R.Vs as a joint pdf of (absolute value squared) random variables The joint pdf of $|g_1|,|g_2|,\ldots,|g_N|$ is given by
$P_{|g_1|,|g_2|,\ldots\ldots,|g_N|}(r_1,r_2,\ldots,r_N)=$ $$\prod_{
\substack{k=1 \\ \mu_1 \triangleq 0}}^N \frac{2r_k}{\sigma^2(1-\mu_k^2)} \exp\left(-\frac{r_k^2 + \mu_k^2r_1^2}{\sigma^2(1-\mu_k^2)} \right) I_0 \left(\frac{2\mu_kr_1r_2}{\sigma^2(1-\mu_k^2)}\right) $$
where the $k$th correlated fading random variable can be parametrized as
$g_{kl} = \sqrt{1-\mu_k^2} x_{kl} + \mu_{kl} x_{0l} + j\left(\sqrt{1-\mu_k^2} y_{kl} + \mu_{kl}y_{0l} \right), l = 2,\ldots,m,$
where $x_{kl}$ and $y_{kl}$ are independent Gaussian RVs with zero mean and variance $1/2.$
The correlation coefficient, denoted by $0\lt\mu_k\lt1$ is found as $\mu_k = J_0\left(\frac{2\pi(k-1)\alpha}{N-1}\right)$ for $k = 2,\ldots,N,$ where $J_0(\cdot)$ is the zero order Bessel function of the first kind.
I would like to get the joint pdf of $|g_1|^2,|g_2|^2,\ldots,|g_N|^2$.

The above is my original posted question. below is my more precise question
I have $N$ flat fading channels experienced by $N$ ports on a user equipment that are placed closer than the $\lambda/2$ distance meaning there is spatial correlation among them (as given above by the zero order Bessel function) and the amplitudes of the channels are Rayleigh distributed each with pdf $P_{|g_k|}(r) = \frac{2r}{\sigma^2}e^{-\frac{r^2}{\sigma^2}},$ for $r \geq 0$ with $E[|g_k|^2]= \sigma^2$. Also the average SNR is defined as $\Gamma = \frac{E[|g_k|^2] E[|x|^2]}{\sigma_\eta^2}.$
The user equipment is such that it can select the best port with the strongest signal, i.e., $|g| = max{(|g_1|,|g_2|,...,|g_N|)}.$ This is equivalent to selection combining technique.
If the instantaneous SNR at any of the $N$ branches is given as $\gamma_i = \frac{A^2|g_i|^2}{N_i}$ where $A$ is the amplitude of the transmitted signal and $N_i$ is the noise spectral density, the selection combiner would choose the branch with the maximum instantaneous SNR for a symbol decision. The output SNR would be equal to
$|\gamma| = max{(|\gamma_1|,|\gamma_2|,...,|\gamma_N|)}.$
In my original question i have a joint pdf of $|g_1|,|g_2|,\ldots,|g_N|$. I wanted a (preferably, in closed form) a joint pdf of $|g_1^2|,|g_2^2|,\ldots,|g_N^2|$ and then wanted to transform that to the joint pdf of $|\gamma_1|,|\gamma_2|,\ldots,|\gamma_N|$
The previous answer that I got helped me to cancel out terms that were the reciprocal of the jacobian which just left me with a product of N exponentials and a product of N zero order modified Bessel functions.
So I  am asking for a transformation of my joint pdf in terms of channel coefficients into joint pdf of snrs.
I Hope the way I posed my question is not confusing for the experts out there whose help I seek. Thank you
 A: $\newcommand\p\partial\newcommand\ga\gamma$By the substitution formula for multiple integrals, the joint pdf of $|g_1|^2,\dots,|g_N|^2$ is given by
$$P_{|g_1|^2,\dots,|g_N|^2}(s_1,\dots,s_N)=P_{|g_1|,\dots,|g_N|}(s_1^{1/2},\dots,s_N^{1/2})J(s_1,\dots,s_N)$$
for $s_1>0,\dots,s_N>0$,
where
$$J(s_1,\dots,s_N)=\frac1{2^N}\prod_{j=1}^N s_j^{-1/2}$$
is the (here positive) determinant of the (here diagonal) Jacobian matrix
$$\frac{\p(s_1^{1/2},\dots,s_N^{1/2})}{\p(s_1,\dots,s_N)}
:=\Big(\frac{\p(s_i^{1/2})}{\p s_j}\Big)_{i,j=1}^N$$
of the transformation $(s_1,\dots,s_N)\mapsto(s_1^{1/2},\dots,s_N^{1/2})$, which latter is inverse to the transformation $(t_1,\dots,t_N)\mapsto(t_1^2,\dots,t_N^2)$.
So, if you know $P_{|g_1|,\dots,|g_N|}$, it is very easy to get $P_{|g_1|^2,\dots,|g_N|^2}$.

The OP has changed the question -- now wanting the joint pdf of $\ga_1,\dots,\ga_N$ instead of the joint pdf of $|g_i|^2,\dots,|g_N|^2$, where $\ga_i=\frac{A^2|g_i|^2}{N_i}$.
By the same substitution formula for multiple integrals,
the joint pdf of $\ga_1,\dots,\ga_N$ is given by
$$P_{\ga_1,\dots,\ga_N}(r_1,\dots,r_N)
=P_{|g_i|^2,\dots,|g_N|^2}(\tfrac{N_1}{A^2}\,r_1,\dots,\tfrac{N_N}{A^2}\,r_n)K=P_{|g_1|,\dots,|g_N|}(\tfrac{N_1^{1/2}}A\,r_i^{1/2},\dots,\tfrac{N_1^{1/2}}A\,r_i^{1/2})J(N_1 r_1/A^2,\dots,N_N r_N/A^2) K$$
for $r_1>0,\dots,r_N>0$,
where
$$
K=\frac1{A^{2N}}\prod_{j=1}^N N_j$$
is the (here positive) determinant of the (here diagonal) Jacobian matrix
$$\frac{N_1 r_1/A^2,\dots,N_N r_N/A^2}{\p(r_1,\dots,r_N)}
$$
of the transformation $(r_1,\dots,r_N)\mapsto(N_1 r_1/A^2,\dots,N_N r_N/A^2)$, which latter is inverse to the transformation $(r_1,\dots,r_N)\mapsto
(A^2s_1/N_1,\dots,A^2s_N/N_N)$.
So,
$P_{\ga_1,\dots,\ga_N}(r_1,\dots,r_N)=$
$$\frac1{(2A)^N}
\prod_{k=1}^N \left(
\frac{N_k}{A\sigma^2(1-\mu_k^2)} 
\exp\left(-\frac{N_k r_k + \mu_k^2 N_1r_1}{A^2\sigma^2(1-\mu_k^2)} \right)\right.$$
$$\left.\times I_0 \left(\frac{2\mu_k(N_1N_2r_1r_2)^{1/2}}
{A^2\sigma^2(1-\mu_k^2)}\right) 
\right) $$
for $r_1>0,\dots,r_N>0$.
