When looking for a copula I could simulate from (rather than the Gaussian copula), I found an algorithm for the Ali–Mikhail–Haq copula

that was proposed by Kumar (2010) as reproduced above. But the method seriously fails in that the range of (U,V) resulting from the simulation does not even cover the (0,1)² square! Unless I made an R coding mistake (which is always a possibility).

sim=function(T=1e3,h=.5){ 
  o=matrix(runif(2*T),T,2)
  a=1-o[,1];b=1-h*(1+2*a*o[,2])+2*h^2*a^2*o[,2]
  d=1+h*(2-4*a+4*a*o[,2])+h^2*(1-4*a*o[,2] +4*a^2*o[,2])
  o[,2]=1-2*o[,2]*(a*h-1)^2/(b+sqrt(d))
  return(o)}

There is no explanation in the paper as to why this algorithm is (not!) working, besides the inverse cdf argument—with the reference in R copBasic temporarily worrying until I checked the cdf inversion is completely numerical—, but a correct version can be derived from inverting the conditional cdf of one component, V, given the other, U. Namely, since the conditional cdf is given by

which leads to a second degree polynomial equation (in v) when solving the equation F(v|U=u) = w.

sim=function(T=1e3,h=.5){
  o=matrix(runif(2*T),T,2)
  v1=o[,1];w=o[,2]
  d=(2*h*v1*w-1-h)^2-4*(1-w)*h*(1-h*w*v1^2)
  o[,2]=-2*h*v1*w+1+h-sqrt(d))/(2*h*(1-h*w*v1^2)
  return(1-o)}

And with a more likely outcome (Xed checked by comparing F(u,v) with its empirical version for several pairs (u,v)):