Bivariate ordered logistic models (BOLMs) are appealing to jointly model the marginal distribution of two ordered responses and their association, given a set of covariates. When the number of categories of the responses increases, the number of global odds ratios (or their re-parametrizations) to be estimated also increases and estimating the association structure becomes crucial for this type of data. In fact, such data could be too "rich" to be fully modelled with an ordinary BOLM while, sometimes, the well-known Dale's model could be too parsimonious to provide a good fit. In addition, when the cross-tabulation of the responses contains some zeros, for a number of model configurations, including the bivariate version of the partial proportional odds model (PPOM), estimation of a BOLM by the Fisher-scoring algorithm may either fail or estimate a too "irregular" association structure. In this work, we propose to use a nonparametric approach for the maximum likelihood estimation of a BOLM. We apply penalties to the differences between adjacent row and column effects. As a result, estimation is less demanding than an ordinary BOLM, permitting the fit of PPOMs and/or the smoothing of the marginal and association parameters by polynomial curves and surfaces, with scores chosen by the data. Model selection is based on the penalized log-likelihood ratio, whose limiting distribution has been studied through simulations, and AIC. Our proposal is compared to the Goodman's model and the Dale's model, in terms of goodness-of-fit and parsimony, on a literature data set. Finally, an application on an original data set of liver disease patients is proposed.