cavity_gamma_hj99 Subroutine

   pure subroutine cavity_gamma_hj99(u_star, l_plus, f_cor, const, gamma_t, gamma_s, ierr)
      !! Holland & Jenkins (1999) eqs. (14)-(18) p. 1792:
      !!
      !!   `gamma_{T,S} = u* / (Gamma_Turb + Gamma_Mole^{T,S})`        (14)
      !!   `Gamma_Turb  = (1/k)*ln( u* xi_N eta*^2/(|f| h_nu) )
      !!                  + 1/(2 xi_N eta*) - 1/k`                     (15)
      !!   `Gamma_Mole  = 12.5*(Pr,Sc)^(2/3) - 6`                      (16)
      !!   `h_nu        = 5 nu/u*`                                     (17)
      !!   `eta*        = ( 1 + xi_N u*/(|f| L_O R_c) )^(-1/2) <= 1`   (18)
      !!
      !! Identical to McPhee, Maykut & Morison (1987) eq. (10) p. 7029
      !! with the roughness length replaced by the viscous-sublayer
      !! thickness (the hydraulically-smooth assumption, H&J99 p. 1792).
      !!
      !! Check value: `u* = 1e-2 m/s`, `eta* = 1`, `|f| = 1e-4` gives
      !! `gamma_T = 1.06e-4 m/s` and `gamma_S/gamma_T = 0.041`, against
      !! H&J99 Table 1's `gamma_T ~ 1.0e-4` and their p. 1797 ratio of
      !! 0.04.  (Table 1's companion `gamma_S ~ 5.05e-7` is a Hellmer &
      !! Olbers constant-coefficient value, NOT this formulation — the two
      !! are not a self-consistent check pair.)
      !!
      !! BRANCH (must be reproduced exactly for oracle agreement): a
      !! destabilising or vanishing buoyancy flux sets `eta* := 1`.
      !!
      !! GUARDS: `|f| = 0` is refused outright (`CAVITY_MELT_NO_CORIOLIS`)
      !! — the law has `|f|` inside a logarithm and divides by it, so it
      !! does not exist on the equator.  A non-positive logarithm
      !! argument, a non-positive `eta*` argument and a non-positive
      !! denominator are all refused rather than clamped.
      !$acc routine seq
      real(wp), intent(in) :: u_star
         !! Friction velocity (m/s), > 0.
      real(wp), intent(in) :: l_plus
         !! Trial viscous Obukhov scale, or `CAVITY_L_PLUS_NEUTRAL`.
      real(wp), intent(in) :: f_cor
         !! Coriolis parameter (1/s), used as `|f|`.  A SCALAR, not the
         !! exchange bundle: see `cavity_exchange_velocities_f`.
      type(ocean_cavity_const_t), intent(in) :: const
         !! Constants bundle.
      real(wp), intent(out) :: gamma_t
         !! Heat exchange velocity (m/s).
      real(wp), intent(out) :: gamma_s
         !! Salt exchange velocity (m/s).
      integer, intent(out) :: ierr
         !! `CAVITY_MELT_*` status.
      real(wp) :: fa, eta, arg, l_obukhov, h_nu, g_turb, den_t, den_s

      ierr = CAVITY_MELT_OK
      gamma_t = 0.0_wp
      gamma_s = 0.0_wp

      if (.not. ieee_is_finite(f_cor)) then
         ierr = CAVITY_MELT_NONFINITE_INPUT
         return
      end if
      fa = abs(f_cor)
      if (fa <= 0.0_wp) then
         ierr = CAVITY_MELT_NO_CORIOLIS
         return
      end if

      ! eq. (18), with the H&J99 p. 1792 destabilising branch.  `L+` is
      ! converted back to the DIMENSIONAL Obukhov length the paper uses,
      ! `L_O = L+ * delta_nu = L+ * nu/u*`.
      if (cavity_l_plus_is_neutral(l_plus)) then
         eta = 1.0_wp
      else
         l_obukhov = l_plus*const%nu/u_star
         if ((.not. ieee_is_finite(l_obukhov)) .or. l_obukhov <= 0.0_wp) then
            eta = 1.0_wp
         else
            arg = 1.0_wp + const%xi_N*u_star/(fa*l_obukhov*const%R_c)
            if (.not. ieee_is_finite(arg)) then
               ierr = CAVITY_MELT_NONFINITE_STATE
               return
            end if
            if (arg <= 0.0_wp) then
               ierr = CAVITY_MELT_LAW_DOMAIN
               return
            end if
            eta = arg**(-0.5_wp)
         end if
      end if

      ! eqs. (15)+(17)
      h_nu = HJ99_H_NU_COEFF*const%nu/u_star
      arg = u_star*const%xi_N*eta*eta/(fa*h_nu)
      if (.not. ieee_is_finite(arg)) then
         ierr = CAVITY_MELT_NONFINITE_STATE
         return
      end if
      if (.not. (arg > 0.0_wp)) then
         ierr = CAVITY_MELT_LAW_DOMAIN
         return
      end if
      g_turb = log(arg)/const%kappa_vk + 1.0_wp/(2.0_wp*const%xi_N*eta) &
               - 1.0_wp/const%kappa_vk

      ! eq. (16) + eq. (14)
      den_t = g_turb + (HJ99_MOLE_SLOPE*const%Pr**(2.0_wp/3.0_wp) - HJ99_MOLE_OFFSET)
      den_s = g_turb + (HJ99_MOLE_SLOPE*const%Sc**(2.0_wp/3.0_wp) - HJ99_MOLE_OFFSET)
      if (.not. (ieee_is_finite(den_t) .and. ieee_is_finite(den_s))) then
         ierr = CAVITY_MELT_NONFINITE_STATE
         return
      end if
      if (den_t <= 0.0_wp .or. den_s <= 0.0_wp) then
         ! `Gamma_Turb` can go strongly negative at very small `u*`.  The
         ! molecular terms (65.9 heat, 2254.6 salt) normally dominate;
         ! refuse rather than return a negative exchange velocity.
         ierr = CAVITY_MELT_LAW_DOMAIN
         return
      end if
      gamma_t = u_star/den_t
      gamma_s = u_star/den_s
   end subroutine cavity_gamma_hj99