\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
2006 International Conference in Honor of Jacqueline Fleckinger.
\newline \emph{Electronic Journal of Differential Equations},
Conference 16, 2007, pp. 129--135.
\newline ISSN: 1072-6691.
URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu (login: ftp)}
\thanks{\copyright 2007 Texas State University - San Marcos.}
\vspace{8mm}}

\begin{document} \setcounter{page}{129}

\title[\hfilneg EJDE/Conf/16 \hfil  2D climate EBM coupled with 3D deep
ocean model]
{A 2D climate energy balance model\\
 coupled with a 3D deep ocean model}

\thanks{J.~I. D\'{\i}az was supported by
         project MTM2005-03463 from  DGISGPI, Spain}

\author[J. I. D\'{\i}az,  L. Tello \hfil EJDE/Conf/16 \hfilneg]
{Jes\'us Ildefonso D\'{\i}az, Lourdes Tello} % in alphabetical order

\address{Jes\'us Ildefonso D\'{\i}az \newline
Departamento Matem\'atica Aplicada,
Universidad Complutense de Madrid,
E-28040 Madrid, Spain}
\email{diaz.racefyn@insde.es}

\address{Lourdes Tello \newline
Departamento Matem\'atica Aplicada,
ETS.\ Arquitectura, Universidad  Polit\'ecnica de Madrid,
E--28040 Madrid, Spain}
\email{l.tello@upm.es}

\thanks{Published May 15, 2007.}
\subjclass[2000]{35K55, 35K57, 58J32}
\keywords{Climate atmosphere - ocean models; PDE's on manifolds
 \hfill\break\indent dynamic and diffusive boundary condition}

\dedicatory{Dedicated to  Jacqueline Fleckinger on the occasion \\
of an international conference in her honor}

\begin{abstract}
 We study a three dimensional climate model which represents
 the coupling of the mean surface temperature with the
 ocean temperature.  We prove the existence of a bounded
 weak solution by a  fixed point argument.
\end{abstract}

\maketitle

\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{definition}[theorem]{Definition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{lemma}[theorem]{Lemma}

\section{Introduction}
 In the last decades, several works about global climate energy
balance models (EBM) have appeared which study the evolution of the
mean surface temperature of the Earth (see for example, D\'{i}az
\cite{diaz1993}, D\'{i}az and Tello \cite{dt1999}, Ghil and
Childress \cite{gc}, Hetzer \cite{gh}, North \cite{grn}, etc.). From
the mathematical point of view, two dimensional EBM (latitude --
longitude) has an spatial domain given by a Riemannian manifold
without boundary $\mathcal{M}$ simulating the Earth surface, as
follows
%
\begin{equation}
\begin{gathered}
c(x)u_{t}-\mathop{\rm div}(k(x)|\nabla u|^{p-2}\nabla u)
 +R_{e}(x,u)\in R_{a}(x,u)
\quad (0,T)\times \mathcal{M}, \\
u(x,0)=u_{0}(x) \quad \mathcal{M},
\end{gathered} \label{ebm-2d}
\end{equation}
%
where $u$ represents the mean surface temperature, $R_{e}$ and
$R_{a}$ the emitted and absorbed energy, respectively. $R_{a}$
depends on the planetary coalbedo $\beta $ (which is eventually
discontinuous on $u$). One dimensional models (early proposed)
assume uniform temperature over each parallel. By calling $x$ the
sine of the latitude and doing a change to spherical coordinates, we
obtain
%
\begin{equation}
\begin{gathered}
c(x)u_{t}-(k(x)(1-x^{2})^{p/2}|u_{x}|^{p-2}u_{x})_{x}+R_{e}(x,u)\in
R_{a}(x,u) \quad (0,T)\times (-1,1), \\
(1-x^{2})^{p/2}|u_{x}|^{p-2}u_{x}=0 \quad x\in \{-1,1\}, \\
u(x,0)=u_{0}(x) \quad (-1,1),
\end{gathered}  \label{ebm-1d}
\end{equation}
%
In these models, the effect of the oceans is only considered in a
implicit and empirical way in the spatial dependence of the
coefficients. However, some works about the rapid climatic change in
Glacial-Holocene transition (see p.e. Berger et al \cite{ber}) show
that the cause could be the past changes in deep water formation. In
this work we study a model including the effect of the deep Ocean
based in the model proposed by Watts - Morantine \cite{ww}. In such
3D model, the atmosphere temperature comes to an energy balance
model closed to (\ref{ebm-2d}).

\section{The Mathematical Model}

The model represents the evolution of the temperature in a global
ocean $\Omega$ with constant depth $H$. The upper boundary of
$\Omega$ is a manifold $\mathcal{M}$ simulating the Earth surface.
The bottom of the ocean is a manifold $\mathcal{N}$. For the
mathematical treatment, we assume that $\mathcal{M}$ and $\mathcal{N}$
are $C^{\infty} $ two dimensional compact connected oriented
Riemannian Manifold of $\mathbb{R}^3$ without boundary and
dist$(\mathcal{M},\mathcal{N})=H$. For example $\mathcal{M}$ and
$\mathcal{N}$ can be two spheres with the same center and radius
$R$ and $R-H$, respectively.

The shape of the spatial domain $\Omega$ suggest us to use a
suitable coordinate system $(\theta , \varphi,z)$ where $z\in
(-H,0)$ is measured from $\mathcal{M}$ to $\mathcal{N}$ in the
orthogonal direction to the tangent space $T_p\mathcal{M}$.

The governing equation for the ocean interior is a heat equation
with transport,
%
\begin{equation} \label{eq-oceano}
\frac{\partial U}{\partial t}-\mathop{\rm div}(\nabla U)+w\frac{\partial
^{{}}U}{\partial z^{{}}} =0 \quad  (t,x) \in  (0,T)\times \Omega ,
\end{equation}
%
where $U$ is the temperature, t the time, the variable $z\in (-H,0)$
represents the depth and $w$ is the vertical velocity.

At the bottom of the ocean ($z=-H$), the  advective and diffusive
transports must be equal
%
\begin{equation}
\hat{F}(x,\nabla _{\mathcal{N}} U)+{\frac{\partial U}{
\partial z}}  =0\quad  \mbox{in }(0,T)\times \mathcal{N} ,
\end{equation}
%
where $\hat{F}$ is linear on the gradient $\nabla _{\mathcal{N}}U$
and the gradient is understood in the sense  of the Riemannian
metric of $\mathcal{N}$. Analogously for the $F$ and  $\mathcal{M}$
below.

The boundary condition at $z=0$ comes from an energy balance
%
\begin{equation}
{\frac{\partial u}{\partial t}}-\mathop{\rm div}( |\nabla _{\mathcal{M}}
u|^{p-2} \nabla _{\mathcal{M}}u) +{\frac{\partial U}{
\partial n}}+F(x,\nabla _{\mathcal{M}}u)+R_e(u)
\in \frac{1}{\rho c}QS(x)\beta (U)
\end{equation}
%
where $u:\mathcal{M}\to \mathbb{R}$ and
\begin{equation}
 U\big|_{\mathcal{M}}=u.
\end{equation}
%
Here, $u$ represents the temperature on $\mathcal{M}$. $R_{e} $ is
the energy emitted by the surface $\mathcal{M}$ to the outer space
 by cooling. The diffusion operator for
$p=2$ is the Laplace - Beltrami operator on the manifold $\mathcal{M}$. The pioneering climate energy balance models considered the
linear case $p=2$, but, later,  Stone \cite{phs} proposed the
case $p=3$ in order to consider the negative feedback in the eddy
fluxes. In order to include both cases we consider $p\geq 2$. $\beta
(u) $ the coalbedo function (a nondecreasing function of $u$ of the
type $\beta (u)=0.7$ if $u>-10+\epsilon $, $\beta (u)=0.4$ if
$u<-10-\epsilon$ with $\epsilon \geq 0$). The coalbedo function will
be treated in the general class of multivalued graphs. More
precisely, we shall assume that $ \beta $ is a maximal monotone
graph of $ \mathbb{R}^2$ such that $\beta$ is bounded (i.e. $ m\leq
b \leq M$ for any $b\in \beta (s)$  for any $ s\in \mathbb{R} $). We
recall that $ \beta $ was assumed to be $ multivalued$ at $u=-10$ by
Budyko \cite{budyko1969} and $ \beta \ locally \ Lipschitz $ by
Sellers \cite{sellers1969}. The function $S:\mathcal{M}\to
\mathbb{R}$ is the insolation function and $Q>0$ the Solar constant.
This kind of models are very sensitive to small fluctuations of
Solar and terrestrial parameters, see Hetzer \cite{gh},
Arcoya, Diaz, Tello \cite{adt}, Diaz, Tello \cite{dt2006} and the references
therein. The equation describing the energy balance at $\mathcal{M}$ is similar to the equation studied in Diaz, Tello \cite{dt1999}
except for the two terms describing the thermal communication with
the deep ocean through the advective and diffusive transports of
heat.

We also need to add the initial conditions at $t=0$,
%
\begin{equation} \label{initialdata}
U(0,\cdot)=U_0(\cdot) \quad \mbox{and} \quad
u(0,\cdot)=u_0(\cdot).
\end{equation}
We notice that $u$ is also an unknown in this model.

The model can be simplified by considering that $\mathcal{M}$ is
the sphere of radius $R$ and the temperature is constant over each
parallel. In this case, the spatial domain $\Omega $ is reduced to
a rectangle with boundary $\Gamma _{0}\cup \Gamma _{H}\cup \Gamma
_{1}$ and the model obtained $(P_{2D}$) is
%
\begin{gather*}
{\frac{\partial U}{\partial
t}}-{\frac{K_{H}}{R^{2}}}\frac{
\partial }{\partial x}((1-x^{2}){\frac{\partial U}{\partial x})-}K_{v}{\frac{
\partial ^{2}U}{\partial z^{2}}}+w\frac{\partial ^{{}}U}{\partial z^{{}}}
=0 \quad \text{in }(0,T)\times \Omega ,
\\
 wx{\frac{\partial U}{\partial x}}+K_{V}{\frac{\partial U}{
\partial z}} =0 \quad \text{in }(0,T)\times \Gamma _{H},
\\
\begin{aligned}
&D{\frac{\partial U}{\partial t}}-{\frac{DK_{H_{0}}}{R^{2}}}
{\frac{\partial }{\partial x}}\left( (1-x^{2})^{p/2}
|{\frac{\partial U}{\partial x}}|^{p-2}{\frac{\partial U}{\partial x}}\right)
+ K_{V}{\frac{\partial U}{\partial n}}\\
& +wx {\frac{\partial U}{\partial x}}+\mathcal{G}(U) \in \frac{1}{
\rho c}QS(x)\beta (x,U) \quad \text{in }(0,T)\times \Gamma _{0},
\end{aligned}
\\
(1-x^{2})^{p/2}|{\frac{\partial U}{\partial
x}}|^{p-2}{\frac{\partial U}{\partial x}} {=0} \quad
 \text{in }(0,T)\times \Gamma _{1} , \\
U(0,x,z)=U_{0}(x,z)   \quad \text{in }\Omega , \\
U(0,x,0)=u_{0}(x)  \quad \text{in }\Gamma _{0},
\end{gather*}
%
where $x$ represents the sine of the latitude. The physical
description and the numerical approximation if $p=2$ and where
$\beta$ depends only on $x$ is in Watts Morantine \cite{ww}. The
proof of existence of bounded weak solution to this 2-D model is in
Diaz and Tello \cite{cedya-salamanca}. The number of steady states
of $(P_{2D})$ was studied in \cite{dt2006}.

We notice that in the 3D model some physical constants were assumed
equal to one.

The goal of the present work is to prove the existence of weak
solutions to the 3D evolution model with dynamical and diffusive
boundary condition described in
(\ref{eq-oceano})-(\ref{initialdata}).


\section{Existence of weak solutions}

The model represents the temperature evolution in a global ocean
$\Omega$.  The unknown are $U:\Omega \to \mathbb{R}$ and
$u:\mathcal{M} \to \mathbb{R}$. The independent spatial
variables are $(\theta, \varphi , z)$  where $z\in (-H,0)$. See
\cite{dt1999} and \cite{ltw1992} to the expression of the operators
$\nabla$ and $\mathop{\rm div}$ in the new coordinates.
Let (P) be the problem
\begin{equation} \label{P}
\begin{gathered}
{\frac{\partial U}{\partial t}} -\mathop{\rm div}(\nabla
U)+w\frac{\partial ^{{}}U}{\partial z^{{}}}
=0 \quad
\text{in }(0,T)\times \Omega ,
\\
\begin{aligned}
& {\frac{\partial u}{\partial t}}-\mathop{\rm div}(|\nabla _{\mathcal{M}}
u|^{p-2}  \nabla _{\mathcal{M}}u) +{\frac{\partial U}{\partial n}}
 +F(x,\nabla _{\mathcal{M}}u)+\mathcal{G}(u)\\
&\in \frac{1}{\rho c}QS(x)\beta (U)\quad \text{in }(0,T)\times \mathcal{M}
\end{aligned} \\
U\big|_{\mathcal{M}}=u
\\
 \hat{F}(x,\nabla _{\mathcal{N}} U)+{\frac{\partial
U}{\partial z}}  =0 \quad  \text{in }(0,T)\times  \mathcal{N}
\\
U(0,x,z)=U_{0}(x,z) \quad   \text{in }\Omega , \\
u(0,\cdot)=u_{0}(\cdot)  \quad \text{on } \mathcal{M},
\end{gathered}
\end{equation}
%
where
\begin{itemize}
\item[(H1)] $\Omega$ is a bounded and open set of $\mathbb{R} ^3$
with constant depth $H$ and $\partial \Omega =\mathcal{M} \cup
\mathcal{N}$. $\mathcal{M}$ and $\mathcal{N}$ are $C^{\infty} $
two dimensional compact connected oriented Riemannian Manifold of
$\mathbb{R}^3$ without boundary and $\mathop{\rm
dist}(\mathcal{M},\mathcal{N})=H$.
\end{itemize}

Here $\nabla_\mathcal{M}$ and $\mathop{\rm div}$ are  understood
in the sense of the Riemannian metric on $\mathcal{M}$. We study
the existence of solutions of \eqref{P} under the following
structure hypotheses:

\begin{itemize}
%
\item[(H2)]  $\beta $ is a bounded maximal
monotone graph, i.e. $|v|\leq M$  for all $v\in \beta (s)$,
and all $s\in D(\beta )=\mathbb{R}$.

\item[(H3)]  $\mathcal{G}:\mathbb{R}\to \mathbb{R}$
is a continuous strictly increasing function such that
$\mathcal{G}(0)=0$ and $|\mathcal{G}(\sigma )|\geq C|\sigma |^{r}$
for some $r>0$.

\item[(H4)]  $S:\mathcal{M} \to \mathbb{R}$,
$s_{1}\geq S(x)\geq s_{0}>0$ a.e. $x\in \mathcal{M}$.

\item[(H5)]  $f\in L^{\infty }((0,T)\times \Omega )$,

\item[(H6)]  $F:\mathcal{M}\times T\mathcal{M}\to
\mathbb{R}$ and $\hat{F}:\mathcal{N}\times T\mathcal{N}\to \mathbb{R}$
 are linear on the tangent bundle spaces
$T\mathcal{M}$ and $T\mathcal{N}$ with bounded coefficients.

\item[(H7)]  $w\in C^{1}(\overline{\Omega })$.
\end{itemize}


We define the functional space on the Riemannian manifold $\mathcal{M}$, as follows,
%
$$
V:= \{ u \in L^2(\mathcal{M}): \nabla_{\mathcal{M}} u \in
L^p(T\mathcal{M})\}
$$
where $T\mathcal{M}=\cup_{p\in \mathcal{M}} T_pM$ is the tangent
bundle space (see Aubin \cite{aubin}).



\begin{definition} \rm
We say that the pair $(U,u)\in C([0,T];L^2(\Omega)\times
L^2(\mathcal{M}))$  is a  bounded weak solution of \eqref{P} if

\begin{itemize}
\item[(i)] $(U,u)\in L^{\infty }((0,T)\times \Omega) \times L^{\infty }((0,T)\times
 \mathcal{M}) \cap
L^2(0,T; H^1(\Omega ))\times L^p(0,T; V)$,

\item[(ii)] there exists
$ Z\in L^{\infty }((0,T) \times \mathcal{M})$ with
$ Z(t,x) \in \beta (u(t,x))$ a.e. $(t,x) \in (0,T)\times \mathcal{M} $
such that
\end{itemize}
\begin{align*}
&\int _{\Omega}U(T,x)\phi (T,x)dA  -\int_0^T
\langle \phi_t(t,x),U(t,x)\rangle_{H^1(\Omega)'\times H^1(\Omega)} dt\\
&+\int_0^T\int_{\Omega} \nabla U  \nabla \phi \,dA\,dt +
 \int_0^T\int_{\Omega}w
\frac{\partial U}{\partial z}\phi \,dA\,dt \\
&-\int_0^T\int_{\mathcal{M}}
\frac{\partial U}{\partial n}\phi \,dS\,dt
+ \int_0^T\int_{\mathcal{N}} \hat{F}(x,\nabla_\mathcal{N})\phi \,dS\,dt \\
&=\int _{\Omega}U_0(x)\phi (0,x)dA,
\end{align*}
and
\begin{align*}
&\int _{\mathcal{M}}u(T,x)\psi (T,x)dA  -\int_0^T
\langle \psi_t(t,x),u(t,x)\rangle_{V'\times V} dt \\
&+ \int_0^T\int_{\mathcal{M}} |\nabla u |^{p-2}\nabla u \nabla \psi
\,dS\,dt + \int_0^T\int_{\mathcal{M}} {\mathcal G}(u)\psi \,dS\,dt\\
&+\int_0^T\int_{\mathcal{M}} \frac{\partial U}{\partial n}\psi
\,dS\,dt+\int_0^T\int_{\mathcal{M}} F(x,\nabla_\mathcal{M})\psi \,dS\,dt \\
&= \int_0^T\int_{\mathcal{M}} QS(x)Z(t,x)\psi
\,dA\,dt+\int_0^T\int_{\mathcal{M}} f\psi \,dA\,dt
+\int_{\mathcal{M}}u_0(x) \psi (0,x)dS
\end{align*}
%
for every test function $( \phi,\psi ) \in L^2(0,T;H^1(\Omega))
\times L^p((0,T);W^{1,p} (\mathcal{M}))$ such that  $(\phi_t,\psi
_t) \in L^{2}(0,T; H^1(\Omega)') \times L^{p'}(0,T; V')$. Here
$<,> _{V'\times V} $ denotes the duality product in $V'\times V$.
\end{definition}

\begin{theorem} \label{th.existence}
Let $U_{0}\in L^{\infty }(\Omega)$ and $u_{0}\in L^{\infty }(\mathcal{M})$.
Then there exists at
least a bounded weak global solution of \eqref{P}.
\end{theorem}

The main idea for the proof is to construct  an operator
$\mathcal{T}$ and to find a fixed point which is the solution of \eqref{P}.
The proof consist of several steps.


\noindent{\bf Step 1}. For every $h\in L^{\infty }((0,T)\times
\mathcal{M})$ we consider the problem $(P_{h})$ by replacing the
coalbedo term in \eqref{P} by $h$. The proof of the existence of
solution of $(P_{h})$ is inspired in Diaz and Jimenez
\cite{dj1984} and Bejaranu, Diaz and Vrabie \cite{bejenaru2001}.
We define the vectorial operator $\mathcal{A}$
%
\begin{equation*}
\mathcal{A}(U,u)\longmapsto (AU,Bu)
\end{equation*}
%
on the domain
%
\begin{equation*}
D(\mathcal{A})=\{(U,u)\in L^{2}(\Omega )\times L^{2}(\mathcal{M})
:AU\in L^{2}(\Omega ),Bu\in L^{2}(\mathcal{M}),U_{|\mathcal{M}}=u\},
\end{equation*}
%
where
%
\begin{gather*}
AU=-\mathop{\rm div}(\nabla U)+w\frac{\partial ^{{}}U}{\partial z},
\\
Bu=-\mathop{\rm div}( |\nabla _{\mathcal{M}} u|^{p-2} \nabla _{\mathcal{M}}u)
+{\frac{\partial U}{
\partial n}}+F(x,\nabla _{\mathcal{M}}u)+\mathcal{G}(u).
\end{gather*}
%
The existence of solution of $(P_h)$ is a consequence of the
following properties of $\mathcal{A}$.

\begin{lemma} \label{A-properties}
There exists $\lambda_0>0$ such that for every $\lambda>\lambda_0$,
we have:
\begin{itemize}
\item[(i)] $\mathcal{A}+\lambda I$ is $T$-accretive in $L^{2}(\Omega )\times
L^{2}(\mathcal{M})$.

\item[(ii)] $R(\mathcal{A}+\lambda I)=L^{2}(\Omega )\times
L^{2}(\mathcal{M})$.
\end{itemize}
\end{lemma}

Note that (i) allows us to prove a comparison principle for the
system
\begin{equation} \label{Pfg}
\begin{gathered}
\lambda U+AU=f \quad \text{in }L^{2}(\Omega ) \\
\lambda u+Bu=g \quad \text{in }L^{2}(\mathcal{M}) \\
U_{|\mathcal{M}}=u   \\
\widehat{F} (x,\nabla_{\mathcal{N}}U)+{\frac{\partial U}{\partial z}}=0
\quad \mathcal{N}.
\end{gathered}
\end{equation}

In fact, if $f_{1}\leq f_{2}$ and $g_{1}\leq g_{2}$ then the
solutions of \eqref{Pfg} with $f=f_{1}$, $g=g_{1}$ and
of \eqref{Pfg} with $f=f_{2}$, $g=g_{2}$ satisfy
\begin{equation*}
U_{1}\leq U_{2}, \quad\text{and}\quad u_{1}\leq u_{2}.
\end{equation*}

\begin{remark} \rm
Lemma \ref{A-properties} and similar arguments to those
in \cite[Lemma 3]{dt1999} allow us to prove the existence
of a maximal and minimal solutions.
\end{remark}

To prove (ii) in Lemma \ref{A-properties} we notice that the
operator $B$ can be expressed as $B_{1}+B_{2}+B_{3}$, where $B_{1}$
and $B_{2}$\ are maximal monotone operators in $L^{2}(\mathcal{M})$,
%
\begin{equation*}
B_{1}u= -\mathop{\rm div}(|\nabla _{\mathcal{M}} u|^{p-2}
\nabla _{\mathcal{M}}u) +\mathcal{G}(u)
\end{equation*}
%
and the pseudo-differential operator
\begin{equation*}
B_{2}u={\frac{\partial U}{\partial n},}
\end{equation*}
where $U$ is the solution of the problem
\begin{gather*}
\lambda U+AU=f \quad \text{in }L^{2}(\Omega ) \\
U_{|\mathcal{M}}=u.
\end{gather*}
The operator $B_{3}$ is defined by
\begin{equation*}
B_{3}u=F(\nabla _{\mathcal{M}}u).
\end{equation*}
This operator is not necessarily monotone but it is dominated
(in some sense) by the operators $B_{1}$ and $B_{2}$. Consequently, it is
possible to apply the abstracts results of perturbation of maximal
monotone operators (see e.g. \cite[Proposition 2.10]{brrr}),
and we arrive to the conclusion.

\noindent {\bf Step 2.} We follow the proof in \cite[Theorem 3]{dt1999}.
 We define the operator $\mathcal{T}:h\to g$ where $g\in \beta (u_{h})$
and $u_{h}$ is the solution of $(P_{h})$. It is easy to see that
every fixed point of $\mathcal{T}$ is a solution of \eqref{P}.

We prove that $\mathcal{T}$ satisfies the hypotheses of
Kakutani fixed point Theorem (see p.e. Vrabie \cite{vrabie}). We
denote $X=L^{p}((0,T),L^{2}(\mathcal{M}))$ then

\begin{itemize}
\item[(i)]  $K=\{h\in L^{p}((0,T),L^{\infty }(\Omega )):||h(t)||\leq C_{0}$
a.e. $t\in (0,T)\}$ is a nonempty, convex and weakly compact set of
$X;$

\item[(ii)]  $\mathcal{T}:K\mapsto 2^{X}$ with nonempty, convex and closed
values such that $\mathcal{T}(g)\subset K$, $\forall g\in K;$

\item[(iii)]  graph($\mathcal{T}$) is weakly$\times $weakly sequentially
closed.
\end{itemize}

 Consequently, $\mathcal{T}$ has at least one fixed point
in $K$ which is a solution of \eqref{P}.
This completes the proof of Theorem \ref{th.existence}.

\begin{remark} \rm
The boundedness of the solutions is a consequence of the existence of
 super-solutions and sub-solutions which are constants.
\end{remark}

\begin{thebibliography}{00}

%
\bibitem{adt} D. Arcoya, J.I. D\'{\i}az, L. Tello,
     \emph{S-shaped bifurcation branch in a quasilinear multivalued model
     arising in Climatology},
     Journal of Differential Equations, 150 (1998) 215-225.
%
\bibitem{aubin} T. Aubin, \emph{Nonlinear Analysis on Manifolds. Monge-Ampere
Equations}. Springer-Verlag, New York, 1982.
%
\bibitem{bejenaru2001}
I. Bejenaru, J. I. Diaz, I. I. Vrabie;
\emph{An abstract approximate
controllability result and applications to elliptic and parabolic
systems with dynamic boundary conditions},
Electronic J. Diff. Eqns., {\bf 2001}, 50, (2001), 1-19.
%
\bibitem{ber}
W. H. Berger, S. Burker, E. Vincent, \emph{Glacial-Holocene
transition: Climate Pulsations and Sporadic Shutdown of NADW
production}, in Abrupt Climatic Change - Evidence and
Implications, (eds. W.H. Berger, L.D. Labeyrie), Reidel Publishing
Co. Dordrecht Holland (1987).
%
\bibitem{brrr} H. Brezis, \emph{
Operateurs maximaux monotones et semigroupes de contractions dans
les espaces de Hilbert}. North Holland, Amsterdam (1973).
%
\bibitem{budyko1969}
M. I. Budyko, \emph{The effects of solar radiation variations on the
climate of the Earth}, Tellus, {\bf 21} (1969), 611-619.
%
\bibitem{diaz1993}  J. I. Diaz, \emph{Mathematical analysis of some diffusive
energy balance climate models}, in the book Mathematics,
Climate and Environment, (J. I. Diaz and J. L. Lions, eds. Masson,
Paris, 28-56 (1993).
%
\bibitem{dj1984}%
J. I. Diaz,  R. Jimenez, \emph{Aplicaci\'on a la teoria no lineal de
semigrupos a un operador pseudodiferencial}, Actas VII CEDYA,
Univ. Granada (1984) 137-142.
%
\bibitem{dt1999}  J. I. Diaz, L. Tello,
\emph{A nonlinear parabolic problem on a Riemannian manifold without
boundary arising in Climatology}, Collectanea Mathematica
{\bf50},1 (1999), 19-51.
%
\bibitem{cedya-salamanca} J. I. D\'{\i}az, L. Tello,
\emph{Sobre un modelo climatico
de balance de energia superficial acoplado con un oceano profundo},
 Actas del XVII CEDYA/ VI CMA, (2001).
%
\bibitem{dt2006}  J. I. Diaz, L. Tello,
\emph{On a climate model with a dynamic nonlinear diffusive boundary
condition}, submitted (2006).
%
\bibitem{gc} M. Ghil, S. Childress,
\emph{Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics
Dynamo Theory and Climate Dynamics}.  Springer Verlag. Applied
Mathematical  Sciences. 1987.
%
\bibitem{gh} G. Hetzer, \emph{The structure of the principal component
for semilinear diffusion equations from energy balance climate
models}, Houston Journal of Math. {\bf16}, 203-216 (1990).
%
\bibitem{ltw1992} J. L. Lions, R. Temam, S. Wang,
\emph{New formulations of the primitive
equations of atmosphere and applications}, Nonlinearity {\bf
5}, 237-288 (1992).
%
\bibitem{grn}  G. R. North,  \emph{Multiple solutions in energy balance
climate models}. Paleogeography, Paleoclimatology, Paleoecology
{\bf82}, Elsevier Science Publishers B.V. Amsterdam, 225-235 (1990).
%
 \bibitem{sellers1969} W. D. Sellers, \emph{A global climatic model
based on the energy balance of the earth-atmosphere system}, J.
Appl. Meteorol. {\bf 8} (1969),  392-400.
%
\bibitem{phs}  P. H. Stone, \emph{A simplified radiative - dynamical model
for the static stability of rotating atmospheres}, Journal of
the Atmospheric Sciences, {\bf29},  No. 3, 405-418 (1972).
%
\bibitem{vrabie}  I. I. Vrabie,  \emph{Compactness methods for
nonlinear evolutions}, Pitman Longman. London. 1987.
%
\bibitem{ww} R. G. Watts, M. Morantine,
\emph{Rapid climatic change and the deep ocean},
Climatic Change, {\bf16}, (1990) 83-97.

\end{thebibliography}

\end{document}