invert_density_targets Subroutine

public pure subroutine invert_density_targets(nk, hc, rhoc, n_int, rho_tgt, z_new)

Density-space interface inversion — the single source of truth for the RHO vcoord regrid (ocean_vcoord_compute_target_h_rho_impl) AND the DENSITY diagnostic remap (rdb_ocean_diag_fills).

Given a TOP-DOWN compacted column of nk layers with thicknesses hc(1:nk) (index 1 = surface, depth positive-down) and layer-mean potential densities rhoc(1:nk), plus n_int monotone-increasing interior target densities rho_tgt(1:n_int), return the n_int + 2 interface depths z_new(1:n_int+2) (top-down, 0 .. column total), monotone non-decreasing. z_new(1) = 0 (surface), z_new(n_int+2) = Σ hc (bed); interior interfaces 2..n_int+1 invert each target.

Algorithm (clean-room from Bleck 2002 / White & Adcroft 2008; MOM6 coord_rho is the behavioural oracle): 1. PPM (Colella-Woodward monotone) reconstruction of the density profile over hc. 2. Per interior target: light boundary → surface; discontinuous- jump sweep; dense boundary → bed; else fixed-NR_ITERS-iter Newton on xi in [0,1] (convergence on |delta| < NR_TOL AFTER xi += delta; zero-gradient NR_OFFSET escape at both ends; masked fallback to the previous interface on no-bracket). 3. Monotone non-decreasing clamp on the interfaces.

Caller supplies nk >= 2. The RHO regrid kernel pre-compacts vanished layers (so nk is the surviving count) and fast-paths nk <= 1 upstream; the DENSITY diagnostic remap passes the full nz column, with every vanished layer given zero thickness and the density of its nearest live neighbour, so the PPM edges it touches are the live layer’s own value (no compaction, same effect on the inversion). Lightest target maps to the surface (index 2), densest to the bed (the surface→bed ordering the callers FLIP into the bottom-up state).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nk
real(kind=wp), intent(in) :: hc(nk)
real(kind=wp), intent(in) :: rhoc(nk)
integer, intent(in) :: n_int
real(kind=wp), intent(in) :: rho_tgt(n_int)
real(kind=wp), intent(out) :: z_new(n_int+2)

Called by

proc~~invert_density_targets~~CalledByGraph proc~invert_density_targets invert_density_targets proc~ocean_vcoord_rho_target_column ocean_vcoord_rho_target_column proc~ocean_vcoord_rho_target_column->proc~invert_density_targets proc~remap_layer_to_density_impl remap_layer_to_density_impl proc~remap_layer_to_density_impl->proc~invert_density_targets proc~ocean_vcoord_rho_target ocean_vcoord_rho_target proc~ocean_vcoord_rho_target->proc~ocean_vcoord_rho_target_column proc~remap_layer_to_density remap_layer_to_density proc~remap_layer_to_density->proc~remap_layer_to_density_impl proc~ocean_vcoord_compute_target_h_rho_impl ocean_vcoord_compute_target_h_rho_impl proc~ocean_vcoord_compute_target_h_rho_impl->proc~ocean_vcoord_rho_target proc~ocean_vcoord_compute_target_h_rho ocean_vcoord_t%ocean_vcoord_compute_target_h_rho proc~ocean_vcoord_compute_target_h_rho->proc~ocean_vcoord_compute_target_h_rho_impl proc~ocean_apply_ale_remap_centres ocean_apply_ale_remap_centres proc~ocean_apply_ale_remap_centres->proc~ocean_vcoord_compute_target_h_rho proc~ocean_apply_ale_remap_step ocean_apply_ale_remap_step proc~ocean_apply_ale_remap_step->proc~ocean_vcoord_compute_target_h_rho

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: dd
real(kind=wp), private :: delta
real(kind=wp), private :: df
real(kind=wp), private :: edge(NZ_STACK_MAX+1)
real(kind=wp), private :: fval
real(kind=wp), private :: grad
real(kind=wp), private :: hi
integer, private :: ii
integer, private :: k
integer, private :: kk
real(kind=wp), private :: lo
logical, private :: placed
real(kind=wp), private :: q6
real(kind=wp), private :: rhoL(NZ_STACK_MAX)
real(kind=wp), private :: rhoR(NZ_STACK_MAX)
real(kind=wp), private :: rho_dense
real(kind=wp), private :: rho_light
real(kind=wp), private :: six
real(kind=wp), private :: tgt
real(kind=wp), private :: ww
real(kind=wp), private :: xi
real(kind=wp), private :: z_old(NZ_STACK_MAX+1)

Source Code

   pure subroutine invert_density_targets(nk, hc, rhoc, n_int, rho_tgt, z_new)
      !! Density-space interface inversion — the single source of truth
      !! for the RHO vcoord regrid (`ocean_vcoord_compute_target_h_rho_impl`)
      !! AND the DENSITY diagnostic remap (`rdb_ocean_diag_fills`).
      !!
      !! Given a TOP-DOWN compacted column of `nk` layers with thicknesses
      !! `hc(1:nk)` (index 1 = surface, depth positive-down) and layer-mean
      !! potential densities `rhoc(1:nk)`, plus `n_int` monotone-increasing
      !! interior target densities `rho_tgt(1:n_int)`, return the `n_int + 2`
      !! interface depths `z_new(1:n_int+2)` (top-down, 0 .. column total),
      !! monotone non-decreasing.  `z_new(1) = 0` (surface), `z_new(n_int+2)`
      !! = Σ hc (bed); interior interfaces 2..n_int+1 invert each target.
      !!
      !! Algorithm (clean-room from Bleck 2002 / White & Adcroft 2008;
      !! MOM6 `coord_rho` is the behavioural oracle):
      !!   1. PPM (Colella-Woodward monotone) reconstruction of the density
      !!      profile over `hc`.
      !!   2. Per interior target: light boundary → surface; discontinuous-
      !!      jump sweep; dense boundary → bed; else fixed-`NR_ITERS`-iter
      !!      Newton on `xi in [0,1]` (convergence on `|delta| < NR_TOL`
      !!      AFTER `xi += delta`; zero-gradient `NR_OFFSET` escape at both
      !!      ends; masked fallback to the previous interface on no-bracket).
      !!   3. Monotone non-decreasing clamp on the interfaces.
      !!
      !! Caller supplies `nk >= 2`.  The RHO regrid kernel pre-compacts
      !! vanished layers (so `nk` is the surviving count) and fast-paths
      !! `nk <= 1` upstream; the DENSITY diagnostic remap passes the full
      !! `nz` column, with every vanished layer given zero thickness and the
      !! density of its nearest live neighbour, so the PPM edges it touches
      !! are the live layer's own value (no compaction, same effect on the
      !! inversion).  Lightest
      !! target maps to the surface (index 2), densest to the bed (the
      !! surface→bed ordering the callers FLIP into the bottom-up state).
      !$acc routine seq
      integer, intent(in)  :: nk
      integer, intent(in)  :: n_int
      real(wp), intent(in)  :: hc(nk)
      real(wp), intent(in)  :: rhoc(nk)
      real(wp), intent(in)  :: rho_tgt(n_int)
      real(wp), intent(out) :: z_new(n_int + 2)
      integer  :: kk, ii, k
      real(wp) :: rhoL(NZ_STACK_MAX), rhoR(NZ_STACK_MAX)
      real(wp) :: edge(NZ_STACK_MAX + 1), z_old(NZ_STACK_MAX + 1)
      real(wp) :: tgt, lo, hi, q6, xi, fval, df, delta, grad, ww
      real(wp) :: rho_light, rho_dense, dd, six
      logical  :: placed

      ! --- step 1: PPM (Colella-Woodward) edges on the compacted column ---
      edge(1) = rhoc(1)
      edge(nk + 1) = rhoc(nk)
      do kk = 2, nk
         ww = hc(kk - 1) + hc(kk)
         ! Floor guards an uncompacted caller (the DENSITY diagnostic remap
         ! passes the raw column) where a vanished layer pair sums to ~0;
         ! the RHO regrid pre-compacts so ww is always >> the floor there
         ! (this branch is inert for it — bit-identical).
         if (ww < 1.0e-30_wp) ww = 1.0e-30_wp
         edge(kk) = (hc(kk)*rhoc(kk - 1) + hc(kk - 1)*rhoc(kk))/ww
      end do
      do kk = 1, nk
         rhoL(kk) = edge(kk)
         rhoR(kk) = edge(kk + 1)
         ! CW monotonic limiter
         if ((rhoR(kk) - rhoc(kk))*(rhoc(kk) - rhoL(kk)) <= 0.0_wp) then
            rhoL(kk) = rhoc(kk)
            rhoR(kk) = rhoc(kk)
         else
            dd = rhoR(kk) - rhoL(kk)
            six = 6.0_wp*(rhoc(kk) - 0.5_wp*(rhoL(kk) + rhoR(kk)))
            if (dd*six > dd*dd) then
               rhoL(kk) = 3.0_wp*rhoc(kk) - 2.0_wp*rhoR(kk)
            else if (dd*six < -dd*dd) then
               rhoR(kk) = 3.0_wp*rhoc(kk) - 2.0_wp*rhoL(kk)
            end if
         end if
      end do

      ! Cumulative compacted-grid interface depths (top-down, 0..H).
      z_old(1) = 0.0_wp
      do kk = 1, nk
         z_old(kk + 1) = z_old(kk) + hc(kk)
      end do
      rho_light = rhoL(1)
      rho_dense = rhoR(nk)

      ! --- step 2: invert each interior target interface (top-down) ---
      z_new(1) = 0.0_wp
      z_new(n_int + 2) = z_old(nk + 1)     ! total compacted depth (= col extent)
      do kk = 2, n_int + 1                  ! interior interfaces
         tgt = rho_tgt(kk - 1)
         if (tgt <= rho_light) then
            z_new(kk) = 0.0_wp              ! lighter than column -> surface
         else if (tgt >= rho_dense) then
            z_new(kk) = z_old(nk + 1)       ! denser than column -> bed
         else
            placed = .false.
            do ii = 1, nk
               lo = rhoL(ii)
               hi = rhoR(ii)
               ! Discontinuous jump at the TOP interface of cell ii
               ! (between cell ii-1's right edge and cell ii's left
               ! edge).  Checked FIRST and independent of whether the
               ! cells are limiter-flattened — a 2-layer column flattens
               ! both cells, and a target inside the jump must still
               ! land on the interface (MOM6 coord_rho behaviour).
               if (ii > 1) then
                  if (rhoR(ii - 1) <= tgt .and. tgt <= rhoL(ii)) then
                     z_new(kk) = z_old(ii)
                     placed = .true.
                     exit
                  end if
               end if
               if (lo == hi) then
                  ! Flat (limiter-collapsed) cell: only an exact match
                  ! places here; otherwise advance to the next cell.
                  if (abs(tgt - lo) < NR_TOL) then
                     z_new(kk) = z_old(ii)
                     placed = .true.
                     exit
                  end if
                  cycle
               end if
               if ((lo - tgt)*(hi - tgt) <= 0.0_wp) then
                  ! Newton on xi in [0,1]: ppm(xi) - tgt = 0
                  q6 = 6.0_wp*rhoc(ii) - 3.0_wp*(lo + hi)
                  xi = 0.5_wp
                  do k = 1, NR_ITERS
                     fval = lo + xi*((hi - lo) + q6*(1.0_wp - xi)) - tgt
                     df = (hi - lo) + q6*(1.0_wp - 2.0_wp*xi)
                     if (abs(df) > 1.0e-30_wp) then
                        delta = -fval/df
                     else
                        delta = 0.0_wp
                     end if
                     xi = xi + delta
                     ! clamp inside the iteration; zero-gradient nudge
                     if (xi < 0.0_wp) then
                        xi = 0.0_wp
                        grad = (hi - lo) + q6           ! d(ppm)/dxi at xi=0
                        if (abs(grad) < 1.0e-30_wp) xi = NR_OFFSET
                     else if (xi > 1.0_wp) then
                        xi = 1.0_wp
                        grad = (hi - lo) - q6           ! d(ppm)/dxi at xi=1
                        if (abs(grad) < 1.0e-30_wp) xi = 1.0_wp - NR_OFFSET
                     end if
                     if (abs(delta) < NR_TOL) exit
                  end do
                  z_new(kk) = z_old(ii) + xi*hc(ii)
                  placed = .true.
                  exit
               end if
            end do
            if (.not. placed) then
               ! masked fallback: previous interface (no FATAL in a DC)
               z_new(kk) = z_new(kk - 1)
            end if
         end if
      end do

      ! --- step 3: monotone non-decreasing interfaces ---
      do kk = 2, n_int + 2
         if (z_new(kk) < z_new(kk - 1)) z_new(kk) = z_new(kk - 1)
      end do
   end subroutine invert_density_targets