A6 Elevation Class Decomposition Design Document
The Design Document page provides a description of the algorithms, implementation and planned testing including unit, verification, validation and performance testing. Please read Step 1.3 Performance Expectations that explains feature documentation requirements from the performance group point of view.
Design Document
In the table below, 4.Equ means Equations and Algorithms, 5.Ver means Verification, 6.Perf - Performance, 7. Val - Validation, 
- completed, 
- in progress, 
- not done
Overview table for the owner and an approver of this feature
1.Description | Elevation Class Decomposition |
| 2.Owner | |
| 3.Created | |
| 4.Equ | |
| 5.Ver | |
| 6.Perf | |
| 7.Val | |
| 8.Approver | |
| 9.Approved Date | |
| V2.0 | Pending |
| V1.0 | Declined |
Title: Elevation Class Decomposition
Requirements and Design
ACME Atmosphere Group
Date: 2015-9-10
Summary
This design is for Phase 2 of Elevation Class development in ACME. Phase 1 only involves the land tasks described here, plus a lapse-rate downscaling of the atmosphere fields for the elevation classes in the land model.
This activity introduces elevation classes to the atmosphere and land physics modules in the ACME model. All of the atmospheric and land physics is applied to each of a modest set of elevation classes within each model grid cell. The number of elevation classes and their fractional area and mean elevation in each grid cell is read into the land and atmosphere models. The surface elevation is used to diagnose the vertical displacement of air parcels flowing through the grid cell, and conservation of energy and moisture is used to diagnose an orographic tendency for temperature and the tracers. Initially the land and atmosphere are assumed to use the same grid, which greatly simplifies coupling. The atmosphere and land history is written for each elevation class; this history can then be distributed according to a high-resolution surface elevation dataset, producing high resolution distributions of the simulated climate.
Implementation requires breaking the assumption that that atmospheric physics and dynamics are calculated on the same grid.
Requirements
Requirement 1: Apply all model physics to multiple subcolumns (elevation classes) within each grid cell
Date last modified: 2015-09-10
Contributors: Steve Ghan (Unlicensed) Steve Goldhaber (Unlicensed) Teklu Tesfa Ruby Leung
Requirement 2: Subgrid surface elevation influences atmospheric physics
Date last modified: 2015-09-10
Contributors: Steve Ghan (Unlicensed) Steve Goldhaber (Unlicensed) Teklu Tesfa Ruby Leung
Requirement 3: Atmosphere communicates with land and ocean through all elevation classes in coupler
Date last modified: 2015-09-10
Contributors: Steve Ghan (Unlicensed) Steve Goldhaber (Unlicensed) Teklu Tesfa Ruby Leung
In the simplest implementation the land and atmosphere share the same grid and elevation classes, so there is a one-to-one mapping between land and atmosphere columns. Coupler mods are required to handle cells that are partial land and partial ocean.
In the more general case of different land and atmosphere grids (such as when the land model uses watersheds), vertical interpolation between elevation classses is necessary. A 19-page document on this is available at /wiki/spaces/ATM/pages/86442139.
Some complications arise near ocean shore, when an elevation class spans both land and ocean.
Requirement 4: Write all model history for all elevation classes.
Date last modified: 2015-09-10
Contributors: Steve Ghan (Unlicensed) Steve Goldhaber (Unlicensed) Teklu Tesfa Ruby Leung
History for each elevation class is needed for postprocessing to distribute fields at high spatial resolution.
Each requirement is to be listed under a ”section” heading, as there will be a one-to-one correspondence between requirements, design, proposed implementation and testing. Requirements should not discuss technical software issues, but rather focus on model capability. To the extent possible, requirements should be relatively independent of each other, thus allowing a clean design solution, implementation and testing plan.
Algorithmic Formulations
Design solution 1:
To address Requirement 1, use Goldhaber's atmospheric physics-dynamics coupling generalization to allow multiple physics column for each dynamics column. Expand physics columns to span elevation classes, with all classes in a cell on the same PE.
Date last modified: 2017-03-21
Contributors: Steve Ghan (Unlicensed) Steve Goldhaber (Unlicensed) Teklu Tesfa Ruby Leung
Design solution 2:
To address Requirement 2, an orographic forcing tendency is diagnosed from the difference between the grid cell mean surface elevation z and the surface elevation zj of the elevation class j. First, the surface pressure pj of the elevation class is diagnosed from the hydrostatic relation. This will be used to determine the air mass mj of layers in each elevation class, as well as the pressures used for thermodynamics.
Then, the height rise (or descent) of air parcels traversing a grid cell is diagnosed h=min(hmax,z-zmean) where hmax=U/N , U is the grid cell mean wind speed at the model level, and N is the grid cell mean Brunt-Vaisalla frequency. Once the vertical displacement is known, the orographic profile X’ can be diagnosed, X'(z)=Xmean(z-h)
Note that, since T is not conserved by vertical displacement, potential temperature is used for the translation, which is then converted back to temperature. Furthermore, to ensure conservation of energy and tracer mass, the orographic profile is normalized by the area-weighted atmosphere mass. This profile is not applied directly to T and q, because that produces an excessively strong orographic profile. Instead, it is used to determine an orographic tendency: O=(X'-X)/tau where tau is an orography time scale. The tendency is applied to the conservation equations for T and q
Date last modified: 2015-09-10
Contributors: Steve Ghan (Unlicensed) Steve Goldhaber (Unlicensed) Teklu Tesfa Ruby Leung
Design solution 3: Read land-atmosphere mapping indices from preprocessor to address Requirement 3
Date last modified: 2017-03-21
Contributors: Steve Ghan (Unlicensed) Steve Goldhaber (Unlicensed) Teklu Tesfa Ruby Leung Robert Jacob
Coupler can work with elevation classes using enhanced attribute vectors. Need to know what coupler interface looks like. one long vector per processor. define vector indices in preprocessor. data layout done at runtime. uniform memory, but spatially variable computation. data layout in atmosphere and land must be done in a consistent way. required input: global grid id, topo unit id.
Design solution 4: Use Goldhaber's atmospheric physics-dynamics coupling generalization to address Requirement 4
Date last modified: 2017-03-21
Contributors: Steve Ghan (Unlicensed) Steve Goldhaber (Unlicensed) Teklu Tesfa Ruby Leung
Coupling between the dycore and the physics state will be used mostly unchanged.
At the end of d_p_coupling, a call will be made to a routine which derives tendencies from the newly-created grid-cell mean physics state and the grid-cell mean physics state saved from the last time step. These tendencies are used to update the elevation class states. This routine may also call the orographic tendency since it will have the correct grid-cell mean state information (saving a computation which would have to be done if called later in the physics timestep).
At the beginning of p_d_coupling, the physics grid-cell mean state will be saved as the previous timestep state. This will be used on the next timestep to compute elevation-class tendencies.
Design and Implementation
Coupler Implementation: Components of implementation of elevation classes in the coupler
Date last modified: 2017-05-18
Contributors: Robert JacobSteve Goldhaber (Unlicensed)
Because there is no SCRIP grid for elevation classes, coupling between the atmosphere and surface components will be via the standard SE grid (e.g., ne30np4). To pass elevation class information, elevation class state and flux data will be added to the CAM attribute vector. For each field which will be passed as elevation class data, the current scalar field (e.g., Sa_dens which is the air density at the lowest model state) will represent the ocean/ice portion of the grid cell (if any) and a new vector field will be added to the attribute vector to pass all land elevation class data (ordered from lowest to highest elevation). For instance, in the case of Sa_dens, there will be attribute vector entries, Sa_dens_ec01, Sa_dens_ec02, ... Sa_dens_ecxx where xx is MaxNoClass. A given AV must have the same number of attributes for every grid point so MaxNoClass attributes will be declared even if an SE grid point has fewer then MaxNoClass elevation classes.
Some terminology: We call the set of attribute vector fields (states and fluxes) for a non-elevation class run the 'original attribute set' or the 'ocean/ice attribute set'. There are MaxNoClass copies of the original attribute set which we call the elevation class attribute sets. This is just terminology as all attribute fields in all of these sets are members of the singleton attribute vector.
The elevation class file has separate elevation classes for grid cells that contain both ocean and sea-level land. If an atmosphere grid cell contains any ocean in its surface, then the first elevation class is an ocean elevation class. For fractional grid cells (part land, part ocean), there may also be a second sea-level elevation class which will be land (if there is sea-level land under that grid cell). A mask indicating whether the lowest elevation class is land or ocean will be added to the elevation class file if required.
The land will find its elevation class information via its initial dataset while the atmosphere will read the data in a special-purpose file (see above). In order for the land and atmosphere to verify that they are running a compatible simulation, the land and atmosphere will create a new attribute describing how many active elevation classes they have on each cell. The land will use the number of each grid cells elevation classes while the atmosphere will only send the number of land elevation classes (so grid cells entirely over ocean will send zero). The coupler will check that these match at each grid point during the initialization phase of the coupler and throw an error if they do not.
The coupler will no longer do the merge of all the surface → atmosphere fluxes. For land-atmosphere, this will happen on the atmosphere side by copying the right flux to the right flux in the right ec field.
If an elevation class is all-ocean, use the (merged ocean-sea-ice) - atmosphere fluxes. If an elevation class is all-land, use the land-atmosphere fluxes.
"if ec" conditions will be used in prep_atm_merge in prep_atm_mod.F90 to adjust the merge for the presence of elevation classes. Away from the coastline, ocean, sea ice fluxes will be merged in the coupler and passed to atmosphere as the "non-land sea level flux".
Atmosphere Implementation: Components of implementation of elevation classes in the atmosphere component
Date last modified: 2017-05-18
Contributors: Steve Goldhaber (Unlicensed) Steve Ghan (Unlicensed)
Storage of state for each elevation class
Since the atmospheric state variables in elevation classes are treated prognostically, the atmospheric state for each elevation class needs to be stored across time steps.
pbuf cannot be used for this, but the state can be stored as a module variable which is also written to the restart file
The grid cell mean state is also needed to diagnose the orographic tendency; this can be done with a separate stateg module variable
Elevation Class Initialization
Elevation classes will be initialized with an additional input file. This file specifies all the information needed to create elevation classes although it does not contain subgrid-scale initial data (i.e., the model is still initialized with standard grid-scale data). The specification for the file is:
Dimensions:
- MaxNoClass: The maximum number of elevation classes in any atmosphere physics column
- grid_size: The number of atmosphere physics columns
Variables:
- GridID(grid_size): The ID corresponding to the column ID in standard atmosphere input files
- NumOfSubgrid(grid_size): Number of elevation classes in each atmosphere physics column. NumOfSubgrid(i) <= MaxNoClass
- SubgridAreaFrac(MaxNoClass, grid_size): The fraction of the atmosphere column taken up by each elevation class
- AveSubgridElv(MaxNoClass, grid_size): The average elevation of each elevation class
Physics (Elevation Class) <==> Dynamics coupling
The atmospheric physics is influenced by the dynamics through the advective and adiabatic compression tendencies. For temperature, energy is conserved if those grid cell mean tendencies A are directly applied to the physics in each elevation class. For tracers greater care is required to prevent negative values. Following Ghan et al. (Climate Dynamics, 2002), negative values can be prevented without affecting the grid cell mean impact on physics by weighting A by the ratio of the tracer mixing ratio in each class / the grid cell mean tracer mixing ratio.
The dynamics is influenced by the physics through the grid cell mean heating and tracer tendencies. To conserve energy and tracer mass, the grid cell means are formed by weighting the terms in each elevation class by the fractional area and the layer mass for each class.
Physics (Elevation Class) <==> Surface coupling
When preparing the attribute vector for coupling to the surface, the atmosphere will export its state and fluxes from each of its land elevation classes to the corresponding attribute-vector field (e.g., atm EC1 cam_out(c)%rho ==> Sa_dens_ec01, atm EC2 cam_out(c)%rho ==> Sa_dens_ec02, etc. See the coupler section above). It will also copy the state and fluxes from its first elevation class to the original attribute vector field (e.g., atm EC1 cam_out(c)%rho ==> Sa_dens). This behavior is so that the ocean/ice will always find the correct information in the original attribute vector field (e.g., Sa_dens) and requires no modification for elevation class runs. Note that when the first elevation class is an ocean/ice EC, the EC attribute vector sets will be copied starting at atm EC 2 (e.g., atm EC2 cam_out(c)%rho ==> Sa_dens_ec01, atm EC3 cam_out(c)%rho ==> Sa_dens_ec02, etc.).
When processing the attribute vector from the coupler (after coupling), the atmosphere will only look at the elevation class attribute vector fields (e.g., Faxx_taux_ec01 ==> atm EC1 cam_in(c)%wsx, Faxx_taux_ec02 ==> atm EC2 cam_in(c)%wsx) for land elevation classes. For any elevation class which is over ocean, it will use the values from the original attribute vector fields (e.g., Faxx_taux ==> atm EC1 cam_in(c)%wsx). Also, if there is an ocean fraction, the land copying will skip atm EC1 (e.g., Faxx_taux_ec01 ==> atm EC2 cam_in(c)%wsx, Faxx_taux_ec02 ==> atm EC3 cam_in(c)%wsx).
Land Implementation: Components of implementation of elevation classes in the land component
Date last modified: 2017-05-18
Contributors: Steve Goldhaber (Unlicensed)
Land (Elevation Class) <==> Atmosphere coupling
When preparing the attribute vector for coupling to the atmosphere, the land will export its state and fluxes from each of its elevation classes to the corresponding attribute-vector field (e.g., lnd EC1 lnd2atm_vars%t_rad_grc ==> Sl_t_ec01, lnd EC2 lnd2atm_vars%t_rad_grc ==> Sl_t_ec02, etc. See the coupler section above).
When processing the attribute vector from the coupler (after coupling), the land will simply reverse the process outlined above to handle states and fluxes from the atmosphere (e.g., Sa_u_ec01 ==> land EC1 atm2lnd_vars%forc_u_grc, Sa_u_ec02 ==> land EC2 atm2lnd_vars%forc_u_grc).
Post-Processing ImplementationDate last modified: 2017-03-21
Contributors: Steve Goldhaber (Unlicensed)
Given the history written to each elevation class, postprocessors are needed to (a) form a gridcell mean for plotting at the scale of the grid cells, and (b) mapping to much higher resolution in regions with complex terrain. The gridcell mean is easily calculated using the fractional area of each elevation class. Mapping to high resolution is best done by vertically interpolating values in adjacent grid cells to the surface elevation of a high-resolution grid, and then horizontally interpolating those values from the grid cell locations to the location of the high-resolution grid. This can be done in a way that does not yield apparent discontinuities in the distributions. Grid-cell mean history output will also be available (using a separate addfld call).
Planned Verification and Unit Testing
Verification and Unit Testing:
Date last modified: 2017-02-09
Contributors: Steve Ghan (Unlicensed) Steve Goldhaber (Unlicensed)
- One elevation class per grid cell (roundoff changes to no elevation class run)
- Elevation classes in atmosphere only (smoke test / science test for validation)
- Use class numbers from EC file but all classes have standard CAM elevation (roundoff changes to no elevation class run)
- Identical elevation classes in land and atmosphere (science test for validation)
- restart bit-for-bit
Configuration
To set up experiments and tests for elevation classes, compsets will be defined for each. For each compset, the land will set an initial data set in its config_components.xml file while the atmosphere will set the elevation class data file in its config_components.xml file.
Planned Validation Testing
Validation Testing: Multiscale evaluation
Date last modified: 2017-2-09
Contributors: Steve Ghan (Unlicensed)
Validation will be performed with an AMIP simulation, using high-resolution observation datasets such as PRISM (2.5 minute), NOHRSC (2.5 minute), and Climate Research Unit (10 minute) for validation in regions with complex terrain. Key variables such as surface air temperature, precipitation rate, and snow water equivalent will be mapped to the high-resolution grid of the observations according to the surface elevation of the datasets, and then compared with the data at high resolution. Snow water validation in regions lacking high resolution datasets will be achieved by comparing summertime snow distributions with observed glacier distributions. See, e.g., Ghan, S. J., T. Shippert, and J. Fox, 2006: Physically-based global downscaling: Regional evaluation. J. Climate, 19, 429–445, http://dx.doi.org/10.1175/JCLI3622.1. Scatterplots showing bias vs elevation will be used to determine whether the orographic signature is too strong or too weak. Transects across mountain ranges will be used to evaluate resolution of rainshadows. Global validation will be performed by forming grid cell means of all history for comparison with global coarse resolution reanalyses and satellite data.
Planned Performance Testing
Performance Testing: Dependence on elevation classification and load balancing
Date last modified: 2017-2-09
Contributors: Steve Ghan (Unlicensed)
Performance expectations are laid out in
Ghan, S.J., X. Bian, A. G. Hunt, and A. Coleman, 2002: The thermodynamic influence of subgrid orography in a global climate model, Climate Dynamics, 20, 31-44, 10.1007/s00382-002-0257-5.
Ghan, S.J., and T. Shippert, 2005: Load balancing and scalability of a subgrid orography scheme in a global climate model. Int. J. High Performance Comput. Appl., 19, 237-245, doi: 10.1177/1094342005056112.
Since the number of elevation classes per grid cell is quite heterogeneous, we expect performance to depend on the domain decomposition of the atmospheric physics columns and land subunits. If the load can be distributed evenly, we expect that elevation classes will increase the cost of atmospheric and land physics by a factor of less than two, depending on the resolution of the grid (which explicitly resolves surface elevation) and on the elevation classification scheme. There is a cost to domain decomposition of the physics, which will require communication with the dynamics cells. The cost factor should go down as explicit resolution is refined.
For the ACME implementation we will perform atmosphere-land timing tests with several potential load balancing schemes, a variety of MPI tasks and levels of threading, and on several hardware architectures. Note that implementation requires all elevation classes in cell to be allocated to the same physics chunk, so that limits how small chunks can be.
This testing will yield optimal domain decompositions and processor counts for the targeted machines.