What about $\int_0^1\int_0^1\frac{\text{gd}(\log(xy))}{1-xy}dxdy,$ where $\text{gd}(z)$ is the Gudermannian function? This morning I was thinking in integrals of the kind $$\int_0^1\int_0^1\frac{\text{Numerator}(x,y)}{1-xy}\,dxdy,$$ where $\text{Numerator}(x,y)$ is a function of $x$ and $y$ involving special functions. Then as a particular case of these I am trying to think if 

Question. Is it possible to calculate an expression (a series, a closed-form involving special functions or a very good approximation) for $$\int_0^1\int_0^1\frac{\text{gd}(\log(xy))}{1-xy}\,dxdy,$$ where $\text{gd}(z)$ is the Gudermannian function? What are your calculations and approach? If you don't know this special function see this MathWorld. Thanks in advance.

Notice that I've choose this example with the purpose to combine with the expression involving the Guermannian function and the inverse tangent function, $(4)$ in previous reference. Then using a Cauchy product I wrote 
$$\int_0^1\int_0^1\frac{\text{gd}(\log(xy))}{1-xy}\,dxdy=2\sum_{k=0}^\infty\sum_{l=0}^k\frac{(-1)^l}{(2l+1)(l+k+2)^2}-\frac{\pi^3}{12}.$$ On the other hand Wolfram Alpha online calculator knows how calculate indefinite integrals as  
int arctan(xy)/(1-xy)dx
but with my standard time of computation and code I am not able to calculate a closed-form for our integral in the Question.
 A: We have $\text{gd}\log z=2\arctan z-\frac{\pi}{2}$, hence
$$ \iint_{(0,1)^2}\frac{\text{gd}\log(xy)}{1-xy}\,dx\,dy=2\iint_{(0,1)^2}\frac{\arctan(xy)}{1-xy}\,dx\,dy-\frac{\pi}{2}\sum_{n\geq 0}\iint_{(0,1)^2}(xy)^n\,dx\,dy$$
where the last integral clearly equals $\zeta(2)=\frac{\pi^2}{6}$. Since
$$ \arctan(z)=\sum_{n\geq 0}\frac{(-1)^n z^{2n+1}}{2n+1} $$
we have
$$ 2\iint_{(0,1)^2}\frac{\arctan(xy)}{1-xy}\,dx\,dy = 2\sum_{n\geq 0}\frac{(-1)^n}{2n+1}\sum_{m\geq 2n+2}\frac{1}{m^2}=\frac{\pi^3}{12}-\color{blue}{2\sum_{n\geq 0}\frac{(-1)^n H_{2n+1}^{(2)}}{2n+1}} $$
and the whole problem boils down to evaluating the blue series. We have
$$ \sum_{n\geq 1} H_{n}^{(2)} z^n = \frac{\text{Li}_2(z)}{1-z},\qquad \sum_{n\geq 0} 2\,H_{2n+1}^{(2)} z^{2n+1} = \frac{\text{Li}_2(z)}{1-z}-\frac{\text{Li}_2(-z)}{1+z} $$
hence the evaluation of the blue series is equivalent to the computation of the integral
$$ \frac{1}{i}\int_{0}^{i}\left(\frac{\text{Li}_2(z)}{z(1-z)}-\frac{\text{Li}_2(-z)}{z(1+z)}\right)\,dz$$
which by partial fraction decomposition and the functional identities for the dilogarithm and trilogarithm equals:
$$ 2\sum_{n\geq 0}\frac{(-1)^n H_{2n+1}^{(2)}}{2n+1}=\color{blue}{-\frac{\pi^3}{24}+G \log 2-\frac{\pi}{8}\log^2(2)+4\,\text{Im }\text{Li}_3\left(\frac{1+i}{2}\right)} $$
where $G$ is Catalan's constant.
