tripolar_supergrid_arrays Subroutine

   subroutine tripolar_supergrid_arrays(grid, lon_west, lat_south, dlat_deg, rad_earth, &
                                        phi_join, lon_pole, sg_x, sg_y, sg_dx, sg_dy, sg_area)
      !! The in-memory MOM6-style supergrid of the analytic tripolar grid
      !! (node geography, great-circle edge lengths, spherical sub-cell
      !! areas) that `metrics_fill_tripolar` assembles.  Public so a test
      !! can write the very same grid as a mosaic file and check that the
      !! NetCDF reader reproduces the generator's metrics, ghosts included.
      type(hgrid_t), intent(in) :: grid
      real(wp), intent(in) :: lon_west, lat_south, dlat_deg
      real(wp), intent(in) :: rad_earth, phi_join, lon_pole
      real(wp), allocatable, intent(out) :: sg_x(:, :), sg_y(:, :)
         !! Node longitude / latitude (degrees), `(2ni+1, 2nj+1)`.
      real(wp), allocatable, intent(out) :: sg_dx(:, :)
         !! Along-i segment lengths (m), `(2ni, 2nj+1)`.
      real(wp), allocatable, intent(out) :: sg_dy(:, :)
         !! Along-j segment lengths (m), `(2ni+1, 2nj)`.
      real(wp), allocatable, intent(out) :: sg_area(:, :)
         !! Sub-cell areas (m^2), `(2ni, 2nj)`.

      integer :: ni, nj, sg_nxp, sg_nyp, m, n
      real(wp) :: lat_top, dlam, dlat_sg

      ni = grid%nx_phys
      nj = grid%ny_phys
      sg_nxp = 2*ni + 1
      sg_nyp = 2*nj + 1

      ! Supergrid node spacing: half a model cell in each logical direction.
      ! The i-direction must wrap exactly 360 deg of pseudo-longitude over
      ! ni model cells, so dlam per supergrid i-step = 360/(2*ni).
      dlam = 360.0_wp/real(2*ni, wp)
      dlat_sg = dlat_deg*0.5_wp

      ! Top corner latitude of the lon-lat ladder (the cap is mapped between
      ! phi_join and lat_top).  Bu corner j runs 1..nj+1; supergrid corner
      ! row n = 2*nj+1 is the top corner = lat_south + nj*dlat.
      lat_top = lat_south + real(nj, wp)*dlat_deg

      allocate (sg_x(sg_nxp, sg_nyp), source=0.0_wp)
      allocate (sg_y(sg_nxp, sg_nyp), source=0.0_wp)
      allocate (sg_dx(2*ni, sg_nyp), source=0.0_wp)
      allocate (sg_dy(sg_nxp, 2*nj), source=0.0_wp)
      allocate (sg_area(2*ni, 2*nj), source=0.0_wp)

      ! ---- Supergrid node geography (lon-lat below join, bipolar above) ----
      ! Supergrid node (m,n), m=1..2ni+1, n=1..2nj+1.  Corner Bu(1,1) is at
      ! node (1,1) = (lon_west, lat_south).  A T-centre lies at even (m,n).
      ! lon-lat geographic at node (m,n):
      !   lon0 = lon_west + (m-1)*dlam
      !   lat0 = lat_south + (n-1)*dlat_sg
      do n = 1, sg_nyp
         do m = 1, sg_nxp
            call tripolar_node_latlon(m, n, lon_west, lat_south, dlam, dlat_sg, &
                                      phi_join, lat_top, lon_pole, &
                                      sg_y(m, n), sg_x(m, n))
         end do
      end do

      ! ---- Edge lengths by great-circle distance between adjacent nodes ----
      do n = 1, sg_nyp
         do m = 1, 2*ni
            sg_dx(m, n) = great_circle(rad_earth, sg_y(m, n), sg_x(m, n), &
                                       sg_y(m + 1, n), sg_x(m + 1, n))
         end do
      end do
      do n = 1, 2*nj
         do m = 1, sg_nxp
            sg_dy(m, n) = great_circle(rad_earth, sg_y(m, n), sg_x(m, n), &
                                       sg_y(m, n + 1), sg_x(m, n + 1))
         end do
      end do

      ! ---- Sub-cell areas: spherical quad from the four corner nodes ----
      do n = 1, 2*nj
         do m = 1, 2*ni
            sg_area(m, n) = spherical_quad_area(rad_earth, &
                                                sg_y(m, n), sg_x(m, n), &
                                                sg_y(m + 1, n), sg_x(m + 1, n), &
                                                sg_y(m + 1, n + 1), sg_x(m + 1, n + 1), &
                                                sg_y(m, n + 1), sg_x(m, n + 1))
         end do
      end do
   end subroutine tripolar_supergrid_arrays