Time-averaging of barotropic fields and Power Law Filter

[shchepet]

It is known that the power function weights for time filtering of the vertically integrated fields (POWER_LAW CPP option) result in unstable code when the splitting ratio (ndtfast) becomes too small. For a natural setting of ndtfast=60 to ndtfast=70(most of the realistic basin-scale problems) works just fine, but if one wants to use ndtfast=20, the power function weights become unstable. Normally users tend to revert back to the squared-cosine (COSINE2 CPP option) or even flat-shaped filters in this situation.

The problem is due to restart of the barotropic mode at every 3D time step: it causes weak instability and weights should provide enough dissipation to suppress it. The effect is of 1/ndtfast --> 0 order, when ndtfast becomes large.

Consequently, I introduced a 1-line change into set_weights.F, to compensate for it. After that I can run applications with splitting ratio as low as ndtfast=12 ... 15, which is more robust in this sense than even Hamming window weights. I highly recommend to everybody to get the updated set_weights.F file.

Sasha

[arango]

For those of you using ROMS 2.1, please download the following tar file containing an update to mod_scalars.F and set_weights.F:

http://www.myroms.org/links/weights.tar.gz

For the few of you testing ROMS 2.2 beta, please contact me directly since mod_scalars.F has several additional changes.

For those of you using UCLA version of ROMS, please download the following tar containing an update to set_weights.F:

http://www.myroms.org/links/weights_ucla.tar.gz

We advice you to update these routines. This change will be available in the next release of ROMS.

Good luck

[m.hadfield]

For those of you using ROMS 2.1, please download the following tar file containing an update to mod_scalars.F and set_weights.F:

I suggest in addition the following modification to mod_scalars.F: after the following

Code: Select all

      allocate ( Tobc_in(MT,Ngrids,4) )
      allocate ( Tobc_out(MT,Ngrids,4) )

add

Code: Select all

      Tobc_in = IniVal
      Tobc_out = IniVal

Why? To bring the model down when Tobc_in and Tobc_out are used without being initialised. (Some compilers have an option to prevent the use of uninitialised variables, but in my experience this can't always be relied on.) I've made this change in my own code. I can't remember exactly why, but I'm sure I must have had a good reason!

[arango]

Please notice that:

1) We are trying to provide a foolproof code for those users that do not understand well the time split-explicit algorithm of ROMS. It requires solving the vertically integrated equations with smaller time-step to resolve fast gravity-wave dynamics. Recall that this time step (dtfast) is limited by dx/SQRT(g*h). In ROMS, dtfast=dt/ndtfast where dt is the time-step for the slow, Rossby-like dynamics. We are addressing above the issue of ndtfast being too small. We have mentioned in the past that ROMS is designed to have a big value for ndtfast and not a small one, say less than 15. What to use for ndtfast depends on the application and we expect the user to apply his/her knowledge to determine its appropriate value.

2) There are few papers in the literature that address the issue of stability in time splitting algorithms used in ocean modeling and that is what is implied above by the It is know .... quote above. Users are welcome to examine the literature on this subject.

[shchepet]

Dear All,

First remark is that this set_weight modification is not a bug fix, but it would better fall into category of optimization of code usage. The purpose of filter is to remove fast frequences of the barotropic motions to prevent their aliasing and instability of the model. The downside of it is that it also introduces numerical inaccuracy, at least in formal sense (i.e., formal numerical order of accuracy of time discretization is just first). The overall methodology of designing filter shapes is described in the ROMS paper, which was hanging around for several years and now (at last) has been finally published. That includes the description of S-shaped filters (a.k.a POWER_LAW).

A publishable paper cannot describe each and every detail because of volume constaints (good luck to push it through reviewers, and besides stylistically it would loose focus on what is actually new, if one attempts to do that).

Besides this paper, the theory of digital filtering is quite well understood in general literature, starting with Richard Hamming bell-shaped filter back in 195x.

The original recommended parameter setting in set_weigts (p=2, q=4, r=0.25) is already biased slightly away from the ideal value of r=0.284. It is true that the bias was never fully explained in literature, but the implication is quite clear: to privide a residual dissipation in second-order in order to ensure code stability. In this sense, the magic number of r=r(p,q) which eliminates the leading second-order dissipative truncation term is to be understood as maximum bound beyond which the code becomes unstable, even if 2D time stepping is exact. The value of 0.25 came from experience with regional to basin-scale configurations which are characterised by 4km depth resulting in splitting ratio of 60 to 80: the choice of 0.25 is accurate, and yet provides sufficient margine of safety (in practice ndtfast can be reduced to 30 by reducing 3D time step, and the code still runs). But it is not appropriate for those who primarely interested in estuaries and model setups with maximum depth of 100 meters or so: their splitting ratio tend to be 20 or less. Nobody on the West Coasts [that includes US West Cost (UCLA) , as West Costs of France, Portugal, Pery, Chili, North and South Africa (Bangella)] ever bothered with model setups shallow enough so that natural splitting ratio become too small. So that the matter was naturally deferred.

Meanwhile ROMS community "cooking experience" is also well established: FLAT weigths settings yields the most robust code, then goes hamming window, then S-shaped weights. Accuracy rating goes other way around.

The goal of the last change is to unify, and ultimately construct an algorithm which optimally works for both North Pacific model (splitting ratio ndtfast=78), for upwelling example (ndtfast=20), and for gravitational adjustment test [splitting ratio is just 7 (!)].

Regarding policy stated in globaldefs.h [UCLA equivalent is set_global_definitions.h]. Very early in ROMS evolution in became clear that there are two kind of CPP switches in ROMS:

--- those which define physical configuration, and

--- all others: machine/compiler dependent switches, those, which facilitate code evolutiuon; unsettled numerical issues, to control various model infrastrictures, etc...

It is natural to separate them. The ultimate goal is that switches defined in globaldefs.h would not need to be changed by a general user. In practice is does not always work this way because ROMS model covers very wide range of applications, so algorithms need to be changed sometimes.

[debcox]

I am using ROMS to model shallow water (generally < 30-50m) island and headland wakes, with both idealised and real bathymetry. Most of my work to date has been for unstratified flows. The resulting flows are highly three-dimensional with significant upwelling and downwelling as a result of interactions with the topography. I guess this is a bit different to the type of larger scale flows that ROMS was originally designed for.

I find that the model generally crashes if I set NDTFAST > 1, even using the COSINE2 filter. Is this what would be expected in this type of situation? It would be great to get the model running faster if possible.

Deborah