ocean_obc_update_reservoirs Subroutine

   subroutine ocean_obc_update_reservoirs(grid, bc, ms, dt)
      !! Evolve per-edge reservoir concentrations `tres` one timestep.
      !! Called after `continuity_tracer_step_split` while
      !! `ms%mass_flux_{x,y}_layer` still hold the stage's wall-face fluxes.
      !! Per open-ish edge (with allocated `tres_*`), per (j|i, k, tracer):
      !!   u_n = sign_edge * mass_flux(wall) / max(h_int, h_min)  (outward normal)
      !! Implicit backward-Euler (Marchesiello et al. 2001):
      !!   c_out = max(0,u_n)*dt/L_out, c_in = max(0,-u_n)*dt/L_in  (0 if L==0)
      !!   tres = (tres + c_out*T_int + c_in*T_data)/(1 + c_out + c_in)
      !! Degenerate L_out==0 & u_n>0 ⇒ tres = T_int (L_in==0 & u_n<0 ⇒ T_data).
      !! T_int = hTr(int)/max(h(int),h_min); T_data = clamped_tracer(it).
      !! No-op when both length scales are zero (default path).
      type(hgrid_t), intent(in) :: grid
      type(ocean_bc_state_t), intent(inout) :: bc
      type(multilayer_state_t), intent(in) :: ms
      real(wp), intent(in) :: dt

      logical :: do_res
      integer :: nxt, nyt, nz, nx, ny
      integer :: i_w, i_e, j_s, j_n
      integer :: it, n_tr

      do_res = (bc%res_lscale_out > 0.0_wp .or. bc%res_lscale_in > 0.0_wp)
      if (.not. do_res) return

      nxt = grid%nx_total
      nyt = grid%ny_total
      nz = ms%nz_ml
      nx = grid%nx_phys
      ny = grid%ny_phys
      i_w = grid%nghost + 1
      i_e = grid%nghost + nx
      j_s = grid%nghost + 1
      j_n = grid%nghost + ny

      n_tr = bc%n_tracers
      if (n_tr <= 0) return
      if (.not. allocated(ms%tracers)) return

      ! Per-tracer outer shim: pass full explicit-shape arrays to kernels.
      do it = 1, n_tr
         if (.not. allocated(ms%tracers(it)%hTr)) cycle

         ! ---- West ----
         ! MPI-seam neutralisation (O0): dispatch here is on allocated(tres_*),
         ! not bc_type, so the remap trick cannot be used — explicit has_* guard.
         if (allocated(bc%tres_west) .and. bc%has_west) then
            call update_reservoir_zonal_west( &
               ms%tracers(it)%hTr, ms%h_layer, ms%mass_flux_x_layer, &
               bc%tres_west, &
               get_clamped_tracer(bc%west%clamped_tracer, it), &
               nxt, nyt, nz, i_w, j_s, j_n, &
               bc%res_lscale_out, bc%res_lscale_in, dt, it, n_tr)
         end if

         ! ---- East ----
         if (allocated(bc%tres_east) .and. bc%has_east) then
            call update_reservoir_zonal_east( &
               ms%tracers(it)%hTr, ms%h_layer, ms%mass_flux_x_layer, &
               bc%tres_east, &
               get_clamped_tracer(bc%east%clamped_tracer, it), &
               nxt, nyt, nz, i_e, j_s, j_n, &
               bc%res_lscale_out, bc%res_lscale_in, dt, it, n_tr)
         end if

         ! ---- South ----
         if (allocated(bc%tres_south) .and. bc%has_south) then
            call update_reservoir_meridional_south( &
               ms%tracers(it)%hTr, ms%h_layer, ms%mass_flux_y_layer, &
               bc%tres_south, &
               get_clamped_tracer(bc%south%clamped_tracer, it), &
               nxt, nyt, nz, j_s, i_w, i_e, &
               bc%res_lscale_out, bc%res_lscale_in, dt, it, n_tr)
         end if

         ! ---- North ----
         if (allocated(bc%tres_north) .and. bc%has_north) then
            call update_reservoir_meridional_north( &
               ms%tracers(it)%hTr, ms%h_layer, ms%mass_flux_y_layer, &
               bc%tres_north, &
               get_clamped_tracer(bc%north%clamped_tracer, it), &
               nxt, nyt, nz, j_n, i_w, i_e, &
               bc%res_lscale_out, bc%res_lscale_in, dt, it, n_tr)
         end if

      end do
   end subroutine ocean_obc_update_reservoirs