ocean_pressure_force_init Subroutine

   subroutine ocean_pressure_force_init(this, grid, nz_ml)
      !! Allocate the scratch buffers.  Default nz_ml=1 keeps the
      !! barotropic-only path constructible; passing nz_ml sizes them
      !! for the multilayer kernel.
      !!
      !! ALLOCATION GATE (`this%scratch_gated`, see the type docstring):
      !! when `.true.` only the buffers the active `variant` /
      !! `reconstruct_for_pressure` can actually reach are allocated.  The
      !! gate defaults `.false.` so a bare `pgf%init(...)` (every direct
      !! test/benchmark call site) keeps the historical allocate-everything
      !! behaviour.  Each gate's unreachability proof is stated inline.
      class(ocean_pressure_force_t), intent(inout) :: this
      type(hgrid_t), intent(in) :: grid
      integer, intent(in), optional :: nz_ml
      integer :: nx, ny, nz
      logical :: gate, need_p_edge, need_z_centre, need_rho_insitu
      logical :: need_mont_M, need_fv_mom6, need_recon

      nx = grid%nx_total
      ny = grid%ny_total
      nz = 1
      if (present(nz_ml)) nz = nz_ml

      gate = this%scratch_gated
      ! `p_edge`: filled by Pass 1 and read by Pass 2/3.  The MONT, GPRIME and
      ! FV_MOM6 branches of `ocean_pressure_force_compute` `return` or branch
      ! away BEFORE Pass 1, so none of them ever touches it.  MONT builds the
      ! Montgomery potential straight from `rho_layer` + interface heights and
      ! needs no pressure stack at all.  No other module references it.
      need_p_edge = .not. gate .or. &
                    (this%variant == OPGF_VARIANT_FV_LITE .or. &
                     this%variant == OPGF_VARIANT_FV_WRIGHT)
      ! `z_centre`: written by Pass 1b and read by the MONT / FV_LITE /
      ! FV_WRIGHT face passes — every one of those sites sits inside a branch
      ! testing for exactly those three variants.  The Pass-4 grounded-layer
      ! gate reads it too, and that pass runs for FV_MOM6 as well — but ONLY
      ! when `skip_nonoverlap` is on, which is a pure config-time decision
      ! (vcoord + namelist) that the caller latches before `init`.  So FV_MOM6
      ! gets the buffer exactly when something will read it, and an ungated
      ! FV_MOM6 run still pays nothing.  No other module references it.
      need_z_centre = .not. gate .or. &
                      (this%variant == OPGF_VARIANT_MONT .or. &
                       this%variant == OPGF_VARIANT_FV_LITE .or. &
                       this%variant == OPGF_VARIANT_FV_WRIGHT .or. &
                       (this%variant == OPGF_VARIANT_FV_MOM6 .and. &
                        this%skip_nonoverlap))
      ! `mont_M`: written and read by the MONT branch alone.  No other
      ! variant, and no other module, references it.
      need_mont_M = .not. gate .or. (this%variant == OPGF_VARIANT_MONT)
      ! `rho_insitu`: written by the Pass-1 FV_WRIGHT branch (both the EOS
      ! column sweep and its no-tracer fallback) and read by the Pass 2/3
      ! FV_WRIGHT branches only.  No other module references it.
      need_rho_insitu = .not. gate .or. (this%variant == OPGF_VARIANT_FV_WRIGHT)
      ! FV_MOM6 stack: passed only to `compute_fv_mom6_impl` /
      ! `compute_fv_mom6_reconstruct_impl`, both inside the
      ! `if (variant == FV_MOM6)` branch that `return`s.  The one external
      ! reader, `compute_pbce` (rdb_barotropic_coupling), opens with an
      ! `error stop` unless `variant == OPGF_VARIANT_FV_MOM6`.
      need_fv_mom6 = .not. gate .or. (this%variant == OPGF_VARIANT_FV_MOM6)
      ! Reconstruction edge scratch: passed only to
      ! `compute_fv_mom6_reconstruct_impl`, whose call site additionally
      ! requires `reconstruct_for_pressure` (configure_ocean_pgf `error
      ! stop`s if that knob is set with any variant other than FV_MOM6).
      need_recon = .not. gate .or. this%reconstruct_for_pressure

      ! Cell-centred edge stack: (nx, ny, nz+1)
      if (need_p_edge) call this%p_edge%init(nx, ny, nz + 1, "ocean_pgf_p_edge")
      ! Layer-centre z-coordinate: (nx, ny, nz)
      if (need_z_centre) call this%z_centre%init(nx, ny, nz, "ocean_pgf_z_centre")
      ! Montgomery potential at layer centres: (nx, ny, nz)
      if (need_mont_M) call this%mont_M%init(nx, ny, nz, "ocean_pgf_mont_M")
      ! In-situ density at layer centres: (nx, ny, nz)
      if (need_rho_insitu) call this%rho_insitu%init(nx, ny, nz, "ocean_pgf_rho_insitu")
      ! East-face: (nx+1, ny, nz) — same shape as u_face_x_layer.  Every
      ! variant writes these (they ARE the PGF output), so never gated.
      call this%dpdx_face%init(nx + 1, ny, nz, "ocean_pgf_dpdx_face")
      ! North-face: (nx, ny+1, nz)
      call this%dpdy_face%init(nx, ny + 1, nz, "ocean_pgf_dpdy_face")

      ! FV_MOM6 scratch (interface heights, pa stack, per-layer
      ! integrals, per-face horizontal integrals).
      if (need_fv_mom6) then
         call this%e_face%init(nx, ny, nz + 1, "ocean_pgf_fv_mom6_e_face")
         call this%pa%init(nx, ny, nz + 1, "ocean_pgf_fv_mom6_pa")
         call this%intz_dpa%init(nx, ny, nz, "ocean_pgf_fv_mom6_intz_dpa")
         call this%intx_pa%init(nx + 1, ny, nz + 1, "ocean_pgf_fv_mom6_intx_pa")
         call this%inty_pa%init(nx, ny + 1, nz + 1, "ocean_pgf_fv_mom6_inty_pa")
         call this%intx_dpa%init(nx + 1, ny, nz, "ocean_pgf_fv_mom6_intx_dpa")
         call this%inty_dpa%init(nx, ny + 1, nz, "ocean_pgf_fv_mom6_inty_dpa")
         call this%conc_T%init(nx, ny, nz, "ocean_pgf_fv_mom6_conc_T")
         call this%conc_S%init(nx, ny, nz, "ocean_pgf_fv_mom6_conc_S")
      end if

      ! In-layer reconstruction edge scratch (nx, ny, nz).
      if (need_recon) then
         call this%recon_T_t%init(nx, ny, nz, "ocean_pgf_recon_T_t")
         call this%recon_T_b%init(nx, ny, nz, "ocean_pgf_recon_T_b")
         call this%recon_S_t%init(nx, ny, nz, "ocean_pgf_recon_S_t")
         call this%recon_S_b%init(nx, ny, nz, "ocean_pgf_recon_S_b")
      end if

      ! Bathymetry copy.  Default zero = flat bed.  Driver overwrites
      ! via `set_bathymetry` after `state%barotropic%b` is populated.
      allocate (this%b(nx, ny), source=0.0_wp)

      this%is_init = .true.
   end subroutine ocean_pressure_force_init