Hello,
We just fixed a quite important bug (we take it as a bug) in the Roms_Agrif version of the code, that is also present in the Rutgers version. This concerns the computation of the pressure gradient term in the corrector step of the barotropic mode. This computation involves zeta at time n and n-1 which should come from the preceding corrector step, but this is not the case. They actually come via the storage in rubar, rvbar that has been done in the preceding predictor step where the values of zeta have been "filtered". So the algorithm do not match eqs 2.27-2.30 of Sasha and Jim's paper.
More problematic, these "wrong" values of rubar, rvbar computed in the first predictor step are used for the coupling with the 3D part. So that the free surface gradient does not cancel with the one computed in the 3D part. (Additionaly, computation in the first predictor is not centered at time n, since in this first step a trapezoidal scheme is used, this makes the coupling really wrong and the second time step of the barotropic mode is only first order accurate).
A quick fix is to recompute at every predictor step the right value of the pressure gradient term and to use this term both for the coupling and for
the storage in rubar(:,:,krhs). This has been done in the Roms_Agrif version.(you can check the website)
Configurations where a large ndtfast is used will probably be only slightly affected.
We have done several experiments that tend to show that configurations with small ndtfast (ndtfast <20,30) that were blowing up before the corrections are now running ok.
Laurent
2D time stepping and 2D/3D coupling
RE: 2D time stepping and 2D/3D coupling
This observation is generally correct: the timestepping algorithm of barotropic modes of both Agrif and Rutgers code (all versions of both) cannot be accurately mapped into Eqs. (2.27)---(2.30) of the time stepping paper because r.h.s. terms saved from previous time step contain barotropic pressure gradient computed using not zeta[n] and zeta[n-1], but rather their averaged versions (which are averaged during prediction),
0.16*zeta[n-1] + 0.68*zeta[n] + 0.16*zeta[n+1,*]
and
0.16*zeta[n-2] + 0.68*zeta[n-1] + 0.16*zeta[n,*]
(with asterisk (*) means that these values come from predictor sub-step) consequently, computation of barotropic pressure gradient terms during corrector state does not map onto (2.30), but rather onto a version of that equation with zeta[n] and zeta[n-1] replaced as stated above.
This is not a bug, however, but is rather a different algorithm.
In reality, this particular version of barotropic mode was created back in 1998, long before the paper was written (the paper was finalized in early 2003, but it took another two years to publish it).
Back in 1998 the whole apparatus of studying stability of time-stepping algorithms was not developed, and the coefficients 0.16 (see above) and 0.4 (weighting between predicted and corrected values of zeta during corrector stage) was chosen empirically by playing with oscillator system. The 2005 paper actually DERIVES a time stepping algorithm with the desired properties by optimising coefficients of characteristic equations in order to make them produce the desired eigen values. There was no such procedure back in 1998.
HOWEVER, THE PROPOSED FIX DESIGNED TO BRING THIS CODE IN CONSISTENCY WITH THE 2005 PAPER DOES NOT MAKE ANY SENSE WHATSOEVER.
THE WHOLE BAROTROPIC MODE SHOULD BE REPLACED INSTEAD.
The existing step2D codes (both in Rutgers ROMS and Agrif, all versions of both) are now helplessy obsolete.
What is proposed here is to compute barotropic pressure gradient terms twice during predictor substep, which simply adds more computational load into already inefficient code.
A new version of step2D, which is consistent with Eqs. (2.38)---(2.41) of 2005 paper was available in the middle of 2001 (in fact, it broadcasted to ROMS at the end of August 2001, soon after the Boulder meeting) and a newer Forward-Backward (now standard in UCLA codes) scheme followed in 2002. They are approximately 1.5 and 2.5 times more efficient than the original LF-AM3 code. (Due to combined effect of better stability limits and more efficient codes, for example, not storing r.h.s. terms between time steps).
Aparently, both Hernan and Agrif people (with some help from the community) did a very good job in consealing and, in fact, supressing their sole existence.
0.16*zeta[n-1] + 0.68*zeta[n] + 0.16*zeta[n+1,*]
and
0.16*zeta[n-2] + 0.68*zeta[n-1] + 0.16*zeta[n,*]
(with asterisk (*) means that these values come from predictor sub-step) consequently, computation of barotropic pressure gradient terms during corrector state does not map onto (2.30), but rather onto a version of that equation with zeta[n] and zeta[n-1] replaced as stated above.
This is not a bug, however, but is rather a different algorithm.
In reality, this particular version of barotropic mode was created back in 1998, long before the paper was written (the paper was finalized in early 2003, but it took another two years to publish it).
Back in 1998 the whole apparatus of studying stability of time-stepping algorithms was not developed, and the coefficients 0.16 (see above) and 0.4 (weighting between predicted and corrected values of zeta during corrector stage) was chosen empirically by playing with oscillator system. The 2005 paper actually DERIVES a time stepping algorithm with the desired properties by optimising coefficients of characteristic equations in order to make them produce the desired eigen values. There was no such procedure back in 1998.
HOWEVER, THE PROPOSED FIX DESIGNED TO BRING THIS CODE IN CONSISTENCY WITH THE 2005 PAPER DOES NOT MAKE ANY SENSE WHATSOEVER.
THE WHOLE BAROTROPIC MODE SHOULD BE REPLACED INSTEAD.
The existing step2D codes (both in Rutgers ROMS and Agrif, all versions of both) are now helplessy obsolete.
What is proposed here is to compute barotropic pressure gradient terms twice during predictor substep, which simply adds more computational load into already inefficient code.
A new version of step2D, which is consistent with Eqs. (2.38)---(2.41) of 2005 paper was available in the middle of 2001 (in fact, it broadcasted to ROMS at the end of August 2001, soon after the Boulder meeting) and a newer Forward-Backward (now standard in UCLA codes) scheme followed in 2002. They are approximately 1.5 and 2.5 times more efficient than the original LF-AM3 code. (Due to combined effect of better stability limits and more efficient codes, for example, not storing r.h.s. terms between time steps).
Aparently, both Hernan and Agrif people (with some help from the community) did a very good job in consealing and, in fact, supressing their sole existence.
RE: 2D time stepping and 2D/3D coupling
The second aspect of this problem --- interference between computation of 3D--> baroclinic mode forcing term and the restart of barotropic mode --- is a more serious issue because all ROMS codes (UCLA, Rutgers, Agrif) are affected by this issue.
The problem comes from the fact that the full 3D pressure gradient term is computed assuming that free-surface elevation is exactly at time step "n", while the r.h.s. of 2D at the first barotropic time step is computed from free-surface updated by the first barotropic time step. As a result, when the fully computed r.h.s. of 2D mode is subtracted from vertically integrated r.h.s. of 3D,
the resultant forcing term implicitly contains the difference of barotropic pressure gradient terms computed by 3D and 2D modes, which contains error due to different time-placement of free surface. Specific errors in time-placement are:
The motivation to use updated free-surface during the first barotropic step comes from the desire to construct a stable startup algorithm, for example, in FB version a backward step for pressure gradient is the only alternative because the first step for free surface is always forward Euler. Changing backward step to forward would completely eliminate the error in computation of the coupling term, but it brings another source of instability on its own. Overall stability of the whole code depends on total balance between this startup instability and numerical dissipation during subsequent stable steps and dissipation by fast-time averaging filter. One way or the other (admitting unstable startup or coupling error), this leads to unstable code in the case when the mode-splitting ratio "ndtfast" becomes small. Choosing more dissipative weights always helps, but degrades accuracy (this matter was discussed previously on this board under "Time-averaging of barotropic fields and Power Law Filter").
This is a lesser problem in the case of LF-AM3: it is just a matter of re-tuning coefficients, because a perfectly stable algorithm can be obtained from a modified Runge Kutta (RK2) scheme, Eqs. (2.18 ) --- (2.20) by setting beta=0, epsilon=1, which results in a stability limit of 2 (with a characteristic equation exactly the same as classical FB --- not bad at all), and no coupling term error. [It is still a bit of a problem in the case of UCLA LF-AM3 version because the starting barotropic step ends up being run slightly outside the stability limit set by the regular time-stepping optimized for the largest stability range, say LF-TR with gamma=0., beta=0.14, epsil=0.74 recommended there].
A better, and more regular solution to eliminate the conflict between start up and coupling error is to compute barotropic pressure gradient term during the first barotropic step in two stages:
http://www.atmos.ucla.edu/~alex/ROMS/st ... _patch.tar
To test the effect of this change, I run my Pacific model using ndtfast = 61 (which is considered rather large), but I tuned my time filter to very aggressive settings:
p = 2, q = 4, r=0.28
where 0.28 is very close to the ideal value of 0.284. Note that the difference 0.284-r sets the magnitude of second-order dissipative truncation error; which is now nearly 8 times less than that of more typical recommended setting
p = 2, q = 4, r=0.25
The code is still stable with r=0.28, and I do not see abnormal oscillations in free surface. This was not achievable before.
See also "Time-averaging of barotropic fields and Power Law Filter" in this section below.
The problem comes from the fact that the full 3D pressure gradient term is computed assuming that free-surface elevation is exactly at time step "n", while the r.h.s. of 2D at the first barotropic time step is computed from free-surface updated by the first barotropic time step. As a result, when the fully computed r.h.s. of 2D mode is subtracted from vertically integrated r.h.s. of 3D,
Code: Select all
rufrc(i,j) = rufrc(i,j)-rubar(i,j)
rvfrc(i,j) = rvfrc(i,j)-rvbar(i,j)
Code: Select all
0.5*dt/ndtfast --- Rutgers codes (all versions); Agrif (except corrected by Laurent)
(1./3.)*dt/ndtfast --- UCLA step2D_LF_AM3.F
dt/ndtfast --- UCLA step2D_FB.FThis is a lesser problem in the case of LF-AM3: it is just a matter of re-tuning coefficients, because a perfectly stable algorithm can be obtained from a modified Runge Kutta (RK2) scheme, Eqs. (2.18 ) --- (2.20) by setting beta=0, epsilon=1, which results in a stability limit of 2 (with a characteristic equation exactly the same as classical FB --- not bad at all), and no coupling term error. [It is still a bit of a problem in the case of UCLA LF-AM3 version because the starting barotropic step ends up being run slightly outside the stability limit set by the regular time-stepping optimized for the largest stability range, say LF-TR with gamma=0., beta=0.14, epsil=0.74 recommended there].
A better, and more regular solution to eliminate the conflict between start up and coupling error is to compute barotropic pressure gradient term during the first barotropic step in two stages:
- at first use just zeta(:,:,kstp) to ensure exact consistency with 3D mode;
- then, after "rufrc,rvfrc" , which initially contain vertically integrated r.h.s. of 3D mode are converted into coupling terms by subtracting "rubar,rvbar", add correction based on the difference zeta_new(...) - zeta(... , kstp) to "rubar,rvbar" to make them consistent with the desired Backward step for pressure gradient terms (FB version), or 1/3 of that change (LF-AM3).
http://www.atmos.ucla.edu/~alex/ROMS/st ... _patch.tar
To test the effect of this change, I run my Pacific model using ndtfast = 61 (which is considered rather large), but I tuned my time filter to very aggressive settings:
p = 2, q = 4, r=0.28
where 0.28 is very close to the ideal value of 0.284. Note that the difference 0.284-r sets the magnitude of second-order dissipative truncation error; which is now nearly 8 times less than that of more typical recommended setting
p = 2, q = 4, r=0.25
The code is still stable with r=0.28, and I do not see abnormal oscillations in free surface. This was not achievable before.
See also "Time-averaging of barotropic fields and Power Law Filter" in this section below.