SUBROUTINE unoh(otf, octh, logdz, logdzh, logzh, sqdz, utangInt, dT, Ribl0, fkar2) !for heat, the bit that does not change no matter what wall
use modglobal, only : prandtlturb
!flux in Km/s
IMPLICIT NONE
REAL, INTENT(in) :: logdz, logdzh, logzh, sqdz, utangInt, dT, Ribl0, fkar2
REAL, INTENT(out) :: octh, otf
REAL :: Ribl1, Fm, Fh, cm, ch, M, dTrough, cth
REAL, PARAMETER :: b1 = 9.4 !parameters from Uno1995
REAL, PARAMETER :: b2 = 4.7
REAL, PARAMETER :: dm = 7.4
REAL, PARAMETER :: dh = 5.3
octh = 0.
otf = 0.
IF (Ribl0 > 0.) THEN !0.25 approx critical for bulk Richardson number => stable
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 !0.25 approx critical for bulk Richardson number => stable
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
octh = SQRT(utangInt)*fkar2/(logdz*logdz)*Fh/prandtlturb !Eq. 8
otf = octh*dTrough !Eq. 2, Eq. 8
! ! Uno (8)
! octh = SQRT(utangInt)*fkar2/(logdz*logdzh)*Fh/prandtlturb !Eq. 8
! otf = octh*dT !Eq. 2, Eq. 8
END SUBROUTINE unoh