metrics_fill_coriolis Subroutine

   subroutine metrics_fill_coriolis(this, scheme, f_0, beta, y_ref, omega, &
                                    grid, f_corner, f_centre)
      !! Fill a corner array AND a centre array with the Coriolis
      !! parameter, from one of two schemes (D7).  Does NOT touch any
      !! existing fill sites in coriolis_adv / EPBL / kappa-shear (that
      !! re-routing is M2d); this routine just exists + is tested.
      !!
      !!   `beta_plane` : f = f_0 + beta*(y - y_ref), y from the
      !!      GLOBAL CARTESIAN coordinate.  On an undecomposed grid
      !!      (`grid%j_offset_global == 0`) the corner fill is
      !!      BIT-IDENTICAL to `coriolis_adv_set_beta_plane`
      !!      (y = (j-1-nghost)*dy, raw f) and the centre fill is
      !!      BIT-IDENTICAL to the EPBL / kappa-shear `set_f_centre`
      !!      (y = (j-nghost-0.5)*dy, abs(f)).  Under MPI y-decomposition
      !!      `grid%j_offset_global` shifts the local index to the GLOBAL
      !!      row so each rank's beta-plane y is correct.
      !!   `planetary` : f = 2*omega*sin(geolat) at the respective
      !!      stagger — uses `this%geolatBu` (corner) and `this%geolatT`
      !!      (centre), so a spherical generator must have run first.
      !!
      !! `f_corner` is shaped `(nx+1,ny+1)` (like the metric corner
      !! arrays / `coriolis_adv%f_corner`); `f_centre` is `(nx,ny)`.
      type(ocean_metrics_t), intent(in) :: this
      integer, intent(in) :: scheme
      real(wp), intent(in) :: f_0, beta, y_ref, omega
      type(hgrid_t), intent(in) :: grid
      real(wp), intent(out) :: f_corner(:, :)
         !! Coriolis at C-grid corners (1/s).
      real(wp), intent(out) :: f_centre(:, :)
         !! |Coriolis| at cell centres (1/s).

      integer :: i, j, ng
      real(wp) :: y

      ng = grid%nghost

      select case (scheme)
      case (CORIOLIS_SCHEME_PLANETARY)
         ! 2*omega*sin(geolat) at the respective stagger.
         do j = 1, size(f_corner, 2)
            do i = 1, size(f_corner, 1)
               f_corner(i, j) = 2.0_wp*omega*sin(this%geolatBu(i, j)*DEG2RAD)
            end do
         end do
         do j = 1, size(f_centre, 2)
            do i = 1, size(f_centre, 1)
               f_centre(i, j) = abs(2.0_wp*omega*sin(this%geolatT(i, j)*DEG2RAD))
            end do
         end do
      case default   ! CORIOLIS_SCHEME_BETA_PLANE
         ! Corner: BIT-IDENTICAL to coriolis_adv_set_beta_plane on an
         ! undecomposed grid (j_offset_global == 0).  Under MPI y-split
         ! the global row index is j + j_offset_global, giving the correct
         ! physical y coordinate on every rank.
         do j = 1, size(f_corner, 2)
            y = real(j + grid%j_offset_global - 1 - ng, wp)*grid%dy
            do i = 1, size(f_corner, 1)
               f_corner(i, j) = f_0 + beta*(y - y_ref)
            end do
         end do
         ! Centre: BIT-IDENTICAL to EPBL / kappa-shear set_f_centre on an
         ! undecomposed grid (j_offset_global == 0).
         do j = 1, size(f_centre, 2)
            y = (real(j + grid%j_offset_global - ng, wp) - 0.5_wp)*grid%dy
            do i = 1, size(f_centre, 1)
               f_centre(i, j) = abs(f_0 + beta*(y - y_ref))
            end do
         end do
      end select
   end subroutine metrics_fill_coriolis