Some time ago I found that Computation of Courant number for the purpose of

diagnostics in done by averaging velocities toward the tracer point (see around

line 220 of diag.F),

Code: Select all

```
my_Cu=0.5_r8*ABS(u(i,j,k,nstp)+u(i+1,j,k,nstp))* &
& dt(ng)*pm(i,j)
my_Cv=0.5_r8*ABS(v(i,j,k,nstp)+v(i,j+1,k,nstp))* &
& dt(ng)*pn(i,j)
my_Cw=0.5_r8*ABS(wvel(i,j,k-1)+wvel(i,j,k))* &
& dt(ng)/Hz(i,j,k)
my_C=my_Cu+my_Cv+my_Cw
IF (my_C.gt.my_max_C) THEN
my_max_C =my_C
my_max_Cu=my_Cu
my_max_Cv=my_Cv
my_max_Cw=my_Cw
my_max_Ci=i
my_max_Cj=j
my_max_Ck=k
END IF
```

Courant number consistently with finite-volume interpretation of grid-point

discrete data.

In the code segment below it advective Courant number is defined as the sum

of fluxes directed OUTWARD from grid box Hz divided by its control volume.

Essentially it is the fraction of water replaced within the grid box "Hz" during

one time step.

Why OUTWARD? Because under this definition, Cu=1 is the condition when

flux-split first-order upstream advection scheme looses its positive definiteness

property (hence stability). Of course, other combinations of spatial discrete

advection schemes and time stepping algorithm should become unstable when

Cu exceeds a different critical value, but it is also comparable to 1.

Note that in the case of "Hz"s fixed in time it does not matter whether it is

incoming or outgoing fluxes, because of the sum of incoming is equal to the

sum of outgoing. If "Hz"s are allowed to change, then it matters and it should

be the sum of outgoing. In principle, this definition is applicable to a code

which allows "drying" cells (i.e., one cannot take out more volume from an

about-to-dry cell that this cell has at the beginning of the time step), or,

in fact, any weird/unstructured grids with not necessarily six facets for

each grid box -- it can be any.

This also solves the dilemma which numbers to report: separate in each

direction, Cx,Cy,Cz or some combination of them, Cx+Cy+Cz. As the practical

utility of computing advective Courant number is to diagnose "hot spots" of

ROMS grids, which cause model blow up, and identify the primary reason -- either

exceeding horizontal or vertical velocity (say due to combination of a large

topographic slope and very fine local vertical resolution) , so it makes sense

to separate the two, but reporting two horizontal components is hardly useful:

ROMS grid resolutions are more or less horizontally isotropic most of the time.

Reported values are: "Cu_Adv" is full (tri-dimensional) Courant number;

"i,j,k" are coordinates where its maximum occurs; "Cu_W" is contribution

into "Cu_Adv" due to vertical velocity.

Code: Select all

```
# ifdef MAX_ADV_CFL
# ifdef MASKING
ciV=dt*rmask(i,j)*pm(i,j)*pn(i,j)/Hz(i,j,k)
# else
ciV=dt*pm(i,j)*pn(i,j)/Hz(i,j,k)
# endif
cw=ciV*( max(W(i,j,k), 0.) - min(W(i,j,k-1), 0.))
cx=cw+ciV*( max(FlxU(i+1,j,k), 0.) -min(FlxU(i,j,k), 0.)
& +max(FlxV(i,j+1,k), 0.) -min(FlxV(i,j,k), 0.))
if (cx .gt. my_Cu_Adv) then
my_Cu_Adv=cx
my_Cu_W=cw
my_i_cmax=i
my_j_cmax=j
my_k_cmax=k
endif
# else
...........
```

In theory, the most dramatic difference between the two diagnostics is when

velocity field is a tri-dimensional checkerboard pattern. Then averaging of

velocity components leads to zero, hence zero Courant number(?), while

computing it via sum of outgoing fluxes may lead to a large number.

Hernan, adapt it from http://www.atmos.ucla.edu/~alex/ROMS/tmp/diag.F.