2 Dynamics Package

2.1 Overview

Three dynamical cores are currently suppported in the CAM code: the Eulerian-Spectral, the Semi-Lagrangian Spectral and the Lin-Rood Semi-Lagrangian Finite Volume. Each dynamical core is configured as a package that is independently testable using idealised forcings.

The dynamical core actually performs 2 nearly separate functions

• time stepping the dynamics
• advecting tracers

Although these functions are closely related, they are performed separately in most dynamical cores. A chemical transport model needs to replace the time stepping of the dynamics with reading of datasets, time interpolation, and a variety of consistency requirements. The advection of tracers is unchanged. The interfaces inside the core are specified so that only the dynamics could be replaced, or the dynamical core could be an offline advection module.

Dynamical cores provide a public type called state which holds (or returns upon query) the current values of the prognostic, derived, and boundary fields that are used to set the state communicated to the physical parameterizations. An interface routine converts between the dynamical state and the physical parameterization state. The state will hold information on distributed data decomposition which can be used for transpose and redistribution functions.

The variables needed for the physical parameterizations are $u, v, T, p, z, \omega, q$, and $C$. The dynamics state consists of the following prognostic, derived, and boundary variables for each of the dynamical cores.

Eulerian spectral	Semi-Lagrangian spectral	Lin-Rood
$\zeta, \delta, u,v, T, ln(p_s), q, \phi_s$	$u, v, \zeta, \delta, T, \nabla T, ln(p_s), q, \phi_s, \nabla \phi_s$	$u, v, \theta_v, \Delta p, q, \phi_s$

The dynamical core takes as input the current state (which must be consistent with the state history of the dynamics or be an initial or restart condition). On output, the dynamical core returns the updated state (and has updated the internal state history). Thus the time stepping method, leap frog explicit, semi-implicit, etc, are considered part of the dynamical core.

In addition, to initialization routines, the dynamical core must also interface with restart logic. Depending on the time stepping method, this may require saving state history and time integration information to the restart file.

2.2 Requirements

The functional requirements for the atmospheric package software are given in the Requirements Document. The most important for this section is that the dynamics be swapable (Requirement DYN04) in the sense that multiple dynamical cores work within the same interface and package structure. The physics and dynamics aspects are then represented in modular code. Interfaces with the history file package and the restart packages also comply with the package interface structures and requirments.

To support different dynamics formulations the CAM interface must support both a timesplit and a process split formulation. This implies that the physics package return tendencies for the dynamics. It also implies that the function of the coupling step between physics and dynamics is different for different dycores and physics packages.

Requirements are also specified for the support of parallelism and memory hierarchies so that the dynamics packages work effectively on target parallel architectures. These issues are discussed in the Architecture Document. The design presented satisfies these requirements and supports the approach of the Architecture Document.

2.3 Array Dimensions

The array dimensions are managed by the PMGRID module where some parameters are set at compile time. Many dimensions are determined at run time with arrays allocated based on runtime parameters such as the number of processors. To support a general parallel decomposition satisfying Requirements EUL02 and SLD02 as well as providing support for reduced grids in both the spectral dynamical cores, Requirement EUL03 and SLD03, the number of blocks for the dynamics are determined in DYN_GRID. However, static array dimensions are set in the PMGRID module. The parameters PLON, PLAT and PLEV are provided at compile time and determine the global grid resolution.

!-----------------------------------------------------------------------
!
! Purpose: Parameters and variables related to the dynamics grid
!
!-----------------------------------------------------------------------

   integer, parameter :: plon   = PLON                     ! number of longitude s
   integer, parameter :: plev   = PLEV                     ! number of vertical levels
   integer, parameter :: plat   = PLAT                     ! number of latitudes
   integer, parameter :: plevp  = plev + 1                 ! plev + 1
   integer, parameter :: nxpt   = 1                        ! no. of pts outside active domain of interpolant
   integer, parameter :: jintmx = 2                        ! number of extra latitudes in polar region
   integer, parameter :: i1     = 1 + nxpt                 ! model starting longitude index
   integer, parameter :: j1     = 1 + nxpt + jintmx        ! model starting latitude index
   integer, parameter :: plond  = plon + 1 + 2*nxpt        ! slt extended domain longitude
   integer, parameter :: plond1 = plond - i1 +1            ! slt extended domain longitude starting at i1
   integer, parameter :: platd  = plat + 2*nxpt + 2*jintmx ! slt extended domain lat.
   integer, parameter :: numbnd = nxpt + jintmx            ! no.of lats passed N and S of forecast lat
   integer, parameter :: plnlv  = plon*plev                ! Length of multilevel field slice
   integer, parameter :: plndlv = plond*plev               ! Length of multilevel 3-d field slice
!
id
   integer :: nblks      ! number of blocks owned by the processor
   integer :: begblk     ! beg. index for blocks owned by a given proc
   integer :: endblk     ! end. index for blocks owned by a given proc
!

2.4 Precision of Reals

The precision of reals is set by assigning r8 token the appropriate value in in shr_kind_mod. Each CAM module then uses the following statement

  use shr_kind_mod, only: r8 => shr_kind_r8      ! atmospheric model precision

Precision is then declared with the statement

 real(r8) coslat(plond) .

2.5 Dynamics Interface

The interface to the dynamics drivers occurs in subroutine CAM_run.

All dycore drivers use the call sequence

   call dyn_run(dyn_state, dyn_tend, dyn_buffer, dtime )

2.6 Derived Types

The derived data types are specified in the DYN_GRID module and the DYNAMICS_TYPES module. Basic decomposition information and the blocked grid data structures are in DYN_GRID where the derived type BLOCK is defined. The dynamics grid consists of multiple blocks per MPI process. The public state of the dynamics package is given in the DYNAMICS_TYPES module. Public state and source information are carried in types dyn_state and dyn_tend.

2.6.1 Dynamics Grid

The resolution parameters are provided through the PMGRID and RGRID modules. Precision is specified using the r8 token. The distributed grid is specified using 3-D grid blocks. It is assumed that these blocks do not overlap (except in a halo region). Since a variety of discretizations must be supported with this structure several special cases will be supported. For the three dycores presently in the model, the following indexing specification for blocks is sufficient.

   use pmgrid
   use rgrid
 
   type block
     integer  :: beglon,endlon         ! beginning(and ending) global lon index     
     integer  :: beglonex,endlonex     ! beginning(and ending) global lon index     
     integer  :: beglev,endlev         ! beginning(and ending) global level index     
     integer  :: beglevex,endlevex     ! beginning(and ending) global level index     
     integer  :: beglat,endlat         ! begining (and ending) global lat index
     integer  :: beglatex,endlatex     ! begining (and ending) global lat extensions
     integer  :: owner                 ! id of MPI process where block assigned
   end type block

   type dynamics_grid
     type (block), dimension(:), allocatable, private :: local blocks of the
                                       ! global computational grid
     integer, private :: nblocks       ! global block count
     integer          :: begblk        ! local block beginning index
     integer          :: endblk        ! local block ending index          
   end type dynamics_grid

The indices are all global referring to a logically rectangular lon-lat-lev base grid or the reduced grid indexing. For the SEAM grid, a different base grid is used with an indexing scheme appropriate for finite elements. We expect that a block will consist of a collection of elements. Points of the grid may be nodal points or quadrature points. If a 3-d decomposition is used it will be assumed that grid points in a block can be associated with a base surface grid and a public method will be provided that identifies the base grids global index. We need this in order to assemble atmospheric columns for the physics package. The global ordering of the base grid may be given in a linear array.

For adaptive grids or more general 3-D unstructured grids, eg. tetrahedral, the identification of points of a block with columns will pass through an intermediary interpolation.

The ESMF framework provides types for distributed grids and for physical grids that (roughly) correspond to the grids defined here. These types will be used, replacing the definitions here, when the framework is ready.

2.6.2 Dynamics Fields

Fields are defined inheriting properties of the block decomposition. They may be either surface or 3-D fields with appropriate block structure. The declaration and allocation of the public fields of the dynamics interface are handled in the DYNAMICS_TYPES module.

2.6.3 Dynamics Types

For the Eulerian and Semi-Lagrangian spectral dynamical cores the dynamics state is a bundle of fields of the following form

type dyn_state
    real(r8) :: ps   ! surface pressure
    real(r8) :: t3   ! temperature
    real(r8) :: u3   ! u-wind component
    real(r8) :: v3   ! v-wind component
    real(r8) :: q3   ! constituents
    real(r8) :: omga ! vertical velocity
    real(r8) :: phis ! Surface geopotential                                               
end type dyn_state

For the Lin-Rood dynamical core the dynamics state is a bundle of fields,

type dyn_state
    real(r8) :: ps           ! surface pressure
    real(r8) :: u3s          ! Staggered grid winds, latitude
    real(r8) :: v3s          ! Staggered grid winds, longitude
    real(r8) :: delp         ! delta pressure
    real(r8) :: pt           ! virtual potential temperature
    real(r8) :: q3           ! Moisture and constituents
    real(r8) :: omga         ! vertical velocity
    real(r8) :: phis         ! Surface geopotential
    real(r8) :: pe           ! edge pressure
    real(r8) :: pk           ! pe**cappa
    real(r8) :: piln         ! ln(pe)
    real(r8) :: pkz          ! finite-volume mean pk
    pointer  :: xy_dyn_state ! Pointer to dynamics state for 2D decomposition
end type dyn_state

The basic state information is needed for the CAM interface and the regridding through DP_COUPLING. Other dynamics fields, for example, diagnositic variables as well as prognostic fields such as the potential vorticity for the Eulerian and Semi-Lagrangian arrays are held in private variables of the dynamics buffer module.

For Lin-Rood dynamics pe, pk, piln, pkz, and the xy_dyn_state should be moved into the dynamics. They are currently used at the dp_coupling layer and hence appear here.

These dynamics types will become instances of the ESMF_Bundle with individual fields being typed using ESMF_Field.

2.6.4 Dynamics Source Terms

The Eulerian and Semi-Lagrangian Spectral dynamical cores are process split and require tendency terms. These come from the physics tendencies.

type dyn_tend
   real(r8) :: fu      ! u momentum source term
   real(r8) :: fv      ! v momentum source term
   real(r8) :: t2      ! energy equation source term
   real(r8) :: qminus  ! advected contituent source terms
   real(r8) :: qnats   ! non-advected contituent source terms
end type dyn_tend

The Lin-Rood finite volume dynamical core does not require source terms as states are updated from DP_COUPLING using the dynamics state update method. However, currently the following structure is used.

type dyn_tend
   real(r8) :: dudt        ! u momentum source term
   real(r8) :: dvdt        ! v momentum source term
   pointer  :: xy_dyn_tend ! Pointer to dynamics tend for 2D decomposition
end type dyn_tend

2.6.5 Dynamics Buffer

The dyn_buffer holds the time level history, the private copy of the state and the derived variables. This is required in the interface of the dynamics package with the atmospheric driver only to allow multiple instatiations of dynamics packages. It is only used by the dynamics and restart packages.

2.7 Public Methods

      dyn_grid_init       ! initialize block'ed data structure
      dyn_run             ! step the dynamics
      get_block_owner_d   ! get MPI process ID of block
      get_block_indices_d ! get local block index range
      get_block_xxx_d     ! get global base indices for a block
      where xxx is
       lat                ! for global latitude index
       lon                ! for global longitude index
       lev                ! for global level index
       unstruct           ! for global unstructured index

2.7.1 Initialization

The blocked data structure is initialized by calling the _init routines. This must be done before the fields of the dynamics state are initialized since they will use the same blocking information and indexing. (The _init routine calls GRDINI and GAUAW setting the coordinates of the horizontal Gaussian grid as well as the vertical levels.)

      dyn_grid_init      initialize block'ed data structure

2.7.2 Advective Transport Modules

To implement Requirement DYN05 and DYN06, the transport modules have been isolated and support either the hosting dynamical core advection or advection with a "data atmosphere model".

2.7.3 Unit Testing

The dynamics package can be tested in a stand alone mode by setting the adiabatic and or ideal_physics flag to true on the namelist input. This results in a Held-Suarez forcing from the physics package. Though physics is still compiled, the tendency terms provided to the dynamics package through DP_COUPLING are filled with other forcings and radiation, etc are not computed.

2.8 Internal Design

The internal design of data structures is highly dependent on the specifics of the dynamical core algorithms. The block structures are designed to provide a general way to interface with the internals that is also efficient. The decompositions of particular parts of the intenal methods may or may not utilize the public decompositions. For example, in the spectral dynamical cores, the code that calls the FFT's reverts to a latitude line decomposition no matter what the previous decomposition. The data transpositions required to support these reallocations for efficient computation is not constrained by this design.

The block structure may be configured to support efficient cache utilization in the dynamics. The block structure may also be used in a vector friendly fashion to provide performance through vector libraries, compilers and machine architectures. The cache based considerations suggest many small blocks while the vector considerations may show better performance with few large blocks. Since the indexing through the blocks is an internal matter of the dynamics code, specific optimizations can be achieved with out changing the interface specification.

The EUL and SLD dycores are organized to pass through the memory in a series of "scans". The strategy is inherited from low memory and out-of-core solvers and is advantageous on the non-uniform memory machines that constitute our target platforms. The CAM codes processed memory one (or two) latitude slice at a time, streaming the next latitude slice from slow memory. The CAM scans memory one block at a time for the grid point calculations. A transpose of the data arranges the relevant fields for transformation to spectral space and a parallel decomposition, by wave number, of the spectral coeffients. After processing in spectral space, the inverse transformation (synthesis) leaves the grid points in the block decomposition.