Handling of horizontal gradients near rough bathymetry

Versions v4.0.x
1

This may be a basic question but I am curious to know how nemo handles the computation of horizontal gradients next to a boundary edge, such as in regions of rough ocean bathymetry. In a routine I am implementing I need to calculate horizontal gradients of buoyancy, which are done like

DO jk = 1, jpkm1
   DO jj = 1, jpjm1
      DO ji = 1, jpim1
         !== gradients of buoyancy on U and V grid ==!
         zbudxup(ji,jj,jk) = r1_e1t(ji,jj) * ( zbu(ji+1,jj,jk) - zbu(ji,jj,jk) )
         zbudyvp(ji,jj,jk) = r1_e2t(ji,jj) * ( zbu(ji,jj+1,jk) - zbu(ji,jj,jk) )
      END DO
   END DO
END DO

where zbu is on the T-point. Several other horizontal gradients are required (relative vorticity, divergence).

For example, what happens when ji+1 is outside the ocean domain due to, say, an ocean ridge. Is zbu(ji+1,jj,jk) equal to NaN or zero? And how does this impact the horizontal gradient? Is there a technique that I need to use to avoid issues at the boundary edges?

2

Over a C-grid, NEMO basically relies on mask arrays (=0/1) to handle such boundary conditions.
In your particular case, such gradients generally refer to fluxes, which are zero across solid boundaries (so that this does not create any source/sink globally. Practically speaking, zbudxup and zbvdxvp would have been multiplied by umask and vmask respectively, which translates into a 0-gradient boundary condition also.
Apply boundary conditions for dynamical conditions can still make use of masking, but this may get more complicated depending on the physical hypothesis you make at the boundary.