JULY 13, 2021

principles of numerical weather prediction and modelling

By AYAPILLA MURTY
 Prof. Ayapilla Murty

CHAPTER-7

PRINCIPLES OF NUMERICAL WEATHER PREDICTION AND MODELING

1. An Overview of Numerical Weather Prediction:

The main goal of a Meteorological office is to issue weather forecasts
in different time scales to various user agencies like aviation, marine,
agriculture, water management, builders, tourism industry, planners
etc. and to general public. The forecast requirement varies from
detailed weather forecasts in time scales of a few hours to days and to
more general indication of the broad weather patterns of succeeding
months, seasons or even beyond. For example, aviation industry needs
weather information in time scales of a few hours or a day, whereas
agriculture sector demands weather forecasts in time scales of a
week. Weather forecasts of a month or a season in advance are
required mostly by planners. Thus, the weather forecasts are broadly
classified as (i) short-range (forecast validity upto 3 days), (ii) medium
range (forecast validity beyond 3 days upto 10 days) and (iii) long-
range (forecast validity beyond 10 days or a few weeks or a month or a
season or even beyond).
Weather forecasting basically consists of two steps. The first step is
to have an accurate assessment of the present/initial state of the
atmosphere. This helps in identifying the different weather systems
and their horizontal and vertical state. As we know, there are a variety
of phenomena occurring in the atmosphere having different space and
time scales. The characteristic sizes of these motions vary from a
fraction of a centimetre to several thousands of kilometres, with time
scales of a fraction of a section to several weeks. Each of the various
scales of motions has a varying degree of influence upon all the others
and it is important to properly observe, analyze and account them in
atmospheric studies and weather forecasting. As weather has no
political boundary, one needs weather data from a fairly large region;
the area from which data required increases with the duration of
forecast made, since weather systems from one part of a region may
travel and affect the weather condition over a far off region in course
of time.

1. Historical Back ground:

1
The roots of numerical weather prediction can be traced back to the work of Vilhelm Bjerknes, a
Norwegian physicist who has been called the father of modern meteorology. In 1904, he published a
paper suggesting that it would be possible to forecast the weather by solving a system of nonlinear
partial differential equations. The first attempt to predict the weather numerically was by a British
mathematician named Lewis Fry Richardson spent three years developing Bjerknes’s techniques and
procedures to solve these equations. Armed with no more than a slide rule and a table of logarithms,
and working among the World War I battlefields of France where he was a member of an ambulance
unit, Richardson computed a prediction for the change in pressure at a single point over a six-hour
period. The calculation took him six weeks, and the prediction turned out to be completely unrealistic,
but his efforts were a glimpse into the future of weather forecasting. Richardson’s book Weather
Prediction by Numerical Process was published in 1922. Richardson showed how the differential
equations governing atmospheric motions could be written approximately as a set of algebraic
difference equations for values of the tendencies of various field variables at a finite number of points in
space. Given the observed values of these field variables at these grid points the tendencies could be
calculated numerically by solving the algebraic difference equations.
The new values of the field variables could then be used to re-compute the tendencies which in
turn could be used to extrapolate further ahead in time. By extrapolating the computed tendencies
ahead a small increment in time, an estimate of the fields at a short time in the future could be
obtained. Even for a short-range forecast over a small area of the earth this procedure requires an
enormous number of arithmetic calculations. Richardson did not foresee the development of high speed
computers at that time. Estimated that a work force of 64 000 people would be required simply to keep
up with the weather on a global basis.
Richardson tried a “forecast” by hand –unfortunately the results were poor due to poor initial
data and the inclusion of fast waves like sound and gravity waves. Numerical weather prediction was not
attempted again until after World War 2 –interest grew due to improved meteorological observational
network and the development of digital computers. J. G. Charney showed in 1948 how the dynamical
equations could be simplified using the geostrophic and hydrostatic approximations so that sound and
gravity waves were filtered out (essentially the quasi-geostrophic model). A special case of his model the
Equivalent barotropic model was used in 1950 to make the first numerical forecast.
The model provided forecasts of the geopotential near 500 mb. Thus it did not predict weather
in the usual sense. However, it could be used to assist forecasters in predicting local weather as a result
of large-scale circulations. Later versions of the quasi-geostrophic model provided surface pressure and
temperature predictions, however, the accuracy of the prediction was limited by the quasi-geostrophic
assumptions. With the development of more powerful computers and better modeling techniques
numerical weather prediction has returned to models that are quite similar to Richardson’s model and
are more accurate. Current operational models are generally hydrostatic. The exception to this are the
research prototype forecast models such as MM5 and RAMS.
2. Continued Evolution of Numerical Weather Prediction:
The evolution of numerical weather prediction throughout the latter part of the 20th century
proceeded at a similar pace at many operational numerical weather prediction centers around the
globe. The restriction on NWP was primarily due to the amount of computer power as well as the
amount of data available to initialize the model. However, a progression of more and more powerful
computers in the 1960s and 1970s as well as increasing sources of data—particularly from weather
satellites—allowed the expansion of both the domains and the number of models run. Increases were
also made in the number of vertical levels and the horizontal resolution of the models. A three-layer
hemispheric model was introduced in 1962 and a six-layer primitive equation model appeared in 1966.

2
Additional atmospheric layers allowed more accurate forecasts of winds and temperature, resulting in
better prediction of storm motion.
The late 1970s and the 1980s saw the introduction of the use of spectral coordinates and
advances in methods of data assimilation and more accurate representation of physical processes in the
atmosphere such as the formation of clouds and precipitation. During the 1990s, models of finer and
finer scales were developed.
A more recent addition to the numerical weather prediction suite of products is “ensemble
modeling.” With ensemble modeling, many forecasts are run with slightly varying initial conditions from
different models and an average, or "ensemble mean," of the different forecasts is created. This
ensemble mean can extend forecast skill to the two-week range because it averages over the many
possible initial states, to essentially “smooth” the chaotic nature of the atmosphere. In addition, this
technique provides information about the level of uncertainty for different conditions because of the
large ensemble of forecasts available.
As soon as new computer resources became available, they were fully utilized with higher-
resolution models with improved representation of physical processes, using more and more
sophisticated and complex modeling techniques and an increasingly larger number of observations from
around the globe, especially from satellites. Advances in predictive skill could be shown to be directly
proportional to advances in computing and modeling.
However, the computer-based forecast is not the only reason for the improvements in
forecasting accuracy the world enjoys today. The interpretation of the model forecasts by human
forecasters has been shown to provide a consistent amount of additional forecast improvement
throughout the evolution of numerical weather prediction. Our ability to interpret the many “moods” of
the atmosphere has not yet been matched by equations and computer chips.
Though computing power increased in the 1960s, Edward Lorenz found that these equations are
extremely sensitive to initial conditions and that very small perturbations to the initial conditions (or
boundary conditions or other forcings) could cause large deviations in the computed flow conditions at a
later time (Lorenz 1963). This sensitivity to initial conditions has become known as “chaos” and has
altered meteorologists’ views of predictability. The sensitivity is caused by an inability to accurately
represent the current atmospheric state. Thus, resulting errors grow and saturate the energy spectrum.
This inherent sensitivity means that because we cannot precisely know the conditions at all places at the
initial time, we are doomed to a loss of predictability at future times. Because of this limitation
meteorologists have embraced probabilistic forecasting to quantify uncertainties in the forecast.

3. Prediction across scales

It is important to understand the wide range of scales in weather forecasting and their impact on
providing an accurate forecast. One can think of the atmospheric motion as being composed of a large
number of nonlinearly interacting waves. The largest scale of wave is those weather systems (lows
followed by highs) that circle the globe in the mid-latitudes (Fig. 1), the Rossby waves that form the
interesting daily weather in mid-latitudes. These weather systems have wavelengths on the order of
about ten thousand kilometers and their influence persists over a given region for a few days. Other
circulations develop at smaller scales. The mesoscale is responsible for diurnal circulations usually
originated by heterogeneities at the surface. The sea breeze associated with a differential heating of the
land and the ocean, like monsoon, is a clear example. Other types of waves exist are gravity waves (like
those waves observed on the surface of a water body), inertial waves that travel through the interior of
the rotating atmosphere, and sound waves due to the compressibility of the air. Some of the smallest

3
eddies in the boundary layer have wavelengths on the order of centimeters. This implies that the range
of scales of atmospheric motion is on the order of eight decades. Given current computing power, one
cannot resolve all these scales with a numerical model. Energy from variations in the smallest scales also
can backscatter, or propagate to larger scales and can cause a numerical integration to “blow up” due to
numerical instability.

Figure 1. Depiction of the general circulation of the atmosphere and associated Rossby waves.

At the very short range, on the order of seconds or minutes to a couple of hours, it is best to
have situational awareness through measurements at the site, and also at short distances upwind from
the site. It is difficult to improve upon persistence forecasts at these temporal scales; that is, a forecast
would assume that the current value of wind speed or irradiance persists to the next time period.
Methods that improve upon persistence tend to be statistical learning methods (Sharma et al. 2011).
Beyond an hour or two, NWP becomes important. This is where we leverage the equations of
motion and numerically integrate them forward in time. It is important to initialize those equations with
the best estimate of the atmospheric state (also known as an analysis) via data assimilation. Integrating
the equations forward from the current analysis state provides an estimate at the future time. NWP
continues to provide value until about 2 weeks, when predictability approaches climatology.
Beyond a couple of weeks, the specified integrations are washed out, but one can still use
general trends associated with long-range atmospheric oscillations, such as the ones associated with El
Niño/Southern Oscillation, which is associated with heating in the tropical Pacific that produces standing
waves in the atmosphere that may persist for a period of time, affecting global weather patterns for
months at a time. Understanding and predicting such oscillations leads to seasonal predictions that may

4
be better than climatology. On the very longest (climatic) scales, there is a potential for long-term trends
due to changes in atmospheric composition.
4. Data assimilation (DA):
As mentioned earlier, the atmosphere is chaotic. Therefore, a major key to reduce errors in NWP
forecasts is to specify the initial state of the atmosphere as accurately as possible. Some form of data
assimilation (DA) is used for this task in NWP. Data assimilation is the “process in which observations
distributed with time are merged together with a dynamical numerical model of the flow in order to
determine the state of the atmosphere as accurately as possible.” There are a few basic classes of DA
algorithms in use today, all of which use observations from the numerous platforms mentioned above to
increment from a background (prior) state to an analysis (posterior) state, and a few of these are
highlighted below. The analysis state serves as the initial condition for a subsequent forecast.
4.1 Nudging

Observation nudging, or Newtonian relaxation, is an empirical analysis scheme in which the

NWP model state is “nudged” toward an observation by means of adding a small extra term to the
model prognostic equations. In modern four-dimensional data assimilation (FDDA) nudging schemes,
this nudging term is weighted in both time and space around the observation to avoid introducing
artificial gravity waves. FDDA is a continuous data assimilation technique, as it assimilates observations
sequentially through time, rather than treating all observations within a given time window as being
valid at a single analysis time. The limitation of FDDA is that it can only assimilate observations that can
be converted to model prognostic variables.

4.2 Variational assimilation

Variational assimilation techniques are a form of statistical interpolation. All statistical

interpolation techniques require an estimation of the error covariances between variables in the
background state, as well as error covariances between the observed variables. These techniques find
the optimal analysis by globally minimizing a cost function that incorporates the distance between the
analysis and observations within the assimilation window. This method also requires the observation
error covariance matrix and the background error covariance matrix which at times may be difficult to
compute accurately. In three-dimensional variational Data assimilation (3D-Var) schemes, these error
covariance matrices come from a static climatology, and all observations within a given assimilation
window are assumed to be valid at the analysis time. These assumptions reduce the computational
burden. In contrast, four-dimensional variational Data assimilation (4D-Var) schemes seek to minimize
the cost function, subject to the NWP model equations, to find the best model trajectory through the
entire assimilation window, rather than just at the analysis time. In addition, the error covariance
matrices are flow-dependent in 4D-Var. These differences make 4D-Var significantly more
computationally intensive than 3D-Var, but also more accurate.

4.3 Coupled models

In addition to Data Assimilation systems, another way to further reduce forecast errors from NWP
models is to couple them to models that simulate other parts of the Earth system. For instance, it is

5
known that both the wind and the waves mutually impact each other in the marine atmospheric
boundary layer. Efforts to build two-way coupled atmosphere-wave models are crucial to improve
predictions for a number of applications, including offshore wind energy. Additionally, coupled
atmosphere-hydrology models reduce forecast errors for several applications, including hydro energy
and even wind energy, thanks to the impact of soil moisture on boundary layer winds.

5. Fundamentals of NWP

The rational basis for NWP was proposed at the beginning of the 20th century. Bjerknes (1904)
suggested that the weather could be predicted using the fundamental conservation laws of nature.
These equations, when expressed in mathematical form, form a system of equations with the same
number of unknowns and equations. In principle a solution can be found, allowing a forecast of the
evolution of the atmosphere. He anticipated that in order to solve the system of atmospheric equations
of motion, one would first need to develop an analysis of the atmospheric state (i.e. temperature,
humidity, winds, etc.) and then use this information for prognostic analysis of subsequent states. The
first step of the analysis is facilitated by data assimilation. The second step, the prognosis or forecast, is
challenging because the analytical solution of the atmospheric equations of motion is unknown. Hence,
we need to numerically solve the atmospheric equations of motion. The method used to solve and get
the solution is the NWP model.
The number of calculations necessary to approximate a solution is enormous and NWP was not
feasible until the advent of digital computers. In a visionary work, Richardson (1922) approximated a
solution to the atmospheric equations of motion using the basic ideas that we still use in modern NWP
models. He discretized the model equations onto a grid, beginning with observed conditions (which at
that time required wiring weather conditions from upstream locations), and envisioned that a team of
human “calculators” could perform the numerical calculations. Although his method failed (he did not
fully understand the impact of the small scales on the larger scales, and the unfiltered equations “blew
up”), thus was born the concept of numerical weather prediction. Based on this work, the development
of digital computers, and with improved understanding of atmospheric dynamics, Charney et al. (1950)
made the first successful numerical weather prediction. A rapid increase in NWP model developments
took place after this pioneering work, leading to modern NWP models (Lynch 2006).
Current NWP models have two main components: a dynamic solver and a set of physical
parameterization schemes. These components are described below.

5.1.Dynamic solver

The dynamic solver uses algebraic approximations to represent the differential terms in the atmospheric
equations of motion. This requires the selection of appropriate spatial and temporal discretization.
Typical terms approximated by the dynamic solver include advection, the pressure gradient force, and
the Coriolis force resulting from Earth’s rotation.
Different approaches have been proposed to represent the differential terms. The most
common approach is to use finite difference schemes. These methods expand the spatial and temporal
differential terms into a polynomial series (a truncated Taylor series). This type of dynamic solver is

6
simple, relatively easy to code, and widely used in limited area models. An alternative approach is to use
spectral methods. The fundamental idea of spectral methods consists of expanding atmospheric
variables (temperature, winds, moisture, etc.) in frequency (or wavenumber) space. Spectral methods
provide the highest possible order of accuracy for a discretized set of partial differential equations, and
therefore result in reduced errors in representing the transfer of energy towards the finest atmospheric
scales. Spectral methods are most often used in global NWP models and in fine-scale large eddy
simulations (LES).

5.2. Parameterizations

Parameterization is the process of Inclusion of the effects of physical processes implicitly

(indirectly) when we cannot include them explicitly (directly). In other words each important physical
process that cannot be directly predicted requires a parameterization scheme. Thus it is necessary due
to meagre computation facilities and some of the physical processes cannot be modelled or expressed in
numerical means that computer can do.
Regardless of the use of dynamic solver, there are additional processes that need to be
represented in an NWP model. These include physical processes that we cannot resolve at the selected
grid spacing, processes that we do not know how to represent analytically and processes that require
large amounts of calculations to represent explicitly (directly). In these situations, the model includes
the effects of these terms in the form of physical parameterization schemes. The model physics,
therefore, can introduce large approximations than those introduced in the dynamic solver. All NWP
models parameterize a number of physical processes. These include atmospheric radiation, land-surface
interactions, turbulent mixing, convective clouds, and cloud microphysics.
Radiation parameterizations provide atmospheric heating rates and downward irradiance at the
surface. The radiative transfer calculations are divided into shortwave and longwave components.
Longwave radiation accounts for the infrared radiation absorbed and emitted by atmospheric gases,
clouds, and the earth’s surface. On the other hand, shortwave radiation accounts for the absorption,
reflection and scattering by gases, clouds and the earth's surface of the solar spectrum centered in the
visible region. Both shortwave and long wave radiative transfers are affected by the modeled water
vapor, clouds, and gases such as ozone, oxygen, and carbon dioxide. Some advanced radiation schemes
also model the interaction with various species of modeled atmospheric aerosols.
PBL parameterizations account for the vertical mixing associated with atmospheric turbulence in
the PBL, which grows and decays with diurnal heating and cooling of the underlying surface. The
turbulent mixing usually includes temperature, water vapor, horizontal momentum, and trace gases.
The atmospheric turbulent fluxes are usually modeled as a combination of local and non-local mixing.
The parameterizations require a quantification of the surface turbulent fluxes, which are usually
computed via Monin-Obukhov (1954) similarity theory in a surface layer parameterization, which
models these fluxes below the lowest model layer of the PBL scheme, and in the Land Surface Model
directly.
Cumulus parameterizations account for the effects on temperature and moisture profiles
associated with unresolved deep convective and shallow convective clouds. In addition, the schemes
provide surface precipitation that feeds back into the LSM. The deep convective schemes are usually

7
based on closure assumptions valid for horizontal grid spacing higher than 10 km. These approximations
are less accurate at finer grid spacing because the NWP starts to resolve the convective
updrafts/downdrafts that are parameterized. The deep convective scheme is usually turned off at a grid
spacing finer than about 5 km. However, shallow convection may still need to be parameterized at a grid
spacing of about 1 km.
Finally, microphysics parameterizations represent processes involving cloud droplets and ice
crystals and the effects of their life cycle on the temperature and moisture profiles, precipitation, and
cloud radiative properties. Microphysics can account for processes in liquid water clouds (warm clouds)
or clouds consisting in water and ice hydrometeors (mixed phase clouds). Hence, the parameterizations
differ in the number of hydrometeors that are predicted and advected by the model dynamics.
Parameterizations in microphysics are usually based on a bulk approach wherein the size distribution of
hydrometeors follows a functional form, and one or more parameters of the distribution are predicted
by the parameterization (e.g. mixing ratios).

6. Methods of Modelling:

Since the equations governing the meteorological variables change with time, if we know the initial
condition of the atmosphere, we can solve the equations to obtain new values of those variables at a
later time (i.e., make a forecast). The model forecast equations are simplified versions of the actual
physical laws governing atmospheric processes. Due to their complexity, the primitive equations must
be solved numerically using algebraic approximations, rather than calculating complete analytic
solutions as computers can only estimate the physics if expressed in terms of mathematics .
For example, to represent an NWP model in its simplest form, we can write:

Where ΔA gives the change in the forecast variable at a particular point in space, Δt gives the change in
time (how far into the future we are forecasting) and F(A) represents the terms that can cause changes
in the value of A . This equation means that the change in forecast variable ‘A’ during the time period ‘t’
is equal to the cumulative effects of all processes that force ‘A’ to change.

Future values of meteorological variables are solved by finding their initial values and then adding the
physical forcing that acts on the variables over the time period of the forecast. This is stated as:

8
where F(A) stands for the combination of all kinds of forcings that can occur. The process is
then repeated for the duration of the simulation. This stepwise process represents the
configuration of the prediction equations used in NWP.
The actual equations used in ECMWF model are:

6.1.Expressing the equations in to finite difference form:

As computers can solve only mathematics we have to express the equations of motions into finite
difference form and neglecting higher order (nonlinear) terms called truncation.

9
Similarly For example, the moisture equation can be expressed in finite difference form as:

10
Fig.2.For writing the past, present and future quantities

6.2. Types of Models:

There are four different types of models: .Grid point models, Spectral models, Hydrostatic models, and
Non-hydrostatic models.

6.2.1. Grid point Models:

11
In the real atmosphere, wind, pressure, temperature, and moisture vary from location to location in a
smooth and continuous way. Grid point models, however, perform their calculations on a fixed array of
spatially disconnected grid points. The values at the grid points actually represent an area average over
a grid box. For example, the continuous temperature field shown in the Fig.3, must be represented at
each grid point as shown by the black numbers in the right panel of the graphic. The temperature value
at the grid point represents the grid box volume average.

Fig.3.Grid point model (right panel) representation of temperature of the isotherms in the left panel

The Grid point models actually represent the atmosphere in three-dimensional grid cubes, as shown in
Fig.4 below. The temperature, pressure, and moisture (T, p, and q), shown in the center of the cube,
represent the average conditions throughout the cube. The east-west winds (u) and the north-south
winds (v), located at the sides of the cube, represent the average of the wind components between the
center of this cube and the center of the adjacent cubes. Similarly, the vertical motion (w) is represented
on the upper and lower faces of the cube. This arrangement of variables within and around the grid
cube (called a staggered grid) has advantages when calculating derivatives. It is also physically intuitive;
average thermodynamic properties inside the grid cube are represented at the center, while the winds
on the faces are associated with fluxes into and out of the cube. Grid points are generally evenly spaced
(or nearly so) except latitude –longitude grids and adaptive grids. The grid point models use finite
difference techniques to solve the forecast equations. Deficiencies in the ability of the finite difference
approximations to calculate exactly the gradients and higher order derivatives are called truncation
errors. If the grid spacing is large we may miss the short period variations and increase the magnitude of
errors introduced into the solutions of the model forecast equations. For example, as shown in Fig.5, in
the forecast of an urban smog of 40 km distance with only two grid points, the peaks are missing and
nearly 25% of the information is lost. So the distance between the grid points is selected in such a way
that there are sufficient number of grid points to adequately represent the smallest feature of interest.

12
Fig.4.Staggered grid showing winds and other parameters in three dimensional form

Fig.5. Two grid points missed the peak value forecast of Urban Smog of 40 km radius

6.2.2. Spectral Model:

Spectral Methods Involves the forward transformation of the standard dependent variables (u, v, T, etc.)
into transform space utilizing Fourier transforms. The resultant Fourier series and Fourier-Legendre
functions represent horizontal variability. Temporal and vertical derivatives are handled with
conventional methods. The dependent variables are subsequently transformed back into physical space
to get interpretable forecast fields.
Spectral models represent the spatial variations of meteorological variables as a finite series of waves of
differing wavelengths. For example, consider the hemispheric 500-hPa height field in the top portion of
the Fig.6. If the height data are tabulated at 40°N latitude every 10 0 of longitude (represented at each
yellow dot on the chart), there are 36 points around the globe. It takes a minimum of 5 to 7 points to
reasonably represent a wave and, in this case, five or six waves can be defined with this data. The
locations of the wave troughs are shown in the top part as solid red lines.

13
 Fig.6. Spectral wave representation of 500-hPa height field in N.H.

When the data are plotted on the graph (Lower panel of Fig.6), the five wave troughs are definable by
the blue dots but are unequally spaced. This indicates the presence of more than one wavelength of
small-scale variations. In this case, the shorter waves represent the synoptic-scale features, while the
longer waves represent planetary features. Thus by spectral method we are able to get the finer details
of physics.
Recall that in a grid point model, truncation error is associated with the finite difference
approximations used to evaluate the derivatives of the model forecast equations. One of the nice
features of the spectral formulation is that most horizontal derivatives are calculated directly from the
wave and so relatively more accurate. However, the degree of truncation error in a given spectral model
depends on the scale of the smallest wave represented by the model.
Though several types of wave shapes are possible in spectral models, triangular (T, as in T170)
configuration is the most common in operational models since it has roughly the same resolution in the

14
zonal and meridional directions around the globe. However, if the number of waves in the model is
small (for example, T80), only larger wave features (A 1 & A2 in Fig.7)) can be represented and smaller
scale wave features (A3 with red Zig-zagged line in Fig.7) observed in the atmosphere will be entirely
eliminated from the forecast model as shown in Fig.7.

Fig.7. Smaller scale waves (A3) in spectral method cannot be represented and so ignored

6.2.3. Hydrostatic models:

Most grid point models and all spectral models in the current operational NWP models are hydrostatic.
The hydrostatic assumption is taken in this modelling. This means that no vertical accelerations are
calculated explicitly. So these models cannot predict details of small-scale processes associated with
buoyancy. The hydrostatic assumption is valid for synoptic-and planetary-scale systems and for some
mesoscale phenomena. An important exception is deep convection, where buoyancy becomes an
important force. Hydrostatic models account for the effects of convection using statistical
parameterizations approximating the larger-scale changes in temperature and moisture caused by non-
hydrostatic processes. These models are used for forecasting synoptic-scale phenomena, can forecast
some mesoscale phenomena.

6.2.4. Non-hydrostatic models:

Currently, most non-hydrostatic models are grid point models. They are generally used in forecast or
research problems requiring very high horizontal resolution (from tens of meters to a few kilometers)
and cover relatively small domains. Non-hydrostatic models can explicitly (directly) forecast the release

15
of buoyancy in the atmosphere and its effects on the development of deep convection. To do this, non-
hydrostatic models must include an additional forecast equation that accounts for vertical accelerations
and vertical motions directly, rather than determining the vertical motion diagnostically from horizontal
divergence. The equation is used is as given in the box.

6.3.Vertical Coordinates:

To represent the vertical structure of the atmosphere properly requires selection of a suitable vertical
coordinate and sufficient vertical resolution. Unlike the horizontal structure of models where discrete or
continuous (grid point or spectral) configurations can be used, virtually all operational models use
discrete vertical structures. As such, they produce forecasts for the average over an atmospheric layer
between the vertical-coordinate surfaces, not on the surfaces themselves.
The different vertical coordinate systems used are by height, pressure, Sigma, potential
temperature, Eta and hybrid. Height and pressure normally cannot be directly used as vertical
coordinate as it is difficult to represent the topography correctly as shown in Fig.8A

Fig.8A. Schematic illustration of a rectangular, (b) isobaric, (c) isentropic and (d) sigma coordinate
representation

16
6.3.1. Sigma-vertical coordinate system:

The equations of motion have their simplest form in pressure coordinates. Unfortunately, pressure
coordinate systems are not particularly suited to solving the forecast equations because, like height
surfaces, they can intersect mountains and consequently 'disappear' over parts of the forecast domain.
To deal with this problem Phillips (1957) developed a terrain-following coordinate called the sigma (σ)
coordinate. The sigma coordinate or variants are used in the NGM, GFS, ECMWF, NOGAPS, and UKMET
models and appear in some mesoscale models, such as MM5, COAMPS, and RAMS.

The sigma coordinate is defined by σ = p/ps

where p is the pressure at the forecast level within the model and ‘ps’ is the pressure at the earth's
surface, not mean sea level pressure. As shown in Fig. 8B, the lowest coordinate surface (usually labeled
σ = 1) follows a smoothed version of the actual terrain. Note that the terrain slopes used in sigma
models are always smoothed to some degree. The other sigma surfaces gradually change from nearly
parallel to the smoothed terrain at the bottom of the model (σ = 1) to being nearly horizontal to the
constant pressure surface at the top of the model (σ = 0). The sigma vertical coordinate can also be
formulated with respect to height (z), rather than pressure. The advantage of sigma coordinate system
is easy to represent top and bottom of the atmosphere and the disadvantage is equations need to be
transformed as errors occur in horizontal pressure gradient force when terrain slope is steep.

Fig.8B. Representation of Sigma coordinate on a terrain

17
Schematic illustration of a rectangular, (b) isobaric, (c) isentropic and (d) sigma coordinate
representation

6.3.2. Eta Coordinate system:

The Eta (η) coordinate was created in the early 1980s in an effort to reduce the errors incurred in
calculating the pressure gradient force using sigma coordinate models. The Eta coordinate is, in fact,
another form of the sigma coordinate, but uses mean sea level pressure instead of surface pressure as a
bottom reference level. As such, Eta is defined as:
Pr ( Z S )−P t
η S=
P r ( Z=0 )−Pt
where pt is pressure at the model top; pr(z=0) is the standard atmosphere MSL pressure (1013 hPa) and
pr(zs) is the standard atmospheric pressure at the model terrain level ZS..
The Fig.9 shows the representation of Eta (η) at different pressure levels as per the above equation.

6.3.2.1. Calculation of ‘η’ surfaces:

Let suppose the terrain is as shown in Fig. 9 and distributed evenly with respect to pressure from sea
level to the top of the atmosphere. The standard atmospheric pressures are then determined at each of
the heights at points 1,2, 3 etc. The calculation of Eta (η) level at point 1 is as follows: the actual terrain
elevation is 848m which is close to the standard atmospheric pressure at that height is 900 hPa. The sea
level pressure is closer to the 1000-m height. Using the equation of Eta (η) can be calculated as:

Pr ( Z S )−P t 900−0
η S= = =0. 9
P r ( Z=0 )−Pt 1000−0

Fig.9. Calculation of Eta (η) levels at different pressure levels on an irregular terrain

Similarly if we go to point 2, the actual terrain height is 1126 m and is also closer to 1000-m height.
Therefore, the eta level closer to this point is again 0.9. But, if we go to point 3 (1832 m), the nearest eta

18
surface in the model is at 2000 m. Here the standard atmospheric pressure is 800 hPa. The nearest eta
(η) level is therefore = 1 x (800-0 / 1000-0) = 0.8.
Note that the eta (η) levels are predefined and the model topography is set to the nearest eta (η)
surface even if it does not quite match the average or smoothed terrain height in the grid box. This has
been a simplified example. In reality, it is necessary to choose the intervals between eta (η) levels in a
way that both depicts the planetary boundary layer (PBL) with sufficient detail and yet represents the
average changes in elevation over the entire forecast domain.
Eta is usually labeled from 0 to 1 from the top of the model domain to mean sea level. Some of
the model's grid cubes are located underground in areas where the surface elevation is notably above
sea level. This requires special numerical formulations to model flow near the earth's surface. The Eta
coordinate systems allows the bottom atmospheric layer of the model to be represented within each
grid box as a flat "step," rather than sloping like sigma in steep terrain. This configuration eliminates
nearly all errors in the Pressure Gradient Force (PGF) calculation and allows models using the eta
coordinate to have extreme differences in elevation from one grid point to its neighbor.
Eta coordinate models can therefore develop strong vertical motions in areas of steep terrain
and thus more accurately represent many of the blocking effects that mountains can have on stable air
masses. Even when the step-like Eta is used as the vertical coordinate, model terrain is still much coarser
than real terrain, but the topographic gradients are less smoothed than in sigma models. Note that the
Eta levels are predefined and the model topography is set to the nearest eta surface even if it does not
quite match the average or smoothed terrain height in the grid box.

6.4.Sigma verses Eta coordinate system:

Eta models do not need to perform the vertical interpolations that are necessary to calculate the PGF in
sigma models. Although the numerical formulation near the surface is more complex, the low-level
convergence in areas of steep terrain are more representative of real atmospheric conditions than in the
simpler formulations in sigma models. Compared with sigma models, Eta models can often improve
forecasts of cold air outbreaks, damming events, and leeside cyclogenesis. For example, in cold-air
damming events, the inversion in the real atmosphere above the cold air mass are preserved almost
exactly in an eta model. As a result, there is little or no contribution to erroneous horizontal
temperature gradients in the lee of the mountains

19
Fig. 10 Sigma versus Eta model representation

6.5.Isentropic Vertical Coordinate system:

Since the flow in the free atmosphere is predominantly isentropic, potential temperature (θ) can be very
useful as a vertical coordinate system. However, non-adiabatic processes dominate in the boundary
layer and isentropic surfaces intersect the earth's surface. For these reasons, potential temperature
alone is not currently used as a vertical coordinate in any operational numerical model system.
However, isentropic coordinates do form an essential part of many hybrid vertical coordinate systems.
The theta coordinate allows for more vertical resolution in the vicinity of baroclinic regions, such
as fronts, and near the tropopause. For adiabatic motion, air flows along constant theta (isentropic)
surfaces and implicitly includes both horizontal and vertical displacement. Vertical motion through
isentropic surfaces is caused almost exclusively by diabatic heating. The vertical component of the
isentropic forecast equations is related entirely to diabatic processes. Isentropic coordinate models
conserve important dynamical quantities, such as potential vorticity.

20
 Fig.11. Isentropic vertical coordinate system

6.5.1. Hybrid Vertical Coordinates:

Different hybrid combinations are currently in use, eg: hybrid isentropic-sigma vertical
coordinates, hybrid sigma-pressure coordinate, hybrid isentropic-sigma coordinate
Hybrid isentropic-sigma coordinate models have a combination of sigma layers at the bottom
that shift to isentropic layers above. The union of theta and sigma into one vertical coordinate system
combines the terrain-following advantages of sigma and the increased vertical resolution in key
baroclinic areas due to the adaptive nature of isentropic surfaces.
The table 1 gives a summary of advantages and limitations of different vertical coordinate
systems.

21
Fig.12. Hybrid (Isentropic-Sigma) vertical coordinate system

Table1: Summary of vertical coordinate system

22
6.6. Horizontal resolution:

Horizontal resolution is defined as the distance between grid points for grid point models and the
number of waves that can be resolved for spectral models as shown in Fig.13.

Fig.13. Horizontal resolution in case of grid point and spectral models

6.6.1. Spectral resolution and T number:

In spectral models, the horizontal resolution is designated by a "T" number (for example, T80), which
indicates the number of waves used to represent the data. The "T" stands for triangular truncation,
which indicates the particular set of waves used by a spectral model. Spectral models represent data
precisely out to a maximum number of waves, but omit all, more detailed information contained in
smaller waves. The wavelength of the smallest wave in a spectral model is represented as:

3600
minimum wavelength = N

where N is the total number of waves (the "T" number).

Terrain representation in numerical models is usually much smoother than in reality, especially
in coarse resolution models. Terrain smoothing is largely a function of a model's horizontal resolution
and the detail of the topographic dataset. Terrain smoothing can be a large source of error in regions
affected by small-scale orographic features. Several methods are used to convert terrain height data to
model orography in order to get a representative terrain value for each model grid box.
Three different types of orography is in vogue as shown in Fig.14 . Envelope orography is like a
blanket over the terrain that clings to plateaus but rises smoothly to cover all but the very sharpest
peaks. This ensures that the full effect of terrain blocking and lee side cyclogenesis can be simulated.
However, the valleys are filled and the terrain is considerably smoothed. Mean orography uses the

23
average of the terrain data inside a model grid box. This trims the tops of mountains considerably,
especially if the grid box is large, greatly diminishing the blocking effect on cross-mountain flow.
Silhouette orography averages only the tallest features in each grid box. Thus, it lies below the
mountaintops but above the valleys. The NCEP Eta Model uses a variation of silhouette orography that
enables more detail of valleys.

Fig.14. Three different types of Orographic averaging used in NWP

7. Domain and Boundary conditions:

Model domain refers to a model's area of coverage (Fig.15). Limited-area models (LAMs) have
horizontal (lateral) and top and bottom (vertical) boundaries. Global models cover the entire earth
and have only vertical boundaries. For limited-area models, larger-domain models supply the data
for the lateral boundary conditions.

24
 Fig.15. Model domain and associated boundary conditions

For all models, accurate information must be provided for all forecast variables and along each model
boundary (lateral, top, and bottom) in order to solve the forecast equations. Boundary values can be
obtained from a variety of sources, including: Data assimilation systems, Forecast values from a current
or previous cycle of a large-scale model (as is the case for lateral boundary conditions used in LAMs).
Some type of climatological or fixed value (for specifying certain surface characteristics, such as soil
moisture, sea surface temperature, and vegetation type).

8. Forecasting types:
Usually a weather forecast is for a given region and has a definite period of validity. According to the
period of validity one can classify as:
– Nowcasting : 0-6hrs , – Short range : up to 3 days , – Medium range : 3-10 days, – Long/Extended
range : more than 10 days ,Monthly and seasonal means ( -- Climate ): 30yrs or more

There are three weather forecasting methods. They are i) synoptic weather forecasting ii) Numerical
methods and iii) statistical methods.
i) Synoptic weather forecasting:
Synoptic means the observation of different weather elements refers to a specific time of observation.
From the careful study of weather charts over many years, certain empirical rules can be formulated.
These rules will help to estimate the rate and direction of movement of weather systems like a trough
line, hurricane or thunderstorm.
In the synoptic method, a forecaster attempts to predict the future changes in the state of atmosphere
from its initial state using his theoretical knowledge and experience. Here, various weather charts,
commonly known as synoptic charts are analyzed at a fixed time to understand the three dimensions of
the atmospheric state. This method can generate forecast for a broad period over a broad region, can’t
generate time and location specific objective forecast. That way, this method is subjective and the
success of the forecast depends heavily on the knowledge, experience and skill of the forecaster. The
synoptic method is widely used in short range weather forecasting, especially in tropics.
ii) Numerical methods:
Numerical method is based on the fact that gases of the atmosphere follow a number of
physical principles. If the current conditions are known, these physical laws help to forecast the future
condition.
A series of mathematical equations are used to develop theoretical models of the general
circulation of the atmosphere. These equations are used to specify changes in the atmosphere as time
changes using the weather elements like wind, temperature, humidity, clouds, rain, snow, evaporation
etc.
The atmosphere is divided into 6 or 11 layers for this purpose. As the initial conditions (initial
state of the atmosphere) are known through observations using surface and upper air data, it is possible
to predict using numerical models for future conditions.
Valuable contribution in the development of numerical methods was made by Charney (1950)
and Obukov (1954).
As there is large number of weather variables and only a less number of equations are available,
some weather variables are ignored on the assumption that they do not change with time. As large
numbers of calculations are involved even with simplified models, a high speed super computer is
necessary.
For carrying out numerical models the following considerations are followed:

25
The equations are simplified. The equations are designed in such a way that they guarantee the
conservation of air mass, momentum, water vapor and total energy for all time. The region to be
forecasted is divided into number of nodal points in a rectangular shape. The equations are solved at
each nodal point for a short period of 10 minutes. By repetitive calculations for every next 10 minutes
forecast is obtained for 24, 48 or 72 hours.
L.F.Richardson in 1922 attempted to solve the differential equations of momentum, energy, state and
conservation using finite difference method with the help of a desk calculator.
After World War II and the development of digital computer, it was realized that Richardson’s
scheme was not simple as the governing equations contains not only slow moving meteorologically
important motions but also high speed sound and gravity waves. Though these waves are very weak,
they amplify surprisingly in numerical methods and thus generate ‘noise’.
So J.G.Charney simplified the equations by introducing the scale analysis, geostrophic and
hydrostatic assumptions. By these techniques sound and gravity waves are filtered. Then these
equations are called quasi-geostrophic models. A special case of this model, called equivalent barotropic
model, was used in 1950 which gave the first successful numerical forecast. The models and the forecast
process is given in Fig.16.

Fig.16. Models and the forecast process

A. Objective Analysis
The automatic process of transforming of irregularly distributed observed data from different
geographic locations at various times to numerical values at regularly spaced grid points at fixed times is
called objective analysis. In the earlier days of numerical weather prediction, the two-dimensional fitting
of polynomial was applied to observational data in an area surrounding a grid point at which the value

26
of analysis is required. Polynomial fitting is suitable for processing relatively dense and redundant data,
but it often gives unreasonable results in data of sparse areas.
More recent procedures of objective analysis practiced at most of the operational forecasting
centers consist of the steps schematically described in Figure 17 below. Most of the observed data are
collected worldwide in every 6 hours interval, say 00, 06, 12 etc. UTC (Coordinated Universal Time), and
are blended into forecast values appropriate to the map time by the process referred to as data
assimilation. The basic process of data assimilation, indicated by the box enclosed by dashed lines in
Figure 2 below, consists of two steps: One is objective analysis and the other is called initialization. The
step of objective analysis is carried out by statistical interpolation by taking into account all of the
available observations plus other prior information. Short-term forecasts from the previous analysis
cycle are used as prior information and are referred to as background fields.

B. Concept of Initialization
The process of providing initial value data to a model is known as initialization. The quality of a
numerical forecast strongly depends on the quality of the initial conditions.
The solutions of primitive equation models correspond to two distinct types of motion. One type
has low frequency. Its motion is quasi-geostrophic and meteorologically dominant. The other
corresponds to high-frequency gravity-inertia modes. The amplitude of the latter type of motion is small
in the atmosphere. Hence, it is important to ensure that the amplitudes of high-frequency motions are
small initially and remain small during the time integration when the primitive equations are solved as
an initial value problem. The process of adjusting the input data to ensure this dynamical balance is
called initialization.
Solution to the initialization problem was central to the successful transition in forecast practice
from the use of quasi-geostrophic models to primitive equation systems during the 1970s. Since gravity-
inertial motions are filtered out in quasi-geostrophic models, no special procedure was necessary and
the objectively analyzed data could be used immediately as the input data to quasi-geostrophic models.
Actually, there is an interesting history that led to solution of the initialization problem that is equivalent
to the quest of understanding what the primitive equation forecast system really is. A promising break
came with the advent of so-called nonlinear normal mode initialization (NNMI) during the latter part of
the 1970s.

27
Fig.17.Schematic diagram for atmospheric forecast-analysis system. The two processes in the square enclosed by
dashed lines can now be combined in the four-dimensional variational data assimilation using a continuous stream
of observed data.

There are two general types of initializations: a) Static initialization (cold start) and b) Dynamic
initialization (warm / hot start). In the Static initializations observations are interpolated at t=0 to the
model grid, one should ensure that the initial conditions are appropriately balanced and begin model
forward integration. Static initializations require a “spin up” period and no data are provided on scales
smaller (e.g., convective or topographic scale) than that of the available observations.
Dynamic initializations utilize some means of model “spin up” to ensure that local circulations
are represented at the initial forecast time. Typically accomplished via forecast cycling. A short-term
forecast from an earlier model simulation is used as a “first guess” for the desired simulation. Available
observations are assimilated, whether at specified times (three-dimensional) or through time (four-
dimensional), to improve the “first guess” analysis. The model is then integrated forward in time. Fig.18
shows the Illustration of forecast cycling.

28
Fig.18. Illustration of forecast cycling, where the model is initialized and a new forecast is launched every
6h. The vertical dashed lines show that output from a previously initialized forecast is used as a first
guess (FG) for an objective analysis that employs observations (Obs) that are available within a time
window around the initialization time. The cycle number is shown at the initial time of each forecast.

The need of Parameterization:

Expressing an unknown parameter into a known parameter is known as parameterization. NWP models
cannot resolve weather features and processes that occur within a single model grid box. This is due to
various surface features that cause: a) Friction that is large over tall trees, b) Turbulent eddies created
around buildings or other obstacles, c) Much less surface friction over open areas.
A model cannot resolve any of these local flows, swirls, or obstacles if they exist within a grid box.
However, the model must include or take into account the effects of these surfaces on the low-level
wind flow with a single number or expression that goes into the friction (F) term in the forecast wind
equation. Taking into account these effects without actually simulating them is called a
PARAMETERIZATION. There are many complex processes in the atmosphere that need to be
parameterized in NWP models for example, radiative processes, cloud processes, turbulence etc. The
number and type of parameterizations that are used depend on the model resolution and the type of
model used.

The processes shown in Fig.19 are often parameterized as they cannot be explicitly predicted in full
detail by model forecast equations in spite of the model resolution and their effects are important to the
simulation.

Fig.19. The processes that require parameterization

Numerical solution of equations

29
 The five prognostic equations used for forecasting are:

du ∂p
=−α +fv …………………………..(1)
dt ∂x

dv ∂p
=−α −fu …………………………(2)
dt ∂y

dw ∂p
=−α −g …………………………..(3)
dt ∂z

1 dρ ∂ u ∂ v ∂ w
= + + …………………………..(4)
ρ dt ∂ x ∂ y ∂ z

1 dp
p dt
=−γ +
[
∂u ∂ v ∂w
+
∂x ∂ y ∂z ] …………………………..(5)

The first three equations are equation of motion and the other two are equations of continuity and law
of conservation of energy for an adiabatic process. There are five dependent variables u,v.w,p and ρ
and five equations. So they can be solved.

Finite difference scheme to solve the equations

Using the Taylor’s equation we can write

❑
∆ x ∂ (f ( x )) ∆ x
2
∂
f ( x + ∆ x ) =f ( x )+ f (x ) + +… .. ❑
∂x 1! ∂x 2
2! ❑

Taking first two terms in the right hand side, we can write

∂ f ( x+ ∆ x )−f ( x )
f ( x )=
∂x ∆x

The finite difference grid is as given Fig.20.The atmosphere is represented as a two- dimensional spatial
grid defined on a specified map projection. Grid points are generally evenly spaced (or nearly so).But
exceptions are in the case of latitude-longitude grids and adaptive grids. Spatial derivatives are solved
using Taylor series- derived finite difference approximations.

30
 Fig.20. Finite difference grid

∂p
For example, to get at middle point ‘O’ we can write:
∂y

∂ p p 2− p 4
=
∂y 2d
Similar expressions can be written for all other five variables. Thus incorporating all these values at each
grid point, future time value is known.

Filtering of Sound and Gravity waves

Sound waves are longitudinal waves which are generated because of alternating adiabatic
compression and expansion of the medium. This arises because of pressure gradient. In the case of
hydrostatic atmosphere, the vertical pressure gradient cannot be influenced by adiabatic compression.
Hence the vertically propagated sound waves are filtered under the hydrostatic approximation.
However, horizontally propagating sound waves are eliminated by the assumption that at the lower
dp
boundaryω= =0 .
dt
After the sound waves are filtered, the prediction equations, such as momentum equation, continuity
equation and hydrostatic thermodynamic equation are written as:

( ∂∂t +V . ∇) V + ω ∂V∂p + fK ×V =−∇ φ……………………(1)

∂ω
∇ .V + =0 ……………………(2)
∂p

( ∂∂t +V . ∇) ∂∂ Vp +σω=0 ……………………(3)

31
 −α ∂ θ
Here σ = and the lower boundary condition is ω=0 . The above set of equations along with the
θ ∂p
hydrostatic and equation of state are called primitive equations. They contain the gravity waves as
solutions.
The gravity waves are waves in which buoyancy acts as the restoring force on air parcels
displaced vertically. A disturbance is propagated because of alternate horizontal convergence and
divergence in individual columns of fluid. A divergent horizontal velocity field which changes in time is
essential for the propagation of gravity waves. Neglecting the local rate of change of the horizontal
divergence in computing the balance between mass and velocity fields is sufficient to filter out the time
dependent gravity waves. As the prediction equations of 1 to 3 do not contain the local derivatives ∇ .V
explicitly, we use the vorticity and divergence equations instead of horizontal momentum equations.
The vorticity and divergence equations are as below:

∂ζ
∂t
=−V . ∇ ( ζ + f )−ω
∂ζ
∂p
−( ζ +f ) ∇ .V + K . (
∂V
∂p )
× ∇ ω ………………….(4)

The divergence equation is

∂V
∂t
=−∇ φ+ (
V .V
2 )
−K ×V ( ζ + f )−ω .
∂V
∂p ( )
By operating ∇ .() on this equation, the divergence equation can be written as :

∂
∂t (
( ∇ . V )=−∇2 φ +
V .V
2 )
−∇ . [ K × V ( ζ + f ) ] −ω
∂
∂p
( ∇ .V )−
∂V
∂p
. ∇ ω …………..(5)

If the Left hand side of equation (5) is set to zero, all solutions corresponding to time dependent gravity
waves can be eliminated. But the equation (5) becomes an intractable diagnostic relationship between
, V ∧φ . Hence we do the scale analysis for systematic elimination of unimportant terms. Usually
neglecting the second order terms is sufficient to filter the gravity waves.
In order to exhibit the connection between scaling and filtering, it is convenient to partition the
wind as
V =V ψ +V e ………………………(6)

Where Vψ is the non divergent part and Ve is the divergent part of the wind.
If the velocity field is two dimensional, the non-divergent part can be expressed in terms of
stream function as
V ψ =K × ∇ ψ ……………………………..(7)

As ∇ .V ψ = 0, ∇ × V e = 0 and

2
ζ =K . ∇ ×V =∇ ψ ……………………………………..(8)

The horizontal divergence is small compared to the vorticity (by scale analysis) for synoptic scale
systems in mid- latitudes. That is V is quasi- nondivergent

|V ψ|≫|V e|
32
Thus to a first approximation we can replace V by V ψ everywhere in equations 4 and 5 except in terms
involving horizontal divergence.
So the resulting vorticity and divergence equations after scaling are
∂ζ
=−V ψ . ∇ ( ζ + f )−f 0 ∇ . V e ……………………………..(9)
∂t
2 2
∇ φ=−f 0 ∇ . ( K ×V ψ )=f 0 ∇ ψ ……………………………..(10)

In the above equations gravity wave terms were filtered. From equation (10), we can interpret as

ϕ K×∇ϕ
ψ= and V ψ =
f0 f0

This means the geopotential field on a constant pressure chart is approximately proportional to the
stream function field. Equation 10 also says that to a first approximation, the vorticity is the vorticity of
geostrophic wind calculated using the constant mid latitude value of the coriolis parameter.
The geostrophic vorticity equation and the hydrostatic thermodynamic equation can be written
as:
∂ 2
( ∇ ψ )=−∇ψ . ∇❑ ( ∇ 2 ψ + f )−f 0 ∂ω ………………………(11)
∂t ∂p

∂ ∂ψ
( )
∂t ∂ p
=−∇ψ . ∇ ❑
∂ψ σ
( )
− ω ………………………(12)
∂ p f0 ❑

By eliminating ω from these equations, the quasi-geostrophic potential vorticity equation as:

❑
( ∂∂t +V . ∇) q=0 …………………………….. (13)
ψ

Where q=( ∇2 ψ+ f ) + f 02 (
∂ 1 ∂ψ
∂ p σ ∂p )
We have assumed that σ is a function of pressure only. This relation states that the geostrophic
potential vorticity is conserved following the non-divergent wind in pressure coordinates.
In practice solving equation (13) is extremely difficult and so one has to go for numerical schemes.
Models have been developed which simplify the numerical structure by referring to one, two or few
layers.
One Parameter models

There are some simplest models with only one cluster level in the vertical.

The barotropic model:

In this model one assumes that there exists a level of non divergence in the atmosphere. Usually this
level is assumed as 500mb level.

33
 −∂ ω
Thus ∇ .V = =0 at 500 mb level
∂p

∂ 2
Hence equation (9) becomes ( ∇ ψ )=−∇ψ . ∇❑ ( ∇ 2 ψ + f ) …………….(14)
∂t

This equation allows to compute the evolution of the fluid at the level of non- divergence only.

b) The equivalent barotropic model:

Same type of variation of wind with height is taken in equivalent barotropic model by averaging in
the vertical. But one level is only considered in the model. In this model the non-divergent wind is
expressed as:

V ψ ( x , y , p )= A ( p ) ⟨ V ψ ( x , y ) ⟩

The angle bracket denotes the vertical average and A (p) is a weight function. The vorticity equation
for this model is written at the stairred level where

A ( p )= ⟨ A 2 ( p ) ⟩

ζ ∧ψ are also expressed as,

ζ =A ( p ) ⟨ ζ ⟩∧ψ =A ( p ) ⟨ ψ ⟩ …………..(15)

To the extent that equation (15) is correct, the results from the solution of vorticity equation can be
utilized to deduce the evolution of ψ at other levels.

Two parameter baroclinic model

The one parameter models do not allow temperature advection. But the thermal advection processes
are essential to the development of synoptic systems. The barotropic models are quite effective in
predicting the evolution of mid tropospheric flow for periods upto a couple of days. For consistently
reliable forecast one should predict the development of new systems also. Hence there is a need of
baroclinic models. The simplest baroclinic model is a two layer model. The atmosphere is divided into
four discrete layers and the disposition of variable is given in the following schematic diagram (Fig.21).

34
 Fig.21.Division of atmosphere into discrete layers

Here ‘l’ is level and p is pressure surface.

For a flat ground the lower boundary condition is ω 4 = 0. Applying the vorticity equation (11) at
levels 1 and 3 and utilizing the finite difference technique, we can write:

ω2−ω 0 ω 2
( )
∂ω
∂p 1
=
∆P
=
∆P
( ∵ ω 2=0)

ω4 −ω 2 −ω 2
( )
∂ω
∂p 3
=
∆P
=
∆P
( ∵ ω 4=0)

The vorticity equations at the two levels are:

∂ f0
∂t
( ∇ ψ 1 ) =−( K ×∇ ψ 1 ) . ∇ ( ∇ ψ 1 +f ) +
2 2
ω ………………………(16)
∆P 2

∂ f0
∂t
( ∇ ψ 3 ) =−( K ×∇ ψ 3 ) . ∇ ( ∇ ψ 3 + f )−
2 2
ω ………………………(17)
∆P 2

Writing the thermodynamic equation (12) at level (2) by utilizing the approximation

ψ 3−ψ 1
( ∂∂ ψp ) =
2 ∆P
we get

∂ σ∆P
( ψ −ψ 3 )=−( K × ∇ ψ 2 ) . ∇ ( ψ 1−ψ 3 ) + f ω2 ………………………(18)
∂t 1 0

Here

1
ψ 2= ( ψ 1 +ψ 3 )
2

35
 Now equations 16, 17 and 18 constitute a close set of prediction equations.

First eliminate ω 2 between 16 and 17 by adding both equations to get

∂ 2
∇ ( ψ 1+ ψ 3 )=−( K × ∇ ψ 1 ) . ∇ ( ∇ ψ 1 + f ) −( K ×∇ ψ 3 ) . ∇ ( ∇ ψ 3 + f ) …………(19)
2 2
∂t

This equation governs the barotropic part of the flow.

2
−2 f 0
Next subtract ( 17) from (16 )and add the results to 2
times (18) to get
σ∆P

∂
∂t
[ ]
( ∇ 2−2 λ2 ) ( ψ 1−ψ 3 ) =−( K × ∇ ψ 1 ) . ∇ ( ∇ 2 ψ 1 + f ) + ( K × ∇ ψ 3 ) . ∇ ( ∇2 ψ 3+ f ) + 2 λ 2 ( K × ∇ ψ 2 ) . ∇ ( ψ 1−ψ 3 )
……………….(20)
2
2 β
Here λ = 2 is anon dimensional parameter. Equation (19) states that the local rate of
σ (∆ P)
change of vertically averaged vorticity is equal to the average of vorticity advections at the two
levels. Equation (20) gives the local rate of change of thickness of 250-750 mb layer.

The time derivatives between equations 16, 17 and 18 can be eliminated to get the omega equation
as

f0
[ ( ∇ −2 λ ) ( ω ) ]= σ ∆ P { [ ∇ ( K × ∇ ψ ) . ∇ ( ψ −ψ ) ]−( K × ∇ ψ ) . ∇ ( ∇ ψ + f )+( K × ∇ ψ ) . ∇ (∇ ψ + f )
2 2
2
2
2 1 3 1
2
1 3
2
3

} …………….(21)

Since equations 19, 20 and 21 are derived from the quasi-geostrophic system of equations, in the two
parameter model, the vorticity changes are geostrophic and the temperature changes are hydrostatic.
Due to these constraints, the omega field is to be adjusted at every instant.
Primitive Equation models

The filtered models discussed above, involves various approximations of the quasi-geostrophic system.
So there are some limitations on the accuracy of forecasting. In low latitudes, the stream function ψ can
not be obtained from the geostrophic vorticity field and can be obtained from more complicated
balance equation which is a nonlinear equation. It is difficult and time consuming to solve the balance
equation at each time step.
The following balance equation is derived from the divergence equation by putting initially
∂
( ∇ . V )=0 and ∇ .V =0 .
∂t
∇ 2 [ ϕ+ ½ ( ∇ ψ )2 ]=∇ . [( f +∇ 2 ψ )∇ ψ ] …………………….(22)

Therefore, it is worthwhile to investigate the possibility of using less restrictive modeling assumptions.
The scale analysis indicated that the hydrostatic approximation is far more accurate than other filtering
approximations. Hence we continue to assume that the motions are hydrostatic. The resulting system of

36
equations is isobaric coordinates consists of three prognostic equations (x and y components of
momentum equations and thermodynamic equation) and three diagnostic equations (the continuity,
hydrostatic ad the equation of state).
These six equations constitute a closed set with the dependent variables u,v,w, α ,θ ,∧ϕ which
are referred to as primitive equations.
Let us write the basic primitive equations in isobaric coordinate system as:

dV
+ fK × V =−∇ ϕ
dt

∂ω
∇ .V + =0
∂p

∂ϕ
+ α=0
∂p

d
Cp ( lnθ )=q̇
dt

( ) ( pαR )
R
p0 Cp
Where θ=
p

d ∂ ∂
= +V . ∇+ ω …………………(23)
dt ∂ t ∂p

The transformation from p to σ coordinates can be derived in a similar manner of transformation

from Z to P coordinates.

Isobaric and sigma coordinates derivations

we know the equations of motion in Cartesian and spherical polar coordinate systems taking z as
vertical coordinate. But modelers find it convenient to use the vertical coordinate as pressure, potential
temperature (entropy) or the vertical stability parameter (σ) and derive the equations of momentum,
continuity, vorticity etc. As density does not appear in these equations, the complicacy of its variation
can easily be eliminated by expressing these equations in a modified vertical coordinate.

ISOBARIC COORDINATE SYSTEM:

TRANSFORMATION EQUATIONS:

Since pressure and height are directly related in the atmosphere, it is common to express the
variables in terms of pressure surfaces rather than geometrical heights in meteorology. So it is essential
to replace the vertical coordinate (z) by pressure (p). Once pressure is taken as vertical coordinate
instead of z, it is called isobaric coordinate system.

37
 When ‘p’ is taken along ‘z’ direction, the pressure coordinate is not normal to the isobaric surfaces
but is exactly vertical upward on the plumb line.

Consider a cross section in x-z plane and let some arbitrary property ‘q’ is constant along the small
section of the surface at an instant of time as shown in Fig.22. Let the variable is ζ (x, z, t).

Fig.22. section in x-z plane

Then we can write

(δζ )1−2=(δζ )1−3 +(δζ )2−3

for small increments we can write

ζ 2 −ζ 1 ζ 3−ζ 1 ζ 2 −ζ 3 Δz
= +
Δx Δx Δz Δx which can be written for small increments as

( ) ( ) ( )( ) ( ) ( ) ( )( )
∂ζ
=
∂ζ
+
∂ζ
∂x q ∂x z ∂z
∂z
∂x q or
∂ζ
=
∂ζ
−
∂ζ
∂x z ∂ x q ∂ z
∂z
∂x q ……….(24)

As we can use any variable in place of ζ we can write (24) as

( ∂∂x ) =( ∂∂x ) −( ∂∂z )( ∂ x ) ……………………………………………….(25)

∂z
z q q

Suffix ‘q’ implies ‘holding q as constant’ in the operation.

For isobaric coordinate system q = p (pressure). A derivative with respect to z can be written as
∂ =∂ p ∂
∂z ∂z ∂ p and if we make the hydrostatic assumption this becomes

38
 ∂ =−ρg ∂
∂z ∂p ………………………………………………………(26)

So using (26), equation (25) may be written as

( ∂∂x ) =( ∂∂x ) +ρg ( ∂ x )

∂z ∂
z p p ∂ p …………………………………..(27)

Similarly we can write for ‘y’ and ‘t’ as

( ) ( ) ( )
∂ = ∂ + ρg ∂ z ∂
∂y z ∂y p ∂ y p∂ p ……………………………….(28)

∂ = ∂ +ρg ∂ z
( ) ( )
∂t z ∂t p ∂t ( ) ∂∂p
p ……………………………… (29)

The last term of equation (29) denotes the rate at which the height of an isobaric surface
changes with time. Equations (27) to (29) denote the inputs to transform the (x,y,z,t) Cartesian
coordinate system to (x,y,p,t) isobaric coordinate system.

ISENTROPIC COORDINATE SYSTEM:

In many meteorological phenomena particularly in the free atmosphere, the air motions are
assumed to be isentropic. So an air parcel will always have the same potential temperature ( θ ) in
order to remain on the same isentropic surface (equal entropy line). Under these circumstances it is
convenient to use θ as the vertical coordinate in place of pressure (p) or ‘z’. In other words in an
adiabatic flow, potential temperature is conserved, so isentropic coordinates are useful, for tracing the
actual paths of travel of individual air parcels.

So taking equation (27) and considering q = p (pressure) can be written as

∂ = ∂ + ∂ ∂z
( ) ( )
∂x p ∂x z ∂z ∂ x ( ) p ………..(30)

Similarly for ‘y’ coordinate can be written as

∂ = ∂ + ∂ ∂z
( ) ( )
∂y p ∂y z ∂z ∂y ( ) p ……………(31)

Multiplying equation (30) by i and (31) by j and then on adding we get

39
( ) ( ) (
i ∂ + j ∂ = i ∂ +j ∂ + ∂ i ∂ + j ∂ z
∂x ∂ y p ∂x ∂ y z ∂z ∂x ∂ y p )
∇ q =∇ z + ∂ ∇ q z
Which means ∂z where the suffix ‘q’ stands for any constant surface.

Let the constant surface ‘q’ = θ . Then the equation turns out as

∇ θ =∇ z + ∂ ∇ θ z
∂z
Let the variable is pressure ‘p’ then

∂p
∇ θ p=∇ z p+ ∇ z
∂ z θ ………..(32)

On a horizontal surface (z = H) we can write

∂p
∇ H p=∇ H p + ∇ z
∂ z θ …………..(33)

EQUATION OF MOTION UNDER ISENTROPIC COORDINATE SYSTEM

We know the equation of motion in Cartesian coordinate system as

du ∂p
=−α + fv
dt ∂x

dv ∂p
=−α −fu
dt ∂y
Multiplying first one by ‘i and second one by ‘j’ and adding both and considering

u+iv = V we can write the vector form as

dV
=−α ∇ H p−fk×V
dt ………(34)

SIGMA COORDINATE SYSTEM:

To avoid the serious disadvantage of the isobaric coordinate system, σ coordinate system has
been introduced in modeling work. The disadvantage of isobaric-coordinate system is the coincidence of

surface pressure at the lower boundary with the ground (

z0 )
which is not true.

40
 In the sigma (σ ) system, the vertical coordinate is the pressure normalized with the surface
pressure (ps).

p
σ=
∴ p s …………..(35)

This is called the static stability parameter. Thus σ here is a non- dimensional independent
vertical coordinate which decreases upward from a value σ = 1 at the ground to σ = 0 at the top of the
atmosphere. Here in σ coordinates the lower boundary condition will always apply exactly at σ =1.
Further the vertical sigma velocity is given as

.
dσ
σ=
dt

This will always be zero at the ground even in the presence of sloping terrain. So the lower
.
boundary condition in the σ system is σ = 0 at σ = 1. Thus this σ coordinate system is particularly
useful in the regions of strong topographic variations.

TRANSFORMATION OF SIGMA COORDINATES:

Since in general σ and p surfaces will not coincide, partial derivatives evaluated at constant σ
will not be equal to the analogous derivatives evaluated at constant p. The transformation of
coordinates to σ can be achieved in the same previous lines.

Take eqn (32) and replace θ by σ , then the transformation equation can be written as:

∂p
∇ z p=∇ σ p+ ∇ z
∂z σ
We know from (35) p = σ ps. So replace p by σ ps) and z = H, then

∂p
∇ H p=∇ σ ( p s σ ) + ∇ z
∂z σ

g ∂p
∇ H p=∇ σ ( p s σ ) − ∇ σ z (∵ =−ρg=−α )
α ∂z
This can be further written as

α ∇ H p=α ∇ σ ( ps σ ) −∇ σ ( gz )

α ∇ H p=α ∇ σ ( ps σ ) −∇ σ (φ )
Or …………..(36)

41
The equation (36) is the transformation equation for pressure gradient term of momentum equation
(34). So on substituting eqn (36) in (34) we get

dV
=−α ∇ σ ( pS σ )+∇ σ φ−fk×V )
dt

RT α=
RT
p S =RρT =
Since α and so replacing PS we can get

dV RT
=− ∇ ( p σ )+∇ σ φ−fk×V )
dt PS σ S ……….(37)

This is the momentum equation in sigma coordinates.

Please recall in isobaric coordinates we put w ( p 0 )=−ρ 0 gw ( Z 0 ). We assume that the height of
the ground Z 0 is coincident with the pressure surface p 0 (= 1000mb). However, this is not strictly true, as
pressure surface will not coincide with the ground even when the ground is level. So there is a need of σ
p
system of coordinate is used in the vertical in numerical modeling such that σ = where ps denotes the
ps
pressure at the surface. In sigma coordinate system the lower boundary condition is applied at σ =1 .
dσ
The vertical velocity in σ coordinates is σ̇ = will always be zero at the ground as σ =1 is constant
dt
even in the presence of sloping terrain. Hence in σ coordinates the lower boundary condition is σ̇ = 0
at σ =1.
σ ∂
Using the definitions of equations 25 & 32 , one can write: ∇ p ()=∇ σ ()− ∇ p () ,
ps ∂σs

where ( )denotes for the use of any parameter.

So the system of equations in σ coordinates similar to equation 37 can be written as:

dV σ ∂ϕ
+ fK × V =−∇ ϕ + ∇ p …………(i)
dt ps ∂σ s

∂ σ̇ ∂ p s
∇ . ( p s V )+ p s + =0 ………………….(ii)
∂σ ∂t

∂ ϕ −RT
= ………………(iii) ………………(38)
∂σ σ

d
Cp ( lnθ )=q̇ ……………………(iv)
dt

( )
R
p0 Cp
Where θ=T …………(v)
σ ps

42
 d ∂ ∂
= +V . ∇+ σ ……………..(vi)
dt ∂ t ∂σ

The set of six equations in (38) contain seven dependent variables u, v, σ̇ , θ , T ϕ and ps . So another
equation is necessary to solve these variables. This equation is surface pressure tendency equation. It
can be obtained by integrating the continuity equation vertically by using the boundary conditions σ̇ = 0
at σ =0∧1.
1
∂ ps
The surface pressure tendency equation is =−∫ ∇ . ( p s V ) dσ
∂t 0
This equation states that the rate of increase of the surface pressure at a given point is equal to the
mass convergence in a unit cross section of column above that point.

Two layer primitive equation model (PE Model)

In the primitive equation model the static stability parameter σ is a function of space and time
∂θ
whereas in earlier models, σ (p) is usually a specified parameter. It can be shown that if is
∂p
specified as a function of pressure alone, that is σ =σ (p), it is not possible to formulate a PE model
∂θ
which satisfies the requirements of energy conservation. Hence it is necessary to compute explicitly
∂p
at each time step. Hence temperature must be predicted at a minimum of two levels rather than at one
level as in the case of Z-level quasi-geostrophic model.
The schematic diagram of Z-level PE model with variables is as below:

The momentum equation and the thermodynamic equations are applied at levels 1 and 3. The
vertical differencing of these equations requires knowledge of sigma vertical motion at level 2. σ˙ 2 is
obtained diagnostically through the use of continuity equation as follows:

43
 σ̇ 2−σ̇ 0
( )
∂ σ̇
∂σ 1
=
1
=2 σ̇ 2
2

σ̇ 4 −σ̇ 2
( )
∂ σ̇
∂σ 3
=
1
=−2 σ̇ 2
2

From the continuity equation in (38 vi) we can write:

∂ ps
+ ∇ . ( ps V 1 ) +2 p s σ̇ 2=0 ………………..(39)
∂t

∂ ps
+ ∇ . ( ps V 3 )−2 ps σ̇ 2=0 …………………(40)
∂t

Adding (39) and (40) we can write:

∂ ps 1
+ ∇ . [ p s ( V 1 +V 3 ) ]=0
∂t 2

Subtracting (40) from (39) we can get:

−1
σ˙ 2= ∇ . [ ps ( V 1−V 3 ) ]
4 ps

In a similar manner ϕ 1 and ϕ 3 can be obtained diagnostically by vertical integration of hydrostatic

relation in equation (38).

In short the following steps may be followed for the forecast using PE model.

1. Suitable finite difference analogues of the momentum and thermodynamic equations at levels 1
and 3 and the surface pressure tendency equation are written.
2. These prognostic equations are utilized to obtain tendencies of V 1 ,V 3 , θ1 ,θ 3 , and ps .
3. Extrapolating the tendencies ahead using a suitable time differene scheme.
4. The new variables of V 1 ,V 3 , θ1 ,θ 3 , and ps are used to determine σ˙ 2, ϕ 1 and ϕ 3 diagnostically.
5. Repeat steps 2 to 4 until the desired forecast time is reached.

It must be noted that the PE model contains the mechanism for both gravity wave and horizontal
acoustic wave propagation. It is necessary that the time increment be kept small enough to satisfy the
C Δt 1
condition: < where C is horizontally propagating sound wave which is the fastest moving wave
d 2
solution, d is horizontal resolution and Δ t is time increment. Therefore the time steps in PE models are
considerably small ( ≈ 10 minutes).

Initialization Scheme

44
Initialization scheme involves the use of accurate initial data. One reason for the failure of Richardson is
due to lack of accurate initial data. Though the data coverage is vastly improved since then actual
measured wind and pressure data is still lacking as most of the data is derived from remote sensing
satellite data.

To introduce the initialization method, consider the relative magnitudes of the various terms in
the momentum equation in isobaric coordinates:

∂V ∂V
+( V . ∇ ) V + ω =−fK ×V + ∇ φ
∂t ∂p

For synoptic scale motions in mid-latitudes, we have seen that the wind and pressure fields are in
approximate geostrophic balance. Thus the acceleration following the motion is measured by the small
difference between the two nearly equal terms Coriolis and geopotential ( −fK ×V ∧∇ φ¿ . Although φ
can be determined with good accuracy using observed data, observed winds are often 10 to 20% in
error. Hence if wind is used as initial data, the Coriolis force works out to 10 to 20% in error. Since the
acceleration is normally 10% of the Coriolis force, the acceleration comes to 100% in error (10x10). Such
spurious acceleration leads to poor estimates of initial pressure and velocity tendencies. This further
produces large amplitude gravity wave oscillations.

To avoid this problem, let us derive the balanced equation as follows. In the divergence
∂
equation at t = 0, put +V . ∇=0 and V . ∇=0
∂t

[
∇2 ϕ +
V .V
2 ]
=∇ . [ K ×V ( ζ + f ) ]

Since we have assumed that the initial velocity is non divergent, we can replace V by ψ by letting

V ψ = K × ∇ ψ then we get

∇ 2 [ ϕ+ ½ ( ∇ ψ )2 ]=∇ . [( f +∇ 2 ψ )∇ ψ ]……………. (41)

This equation is called balance equation. Equation (41) is linear in ϕ and non linear in ψ . In mid-latitudes
ϕ is known and hence ψ is solved with the help of equation 41. The ψ has two solutions. One solution
corresponds to the gradient wind balance and the other one is anomalous. Thus as one ψ is known, the
initial velocity field V ψ can be computed. This velocity is in gradient wind balance with the geopotential
field. The initial local accelerations are thus very small. Since the initial rate of change of divergence is
zero, the large amplitude gravity wave oscillations are not excited by the initial conditions.

Present Status of NWP

So far, we have traced only the early developments of the subject of NWP before, say, 1970. Since then,
the field has expanded so much that it is practically impossible to cover all aspects of the later

45
developments in this brief survey. Computing machines are now millions of times faster; a development
which has allowed meteorologists to take up models with much higher spatial resolution over a global
domain both in the horizontal and the vertical than before. Data coverage has improved a great deal.
Apart from conventional observational network, geostationary and polar -orbiting satellites now
regularly report data from different parts of the global atmosphere. Four-dimensional data assimilation
schemes have been devised which prepare initial data directly for integration through a 6-hr data
assimilation prediction scheme. With improved computational schemes in place (at many centers, finite-
difference models have been replaced by spectral models), the global forecast centers are now able to
issue extended-range forecasts of mass and flow fields over the globe for periods ranging from a day to
a week or more.
Verification statistics of present-day numerical forecasts appear to suggest that by and large these
improved numerical models perform better than those based on persistence or other statistical or
subjective methods, especially at extended ranges. For example, according to a study by Leith (1978),
quoted by Holton (1979), the root-mean-square error in a 24-hr forecast of 500 mb height was 40m by
the primitive equation model as against 65m by persistence, the initial height error being 20m (now less
than 10 m). With extension of the forecast period to 48 h, the errors in the two methods grew to about
60m and 95m respectively. Experiments on how the growth of errors limits the inherent predictability of
the atmosphere have been carried out by a number of workers using primitive equation forecast
models. Results of these experiments appear to indicate that the theoretical limit for any useful
forecasts of synoptic-scale motion is probably one to two weeks. Actual predictive skill of present-day
models falls far short of this theoretical limit due to various sources of error, some of which are
computational and some physical. Among the main sources of error are the replacement of derivatives
of the dependent variables by finite-differences, inadequate representation of mountain effects,
boundary layer processes, frictional dissipation, cloud and precipitation, and horizontal and vertical
resolution of different scales of motion. Further, it is to be understood that the above error statistics
apply only to forecasts of variables of the atmosphere which are continuous functions of space and time,
such as pressure, temperature, geopotential height, wind components, etc. They do not apply to
forecasts of actual weather phenomena experienced in daily life, such as fog, cloud, rain,
thunderstorms, tornadoes, etc, which may or may not occur during the forecast period. However,
practicing forecasters give due weightage to NWP model output statistics and interpret them suitably
along with other guiding materials available to them in assessing the probability of occurrence or non-
occurrence of any of these weather elements.

That is how the art or science of weather forecasting stands to-day. However, there is no doubt that
attempts to develop numerical models with better coverage of initial data over land and ocean,
improved numerical techniques, and higher resolution in space and time have led to better
understanding of the dynamics of the atmosphere during the last several decades and the trend will
undoubtedly continue. So far as long-range forecasting is concerned, there has been increasing
realization in recent years of the role played by oceans in atmospheric circulation and the atmosphere in
ocean circulation. This has led to development of interacting coupled ocean-atmosphere models at
many of the world meteorological centers (e.g. Saha et al., 2006) and the trend holds out a great
promise for future climate forecasts.

References:

46
Bjerknes, V., 1904: Das Problem der Wettervorhersage, betrachtet vom Standpumkte der Mechanik und
der Physik. Met. Zeit., 21, 1-7. [Translated and edited by Volken, E. and S. Brönnimann, 2009: Met.
Zeit., 18, 663-667, doi:10.1127/0941-2948/2009/416.]

Charney, J., R. Fjörtoft and J. von Neumann, 1950: Numerical integration of the barotropic vorticity
equation. Tellus, 2, 237 – 254, doi:10.1111/j.2153-3490.1950.tb00336.x.

Holton JR (1979): An Introduction to Dynamic Meteorology, 2nd edn. Academic Press, New York,
p 391
Leith C.E (1978): Objective methods of weather prediction. Ann Rev Fluid Mechan 10.

Lorenz, E.N., 1963: Deterministic non-periodic flow. J. Atmos. Sci., 20, 130-141.

Lynch, P., 2006: The emergence on numerical weather prediction: Richardson’s dream. Cambridge
University Press. 290 pp.

Monin, A. S., and A. M. Obukhov, 1954: Basic laws of turbulent mixing in the surface layer of the
atmosphere. Tr. Akad. Nauk. SSSR Geophiz. Inst. 24, 163-187. [English translation available online at
http://mcnaughty.com/keith/papers/Monin_and_Obukhov_1954.pdf.]

Richardson, L. F., 1922: Weather prediction by numerical process. Cambridge University Press. 262 pp.

Saha S, Nadiga S, Thiaw C, Wang J, Wang W, Zhang Q, van den Dool HM, Pan H.-L, Moorthi S, Behringer
D, Stokes D, PenaM, Lord S,White G, EbisuzakiW, Peng P, Arkin P, Xie PP (2006): The NCEP climate
forecast system. J Climate, 19:3483–3517.

Sharma, N., P. Sharma, D. Irwin, and P. Shenoy, 2011: Predicting Solar Generation from Weather
Forecasts Using Machine Learning: Proceedings of the 2nd IEEE International Conference on Smart
Grid Communications, Brussels, 17-20 October, pp. 32-37.

47