Home Contact Us Site Map  
 
       
    next up previous contents
Next: 2.9.1 Crank-Nickelson barotropic time Up: 2. Discretization and Algorithm Previous: 2.8 Non-hydrostatic formulation   Contents


2.9 Variants on the Free Surface

We now describe the various formulations of the free-surface that include non-linear forms, implicit in time using Crank-Nicholson, explicit and [one day] split-explicit. First, we'll reiterate the underlying algorithm but this time using the notation consistent with the more general vertical coordinate $ r$. The elliptic equation for free-surface coordinate (units of $ r$), corresponding to 2.16, and assuming no non-hydrostatic effects ( $ \epsilon _{nh}=0$) is:

$\displaystyle \epsilon_{fs} {\eta}^{n+1} - {\bf\nabla}_h \cdot \Delta t^2 (R_o-R_{fixed}) {\bf\nabla}_h b_s {\eta}^{n+1} = {\eta}^*$     (2.72)

where
$\displaystyle {\eta}^* = \epsilon_{fs} \: {\eta}^{n} - \Delta t {\bf\nabla}_h \... ...int_{R_{fixed}}^{R_o} \vec{\bf v}^* dr \: + \: \epsilon_{fw} \Delta_t (P-E)^{n}$     (2.73)

\fbox{ \begin{minipage}{4.75in} {\em S/R SOLVE\_FOR\_PRESSURE} ({\em solve\_for\... ..._PRESSURE.h) \par $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h) \par \end{minipage} }

Once $ {\eta}^{n+1}$ has been found, substituting into [*] yields $ \vec{\bf v}^{n+1}$ if the model is hydrostatic ( $ \epsilon _{nh}=0$):

$\displaystyle \vec{\bf v}^{n+1} = \vec{\bf v}^{*} - \Delta t {\bf\nabla}_h b_s {\eta}^{n+1} $

This is known as the correction step. However, when the model is non-hydrostatic ( $ \epsilon_{nh}=1$) we need an additional step and an additional equation for $ \phi'_{nh}$. This is obtained by substituting 2.54, 2.55 and 2.56 into continuity:

$\displaystyle \left[ {\bf\nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}... ...a t} \left( {\bf\nabla}_h \cdot \vec{\bf v}^{**} + \partial_r \dot{r}^* \right)$ (2.74)

where

$\displaystyle \vec{\bf v}^{**} = \vec{\bf v}^* - \Delta t {\bf\nabla}_h b_s {\eta}^{n+1} $

Note that $ \eta^{n+1}$ is also used to update the second RHS term $ \partial_r \dot{r}^* $ since the vertical velocity at the surface ( $ \dot{r}_{surf}$) is evaluated as $ (\eta^{n+1} - \eta^n) / \Delta t$.

Finally, the horizontal velocities at the new time level are found by:

$\displaystyle \vec{\bf v}^{n+1} = \vec{\bf v}^{**} - \epsilon_{nh} \Delta t {\bf\nabla}_h {\phi'_{nh}}^{n+1}$ (2.75)

and the vertical velocity is found by integrating the continuity equation vertically. Note that, for the convenience of the restart procedure, the vertical integration of the continuity equation has been moved to the beginning of the time step (instead of at the end), without any consequence on the solution.

\fbox{ \begin{minipage}{4.75in} {\em S/R CORRECTION\_STEP} ({\em correction\_ste... ...m DYNVARS.h}) \par $v^{n+1}$: {\bf vVel} ({\em DYNVARS.h}) \par \end{minipage} }

Regarding the implementation of the surface pressure solver, all computation are done within the routine SOLVE_FOR_PRESSURE and its dependent calls. The standard method to solve the 2D elliptic problem (2.74) uses the conjugate gradient method (routine CG2D); the solver matrix and conjugate gradient operator are only function of the discretized domain and are therefore evaluated separately, before the time iteration loop, within INI_CG2D. The computation of the RHS $ \eta^*$ is partly done in CALC_DIV_GHAT and in SOLVE_FOR_PRESSURE.

The same method is applied for the non hydrostatic part, using a conjugate gradient 3D solver (CG3D) that is initialized in INI_CG3D. The RHS terms of 2D and 3D problems are computed together at the same point in the code.


Subsections
next up previous contents
Next: 2.9.1 Crank-Nickelson barotropic time Up: 2. Discretization and Algorithm Previous: 2.8 Non-hydrostatic formulation   Contents
mitgcm-support@dev.mitgcm.org
Copyright © 2002 Massachusetts Institute of Technology