[R] heteroskedastic ordered probit

[Rendas-Baum, Regina]

Dear all,
 
is there existing code to fit a ordinal regression model with probit link where the variance can be modeled? For example, the polr function from the MASS package fits:
 
probitP(Y <= k | x) = zeta_k - eta

Instead of,
 
probitP(Y <= k | x) = (zeta_k - eta)/scale, where scale = exp(b*z) 
 
which is what I need.
 
I've seen these models go by the name of heteroskedastic ordered probit, location-scale glms, generalized nonlinear models, ... .  I know the package gnlm (Jim Lindsey) has a function to fit this model using the logit link, so I'm wondering if this is already done in R using probit instead.
 
Thanks.
 
Regina