metrics_fill_spherical Subroutine

   subroutine metrics_fill_spherical(this, grid, lon_west, lat_south, &
                                     dlon_deg, dlat_deg, rad_earth)
      !! Spherical lon-lat sector.  For each stagger, geolat/geolon are
      !! evaluated at THAT point's own location; the metric lengths use
      !! the cos of that stagger's own latitude (the consistency trick
      !! that keeps the C-grid metrics compatible, D6):
      !!   dx = rad_earth * cos(lat) * dlon_rad
      !!   dy = rad_earth * dlat_rad
      !!   area = dx * dy   (analytic-derivative form, NOT great-circle).
      !!
      !! Indexing: the first INTERIOR T cell is `(1+nghost, 1+nghost)`,
      !! centred at `(lon_west + (i+i_offset_global-nghost-0.5)*dlon,
      !! lat_south + (j+j_offset_global-nghost-0.5)*dlat)`.  On an
      !! undecomposed grid the offsets are 0 and the formula reduces to
      !! the original single-rank form.  Under MPI decomposition the
      !! offsets shift the local (i,j) to the correct GLOBAL coordinate
      !! so every rank computes the right geolat/geolon.  Corners (Bu)
      !! sit half a cell up/right of their cell centre.  u-faces share
      !! the T latitude, v-faces / corners use the corner latitude.
      !! ALL ghost rows/columns are filled (the formula extends naturally
      !! past the physical sector).
      type(ocean_metrics_t), intent(inout) :: this
      type(hgrid_t), intent(in) :: grid
      real(wp), intent(in) :: lon_west, lat_south, dlon_deg, dlat_deg, rad_earth

      integer :: i, j, nx, ny, ng
      real(wp) :: dlon_rad, dlat_rad, dy_len
      real(wp) :: lat_t, lon_t, lat_b, lon_b

      nx = grid%nx_total
      ny = grid%ny_total
      ng = grid%nghost
      dlon_rad = dlon_deg*DEG2RAD
      dlat_rad = dlat_deg*DEG2RAD
      dy_len = rad_earth*dlat_rad   ! meridional length is lat-independent

      ! ---- T points + u-faces (share the T-row latitude) ----
      do j = 1, ny
         lat_t = lat_south + (real(j + grid%j_offset_global - ng, wp) - 0.5_wp)*dlat_deg
         do i = 1, nx
            lon_t = lon_west + (real(i + grid%i_offset_global - ng, wp) - 0.5_wp)*dlon_deg
            this%geolatT(i, j) = lat_t
            this%geolonT(i, j) = lon_t
            this%dxT(i, j) = rad_earth*cos(lat_t*DEG2RAD)*dlon_rad
            this%dyT(i, j) = dy_len
            this%areaT(i, j) = this%dxT(i, j)*dy_len
         end do
      end do

      ! Cu (u-face): x-face on the T-row latitude.  Cu(i,j) sits on the
      ! west edge of T(i,j); its length uses lat_t (same row).
      !
      ! Uniform-dlon coincidence note: the correct T-to-T dxCu definition is
      ! the sum of the two half-segments straddling the face node (one from
      ! the cell to the west, one from the cell to the east).  On this analytic
      ! generator DLON is constant, so both half-segments are equal and the
      ! two-half-segment sum = R·cos(lat_face)·dlon_rad — which is exactly
      ! what this formula computes (lat_face = lat_t for Cu on the T row).
      ! For a variable-resolution supergrid the two definitions diverge; that
      ! path is corrected in `metrics_fill_from_supergrid` (see comment there).
      do j = 1, ny
         lat_t = lat_south + (real(j + grid%j_offset_global - ng, wp) - 0.5_wp)*dlat_deg
         do i = 1, nx + 1
            this%dxCu(i, j) = rad_earth*cos(lat_t*DEG2RAD)*dlon_rad
            this%dyCu(i, j) = dy_len
            this%dy_cu(i, j) = dy_len
            this%areaCu(i, j) = this%dxCu(i, j)*dy_len
         end do
      end do

      ! Cv (v-face) + Bu (corner): on the corner latitude row.
      ! dyCv uniform-dlat coincidence note: the T-to-T dyCv definition is the
      ! sum of two half-segments straddling the face node row.  On this analytic
      ! generator DLAT is constant so the sum = R·dlat_rad = dy_len — identical
      ! to what is coded here.  Variable-resolution corrected in supergrid reader.
      do j = 1, ny + 1
         lat_b = lat_south + real(j + grid%j_offset_global - ng - 1, wp)*dlat_deg
         do i = 1, nx
            this%dxCv(i, j) = rad_earth*cos(lat_b*DEG2RAD)*dlon_rad
            this%dyCv(i, j) = dy_len
            this%dx_cv(i, j) = this%dxCv(i, j)
            this%areaCv(i, j) = this%dxCv(i, j)*dy_len
         end do
      end do

      do j = 1, ny + 1
         lat_b = lat_south + real(j + grid%j_offset_global - ng - 1, wp)*dlat_deg
         do i = 1, nx + 1
            lon_b = lon_west + real(i + grid%i_offset_global - ng - 1, wp)*dlon_deg
            this%geolatBu(i, j) = lat_b
            this%geolonBu(i, j) = lon_b
            this%dxBu(i, j) = rad_earth*cos(lat_b*DEG2RAD)*dlon_rad
            this%dyBu(i, j) = dy_len
            this%areaBu(i, j) = this%dxBu(i, j)*dy_len
         end do
      end do
   end subroutine metrics_fill_spherical