<p><b>mpetersen@lanl.gov</b> 2012-07-26 16:49:55 -0600 (Thu, 26 Jul 2012)</p><p>Updating timesplitting design document to reflect state of code.  This includes updated details on namelist variables, and expands stage 2 to be more clear.<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/documents/ocean/current_design_doc/time_splitting_design/time_split_design.pdf
===================================================================
(Binary files differ)

Modified: trunk/documents/ocean/current_design_doc/time_splitting_design/time_split_design.tex
===================================================================
--- trunk/documents/ocean/current_design_doc/time_splitting_design/time_split_design.tex        2012-07-26 20:02:51 UTC (rev 2065)
+++ trunk/documents/ocean/current_design_doc/time_splitting_design/time_split_design.tex        2012-07-26 22:49:55 UTC (rev 2066)
@@ -12,6 +12,9 @@
 </font>
<font color="blue">ewcommand{\ds}{\displaystyle}
 \setlength{\parskip}{1.2ex}
 \setlength{\parindent}{0mm}
+</font>
<font color="blue">ewcommand{\bea}{\begin{eqnarray}}
+</font>
<font color="blue">ewcommand{\eea}{\end{eqnarray}}
+</font>
<font color="black">ewcommand{</font>
<font color="black">n}{</font>
<font color="gray">onumber}
 
 \begin{document}
 
@@ -43,12 +46,12 @@
 method differs from traditional splitting as the barotropic terms are
 subcycled explicitly rather than treated implicitly, as in POP.
 
-This document only addresses Higdon time-splitting in z-level
+This document only addresses split exlicit time stepping in z-level
 coordinates.  This will be applied to isopycnal coordinates at a later
-time.  
+time.  See the ALE vertical coordinate design document.
 
-The Higdon method consists of the following steps: decompose the
-velocity into barotropic and baroclinic components; take a large
+The split explicit method consists of the following steps: decompose
+the velocity into barotropic and baroclinic components; take a large
 timestep with the baroclinic velocities, computing the vertical mean
 forcing ${\overline {\bf G}}$; subcycle the barotropic velocity with
 small explicit timesteps; add velocities, compute other variables.
@@ -62,11 +65,11 @@
 
 \section{Requirement: A split time-stepping in z-level MPAS}
 Date last modified: 2011/05/4 \\
-Contributors: Mark, Chris, Todd \\
+Contributors: Mark, Todd \\
 
 The algorithm should follow Higdon (2005) section 2.3, but with alterations for z-level variables.
-Input variables should be provided for: timestepping type, number of Higdon iterations, number of barotropic subcycles, and number of baroclinic Coriolis iterations.
-A Higdon unsplit version will be provided that is identical to Higdon split but where the full velocity is solved for in the baroclinic stage, and nothing is done in the barotropic stage.
+Input variables should be provided for: timestepping type, number of split explicit iterations, number of barotropic subcycles, and number of baroclinic Coriolis iterations.
+A unsplit version will be provided that is identical to split explicit but where the full velocity is solved for in the baroclinic stage, and nothing is done in the barotropic stage.
 
 
 
@@ -74,9 +77,9 @@
 
 \chapter{Algorithmic Formulations}
 
-%\section{Design Solution: Implement Higdon split time-stepping in z-level mpas}
+%\section{Design Solution: Implement split explicit time-stepping in z-level mpas}
 %Date last modified: 2011/05/4 \\
-%Contributors: Mark, Chris, Todd \\
+%Contributors: Mark, Todd \\
 
 \section{MPAS-Ocean time splitting, z-level}
 The MPAS-ocean z-level formulation solves the following equations for thickness, momentum, and tracers at layer $k$:
@@ -212,19 +215,22 @@
 for the split system.  For the unsplit algorithm, the barotopric equations (\ref{u btr 3}) and (\ref{h btr 3}) are not used, so (\ref{h2}) is used to find the top layer thickness.
 
 </font>
<font color="gray">ewpage
-\section{Higdon time splitting algorithm, z-level}
+\section{Split explicit time stepping algorithm, z-level}
+
 {\bf Prepare variables before first iteration}\\
 Always use most recent available for forcing terms.  The first time, use end of last timestep.
 \begin{eqnarray} &amp;&amp;
 {\bf u}_k^*={\bf u}_{k,n},\;\;
 w_k^*=w_{k,n},\;\;
 p_k^*=p_{k,n}, \;\;
-{\varphi}_k^* = {\varphi}_{k,n},\;\;
-\zeta_{n+1}^{edge} = \zeta_n^{edge} \\&amp;&amp;
+{\varphi}_k^* = {\varphi}_{k,n},\\&amp;&amp;
+h_{k,*} = h_{k,n}, \;\; h_{k,*}^{edge} = h_{k,n}^{edge},
+\zeta_{*} = \zeta_n,\;\;
+{\bf u}_{k,*}^{bolus} = {\bf u}_{k,n}^{bolus} \\&amp;&amp;
 \ubtr_n = \left. \textstyle \sum_{k=1}^{N^{edge}}
-\left(\zeta_{k,n}^{edge}+\Delta z_k\right)  {\bf u}_{k,n}
+h_{k,n}^{edge}  {\bf u}_{k,n}
 \right/\textstyle 
-\sum_{k=1}^{N^{edge}}\left(\zeta_{k,n}^{edge}+\Delta z_k\right) 
+\sum_{k=1}^{N^{edge}}h_{k,n}^{edge} 
 \mbox{, on start-up only.}
 \\&amp;&amp;
 \mbox{Otherwise, }\ubtr_n\mbox{ from previous step.}\\&amp;&amp;
@@ -233,26 +239,25 @@
 {\bf u}_{k,n+1/2}'={\bf u}_{k,n}' 
 \end{eqnarray}
 
-Note (\ref{ubcl10}) is required because implicit vertical mixing at the end of the last step may add a barotropic component to ${\bf u}'_{k,n}$.  That component is then to the barotropic forcing $G$ in stange 1.
+The full algorithm, Stages 1--3 are iterated.  This is typically done using two iterations, like a predictor--corrector timestep.  The flag for the number of these large iterations is {\tt config\_n\_ts\_iter}.
 
 {\bf Stage 1: Baroclinic velocity (3D), explicit with long timestep}\\
 Iterate on linear Coriolis term only.
 \begin{eqnarray} &amp;&amp;
-\mbox{compute } {\bf T}({\bf u}_k^*,w_k^*,p_k^*) +g</font>
<font color="red">abla \zeta^* \\ &amp;&amp;
-\mbox{compute weights, }\omega_k  =  
-\left(\zeta_{k,n+1}^{edge}+\Delta z_k\right) 
-\left/\textstyle\sum_{k=1}^{N^{edge}} 
-\left(\zeta_{k,n+1}^{edge} + \Delta z_k\right)
+\mbox{compute } {\bf T}^u({\bf u}_k^*,w_k^*,p_k^*) +g</font>
<font color="red">abla \zeta^* \\ &amp;&amp;
+\mbox{compute weights, }\omega_k  =  h_{k,*}^{edge}
+\left/\textstyle\sum_{k=1}^{N^{edge}} h_{k,*}^{edge}
 \right.\\ &amp;&amp;
 \left\{\begin{array}{l} \label{u bcl s2} 
 \mbox{compute }{\bf u}_{k,n+1/2}^{'\;\perp}\mbox{ from }{\bf u}_{k,n+1/2}'\\ 
 {\tilde {\bf u}}'_{k,n+1} = {\bf u}'_{k,n} + \Delta t 
-\left( -f{\bf u}_{k,n+1/2}^{'\;\perp} + {\bf T}({\bf u}_k^*,w_k^*,p_k^*) 
+\left( -f{\bf u}_{k,n+1/2}^{'\;\perp} + {\bf T}^u({\bf u}_k^*,w_k^*,p_k^*) 
 +g</font>
<font color="gray">abla \zeta^* \right)
 \\  
 {\overline {\bf G}} = 
 \frac{1}{\Delta t}
 \sum_{k=1}^{N^{edge}} \omega_k {\tilde {\bf u}}'_{k,n+1}
+\mbox{ (unsplit: ${\overline {\bf G}} =0$ )}
 \\ 
 {\bf u}'_{k,n+1} = {\tilde {\bf u}}'_{k,n+1} - \Delta t {\overline {\bf G}}
 \\
@@ -261,113 +266,146 @@
 \mbox{boundary update on }{\bf u}'_{k,n+1/2}
 \end{array}
 \right. \;\;\; l=1\ldots L \\ &amp;&amp;
-{\bf u}^{'*}_k = {\bf u}'_{k,n+1/2}
 \end{eqnarray}
 
+The bracketed computation is iterated $L$ times.  The default method is to use $L=1$ on the first time through Stages 1--3, and $L=2$ the second time through.  This is set by the flags {\tt
+config\_n\_bcl\_iter\_beg = 1, config\_n\_bcl\_iter\_mid = 2, config\_n\_bcl\_iter\_end = 2}.
+In the case of two iterations through Stages 1--3, the flag {\tt config\_n\_bcl\_iter\_mid} is not used.
+
+</font>
<font color="red">ewpage
 {\bf Stage 2: Barotropic velocity (2D), explicitly subcycled}\\
-Advance $\ubtr$ and $\zeta$ as a coupled system.  
+Advance $\ubtr$ and $\zeta$ as a coupled system through $2J$ subcycles, ending at time $t+2\Delta t$. 
 \begin{eqnarray}   
 \label{u btr 5} &amp;&amp;
-\left\{\begin{array}{l}
-\mbox{compute }\ubtr_{n+(j-1)/J}^{\perp}\mbox{ from }\ubtr_{n+(j-1)/J} \\ 
+\mbox{{\bf velocity predictor step:}} </font>
<font color="blue">n \\ &amp;&amp;
+\mbox{compute }\ubtr_{n+(j-1)/J}^{\perp}\mbox{ from }\ubtr_{n+(j-1)/J} \\ &amp;&amp; 
 {\tilde \ubtr}_{n+j/J} = \ubtr_{n+(j-1)/J} + \frac{\Delta t}{J} \left(
 - f\ubtr_{n+(j-1)/J}^{\perp}
 - g</font>
<font color="red">abla \zeta_{n+(j-1)/J} + {\overline {\bf G}_j}
-\right)\\
-\mbox{boundary update on }{\tilde \ubtr}_{n+j/J}\\
-\zeta_{n+(j-1)/J}^{edge} = Interp(\zeta_{n+(j-1)/J}) \\
+\right)\\ &amp;&amp;
+\mbox{boundary update on }{\tilde \ubtr}_{n+j/J}\\ &amp;&amp; 
+</font>
<font color="blue">n \\ &amp;&amp;
+\mbox{{\bf SSH predictor step:}}</font>
<font color="red">n \\ &amp;&amp;
+\zeta_{n+(j-1)/J}^{edge} = Interp(\zeta_{n+(j-1)/J}) \\ &amp;&amp;
 \tilde{\bf F}_j = \left((1-\gamma_1){\ubtr}_{n+(j-1)/J}  + \gamma_1{\tilde \ubtr}_{n+j/J}\right)
-\left(\zeta_{n+(j-1)/J}^{edge} + \sum_{k=1}^{N^{edge}} \Delta z_k \right)\\
+\left(\zeta_{n+(j-1)/J}^{edge} + H^{edge} \right)\\ &amp;&amp;
 \tilde\zeta_{n+j/J} = \zeta_{n+(j-1)/J} + \frac{\Delta t}{J} \left( -
- </font>
<font color="red">abla \cdot \tilde{\bf F}_j \right) \\
-\mbox{boundary update on } \tilde\zeta_{n+j/J} \\
-\mbox{compute }{\tilde \ubtr}_{n+j/J}^{\perp}\mbox{ from }{\tilde \ubtr}_{n+j/J} \\ 
+ </font>
<font color="blue">abla \cdot \tilde{\bf F}_j \right) \\ &amp;&amp;
+\mbox{boundary update on } \tilde\zeta_{n+j/J} \\ &amp;&amp; 
+</font>
<font color="blue">n \\ &amp;&amp;
+\mbox{{\bf velocity corrector step:}}</font>
<font color="blue">n \\ &amp;&amp;
+\mbox{compute }{\tilde \ubtr}_{n+j/J}^{\perp}\mbox{ from }{\tilde \ubtr}_{n+j/J} \\ &amp;&amp; 
 \ubtr_{n+j/J} = \ubtr_{n+(j-1)/J} + \frac{\Delta t}{J} \left(
 - f{\tilde \ubtr}_{n+j/J}^{\perp}
 - g</font>
<font color="red">abla \left((1-\gamma_2)\zeta_{n+(j-1)/J} + \gamma_2\tilde\zeta_{n+j/J}\right) + {\overline {\bf G}_j}
-\right),\\
-\mbox{boundary update on }\ubtr_{n+j/J}\\
+\right),\\ &amp;&amp;
+\mbox{boundary update on }\ubtr_{n+j/J}\\ &amp;&amp; 
+</font>
<font color="blue">n\\ &amp;&amp;
+\mbox{{\bf SSH corrector step:}}</font>
<font color="red">n \\ &amp;&amp;
+\tilde\zeta_{n+j/J}^{\;edge}= Interp\left((1-\gamma_2)\zeta_{n+(j-1)/J} + \gamma_2\tilde\zeta_{n+j/J}\right) \\&amp;&amp;
 {\bf F}_j = \left((1-\gamma_3){\ubtr}_{n+(j-1)/J}  + \gamma_3 \ubtr_{n+j/J}\right)
-\left(\zeta_{n+(j-1)/J}^{edge} + \sum_{k=1}^{N^{edge}} \Delta z_k \right)\\
+\left(\tilde\zeta_{n+j/J}^{\;edge} + H^{edge} \right)\\ &amp;&amp;
 \zeta_{n+j/J} = \zeta_{n+(j-1)/J} + \frac{\Delta t}{J} \left( -
- </font>
<font color="red">abla \cdot {\bf F}_j \right) \\
-\mbox{boundary update on } \zeta_{n+j/J} \\
-\end{array}
-\right. \;\;\; j=1\ldots J \hspace{1cm}
+ </font>
<font color="red">abla \cdot {\bf F}_j \right) \\ &amp;&amp;
+\mbox{boundary update on } \zeta_{n+j/J} 
 \end{eqnarray}   
-Note: The last three lines may be iterated to update the velocity in the Coriolis term.  ${\overline {\bf G}}_j$ may vary over the barotropic subcycles, as long as 
-$\frac{1}{J}\sum_{j=1}^{J}{\overline {\bf G}}_j={\overline {\bf G}}$.
+Repeat $j=1\ldots 2J$ to step through the barotropic subcycles.  There are two predictor and two corrector steps for each subcycle.  At the end of $2J$ subcycles, we have progressed through two baroclinic timesteps, i.e.\ through $2\Delta t$.  In the code, {\tt config\_n\_btr\_subcycles} is $J$, while {\tt config\_btr\_subcycle\_loop\_factor=2} is the ``2'' coefficient in $2J$.  
 
-{\bf Stage 2 continued}
+The input flag {\tt config\_btr\_solve\_SSH2=.true.} runs the algorithm as shown, while {\tt .false.} does not include the SSH corrector step.  Here $H^{edge}$ is the total column depth without SSH perturbations, that is, from the higher cell adjoining an edge to $z=0$.
+
+The velocity corrector step may be iterated to update the velocity in the Coriolis term.  The number of iterations is controlled by {\tt config\_n\_btr\_cor\_iter}, and is usually set to two.
+
+</font>
<font color="red">ewpage
+{\bf Stage 2 continued}\\
+
+The coefficients $(\gamma_1,\gamma_2,\gamma_3)$ control weighting between the old and new variables in the predictor velocity, corrector SSH gradient, and corrector velocity, respectively.  These are typically set as $(\gamma_1,\gamma_2,\gamma_3)=(0.5,1,1)$, but this is open to investigation.  These flags are 
+{\tt config\_btr\_gam1\_uWt1, config\_btr\_gam2\_SSHWt1, config\_btr\_gam3\_uWt2}.
+
+The baroclinic forcing ${\overline {\bf G}}_j$ may vary over the barotropic subcycles, as long as 
+$\frac{1}{2J}\sum_{j=1}^{2J}{\overline {\bf G}}_j={\overline {\bf G}}$.  This option is not currently implemented in the code.
+
 \begin{eqnarray} &amp;&amp;  
-\zeta^* = \frac{1}{J+1} \sum_{j=0}^{J} \zeta_{n+j/J}
-\label{ssh5}
+\ubtr_{avg} = \frac{1}{2J+1} \sum_{j=0}^{2J} \ubtr_{n+j/J} 
+\label{ustar5}\\&amp;&amp;
+{\overline {\bf F}} = \frac{1}{2J} \sum_{j=1}^{2J}{\bf F}_j \\&amp;&amp;
+\mbox{boundary update on }{\overline {\bf F}}
 \\&amp;&amp;
-%\zeta^{*\;edge} = Interp(\zeta^*)\\&amp;&amp;
-\ubtr^* = \frac{1}{J+1} \sum_{j=0}^{J} \ubtr_{n+j/J} 
-\label{ustar5}\\&amp;&amp;
-{\overline {\bf F}} = \frac{1}{J} \sum_{j=1}^{J}{\bf F}_j \\&amp;&amp;
-\mbox{boundary update on }{\overline {\bf F}}\\&amp;&amp;
-\mbox{check: } \zeta_{n} + \Delta t \left( - </font>
<font color="red">abla \cdot {\overline {\bf F}} \right)  
-\mbox{ should match } 
-\zeta_{n+1}
-\mbox{ from subcycling.}  \\&amp;&amp;
 {\bf u}^{corr} = \left( {\overline {\bf F}} 
-  - \sum_{k=1}^{N^{edge}} \left(\zeta_{k,n}^{*\;edge}+\Delta z_k\right)  {\bf u}_k^* \right)
-\left/ \sum_{k=1}^{N^{edge}} \left(\zeta_{k,n}^{*\;edge}+\Delta z_k\right)   \right. 
-\mbox{, correction velocity}\hspace{.5cm}\\&amp;&amp;
-{\bf u}^{tr}_k = \ubtr^* + {\bf u}^{'*}_k + {\bf u}^{corr}
-\mbox{, transport velocity for tracer equation}\\&amp;&amp;
-h_{1}^* = \zeta^*+\Delta z_1,\;\;
-h_{1}^{*\;edge} = interp(\zeta^*) + \Delta z_1
+  - \sum_{k=1}^{N^{edge}} h_{k,*}^{edge} 
+   \left( \ubtr_{avg} + {\bf u}'_{k,n+1/2} + {\bf u}_{k,*}^{bolus}\right) \right)
+\left/ \sum_{k=1}^{N^{edge}} h_{k,*}^{edge}   \right. \label{ucorr}
+\\ &amp;&amp;
+{\bf u}^{tr}_k = \ubtr_{avg} + {\bf u}'_{k,n+1/2} + {\bf u}_{k,*}^{bolus} + {\bf u}^{corr} \label{utr}
 \end{eqnarray}
-At the end of stage 2, we have ${\bf u}^*,\ubtr^*,{\bf u}^{'*}_k, \zeta^*$, where the * indicates a variable that is to be used for tendency terms, and is time-averaged, centered at time $n+1/2$.  We also have variables $\ubtr_{n+1}, \zeta_{n+1}$, the instantaneous variables from the end of the barotropic subcycling.
+where ${\bf u}_{k,*}^{bolus}$ is the GM bolus velocity computed at the end of stage 3, and ${\bf u}^{tr}$ is the transport velocity used in the advection terms for for thickness and tracers.
 
-{\bf Stage 3: Tracer, density, pressure, vertical velocity prediction} 
-\begin{eqnarray}   &amp;&amp;
+For unsplit explicit, skip all computations in stage 2.  Instead, set:
+\bea
+&amp;&amp;
+\ubtr_{avg} = 0 \\ &amp;&amp;
+{\bf u}^{tr}_{k} = {\bf u}'_{k,n+1/2} + {\bf u}_{k,*}^{bolus}
+\eea
+
+</font>
<font color="blue">ewpage
+{\bf Stage 3: update thickness, tracers, density and pressure}
+\bea &amp;&amp;
 \label{w3}
 w^{*\; top}_k 
 = w^{*\; top}_{k+1} 
 - </font>
<font color="red">abla \cdot \left( \Delta z_k {\bf u}^{tr}_k \right), \;\;\; k=N\ldots 2.
 \\&amp;&amp;
-\label{phi3}
-{\varphi}_{k,n+1} = \frac{1}{h_{k,n+1}} \left[
-h_{k,n} \varphi_{k,n} 
-+ \Delta t
-\left(  - </font>
<font color="blue">abla \cdot \left( h_k^{*\; edge}\varphi_k^* {\bf u}_k^{tr} \right)  
+T^h_k = \left(  - </font>
<font color="blue">abla \cdot \left( h_k^{*\; edge} {\bf u}_k^{tr} \right)  
+- \frac{\partial}{\partial z} \left( h_k^* w_k^* \right)  \right)
+%-D_k^* + w_{k+1}^* - w_{k}^* 
+\\ &amp;&amp; 
+T^\varphi %({\bf u}_k^*,w_k^*,\varphi_k^*) 
+= \left(  - </font>
<font color="red">abla \cdot \left( h_k^{*\; edge}\varphi_k^* {\bf u}_k^{tr} \right)  
 - \frac{\partial}{\partial z} \left( h_k^*\varphi_k^* w_k^* \right)  
-\right.\right.
-\\&amp;&amp; </font>
<font color="red">onumber \hspace{5cm}
-\left.\left.
 + </font>
<font color="black">abla\cdot\left(h_k^* \kappa_h </font>
<font color="red">abla\varphi_k^* \right)
 + h_k^* \frac{\partial }{\partial z} 
   \left( \kappa_v \frac{\partial \varphi_k^*}{\partial z} \right)
-\right) \right]
-\\&amp;&amp;
-\mbox{boundary update on }{\varphi}_{k,n+1} \\&amp;&amp;
-\varphi_{k}^* = \textstyle\frac{1}{2}\left(\varphi_{k,n} +\varphi_{k,n+1}\right) 
-\mbox{ (not on last iteration)}
-\\&amp;&amp;
-{\bf u}^*_k = \ubtr^* + {\bf u}^{'*}_k 
-\\&amp;&amp; \label{rho3}
+\right) 
+\\ &amp;&amp;
+\mbox{boundary update on tendencies: } T^h,T^\varphi\\&amp;&amp; 
+h_{k,n+1} = h_{k,n} +  \Delta t T^h_k \\&amp;&amp; 
+\label{phi3}
+{\varphi}_{k,n+1} = \frac{1}{h_{k,n+1}} \left[
+h_{k,n} \varphi_{k,n} 
++ \Delta t T^\varphi_k \right]
+\label{phi}
+\eea
+
+{\bf Reset variables}
+\bea 
+\begin{array}{ll}
+\mbox{\bf if iterating} &amp; \mbox{\bf after final iteration}\smallskip\\
+{\bf u}'_{k,*} = {\bf u}'_{k,n+1/2} &amp; {\bf u}'_{k,n+1} \mbox{ from stage 1}\smallskip\\
+\ubtr_{*} = \ubtr_{avg} \mbox{ from stage 2} &amp; \ubtr_{n+1} = \ubtr_{avg} \mbox{ from stage 2}\smallskip\\
+{\bf u}_{k,*} = \ubtr_{*} + {\bf u}'_{k,*} &amp; {\bf u}_{k,n+1} = \ubtr_{n+1} + {\bf u}'_{k,n+1}\smallskip\\
+h_{k,*} = \frac{1}{2}\left(h_{k,n} +h_{k,n+1}\right)  
+ &amp; h_{k,n+1} \mbox{ from stage 3} \smallskip\\
+\varphi_{k,*} = \frac{1}{2}\left(\varphi_{k,n} +\varphi_{k,n+1}\right)
+ &amp; \varphi_{k,n+1} \mbox{ from stage 3} \smallskip\\
+\mbox{diagnostics:}
+&amp; \mbox{diagnostics:} \\
 {\rho}_k^* = EOS(T^*_k, S^*_k) 
-\\&amp;&amp;
+ &amp; {\rho}_{k,n+1} = EOS(T_{k,n+1}, S_{k,n+1}) \smallskip\\
 \label{p3}
-{p}_k^* =  g {\rho}_1^* \left( \zeta^* + \frac{1}{2}\Delta z_1 \right) + 
-  \frac{g}{2}\sum_{l=2}^{k} 
-\left({\rho}_{l-1}^* \Delta z_{l-1}
-+ {\rho}_{l}^* \Delta z_{l} 
-\right)
-\end{eqnarray}
-{\bf Iteration} \\
-If iterating, return to stage 1.  \\
-If complete, we have ${\bf u}'_{k,n+1}$ from (\ref{u bcl s2}), $\ubtr_{n+1}, \zeta_{n+1}$ from the end of barotropic subcycling (\ref{u btr 5}), and $\varphi_{k,n+1}$ from (\ref{phi3}).  These are instantanious variables, not averaged.  Then
-\begin{eqnarray} &amp;&amp;
-{\bf u}_{k,n+1} =\ubtr_{n+1} + {\bf u}'_{k,n+1}
-\\ &amp;&amp;
-h_{1,n+1} = \zeta_{n+1} + \Delta z_1
-\end{eqnarray}
-and compute the full end-of-step diagnostics, including $w^{top}_{k,n+1}$, $\rho_{k,n+1}$ and $p_{k,n+1}$.
+{p}_k^* =  g\sum_{k'=1}^{k-1} 
+{\rho}_{k'}^* h_{k'}
++ \frac{1}{2}g{\rho}_{k}^* h_k
+ &amp; {p}_{k,n+1} =  g\sum_{k'=1}^{k-1} 
+{\rho}_{k',n+1} h_{k',n+1}
++ \frac{1}{2}g{\rho}_{k,n+1}^* h_{k,n+1} \smallskip \\
+h^{edge}_{k,*} = interp(h_{k,*}) &amp; 
+h^{edge}_{k,n+1} = interp(h_{k,n+1}) \smallskip\\
+\zeta_* = \sum_{k=1}^{kmax}h_{k,*} - H &amp;
+\zeta_{n+1} = \sum_{k=1}^{kmax}h_{k,n+1} - H \\
+\mbox{compute }{\bf u}_{k,*}^{bolus} &amp;
+\mbox{compute }{\bf u}_{k,n+1}^{bolus} 
+\end{array}
+\eea
+where H is the total column height with zero sea surface height.
 
 </font>
<font color="gray">ewpage
 {\bf notes:}
@@ -398,7 +436,7 @@
 \left(  - </font>
<font color="red">abla \cdot \left( \sum_{k=1}^N h_k^{*\; edge} {\bf u}_k^{tr} \right)  
 \right)
 \end{eqnarray}
-For consistency between the barotropic and summed baroclinic thickness flux, (or equivalently, to make \ref{h6} identical to $\zeta^* = \zeta_{n} + \Delta t \left( - </font>
<font color="gray">abla \cdot {\overline {\bf F}} \right)$ in (\ref{ssh5})), we must enforce
+For consistency between the barotropic and summed baroclinic thickness flux, we must enforce
 \begin{eqnarray}
 {\overline {\bf F}} &amp;=&amp; \sum_{k=1}^N h_k^{*\; edge} {\bf u}_k^{tr}.
 \label{BclBtrFluxConsistency}
@@ -412,7 +450,7 @@
 \end{eqnarray}
 
 </font>
<font color="gray">ewpage
-\section{Higdon unsplit algorithm, z-level}
+\section{Unsplit algorithm, z-level}
 {\bf Prep variables before first iteration}\\
 Always use most recent available for forcing terms.  The first time, use end of last timestep.
 \begin{eqnarray} &amp;&amp;
@@ -557,21 +595,21 @@
 
 \chapter{Design and Implementation}
 
-\section{Design Solution: Higdon split time-stepping in z-level mpas}
+\section{Design Solution: split explicit time-stepping in z-level mpas}
 Date last modified: 2011/05/4 \\
-Contributors: Mark, Chris, Todd \\
+Contributors: Mark, Todd \\
 
 \section{New variables}
 
 \begin{verbatim}
 namelist character timestep        config_time_integration 'RK4',
- also: 'higdon_split','higdon_unsplit'
-namelist integer   timestep_higdon config_n_ts_iter     2
-namelist integer   timestep_higdon config_n_bcl_iter_beg   4
-namelist integer   timestep_higdon config_n_bcl_iter_mid   4
-namelist integer   timestep_higdon config_n_bcl_iter_end   4
-namelist integer   timestep_higdon config_n_btr_subcycles  10
-namelist logical   timestep_higdon config_compute_tr_midstage true
+ also: 'unsplit','split_explicit'
+namelist integer   split_explicit_ts config_n_ts_iter     2
+namelist integer   split_explicit_ts config_n_bcl_iter_beg   4
+namelist integer   split_explicit_ts config_n_bcl_iter_mid   4
+namelist integer   split_explicit_ts config_n_bcl_iter_end   4
+namelist integer   split_explicit_ts config_n_btr_subcycles  10
+namelist logical   split_explicit_ts config_compute_tr_midstage true
 \end{verbatim}
 
 
@@ -805,11 +843,11 @@
 
 \chapter{Testing}
 
-\section{Testing and Validation: Higdon split time-stepping in z-level mpas}
+\section{Testing and Validation: split explicit time-stepping in z-level mpas}
 Date last modified: 2011/05/04 \\
 Contributors: (add your name to this list if it does not appear) \\
 
-Testing: Compare Runge-Kutta versus Higdon unsplit and Higdon split.
+Testing: Compare Runge-Kutta versus unsplit and split explicit.
 
 %-----------------------------------------------------------------------
 

</font>
</pre>