## Numerical Issues

- arango
- Site Admin
**Posts:**1091**Joined:**Wed Feb 26, 2003 4:41 pm**Location:**IMCS, Rutgers University-
**Contact:**

### Numerical Issues

Use this thread to discuss any numerical issues associated with the documentation in the ROMS/TOMS manual.

I was wondering if someone could give me hints about where to find documentation of the 3d momentum time-stepping (Adams-Bashforth 3) used in version 2.2. I tried to understand the Shchepetkin-McWilliams paper, and could extract some info about the barotropic LF-AM3, but I think they don't cover the AB3, do they? Is there such a study about the AB3, directed to ocean modelling?

Thanks,

--Stefan

stepping procedure used in ROMS v. 1.8....2.2 in the Shchepetkin

McWilliams,2005 ROMS paper because this algorithm was always considered

as a provisional variant to be replaced with a more refined one.

The main 3D procedure of v. 1.8....2.2 can be classified as an AB3-TR

Generalized forward-backward scheme [AB3 = Adams--Bashforth 3rd order

(open parabolic rule); TR = trapezoidal rule) and is summarized as

follows:

Code: Select all

```
u^{n+1} = u^n +dt*[(23/12)*RHSU^n - (4/3)*RHSU^{n-1} + (5/12)*RHSU^{n-2}]
T^{n+1} = T^n -dt*div[ (1/2)*(U^{n+1}+U^n)*T^{n+1/2} ]
```

interpolation coefficients are recognized as the ones for Adams--Moulton

3rd-order closed parabolic integration rule). In its turn T^{n+1,*} comes

from LF predictor sub-step,

Code: Select all

```
T^{n+1,*}=T^{n-1} -2*dt*div[ U^n * T^n ]
```

step with artificial continuity equation, and, The LF-step is also

combined with AM3 interpolation (this is done in pre_step.F), so that

you actually do not see a naked LF explicitly in the code, but rather

you see

Code: Select all

```
T^{n+1/2} = (5/12)*[ T^{n-1} -2*dt*div(U^n * T^n)]
+ (2/3)*T^n -(1/12)*T^{n-1}
```

Code: Select all

```
T^{n+1/2} = (1/3)*T^{n-1} + (2/3)*T^n -(5/6)*dt*div(U^n * T^n)
```

(5/6) as "1-Gamma"; and Gamma is set to 1/6.

Above for simplicity I have omitted multiplication and division by Hz

and metric factors. Lowercase "u" means velocity; uppercase "U" fluxes

(which map onto Huon and Hvom in the code.

Setting of fluxes for corrector sub-step of tracer equations via

Trapezoidal Rule (TR) occurs in "step3d_uv.F", see

Code: Select all

```
Huon(i,j,k)=0.5_r8*(Huon(i,j,k)+u(i,j,k,nnew)*DC(i,k))
```

(Hz^{n+1}/pn or pm for v-component), and Huon(i,j,k) is previously

computed flux "U^n" (in the actual code Huon(i,j,k) is initially computed

within "set_massflux.F" (in ROMS 2.2; I believe that older codes may have

different file name), and then modified in "step3d_uv.F".

Assuming that the stability limit is governed primarily by the phase speed

of the fastest internal mode (basically the first baroclinic mode), the

details of LF-AM3 stepping for tracer become unimportant because the rate

of change of T is dominated by vertical gradient of T multiplied by

vertical velocity. In this case the whole algorithm maps onto Eq. (2.49)

for Shchepetkin McWilliams, 2005 with beta=5/12, delta=1/2,

gamma=epsilon=0. Its stability limit alpha_max=1.1441551 and its

placement of characteristic roots relatively to unit circle looks similar

to one of Fig. 13, upper left panel (that panel shows AB3-AM3, rather

than AB3-TR version).

The exact figure showing AB3-TR characteristic roots can be seen on page 22 in

http://marine.rutgers.edu/po/Workshops/ ... petkin.pdf

which is presentation on Venice, 2004 workshop.

Thanks for the reply!

I am also searching for information on ROMS finite volume formulations on the c-grid. Specifically, I have troubles understanding how the variables on the boundary of the gridcells are reconstructed to do the integration of the fluxes. It seems like there are many ways to do the reconstruction, and after reading the code and seeing how it is done, I would like to understand why it's done that way. In case anybody is aware of papers on this topic (finite-volumes on a c-grid), it would help me a lot if you could just post the titles here.

Thanks very much,

--Stefan

Sorry, I know this might not be ineresting for most users, but I have to write something into my work, so I would appreciate if somebody could tell me whether the term "finite volume" or "finite differences" ist the more correct one.

Suggestions for literature are welcome!

Thanks,

--Stefan

The main thing is interpretation of the gridded data: "finite-volume" assumes that u(i,j,k) is AVERAGE of field "u" over control volume dV(i,j,k), while finite difference assumes that it is instantaneous value at the location x(i,j,k). This leads to different formulas for things like computing derivatives, although within the second-order accuracy the formulas are mostly agree.

The method of derivation is very different, however. Formally, to derive a finite-volume discrete equation you must integrate you equations of motion over a control volume and, whenever possibly, transform the integral into flux form. Then, technically, it boils down to a three stage procedure: (1) given set of grid-box averages find instantaneous values at grid-box interfaces; (2) compute fluxes (perhaps using nonlinear formulas); and (3) add all fluxes coming in/out each control volume to find change of the amount of quantity there.

Methodologically, operation (3) is always "exact".

Operation (2) is "exact" if flux formulas are linear, which is true only in simplest cases, but generally it requires some approximations to be applied. For example, formally speaking if you have velocity and tracer concentration fields at the interface, then you must multiply them and integrate the product over the interface to get flux. But instead you approximate it by a product of "mean" values of each field averaged over the interface individually, and multiply it by the area of the interface, if you are on a curvilinear grid and areas are changing from one grid element to another. These are not the same, but approximately close to.

Operation (1) is never exact, and is, if fact, root of most errors.

Consider two formulas

a(i+1/2) = -1/16 *a(i-1) + 9/16 *a(i) + 9/16 *a(i+1) -1/16 *a(i+1)

and

a(i+1/2) = -1/12 *a(i-1) + 7/12 *a(i) + 7/12 *a(i+1) -1/12 *a(i+1)

the first one interpolates the value of quantity "a" given at discrete locations x(i), x(i+1), ... to the midpoint x(i+1/2) half-way between x(i) and x(i+1). That is, given instantaneous values, compute instantaneous value.

The second formula assumes that a(i) is the average of the field within the interval x(i-1/2) < x < x(i+1/2), and similarly interprets, a(i+1), a(i+2), and a(i-1). Then it not just interpolates to mid-location, but it performs TRANSLATION from averaged to instantaneous values. As the result, coefficients are different.

What happens if you subtract [a(i+1/2) -a(i-1/2)]/dx ?

If a(i+1/2) is computed using bottom formula, you get the fourth-order accurate formula for approximation for the first derivative. If, on the other hand, you use upper formula, you would not get fourth order accuracy.

If one wants to guarantee things like integral conservation of something, and the grid is curvilinear, then "finite volume" is basically the only way to go.

A good reference to read is

Arakawa and Lamb, 1977 Computational Design of the Basic Dynamical Processes of the UCLA General Circulation Model, Meth. Comput. Phys., vol17, 174-267. Academic Press.

Common usage of terms "finite volume" and "finite difference" is loose and often interchangeable, especially if people are talking about structured grids. SCRUM is a finite volume model as well, and so is POM, but nobody emphasizes that.

In contrast, unstructured grid people, use the term "finite volume" as opposite to "finite element" and tend to make very big deal out of it.

Term "control volume method" is interchangeable with "finite volume".