metrics_assemble_from_supergrid_arrays Subroutine

   subroutine metrics_assemble_from_supergrid_arrays(this, grid, &
                                                     sg_x, sg_y, sg_dx, sg_dy, sg_area, &
                                                     periodic_x, sg_angle_dx)
      !! Fill all model metric arrays from an in-memory MOM6-style
      !! supergrid (2x-refined corner geography + edge segments +
      !! sub-cell areas), using the even/odd index sums.  This is the
      !! battle-tested assembly path the NetCDF reader used inline; the
      !! tripolar generator builds the supergrid analytically and feeds it
      !! here so tripolar metrics flow through identical index logic.
      !!
      !! Supergrid index convention (1-based, node (1,1) = SW corner of
      !! the physical domain; 2*ni+1 nodes in i, 2*nj+1 in j):
      !!     node (2i-1, 2j-1): SW corner of T(i,j)  -- ODD/ODD = Bu
      !!     node (2i,   2j  ): T-cell centre          -- EVEN/EVEN = T
      !!     node (2i-1, 2j  ): west  face of T(i,j)    -- ODD/EVEN = Cu
      !!     node (2i,   2j-1): south face of T(i,j)    -- EVEN/ODD = Cv
      !!   sg_dx(m,n): along-i segment node (m,n)->(m+1,n), shape (2ni,2nj+1).
      !!   sg_dy(m,n): along-j segment node (m,n)->(m,n+1), shape (2ni+1,2nj).
      !!   sg_area(m,n): sub-cell area SW at node (m,n), shape (2ni,2nj).
      !! Boundary (Cu i=1/ni+1, Cv j=1/nj+1, Bu edges) by extrapolation —
      !! EXCEPT the east-west seam when `periodic_x`: there the face lies
      !! between the last and the first column, and its along-i span is
      !! the real one, `sg_dx(2ni) + sg_dx(1)`, not a copy of the
      !! neighbouring face's.
      !! Ghost rows by constant extrapolation of the nearest physical value
      !! (periodic / fold ghost fill is `metrics_fold_periodic_ghosts`).
      type(ocean_metrics_t), intent(inout) :: this
      type(hgrid_t), intent(in) :: grid
      real(wp), intent(in) :: sg_x(:, :), sg_y(:, :)
      real(wp), intent(in) :: sg_dx(:, :), sg_dy(:, :), sg_area(:, :)
      logical, intent(in), optional :: periodic_x
         !! The i-direction is periodic (node column `2ni+1` IS column 1):
         !! build the seam-face Cu / Bu spans across the seam.  Absent ⇒
         !! `.false.` (extrapolated, byte-identical to the previous path).
      real(wp), intent(in), optional :: sg_angle_dx(:, :)
         !! Supergrid-node grid rotation in DEGREES, `(2ni+1, 2nj+1)` (the
         !! mosaic's `angle_dx`); the T value is node `(2i, 2j)`.  Absent ⇒
         !! `angle_dx` stays zero (the axes are taken as east/north).

      integer :: ni, nj, ng, sg_nxp, sg_nyp
      integer :: i, j, si, sj, si1, sj1
      logical :: per_x

      per_x = .false.
      if (present(periodic_x)) per_x = periodic_x
      ni = grid%nx_phys
      nj = grid%ny_phys
      ng = grid%nghost
      sg_nxp = 2*ni + 1
      sg_nyp = 2*nj + 1

      ! ---- T metrics (all physical i,j) ----
      do j = 1, nj
         sj = 2*j       ! T-centre sg j-index (even)
         sj1 = 2*j - 1  ! lower corner / SW sg j-index (odd)
         do i = 1, ni
            si = 2*i       ! T-centre sg i-index (even)
            si1 = 2*i - 1  ! left corner / SW sg i-index (odd)
            this%geolatT(ng + i, ng + j) = sg_y(si, sj)
            this%geolonT(ng + i, ng + j) = sg_x(si, sj)
            this%dxT(ng + i, ng + j) = sg_dx(si1, sj) + sg_dx(si, sj)
            this%dyT(ng + i, ng + j) = sg_dy(si, sj1) + sg_dy(si, sj)
            this%areaT(ng + i, ng + j) = sg_area(si1, sj1) + sg_area(si, sj1) + &
                                         sg_area(si1, sj) + sg_area(si, sj)
         end do
      end do

      ! ---- Cu metrics (u-face: i ∈ [1,ni+1], j ∈ [1,nj]) ----
      do j = 1, nj
         sj = 2*j
         sj1 = 2*j - 1
         do i = 1, ni + 1
            si = min(max(2*i - 1, 1), sg_nxp)  ! col of Cu face node
            if (i >= 2 .and. i <= ni) then
               this%dxCu(ng + i, ng + j) = sg_dx(2*i - 2, sj) + sg_dx(2*i - 1, sj)
            else
               this%dxCu(ng + i, ng + j) = 0.0_wp
            end if
            this%dyCu(ng + i, ng + j) = sg_dy(si, sj1) + sg_dy(si, sj)
            this%areaCu(ng + i, ng + j) = this%dxCu(ng + i, ng + j)* &
                                          this%dyCu(ng + i, ng + j)
            this%dy_cu(ng + i, ng + j) = this%dyCu(ng + i, ng + j)
         end do
      end do
      do j = 1, nj
         this%dxCu(ng + 1, ng + j) = this%dxCu(ng + 2, ng + j)
         this%areaCu(ng + 1, ng + j) = this%dxCu(ng + 1, ng + j)* &
                                       this%dyCu(ng + 1, ng + j)
         this%dxCu(ng + ni + 1, ng + j) = this%dxCu(ng + ni, ng + j)
         this%areaCu(ng + ni + 1, ng + j) = this%dxCu(ng + ni + 1, ng + j)* &
                                            this%dyCu(ng + ni + 1, ng + j)
      end do

      ! ---- Cv metrics (v-face: i ∈ [1,ni], j ∈ [1,nj+1]) ----
      do i = 1, ni
         si = 2*i
         si1 = 2*i - 1
         do j = 1, nj + 1
            sj = min(max(2*j - 1, 1), sg_nyp)  ! row of Cv face node
            this%dxCv(ng + i, ng + j) = sg_dx(si1, sj) + sg_dx(si, sj)
            if (j >= 2 .and. j <= nj) then
               this%dyCv(ng + i, ng + j) = sg_dy(si, 2*j - 2) + sg_dy(si, 2*j - 1)
            else
               this%dyCv(ng + i, ng + j) = 0.0_wp
            end if
            this%areaCv(ng + i, ng + j) = this%dxCv(ng + i, ng + j)* &
                                          this%dyCv(ng + i, ng + j)
            this%dx_cv(ng + i, ng + j) = this%dxCv(ng + i, ng + j)
         end do
      end do
      do i = 1, ni
         this%dyCv(ng + i, ng + 1) = this%dyCv(ng + i, ng + 2)
         this%areaCv(ng + i, ng + 1) = this%dxCv(ng + i, ng + 1)* &
                                       this%dyCv(ng + i, ng + 1)
         this%dx_cv(ng + i, ng + 1) = this%dxCv(ng + i, ng + 1)
         this%dyCv(ng + i, ng + nj + 1) = this%dyCv(ng + i, ng + nj)
         this%areaCv(ng + i, ng + nj + 1) = this%dxCv(ng + i, ng + nj + 1)* &
                                            this%dyCv(ng + i, ng + nj + 1)
         this%dx_cv(ng + i, ng + nj + 1) = this%dxCv(ng + i, ng + nj + 1)
      end do

      ! ---- Bu metrics (corner: i ∈ [1,ni+1], j ∈ [1,nj+1]) ----
      do j = 1, nj + 1
         sj = 2*j - 1   ! ODD: sg j-index of Bu corner
         do i = 1, ni + 1
            si = 2*i - 1  ! ODD: sg i-index of Bu corner
            this%geolatBu(ng + i, ng + j) = sg_y(si, sj)
            this%geolonBu(ng + i, ng + j) = sg_x(si, sj)
            if (i >= 2 .and. i <= ni) then
               this%dxBu(ng + i, ng + j) = sg_dx(2*i - 2, sj) + sg_dx(2*i - 1, sj)
            else
               this%dxBu(ng + i, ng + j) = 0.0_wp
            end if
            if (j >= 2 .and. j <= nj) then
               this%dyBu(ng + i, ng + j) = sg_dy(si, 2*j - 2) + sg_dy(si, 2*j - 1)
            else
               this%dyBu(ng + i, ng + j) = 0.0_wp
            end if
            this%areaBu(ng + i, ng + j) = this%dxBu(ng + i, ng + j)* &
                                          this%dyBu(ng + i, ng + j)
         end do
      end do
      do j = 1, nj + 1
         this%dxBu(ng + 1, ng + j) = this%dxBu(ng + 2, ng + j)
         this%areaBu(ng + 1, ng + j) = this%areaBu(ng + 2, ng + j)
         this%dxBu(ng + ni + 1, ng + j) = this%dxBu(ng + ni, ng + j)
         this%areaBu(ng + ni + 1, ng + j) = this%areaBu(ng + ni, ng + j)
      end do
      do i = 1, ni + 1
         this%dyBu(ng + i, ng + 1) = this%dyBu(ng + i, ng + 2)
         if (i > 1) then
            this%areaBu(ng + i, ng + 1) = this%dxBu(ng + i, ng + 1)* &
                                          this%dyBu(ng + i, ng + 1)
         end if
         this%dyBu(ng + i, ng + nj + 1) = this%dyBu(ng + i, ng + nj)
         this%areaBu(ng + i, ng + nj + 1) = this%dxBu(ng + i, ng + nj + 1)* &
                                            this%dyBu(ng + i, ng + nj + 1)
      end do
      this%areaBu(ng + 1, ng + 1) = this%areaBu(ng + 2, ng + 2)

      ! ---- Periodic east-west seam: the faces at i=1 and i=ni+1 are ONE
      ! face, between column ni and column 1.  Its along-i span is the two
      ! half-cell segments either side of it, which on a periodic grid are
      ! `sg_dx(2ni, .)` (east half of column ni) and `sg_dx(1, .)` (west
      ! half of column 1).  The across-face spans (dyCu, dyBu) already come
      ! from the seam node column itself and need nothing.
      !
      ! The seam CORNER area is the mean of the two corners either side of
      ! the seam (columns 2 and ni), the symmetric form of the one-sided
      ! copy the extrapolated edge uses — NOT `dxBu*dyBu`: on a bipolar cap
      ! whose pole meridian is the seam, the seam node column collapses onto
      ! the pole (`dyBu = 0`), and a zero corner area there is a zero
      ! vorticity cell the Coriolis/viscosity stencils cannot live with
      ! (measured: the analytic-tripolar cross-seam test blows up).
      if (per_x) then
         do j = 1, nj
            sj = 2*j
            this%dxCu(ng + 1, ng + j) = sg_dx(2*ni, sj) + sg_dx(1, sj)
            this%dxCu(ng + ni + 1, ng + j) = this%dxCu(ng + 1, ng + j)
            this%areaCu(ng + 1, ng + j) = this%dxCu(ng + 1, ng + j)*this%dyCu(ng + 1, ng + j)
            this%areaCu(ng + ni + 1, ng + j) = this%dxCu(ng + ni + 1, ng + j)* &
                                               this%dyCu(ng + ni + 1, ng + j)
         end do
         do j = 1, nj + 1
            sj = 2*j - 1
            this%dxBu(ng + 1, ng + j) = sg_dx(2*ni, sj) + sg_dx(1, sj)
            this%dxBu(ng + ni + 1, ng + j) = this%dxBu(ng + 1, ng + j)
            this%areaBu(ng + 1, ng + j) = 0.5_wp*(this%areaBu(ng + 2, ng + j) + &
                                                  this%areaBu(ng + ni, ng + j))
            this%areaBu(ng + ni + 1, ng + j) = this%areaBu(ng + 1, ng + j)
         end do
      end if

      ! ---- Grid rotation at T points (degrees in the mosaic -> radians) ----
      if (present(sg_angle_dx)) then
         do j = 1, nj
            do i = 1, ni
               this%angle_dx(ng + i, ng + j) = sg_angle_dx(2*i, 2*j)*DEG2RAD
            end do
         end do
      end if

      ! ---- Ghost extrapolation: constant copy from nearest physical cell ----
      call supergrid_ghost_fill_2d(this%dxT, grid)
      call supergrid_ghost_fill_2d(this%dyT, grid)
      call supergrid_ghost_fill_2d(this%areaT, grid)
      call supergrid_ghost_fill_2d(this%geolatT, grid)
      call supergrid_ghost_fill_2d(this%geolonT, grid)
      call supergrid_ghost_fill_2d(this%angle_dx, grid)
      call supergrid_ghost_fill_cu(this%dxCu, grid)
      call supergrid_ghost_fill_cu(this%dyCu, grid)
      call supergrid_ghost_fill_cu(this%areaCu, grid)
      call supergrid_ghost_fill_cu(this%dy_cu, grid)
      call supergrid_ghost_fill_cv(this%dxCv, grid)
      call supergrid_ghost_fill_cv(this%dyCv, grid)
      call supergrid_ghost_fill_cv(this%areaCv, grid)
      call supergrid_ghost_fill_cv(this%dx_cv, grid)
      call supergrid_ghost_fill_bu(this%dxBu, grid)
      call supergrid_ghost_fill_bu(this%dyBu, grid)
      call supergrid_ghost_fill_bu(this%areaBu, grid)
      call supergrid_ghost_fill_bu(this%geolatBu, grid)
      call supergrid_ghost_fill_bu(this%geolonBu, grid)
   end subroutine metrics_assemble_from_supergrid_arrays