Slipperiness, psi-masking and analytical grid

[mmagaldi]

Hi everybody!

I am setting up an idealized grid experiment where I want to control the Reynolds number. I also want no slip condition on the boundaries. The cppdefs of interest are:

Code: Select all

#define  ANA_GRID
#define  MASKING
#define  UV_VIS2
#define  TS_DIF2

The psi-mask should be equal to 2 at the boundaries for the slipperiness and this is important for the calculations done in the momentum equations (e.g. in step2d.F). However, maybe, we have a problem.
In both metrics.F and analytical.F the following lines are defined:

Code: Select all

      
     IF (gamma2.lt.0.0_r8) THEN
        DO j=Jstr,Jend
          DO i=Istr,Iend
            pmask(i,j)=2.0_r8-pmask(i,j)
          END DO
        END DO
      END IF

and this means that the operation of subtracting pmask from 2 is performed twice. This results in psi-mask equal to zero instead of 2 on the land points, i.e. exactly as if we had free-slip conditions! Is this correct or am I wrong somewhere? Should we not perform the operation only once?
Thanks in advance to anyone who is going to reply.
M.

[kate]

Good catch! It looks like an honest to goodness bug. I would delete those lines from analytical.F and leave it in metrics.F since you want it to happen for all cases.

[shchepet]

As a side remark here is that setting pmask=2 is acceptable along straight boundaries
to account for the fact that tangential velocity is located half-grid interval from the boundary,
and therefore

Code: Select all

                pmask(i,j)*(u(i,j,k) - u(i,j-1,k))

accurate evaluates derivarive near the boundary, if one of the points, say u(i,j-1,k)=0
because it is masked, while the other, u(i,j,k) is inside the domain, but half grid interval
from the boundary. Then pmask=2 recovers the factor of 2.

However, doing so is wrong near the coastline which is eiter 45-degree to the grid primary
direction (hence the coastline passes through both u- and v-points), or near convex
corners because it leads to an overestimate of tangential stress by a factor of 2.

Therefore, the code mentioned here,

Code: Select all

IF (gamma2.lt.0.0_r8) THEN
        DO j=Jstr,Jend
          DO i=Istr,Iend
            pmask(i,j)=2.0_r8-pmask(i,j)
          END DO
        END DO
      END IF 

should be replaced with

Code: Select all

!
! Set no-slip boundary conditions on land-mask boundaries
! regardless of supplied value of gamma2.
!

          cff1=1.  !<-- computation of off-diagonal nonlinear terms
          cff2=2.

          if (rmask(i-1,j  ).gt.0.5 .and. rmask(i,j  ).gt.0.5 .and.
     &        rmask(i-1,j-1).gt.0.5 .and. rmask(i,j-1).gt.0.5) then
            pmask(i,j)=1.

          elseif(rmask(i-1,j  ).lt.0.5 .and.rmask(i,j  ).gt.0.5 .and.
     &           rmask(i-1,j-1).gt.0.5 .and.rmask(i,j-1).gt.0.5) then
            pmask(i,j)=cff1
          elseif(rmask(i-1,j  ).gt.0.5 .and.rmask(i,j  ).lt.0.5 .and.
     &           rmask(i-1,j-1).gt.0.5 .and.rmask(i,j-1).gt.0.5) then
            pmask(i,j)=cff1
          elseif(rmask(i-1,j  ).gt.0.5 .and.rmask(i,j  ).gt.0.5 .and.
     &           rmask(i-1,j-1).lt.0.5 .and.rmask(i,j-1).gt.0.5) then
            pmask(i,j)=cff1
          elseif(rmask(i-1,j  ).gt.0.5 .and.rmask(i,j  ).gt.0.5 .and.
     &           rmask(i-1,j-1).gt.0.5 .and.rmask(i,j-1).lt.0.5) then
            pmask(i,j)=cff1


          elseif(rmask(i-1,j  ).gt.0.5 .and.rmask(i,j  ).lt.0.5 .and.
     &           rmask(i-1,j-1).gt.0.5 .and.rmask(i,j-1).lt.0.5) then
            pmask(i,j)=cff2
          elseif(rmask(i-1,j  ).lt.0.5 .and.rmask(i,j  ).gt.0.5 .and.
     &           rmask(i-1,j-1).lt.0.5 .and.rmask(i,j-1).gt.0.5) then
            pmask(i,j)=cff2
          elseif(rmask(i-1,j  ).gt.0.5 .and.rmask(i,j  ).gt.0.5 .and.
     &           rmask(i-1,j-1).lt.0.5 .and.rmask(i,j-1).lt.0.5) then
            pmask(i,j)=cff2
          elseif(rmask(i-1,j  ).lt.0.5 .and.rmask(i,j  ).lt.0.5 .and.
     &           rmask(i-1,j-1).gt.0.5 .and.rmask(i,j-1).gt.0.5) then
            pmask(i,j)=cff2

          else
            pmask(i,j)=0.
          endif

which is a part of "setup_grid1.F" (or "metrics.F" in Rutgers codes).

Originally this matter was discovered in late 2004 during study of island wakes which it its
simplest configuration boils down to flow around cylinder made of masked points. It was
found that uniformy setting pmask=2 whenever it touches land makes the cylinger island
behave more like steppy polygon with vorticity generated predominantly at the points
where boundary bends. The above code, with pmask=1 at corner points and 2 at straight
90-degree segments gives smoother flow.

This was reported back then, but somehow slipped through and did not appear in the
officially released codes.

[mmagaldi]

Thank you to all who replied. On the same subject I have another question... What is the effect of land when we run just with implicit diffusion? And what does it mean running just with implicit diffusion? I mean I saw some application (see ADRIA02) just using as cppdefs

Code: Select all

#undef  UV_VIS2         /* turn ON or OFF Laplacian horizontal mixing */
#undef  UV_VIS4         /* turn ON or OFF biharmonic horizontal mixing */
#undef  TS_DIF2         /* turn ON or OFF Laplacian horizontal mixing */
#undef  TS_DIF4         /* turn ON or OFF biharmonic horizontal mixing */

Other instead

Code: Select all

#define  UV_VIS2         /* turn ON or OFF Laplacian horizontal mixing */
#define  TS_DIF2         /* turn ON or OFF Laplacian horizontal mixing */

and putting in the ocean_*.in file

Code: Select all

TNU2 == 0.0d0  0.0d0                    ! m2/s
VISC2 == 0.0d0                           ! m2/s

or similar for the biharmonic viscosity instead of the laplacian one.

Looking at the threads above I think we have always to define the switches for the horizontal mixing and, in the case we want implicit diffusion, put the viscosity and diffusivity coefficients equal to zero. If this is true, according to what Shchepetkin said above, this psi-mask bug always affected our simulations, even in NON ANALYTICAL GRIDS, because instead of 2, along not straight boundaries, pmask should have been less to generate the right amount of vorticity!
If, instead I am wrong (more likely!), could someone explain to me how ROMS takes into account the land masking and the slipperiness parameter if all the viscosity/diffusivity switches are off? It seems to me, in fact, that this happens only when these switches are on...

Thanks again,
M.

[shchepet]

pmask settings have any effect only if viscosity is nonzero.

The code may be run with or without explicit viscosity, if the upstream-biased advection is
used, because the advection scheme is already dissipative enough on grid scale to avoid
numerical issues.

If, on the other hand, centered advection is selected, some viscosity is required.

In Navier-Stokes equation tangential boundary condition appear only in context of viscous
stress at the boundaries and nowhere else. Numerically they are relevant only if viscosity
is sufficiently large to make viscous boundary layer next to the boundary resolved, that
is grid-scale-based Reynolds number should not be very large. This may or may not be the
case in ocean modeling, and typically in ROMS practice the value for explicit viscosity
coefficient is selected to be small enough so that viscous boundary is unresolved.
Setting of pmask in this case does not make much difference. Instead, code relies
on implicit boundary conditions applied to high-order derivatives via masking of
intermediate fields computed in order to compute advection terms. Qualitatively this
acts as no-slip boundary conditions in sense you see vortex sheet generation off the
boundary, but actually these are highly numerically controlled and should not be viewed
as physical, although they produce physically reasonable results all the time. It reality,
when going close to the shore, one should expect that bottom drag takes over control,
making side boundary conditions essentially irrelevant.

On the other hand, if you want to pose a problem which has well-defined viscous boundary,
for example, flow around cylinder (starting with Stokes highly viscous flow), and you want
to observe regime change with Reynolds number change (RE is based on physical scale
this time, not grid scale), and your cylinder is implemented as approximate area of masked
points, then pmask matter, if you want your results more or less agree with classical (say
transition from statuinary flow to Karmann street occurs at RE=28 more or less. This is,
of course, a different class of problems.