  Numerical models of Photodissociation Regions (PDRs) are an essential tool to
quantitatively understand observations of massive star forming regions through
simulations. Few mature PDR models are available and the Cologne KOSMA-$\tau$
PDR model is the only sophisticated model that uses a spherical cloud geometry
thereby allowing us to simulate clumpy PDRs. We present the current status of
the code as reference for modelers and for observers that plan to apply
KOSMA-$\tau$ to interpret their data. For the numerical solution of the
chemical problem we present a superior Newton-Raphson stepping algorithm and
discuss strategies to numerically stabilize the problem and speed up the
iterations. The chemistry in KOSMA-$\tau$ is upgraded to include the full
surface chemistry in an up-to-date formulation and we discuss a novel
computation of branching ratios in chemical desorption reactions. The high dust
temperature in PDRs leads to a selective freeze-out of oxygen-bearing ice
species due to their higher condensation temperatures and we study changes in
the ice mantle structures depending on the PDR parameters, in particular the
impinging UV field. Selective freeze-out can produce enhanced C abundances and
higher gas temperatures resulting in a fine-structure line emission of atomic
carbon [C] enhanced by up to 50% if surface reactions are considered. We show
how recent ALMA observations of HCO$^+$ emission in the Orion Bar with high
spatial resolution on the scale of individual clumps can be interpreted in the
context of non-stationary, clumpy PDR ensembles. Additionally, we introduce
WL-PDR, a simple plane-parallel PDR model written in Mathematica to act as
numerical testing environment of PDR modeling aspects.
