subroutine heat_transfer_coef_flux(utan, dist, z0, z0h, Tair, Tsurf, cth, flux, htc)
use modglobal, only : grav, fkar, prandtlturb
implicit none
real, intent(in) :: dist, z0, z0h, Tsurf, Tair, utan
real, intent(out) :: cth, flux, htc
real, parameter :: b1 = 9.4 !parameters from Uno1995
real, parameter :: b2 = 4.7
real, parameter :: dm = 7.4
real, parameter :: dh = 5.3
!real :: Pr
real :: dT, Ribl0, logdz, logdzh, logzh, sqdz, fkar2, Ribl1, Fm, Fh, cm, ch, M, dTrough
!Pr = 1.
!Pr = prandtlmol
dT = Tair - Tsurf
Ribl0 = grav * dist * dT / (Tsurf * utan**2) !Eq. 6, guess initial Ri
logdz = log(dist/z0)
logdzh = log(dist/z0h)
logzh = log(z0/z0h)
sqdz = sqrt(dist/z0)
fkar2 = fkar**2
cth = 0.
flux = 0.
if (Ribl0 > 0.) then
Fm = 1./(1. + b2*Ribl0)**2 !Eq. 4
Fh = Fm !Eq. 4
else ! => unstable
cm = (dm*fkar2)/(logdz**2)*b1*sqdz !Eq. 5
ch = (dh*fkar2)/(logdz**2)*b1*sqdz !Eq. 5
Fm = 1. - (b1*Ribl0)/(1. + cm*sqrt(abs(Ribl0))) !Eq. 3
Fh = 1. - (b1*Ribl0)/(1. + ch*sqrt(abs(Ribl0))) !Eq. 3
end if
M = prandtlturb*logdz*sqrt(Fm)/Fh !Eq. 14
Ribl1 = Ribl0 - Ribl0*prandtlturb*logzh/(prandtlturb*logzh + M) !Eq. 17
!interate to get new Richardson number
if (Ribl1 > 0.) then
Fm = 1./(1. + b2*Ribl1)**2 !Eq. 4
Fh = Fm !Eq. 4
else ! => unstable
cm = (dm*fkar2)/(logdz**2)*b1*sqdz !Eq. 5
ch = (dh*fkar2)/(logdz**2)*b1*sqdz !Eq. 5
Fm = 1. - (b1*Ribl1)/(1. + cm*sqrt(abs(Ribl1))) !Eq. 3
Fh = 1. - (b1*Ribl1)/(1. + ch*sqrt(abs(Ribl1))) !Eq. 3
end if
! Uno (2)
M = prandtlturb*logdz*sqrt(Fm)/Fh !Eq. 14
dTrough = dT*1./(prandtlturb*logzh/M + 1.) !Eq. 13a
cth = fkar2/(logdz*logdz)*Fh/prandtlturb ! Ivo's heat transfer coefficient
flux = abs(utan)*cth*dTrough
if (abs(abs(utan)*dT) > 0.) then
htc = flux / (abs(utan)*dT)
else
htc = 0.
end if
! ! Uno (8)
! cth = abs(utan)*fkar2/(logdz*logdzh)*Fh/prandtlturb !Eq. 8
! flux = cth*dT !Eq. 2, Eq. 8
end subroutine heat_transfer_coef_flux