Goddard Cumulus Ensemble Model. Part I: Model Description

During the past two decades, convective scale models have advanced sufficiently to study the dynamic and microphysical pro • êu1 • ê associated with mesoscale convective systems. The basic features of these models are that they are non-hydrostatic and include a good representation of microphysical processes. The Goddard Cumulus Ensemble (GCE) model has been extensively applied to study cloud-environment interactions, cloud interaction and mergers , air-sea interaction, cloud draft structure and trace gas transport. The GCE model has improved significantly during the past decade. For example, ice-microphysical processes, and solar and infrared radiative transfer pro-C • ü ses have been included. These model improvements allow the GCE model to study cloud-radiation interaction, cloud-radiation-climate relations and to develop rain retrieval algorithms for Tropical Rainfall Measuring Mission (TRMM). In Part I , a full description of the GCE model is presented, PP well as several sensitivity tests associated with its assumptions. In Part II (Simpson and Tao, 1993), we will review GCE model applications to cloud precipitating proce bmfòand to the Tropical Rainfall Measuring Mission (TRMM), a joint U.S.-Japan satellite project to measure rain and latent heat release over the global tropi VÛ .


INTRODUCTION .
During the past 20 years, observational data on atmospheric convection has been accu mulated from measurements by various means, including radars, instrumented aircraft, satel lites, and rawinsondes in special field observations (e.g., GATE, PRE-STORM, COHMEX, TAMEX, EMEX and several others).This has made it possible for convective cloud modelers to test their simulations against observations, and thereby improve their models.In turn, the models have provided a necessary framework to relate the fragmentary observations to help in understanding the complex physical processes interacting in atmospheric convective systems, for which observations alone still cannot provide a dynamically consistent four-dimensional TAO, Vo l.4, No .l, March 1993 picture.The past decades have also seen substantial advances in the numerical modeling of convective clouds and mesoscale convective systems (e.g., squall-type and non-squall type convective systems), which have substantially elucidated complex dynamical cloud environment interactions in the presence of varying vertical wind shear.With the advent of powerful scientific computers such as the CRAY-YMP/C90, many imponant and complex processes (which require extensive computations), such as ice-microphysics and radiative transfer, can now be simulated to a useful (but still oversimplified) degree in these numerical cloud models.
A basic characteristic of the convective cloud models is that their governing equations are non-hydrostatic since the vertical and horizontal scales of convection are similar.Such models are also necessary in order to allow gravity waves, such as those triggered by clouds, to be resolved explicitly.Another requirement of these models is that they need increasingly sophisticated representation of microphysical processes, to which the dynamic processes have shown to be highly sensitive (Tao et al., 1991;Ferrier et aL, 1992).The Goddard Cumulus Ensemble (GCE) model is one such model.Others have been developed by Klemp and Wilhehnson (1978), Cotton and Tripoli (1978), Schlesinger (1978), Clark (1979), Chen and Orville (1980), and Chen (1991).
The GCE model has been used to provide essential insights into the interactions of clouds with each other (Tao and Simpson, 1984;Tao and Simpson, 1989a).with their surroundings (Soong and Tao, 1980;1984;Tao and Soong, 1986;Tao and Simpson, 1989b), with long wave radiative transfer processes (Tao et al., 1991(Tao et al., , 1993)), with ocean surfaces (Tao et al., 1991;Lau et al., 1993;Sui et al., 1993a and b), and with trace gas distributions (Scala et al., 1990;Pickering et al., 1991Pickering et al., , 1992a, b, c), b, c).The GCE model has also •been used to convert the radiances received by cloud-observing microwave radiometers into predicted rainfall rates (Simpson et aL, 1988).Remote sensing of cloud-top properties by high-flying aircraft bearing microwave and other TRMM instruments have already begun to provide a new and powerful means of testing the GCE model, particularly when such observations are augmented by simultaneous ground-based radar measurements (Adler et al., 1988(Adler et al., , 1991;;Yeh et al., 1992).
The Goddard Cumulus Ensemble Model's framework was originally based on the cloud ensemble model developed by Soong and Ogura (1980), Soong and Tao (1980), Tao (1983) and Tao and Soong (1986).However, many significant improvements have been made during the past decade at Goddard (it was named the GCE model in 1989).In this paper, we will present the full formulations (including the improvements) used in the GCE model, as well as several sensitivity tests associated with the assumptions.In Part II (Simpson and Tao, 1993), we will review GCE model applications to the study of cloud precipitating processes and to their measurement from space.

Equations of Motion
The derivation of the momentum equations requires the use of the moist equation of state, which is commonly written in the form p = pRT(l + 0.61qv) (1 where p is the pressure, p the density of moist air, R the gas constant for dry air, T the temperature and q v the mixing ratio of water vapor.The non-dimensional function of pressure is defined as (2 where P0 is a reference pressure (1000 mb) and Gp the specific heat of dry air at constant pressure.The virtual potential temperature ( 811) is defined as The potential temperature () is related to 7r by (} = T / 7r.where u and v are the horizontal and w the vertical velocity components, g the acceleration of gravity and qi the sum of the mixing ratios of liquid and ice water.A prime is defined as a deviation from the horizontal average.The subgrid-scale diffusion effects are represented by Du.D v and Dw which will be defined in Section 2.4.
The equations for(} and qv are written in the form --uB --v{} --:: -p w B + D9 + -( c -e c -er)+ -(f-m) Ls( aq., a a 1 a _ -= --uq --vq ---pwq + D + ( c-e -e) + (d-s ) ( 8) where L v. L f and L 8 are the latent heats of condensation.fusion and sublimation, respec tively.The variables c, Ee, Er, f, m, d ands stand for the rates of condensation, evaporation of cloud droplets, evaporation of raindrops, freezing of raindrops, melting of snow and grau pel/hail, deposition of ice particles and sublimation of ice particles, respectively.The term QR is the cooling/heating rate associated with radiative processes and will be discussed in Section 2.3.The subgrid-scale diffusion terms D9 and Dq., will be defined in section 2.4.

Cloud Microphysical Processes
The cloud microphysics include a parameterized Kessler-type two-category liquid water (0.0004 cm-4 ), respectively.The density for rain, snow and graupel (hail) are 1 g cm-3 , 0.1 g cm-3, and 0.4 g cm-3 (0.917 g cm-3), respectively.The cloud ice has a single size (mono-disperse) where its diameter and density are assumed to be 2 X 10-3 cm and 0.917 g cm-3, respectively.The values of these parameters are based on limited observational data.
The predictive equations for cloud water (qc).rain (qr).cloud ice (qi).cloud water and cloud ice will move with the surr ounding air.The terms T� c • T q r • Tq i > T q8 and Tq g are microphysical transfer rates between hydrometeor species, and their sum is zero.They are defined as: T q i = -(Psaut + Paaci + Praci + Psji + D9aci + Wgaci) + Pihom -Pimlt + Pidw (16) For T > 273.16°K (20) h2 is defined as 1 for a grid box which has qr and q8 < 1 X 10-4 9 9 -1 , and otherwise is defined as zero.63 is defined as 1 for a grid box which has qr < 1 X 10 -4 9 9 -1 • and otherwise is defined as zero (see Lin et al., 1983).Dgaci • D9acr and Dgacs ( Wgaci . W9acr and W9ac s) are production rates for diy (wet) growth of hail.These microphysical processes are demonstrated in Figure 1 and explained in    A saturation adjustment scheme that calculates the amount of condensation (and/or deposition) necessary to remove any supersaturated vapor, or the amount of evaporation (and/or sublimation) necessary to remove any subsaturation in the presence of cloud water (cloud ice) is needed for a non-hydrostatic cloud model (e.g., Soong and Ogura, 1973;Klemp and Wilhelmson, 1978;Cotton and Tr ipoli, 1978;Clark, 1979;Lin et aL, 1983).One approach frequently used is a relaxation technique (e.g., Newton-Raphson method) to iteratively balance the heat exchange and change of phase of water substance.A water phase only saturation adjustment without a need for iterative computation was fir st proposed by Soong and Ogura (1973).A modified version of their scheme with the inclusion of an ice-phase was presented by Tao et al. (1989).Two major assumptions involved in this water-ice saturation adjustment scheme are the following: 1) The saturation vapor mixing ratio is defined as a mass-weighted combination of the saturation values over liquid water and ice between 0°C and -40°C.into cloud water/cloud ice instantaneously when mixed-phase supersaturation occurs.The water/cloud ice is asswned to evaporate/sublimate instantaneously if the air is subsaturated.
Precipitating rain drops (snow and graupel) evaporate (sublimate) only if the cloud water (cloud ice) is exhausted before saturation is attained.In application of this saturation ad justment, the terms, Pihom• Pidw and Pimlt• will be activated first.In addition, initiation of cloud ice (Pint) and depositional growth of cloud ice (Pde p i) discussed in Rutledge and Hobbs (1984) will be used to initiate the cloud ice in a saturated environment.This proce dure weighs the saturation mixing ratio in favor of ice at levels above the freezing level 0°C.This adjustment scheme will almost guarantee that the cloudy region (defined as the area which contains cloud water and/or cloud ice) is always saturated (100% relative hwnidity).
This permits subsaturated downdrafts with rain and haiVgraupel particles but not cloud-sized particles.Tao et aL (1989) presented the derivation (sensitivity) as well as the performance of the ice-water saturation adjustment scheme.

Radiative Transfer Processes
Short-wave (solar) and long wave (infrared) radiation parameterizations have recently been included in the GCE model.The long wave and short-wave radiative transfer calcula tions implemented in the GCE model were developed by Chou (1984Chou ( , 1986) ) and Chou and Kouvaris (1991).This radiation scheme has also been applied in several General Circulation Models (e.g., the UCLA-GLA GCM and the FSU global model).In computing the absorption of solar radiation by atmospheric gases, the spectrum is divided into two regions: short-wave ( < 0.69 µm) and near infrared ( > 0.69 µm ).In the near infrared region, the computation of absorption of solar radiation by water vapor follows that of Chou and Arking (1981) and Chou (1986), which uses a broad band transmission parameterization for clear atmospheres and k-distribution approach for cloudy atmospheres.In the short-wave region, the compu tation of the absorption by ozone follows Lacis and Hansen (1974).The absorption by 02 and C0 2 near the visible and near IR regions is included using the parameteriz<1-tion of Chou (1992).In cloudy atmospheres, the four-stream discrete ordinate of Stamnes et al. (1988) is used to compute transmittance and reflectance of cloudy layers.The addirlg method is then used to compute fluxes in a composite of clear and cloudy layers.The long-wave scheme takes into account the effects of both gaseous and hydrometeor absorption.For gaseous absorption, the effects of water vapor, carbon dioxide and ozone are computed.
It is necessary to parameterize the emissivity of gases and hydrometeors in order to simulate long-wave and short-wave radiative transfer processes.The parameterizations for water vapor and cloud effects are based on the commonly used broad-band emissivity method (see Stephens, 1978Stephens, , 1984)).Both liquid and solid phases of water as well as their size distributions are considered.In the water phase, the optical depth ( r) and effective emissivity J.00 n(r)r3dr r e = "'-0"-oo�--- fo n(r)r 2 dr (24) (25) where W1 is the liquid water content (g m -3 ).For the ice-cloud, the optical depth and effec tive emissivity of cloud are parameterized following Starr and Cox (1985) and Ramaswamy and Ramanathan (1989), giving the relations: and, e!}1 = 1 -exp( .... .bo ttwi) b0 l = 0.05; bo where Wi is ice water path (g m-2); is the ice water content (g m-3); RR is 1.6, 1.4 and 1.2 for the tropics, mid latitude day and mid latitude night, respectively.The effective radius of cloud water and cloud ice are specified to be 0.0003 and 0.002 cm, respectively.Model-predicted radiative cooling and heating rates at cloud top are on the order of 30 to 50°K/day, which is in good agreement with Stephens (1978).A quantitative estimate of the impact of radiative processes on the total surface precipitation as well as on the development and structure of squall lines will be presented in Part II.

Subgrid-scale Turbulence
The subgrid-scale turbulence used in the GCE model is based on work by Deardorff (1975), Klemp and Wilhelmson (1978), and Soong and Ogura (1980).In their approach, one prognostic equation is solved for subgrid kinetic energy (E), which is then used to specify the eddy coefficients (J( m).The effect of condensation on the generation of subgrid-scale kinetic energy is also incorporated in the model.The equation for the subgrid kinetic energy is written as 1 ax j axj m axi The terms on the right-hand side of equation ( 28) are stability, shear, diffu sion and dissipation effects, respectively.The double overbar in the equation denotes an average over each grid point in the cloud model and double primed variables are the deviation from the grid values.
The variables .6.x, .6.y and .6.z are the grid int ervals .The higher value of Cm and the lower value of Ce imply that stronger mixing will occur.The values Cm=0.2 and Ce=0.7 are used in the GCE model (Deardorff, 1975).
All the subgrid-scale terms have been defined except for the vertical buoyancy fluxes in a saturated area.These terms in a saturated area can be expressed as a function of the nearly conservative equivalent potential temperature 8 e and the total mixing ratio of water substance qt = qv + qc + qi + qr + qs + q9'.Following Klemp and Wilhelmson (1978), In an unsaturated area, 0 11 1 8 0 a w "(-+ 0 61q" -q f ' ) = -Kh("""' -+ 0 61 q 11 ) (J • 11 Equation ( 31) is applied to both the cloud and subcloud layers wherever the air is not saturated.
This term creates turbulent kinetic energy in the subcloud layer where the surface heat flux and evaporation maintain a slight instability.In the unsaturated portion of the cloud layer, this term produces a strong damping effect and destroys turbulent kinetic energy rapidly.The ratio K h /Km =3 was used in Klemp and Wilhelmsen (1978).They also used the same K h in both the horizontal and the vertical directions, implying an isotropic diffusion process.
However, the ratio of /( h / /( m should actually decrease with increasing stability (Deardorff, 1973 ).Soong and Ogura (1980) incorporated the stability dependence of the ratio /( h /Km and the anisotropy of turbulence by setting K h =2K m everywhere except [( h =3K m in the vertical direction only at grid points where the sounding is unstable.Once the E, K h and Km are determined, the subgrid-scale terms can be computed from oz q c/i :::: :: : oz w q + fJ z h fJ z The derivation of these equations can be found in Klemp and Wilhelmson (1978).
2.5 Heat and Moisture Fluxes from the Ocean surface ( 33 ) Surface heat and moisture are included in the diffusion tenns Da and Dqv at the lower boundary of the model obtained from bulk formulations.The vertical flux of sensible heat from the sea surface is assumed to take the form ( 34 ) where C v is the drag coefficient.V0 the surface wind speed, Ts the sea surface temperature.and T0 the air temperature at the anemometer level.The GCE model uses the formulation by Roll (1965) for the drag coefficient, namely: where V0 is in meters per second.The moisture flux from the sea surface has a similar form ( 3 6) where q 0 is the mixing ratio at the anemometer level and q8 is the saturation mixing ratio at the sea surface temperature.fu the model computations, the values of V0• T0 and q 0 are taken from the lowest grid level.

Model Domain, Numerical Technique, and Boundary Conditions
A stretched vertical coordinate (height increments from 220 to 1050 m) with 31 grid points is used in order to maximize resolution in the lowest levels.The model top can be 20-25 km.For a typical two-dimensional simulation, 612 grid points were used in the horizontal, the central 504 of which comprise the fine-grid area with a constant 750 m resolution.Simulated cloud activity is confined to the fine-grid region through a Galilean transformation, namely by subtracting the storm propagation speed from the initial wind field.Outside of this region, the grid is horizontally stretched with a ratio of 1.0625: 1 between adjacent grid points.This results in a domain that is about 1025 km wide.A staggered grid arr angement (Arakawa C grid) is used in the model.A leapfrog time integration and second-order space derivative scheme in the vertical direction are used.Either a second er a fourth-order space derivative scheme (Chen, 1980) is used in the horizontal direction.
To avoid the problem of time�splitting.a time smoother is adopted.The time smoothing coefficient is set to 0.1 with a time interval of 7.5 s.
Open-type lateral boundary conditions are most frequently used (Klemp and Wilhelm sen.1978).The goal was to allow the advection of the normal velocity component out through the boundary at a speed given by the sum of the nonnal velocity component and an approximated inttinsic phase velocity c* of the dominant gravity wave modes with mini mum reflection.An approximation for the nonnal velocity equation at an outflow boundary location is given by &u &u where c* must be selected in some way.Klemp and Wilhelmson (1978) suggested the c* might be about 30 m s-1.Chen (1982), however.found that in some cases a c* less than 30 m s-1 leads to a better comparison (small domain versus larger ones), particularly in some three-dimensional tests.The use of stretched horizontal coordinates causes the model results to be less sensitive to the choice of gravity wave speed associated with the open lateral boundary conditions (Fovell and Ogura, 1988).A free-slip boundary condition is used for u, v, (} , qv and all five water hydrometeors.The w-velocity is set to zero.A 5 km deep Rayleigh relaxation (absorbing) layer is also used at the top of the model.
In the three-dimensional model, however, the current computer limits us to about 218 X 54 grid points (using 1 or 1.5 km resolution) in the horizontal and 32 grid points (to a depth of 22 km) in the vertical.For line-type convection, an open-type lateral boundary condition is used along the x-axis in the direction perpendicular to the line convection and a periodic lateral boundary condition is applied along the y-axis.However, open lateral boundary conditions can also be used on the y-axis for simulating isolated stonns (e.g., super-cells).In actual computation, the number of grid points in each simulation can be adjusted for specific purposes and economy of computer time.

SENSITIVITY TESTS
3.1 Microphysical Processes

Ice versus no-ice runs
To investigate the role of ice-phase microphysics in stratiform precipitation, a no-ice phase version of the GCE model, with a Kessler-type two-category liquid water microphysics.was utilized for a GATE squall line simulation (Tao and Simpson, 1989b).The system's life cycle and total precipitation do not differ significantly from the ice run (see Table 2).The modeled propagation speed of about 8.2 m s-1 differs only slightly from that of the ice run.These results suggest that the ice-phase microphysical processes do not significantly affect the propagation speed of the simulated systems.From a close examination of the model results.it was found that the propagation speed of the modeled systems is closely related to the proximity and vigor of simulated downdrafts and their associated cold outflow structure in the lower troposphere.The total amount of evaporative cooling in the stratiform and convective regions in the ice-free-run does not change signifi cantly from the run with ice in these particular GATE simulations.The small difference (< 1 to 2 m s-1) in simulated squall line propagation speeds for the runs with and without ice-phase processes were also reported in Yoshzaki (1986), Nicholls (1987), Fovell and Ogura (1988) and Chen (1991).In each study.the surface wind gust, temperature drop, and pressure rise, which are all usually associated with cold outflow, are quite similar between the runs with and without ice-phase processes.
Table 2.Estimated surface rainfall and rainfall area at the surface.Rain is given in relative units (millimeters per grid point accumulated over a 12 h period).The total rain area is given as a percentage of the total domain area Anvil portion and anvil rain area are given as a percentage of the total rain and total rain area, respectively.

47
Several differences, however, are noted between the ice run and the ice-free run.The first difference concerns the precipitation statistics of simulated squall MCSs .It was found that light rain (< 10 mm h-1) only accounts for about 26.5 % of the total rain, but it covers 90% of the total rain area for the ice run.By contrast, heavy precipitation (> 30 mm h-1) accounts for a large portion of the total rain, but it only occupies a very small portion of the total rain area.Without the inclusion of ice-phase microphysics, the role of heavy precipitation is increased significantly.Also, only 12% of the rain is characterized as stratiform (Figure 2).The model results with ice-phase microphysics are generally in good agreement with GATE observations (see Table 1 in Ta o and Simpson, 1989b).Another difference is that the depth of the stratiform cloud is reduced with ice-free microphysics.More cellular rather than mesoscale features are also evident for the ice-free run (Figure 3).In a separate modeling study of a subtropical squall line that occurred during TAMEX (Tao et al., 1991), it was found that a slower propagation speed resulted from the run without ice processes.However, the simulated squall system did not last long (decayed about 5 h into the simulation).The rain fell into the front inflow region and cut off the moisture supply in the ice-free run.The updraft titled downshear in relation to the low-level shear.As reported in many modeling and observational studies, an upshear tilt of the updraft is a necessary condition for a long lasting convective system.

Ice scheme comparisons
The GCE model has been used to investigate different ice-phase microphysical pa rameterization schemes as well as their effects on convective motion under different tropical environmental conditions (McCumber et al., 1991).Two specific experiments from Mccum ber et aL will be presented here.The first one adapted the Lin et aL (1983) scheme which includes hail processes.Intercept parameter for snow (Nos), hail (Noh) and density of hail (pg) are 0.03 cm-4 , 0.0004 cm-4 and 0.917 g cm-3, respectively.These parameters were typical for mid latitude convective clouds.In the second experiment.Nos, Noh and pg are 0.04 cm -4, 0.04 cm-4 and 0.4 g cm-3, respectively, which were used for warm frontal cloud simulations (Rutledge and Hobbs, 1984).The graupel has a low density and a large intercept (i.e., high number concentration).By contrast, the hail has a high density and a small intercept (low number concentration).These differences affect not only the description  of the hydrometeor population, but also the relative importance of the microphysical processes (Figure 4).For example, in the hail case, the value selected results in a few large, fast-falling hail stones.These characteristics confine hail to the convective cells.A result of this is that snow becomes the dominant precipitating hydrometeor type in the anvil cloud.The large amount of snow (Figures 4a and 5b) is the reason why melting and deposition of snow are first-order processes for the hail case.In the graupel case, however, the selected parameter values result in a large number of graupel particles that are smaller than the hail particles and thus fall more slowly.These characteristics permit graupel to occur in both the convective cell and the anvil.Also note that in the anvil graupel is much more prevalent than snow (Figure 5a).For this reason, melting and deposition of snow are second-order processes.For the tropical squall-type convective systems, the principal distinction in the heating profiles between the hail and graupel case is the shape and intensity of the heating/cooling component in the anvil region (Figure 6).The cooling in the anvil region from the freezing level at about 4.5 km down to about 2 km in Figure 6a indicates the significance of melting graupel particles.Significantly less cooling due to melting in the anvil region for the hail case (Figure 6b) as opposed to the graupel case suggests that for the tropical squall system, the optimal mix of three-class ice hydrometeors is cloud ice-snow-graupel.A similar conclusion was reached for a TAMEX squall line simulation (Tao et aL, 1991).Recently, Sui et al. \ ... ..... .. .Another characteristic that is sensitive to the parameterization of cloud microphysics is the passive microwave signatures that result from the cloud hydrometeors (e.g., Figure 5).Simpson et al. (1988) and Adler et al. (1991) showed that vertical profiles of water and ice particles predicted by the cloud model can be used in a passive radiative transfer model to calculate the black body temperatures corresponding to signals detected by remote aircraft or satellite-borne sensors, and could be done simultaneously for several important microwave frequencies.Their response is particularly sensitive to the concentration and properties of large ice (Smith and Mugnai, 1989).Since the GCE IIJOdel is used in rain and latent heat retrieval algorithms based on passive microwave (Adler et aL. 1989).it is very important to properly represent the vertical profile of hydrometeors.To appreciate this point, consider again the vertical hydrometeor profiles in Figure 5.The variability that they display as a function of different microphysical formulations provides some hint of the impact of different ice schemes on passive microwave characteristics, which is discussed in more detail in Part II.

Anelastic and Compressible Syste ms
The Goddard Cumulus Ensemble (GCE) modeled flow can be either anelastic (Ogura and Phillips, 1962), filtering out sound waves, or compressible (Klemp and Wilhelmson, 1978), which allows the presence of sound waves.The sound waves are not important in thermal convection, their processes can place severe restrictions on• the time step in numerical integrations because of their high propagation speed.For this reason, most cloud modelers use an anelastic system of equations in which sound waves have been removed by eliminating certain terms in the compressible system (e.g., neglecting the local variation of air density with time in the mass continuity equation).
In the anelastic approximation the continuity equation is: (38) A diagnostic pressure equation can be determined for the anelastic set by differentiating (4) with respect to x, (5) with respect to y and (6) with respect to z.Because the boundary conditions for the pressure equation are of the Neumann type (i.e., they involve pressure derivatives), there are an infinite number solutions to the elliptic pressure equation that differ from one another by some constant.This is not a problem because only the derivatives of pressure are needed in the model equations (4-6).The elliptic pressure equation can be solved using direct (e.g., FFT) or iterative methods.However, when terrain is added, the direct methods become more complex and sometimes time consuming.
In the compressible system, the pressure equation is derived by taking the substantial derivative of (2) and using the compressible continuity equation op a -+ -p u ,=0 at axj 1 (39) � eliminate !jtf and the thermodynamic equation (7) to remove �:.The resulting equation where c is the spee d of sowid given by c 2 = cpRdiriiv/cv.Due to the presence of sowid waves, a very small time step, 2 s, for a 1000 m spatial resolution is needed for time integration of the entire model equation set However, Klemp and Wilhelmson (1978) developed a semi-implicit time-splitting scheme, in which the equations are split into sowid wave and gravity-wave components, to achieve computational efficiency.The small time step integration is semi-implicit in the vertical using a 2 s time step.The remaining time integration can use a 10 s time-step.One advantage of the compressible system is its computational simplicity and flexibility.The numerical code then remains a set of explicit prognostic equations and alterations involving such things as stretched or nested grids, surface terrain and bowidary conditions (e.g., radiative upper boundary) can be incorporated into the numerical model without complicating the solution procedure (Klemp and Wilhelmson, 1978).Anderson et al. (1985) have tested an anelastic system using a 4 s time step against the results from a fully compressible system without using the time-splitting technique (which needs a 0.3 s time step).They have found that the anelastic system produced essentially identical results to those of the compressible system for 2-D cool pool experiments lasting 500 s, where the convection was initiated by a cool pool.Ikawa (1988) has also compared the anelastic and the compressible system with a 2-D case involving orography.His results indicated that both simulated systems are similar if sound waves are damped enough in the compressible system.Due to the stretched grid arr angement in the horizontal (Section 2.6), the three-dimensional version of the GCE model (which has open lateral boundary conditions and a stretched grid) adapted the compressible system.All of the two-dimensional GCE simulations (including published results) used the anelastic system.A pair of sensitivity tests, anelastic versus compressible, were perfonned with an initial condition ass ociated with a PRE-STORM squall line (Figure 1 in Tao et al., 1993).All model physics are activated in these two runs for a 12 h time simulation in order to maximize the possibility of differences.Figure 7 shows that the GCE model does not produce identical results.The differences between the anelastic and compressible systems are much smaller, however, than those obtained by changing microphysical processes and advection schemes (Section 4.2).
The precipitation statistics (e.g., total precipitation and its stratiform component) from the anelastic and the compressible systems are also quite similar especially during the first 4 h (Table 3).Several observed features for this PRE-STORM squall system, such as the propagation speed of the squall system, the weakly evolving multicellular structure, the mesolow aloft, the squall pressure high, and the rear inflow from the midtroposphere, were well-simulated in both systems.No significant differences were found in simulated cloud updraft/downdraft structure and cloud top height between the two systems.It is difficult to tell which method is the most promising or will prevail in the future.and flow of the environmental wind around a cloud.However, doubling the number of grid points in a three-dimensional domain implies an eight fold increase in computer time, which exceeds the computing power currently available.Mesoscale convective systems (e.g., squall lines) can cover a very large area (300-400 km in width and 1000 km in length).Both two and three-dimensional simulations thus need to be run in combination.
Tao et aL (1987) used 2-D and 3-D versions of the GCE model to study the statistical properties of cloud ensembles for a well-organized intertropical convergence zone (ITCZ) rainband that occurred during GATE.The statistical properties of clouds, such as mass flux by cloud drafts and vertical velocity as well as condensation and evaporation have been examined.Figure 8 shows heating rates by condensation (c) and evaporation (e) in the 2-D model and the 3-D model.The heating rate estimated from large-scale observations, Q 1 -QR• and the total cloud heating rate in the 3-D model are also included.The rate of condensation in the 3-D case is slightly larger than its 2-D counterpart, but so is the rate of evaporation.As consequence, the net heating effect of clouds in the 3-D model is nearly equal to the 2-D counterpart, and both of them agree with Qi -QR as estimated from the large-scale heat budget (Note a minimum in the observed heating profile near 700 mb as discussed in Section 3.1).100 r---------.------------Zipser and LeMone (1980) and LeMone and Zipser (1980) presented the results of statistical analyses of convective updrafts and downdrafts.Their analysis is based on aircraft data gathered from cumulonimbus cloud penetrations for six days during GATE.In order to facilitate a comparison between our model results and their analysis, we have subdivided the updrafts and downdrafts into active or inactive updrafts and downdrafts (Table 4).For example, a grid point in the model is designated as an active updraft region if (a) the total liquid water content exceeds 0.01 g kg-1 and (b) the vertical velocity is larger than 1 m s-1 (or 2 m s-1 , depending upon how we defi ne "active") at that grid point and at that integration time.The inactive and active cloud fractional area coverage's as determined by aircraft measurement and model results are shown in Table 5.These ratios between active cloud updrafts and downdrafts indicate an excellent agreement among results between 2-D and 3-D models as well as the cores as measured by Zipser and LeMone (1980).Both 2-D and 3-D model results also showed a similar feature in that the active updrafts account for approximately 75% of the upward mass flux due to clouds and yet they only cover about 12-14% of the total area.This result is consistent with the concept, first proposed by Riehl and Malkus (1958; see also Riehl and Simpson,1979), that hot towers play a critical role in the heat and moisture budgets in the tropics, even though they occupy a small fraction of the area.Overall, our comparison study has indicated that the statistical properties of the clouds obtained in the 2-D model are essentially the same as the 3-D counterpart given an identical large-scale environment (see Tao , 1983).
In another GATE squall line simulation, Tao and Simpson (1989b) found that the three dimensional model results could not adequately reproduce many of the features associated with the stratiforrn region due to a small domain as well as a shorter time integration (4h).
The orientation of the three-dimensional simulated convective line is perpendicular to the environmental wind shear as observed during GATE.Both 2-D and 3-D modeled propagation speeds are in fa ir agreement with observational studies.Based on the trajectory analyses, it was also found that the air circulation near the leading edge of a squall-type system is very similar in both the two and three-dimensional model simulations.This conclusion is consistent with the results of Rotunno et al. (1988) and Moncrieff (1992).However, all of these conclusions were based on simulations of well-organized convective lines.A new microphysical parameterization has been developed recently by Ferrier (1993) and Ferrier et al. (1993), which combines the main features of previous three-class ice schemes by calculating the mixing ratios of both graupel and frozen.drops/hail.Additional model variables include the number concentrations of all ice particles (small ice crystals, snow, graupel and frozen drops), as well as the mixing ratios of liquid water in each of the precipitation ice species during wet growth and melting for purposes of accurate active and passive radiometric calculations.The parameterization also includes the following: (1) more accurate calculation of accretion processes, including partitioning the freezing of raindrops as sources of snow, graupel and fr ozen drops/hail; (2) consideration of rime densities and riming rates in converting between ice species due to rapid cloud water riming; (3) incorporation of new parameterizations of ice nucleation processes, the rime splintering mechanism using laboratory data, and the aircraft observations of high ice particle concentrations; (4) shedding of liquid water from melting ice and from excessive amounts of water accumulated on su percooled frozen drops/hail; (5) preventing unrealistically large glaciation rates immediately above the freezing level by explicitly calculating freezing rates of raindrops and freezing rates of liquid water accreted onto supercooled ice; (6) introducing fall speeds and size distribu tions for small ice crystals; (7) calculating radar reflectivities of particles with variable size distributions and liquid water coatings from Rayleigh theory; (8) basing conversion of particle number concentrations between hydrometeor species on preserving spectral characteristics of particle distributions rather than conserving their number concentrations (important).Details of these parameterized processes can be found in Ferrier (1993).
The three-class ice formulations perform well only in specific large-scale environments, such that tropical-maritime storms are better represented using the Lin et al. (1983) scheme as modified by Rutledge and Hobbs (1984) which substitutes slower falling graupel for hail (McCumber et aL, 199 1;Tao et aL, 1991).whereas the unmodified Lin et aL scheme is preferred in simulating midlatitude-continental storms (Ferrier et al. , 1991;Tao et al., 1993).
Furthermore, reasonable agreement with the observations may be obtained only after trial and error adjustment of several tunable coefficients in these three-class ice schemes.Figure 9 shows the radar reflectivities of an intense mid latitude-continental COHMEX (COoperative Huntsville Meteorological Experiment) squall line and a mature tropical-maritime GATE squall line simulated using the 4ICE microphysics.The simulated reflectivities compare well against radar observations of the COHMEX storm (Adler et al., 1990;Yeh et al., 1990) and the GATE squall line (Zipser et al., 1983;Szoke and Zipser, 1986; not shown here).
The 4ICE parameterization is much improved over the three-class ice schemes in the ability to simulate the radar and microphysical structures of storms that develop in vastly different environments with minimal tuning.

Numerical Advection Schemes
A second-or fourth-order horizontal advection scheme was generally used in the model with a time step of 7.5 or 10 s.However, using second-order or higher-order accuracy ad vection schemes can introduce some difficulties because negative values arise in the solution (Soong and Ogura, 1973).This effect can be especially important in cases where the solution of the advection is used as input to nonlinear equations describing microphysical phenomena or inert tracers which can eventually lead to instability of the whole system (Smolarkiewicz, 1983).The use of upstream differencing or other low-order schemes (Soong and Ogura, 1973) would not produce dispersive "ripples" but would suffer from excessive numeric al diffusion.
Smolarkiewicz (1983) has reduced the implicit diffusion by using a second "upstream" step where a specifically defi ned velocity field leads to a new form of a positive defi nite advec tion scheme with small implicit diffusion.The positive definite advection scheme involves iterations and needs more computational resources.Th is scheme has been improved with mul tidimensional application (Smolarkiewicz, 1984) and a non-oscillatory option (Smolarkiewicz and Grabowski, 1990).
The 2-D version of the GCE model has tested this Multi-dimensional Positive Definite Advection Transport Algorithm (MPDATA) for a PRE-STORM squall line case.Three spe cific tests have been performed.The first one was that all scalar variables (8 , qv, J( m and all five hydrometeor classes) used forward time differencing and the MPDATA for advection.The dynamic variables, u, v and w, used a second-order accurate advection scheme and a leapfrog time integration (kinetic energy semi-conserving method).In the GCE model code, the perturbed potential temperature was predicted and it could be negative.In application of MPDATA, a constant (large number) needs to be added on to 8'.This can cause some prob lems (e.g., invariance and stability) as discussed in Smolarkiewicz and Clark (1986).Thus, the second test applied the MPDATA method to the water vapor and hydrometeors only.The third test used the leapfrog time integration and fourth-order space derivative scheme for scalar and dynamic variables.Table 6 shows total surface precipitation and its anvil portion fr om these three runs.Several observed fe atures for the PRE-STORM squall system, such as the weakly evolving multicellular structure, the mesolow aloft, the squall pressure high, and the rear infl ow from the midtroposphere, were simulated in all three runs.All three runs also indicate a gradual growth of anvil precipitation with time (Table 6).However, the runs with the positive definite advection scheme generated more stratiform rain as well as broadening The COHMEX storm propagated eastward (to the right) and had a large leading (pre-line) anvil, whereas the GATE system propagated westward (to the left) and had an extensive area of trailing stratiform precipitation.its horizontal extent (Figure 10).These features are in better agreement with observational studies (Johnson and Hamilton, 1988;and Ruwedge et al., 1988).The non-oscillatory option was used in the first runs, but oscillation in temperature still occurr ed in the first two runs.MPDATA as well as other positive-defi nite advection schemes will be tested for the new double-moment ice formulation in the near future.

Modeling the Interactions . between Cloud Ensembles and Regional-scale Circulations
Observations of cumulus clouds generally indicate that clouds do not exist as isolated entities, but occur in groups or clusters.Observations also show that low-level mesoscale convergence (on the scale of 30-300 km) in a conditionally unstable environment is gener ally necessary to generate ensembles of deep convection (Ogura et al, 1979; Simpson et al.,   1980).In general , numerical cloud model simulations have been initialized with an atmo sphere having zero large-scale vertical velocity and with thermodynamic and horizontal wind profiles that are typical of the conditions upstream of severe convection (e.g., Klemp and Wilhelmson, 1978;Schlesinger, 1978).Then, a perturbation is introduced by either a positive temperature anomaly at low levels (the "warm bubble" method) or by cooling the low�t lay ers (the "cold pool" method).Several different methods of inducing large-scale convergence have been described in the literature (Soong and Tao, 1980;Chen and Orville, 1980;Tr ipoli and Cotton, 1980;Tao , 1983;Dudhia and Moncrieff, 1987;Crook andMoncrieff, 1988, Tao andSimpson, 1989b;and others).One difficulty with these methods is that the large-scale convergence can only be evaluated before hand (a priori condition) from either observational or theoretical studies.It is important to note that in convective conditions, moist-convective processes can feed back (both positively and negatively) into the convergent forcing.For a model to approximate the real atmosphere, it sh ould provide a convergence field that can vary in time and in space and also allow a fe edback between the convection and the large-scale convergence field.In order to include this effect, either an eventual two-way nesting between a convective-scale and a mesoscale model or a much higher resolution capability (through an interactive grid nesting developed by Clark and Farley, 1984;and Chen, 1991) is needed.
The interaction between a mesoscale (GMASS, Goddard Mesoscale Atmospheric Simu lation System) model and our convective scale (GCE) model is underway.We have modified   -� -. .. .,� .� .::: -r->-r--r-....-.--r"'T:""".,...,--r-'l'.::--O-r--r--rb -r-i 1 4 200 -, , / < :"   the GCE model to accept horizontally variable initial conditions provided by GMASS or analytic solutions (e.g., solitary wave or Eady wave) using one-way open lateral bowidary conditions.The buoyancy term in the w-equation (6) has been modified (the prime term is defined as a deviation from its initial value; not fro m its horizontal average).We also ass wne that only the perturbed wind (defined as a deviation from the initial value) will satisfy the continuity equation.Both. the GCE model and GMASS have been initialized with the same analytic solutions to the Hoskins-Bretherton horizontal shear model at 60 h which is based on the semigeostrophic theory (e.g., Hoskins and Bretherton, 1972).The horizontal domain size is equal to the wavelength of the most unstable baroclinic wave and periodic boundary conditions are used.The computed wave length and phase speed of the baroclinic wave are 3507 km (15 km resolution was used in the models) and 14.35 m s-1, respectively.Both models have been tested under a dry environment (Biak et al., 1992).Good qualitative agreement between GCE model and GMASS results have been found (Figure 11 ).The differences were mainly caused by an absorbing layer which was necessary for the non hydrostatic model (in order to reduce the reflection of vertically propagating gravity waves).A much better agreement between the two models was ob tained by applying the absorbing layer to 81 and w ' only.We will examine similarities and differences in the simulated rainbands with full physics in both models.This comparison study will provide some indication as to how well the cumulus parameterizations perform in the mesoscale model.This kind of work is necessary to address the difficult proplems that will be encowitered in an eventual two-way nesting with GMASS.
The GCE model has also been initialized with a solitary wave (derived from nonlinear  Other key GCE model physics modifications are the addition of an energy balance equa tion for the ocean/land surface and incorporation of the high resolution Blackadar planetary boWldary layer (PBL), which are curr ently used in GMASS.The effects of cloud-SST (Sea Surface Temperature)-radiation interaction will be estimated by performing a series of simu lations with various SST's as well as land-ocean distributions.The inclusion of topography into the GCE model is also underway.

CONCLUSIONS AND FUTURE OUTLOOK
The GCE has evolved substantially over the past five years.Its applications to cloud en sembles in several different climatic regimes and to remote sensing of rainfall and associated latent heat release are discussed in Part II.It is clear that observations of cloud processes are greatly needed, particularly of ice-phase microphysics.A microphysical data set obtained in clouds over the Pacific warm ocean pool during TOGNCOARE will be especially valuable and will help to improve the formulations in the model.
It is recognized that bulk formulations of microphysics leave out important processes (such as particle sorting by differing fall speeds) which impact both the dynamical and radiative transfer aspects of convective systems.A close collaboration has been established with the groups of Kogan in Oklahoma U.S. and Rosenfeld in Israel who are working on explicit microphysics of both water and ice.The aim is to improve the bulk formulas by model-to-model comparison of results with detailed microphysics, since present computing capability does not permit many runs with a large number of size classes of hydrometeors.Since computing capability also limits simulation of scale interactions, improved boundary layer physics, and the allowable three-dimensional model domain, it is likely that a massive parallel processing approach will be required within the next few years.
2) Under super-(or sub-) saturated conditions condensation (evaporation) •and deposition (sublimation) occur in proportions that depend linearly on the temperature in the range 0°C and -40°C; Excess water vapor is assumed to condense/deposit Kuo Ta.o & Joanne Simpson

Table 1 .
Key to Figure 1.Meaning Depositional growth of cloud ice.Initiation of cloud ice.Melting of cloud ice to form cloud water.Depositional growth of cloud ice at the expense of cloud water.Homogeneous freezing of cloud water to form cloud ice.Accretion of rain by cloud ice; producing snow or graupel depending on the amount of rain.Accretion of cloud ice by rain; producing snow or graupel depending� the amount of rain.Autoconversion of cloud water to form rain. Accretion of cloud water by rain.Evaporation of rain.

2 )
The expressions for Dqv , Dqc• Dqr• Dqit Dq8, and Dqg are similar to De.Where a grid point is saturated, the following vertical turbulence diffusion equations are used: We i-Kuo Ta o & Joanne Simpson 8 11011 8 T.( fJBe (l + eL 2 qv )

Fig. 2 .
Fig. 2. Total rain intensity integrated over the grid points designated as the con vective and stratifonn (anvil) regions.(a) Ice run, and (b) ice-free run.

Fig. 3 .Fig. 5 .
Fig. 3. Vertical cross section of a simulated tropical squall-type convective line at its mature stage.The heavy solid lines indicate the 10, 30 and 40 dBZ radar reflectivity contours.In the dark shaded area, w > 0.5 m s-1 , and in the light shaded area, w. < --0.5 m s-1 • The convective, st ratifonn and no-surface precipitation regions are also indicated.(a) Ice run, and (b) ice-free run.
(1993a)    have found a distinct minimum in the temperature lapse rate near 4.5 km level in the GATE and Western Pacific regions.They suggested that this minimum could be caused by melting processes.

Fig. 6 .
Fig. 6.Convective heating profiles (K h-1 ) computed for the last 5 h of 2D simulations for tropical squall-type convective lines.Shown are the heat ing in the convective region (solid), heating in the anvil region (double dot-dash), and the total heating (dot-dash).(a) Graupel case, and (b) hail case.

Fig. 7 .
Fig. 7. Ve rtical cross section of the pressure gradient force term in the w-momenturn: equation over part of the model domain at (a) 60 min and (b) 480 min for a simulated mid latitude squall-type convective line using the anelastic assumption.(c) and (d) are the same as (a) and (b).respectively, except for the compressible system.
-dimensional model simulations are necessary for studying the rel ationship between vertically varying environmental winds and such features as rotation within and orientation of updrafts and downdrafts, cloud movement relative to the mean wind direction.

Fig. 8 .
Fig. 8. Heating rates for condensation, c, and evaporation, e, in the two-dimensional model (dashed line) and the three-dimensional model (solid line) .The heating rate estimated from the large-scale observations, Q 1-Q R• and the total cloud heating rate in the three-dimensional model are also included.
Ji: Table 4. Calibration of clear, active and inactive cloud drafts .b is either 1 or 2 Ratio of fractional cloud coverage (R = cloud updraft coverage I cloud downdraft coverage).Fr actional coverage occupied by cloud drafts and active cloud drafts over the domain are also shown within the parentheses.

Fig. 9 .
Fig. 9. Simulated radar reflectivities using the 41CE scheme for (a) the 29 June 1986 COHMEX squall line at 150 min of the simulation time and (b) the 12 September 1974 (Day 255) GATE squall line 360 min into the simulation.A grid increment of Ll.x=l km is used in both simulations.

Fig. 10 .
Fig. 10.Ve rtical cross section of (a) the simulated radar refiectivities, (b) the pres sure deviation and (c) the wind speed relative to the squall line movement for the PRE-STORM run using the positive definite advection scheme (840 minutes into the simulation).The contour interval is 10 dBZ , 30 Pa and 2 m s-1 for (a), (b) and (c), respectively.Figs. 10 (dHf) are the same as Fig. 10 (a)-(c), respectively, except for a fourth-order advection scheme.
analytical theory by Dr. F. Einaudi at Goddard1) as a test to see if convection could be triggered by such a mechanism as opposed to other methods.Preliminary results show that the model produced a realistic long-lasting cloud system compared with limited observations, such as echo top heights and propagation speeds derived from CHILL radar data and NMC

Table 3 .
Surface rainfal l accumulated over 4 h, 8 h and 12 h simulations for the anelastic and compressible systems, and normalized with respect to the total model domain (P 0 in mm/grid/4 h).Anvil portion is given as per centage of the total rain.

Table 6 .
Same as Table3except for the experiments with the positive definite ad vection(Smolarkiewicz)schemes and the fourth-order advection scheme.