compute_w_from_continuity Subroutine

public pure subroutine compute_w_from_continuity(grid, this, ms)

Fill ms%w_interface by integrating the horizontal- continuity residual upward from the bed. Eulerian z:

w_interface(:, :, 1) = 0 (bed BC) w_interface(:, :, k+1) = w_interface(:, :, k) - flux_h_layer(k)

With this w, ∂h/∂t = -horizontal_div + (w(k) - w(k+1)) = 0, so layer thicknesses stay at their initial Eulerian z positions. continuity_compute_fluxes must run first — flux_h_layer is the input.

Recurrence is data-dependent in k, so the outer loop is serial in k inside each column. Columns are parallelised over (j, i).

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_vertical_advection_t), intent(in) :: this
type(multilayer_state_t), intent(inout) :: ms

Called by

proc~~compute_w_from_continuity~~CalledByGraph proc~compute_w_from_continuity compute_w_from_continuity proc~run_continuity_chain run_continuity_chain proc~run_continuity_chain->proc~compute_w_from_continuity proc~run_gm_step run_gm_step proc~run_gm_step->proc~compute_w_from_continuity proc~run_stage run_stage proc~run_stage->proc~compute_w_from_continuity proc~ocean_dyn_step ocean_dyn_step proc~ocean_dyn_step->proc~run_stage proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~run_gm_step proc~run_stage_split run_stage_split proc~ocean_dyn_step_split->proc~run_stage_split proc~run_stage_split->proc~run_continuity_chain proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step proc~engine_step->proc~ocean_dyn_step_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 proc~driver_run driver_run proc~driver_run->proc~driver_run_ocean

Variables

Type Visibility Attributes Name Initial
logical, private :: enforce_bed
integer, private :: i
integer, private :: j
integer, private :: k
integer, private :: nx
integer, private :: ny
integer, private :: nz

Source Code

   pure subroutine compute_w_from_continuity(grid, this, ms)
      !! Fill `ms%w_interface` by integrating the horizontal-
      !! continuity residual upward from the bed.  Eulerian z:
      !!
      !!   w_interface(:, :, 1)    = 0  (bed BC)
      !!   w_interface(:, :, k+1)  = w_interface(:, :, k) - flux_h_layer(k)
      !!
      !! With this w, ∂h/∂t = -horizontal_div + (w(k) - w(k+1)) = 0,
      !! so layer thicknesses stay at their initial Eulerian z
      !! positions.  `continuity_compute_fluxes` must run
      !! first — `flux_h_layer` is the input.
      !!
      !! Recurrence is data-dependent in k, so the outer loop is
      !! serial in k inside each column.  Columns are parallelised
      !! over `(j, i)`.
      type(hgrid_t), intent(in) :: grid
      type(ocean_vertical_advection_t), intent(in) :: this
      type(multilayer_state_t), intent(inout) :: ms

      integer :: i, j, k, nx, ny, nz
      logical :: enforce_bed

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      enforce_bed = this%enforce_bed_bc

      do concurrent(j=1:ny, i=1:nx)
         if (enforce_bed) then
            ms%w_interface(i, j, 1) = 0.0_wp
         end if
         do k = 1, nz
            ms%w_interface(i, j, k + 1) = ms%w_interface(i, j, k) - &
                                          ms%flux_h_layer(i, j, k)
         end do
      end do
   end subroutine compute_w_from_continuity