REAL FUNCTION unom(logdz, logdzh, logzh, sqdz, utangInt, dT, Ribl0, fkar2) !for momentum, this bit is not depended on orientation etc
!momentum flux in m2/s2
!dT,utang and logdzh are unused and could be removed
IMPLICIT NONE
REAL, INTENT(in) :: logdz, logdzh, logzh, sqdz, utangInt, dT, Ribl0, fkar2
REAL :: Ribl1, Fm, Fh, cm, ch, Ctm, M
REAL, PARAMETER :: b1 = 9.4 !parameters from Uno1995
REAL, PARAMETER :: b2 = 4.7
REAL, PARAMETER :: dm = 7.4
REAL, PARAMETER :: dh = 5.3
REAL, PARAMETER :: prandtlmol = 0.71
IF (Ribl0 > 0.21) 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 = prandtlmol*logdz*SQRT(Fm)/Fh !Eq. 14
Ribl1 = Ribl0 - Ribl0*prandtlmol*logzh/(prandtlmol*logzh + M) !Eq. 17
!interate to get new Richardson number
IF (Ribl1 > 0.21) THEN !0.25 approx critical for bulk Richardson number => stable
Fm = 1./(1. + b2*Ribl1)**2 !Eq. 4
ELSE ! => unstable
cm = (dm*fkar2)/(logdz**2)*b1*sqdz !Eq. 5
Fm = 1. - (b1*Ribl1)/(1. + cm*SQRT(ABS(Ribl1))) !Eq. 3
END IF
Ctm = fkar2/(logdz**2)*Fm !Eq. 7
unom = Ctm !Eq. 2, Eq. 8
END FUNCTION unom