We present a computationally-efficient method for obtaining the fully spun-up state of a seasonally-varying global ocean biogeochemistry model. The solver uses a Newton–Krylov method to find the fixed points of the map that assigns to an initial state the value of the model state at the end of a one-year run. Apart from the preconditioner, which we describe in the paper, the method relies on a black-box public-domain Newton–Krylov solver that does not require the explicit construction of the model’s Jacobian matrix. Applied to the PO4 plus dissolved organic phosphorus (DOP) cycle of an Ocean Carbon Model Intercomparison Project II (OCMIP-2) type model, the solver is more than two orders of magnitude faster than the traditional time-stepping method for spinning up the model. The efficiency of the solver is illustrated by using the seasonally varying globally-gridded PO4 climatology to objectively optimize the parameters that control the mean lifetime of semi-labile DOP and the fraction of new production allocated to DOP. The optimization study demonstrates that the information in the seasonal variations of PO4 do not provide a significantly stronger constraint than the annually averaged data used in previous optimization studies.