See the FORTRAN code for a nonlinear 4box THC model used by

Y. Kaspi, R. Sayag, and E. Tziperman. A ‘triple sea-ice state’ mechanism for the abrupt warming and synchronous ice sheet collapses during heinrich events. Paleoceanography, 19(3):PA3004, 10.1029/2004PA001009, 2004.
http://www.seas.harvard.edu/climate/eli/reprints/Kaspi-Sayag-Tziperman-2004.pdf

here:
http://www.seas.harvard.edu/climate/eli/Downloads/Kaspi-etal-2003.tar