continuity_init Subroutine

   subroutine continuity_init(this, grid, nz_ml)
      !! Allocate the 4 face-reconstruction scratch buffers sized at
      !! (nx_face, ny_face, nz).  Default nz=1 covers the barotropic
      !! kernel; passing `nz_ml` sizes them for the multilayer
      !! kernel without forcing a separate init routine.  Ocean
      !! init passes `state%multilayer%nz_ml` when the multilayer
      !! state is in play.
      class(continuity_t), intent(inout) :: this
      type(hgrid_t), intent(in) :: grid
      integer, intent(in), optional :: nz_ml
      integer :: nx, ny, nz

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

      ! East-face shapes: (nx+1, ny, nz)
      call this%h_face_left_x%init(nx + 1, ny, nz, "continuity_h_face_left_x")
      call this%h_face_right_x%init(nx + 1, ny, nz, "continuity_h_face_right_x")
      ! North-face shapes: (nx, ny+1, nz)
      call this%h_face_left_y%init(nx, ny + 1, nz, "continuity_h_face_left_y")
      call this%h_face_right_y%init(nx, ny + 1, nz, "continuity_h_face_right_y")
      ! Cell-centred conservative min-thickness target scratch: (nx, ny, nz)
      call this%mt_h_new%init(nx, ny, nz, "continuity_mt_h_new")
      call this%mt_grounded%init(nx, ny, 1, "continuity_mt_grounded")
      ! P2 positive-definite outflux-limiter θ scratch (unconditional).
      call this%pd_theta%init(nx, ny, nz, "continuity_pd_theta")

      this%t_dyn_rel_adv = 0.0_wp
      ! Phase 2 flux accumulators + (6b) windowed-drain workspace — 13
      ! 3D arrays (~3.7 GiB at 600²x100) consumed ONLY by the
      ! TR_MODE_ACCUMULATE path (`dt_tracer_advect_ratio > 1`), so their
      ! allocation is gated on `windowed_advection` (latched from cfg
      ! BEFORE init by `ocean_state_init_from_config`, the same
      ! conditional-allocation contract as the default-off closures).
      ! Default .false. — direct `ct%init(...)` call sites that DO drive
      ! the drain opt in explicitly.  Plain host allocation + zero (no `do
      ! concurrent` before enter_data — that would force per-loop H<->D
      ! round-trips; setup code is host-side); mapped onto the device in
      ! continuity_enter_data_impl behind its existing `allocated()`
      ! guards.
      if (this%windowed_advection) then
         allocate (this%uhtr(nx + 1, ny, nz))
         allocate (this%vhtr(nx, ny + 1, nz))
         this%uhtr = 0.0_wp
         this%vhtr = 0.0_wp
         allocate (this%hprev_work(nx, ny, nz), source=0.0_wp)
         allocate (this%h_win_start(nx, ny, nz), source=0.0_wp)
         allocate (this%uhr_x(nx + 1, ny, nz), source=0.0_wp)
         allocate (this%uhr_y(nx, ny + 1, nz), source=0.0_wp)
         allocate (this%uhh_x(nx + 1, ny, nz), source=0.0_wp)
         allocate (this%uhh_y(nx, ny + 1, nz), source=0.0_wp)
         allocate (this%tr_flux_x(nx + 1, ny, nz), source=0.0_wp)
         allocate (this%tr_flux_y(nx, ny + 1, nz), source=0.0_wp)
         allocate (this%tr_work(nx, ny, nz), source=0.0_wp)
         allocate (this%pal(nx, ny, nz), source=0.0_wp)
         allocate (this%par(nx, ny, nz), source=0.0_wp)
         allocate (this%pa6(nx, ny, nz), source=0.0_wp)
      end if

      this%is_init = .true.
   end subroutine continuity_init