gm_compute_impl Subroutine

   subroutine gm_compute_impl(nx, ny, nz, dt, khth, khth_max_cfl, slope_max, &
                              rho0, use_ext, use_open, dy_cu, dx_cv, idxCu, &
                              idyCv, idyCu, idxCv, areaT, wet_u, wet_v, h_layer, &
                              bathy, slope_x, slope_y, &
                              n2_u, n2_v, khth_ext_u, khth_ext_v, open_u, open_v, &
                              khth_u, khth_v, uhD, vhD, gm_src)
      !! Flat-impl GM kernel.  Passes: CFL-clamp the 2D face KhTh, the
      !! u-face and v-face column recurrences into `uhD`/`vhD`, the
      !! `gm_src` PE release.  Each face's column sweep is serial in k (the
      !! `uhtot` recurrence) but faces parallelise over (i,j).
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in) :: dt, khth, khth_max_cfl, slope_max, rho0
      logical, intent(in) :: use_ext
      logical, intent(in) :: use_open
         !! z-level closed faces active (`metrics%use_closed_faces`).
         !! `.false.` ⇒ `open_u`/`open_v` are never indexed.
      real(wp), intent(in) :: dy_cu(nx + 1, ny)
      real(wp), intent(in) :: dx_cv(nx, ny + 1)
      real(wp), intent(in) :: idxCu(nx + 1, ny)
      real(wp), intent(in) :: idyCv(nx, ny + 1)
      real(wp), intent(in) :: idyCu(nx + 1, ny)
      real(wp), intent(in) :: idxCv(nx, ny + 1)
      real(wp), intent(in) :: areaT(nx, ny)
      real(wp), intent(in) :: wet_u(nx + 1, ny)
      real(wp), intent(in) :: wet_v(nx, ny + 1)
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(in) :: bathy(nx, ny)
         !! Bed depth `D` (m, positive down; `slopes%bathy`): the bed of each
         !! column is at `z = −D` for the bottom-blocking limiter.
      real(wp), intent(in) :: slope_x(nx + 1, ny, nz + 1)
      real(wp), intent(in) :: slope_y(nx, ny + 1, nz + 1)
      real(wp), intent(in) :: n2_u(nx + 1, ny, nz + 1)
      real(wp), intent(in) :: n2_v(nx, ny + 1, nz + 1)
      real(wp), intent(in) :: khth_ext_u(nx + 1, ny)
      real(wp), intent(in) :: khth_ext_v(nx, ny + 1)
      real(wp), intent(in) :: open_u(nx + 1, ny, nz)
         !! Per-layer 0/1 u-face open mask (`metrics%open_u`), read only
         !! when `use_open`; an inert device-present stand-in otherwise.
      real(wp), intent(in) :: open_v(nx, ny + 1, nz)
         !! v-face twin of `open_u`.
      real(wp), intent(out) :: khth_u(nx + 1, ny)
      real(wp), intent(out) :: khth_v(nx, ny + 1)
      real(wp), intent(out) :: uhD(nx + 1, ny, nz)
      real(wp), intent(out) :: vhD(nx, ny + 1, nz)
      real(wp), intent(out) :: gm_src(nx, ny)

      real(wp) :: i_smax2, i4dt
      integer :: i, j, k

      i_smax2 = 1.0_wp/(slope_max*slope_max)
      i4dt = 1.0_wp/(4.0_wp*dt)

      ! ---- 1. CFL-clamp the constant KhTh onto the 2D face fields, using
      ! native face metrics (idxCu/idyCu at u-faces, idxCv/idyCv at
      ! v-faces) — exact on curvilinear grids.  Wall faces -> 0.
      call gm_clamp_khth(nx, ny, dt, khth, khth_max_cfl, use_ext, &
                         khth_ext_u, khth_ext_v, idxCu, idyCu, &
                         idxCv, idyCv, wet_u, wet_v, khth_u, khth_v)

      ! ---- 2. u-face bolus transport (column recurrence).
      do concurrent(k=1:nz, j=1:ny)
         uhD(1, j, k) = 0.0_wp
         uhD(nx + 1, j, k) = 0.0_wp
      end do
      call gm_column_x(nx, ny, nz, i_smax2, i4dt, use_open, &
                       dy_cu, areaT, h_layer, bathy, slope_x, khth_u, open_u, uhD)

      ! ---- 3. v-face bolus transport (mirror).
      do concurrent(k=1:nz, i=1:nx)
         vhD(i, 1, k) = 0.0_wp
         vhD(i, ny + 1, k) = 0.0_wp
      end do
      call gm_column_y(nx, ny, nz, i_smax2, i4dt, use_open, &
                       dx_cv, areaT, h_layer, bathy, slope_y, khth_v, open_v, vhD)

      ! ---- 4. gm_src PE release (cell centres).
      call gm_pe_release(nx, ny, nz, slope_max, rho0, h_layer, &
                         slope_x, slope_y, n2_u, n2_v, khth_u, khth_v, gm_src)
   end subroutine gm_compute_impl