compute_pbce Subroutine

public pure subroutine compute_pbce(grid, bt_work, pgf, ms)

Per-layer pressure-anomaly gravity coefficient (m/s²): the response of layer k’s pressure to a unit change in η. Montgomery form, bottom-up convention (k=1 bed, k=nz surface): pbce(:,:,nz) = g·ρ_ref/ρ_0 do k = nz-1, 1, -1 g_prime_K = g·(rho_layer(k+1) − rho_layer(k))/ρ_0 pbce(:,:,k) = pbce(:,:,k+1) + g_prime_K·(e_top_of_k − e_bed)/H Uniform-density column ⇒ pbce−gtot ≡ 0 ⇒ bc-PGF correction a no-op. Reads pgf%e_face; requires pgf%variant == OPGF_VARIANT_FV_MOM6 (other variants don’t fill e_face). validate_config and configure_ocean_pgf refuse correction_bc_pgf with any other form, so the error stop below is a backstop for direct callers.

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(barotropic_workstate_t), intent(inout) :: bt_work
type(ocean_pressure_force_t), intent(in) :: pgf
type(multilayer_state_t), intent(in) :: ms

Calls

proc~~compute_pbce~~CallsGraph proc~compute_pbce compute_pbce local local proc~compute_pbce->local

Called by

proc~~compute_pbce~~CalledByGraph proc~compute_pbce compute_pbce proc~run_stage_split run_stage_split proc~run_stage_split->proc~compute_pbce proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~run_stage_split proc~engine_step engine_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
real(kind=wp), private :: e_above
real(kind=wp), private :: e_bed
real(kind=wp), private :: g_prime_K
real(kind=wp), private :: g_surf
real(kind=wp), private :: h_col
integer, private :: i
real(kind=wp), private :: inv_rho0
integer, private :: j
integer, private :: k
integer, private :: nx
integer, private :: ny
integer, private :: nz

Source Code

   pure subroutine compute_pbce(grid, bt_work, pgf, ms)
      !! Per-layer pressure-anomaly gravity coefficient (m/s²): the response of
      !! layer k's pressure to a unit change in η. Montgomery form, bottom-up
      !! convention (k=1 bed, k=nz surface):
      !!     pbce(:,:,nz) = g·ρ_ref/ρ_0
      !!     do k = nz-1, 1, -1
      !!        g_prime_K = g·(rho_layer(k+1) − rho_layer(k))/ρ_0
      !!        pbce(:,:,k) = pbce(:,:,k+1) + g_prime_K·(e_top_of_k − e_bed)/H
      !! Uniform-density column ⇒ pbce−gtot ≡ 0 ⇒ bc-PGF correction a no-op.
      !! Reads `pgf%e_face`; requires `pgf%variant == OPGF_VARIANT_FV_MOM6`
      !! (other variants don't fill e_face).  `validate_config` and
      !! `configure_ocean_pgf` refuse `correction_bc_pgf` with any other
      !! form, so the `error stop` below is a backstop for direct callers.
      type(hgrid_t), intent(in) :: grid
      type(barotropic_workstate_t), intent(inout) :: bt_work
      type(ocean_pressure_force_t), intent(in) :: pgf
      type(multilayer_state_t), intent(in) :: ms

      integer :: i, j, k, nx, ny, nz
      real(wp) :: g_surf, inv_rho0, h_col, g_prime_K, e_above, e_bed

      if (pgf%variant /= OPGF_VARIANT_FV_MOM6) then
         error stop "compute_pbce: requires ocean_pgf_form = 'fv_mom6' "// &
            "(pgf%e_face is not populated by other variants)."
      end if

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      inv_rho0 = 1.0_wp/pgf%rho0
      g_surf = GRAVITY*pgf%rho_ref*inv_rho0

      do concurrent(j=1:ny, i=1:nx) local(k, h_col, g_prime_K, e_above, e_bed)
         h_col = 0.0_wp
         do k = 1, nz
            h_col = h_col + ms%h_layer(i, j, k)
         end do
         if (h_col <= 0.0_wp) h_col = 1.0_wp     ! dry cell — pbce won't be used
         e_bed = pgf%e_face%data(i, j, 1)
         bt_work%pbce(i, j, nz) = g_surf
         do k = nz - 1, 1, -1
            g_prime_K = GRAVITY*(ms%rho_layer(i, j, k + 1) - ms%rho_layer(i, j, k))*inv_rho0
            e_above = pgf%e_face%data(i, j, k + 1)
            bt_work%pbce(i, j, k) = bt_work%pbce(i, j, k + 1) &
                                    + g_prime_K*(e_above - e_bed)/h_col
         end do
      end do
   end subroutine compute_pbce