continuity_apply_meridional Subroutine

private pure subroutine continuity_apply_meridional(grid, metrics, ms, dt, h_min)

Apply the meridional (y-flux) thickness update on top of the zonally-updated state: h(i,j,k) ← h(i,j,k) - dt · (Φy(i,j+1,k) - Φy(i,j,k)) · iareaT Adds the y-divergence to flux_h_layer so the field ends the split step holding the total horizontal divergence that the vertical-advection kernel consumes (w_interface(k+1) = w(k) - flux_h_layer(k)). Φy carries dx_cv; iareaT = inv_dy on uniform metrics.

Optional h_min (m): when > 0, applies max(h_new, h_min) on the h-update (Phase-1 Lagrangian floor). mass_budget_continuity records the RAW divergence regardless — the floor injection shows up in the mass Error diagnostic (R2). 0 or absent ⇒ bit-identical.

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in), optional :: h_min

Minimum-thickness floor (m). 0 or absent ⇒ off ⇒ bit-identical.

Called by

proc~~continuity_apply_meridional~~CalledByGraph proc~continuity_apply_meridional continuity_apply_meridional proc~continuity_step_split continuity_step_split proc~continuity_step_split->proc~continuity_apply_meridional proc~continuity_tracer_step_split continuity_tracer_step_split proc~continuity_tracer_step_split->proc~continuity_apply_meridional proc~run_continuity_chain run_continuity_chain proc~run_continuity_chain->proc~continuity_tracer_step_split proc~run_stage run_stage proc~run_stage->proc~continuity_tracer_step_split proc~ocean_dyn_step ocean_dyn_step proc~ocean_dyn_step->proc~run_stage proc~run_stage_split run_stage_split proc~run_stage_split->proc~run_continuity_chain proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step proc~ocean_dyn_step_split ocean_dyn_step_split proc~engine_step->proc~ocean_dyn_step_split proc~ocean_dyn_step_split->proc~run_stage_split proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~engine_step proc~rdb_ocean_step rdb_ocean_step proc~rdb_ocean_step->proc~engine_step

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: div_y
real(kind=wp), private :: h_min_use
integer, private :: i
integer, private :: j
integer, private :: k
integer, private :: nx
integer, private :: ny
integer, private :: nz

Source Code

   pure subroutine continuity_apply_meridional(grid, metrics, ms, dt, h_min)
      !! Apply the meridional (y-flux) thickness update on top of
      !! the zonally-updated state:
      !!   h(i,j,k) ← h(i,j,k) - dt · (Φy(i,j+1,k) - Φy(i,j,k)) · iareaT
      !! Adds the y-divergence to `flux_h_layer` so the field ends
      !! the split step holding the *total* horizontal divergence
      !! that the vertical-advection kernel consumes
      !! (`w_interface(k+1) = w(k) - flux_h_layer(k)`).  Φy carries
      !! `dx_cv`; `iareaT` = `inv_dy` on uniform metrics.
      !!
      !! Optional `h_min` (m): when > 0, applies max(h_new, h_min) on the
      !! h-update (Phase-1 Lagrangian floor). `mass_budget_continuity` records
      !! the RAW divergence regardless — the floor injection shows up in the
      !! mass Error diagnostic (R2). 0 or absent ⇒ bit-identical.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: dt
      real(wp), intent(in), optional :: h_min
         !! Minimum-thickness floor (m). 0 or absent ⇒ off ⇒ bit-identical.

      integer :: i, j, k, nx, ny, nz
      real(wp) :: div_y, h_min_use

      h_min_use = 0.0_wp
      if (present(h_min)) h_min_use = h_min

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml

      if (h_min_use > 0.0_wp) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx)
            div_y = (ms%mass_flux_y_layer(i, j + 1, k) - &
                     ms%mass_flux_y_layer(i, j, k))*metrics%iareaT(i, j)
            ms%flux_h_layer(i, j, k) = ms%flux_h_layer(i, j, k) + div_y
            ms%h_layer(i, j, k) = max(ms%h_layer(i, j, k) - dt*div_y, h_min_use)
            ms%mass_budget_continuity(i, j, k) = &
               ms%mass_budget_continuity(i, j, k) - dt*div_y
         end do
      else
         do concurrent(k=1:nz, j=1:ny, i=1:nx)
            div_y = (ms%mass_flux_y_layer(i, j + 1, k) - &
                     ms%mass_flux_y_layer(i, j, k))*metrics%iareaT(i, j)
            ms%flux_h_layer(i, j, k) = ms%flux_h_layer(i, j, k) + div_y
            ms%h_layer(i, j, k) = ms%h_layer(i, j, k) - dt*div_y
            ms%mass_budget_continuity(i, j, k) = &
               ms%mass_budget_continuity(i, j, k) - dt*div_y
         end do
      end if
   end subroutine continuity_apply_meridional