Fully coupled ocean-atmosphere models are needed to represent and understand the complicated interactions, becoming increasingly important in climate change assessment in recent years. Numerical stability issues may arise because of cost-effective time integration. In particular, the contributing factors include using large time stepsize, lack of accurate interface flux, and singe-iteration coupling. We investigate the stability of the coupled ocean-atmosphere models for a variety of interface conditions such as DIrichlet-Neumman condition and bulk condition which is unique to climate modelling. We will also discuss the use of Schwarz-in-time iterative schemes that can add implicitness to the interface points for better stability. The efficiency can be achieved by using modified interface conditions that lead to faster convergence and by exploiting task-level parallelism.