ocean_slopes_init Subroutine

   subroutine ocean_slopes_init(this, grid, nz_ml)
      !! Allocate the slope / N² outputs + the vert-fill T/S scratch +
      !! the interface-height buffer.  Default `nz_ml = 1` preserves the
      !! barotropic-only constructor; pass `nz_ml = ms%nz_ml` for the
      !! multilayer driver.  Setup code uses plain host allocation (no
      !! `do concurrent` before `enter_data`).
      class(ocean_slopes_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
      if (nz < 1) nz = 1
      ! Fail loud at configure: vert_fill uses NZ_STACK_MAX-sized column
      ! locals; nz beyond that silently overruns the GPU stack.
      if (nz > NZ_STACK_MAX) then
         error stop "ocean_slopes_init: nz_ml exceeds NZ_STACK_MAX "// &
            "(raise NZ_STACK_MAX in rdb_constants)"
      end if
      this%nx_total = nx
      this%ny_total = ny
      this%nz_ml = nz

      allocate (this%slope_x(nx + 1, ny, nz + 1), source=0.0_wp)
      allocate (this%slope_y(nx, ny + 1, nz + 1), source=0.0_wp)
      allocate (this%n2_u(nx + 1, ny, nz + 1), source=0.0_wp)
      allocate (this%n2_v(nx, ny + 1, nz + 1), source=0.0_wp)
      allocate (this%t_fill(nx, ny, nz), source=0.0_wp)
      allocate (this%s_fill(nx, ny, nz), source=0.0_wp)
      allocate (this%e_int(nx, ny, nz + 1), source=0.0_wp)
      allocate (this%bathy(nx, ny), source=0.0_wp)
      this%bathy_set = .false.
      this%is_init = .true.
   end subroutine ocean_slopes_init