1. PRECIPITATION AND DIFFUSION PROCESSES
The parameterization of non-adiabatic atmospheric physics in the model (excluding radiation) consists of vertical diffusion, deep and shallow convection, and grid-scale condensation/precipitation. It is performed in two separate sections; vertical turbulent diffusion of momentum, heat, and water vapor is carried out just after the non-linear dynamics, and the tendencies due to all of these processes are added together.

The convection and condensation calculations are performed after each leap-frog forward time step and can be considered as adjustments. Deep convection is done first, followed by shallow convection, and finally grid-scale precipitation. No dry adiabatic adjustment is included

1. Vertical Diffusion
Subroutine MONIN3 calculates the effects of vertical turbulent eddy transfer of momentum, heat, and moisture throughout the atmosphere. It uses diffusion coefficients that are dependent upon stability (via the Richardson number) and upon height (via the mixing length). It incorporates lower boundary fluxes by the use of drag coefficients that have been calculated in subroutine PROGTN. MONIN3 takes as input the conditions at time -l and returns solutions in the form of time tendencies of u, v, T, and q, to which will be added the non-linear tendencies due to adiabatic dynamics that have already been calculated in subroutine FIDI. A split implicit numerical approximation is used to set up discretized versions of the basic equations, which are reduced to tri-diagonal coefficient matrices under the boundary conditions

1. Basic Equations

KM (m2/sec) is the diffusion coefficient for momentum and KQT the coefficient for heat and moisture; is the dry adiabatic lapse rate (lK/l00 m), and A, B, C, and D are the non-linear tendencies.

1. Calculations of the Diffusion Coefficients
In general, we use

where is the mixing length, taken as =kz/(1+kz/ ), with k=0.4, =250 m, and S a set of semi-empirical stability functions depending on the Richardson number, R:

For statically stable conditions, R > 0, and K is the same for heat momentum, and moisture with

For unstable conditions, R < 0, we use the expressions

for momentum, and

for heat and moisture.

When S is expanded in a Taylor series to lowest order R, it yields the stable surface layer relations used in Section 4. For large R, S diminishes slowly with increasing R, but never reaches zero (which would happen for critical R). It has been shown by P. Long (unpublished) that the unstable relations for S represent the lowest order asymptotically correct Pade approximations (see Blackadar, 1979). For strong instability approaching free convection, K becomes independent of wind shear

The wind shear is constrained to be at least l m/sec across each layer. The values of K are limited to the range 0.l < K < 300 m2/sec in such a way that if either KM or KQT at a given point is outside this range, both are set back to the appropriate limiting value.

1. Method of Solution
The numerical formulation as applied to the T equation is discussed here. The arrangement and nomenclature of layers is as shown:

The finite-difference centered-in-z form of the T equation is

In the code, the vertical differences are expressed in terms of the model's sigma coordinate through the hydrostatic equation as

Next, define

The general interior equation, for 1 < k < kx has the form

For the top layer, k=kx the boundary condition is K ( )=0, so the equation becomes

Note that water vapor is not carried above layer l2 in the present version of the model.

The lower boundary condition is

where T* is the surface temperature adjusted by the factor (pl/p*)R/Cp to yield a potential temperature gradient, and CH is a stability-dependent drag coefficient passed from subroutine PROGTN.For moisture diffusion, qs(T*) plays the role of T* with an additional factor DPHI multiplied in to account for soil moisture availability. For momentum, a different value, CM, is passed, and the terms in the u and v equations analogous to T* are zero.

The T equation for the lowest layer, k=l, is now

The upper and lower boundary equations now become, respectively,

Recombination of the terms yields, for interior, top, and bottom:

Introduction of and

and further manipulation gives equations for the tendency (Tn+1-Tn-1)/2 t:

If we define

the equations can be written in matrix form

where

and A is a tridiagonal matrix with main diagonal

The upper and lower diagonals are

Finally, collecting expressions for B,

This method of solution has been found to behave well under most conditions, but is occasionally prone to large-amplitude oscillations. A time-averaging scheme used to control these oscillations is described in Section 9.2.

The tendencies resulting from solution of the above system will be added on to the non-linear tendencies originally passed on from PROGTN after one further modification is made; since the model's history variable is virtual temperature, not temperature, it is necessary to convert the temperature tendency calculated here to a virtual temperature tendency. The expression used is

The final step is then the adding in of the non-linear tendencies, A , B, C, and D. These combined tendencies are passed back to the calling program, subroutine FIDI. Also passed back is a quantity for use in adjusting the soil moisture for evaporative loss.

1. Shallow Convection
This subroutine simulates the effects of shallow non-precipitating cumulus clouds by carrying out an enhanced vertical diffusion of specific humidity and heat in a manner similar to that described by Tiedtke (1983). Only model columns that contain conditionally unstable layers near the surface are eligible for this process, but, unlike the model's deep convection, no convergence of any kind is explicitly required here. This allows more vigorous vertical mixing of water vapor that would otherwise tend to accumulate near the surface in synoptically inactive regions.

The shallow convective process is treated in the same way as the model's basic vertical diffusion, acting merely to simulate the increased turbulence and attendant vertical eddy transport known to be a property of shallow cumulus clouds and their local environment. Once calculated, changes in T and q are not carried as tendencies as they are in MONIN3, the vertical diffusion subroutine, but instead are added as direct increments or adjustments.

As in MONIN3, the diffusion equation is used in the form

where T is temperature, q specific humidity, K (m2/sec) the diffusion coefficient, and the dry adiabatic lapse rate.

The enhanced vertical diffusion is brought about by the prescription of somewhat large values for K. The values were selected arbitrarily to produce reasonable-looking profiles of T and q, rather than those calculated from the Richardson number in MONIN3. The calculation follows these guidelines:

1.Cloud base is determined from values of lifting condensation level that have been received from the deep convection subroutine, which has been called just previously. A value for cloud top (the highest unstable layer) is also passed, but stronger restrictions (below) apply.

2.In columns for which deep convection has been performed, shallow convection is not allowed.

3.Shallow convective "clouds" may at most occupy sigma layers two through six, producing relatively strong diffusion within the cloud and weak diffusion above and below it, as governed by an imposed K profile, with -

4.Values of K are tapered in the vertical (to prevent development of unrealistic kinks in the T and q profiles) as follows: At the base of the cloud layer, K = 1.5 m2/sec; At the top, K = 1; For the next-to-top layer, K = 3; For intermediate layers, if any: K = 5.

5.There must be at least one potentially unstable layer within the conditionally unstable column.

SOLUTION

As in the diffusion calculation in MONIN3, the general interior equation is written

where

Solving for T

For the top layer, k = kx, k = 0, so = 0 and writing zero as ,

At the bottom, k = 1 and k = o, so = 0 and

As in MONIN3, these three equations form a tridiagonal system for layers 1 < k < 18, and the system is solved in the same way.

1. Deep Convection
Deep convection is modeled by a fairly basic Kuo-Anthes type of scheme, requiring moisture convergence and deep conditional instability in order to be active (Kuo, 1965, 1974; Anthes, 1977). Not considered are the ice phase , vertical fluxes of heat, momentum, and moisture, the effects of water loading, and the evaporation of falling rain below cloud base. Although the convective portion of a column is referred to here as the "cloud", no account is taken of cloud water.

The moisture convergence between cloud base and cloud top is partitioned into a rain-producing portion and a humidity-increasing portion on the basis of a function related to the column-integrated relative humidity. Finally , the convective heating and moistening of the environment is distributed with height on the basis of the vertical distribution of cloud-environment temperature and specific humidity differences. This adjustment is limited to layers between cloud base and cloud top.

Although in its current version the medium-range model is allowed to carry water vapor only in the lowest 12 sigma layers, this does not impose any ceiling on the deep convection itself; rather the only restriction is that n o convective moistening is allowed above layer 12. Air entering the convective column must come from either layer 2 (sigma = .981) or layer 3 (sigma = .960), or a combination of the two. Layer one is not considered, because with its thickness of only ten millibars it is too subject to the rapid variations in temperature and humidity produced by surface processes.

1. Outline of Calculations
The vertically-integrated moisture convergence is taken to be the mass-weighted integral of the change in specific humidity during the one leap-frog time step,

where pB is the pressure at cloud base and pT is the pressure at the top of the cloud. The change in q can come about through advection, numerical processes, or surface evaporation.

Several checks are performed to put realistic limits on the types of columns that are capable of supporting convection; these include requirements that the low-level temperatures be warm (T(2)>5oC), that no low-level inversions exist (T(2)>T(3)), that the total moisture convergence in layers 1-3 be enough to produce a rainfall rate of at least 2 mm/day, and that the buoyant layers in the column have a combined pressure thickness of at least 0.3 times the surface pressure.

1. Details of the Calculation
Subroutine VCTKUO performs the preliminary tests to determine which columns qualify for convection, then calls MSTADB which provides cloud base, cloud top, and in-cloud profiles of T and q, which are used to calculate the modification of the environment and the rainfall produced each time step.

In MSTADB, the necessary calculations are accomplished primarily by table lookups. First a table is set up to calculate the equivalent potential temperature, e, as a function of parcel temperature T, and pressure p, using the Rossby formula

where es is the saturation vapor pressure at T; R and Rv are the gas constants for dry air and water vapor; cp is the specific heat of dry air at constant pressure, and L is the latent heat of vaporization for water. Pressures are expressed in centibars. The saturation vapor pressure in the formula is calculated from an integration of the Clausius-Clapeyron equation,

where the subscript zero refers to a temperature of 273.16K, and cL is the specific heat of liquid water (which enters through the temperature dependence of L). The last preliminary step in MSTADB is the setting up of a table to recover values of q and T as functions of p and e along a moist adiabat.

The actual calculation of specific parcel ascent trajectories begins with the calculation of dew point at the parcel's origin. This is done by a call to subroutine DEWPOINT which sets up its own eS - vs. -T table and obtain s the dewpoint via an iterative process. Then the temperature and pressure at the lifting condensation are obtained from the formulas

where TD is the dewpoint. The tables in MSTADB are now used to generate a cloud temperature profile above p based on the computed value of e. The resulting vertical profile depends on the layer at which the lifted parcel originates, and so it is done twice, once each for layers 2 and 3. From these two cloud soundings, a composite is selected consisting of the warmest cloud temperatures at each level. The cloud top is considered to be the highest layer at which the parcel is buoyant (in terms of virtual temperature excess). Non-buoyant layers may exist between cloud base and cloud top.

The final cloud sounding is returned to VCTKUO and forms the basis for all the calculations that follow. First the differences between cloud and environment T and q at all levels are calculated. In layers for which the cloud is drier than the environment or cooler than the environment, the differences are zeroed. Then the vertically-integrated moisture convergence W is calculated. This moisture is divided into a precipitating and non-precipitating part (the closure of the scheme) according to the value of a measure of the column-mean capacity for additional water vapor,

where the overbars represent vertical averages in the environment over a depth corresponding to the cloud depth. The vertical distribution of convective heating and moistening are taken to be proportional respectively to the vertical distribution of the cloud buoyancy and the cloud-environment specific humidity difference:

where the subscript c refers to parcel (in-cloud) quantities. The quantities in brackets are not allowed to exceed unity. The increment of convective rainfall in m/sec is for the leapfrog timestep, but since the rain bucket is incremented at regular timesteps, a factor of 1/2 must be introduced in the code. The temperature increments, Tk, are added to the current ( + 1) value of Tk to complete this part of convective adjustment; for the specific humidity, the increments qk are added to the -1 values of qk because the convection has rained out a portion of the newly-converged water vapor, DQk. Up to this point in the calculation, both T and q had been assumed to have their -1 values.

The final adjusted values of Tk and qk are returned to the calling program, along with the values of the lifting condensation levels for all columns, for eventual use in the shallow convection subroutine.

1. Grid-scale Condensation
Subroutine LRGSCL is called after the deep and shallow convection subroutines. It is used to condense excess water vapor into precipitation, some of which is subsequently re-evaporated into unsaturated layers below. Starting with the highest moisture-bearing layer (k=12, in this version), the code proceeds downward a layer at a time checking for supersaturation. Layers having a relative humidity in excess of l00% ar e brought back by an approximate wet-bulb process to a just-saturated state, with all of the resulting liquid water assumed to go into rain and none into cloud or ice. Evaporation of falling precipitation is dependent on a parameterized drop-size distribution and upon relative humidity.

The passage from a supersaturated state at temperature T and specific humidity q, to a just-saturated state, T* and q*, requires the solution of the implicit equation

where q* is the saturation specific humidity at the final temperature. This is done with the help of the approximation

which becomes, on substitution of (5.1),

Solution for the water vapor condensed, , gives

where the Clapeyron equation has been substituted. The required adjustments are thus a decrease in the specific humidity by the amount and a temperature increase of . The condensed water P in units of rainfall depth (m) isfor p in Pa and g in m/sec (in the model code, the pressure is in centibars and the factor is unity).

If the next lower layer in the column is also supersaturated, the process is repeated and the rain P is augmented. If unsaturated, the layer is moistened and cooled by the evaporation of some or all of the falling rain, also according to a wet-bulb process. Here the change of phase is not allowed to occur instantaneously, but rather depends on the drop-size distribution. This must be parameterized as a function of liquid water content. The basis for the needed relations is provided by the work of Kessler (1969) where, starting with a Marshall-Palmer drop-size distribution, integration over drop size provides relations among rainfall rate, water content, and representative drop size, fall speed, and drop-air vapor exchange coefficient. The distribution probability density function for a unit volume of air iswhere N is the number of drops with diameter between D and D + dD, No is 107 m-4 and , after integration over all D, can be expressed as a function of total liquid water content M , asThe evaporation equation for individual drops is

where = the diffusion coefficient for water vapor in air (.226 x l0-4 m2 sec-1 at sea level), the are vapor densities, and G is a dimensionless drop ventilation coefficient. An approximation is used,

where, instead of the true vapor density gradient , between drop surface (a) and environment (b), we substitute a quantity proportional to the saturation deficit of the environment, . Integrated over all drop sizes, it gives an expression for the total evaporation due to all drops as a function of M, which is then expressed in terms of rainfall rate R with the help of an M-weighted representative drop fallspeed. The result is

From this we can derive the total precipitation evaporated during one time step,

where with t in seconds and the rainfall rate P/2 t in m/sec, POTEVP is in units of meters (depth of rain water), as is P. The semi-empirical dimensionless factor 320, contains the diffusivity of water vapor, standard values for pressure and air density, and correction factors to offset both the approximate vapor density gradient used in equation (5.3) and the exponent on the precipitation rate which in Kessler comes out to be 26/45, rather than the square root which we have substituted here for convenience. It is estimated that, as a result of t he vapor density gradient approximation which is most accurate at around 850 mb and 280K, evaporation rates are underestimated by about l5% at 700 mb and overestimated by up to 30% near l000 mb. The approximation to the exponent is most accurate for precipitation rates of about l mm/hr, underestimated by l0% at l0 mm/hr, and overestimated by 25% at 0.l mm/hr.

In the current version of the model, evaporation of falling precipitation is restricted from raising the relative humidity higher than 80%, except in the first layer, where it is halted at a relative humidity of 90%. If the precipitation increment P entering a layer exceeds POTEVP, all of POTEVP is subtracted from it. If not, all of the precipitation must evaporate into the layer. In either case, the evaporative cooling is calculated by use of the wet-bulb equations derived above, with the amount of evaporation now replacing .

The procedure is repeated layer by layer, and the net rainfall emerging from the bottom is halved to take into account the fact that the supersaturation has accumulated over a leapfrog time step, 2 t, while the precipitation is incremented each time step.

# References

Anthes, R. A., 1977: A cumulus parameterization scheme utilizing a one dimensional cloud model. Mon. Wea. Rev., 105, 270-286.

Blackadar, A. K., 1979: High-resolution models of the planetary boundary layer. Advances in Environmental Science and Engineering, Vol. I, J. R. Pfaflin and E. N. Ziegler, eds., Gordon and Breach.

Kessler, E., l969: On the distribution and continuity of water substance in atmospheric circulations, Met. Monographs, Volume l0, No. 32, American Meteorological Society, Boston, 84 pp.

Kuo, H. L., 1965: On formulation and intensification of tropical cyclones through latent heat release by cumulus convection. J. Atmos. Sci., 22, 40-63.

_____, 1974: Further studies of the parameterization of the influence of cumulus convection on large-scale flow. J. Atmos. Sci., 31, 1232-1240.

Tiedtke, M., 1983: The sensitivity of the time-mean large-scale flow to cumulus convection in the ECMWF model. Workshop on Convection in Large-Scale Numerical Models. ECMWF, 28 Nov-1 Dec 1983, pp 297-316.