--- manual/s_examples/baroclinic_gyre/fourlayer.tex 2001/08/08 16:15:41 1.1 +++ manual/s_examples/baroclinic_gyre/fourlayer.tex 2001/10/25 01:15:16 1.7 @@ -1,7 +1,8 @@ -% $Header: /home/ubuntu/mnt/e9_copy/manual/s_examples/baroclinic_gyre/fourlayer.tex,v 1.1 2001/08/08 16:15:41 adcroft Exp $ +% $Header: /home/ubuntu/mnt/e9_copy/manual/s_examples/baroclinic_gyre/fourlayer.tex,v 1.7 2001/10/25 01:15:16 cnh Exp $ % $Name: $ \section{Example: Four layer Baroclinic Ocean Gyre In Spherical Coordinates} +\label{sec:eg-fourlayer} \bodytext{bgcolor="#FFFFFFFF"} @@ -15,10 +16,11 @@ %{\large May 2001} %\end{center} -\subsection{Introduction} - -This document describes the second example MITgcm experiment. The first -example experiment ilustrated how to configure the code for a single layer +This document describes an example experiment using MITgcm +to simulate a baroclinic ocean gyre in spherical +polar coordinates. The barotropic +example experiment in section \ref{sec:eg-baro} +ilustrated how to configure the code for a single layer simulation in a cartesian grid. In this example a similar physical problem is simulated, but the code is now configured for four layers and in a spherical polar coordinate system. @@ -35,14 +37,14 @@ In this experiment the model is configured to represent a mid-latitude enclosed sector of fluid on a sphere, $60^{\circ} \times 60^{\circ}$ in lateral extent. The fluid is $2$~km deep and is forced -by a constant in time zonal wind stress, $\tau_x$, that varies sinusoidally -in the north-south direction. Topologically the simulated +by a constant in time zonal wind stress, $\tau_{\lambda}$, that varies +sinusoidally in the north-south direction. Topologically the simulated domain is a sector on a sphere and the coriolis parameter, $f$, is defined -according to latitude, $\phi$ +according to latitude, $\varphi$ \begin{equation} \label{EQ:fcori} -f(\phi) = 2 \Omega \sin( \phi ) +f(\varphi) = 2 \Omega \sin( \varphi ) \end{equation} \noindent with the rotation rate, $\Omega$ set to $\frac{2 \pi}{86400s}$. @@ -52,20 +54,22 @@ \begin{equation} \label{EQ:taux} -\tau_x(\phi) = \tau_{0}\sin(\pi \frac{\phi}{L_{\phi}}) +\tau_{\lambda}(\varphi) = \tau_{0}\sin(\pi \frac{\varphi}{L_{\varphi}}) \end{equation} -\noindent where $L_{\phi}$ is the lateral domain extent ($60^{\circ}$) and +\noindent where $L_{\varphi}$ is the lateral domain extent ($60^{\circ}$) and $\tau_0$ is set to $0.1N m^{-2}$. \\ Figure \ref{FIG:simulation_config} summarises the configuration simulated. -In contrast to example (1) \cite{baro_gyre_case_study}, the current -experiment simulates a spherical polar domain. However, as indicated +In contrast to the example in section \ref{sec:eg-baro}, the +current experiment simulates a spherical polar domain. As indicated by the axes in the lower left of the figure the model code works internally -in a locally orthoganal coordinate $(x,y,z)$. In the remainder of this -document the local coordinate $(x,y,z)$ will be adopted. +in a locally orthoganal coordinate $(x,y,z)$. For this experiment description +of this document the local orthogonal model coordinate $(x,y,z)$ is synonomous +with the spherical polar coordinate shown in figure +\ref{fig:spherical-polar-coord} \\ The experiment has four levels in the vertical, each of equal thickness, @@ -91,84 +95,165 @@ \noindent with $\rho_{0}=999.8\,{\rm kg\,m}^{-3}$ and $\alpha_{\theta}=2\times10^{-4}\,{\rm degrees}^{-1} $. Integrated forward in -this configuration the model state variable {\bf theta} is synonomous with +this configuration the model state variable {\bf theta} is equivalent to either in-situ temperature, $T$, or potential temperature, $\theta$. For consistency with later examples, in which the equation of state is non-linear, we use $\theta$ to represent temperature here. This is the quantity that is carried in the model core equations. \begin{figure} -\centerline{ +\begin{center} \resizebox{7.5in}{5.5in}{ \includegraphics*[0.2in,0.7in][10.5in,10.5in] {part3/case_studies/fourlayer_gyre/simulation_config.eps} } -} +\end{center} \caption{Schematic of simulation domain and wind-stress forcing function for the four-layer gyre numerical experiment. The domain is enclosed by solid walls at $0^{\circ}$~E, $60^{\circ}$~E, $0^{\circ}$~N and $60^{\circ}$~N. -In the four-layer case an initial temperature stratification is +An initial stratification is imposed by setting the potential temperature, $\theta$, in each layer. The vertical spacing, $\Delta z$, is constant and equal to $500$m. } \label{FIG:simulation_config} \end{figure} -\subsection{Discrete Numerical Configuration} +\subsection{Equations solved} - The model is configured in hydrostatic form. The domain is discretised with -a uniform grid spacing in latitude and longitude - $\Delta x=\Delta y=1^{\circ}$, so -that there are sixty grid cells in the $x$ and $y$ directions. Vertically the -model is configured with a four layers with constant depth, -$\Delta z$, of $500$~m. -The implicit free surface form of the -pressure equation described in Marshall et. al \cite{Marshall97a} is -employed. +The implicit free surface {\bf HPE} form of the +equations described in Marshall et. al \cite{Marshall97a} is +employed. The flow is three-dimensional with just temperature, $\theta$, as +an active tracer. The equation of state is linear. A horizontal laplacian operator $\nabla_{h}^2$ provides viscous -dissipation. The wind-stress momentum input is added to the momentum equation -for the ``zonal flow'', $u$. Other terms in the model +dissipation and provides a diffusive sub-grid scale closure for the +temperature equation. A wind-stress momentum forcing is added to the momentum +equation for the zonal flow, $u$. Other terms in the model are explicitly switched off for this experiement configuration (see section -\ref{SEC:code_config} ), yielding an active set of equations solved in this -configuration as follows +\ref{SEC:eg_fourl_code_config} ). This yields an active set of equations +solved in this configuration, written in spherical polar coordinates as +follows \begin{eqnarray} \label{EQ:model_equations} \frac{Du}{Dt} - fv + - \frac{1}{\rho}\frac{\partial p^{'}}{\partial x} - + \frac{1}{\rho}\frac{\partial p^{\prime}}{\partial \lambda} - A_{h}\nabla_{h}^2u - A_{z}\frac{\partial^{2}u}{\partial z^{2}} & = & -\cal{F} +\cal{F}_{\lambda} \\ \frac{Dv}{Dt} + fu + - \frac{1}{\rho}\frac{\partial p^{'}}{\partial y} - + \frac{1}{\rho}\frac{\partial p^{\prime}}{\partial \varphi} - A_{h}\nabla_{h}^2v - A_{z}\frac{\partial^{2}v}{\partial z^{2}} & = & 0 \\ -\frac{\partial \eta}{\partial t} + \nabla_{h}\cdot \vec{u} +\frac{\partial \eta}{\partial t} + \frac{\partial H \widehat{u}}{\partial \lambda} + +\frac{\partial H \widehat{v}}{\partial \varphi} &=& 0 +\label{eq:fourl_example_continuity} \\ \frac{D\theta}{Dt} - K_{h}\nabla_{h}^2\theta - K_{z}\frac{\partial^{2}\theta}{\partial z^{2}} & = & 0 +\label{eq:eg_fourl_theta} \\ -g\rho_{0} \eta + \int^{0}_{-z}\rho^{'} dz & = & p^{'} +p^{\prime} & = & +g\rho_{0} \eta + \int^{0}_{-z}\rho^{\prime} dz \\ -{\cal F} |_{s} & = & \frac{\tau_{x}}{\rho_{0}\Delta z_{s}} +\rho^{\prime} & = &- \alpha_{\theta}\rho_{0}\theta^{\prime} \\ -{\cal F} |_{i} & = & 0 +{\cal F}_{\lambda} |_{s} & = & \frac{\tau_{\lambda}}{\rho_{0}\Delta z_{s}} +\\ +{\cal F}_{\lambda} |_{i} & = & 0 \end{eqnarray} -\noindent where $u$ and $v$ are the $x$ and $y$ components of the -flow vector $\vec{u}$. The suffices ${s},{i}$ indicate surface and -interior model levels respectively. As described in -MITgcm Numerical Solution Procedure \cite{MITgcm_Numerical_Scheme}, the time -evolution of potential temperature, $\theta$, equation is solved prognostically. -The total pressure, $p$, is diagnosed by summing pressure due to surface -elevation $\eta$ and the hydrostatic pressure. -\\ +\noindent where $u$ and $v$ are the components of the horizontal +flow vector $\vec{u}$ on the sphere ($u=\dot{\lambda},v=\dot{\varphi}$). +The terms $H\widehat{u}$ and $H\widehat{v}$ are the components of the vertical +integral term given in equation \ref{eq:free-surface} and +explained in more detail in section \ref{sect:pressure-method-linear-backward}. +However, for the problem presented here, the continuity relation (equation +\ref{eq:fourl_example_continuity}) differs from the general form given +in section \ref{sect:pressure-method-linear-backward}, +equation \ref{eq:linear-free-surface=P-E+R}, because the source terms +${\cal P}-{\cal E}+{\cal R}$ +are all $0$. + +The pressure field, $p^{\prime}$, is separated into a barotropic part +due to variations in sea-surface height, $\eta$, and a hydrostatic +part due to variations in density, $\rho^{\prime}$, integrated +through the water column. + +The suffices ${s},{i}$ indicate surface and interior of the domain. +The windstress forcing, ${\cal F}_{\lambda}$, is applied in the surface layer +by a source term in the zonal momentum equation. In the ocean interior +this term is zero. + +In the momentum equations +lateral and vertical boundary conditions for the $\nabla_{h}^{2}$ +and $\frac{\partial^{2}}{\partial z^{2}}$ operators are specified +when the numerical simulation is run - see section +\ref{SEC:eg_fourl_code_config}. For temperature +the boundary condition is ``zero-flux'' +e.g. $\frac{\partial \theta}{\partial \varphi}= +\frac{\partial \theta}{\partial \lambda}=\frac{\partial \theta}{\partial z}=0$. + + + +\subsection{Discrete Numerical Configuration} + + The domain is discretised with +a uniform grid spacing in latitude and longitude + $\Delta \lambda=\Delta \varphi=1^{\circ}$, so +that there are sixty grid cells in the zonal and meridional directions. +Vertically the +model is configured with four layers with constant depth, +$\Delta z$, of $500$~m. The internal, locally orthogonal, model coordinate +variables $x$ and $y$ are initialised from the values of +$\lambda$, $\varphi$, $\Delta \lambda$ and $\Delta \varphi$ in +radians according to + +\begin{eqnarray} +x=r\cos(\varphi)\lambda,~\Delta x & = &r\cos(\varphi)\Delta \lambda \\ +y=r\varphi,~\Delta y &= &r\Delta \varphi +\end{eqnarray} + +The procedure for generating a set of internal grid variables from a +spherical polar grid specification is discussed in section +\ref{sec:spatial_discrete_horizontal_grid}. + +\noindent\fbox{ \begin{minipage}{5.5in} +{\em S/R INI\_SPHERICAL\_POLAR\_GRID} ({\em +model/src/ini\_spherical\_polar\_grid.F}) + +$A_c$, $A_\zeta$, $A_w$, $A_s$: {\bf rAc}, {\bf rAz}, {\bf rAw}, {\bf rAs} +({\em GRID.h}) + +$\Delta x_g$, $\Delta y_g$: {\bf DXg}, {\bf DYg} ({\em GRID.h}) + +$\Delta x_c$, $\Delta y_c$: {\bf DXc}, {\bf DYc} ({\em GRID.h}) + +$\Delta x_f$, $\Delta y_f$: {\bf DXf}, {\bf DYf} ({\em GRID.h}) + +$\Delta x_v$, $\Delta y_u$: {\bf DXv}, {\bf DYu} ({\em GRID.h}) + +\end{minipage} }\\ + + + +As described in \ref{sec:tracer_equations}, the time evolution of potential +temperature, +$\theta$, (equation \ref{eq:eg_fourl_theta}) +is evaluated prognostically. The centered second-order scheme with +Adams-Bashforth time stepping described in section +\ref{sec:tracer_equations_abII} is used to step forward the temperature +equation. The pressure forces that drive the fluid motions, ( +$\frac{\partial p^{'}}{\partial \lambda}$ and $\frac{\partial p^{'}}{\partial \varphi}$), are found by summing pressure due to surface +elevation $\eta$ and the hydrostatic pressure. The hydrostatic part of the +pressure is evaluated explicitly by integrating density. The sea-surface +height, $\eta$, is solved for implicitly as described in section +\ref{sect:pressure-method-linear-backward}. \subsubsection{Numerical Stability Criteria} @@ -249,7 +334,7 @@ stability limit of 0.25. \subsection{Code Configuration} -\label{SEC:code_config} +\label{SEC:eg_fourl_code_config} The model configuration for this experiment resides under the directory {\it verification/exp1/}. The experiment files