Build discrete IIR from a transfer function which is not in terms of $s$ Firstly I should state that my maths knowledge is limited.

I'm currently designing an acoustics modelling application which uses a rectilinear FDTD grid for modelling pressure variations - but I've become stuck upon trying to implement a fibrous layer boundary.  From Kowalcyzk 2008, the fibrous boundary model is defined as a frequency domain continuous impedance filter $\xi_w(\omega)$:

\begin{align*}
 \rho(\omega) = 1.2 + [-0.0364({{\rho f} \over {\sigma}})^{-2} - j0.01144({{\rho f } \over {\sigma}})^{-1}]^{{1} \over {2}}
\end{align*}

\begin{align*}
 K(\omega) = 101320{{j29.64 + [2.82({{\rho f} \over {\sigma}})^{-2} + j24.9({{\rho f } \over {\sigma}})^{-1}]^{{1} \over {2}}} \over {j21.17 + [2.82({{\rho f} \over {\sigma}})^{-2} + j24.9({{\rho f } \over {\sigma}})^{-1}]^{{1} \over {2}}}}
\end{align*}

\begin{align*}
 Z _w(\omega) = [\rho(\omega)K(\omega)]^{{1} \over {2}}
\end{align*}

\begin{align*}
 \Gamma(\omega) = j2\pi f[{{\rho(\omega)} \over {K(\omega)}}]^{{1} \over{2}}
\end{align*}

\begin{align*}
 \xi _w(\omega) = {{Z_w(\omega)} \over {\rho c}} \coth [\Gamma(\omega)d]
\end{align*}

Boundaries are implemented in the grid edge nodes as discrete time domain IIR filters of arbitrary order.  The few resources I've read on the subject have transfer functions that consist of a ratio of polynomials in terms of $s$ or $z$ - so I have no idea where to start on converting this filter, any clues?

I read this post, but I don't know if there's simpler solution for my circumstance.
There exist genetic and evolutionary algorithms that attempt to guess and optimize a set of pole/zero placements that approximate an arbitrary frequency response.

There is also a algorithm called Frequency Domain Least Squares : http://www2.units.it/ramponi/teaching/DSP/materiale/ES_5_2.m 

The poles and zeros found can then be used to create an IIR recurance function.