Increasing the Interoperability of an Earth System Model:

Atmospheric-Ocean Dynamics and Tracer Transports

NCC4-624

 

Milestone F

 

ÒUpgrade the parameterization of the planetary boundary layer (PBL) used in the AGCM to a version with multiple layers without loss of performance compared to the one in Milestone E.  Provide scaling curves. 

Documented source code made publicly available via the WebÓ

 

1.              Introduction

 

In the UCLA Earth System Model (ESM) version with which we started this project, the atmospheric component (UCLA atmospheric general circulation model, AGCM) includes a variant of the bulk parameterization for the variable-depth PBL processes proposed by Deardorff (1972).  This is relatively simpler and physically realistic.  In the model, the layer next to the lower boundary is designated as the variable-depth PBL and is assumed to act as a well-mixed layer (Suarez et al. 1983). The PBL top is represented by a coordinate surface, whose height is predicted from the mass budget using a parameterized formulation of entrainment (detrainment). The PBL temperature, moisture and wind fields are predicted using the parameterized surface fluxes and the fluxes associated with entrainment (or detrainment) through the PBL top. A stratocumulus cloud sub-layer can develop along the PBL top if this is higher than the condensation level.  In this case turbulence is primarily driven by the downward buoyancy due to the radiative cooling near the cloud top (Lilly, 1968), and the cloud-topped PBL can be maintained even without positive buoyancy due to the surface heating. The parameterization also incorporates the process that can destroy the cloud sub-layer by entrainment of drier air into the PBL.  The merits of this unique approach to PBL parameterization are demonstrated by the highly successful performance of the UCLA AGCM in the simulation of formation and seasonal variations of marine stratocumulus (Li et al. 1999; Mechoso et al 2000). The AGCM is also fairly successful in capturing the diurnal cycle of precipitation over regions of strong monsoonal circulations (Fig. 1).

 

The approach to modeling the PBL as a single well-mixed layer is not without demerits.  For example, it can provide a poor representation of low-level vertical shear when the PBL is deep. Also, there is no room for multiple cloud layers within the PBL. Of higher relevance to tropical convection simulated by the UCLA AGCM is that using such an approach, cloud base conditions in the Arakawa-Schubert cumulus parameterization (Arakawa and Schubert 1974) are determined by the PBL mean quantities only. Under separate funding we are implementing in the AGCM a new formulation (Konor and Arakawa 2001) that introduces multiple layers within the PBL in the framework of the modified sigma-coordinate (i.e., the upper and lower boundaries of the PBL are still coordinate surfaces). The PBL processes are then formulated following a hybrid approach, in which the effects of large-scale eddies and small-scale eddies are formulated separately. For the large-scale eddies, a relatively well-mixed vertical structure is assumed for conservative thermodynamic variables and a bulk approach is applied to the properties vertically averaged over the entire PBL so that the formulation is non-local. For the small-scale eddies, a K-closure formulation, which tends to produce a well-mixed vertical structure, is applied between the multiple layers within the PBL. Unlike the single-layer approach, the hybrid multi-layer parameterization allows for vertical shears and deviations from well-mixed profiles within the PBL.  These deviations are expected to be small for thermodynamic conservative variables on a convectively active PBL, while they can be significantly large for a convectively inactive PBL. The formulation of processes highly concentrated near the PBL top remains tractable in this approach, while it is more directly applicable to a variety of PBL regimes, particularly the diurnally changing PBL over land.

The work to be performed in order to achieve Milestone F consists of upgrading the single-layer PBL parameterization used in the AGCM to the multiple-layer version without loss of performance.  We have completed this work interpreting ÒupgradingÓ as the changes in the code that allow for a single or multi layer PBL within the same framework.  Namely, the single and multiple layer PBL parameterizations appear as options to be set up at run time.  The component of the code directly affected by the changes is AGCM-Dynamics. The AGCM-Physics part of the code is unchanged. A composite single layer PBL is computed from the multilayer PBL, and computation of physical column processes other than those of the PBL takes place as before.

 

The dynamical equations integrated in this new formulation (Konor and Arakawa 2001) are presented in section 2. Simulation of the effects of turbulent diffusion due to small scale eddies and to convective eddies inside the PBL on the prognostic variables is described in section 3. The diffusion coefficient for the small scale eddies is computed through a local Richardson number formulation, following Louis (1979) and Holstlag and Boville (1993), (section 3.3).  Section 4 presents a selection of preliminary results.

 

 

 

 

 

 

 

 

 

 

 

 

2.      Governing equations

 

2.2    Continuous case.

 

         This section describes the prognostic equations using the sigma vertical coordinate system, defined in Arakawa Suarez (1983).

         We define pB as the pressure at the top of the PBL, and pS as the pressure at the Earth surface. Then, inside the PBL, the sigma vertical coordinate is defined as

 

                                              for                               (2.1.1a)

 

         In the free atmosphere we consider two regions. One, between the PBL top and an intermediate pressure p = pI, and the other, between pI and the atmosphere top, that is considered to be at the level p = pT. The definitions of the sigma coordinate in these regions are

 

                                                 for                                  (2.1.1b)

 

and

 

                                                 for                                  (2.1.1c)

 

         We currently choose pT = 1 hPa and pI = 100 hPa.

         As a consequence of these definitions, s=1 for the PBL top (p = pB), which then corresponds to a coordinate surface. Also the top of the atmosphere and the Earth surface are coordinate surfaces, with s = - 1 and s = 2 respectively. From these definitions the pressure can be obtained as

 

             for , were     (2.1.2a)

 

                               for , were               (2.1.2b)

 

                               for , were               (2.1.2c)

 

Note that pstrat is constant.

 

The continuity equation can be written as

 

                                                                                        (2.1.3)

 

where p is either psurf, ptrop or pstrat.

         The momentum equation, for layers within the PBL, is

 

                                                           (2.1.4a)

 

where the last term is the vertical convergence of the turbulent flux of momentum, which will be discussed in section 3.

         In the free atmosphere (p > pB), the momentum equation is

 

                                                                        (2.1.4b)

 

         The geopotential F = gz is diagnosed from the hydrostatic equation,

 

                                                              ,                                                (2.1.5)

 

where P is the Exner function, defined as

 

                                                                                                         (2.1.6)

 

         The thermodynamic equation in terms of the potential temperature q is, within the PBL,

 

                                              (2.1.7a)

 

where Gq is the contribution of the turbulent fluxes to the tendency of q, and will be discussed in the next section, and Q is the heating rate.

         In the free atmosphere, if ,

 

                                  ,                  (2.1.7b)

 

if ,

 

                                                                                  (2.1.7c)

 

         The continuity equation for the water mixing ratio (r) inside the PBL combines water vapor mixing ratio (q) plus the liquid water mixing ratio (l). This is because we admit that PBL turbulence allows the air to hold both phases of water, and large-scale precipitation processes would not occur there unless the pressure at the condensation level is greater than the surface pressure (condensation level bellow the earth surface), then large scale precipitation is computed, in order to make the condensation level equal to the surface level.

         The continuity equation for r within the PBL is

 

            , with            (2.1.8a)

 

where Gr is the contribution of the vertical convergence of the turbulent flux of r, and C is the condensation rate.

         In the free atmosphere, the equation is formulated in terms of the water vapor mixing ratio (q), since if r > qsat, (qsat is the saturation water vapor contain) liquid water would precipitate (large scale precipitation processes). If

 

                                                    (2.1.8b)

 

and if ,

 

                                                                           (2.1.8c)

 

2.2.      Discretized equations

 

         In this section we discuss the vertical discretization of the equations in section 2.1. A similar discussion for a hybrid vertical coordinate is found in Konor Arakawa 2001.

the atmosphere is vertically divided in layers (Fig. 2), from k = 1 (upper most  layer) to k = kmax-1 (lower most layer). We use half integer indexes for labeling the layer interfaces; k+1/2 is the interface between layer k and layer k+1. The top of the atmosphere is the first layer interface, with index 1/2, and the earth surface, the last layer interface, with index kmax-1/2.

 

         The level  (atmosphere top) will have vertical index equal to 1/2, the level  will correspond to the interface kstrat + 1/2, the level (PBL top) will correspond to the interface ktrop - 1/2, and the level  (earth surface) will correspond to the interface kmax-1/2.

 

 

                                    Atmopshere top

                                    _______________________________          1/2

                  ÉÉÉÉÉÉÉÉÉÉÉÉÉÉÉ..          1                                  first layer

                  _______________________________          3/2                               layer interface

 

                                    _______________________________          k – 1/2

                  ÉÉÉÉÉÉÉÉÉÉÉÉÉÉÉ..          k                                  generic layer

                                    _______________________________          k – 1/2

 

                  _______________________________          ktrop – 1/2

                                    ÉÉÉÉÉÉÉÉÉÉÉÉÉÉÉ..          ktrop                            uppermost PBL layer

                  _______________________________          ktrop + 1/2

                                    ÉÉÉÉÉÉÉÉÉÉÉÉÉÉÉ..

                  _______________________________

                                    ÉÉÉÉÉÉÉÉÉÉÉÉÉÉÉ..

                                    _______________________________          kmax – 3/2

                                    ÉÉÉÉÉÉÉÉÉÉÉÉÉÉÉ..          kmax – 1

                  _______________________________          kmax – 1/2

                                    Earth Surface

 

Figure 2: Scheme of the layers, the layer interfaces and the vertical indexes.

 

 

         Horizontal velocity, temperature and water vapor mixing ratio will be computed at the layers, while vertical velocities (Ds/Dt) will be computed at the layer interfaces, this arrangement consists in the Lorenz vertical grid (Lorenz 1963).

 

We discretize the mass continuity equation by

 

              for    (2.2.1)

 

where p will be either psurf, ptrop or pstrat, according to the vertical interval in which the layer is, and

 

                                                                                                (2.2.2)

 

         In the interface between the PBL and the free atmosphere, we compute the vertical mass flux as

 

                                                                      (2.2.3)

 

where we take E > 0 if there is mass entrainment, or D > 0 if there is detrainment. MB is the cumulus mass detrained from the PBL to the cumulus clouds from their bases. Summation in the vertical of equation (2.2.1) for all the free atmosphere layers yields

 

                       (2.2.4a)

 

where pstrat = constant was used. Summation (2.2.1) for the PBL layers gives

 

                                                        (2.2.4b)

 

and summation for all the atmosphere layers gives

 

.      (2.2.4c)

 

         Equations (2.2.4a-c) allow to predict ptop , psurf, pS and therefore, also pB. and ptop. On the other hand, partial summation gives, within the PBL,

 

           (2.2.4c)

 

which allows to diagnose vertical velocity in the interfaces within the PBL.

         The vertical momentum advection is

 

 (2.2.5a,b,c)

 

         The term vB+ is linearly extrapolated from the layers above the PBL. In the lowest layer of the free atmosphere,

 

   (2.2.5d)

 

and in the other layers, for .

 

                           (2.2.5e)

 

where p is either ptrop or  pstrat, according to the vertical interval where the layer is.

The geopotential, (used in the momentum equation) is computed with the discrete form of the hydrostatic equation,

 

                 (2.2.6)

 

         This discretization maintains the constraint of global total energy conservation (Arakawa Suarez 1983).

 

         Inside the PBL, the pressure gradient force that appears in the momentum equation is computed in a similar way as proposed by Konor and Arakawa (2001), for hybrid vertical coordinate,

 

(2.2.7a)

 

for , with . For the free atmosphere,

 

             for (2.2.7b)

 

                                            for                             (2.2.7c)

 

         In order to obtain the discrete forms of vertical fluxes of potential temperature, we first define qk+1/2, as the value of q interpolated at the interface between layers k and k+1,

 

 for                         (2.2.8)

 

         The values of the Exner function at the layer are computed for its values at the layers interfaces, by

 

                                                                                (2.2.9)

 

with

 

                                                                                            (2.2.10)

 

         Within the PBL, the convergence of potential temperature is discretized as

 

(2.2.11a)

 

         In the lowest layer of the free atmosphere,

 

              (2.2.11b)

 

for ,

 

                             (2.2.11c)

 

and for ,

 

                                                      (2.2.11d)

 

         The convergence of the vertical moisture flux is discretized within the PBL by

 

   (2.2.12a)

 

were rB+ is a value extrapolated form the layers above the PBL, and rk+1/2 is the value of r at the layer interfaces, interpolated from the adjacent layers, with

 

         In the free atmosphere, for ,

 

                               (2.2.11c)

 

and for ,

 

                                                          (2.2.11d)

 

 

 

 

 

3.         Turbulent Fluxes

 

         If y is either v, r, q or h, the moist static energy, defined as , with  for non-saturated air, and  for saturated air, the turbulent flux of y is written as , where is the flux due to convective eddies, and is the flux due to small scale turbulent eddies. For q, the turbulent fluxes due to both convective and small scale eddies are computed from the respective turbulent fluxes of r and h, in a way that depends whether the air is not saturated or saturated. For non saturated air,

 

                                           (3.1)

where L is the latent heat of water and Fq is equal to Fr. For saturated air,

 

, with .                          (3.2)

 

3.1       Turbulent diffusion due to small scale eddies

 

The diffusion of water mixing ratio and moist static energy due to small scale eddies is given by the diffusive equation

 

                                                (3.1.1)

 

were K is the turbulent diffusion coefficient, to be discussed in section 4, and  is the turbulent flux of y due to small scale eddies, . The discrete form of (3.1.1) is

 

                                        (3.1.2)

 

where

 

                                           (3.1.3)

 

Note that as we use the time level n+1 for computing the fluxes of the right side discrete prognostic equations, we have an implicit time scheme.

We consider

                                           (3.1.4)

 

We also must note that, whatever the formulation of F, if the constraint (3.1.4) applies, then

 

                                             (3.1.5)

 

therefore,

 

                                         (3.1.6)

 

and then we conclude that mass weighted average of the diffused quantities must be conserved, that is,

 

                                         (3.1.7)

 

We next develop the discrete form of the equation (3.1.1),

 

                 (3.1.8)

 

These equations can be also written in the following way, that gives a linear system of equations, for computing time level n+1 from n.:

 

(3.1.9)

 

We define A(k) and B(k) as

 

                      (3.1.10)

 

therefore

 

              (3.1.11)

 

in matrix notation,

 

(3.1.12)

 

If we have an exact solution of this system, constraint (7) will apply. After obtaining approximate solutions of the system (3.1.12), we compute a second approximation of the n+1 time level, by

 

                                 (3.1.13)

 

with

 

                                       (3.1.14)

 

and now the constraint of mass weighted average conservation will apply exactly.

In these way we diffuse the water mixing ratio and the moist static energy. Potential temperature is diffused computing the vertical convergence of its turbulent flux, which in turn is computed from the fluxes of h and r given by the last equation, discussing whether air is saturated or not:   

 

                                  (3.1.15)

 

if the air is not saturated. In this case, Fq = Fr. If the air is saturated,

 

                           (3.1.16)

 

3.2       Turbulent diffusion due to convective eddies

 

Fluxes of h and r are computed at the Earth surface using the bulk aerodynamic formulas proposed in Deardorf (1972).

At the top of the PBL, turbulent fluxes are computed considering that there can be a discontinuity or ÒjumpÓ between the values of h and r in a level of the free atmosphere slighty above the PBL top, ÒB+Ó level, and its values in a level slighty bellow the PBL top, ÒB-Ò level.

When we have turbulent mass entrainment E > 0, the turbulent flux of these properties at the top of the PBL will be

 

                                            (3.2.1)

 

with y =  h or r. In the intermediate layer interfaces, we compute the fluxes as linear interpolations, in sigma, between the values at the top and at the surface,

 

                               (3.2.2)

 

For q, the turbulent fluxes are computed in terms of the fluxes of h and r, discussing if the air is saturated or not, in the same way as exposed for the small eddies diffusion. Then, for y (y being q or r), the prognostic equation within the PBL comes to

 

, for           (3.2.3a)

 

in the case k = kmax  - 1, the flux at k + 1/2 is the surface flux. For k = ktrop, (uppermost PBL layer), equation is

 

             (3.2.3b)

-

in this case contribution of the turbulent flux from ktrop – 1/2 is not taken into account, since it was already accounted in the vertical flux of equation (2.2.11a) for q, and (2.2.12a) for r.

 

3.3       Computation of the coefficient for small scale eddies diffusion.

 

The diffusion coefficient due to small scale eddies computed through a local Richardson number formula proposed by Louis (1979), revised by Holstlag and Boville (1993). This formula states that

 

                                          (3.3.1)

 

with

 

                              (3.3.2)

 

                                          (3.3.3)

 

where k = 0.4,  l = 100m.

 

If Ri < 0 (unstable case),

 

                                        (3.3.4)

 

with b  9.4, and

 

                    (3.3.5)

 

with C* = 5.3 for the diffusion of heat and moisture, and C* = 7.4 for the diffusion of momentum.

In the stable case,

 

                                         (3.3.6)

 

with bÕ = b/2 = 4.7.

 

The Richardson number used in these formulas is defined as

 

                                             (3.3.7)

 

and qe is the equivalent potential temperature, defined as

 

                                                                                                              (3.38)

 

where q is the water mixing ratio, L is the specific latent heat of water, and Tc is the temperature at the condensation level, defined as

 

                                                                                                       (3.3.9)

 

Pc is the Exner function computed from the condensation level pc, defined as

 

                                                                                 (3.310)

 

subindex S denotes earth surface conditions, while qB0 is the saturation water vapor mixing ratio at a pressure pB, and at a temperature extrapolated from the subcloud layer to the PBL top (see Fig. 4 at Suarez et alt. 1983). qPBL is computed form the averaged moist static energy, hPBL, and water mixing ratio, rPBL, and the earth surface geopotential FS, by

 

                                                                                (3.3.11)

 

 

4.         Preliminary Results

 

In this section we show some preliminary results of a simulation performed in low resolution (4o latitude, 5o longitude, 14 layers in the free atmosphere and 4 layers in the PBL). The initial conditions were given at October 1st, and the simulation extended through one year.

         Fig. 3 shows a simulated January monthly mean of the global precipitation field. This simulated field shows several features of the observed January climatological precipitation, like the winter storm tracks in the Northern Hemisphere, the Intertropical Convergence Zones (ITCZs) in the tropical Pacific and Atlantic Oceans. In the Southern Hemisphere the simulation is also able to reproduce reasonably the Southern Pacific Convergence Zone (SPCZ) and the South Atlantic Convergence Zone (SACZ) as well as the summer precipitation in the continents.

         Fig.  4 and Fig. 5 shows respectively the simulated January means for the sea level pressure (SLP) field and the zonal mean of the zonal wind, which respectively have patterns reasonably comparable with those observed.

 

 

Figure 3: January mean precipitation simulated using the multilayer PBL. Contour interval is 2mm/day

 

Fig. 4: January mean SLP field simulated with the multilayer PBL. Contour interval is 4 hPa.

 

 

Fig. 5: Zonally averaged January mean the zonal wind simulated with the multilayer PBL. Contour interval is 5 m/s

 

5.         Code Performance

            In order to evaluate the performance of the code while including the new multi layer PBL, several tests were done in the AMES machine Chapman, of the type SGI 0-3000. We used the high resolution version of the model (2o of latitude, 2.5o of longitude, and 28 layers in the free atmosphere), and we considered both a single layer PBL and a four layers PBL. The number of processors used went from 32 to 512. Table I shows the wall clock time per simulated day for the multilayer simulations and the single layer simulations, for each number of processors considered. We found that the multi layer PBL increases wall clock times per simulated day by a percentage between 5 and 7%. Fig. 6 shows the performances in terms of simulated days per wall clock hours.

 

 

 

 

Number of processors

Wall clock time per simulated day.

% of increase

 

Single layer PBL

4 layers PBL

 

32

134.3

142.0

5.8%

64

89.4

94.1

5.2%

128

64.5

67.1

4.1%

256

51.1

54.6

6.8%

384

44.1

47.3

7.2%

512

45.0

47.2

4.9%

 

 Table I. Wall clock time in seconds, per simulated day, for single layer and 4 layers PBL simulations, at AMES SGI 0-3000 computer Chapman. Horizontal resolution was 2ox2.5o, 28 layers were considered in the free atmosphere.

 

 

 

Fig. 6: Simulated days per wall clock hour, for the UCLA AGCM at SGI 3000 Chapman computer. AGCM has 28 layers in the free atmosphere and the horizontal resolution is 2o of latitude and 2.5o of longitude. The higher curve is the single layer PBL case, and the lower, the four layers PBL case.

 

5.         Published Papers

A Multi-Layer Bulk PBL Parameterization for Use in Climate Models. Application to an Atmospheric General Circulation Model

 

6.         References

 

       Arakawa A., M. J. Suarez, 1983: Vertical differencing of the primitive equations for the shallow water equations, Mon. Wea. Rev., 111, 34-45.

 

         Deardorf, J. W., 1972: Parameterization of the planetary boundary layer for use in general circulation models. Mon. Wea. Rev., 100, 93-106.

 

         Holstlag A A., B. A. Boville, (1993) Local Versus Nonlocal Boundary-Layer Diffusion in a Global Climate Model, J. Climate, 6, 1825-1842.

 

         Konor, C. S., A. Arakawa, 2001: Incorporation of moist processes and a PBL parameterization into the generalized vertical coordinate model, Technical Report 102 of the Department of Atmospheric Sciences, University of California, Los Angeles, available at www.atmos.ucla.edu/~csk

 

         Lorenz E. N., 1960: Energy and Numerical Weather Prediction. Tellus, 12, 364-373.

 

         Louis J. F., 1979: A parametric model of vertical eddy fluxes in the atmosphere, B. Layer Met, 17, 187-202.

 

         Suarez, M.J., A. Arakawa and D. A. Randall, 1983: The parameterization of the planetary boundary layer in the UCLA general circulation model: formulation and results. Mon. Wea. Rev., 111, 2224-2243.