redi_neutral_positions_continuous

public pure subroutine redi_neutral_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, Sr, dRdTr, dRdSr, PoL, PoR, KoL, KoR, hEff)

The continuous neutral-surface sweep over a column pair. A single deterministic top→bottom sweep of 2*nk+2 surfaces walking two interface pointers; closed-form linear crossing per step (no inner iteration). Inputs are interface T/S/P + interface dR/dT, dR/dS (nk+1 each). Outputs PoL/PoR (fractional position within layer KoL/KoR) and hEff (harmonic-mean effective thickness between consecutive neutral surfaces; outcrops get hEff=0, not skipped). TOP-DOWN indexing (see module header).

Arguments

Type	Intent		Name
integer,	intent(in)	::	nk
real(kind=wp),	intent(in)	::	Pl(nk+1)	left interface P, T, S
real(kind=wp),	intent(in)	::	Tl(nk+1)	left interface P, T, S
real(kind=wp),	intent(in)	::	Sl(nk+1)	left interface P, T, S
real(kind=wp),	intent(in)	::	dRdTl(nk+1)	left interface dR/dT, dR/dS
real(kind=wp),	intent(in)	::	dRdSl(nk+1)	left interface dR/dT, dR/dS
real(kind=wp),	intent(in)	::	Pr(nk+1)	right interface P, T, S
real(kind=wp),	intent(in)	::	Tr(nk+1)	right interface P, T, S
real(kind=wp),	intent(in)	::	Sr(nk+1)	right interface P, T, S
real(kind=wp),	intent(in)	::	dRdTr(nk+1)	right interface dR/dT, dR/dS
real(kind=wp),	intent(in)	::	dRdSr(nk+1)	right interface dR/dT, dR/dS
real(kind=wp),	intent(out)	::	PoL(2*nk+2)	fractional positions
real(kind=wp),	intent(out)	::	PoR(2*nk+2)	fractional positions
integer,	intent(out)	::	KoL(2*nk+2)	layer indices
integer,	intent(out)	::	KoR(2*nk+2)	layer indices
real(kind=wp),	intent(out)	::	hEff(2*nk+1)	effective thicknesses

Calls

Help

Called by

Help

Variables

Type	Visibility		Name
real(kind=wp),	private	::	dRho
real(kind=wp),	private	::	dRhoBot
real(kind=wp),	private	::	dRhoTop
real(kind=wp),	private	::	hL
real(kind=wp),	private	::	hR
integer,	private	::	k_surface
integer,	private	::	kl
integer,	private	::	klm1
integer,	private	::	kr
integer,	private	::	krm1
integer,	private	::	lastK_left
integer,	private	::	lastK_right
real(kind=wp),	private	::	lastP_left
real(kind=wp),	private	::	lastP_right
integer,	private	::	ns
logical,	private	::	reached_bottom
logical,	private	::	searching_left_column
logical,	private	::	searching_right_column

Source Code

   pure subroutine redi_neutral_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdSl, &
                                                     Pr, Tr, Sr, dRdTr, dRdSr, &
                                                     PoL, PoR, KoL, KoR, hEff)
      !$acc routine seq
      integer, intent(in) :: nk
      real(wp), intent(in) :: Pl(nk + 1), Tl(nk + 1), Sl(nk + 1)     !! left interface P, T, S
      real(wp), intent(in) :: dRdTl(nk + 1), dRdSl(nk + 1)         !! left interface dR/dT, dR/dS
      real(wp), intent(in) :: Pr(nk + 1), Tr(nk + 1), Sr(nk + 1)     !! right interface P, T, S
      real(wp), intent(in) :: dRdTr(nk + 1), dRdSr(nk + 1)         !! right interface dR/dT, dR/dS
      real(wp), intent(out) :: PoL(2*nk + 2), PoR(2*nk + 2)        !! fractional positions
      integer, intent(out) :: KoL(2*nk + 2), KoR(2*nk + 2)         !! layer indices
      real(wp), intent(out) :: hEff(2*nk + 1)                    !! effective thicknesses

      integer :: ns, k_surface, kl, kr, krm1, klm1
      integer :: lastK_left, lastK_right
      real(wp) :: lastP_left, lastP_right
      real(wp) :: dRho, dRhoTop, dRhoBot, hL, hR
      logical :: searching_left_column, searching_right_column, reached_bottom

      ns = 2*nk + 2

      kr = 1
      kl = 1
      lastP_right = 0.0_wp
      lastP_left = 0.0_wp
      lastK_right = 1
      lastK_left = 1
      reached_bottom = .false.
      searching_left_column = .false.
      searching_right_column = .false.

      do k_surface = 1, ns
         klm1 = max(kl - 1, 1)
         krm1 = max(kr - 1, 1)

         ! Cross-column dR/dT, dR/dS-weighted secant density difference:
         ! rho(kr) - rho(kl) (NOT a raw density difference).
         dRho = 0.5_wp*((dRdTr(kr) + dRdTl(kl))*(Tr(kr) - Tl(kl)) &
                        + (dRdSr(kr) + dRdSl(kl))*(Sr(kr) - Sl(kl)))

         if (.not. reached_bottom) then
            if (dRho < 0.0_wp) then
               searching_left_column = .true.
               searching_right_column = .false.
            else if (dRho > 0.0_wp) then
               searching_right_column = .true.
               searching_left_column = .false.
            else  ! dRho == 0: tie-break
               if (kl + kr == 2) then  ! still at surface
                  searching_left_column = .true.
                  searching_right_column = .false.
               else  ! not the surface — change direction
                  searching_left_column = .not. searching_left_column
                  searching_right_column = .not. searching_right_column
               end if
            end if
         end if

         if (searching_left_column) then
            ! rho(kl-1) - rho(kr) (should be negative)
            dRhoTop = 0.5_wp*((dRdTl(klm1) + dRdTr(kr))*(Tl(klm1) - Tr(kr)) &
                              + (dRdSl(klm1) + dRdSr(kr))*(Sl(klm1) - Sr(kr)))
            ! rho(kl) - rho(kr) (will be positive)
            dRhoBot = 0.5_wp*((dRdTl(klm1 + 1) + dRdTr(kr))*(Tl(klm1 + 1) - Tr(kr)) &
                              + (dRdSl(klm1 + 1) + dRdSr(kr))*(Sl(klm1 + 1) - Sr(kr)))

            if (dRhoTop > 0.0_wp .or. kr + kl == 2) then
               PoL(k_surface) = 0.0_wp
            else if (dRhoTop >= dRhoBot) then  ! unstratified left layer
               PoL(k_surface) = 1.0_wp
            else
               PoL(k_surface) = redi_interpolate_position(dRhoTop, Pl(klm1), dRhoBot, Pl(klm1 + 1))
            end if
            if (PoL(k_surface) >= 1.0_wp .and. klm1 < nk) then  ! carry to next layer
               klm1 = klm1 + 1
               PoL(k_surface) = PoL(k_surface) - 1.0_wp
            end if
            ! Monotonic-position backstop
            if (real(klm1 - lastK_left, wp) + (PoL(k_surface) - lastP_left) < 0.0_wp) then
               PoL(k_surface) = lastP_left
               klm1 = lastK_left
            end if
            KoL(k_surface) = klm1
            if (kr <= nk) then
               PoR(k_surface) = 0.0_wp
               KoR(k_surface) = kr
            else
               PoR(k_surface) = 1.0_wp
               KoR(k_surface) = nk
            end if
            if (kr <= nk) then
               kr = kr + 1
            else  ! column exhausted: direction flip
               reached_bottom = .true.
               searching_right_column = .true.
               searching_left_column = .false.
            end if
         else if (searching_right_column) then
            ! rho(kr-1) - rho(kl) (should be negative)
            dRhoTop = 0.5_wp*((dRdTr(krm1) + dRdTl(kl))*(Tr(krm1) - Tl(kl)) &
                              + (dRdSr(krm1) + dRdSl(kl))*(Sr(krm1) - Sl(kl)))
            ! rho(kr) - rho(kl) (will be positive)
            dRhoBot = 0.5_wp*((dRdTr(krm1 + 1) + dRdTl(kl))*(Tr(krm1 + 1) - Tl(kl)) &
                              + (dRdSr(krm1 + 1) + dRdSl(kl))*(Sr(krm1 + 1) - Sl(kl)))

            if (dRhoTop >= 0.0_wp .or. kr + kl == 2) then
               PoR(k_surface) = 0.0_wp
            else if (dRhoTop >= dRhoBot) then  ! unstratified right layer
               PoR(k_surface) = 1.0_wp
            else
               PoR(k_surface) = redi_interpolate_position(dRhoTop, Pr(krm1), dRhoBot, Pr(krm1 + 1))
            end if
            if (PoR(k_surface) >= 1.0_wp .and. krm1 < nk) then  ! carry to next layer
               krm1 = krm1 + 1
               PoR(k_surface) = PoR(k_surface) - 1.0_wp
            end if
            ! Monotonic-position backstop
            if (real(krm1 - lastK_right, wp) + (PoR(k_surface) - lastP_right) < 0.0_wp) then
               PoR(k_surface) = lastP_right
               krm1 = lastK_right
            end if
            KoR(k_surface) = krm1
            if (kl <= nk) then
               PoL(k_surface) = 0.0_wp
               KoL(k_surface) = kl
            else
               PoL(k_surface) = 1.0_wp
               KoL(k_surface) = nk
            end if
            if (kl <= nk) then
               kl = kl + 1
            else  ! column exhausted: direction flip
               reached_bottom = .true.
               searching_right_column = .false.
               searching_left_column = .true.
            end if
         end if

         lastK_left = KoL(k_surface)
         lastP_left = PoL(k_surface)
         lastK_right = KoR(k_surface)
         lastP_right = PoR(k_surface)

         ! Effective thickness between consecutive neutral surfaces (harmonic
         ! mean of the left/right pressure thicknesses; outcrops => hEff=0).
         if (k_surface > 1) then
            hL = redi_absolute_position(nk, Pl, KoL(k_surface), PoL(k_surface)) &
                 - redi_absolute_position(nk, Pl, KoL(k_surface - 1), PoL(k_surface - 1))
            hR = redi_absolute_position(nk, Pr, KoR(k_surface), PoR(k_surface)) &
                 - redi_absolute_position(nk, Pr, KoR(k_surface - 1), PoR(k_surface - 1))
            if (hL + hR > 0.0_wp) then
               hEff(k_surface - 1) = 2.0_wp*hL*hR/(hL + hR)
            else
               hEff(k_surface - 1) = 0.0_wp
            end if
         end if
      end do
   end subroutine redi_neutral_positions_continuous

redi_neutral_positions_continuous Subroutine

Contents

Variables

Source Code

public pure subroutine redi_neutral_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, Sr, dRdTr, dRdSr, PoL, PoR, KoL, KoR, hEff)

Arguments

Calls

Called by

Variables

Source Code