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.
