M, H and bM, H are constants, then the critical gradient Richardson number is given
by
.
As the critical number is approached, FM,H,Q vanish. For example if we take FM,H to be
where M, H and bM, H are constants, then the critical gradient Richardson number is given
by
The critical Richardson number for the profile relations in the MRF is 1.32. This is about 6-7
times higher than the RC for the Dyer (1974) or for the commonly used Businger et al. (1971)
profile relations. It is likely, however, that even RC = 1.32 is too low. When R equals or
exceeds, RC, there is no downward heat flux. The main consequence of a low RC is that too
much of the upward daytime heat flux is trapped in the boundary layer and cannot return
during stable nocturnal conditions. A relaxation of the RC is needed to enhance the downward
flux.
The lack of a precise knowledge of the stable boundary layer's profile laws permits some
constructive tampering, even to the extent of eliminating (moving to +) the critical
Richardson number. This should not be interpreted as a denial of the existence of RC
(probably ~0.20) at a single point, but rather the recognition that an area average of the
nocturnal boundary layer includes portions of small, patchy turbulent eddies mixed with
gravity waves. Flux leakage is realistic and reasonable.
(b) In this section we examine a simple formulation for FM, H that exhibits neither a critical
gradient nor a critical bulk Richardson number. In addition, the formulation does not require
dividing z into separate weakly, mildly, and strongly stable regimes.
The formulation is motivated by the turbulent diffusion equations for R o that the MRF uses
above the surface layer,
To make this relation compatible with the surface layer, we note the following surface layer
relations,
As usual, FH = FQ. For mildly stable cases, experimental evidence suggests FM = FH. For
increasing stability, it is possible that FM and FH begin to differ, but we will not deal with this
possibility at this time. We shall assume that FM = FH for all positive z.
Eqs. (4.4) - (4.6) can also be written in 'K-theory' form:
where
follow by comparison of (4.121) - (4.123) with (4.124) - (4.126). These diffusion coefficients
can be related to (3) provided that in the surface layer we have: the (neutral) mixing length:
1 = kz; the wind shear: 
; the stability function:
since the connecting relation between z and R is
the 'universal' functions FM, H are simply
this means that
or
and
For small values of z (weak stability), we have the standard linear results
but for large , z we have
Eqs. 4.21 - 4.22 mean that for strongly stable conditions the fluxes of momentum and heat are
nonzero and are given by
As 
is increased, the fluxes decrease.
(c) In this section we show how z can be estimated from the bulk Richardson number, RB,
without iteration. We first use the equation that relates z to RB,
in which
From (4.132), we have,
since FM, H is conventionally written in the form
we must change the form of (4.142). The result is
in which
Another standard form for (4.143) is
in which
Now, some manipulation shows that FM,H can be approximated for
= 5 by,
In addition, we can roughly approximate z by
which leads to the quadratic equation
in which
The physical correct solution to (4.152) is
The solution for z can be divided into two portions. Numerical computations show that z can
be adequately approximated by
Substitution of the approximate value of 
into (4.144) leads to momentum and heat
transfer, coefficients CM, H that differ only slightly from the 'exact' values. For 
the approximation 
should not be used and should be replaced instead by 
where is given by (4.153) - (4.154). As with 
should be substituted (4.144) in order to
approximate FM, H which then leads to the transfer coefficients,
where aM, H, 
M, H and pM, H are positive constants. Some examples are: (A) Businger et al.
(1971), {aM, M and pM} = {1, 15, 1/4}; {aH, H and pH} = {0.74, 9, 1/2}; (B) Dyer (1967)
aM=aH=1, M,H=15, pM=0.275, pH=0.55; (C) Carl et al. (1973), {aM, M and pM} = {1, 16,
1/3}; {aH, H and pH} = {.74, 16, 1/4}. The MRF physics formulation uses the results of Dyer
(1974) and Hicks (1976):
(b) In this section we show how (4.157) can be used to compute FM, H = ln (Z/Zo) - YM, H
(z,zo) = ln(Z/Zo) - YM, H + YM, H ( zo). When - zo is small, then (4.160) can be
approximated by
as it commonly done. Now, for 0 
- z 0.5, YM, H (z ) can be approximated by
To determine a0, a1, b1, ( 4.161) is expanded:
Also, we have,
By equating a0=-4, a1-a0 b1 =-20 and requiring YM, (-z=0.5) = 0.7934, we get
for -
0.5. Additional accuracy can be achieved by forcing collocation at -z=0, 0.05, 0.10,
0.25, and 0.50:
Similar calculations for FH give,
and
We now consider the more difficult general problem of determining YM, H from YM, H =
(1-
z ) -p. The task is tedious, but the asymptotic result, Eq. 4.76 is fairly simple. As -z
increases, the accuracy of YM, H increases.
We first express Y as
in which x = - Z/L. We also express (4.168) as Y(x) = I1 (x0)-I2(x)
in which
and
Integral I1 can be broken up into Ja and Jb:
The second term in J a is the exponential integral defined by (Arfken, 1985)
and given by the series expansion
for which 
is the Eulaer-Masheroni constant,
Substitution of (4.174) - (4.175) into (4.172) yields
We are left with a logarithmic divergence in Ja, but this will be shown to cancel with another
term.
We now examine Jb by noting that an integral representation of the logarithm of the gamma
function G(t) can be written as
The integral can be differentiated to get
which, of course, is the same as [1/G(u)] [dG(u)/du], the definition of the Euler 'digamma' or
'psi' function (Gradshteyn and Ryzhik, 1965),
We can change variable in the second term of Jb:
to get
and, therefore,
Values of the digamma function have been tabulated. Some values that are of interest to us
are:
It remains to evaluate I2. We write
in which S (s;p) is the series
For FM (Z/L) with p= 1/4 and FH (Z/L) with p=1/2, we get
and
The presence of ln r in 
cancels the divergence in Ja. The final result is obligingly compact:
When (4.190) is applied to FM and FH, we get
We need to retain only the (-z) -1/4 and (-z) -1/2 terms for -z>0.5 to achieve an adequate
accuracy of 1.8% or greater. As expected, as -z increases, the accuracy increases.
(c) In this section, we delve into the physical and computational difficulties that arise when
h is not many times larger than Z0 and when |L|->Z0. These conditions are not normally
encountered in field conditions whose goal is to measure FM, H (Z/L), but large Z0 and small
|L| must be expected to occur fairly frequently in global modeling and must be dealt with.
We deal with the problem of Zo~O(h) first. It is not expected that the logarithmic wind law
is valid for Z ~Zo. Tennekes (1973) has suggested that the region for validity is Z> 10-100
Zo. Following Garratt (1980), we denote by Z* the lowest Z for which the profile laws are
approximately valid for neutral and unstable lapse rates. The region Z<Z* is Garratt's
'roughness sublayer'.
At heights lower than Z* the individual roughness elements can be 'sensed' as the result of
turbulent wakes created by flow around individual elements. This violates the similarity
conditions by introducing an additional length scale, Z*, into the conventional profiles laws,
FM, H. Garratt's analysis of data from very rough terrain suggests that Z*
35 Zo for
momentum and that Z*100 Zo for heat. This means that if the largest roughness lengths in
the MRF are say, Zo5m, then the standard profile laws do not apply to heights below
175m for momentum and 500m for heat. These depths are as great as or greater than the
heights of the surface layers of most NWP models.
In our calculations, we shall take a more liberal approach and merely require h>10 Zo. For
points where Zo>h/10 for the FM, H calculations.
We now deal with the second question. The Businger-Hicks profile laws are regarded as
reasonably accurate for -z<3. Model computations, however, can significantly exceed this
moderate limitation.
Numerical calculations have uncovered another difficulty. Consider the following situation:
and h are held constant as U(h) increases -RB and -z. This decreases FH and FM,
which, in turn, increases CH. Under normal conditions (that is, |L| >>Zo), the decrease in
U(h) overcomes the increase in CH, and the surface heat flux decreases. For very large Zo and
small |L|, there is an anomaly: As U(h) decreases, the heat flux decreases, reaches a plausible
minimum, and then increases. Eventually, heat flux that equal and exceed by several times the
solar constant can be (mathematically) reached. This is counter-intuitive.
These 'runaway' heat fluxes can be eliminated by noting the following pattern that begins to
emerge after observing the behavior of FH for many examples: Although the examples do not
seem to indicate a minimum -u* Q* for any particular value of U(h), RB, h, or a minimum
occurs when
for m20-50. To help eliminate the problem of runaway heat and humidity fluxes, we choose
This choice has the added advantage of allowing us to neglect Y(Zo/L) compared to Y(Z/L).
This means that YM,H (Z/L; Zo/L) can be accurately approximated by YM, H.
To achieve the restriction imposed by (4.194), U(h) is increased in
That is, we require
for the CM,H,Q and flux calculations. Once these quantities are computed, U(h) is returned
to its original value.
Carl, D. M., T. C. Tarbell, and H. A. Panofsky, 1973: Profiles of wind and temperatures from towers over homogeneous terrain. J. Atmos. Sci., 30, 788-794.
Dyer, A. J., 1967: The turbulent transport of heat and water vapor in an unstable atmosphere. Quart. J. Roy. Meteor. Soc., 87, 401-412.
Garratt, J. R., 1980: Surface influence upon vertical profiles in the atmospheric near-surface layer. Quart. J. Roy. Meteor. Soc., 106, 803-819.
Gradshteyn, I. S. and I. M. Ryzhik, 1965: Table of Integrals, Series and Products, Fourth Edition. Academic Press, Inc., New York. 1086pp.
Hicks, B. B., 1976: Wind profile relationships from the 'Wangara' experiment. Quart. J. Roy. Meteor. Soc., 102, 535-551.
Tennekes, H., 1973: A model for the dynamics of the inversion above a convective boundary layer. J. Atmos. Sci., 30, 558-567.