# Ocean Modeling Discussion

ROMS/TOMS

Search for:
 It is currently Fri Jul 19, 2019 4:39 pm

 All times are UTC
Author Message
 Post subject: 2 possible bugs in Shc boundary conditions Posted: Thu Mar 13, 2014 3:41 am Joined: Mon Oct 17, 2011 1:41 am
Posts: 14
Location: WHOI
I think I found 2 bugs in the code for the new Shchepetkin boundary conditions. In u2dbc_im.F for the eastern boundary, the current code has
Code:
cff=1.0_r8/(0.5_r8*(GRID(ng)%h(Iend  ,j)+                &
&                            GRID(ng)%h(Iend+1,j)))
cff1=SQRT(g*cff)
Cx=dt2d*cff1*0.5_r8*(GRID(ng)%pm(Iend  ,j)+              &
&                             GRID(ng)%pm(Iend+1,j))

cff1 is sqrt(g/h) but Cx (which is c_tilde in Mason et al. (2010) notation) needs sqrt(g*h). However, in the last line we do require sqrt(g/h) when calculating ubar - corresponding line is
Code:
ubar(Iend+1,j,kout)=0.5_r8*                              &
&                            ((1.0_r8-Cx)*ubar(Iend+1,j,know)+    &
&                             Cx*ubar(Iend,j,know)+               &
&                             bry_val+                            &
&                             cff1*(Zx-BOUNDARY(ng)%zeta_east(j)))

So the definition of Cx needs to change to
Code:
Cx=dt2d*SQRT(g*1/cff)*0.5_r8*(GRID(ng)%pm(Iend  ,j)+        &
&                             GRID(ng)%pm(Iend+1,j))

I'm not sure if correcting cff to use h+zeta instead of just h is required for good accuracy.

The second bug is for the case Cx > Co. In the code, Zx is reassigned the value of Zx + cff2*cff3 when this condition is true.
Code:
Zx=(0.5_r8+Cx)*zeta(Iend  ,j,know)+                       &
&           (0.5_r8-Cx)*zeta(Iend+1,j,know)
IF (Cx.gt.Co) THEN
cff2=(1.0_r8-Co/Cx)**2
cff3=zeta(Iend,j,kout)+                                &
&               Cx*zeta(Iend+1,j,know)-                           &
&               (1.0_r8-Cx)*zeta(Iend,j,know)
Zx=Zx+cff2*cff3
END IF

But, Zx+cff2*cff3 gives (0.5 + Cx - (1-Cx)*(1-Co/Cx)^2) * zeta(Iend,j,know) + two other terms. The other two terms have the correct coefficient when comparing against eqn (7) in Mason et al. (2010). But for zeta(Iend,j,know) we need (0.5 + Cx - (1+Cx)*(1-Co/Cx)^2) * zeta(Iend,j,know) for the code to agree with eqn (7). So, the required correction is
Code:
cff3=zeta(Iend,j,kout)+                                &
&               Cx*zeta(Iend+1,j,know)-                           &
&               (1.0_r8+Cx)*zeta(Iend,j,know)

Similar corrections are required for all 4 boundaries.

Deepak

Mason, E., Molemaker, J., Shchepetkin, A. F., Colas, F., McWilliams, J. C., & Sangrà, P. (2010). Procedures for offline grid nesting in regional ocean models. Ocean Modelling, 35(1-2), 1–15. doi:10.1016/j.ocemod.2010.05.007

Top Post subject: Re: 2 possible bugs in Shc boundary conditions Posted: Thu Mar 20, 2014 1:08 am Joined: Mon Oct 17, 2011 1:41 am
Posts: 14
Location: WHOI
I can now confirm that my application behaves a lot better when the two changes are made. Earlier, I had a bunch of waves running around my domain. Now, they exit the open boundaries very nicely.

Deepak

Top Post subject: Re: 2 possible bugs in Shc boundary conditions Posted: Fri Mar 21, 2014 3:05 am  Joined: Fri Nov 14, 2003 4:57 pm
Posts: 185
The original code in u2dbc_im.F is as follows, and it is correct:
Code:
cff=0.5*(h(iend,j)+h(iend+1,j))
hx=sqrt(g/cff)
cx=dtfast*cff*hx*0.5*(pm(iend,j)+pm(iend+1,j))

zx=(0.5+cx)*zeta(iend,j,kstp)+(0.5-cx)*zeta(iend+1,j,kstp)
if (cx > 0.292893218813452) then
zx=zx + ( zeta(iend,j,knew) +cx*zeta(iend+1,j,kstp)
&                               -(1.+cx)*zeta(iend,j,kstp)
&                           )*(1.-0.292893218813452/cx)**2
endif

ubar(iend+1,j,knew)=0.5*( (1.-cx)*ubar(iend+1,j,kstp)
&                                     +cx*ubar(iend,j,kstp)
#    ifdef M2_FRC_BRY
&                       +ubar_east(j) +hx*(zx-zeta_east(j))
#    else
&                   +ubclm(iend+1,j) +hx*(zx-ssh(iend+1,j))
#    endif
&                                                        )

you may use it as the reference.

It also comes with a special test problem WAVE_RAD wich can be used to verify
that all 4 boundaries work correctly and result in very little reflection.

Notes:

1. hx=sqrt(g/cff), that is division by cff, not multiplication, however

2. cx=dtfast*cff*hx*0.5*(pm(iend,j)+pm(iend+1,j)) contains an extra multiplication
by cff, hence cx=dtfast*cff*sqrt(g/cff)*.... = dtfast*sqrt(g*cff)*..., so it contains
sqrt(g*cff)=sqrt(g*h) which is phase speed of barotropic gravity wave. The resultant
overall value of cx is nondimensional, as it should be, because it is Courant number.

3. The above 1. and 2. are constructed in such a way that there is only one division and
one sqrt operation (both division and sqrt are approximately of the theme computational
cost, and are substantially more expensive than multiplicationand addition. This is the
rationale of doing it this way and not the other way around). Note that BOTH hx and cx
are needed by the computations below, and any other variant of their computation
[besides computing sqrt(g/h)] will result in an extra sqrt or an extra division.

4. if (cx > 0.292893218813452) then maps onto cx > Co in Rutgers code, and zx is
modified inside this logical branch in such a way that the additional term is quadratic
with thespect to cx-Co.

5. It can easily be worked with paper and pencil, that after modification 4. to zx
is applied, zx is a weighted sum of zeta(iend,j,kstp), zeta(iend+1,j,kstp), and
zeta(iend,j,knew) with POSITIVE coefficients in front of all three such that the
three coefficients ADD UP TO 1.

Quote:
I'm not sure if correcting cff to use h+zeta instead of just h is required for good accuracy.

Generally speaking the answer is NO, but this judgment is close to neutral than categoric.

These kind boundary conditions (I hate to call them "Flather" because calling them this
way is historically incorrect) are derived under the assumptions of

(i) linearization (left- and right-travelling waves are moving independently from each
other); thought the linearization does not necessarily have to be around zeta=0, it can
be viewed as a small perturbation relative to a final background state (e.g., consider
small short waves on top of large-amplitude tidal elevations and depressions of free
surface).

(ii) The waves are assumed to propagate in the direction normal to the boundary
(which is never true in practice, hence the physical accuracy of these B.C. is
fundamentally limited);

(iii) The Coriolis term is neglected altogether, so the actual wave speed of surface
gravity wave is slightly larger than sqrt(gh), however there is no evident way to
fix this, because in the presense of Coriolis force waves become dispersive and
these is no self-evident way to decompose them into left and right travelling waves;

The limitations (ii) anf (iii) are offset by the fact that the role of these B.C.
condition is to merely provide "radiative cooling" of the solution, or, simply put,
avoid trapping of surface gravity waves inside the domain. The easiest way to
understand this is to compare this with "clamped" boundares. If normal velocity is
specified at the boundary, then waves will be 100% reflected back into the domain,
and the entire solution will have seiche-mode oscillation in free surface which will
never go away. In contrast, a radiation condition (not necessarily perfect
non-reflection) makes seiche mode go away rather quickly -- it cannot sustain itself
as lon as a significant portion of wave is let out of the domain every thime.

Top Post subject: Re: 2 possible bugs in Shc boundary conditions Posted: Wed Mar 26, 2014 3:44 am Joined: Mon Oct 17, 2011 1:41 am
Posts: 14
Location: WHOI
Thanks!

The changes I made agree with the code you posted. I like the "radiation cooling" analogy.

Deepak

Top Post subject: Re: 2 possible bugs in Shc boundary conditions Posted: Wed Mar 26, 2014 8:08 pm  Joined: Fri Nov 14, 2003 4:57 pm
Posts: 185
code looks like
Code:
cff=1.0_r8/(0.5_r8*(GRID(ng)%h(Iend  ,j)+                 &   ! <-- A
&                            GRID(ng)%h(Iend+1,j)))
cff1=SQRT(g*cff)                                                ! <-- B
Cx=dt2d*cff1*0.5_r8*(GRID(ng)%pm(Iend  ,j)+               &    ! <-- C
&                             GRID(ng)%pm(Iend+1,j))
Zx=(0.5_r8+Cx)*zeta(Iend  ,j,know)+                       &
&           (0.5_r8-Cx)*zeta(Iend+1,j,know)
IF (Cx.gt.Co) THEN
cff2=(1.0_r8-Co/Cx)**2
cff3=zeta(Iend,j,kout)+                                 &
&               Cx*zeta(Iend+1,j,know)-                            &
&               (1.0_r8-Cx)*zeta(Iend,j,know)                         ! <-- D
Zx=Zx+cff2*cff3
END IF
ubar(Iend+1,j,kout)=0.5_r8*                               &
&                            ((1.0_r8-Cx)*ubar(Iend+1,j,know)+     &
&                             Cx*ubar(Iend,j,know)+                &
&                             bry_val+                             &
&                             cff1*(Zx-BOUNDARY(ng)%zeta_east(j)))   ! <-- E

where by A,B,C,D,E I highlighted the lines which need to be paid attention to.

There are two errors in the code above: (1) Cx computed by line C is incorrect, and also dimensionally incorrect.
It must be non-dimensional, while it is 1/Length if computed as above; and (2) the three weighting coefficients
for the three free-surface values in computation of cff3 (ending at line D) should add up to 0. Currently they
are not.

The corrected version is
Code:
cff=0.5_r8*(GRID(ng)%h(Iend  ,j)+                        &   ! <-- A
&                    GRID(ng)%h(Iend+1,j))
cff1=SQRT(g/cff)                                                  ! <-- B
Cx=dt2d*cff1*cff* 0.5_r8*(GRID(ng)%pm(Iend  ,j)+         &      ! <-- C
&                                   GRID(ng)%pm(Iend+1,j))
Zx=(0.5_r8+Cx)*zeta(Iend  ,j,know)+                       &
&           (0.5_r8-Cx)*zeta(Iend+1,j,know)
IF (Cx.gt.Co) THEN
cff2=(1.0_r8-Co/Cx)**2
cff3=zeta(Iend,j,kout)+                                 &
&               Cx*zeta(Iend+1,j,know)-                            &
&               (1.0_r8+Cx)*zeta(Iend,j,know)                         ! <-- D
Zx=Zx+cff2*cff3
END IF
ubar(Iend+1,j,kout)=0.5_r8*                               &
&                            ((1.0_r8-Cx)*ubar(Iend+1,j,know)+     &
&                             Cx*ubar(Iend,j,know)+                &
&                             bry_val+                             &
&                             cff1*(Zx-BOUNDARY(ng)%zeta_east(j)))  ! <-- E

A. changed 1.0_r8/(0.5_r8*(...)) into just 0.5_r8*(...)

B. changed SQRT(g*cff) into SQRT(g/cff)

C. introduced an extra multiplication by cff just after cff1. Note: changes A,B,C applied together
do not actually change the computed value of cff1, however Cx is different.

D. changed (1.0_r8-Cx) into (1.0_r8+Cx) Note: this is unrelated to A,B,C,E. The three weighting
coefficients for the three free-surface values should add up to 0.

E. there is no need for any change, however note that the meaning of cff1 must be preserved in changes A,B.

Similar corrections need to be applied to western edge boundary conditions in u2dbc_im.F and in southern
and northern edges in v2dbc_im.F

AGRIF and UCLA codes are not affected.

Top Post subject: Re: 2 possible bugs in Shc boundary conditions Posted: Fri Mar 28, 2014 10:03 pm Site Admin Joined: Wed Feb 26, 2003 4:41 pm
Posts: 1078
Location: IMCS, Rutgers University
Indeed, I fixed those bugs. Please update. Thank you for bringing this to my attention. I also corrected the tangent linear, representer, and adjoint versions of u2dbc_im.F and v2dbc_im.F.

Top Post subject: Re: 2 possible bugs in Shc boundary conditions Posted: Mon Mar 31, 2014 6:16 pm  Joined: Wed Jul 02, 2003 5:29 pm
Posts: 3633
Location: IMS/UAF, USA
The WET_DRY case still needs cff=(h+zeta) for cases in which h goes negative. The following sqrt operation behaves badly if cff is negative.

Top Display posts from previous: All posts1 day7 days2 weeks1 month3 months6 months1 year Sort by AuthorPost timeSubject AscendingDescending

 All times are UTC

#### Who is online

Users browsing this forum: No registered users and 3 guests

 You cannot post new topics in this forumYou cannot reply to topics in this forumYou cannot edit your posts in this forumYou cannot delete your posts in this forumYou cannot post attachments in this forum

Search for:
 Jump to:  Select a forum ------------------ News, Events & Job Opportunities    Meetings/Workshops    Job Opportunities    Ocean News ROMS/TOMS    ROMS Adjoint    ROMS Benchmarks    ROMS Bugs    ROMS Discussion    ROMS Documentation    ROMS Ecosystem    ROMS FAQ    ROMS Ice    ROMS Information    ROMS Installation    ROMS Known Problems    ROMS Messages    ROMS Problems    ROMS Releases    ROMS Results    ROMS Sediment    ROMS Source    ROMS Tips    ROMS Tools and Techniques    ROMS Trivia    ROMS Usage    ROMS Webinar    ROMS Wish List ROMS/TOMS Applications    User Applications    Adriatic Sea Ocean Modeling    Ocean Modeling FAQ