enthalpythickness Subroutine

  subroutine enthalpythickness(output,tinput,uinput)

    use modglobal, only : kb,ke,dzf !,Uinf
    use modsurfdata, only : thls
    implicit none

       real, dimension(kb:ke), intent(in) :: tinput  !< input temperature
       real, dimension(kb:ke), intent(in) :: uinput  !< input velocity
       real, intent(out)                  :: output  !< momentum thickness
       real, dimension(kb:ke)             :: ethick
       real thlsdummy
       integer :: k

       ! An isothermal inlet top (tinput(ke) exactly == thls) makes the enthalpy
       ! -thickness denominator below exactly zero, which -ffpe-trap=zero turns
       ! into a crash. Offsetting thls avoids that. The exact-equality test is
       ! deliberate and must stay: a tolerance test would also capture merely
       ! near-equal values, and clamping those to -0.000001 would flip the sign
       ! of the denominator (and so of ethick) rather than just regularise it.
       ! gfortran's -Wcompare-reals flags this; it is a false positive here.
       thlsdummy = thls
       if (tinput(ke) == thls) then
         thlsdummy = thls -0.000001
       end if
       do k=kb,ke
!         ethick(k) = (uinput(k)/uinput(ke)) * ((tinput(k) - tinput(ke)) /(thls - tinput(ke)) )*dzf(k)
         ethick(k) = (uinput(k)/uinput(ke)) * ((tinput(k) - tinput(ke)) /(thlsdummy - tinput(ke)) )*dzf(k)

       end do
       output   = sum(ethick)  ! enthalpy thickness
       ! Callers divide by the enthalpy thickness, so an exactly-zero sum (an inlet
       ! with no temperature deficit) is regularised rather than allowed to trap.
       ! As above, the exact comparison is intended: only the exactly-zero case is
       ! degenerate, and a small non-zero thickness is physically meaningful.
       if (output==0.) then
         output= 0.000001
       end if

  end subroutine enthalpythickness