uhbt_to_ubt Function

public pure function uhbt_to_ubt(uhbt, BTC) result(ubt)

Invert find_uhbt: recover u from a target transport uhbt. Saturated branches close in one line; cubic branches use Newton + bisection fallback to tol·|uhbt|. Hallberg & Adcroft (2009).

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: uhbt
type(local_BT_cont_u_type), intent(in) :: BTC

Return Value real(kind=wp)


Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: derr_du
integer, private :: itt
integer, private, parameter :: max_itt = 20
real(kind=wp), private, parameter :: tol = 1.0e-10_wp
real(kind=wp), private :: ubt_max
real(kind=wp), private :: ubt_min
real(kind=wp), private :: uhbt_err
real(kind=wp), private :: uherr_max
real(kind=wp), private :: uherr_min

Source Code

   pure function uhbt_to_ubt(uhbt, BTC) result(ubt)
      !! Invert `find_uhbt`: recover `u` from a target transport `uhbt`.
      !! Saturated branches close in one line; cubic branches use Newton +
      !! bisection fallback to tol·|uhbt|. Hallberg & Adcroft (2009).
      real(wp), intent(in) :: uhbt
      type(local_BT_cont_u_type), intent(in) :: BTC
      real(wp) :: ubt

      real(wp) :: ubt_min, ubt_max
      real(wp) :: uhbt_err, derr_du
      real(wp) :: uherr_min, uherr_max
      real(wp), parameter :: tol = 1.0e-10_wp
      integer, parameter :: max_itt = 20
      integer :: itt

      if (uhbt == 0.0_wp) then
         ubt = 0.0_wp
      else if (uhbt < BTC%uh_EE) then
         ubt = BTC%uBT_EE + (uhbt - BTC%uh_EE)/BTC%FA_u_EE
      else if (uhbt < 0.0_wp) then
         ubt_min = BTC%uBT_EE
         uherr_min = BTC%uh_EE - uhbt
         ubt_max = 0.0_wp
         uherr_max = -uhbt
         ubt = BTC%uBT_EE*(uhbt/BTC%uh_EE)
         do itt = 1, max_itt
            uhbt_err = ubt*(BTC%FA_u_E0 + BTC%uh_crvE*ubt**2) - uhbt
            if (abs(uhbt_err) < tol*abs(uhbt)) exit
            if (uhbt_err > 0.0_wp) then
               ubt_max = ubt
               uherr_max = uhbt_err
            end if
            if (uhbt_err < 0.0_wp) then
               ubt_min = ubt
               uherr_min = uhbt_err
            end if
            derr_du = BTC%FA_u_E0 + 3.0_wp*BTC%uh_crvE*ubt**2
            if ((uhbt_err >= derr_du*(ubt - ubt_min)) .or. &
                (-uhbt_err >= derr_du*(ubt_max - ubt)) .or. (derr_du <= 0.0_wp)) then
               ubt = ubt_max + (ubt_min - ubt_max)*(uherr_max/(uherr_max - uherr_min))
            else
               ubt = ubt - uhbt_err/derr_du
               if (abs(uhbt_err) < (0.01_wp*tol)*abs(ubt_min*derr_du)) exit
            end if
         end do
      else if (uhbt <= BTC%uh_WW) then
         ubt_min = 0.0_wp
         uherr_min = -uhbt
         ubt_max = BTC%uBT_WW
         uherr_max = BTC%uh_WW - uhbt
         ubt = BTC%uBT_WW*(uhbt/BTC%uh_WW)
         do itt = 1, max_itt
            uhbt_err = ubt*(BTC%FA_u_W0 + BTC%uh_crvW*ubt**2) - uhbt
            if (abs(uhbt_err) < tol*abs(uhbt)) exit
            if (uhbt_err > 0.0_wp) then
               ubt_max = ubt
               uherr_max = uhbt_err
            end if
            if (uhbt_err < 0.0_wp) then
               ubt_min = ubt
               uherr_min = uhbt_err
            end if
            derr_du = BTC%FA_u_W0 + 3.0_wp*BTC%uh_crvW*ubt**2
            if ((uhbt_err >= derr_du*(ubt - ubt_min)) .or. &
                (-uhbt_err >= derr_du*(ubt_max - ubt)) .or. (derr_du <= 0.0_wp)) then
               ubt = ubt_min + (ubt_max - ubt_min)*(-uherr_min/(uherr_max - uherr_min))
            else
               ubt = ubt - uhbt_err/derr_du
               if (abs(uhbt_err) < (0.01_wp*tol)*(ubt_max*derr_du)) exit
            end if
         end do
      else
         ubt = BTC%uBT_WW + (uhbt - BTC%uh_WW)/BTC%FA_u_WW
      end if
   end function uhbt_to_ubt