Question about Flather 2d boundary condition for ubar

[stef]

I'm trying to write up the boundary condition I use and I wanted to confirm whether the following is correct:

Code: Select all

!
!  Western edge, Flather boundary condition.
!
        ELSE IF (LBC(iwest,isUbar,ng)%Flather) THEN
          DO j=Jstr,Jend
            IF (LBC_apply(ng)%west(j)) THEN
#if defined SSH_TIDES && !defined UV_TIDES
              IF (LBC(iwest,isFsur,ng)%acquire) THEN
                bry_pgr=-g*(zeta(Istr,j,know)-                          &
     &                      BOUNDARY(ng)%zeta_west(j))*                 &
     &                  0.5_r8*GRID(ng)%pm(Istr,j)
              ELSE

Is the factor 0.5 right? Shouldn't there be a factor 2 instead? Unfortunately I haven't read the Mason et al. (2010) paper because I have no access. Maybe someone could send it to me (stefan at sriha.net)?? Thanks for your help and sorry for the inconvenience.

EDIT: Dr. Martyanov kindly sent me the paper via email! Thanks so much!

Regards, Stefan

Mason E., J. Molemaker, A.F. Shchepetkin, F. Colas, J.C. McWilliams, P. Sangrà, 2010: Procedures for offline grid nesting in regional ocean models, Ocean Modelling, 35, 1-15.

[wilkin]

I believe the 1/2 is there because, while the staggered grid formally places zeta_west 1/2 cell to the "left" of the u-face point (this is the code block where UV_TIDES is undefined so we are substituting a dynamical approximation for it), when zeta_west and ubar_west are created by users they are treated as being co-located, not staggered. So, to get d(zeta)/dx, the dx distance is only half of 1/pm.

[stef]

I had the same thought about the staggering of zeta_west, i.e. that it's centered on u-points. This is supported by the comments/code in set_tides.F, see below. zeta_w+=0.5(zeta_{i-1}+zeta_i).

But assuming that's right, a forward difference would be (zeta_i-zeta_w)/(0.5*dx) = 2*pm*(zeta_i-zeta_w).

Code: Select all

!
!  If appropriate, load tidal forcing into boundary arrays.  The "zeta"
!  boundary arrays are important for the Flather or reduced physics
!  boundary conditions for 2D momentum. To avoid having two boundary
!  points for these arrays, the values of "zeta_west" and "zeta_east"
!  are averaged at u-points.  Similarly, the values of "zeta_south"
!  and "zeta_north" is averaged at v-points. Noticed that these
!  arrays are also used for the clamped conditions for free-surface.
!  This averaging is less important for that type ob boundary
!  conditions.
!
        IF (LBC(iwest,isFsur,ng)%acquire.or.                            &
     &      LBC(iwest,isUbar,ng)%acquire.or.                            &
     &      LBC(iwest,isVbar,ng)%acquire) THEN
          update=.FALSE.
          IF (DOMAIN(ng)%Western_Edge(tile)) THEN
            DO j=JstrR,JendR
#  ifdef ADD_FSOBC
              BOUNDARY(ng)%zeta_west(j)=BOUNDARY(ng)%zeta_west(j)+      &
     &                                  0.5_r8*(Etide(Istr-1,j)+        &
     &                                          Etide(Istr  ,j))

[arango]

It can go both ways, but we need to test. The section of the code in question is rarely used since we specify both tidal elevation and tidal currents. For the case that only tidal elevation is activated, we need to approximate the tidal currents with reduced-physics terms evaluated at the open boundary: barotropic pressure gradient, Coriolis, and stresses. Here, bry_pgr (m3/s2) is the barotropic pressure gradient.

Thus, we either leave it as it is:

Code: Select all

#if defined SSH_TIDES && !defined UV_TIDES
              IF (LBC(iwest,isFsur,ng)%acquire) THEN
                bry_pgr=-g*(zeta(Istr,j,know)-                          &
     &                      BOUNDARY(ng)%zeta_west(j))*                 &
     &                  0.5_r8*GRID(ng)%pm(Istr,j)
              ELSE
                bry_pgr=-g*(zeta(Istr  ,j,know)-                        &
     &                      zeta(Istr-1,j,know))*                       &
     &                  0.5_r8*(GRID(ng)%pm(Istr-1,j)+                  &
     &                          GRID(ng)%pm(Istr  ,j))
              END IF

remove the 0.5 factor:

Code: Select all

#if defined SSH_TIDES && !defined UV_TIDES
              IF (LBC(iwest,isFsur,ng)%acquire) THEN
                bry_pgr=-g*(zeta(Istr,j,know)-                          &
     &                      BOUNDARY(ng)%zeta_west(j))*                 &
     &                  GRID(ng)%pm(Istr,j)
              ELSE
                bry_pgr=-g*(zeta(Istr  ,j,know)-                        &
     &                      zeta(Istr-1,j,know))*                       &
     &                  0.5_r8*(GRID(ng)%pm(Istr-1,j)+                  &
     &                          GRID(ng)%pm(Istr  ,j))
              END IF

Or average the grid spacing:

Code: Select all

#if defined SSH_TIDES && !defined UV_TIDES
              IF (LBC(iwest,isFsur,ng)%acquire) THEN
                bry_pgr=-g*(zeta(Istr,j,know)-                          &
     &                      BOUNDARY(ng)%zeta_west(j))*                 &
     &                  0.5_r8*(GRID(ng)%pm(Istr-1,j)+                  &
     &                          GRID(ng)%pm(Istr  ,j))
              ELSE
                bry_pgr=-g*(zeta(Istr  ,j,know)-                        &
     &                      zeta(Istr-1,j,know))*                       &
     &                  0.5_r8*(GRID(ng)%pm(Istr-1,j)+                  &
     &                          GRID(ng)%pm(Istr  ,j))
              END IF

Similar treatment can be done on the ubar eastern boundary, and southern/northern vbar boundaries.

Of course, this needs to be tested on an application. I am swamped now. Perhaps, you can give it a try and let us know.

[stef]

Thanks for the reply! Ok, I'll test it. Along with testing this, I should test prescribing the barotropic velocity. First I wanted to document carefully what I have been doing up to now (reduced physics). I attached an excerpt from a note I'm writing, in case anybody is interested in this topic.

The references in the attachment are not rendered properly:

blayo05: E. Blayo, L. Debreu, Revisiting open boundary conditions from the point of view of characteristic variables, Ocean Modelling, Volume 9, Issue 3, 2005
flather75: Flather, R. A., & Davies, A. M. (1975). The application of numerical models to storm surge prediction.

[wilkin]

Stefan,

I see what you are saying. The reduced physics option to estimate ubar(Istr) (the physical boundary u-point) for use in Flather (when ubar_west is not provided) is solving d(ubar)/dt = -g*d(zeta)/dx and there appears to be a factor of 1/2 multiplying the pressure gradient RHS term.

This is ultimately completed in this line

Code: Select all

              ubar(Istr,j,kout)=ubar(Istr,j,know)+                      &
     &                          dt2d*(bry_pgr+                          &
     &                                bry_cor+                          &
     &                                bry_str)

I think the answer might lie in subtleties of the time stepping.

At the top of u2dbc_im.F there is a block:

Code: Select all

!  Set time-indices
!-----------------------------------------------------------------------
!
      IF (FIRST_2D_STEP) THEN
        know=krhs
        dt2d=dtfast(ng)
      ELSE IF (PREDICTOR_2D_STEP(ng)) THEN
        know=krhs
        dt2d=2.0_r8*dtfast(ng)
      ELSE
        know=kstp
        dt2d=dtfast(ng)
      END IF

where in the PREDICTOR_2D_STEP dt2d is two times the fast DT. This will cancel the factor of 1/2 in bry_pgr assuming kout is the new time step (n+1) and know is krhs if that is (n-1).

I haven't followed through the logic of what happens on the "corrector" step that follows, where dt2d is one times DT fast and know is kstp but I think you have to look at the code in the context of these aspects of the predictor/corrector time stepping of the barotropic equations.

[stef]

Thanks for the tip, yes, I have not looked into the time stepping.

I will write into my note that I don't fully understand it yet. Maybe I can figure it out in spring or summer when it's less busy.

But for the record: I didn't see your first code block in the Flather condition, only in the reduced-physics. Maybe I'm on an older version. For me in Flather it says

Code: Select all

              cff2=GRID(ng)%om_u(Istr,j)*Cx
!!            cff2=dt2d
              bry_val=ubar(Istr+1,j,know)+                              &
     &                cff2*(bry_pgr+                                    &
     &                      bry_cor+                                    &
     &                      bry_str)

But your hypothesis that the time stepping solves the apparent contradiction may still be true. However, wouldn't that make the second option (the ELSE clause) invalid:

Code: Select all

              IF (LBC(iwest,isFsur,ng)%acquire) THEN
                bry_pgr=-g*(zeta(Istr,j,know)-                          &
     &                      BOUNDARY(ng)%zeta_west(j))*                 &
     &                  0.5_r8*GRID(ng)%pm(Istr,j)
              ELSE
                bry_pgr=-g*(zeta(Istr  ,j,know)-                        &
     &                      zeta(Istr-1,j,know))*                       &
     &                  0.5_r8*(GRID(ng)%pm(Istr-1,j)+                  &
     &                          GRID(ng)%pm(Istr  ,j))
              END IF

The two seem inconsistent, regardless of which factor they are later multiplied with.

[arango]

The 0.5 factor is associated with the dx when we evaluate d(zeta)/dx and noting to do with time stepping. The C-grid stencil is:

Code: Select all

          .--------- .-------- vbar(Istr,j+1)-------.
          +          |                              |
          +          |                              |
   zeta(west,j)  ubar(Istr,j)  zeta(Istr,j)     ubar(Istr+1,j)
          +          |                              |
          +          |                              |
          .--------- .-------- vbar(Istr,j)---------.

The distance between zeta(west,j) and zeta(Istr,j) is always a full pm. Thus, for me using:

Code: Select all

#if defined SSH_TIDES && !defined UV_TIDES
              IF (LBC(iwest,isFsur,ng)%acquire) THEN
                bry_pgr=-g*(zeta(Istr,j,know)-                          &
     &                      BOUNDARY(ng)%zeta_west(j))*                 &
     &                  0.5_r8*(GRID(ng)%pm(Istr-1,j)+                  &
     &                          GRID(ng)%pm(Istr  ,j))
              ELSE
                bry_pgr=-g*(zeta(Istr  ,j,know)-                        &
     &                      zeta(Istr-1,j,know))*                       &
     &                  0.5_r8*(GRID(ng)%pm(Istr-1,j)+                  &
     &                          GRID(ng)%pm(Istr  ,j))
              END IF

makes more sense. Therefore, the spatial derivative d(zeta) is centered at the ubar(Istr,j) physical boundary.

[stef]

The distance between zeta(west,j) and zeta(Istr,j) is always a full pm.

Now I'm confused. Isn't this inconsistent with set_tides.F? Isn't zeta(west,j) already set to be centered at u-points:

Code: Select all

322	!  If appropriate, load tidal forcing into boundary arrays.  The "zeta"
323	!  boundary arrays are important for the Flather or reduced physics
324	!  boundary conditions for 2D momentum. To avoid having two boundary
325	!  points for these arrays, the values of "zeta_west" and "zeta_east"
326	!  are averaged at u-points.  

[arango]

Yes, there is an inconsistency in my explanation above. I am trying to show you that 0.5 is associated with x-spacing and not a time-step factor. If the zeta_west is averaged at the u-point, its usage in the zetabc routine for the boundary conditions for zeta is halfway in, but it won't be much trouble. It has been a long time since this was coded. So, it makes sense to use the 0.5 when computing bry_pgr, but the backward spatial-derivative won't be centered (interpolated) at the ubar boundary. After all, it is an approximation! If you are not happy with this approximation, use the tidal currents too! It is the most appropriate way to force ROMS with tides (elevation and currents). Also, if you check step2d, you will notice that the barotropic pressure gradient is more elaborated.

[stef]

Ok, thanks for taking the time to explain it! I will test the various options soon. I just started the thread because I was worried that 1/0.5=2 was mistakenly coded as 0.5, but in hindsight I learned a lot.