<p><b>mpetersen@lanl.gov</b> 2011-03-24 12:11:52 -0600 (Thu, 24 Mar 2011)</p><p>BRANCH COMMIT: Update Requirements and Design document for implicit vertical mixing.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/doc/implicit_vert_mix/implicit_vert_diff_design.pdf
===================================================================
(Binary files differ)
Modified: branches/ocean_projects/doc/implicit_vert_mix/implicit_vert_diff_design.tex
===================================================================
--- branches/ocean_projects/doc/implicit_vert_mix/implicit_vert_diff_design.tex 2011-03-24 18:01:46 UTC (rev 765)
+++ branches/ocean_projects/doc/implicit_vert_mix/implicit_vert_diff_design.tex 2011-03-24 18:11:52 UTC (rev 766)
@@ -97,9 +97,15 @@
The functional forms for vertical viscosity and diffusivity at each layer interface will be identical to POP:
\begin{eqnarray} \label{visc1}
</font>
<font color="black">u_v &=& </font>
<font color="red">u_{bkrd} + Rich\_mix/(1+5Ri)^2\\
-\kappa_v &=& \kappa_{bkrd} + (</font>
<font color="blue">u_{bkrd} + Rich\_mix/(1+5Ri)^2)/(1+5Ri)^2
+\kappa_v &=& \kappa_{bkrd} + (</font>
<font color="blue">u_{bkrd} + Rich\_mix/(1+5Ri)^2)/(1+5Ri)
\end{eqnarray}
+for $Ri>=0$. For unstable stratification, $Ri<0$ and the viscosity and diffusion are set to be very high, based on the vertical diffusive CFL condition,
+\begin{eqnarray} \label{visc1}
+</font>
<font color="blue">u_v =
+\kappa_v = \frac{1}{4} \frac{dz^2}{dt}
+\end{eqnarray}
+
</font>
<font color="gray">ewpage
\section{Design Solution: operator splitting and implicit solve}
@@ -162,17 +168,17 @@
\begin{eqnarray}
\label{tridiag}
&& A_{k-1}\varphi^{n+1}_{k-1} + C_k\varphi^{n+1}_k + A_k\varphi^{n+1}_{k+1}
-= \widetilde{(h\varphi)}^{n+1}_k, \;\; k=1\ldots N\\
-&& A_k = -\frac{2 \Delta t \lambda \kappa^{top}_{k+1}}{h_k + h_{k+1}}, \;\; k=1\ldots N-1 \\
+=\frac{\widetilde{(h\varphi)}^{n+1}}{h^{n+1}}, \;\; k=1\ldots N\\
+&& A_k = -\frac{2 \Delta t \lambda \kappa^{top}_{k+1}}{h_k^{n+1}(h_k^{n+1} + h_{k+1}^{n+1})}, \;\; k=1\ldots N-1 \\
&& A_0=A_N=0, \\
-&& C_k = h_k^{n+1} - A_k - A_{k-1}
+&& C_k = 1 - A_k - A_{k-1}
\end{eqnarray}
where $\kappa^{top}_{k+1}$ is the tracer diffusion on the interface between layers $k$ and $k+1$, and the bottom cell is $N$. This equation set is nearly identical to that in the POP reference manual, section 4.2.3. The boundary condition for the implicit tracer solve is no flux,
\begin{eqnarray}
\label{tridiag bc}
\frac{\partial \varphi}{\partial z} = 0
\end{eqnarray}
-at the top of cell 1 and the bottom of cell N (top of cell $N+1$). When this is implemented in the derivation from (\ref{tridiag}) to (\ref{tridiag bc}), the result is $A_0=A_N=0$ above. T
+at the top of cell 1 and the bottom of cell N (top of cell $N+1$). When this is implemented in the derivation from (\ref{tridiag}) to (\ref{tridiag bc}), the result is $A_0=A_N=0$ above.
We now proceed to operator splitting for the momentum equation. The continuous equation is
\begin{equation}
@@ -213,11 +219,11 @@
where $\delta_z$ is a first order finite difference. The tridiagonal solve is then
\begin{eqnarray}
\label{momtridiag}
-&& A_{k-1}\varphi^{n+1}_{k-1} + C_k\varphi^{n+1}_k + A_k\varphi^{n+1}_{k+1}
+&& A_{k-1} u^{n+1}_{k-1} + C_k u^{n+1}_k + A_k u^{n+1}_{k+1}
= \tilde{u}^{n+1}_k, \;\; k=1\ldots N\\
-&& A_k = -\frac{2 \Delta t \lambda </font>
<font color="blue">u^{top}_{k+1}}{h_k + h_{k+1}}, \;\; k=1\ldots N-1 \\
+&& A_k = -\frac{2 \Delta t \lambda </font>
<font color="gray">u^{top}_{k+1}}{h_k^{n+1}(h_k^{n+1} + h_{k+1}^{n+1})}, \;\; k=1\ldots N-1 \\
&& A_0=A_N=0, \\
-&& C_k = h_k^{n+1} - A_k - A_{k-1}
+&& C_k = 1 - A_k - A_{k-1}
\end{eqnarray}
Again, no flux boundary conditions $\partial u/\partial z=0$ imply that $A_0=A_N=0$.
@@ -252,10 +258,11 @@
Right now, the vertical viscosity variable \verb=vertViscTop(k)= is computed columnwise in subroutine \verb=compute_tend=, and the vertical diffusivity variable \verb=vertDiffTop(k)= is computed columnwise in subroutine \verb=compute_scalar_tend=. We propose to change these to 3D variables computed in \verb=compute_solve_diagnostics=. The revised registry would include:
\begin{verbatim}
-var persistent real RiTopCell ( nVertLevelsP1 nCells Time ) 1 o diagnostics
-var persistent real RiTopEdge ( nVertLevelsP1 nEdges Time ) 1 o diagnostics
-var persistent real vertViscTopEdge ( nVertLevelsP1 nEdges Time ) 1 o diagnostics
-var persistent real vertDiffTopCell ( nVertLevelsP1 nCells Time ) 1 o diagnostics
+var persistent real RiTopOfCell ( nVertLevelsP1 nCells Time ) 1 o diagnostics
+var persistent real RiTopOfEdge ( nVertLevelsP1 nEdges Time ) 1 o diagnostics
+var persistent real vertViscTopOfEdge ( nVertLevelsP1 nEdges Time ) 1 o diagnostics
+var persistent real vertDiffTopOfCell ( nVertLevelsP1 nCells Time ) 1 o diagnostics
+var persistent real rhoDisplaced ( nVertLevels nCells Time ) 2 o state
\end{verbatim}
A note will be added that in the future, arrays could be saved by computing these values columnwise on the fly. We feel that the full 3D allocation allows for clear and consolidated code, so is worthwhile in the prototype version. A new variable class, named diagnostics, with only one time level is used for these variables, in order to avoid the memory required for two time-level state class variables.
@@ -264,14 +271,15 @@
The following steps will be used to ensure that variables are computed at the correct locations:
\begin{enumerate}
-\item drhoTopCell(k) = $\rho^*_{k-1}-\rho^*_k$
-\item interpolate drhoTopCell to drhoTopEdge
-\item duEdgeTop(k) = $u_{k-1}-u_k$
-\item interpolate duTopEdge to duTopCell
-\item compute RiTopCell using drhoTopCell and duTopCell
-\item compute vertDiffTopCell from RiTopCell
-\item compute RiTopEdge using drhoTopEdge and duTopEdge
-\item compute vertViscTopEdge from RiTopEdge
+\item compute density of parcel displaced to surface, $\rho^*$
+\item drhoTopOfCell(k) = $\rho^*_{k-1}-\rho^*_k$
+\item interpolate drhoTopOfCell to drhoTopOfEdge
+\item duTopOfEdge(k) = $u_{k-1}-u_k$
+\item interpolate duTopOfEdge to duTopOfCell
+\item compute RiTopOfCell using drhoTopOfCell and duTopOfCell
+\item compute vertDiffTopOfCell from RiTopOfCell
+\item compute RiTopOfEdge using drhoTopOfEdge and duTopOfEdge
+\item compute vertViscTopOfEdge from RiTopOfEdge
\end{enumerate}
@@ -293,8 +301,6 @@
block => domain % blocklist
do while (associated(block))
- call compute_solve_diagnostics(dt, state, mesh)
-
if (config_implicit_vertical_mix) then
call compute_vertical_mix_coeff(dt, state, mesh, diag)
</font>
</pre>