real function mom_transfer_coef_stability(utan, dist, z0, z0h, Tair, Tsurf)
! By Ivo Suter. calculates the momentum transfer coefficient based on the
! surface tangential velocity 'utan' at a distance 'dist' from the surface,
! for a surface with momentum roughness length z0 and heat roughness length z0h.
! Stability are included using the air temperature Tair and surface temperature Tsurf.
use modglobal, only : grav, fkar, prandtlturb
implicit none
real, intent(in) :: dist, z0, z0h, Tsurf, Tair, utan
real, parameter :: b1 = 9.4 !parameters from uno1995
real, parameter :: b2 = 4.7
real, parameter :: dm = 7.4
real, parameter :: dh = 5.3
real :: dT, Ribl0, logdz, logdzh, logzh, sqdz, fkar2, Ribl1, Fm, Fh, cm, ch, Ctm, M
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
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
ELSE ! => unstable
cm = (dm*fkar2)/(logdz**2)*b1*sqdz !Eq. 5
Fm = 1. - (b1*Ribl1)/(1. + cm*SQRT(ABS(Ribl1))) !Eq. 3
END IF
mom_transfer_coef_stability = fkar2/(logdz**2)*Fm !Eq. 7
end function mom_transfer_coef_stability