ocean_diag_register Subroutine

   subroutine ocean_diag_register(this, name, units, fill, n1, n2, n3, &
                                  long_name, standard_name, time_op, dt_out, &
                                  output_vgrid, remap, mask, is_extensive, has_missing)
      !! Register a new diagnostic variable.  Grows the registry via
      !! capacity doubling on overflow.  Buffer allocation depends on
      !! `output_vgrid`:
      !!   * `LAYER` (default): one `output_buffer(n1, n2, n3)` —
      !!     `fill` writes directly into it.
      !!   * `Z_FIXED` (or any non-LAYER target): two buffers —
      !!     `layer_buffer(n1, n2, n3)` for the fill, plus
      !!     `output_buffer(n1, n2, this%nz_out)` for the remapped
      !!     result.  Caller must have configured `nz_out` via
      !!     `set_output_z_levels` first, and bind a `remap` proc.
      !! Caller binds `fill` to a routine that knows how to populate
      !! the layer-native buffer from the state handle.
      !!
      !! `mask` (optional): restricts accumulation to the region
      !! where `mask%weight > 0`.  See `diag_mask_t` builders in
      !! `rdb_ocean_diag_mask`.
      class(ocean_diag_t), intent(inout) :: this
      character(len=*), intent(in) :: name
      character(len=*), intent(in) :: units
      procedure(diag_fill_proc) :: fill
      integer, intent(in) :: n1, n2, n3
      character(len=*), intent(in), optional :: long_name, standard_name
      integer, intent(in), optional :: time_op
      real(wp), intent(in), optional :: dt_out
      integer, intent(in), optional :: output_vgrid
      procedure(diag_remap_proc), optional :: remap
      type(diag_mask_t), intent(in), optional :: mask
      logical, intent(in), optional :: is_extensive
      logical, intent(in), optional :: has_missing
      type(diag_var_t), allocatable :: tmp(:)
      integer :: i, ovgrid, nzout

      if (.not. this%is_init) return

      if (this%nvars == this%nvars_max) then
         allocate (tmp(2*this%nvars_max))
         do i = 1, this%nvars
            tmp(i) = this%vars(i)
         end do
         call move_alloc(tmp, this%vars)
         this%nvars_max = size(this%vars)
      end if

      ovgrid = DIAG_VGRID_LAYER
      if (present(output_vgrid)) ovgrid = output_vgrid

      this%nvars = this%nvars + 1
      associate (v => this%vars(this%nvars))
         v%name = name
         v%units = units
         v%fill => fill
         v%output_vgrid = ovgrid
         if (present(remap)) v%remap => remap
         if (present(long_name)) v%long_name = long_name
         if (present(standard_name)) v%standard_name = standard_name
         if (present(time_op)) v%time_op = time_op
         if (present(dt_out)) v%dt_out = dt_out
         if (present(has_missing)) v%has_missing = has_missing
         if (present(is_extensive)) v%is_extensive = is_extensive

         if (ovgrid == DIAG_VGRID_LAYER) then
            allocate (v%output_buffer(n1, n2, n3), source=0.0_wp)
         else
            allocate (v%layer_buffer(n1, n2, n3), source=0.0_wp)
            select case (ovgrid)
            case (DIAG_VGRID_DENSITY)
               nzout = this%n_rho_out
            case (DIAG_VGRID_SIGMA)
               nzout = this%n_sigma_out
            case (DIAG_VGRID_ZSTAR)
               nzout = this%n_zstar_out
            case default
               nzout = this%nz_out
            end select
            if (nzout <= 0) nzout = n3
            allocate (v%output_buffer(n1, n2, nzout), source=0.0_wp)
         end if

         if (v%time_op /= DIAG_OP_INSTANT) then
            allocate (v%accumulator(size(v%output_buffer, 1), &
                                    size(v%output_buffer, 2), &
                                    size(v%output_buffer, 3)))
            call reset_accumulator(v)
         end if
         v%n_accum = 0
         v%dt_accum = 0.0_wp

         if (present(mask)) then
            allocate (v%mask, source=mask)
         end if
      end associate
   end subroutine ocean_diag_register