A bug in the 2D barotropic BC routines?

[lanerolle]

I have been trying to reproduce a result from an older version of ROMS (version 92) using a newer version (version 526). The simulation is pretty simple in that it only has tidal forcing with the M2 consituent and has constant T (20C) and S (35 PSU) initial conditions. There is no river, met/atmos or open boundary non-tidal/T/S forcing. So the only things which evolve are water levels and currents and T and S remain fixed. However, the two ROMS versions give different results and I have not been able to reconcile them until I began to look at the code and I found the following :

ROMS version 92 :

Code: Select all

# elif defined SOUTH_M2FLATHER
!
!  Southern edge, Flather boundary condition.
!
        DO i=Istr,Iend
#  if !defined SSH_TIDES && !defined UV_TIDES
#   ifdef FSOBC_REDUCED
          bry_pgr=-g*(zeta(i,Jstr,know)-                                &
     &                BOUNDARY(ng)%zeta_south(i))*                      &
     &            0.5_r8*GRID(ng)%pn(i,Jstr)
#   else

ROMS version 526:

Code: Select all

# elif defined SOUTH_M2FLATHER
!
!  Southern edge, Flather boundary condition.
!
        DO i=Istr,Iend
#  if defined SSH_TIDES && !defined UV_TIDES
#   ifdef FSOBC_REDUCED
          bry_pgr=-g*(zeta(i,Jstr,know)-                                &
     &                BOUNDARY(ng)%zeta_south(i))*                      &
     &            0.5_r8*GRID(ng)%pn(i,Jstr)
#else

As you can see, one has "defined SSH_TIDES" and the other has "!defined SSH_TIDES".

Is this a bug? Which version is correct?

[kate]

I believe the intention is consistent with the new code. The time to use reduced boundary condition is when you have information about surface elevation but not barotropic currents. The new version allows you to set the SSH_TIDES without the UV_TIDES and still approximate a tidal flow.

Still, ROMS is tricksy and there are times when you just have to read the .f90 files to see if you got what you wanted.

P.S. Please don't double post.

[lanerolle]

The issue is the we force tides using an open boundary forcing NetCDF file with water level and barotropic current series in it. So we do NOT use the SSH_TIDES or UV_TIDES CPP options because we do not use a tidal harmonics NetCDF file. So the two code segments in u2dbc_im.f90 become :

ROMS version 92 :

Code: Select all

!
!-----------------------------------------------------------------------
!  Lateral boundary conditions at the eastern edge.
!-----------------------------------------------------------------------
!
      IF (Iend.eq.Lm(ng)) THEN
!
!  Eastern edge, Flather boundary condition.
!
        DO j=Jstr,Jend
          bry_pgr=-g*(zeta(Iend+1,j,know)-                              &
     &                zeta(Iend  ,j,know))*                             &
     &            0.5_r8*(GRID(ng)%pm(Iend  ,j)+                        &
     &                    GRID(ng)%pm(Iend+1,j))
          bry_cor=0.125_r8*(vbar(Iend  ,j  ,know)+                      &
     &                      vbar(Iend  ,j+1,know)+                      &
     &                      vbar(Iend+1,j  ,know)+                      &
     &                      vbar(Iend+1,j+1,know))*                     &
     &                     (GRID(ng)%f(Iend  ,j)+                       &
     &                      GRID(ng)%f(Iend+1,j))
          cff1=1.0_r8/(0.5_r8*(GRID(ng)%h(Iend  ,j)+                    &
     &                         zeta(Iend  ,j,know)+                     &
     &                         GRID(ng)%h(Iend+1,j)+                    &
     &                         zeta(Iend+1,j,know)))
          bry_str=cff1*(FORCES(ng)%sustr(Iend+1,j)-                     &
     &                  FORCES(ng)%bustr(Iend+1,j))
          Cx=1.0_r8/SQRT(g*0.5_r8*(GRID(ng)%h(Iend+1,j)+                &
     &                             zeta(Iend+1,j,know)+                 &
     &                             GRID(ng)%h(Iend  ,j)+                &
     &                             zeta(Iend  ,j,know)))
          cff2=GRID(ng)%om_u(Iend+1,j)*Cx
          bry_val=ubar(Iend,j,know)+                                    &
     &            cff2*(bry_pgr+                                        &
     &                  bry_cor+                                        &
     &                  bry_str)
          cff=1.0_r8/(0.5_r8*(GRID(ng)%h(Iend  ,j)+                     &
     &                        zeta(Iend  ,j,know)+                      &
     &                        GRID(ng)%h(Iend+1,j)+                     &
     &                        zeta(Iend+1,j,know)))
          Cx=SQRT(g*cff)
          ubar(Iend+1,j,kout)=bry_val+                                  &
     &                        Cx*(0.5_r8*(zeta(Iend  ,j,know)+          &
     &                                    zeta(Iend+1,j,know))-         &
     &                            BOUNDARY(ng)%zeta_east(j))
          ubar(Iend+1,j,kout)=ubar(Iend+1,j,kout)*                      &
     &                        GRID(ng)%umask(Iend+1,j)
        END DO
      END IF

ROMS version 526 :

Code: Select all

!
!-----------------------------------------------------------------------
!  Lateral boundary conditions at the eastern edge.
!-----------------------------------------------------------------------
!
      IF (Iend.eq.Lm(ng)) THEN
!
!  Eastern edge, Flather boundary condition.
!
        DO j=Jstr,Jend
          bry_val=BOUNDARY(ng)%ubar_east(j)
          cff=1.0_r8/(0.5_r8*(GRID(ng)%h(Iend  ,j)+                     &
     &                        zeta(Iend  ,j,know)+                      &
     &                        GRID(ng)%h(Iend+1,j)+                     &
     &                        zeta(Iend+1,j,know)))
          Cx=SQRT(g*cff)
          ubar(Iend+1,j,kout)=bry_val+                                  &
     &                        Cx*(0.5_r8*(zeta(Iend  ,j,know)+          &
     &                                    zeta(Iend+1,j,know))-         &
     &                            BOUNDARY(ng)%zeta_east(j))
          ubar(Iend+1,j,kout)=ubar(Iend+1,j,kout)*                      &
     &                        GRID(ng)%umask(Iend+1,j)
        END DO
      END IF

As you can see, the older version of ROMS brings in contributions from Coriolis, pressure gradient, surface and bottom fluxes when imposing the BCs and these contributions are MISSING from the new ROMS.

So how do we reconcile these two versions? What CPP options should we employ if we are forcing using water level and barotropic currents time-series from a NetCDF file as opposed to using a NetCDF file containing tidal harmonic constituents?

[jcwarner]

It seems that the main difference here is:

older version:
bry_val=ubar(Iend,j,know)+ &
& cff2*(bry_pgr+ &
& bry_cor+ &
& bry_str)

newer version:
bry_val=BOUNDARY(ng)%ubar_east(j)

So if you have a netcdf file with ubar_east in it, then the newer version seems to be allowing that option.

[arango]

You will be never be able to reproduce the results because bugs have been corrected in file u2dbc_im.F since version 92. For example, see ticket 118. We update the code to either correct bugs or add new development/improvements. So will always advice the Users to keep the codes up to date.

I highly recommend Users to check trac when questions like the one raised here is asked. The trac web interface has a complete history of the changes of each file in ROMS repository. See Browse Source menu bottom. You will notice that you enter into the ROMS repository directory tree (trunk) and can browse any of the files. The web page for the file in question is here. Then, in the upper right side you will see the link for Last Change in red. If you click on this link you will be able to browse the entire history of changes to this file, one change at the time backwards. You can customize how to display the differences between revisions for this routine. I like to set the side by side menu when displaying the differences. If you want to know why that change was done, click on the trac ticket under the Message heading. Everything that you need to know about the evolution of ROMS is in the trac web site. I wonder how many user actually check this...

I see that you have hacked this routine because we never have the following directive:

Code: Select all

#  if !defined SSH_TIDES && !defined UV_TIDES

As Kate mentioned, there is a lot of flexibility here. If only SSH_TIDES is activated, the tidal forcing velocities are approximated by reduced physics. We actually do not recommend this in realistic applications. There has been a lot of discussions about this in the forum.

If you don't activate both SSH_TIDES and UV_TIDES, why are you setting Flather boundary conditions like this? I will use the values in the open boundary structure arrays. Seems to me that you need instead any of the ****_M2REDUCED options. Please, check this routine carefully and you will see what I am taking about. As John mentioned, there are countless options by manipulating the open boundary structure arrays BOUNDARY(ng)%. We use this extensively in nesting.