find_duhbt_du Function

   pure function find_duhbt_du(u, BTC) result(duhbt_du)
      !! Marginal zonal face area `d(uhbt)/du`. At `u = 0` returns the average
      !! of the two cubic-branch slopes (discontinuity harmless — only consumed
      !! via `max(…, h_neglect)`).
      !$acc routine seq
      real(wp), intent(in) :: u
      type(local_BT_cont_u_type), intent(in) :: BTC
      real(wp) :: duhbt_du

      if (u == 0.0_wp) then
         duhbt_du = 0.5_wp*(BTC%FA_u_E0 + BTC%FA_u_W0)
      else if (u < BTC%uBT_EE) then
         duhbt_du = BTC%FA_u_EE
      else if (u < 0.0_wp) then
         duhbt_du = BTC%FA_u_E0 + 3.0_wp*BTC%uh_crvE*u**2
      else if (u <= BTC%uBT_WW) then
         duhbt_du = BTC%FA_u_W0 + 3.0_wp*BTC%uh_crvW*u**2
      else
         duhbt_du = BTC%FA_u_WW
      end if
   end function find_duhbt_du