16 Calculations Performed in the Coupler

The Coupler performs crucial scientific calculations which, by design of the hub-and-spoke system, are not or can not be handled by the component models but are essential to the coupled integration. This section provides a description of these calculations.

16.1 Mapping

The Coupler is responsible for mapping (also called interpolation or regridding) data from one model's grid to another. The Coupler implements this mapping as a matrix-vector multiply which in case of mapping atmosphere data to the ocean grid would be:

$$\overbrace{\pmatrix{a_1 & a_2 & \ldots & a_n\cr}}^{n\ atmosp... ...pmatrix{o_1 & o_2 & \ldots & o_m\cr}}^{m\ ocean\ grid\ points}$$

The 1x$N$ matrix ${\bf A}$ contains all of the atmosphere grid points in a 2D horizontal plane unrolled into a single vector while the $M$x1 matrix ${\bf O}$ contains all the points in a 2D horizontal slice of the ocean grid. For a T42 atmosphere and x1 ocean grid, the matrix ${\bf W}$ would contain $(64 \times 128)= 8192$ rows and $(320 \times 384) = 122800$ columns! Luckily most of the elements of ${\bf W}$ are zero and this is really a sparse matrix multiply. The Coupler only stores the non-zero elements and their $(i,j)$ locations and performs the multiply with the corresponding $i^{th}$ element from the atmosphere grid and $j^{th}$ element from the ocean using the cpl_map_bun and the underlying MCT method, sMatAvMult.

The mapping weights in ${\bf W}$ are stored in files and pre-calculated using the Spherical Coordinate Remapping and Interpolation Package (SCRIP). See http://climate.lanl.gov/Software/SCRIP/. Two methods for calculating the weights are used in the Coupler. All state data is mapped with weights calculated using SCRIP's bilinear interpolation scheme while all fluxes are mapped with weights calculated using SCRIP's second-order conservative remapping scheme.

Given the assumptions of only 3 grids (Sec. 9.2), there are 5 sets of mapping weights read in by the Coupler: bilinear and conservative mappings for atmosphere to ocean grids, bilinear and conservative mappings for ocean to atmosphere and a conservative remapping for the river to the ocean grid (ocean to river is obviously not needed).

Note that the grid point numbering scheme mentioned in Section 10.2 is also present in the calculation of the mapping weights: the number of a model's grid point corresponds to the number of the row or column of that point in the mapping matrix ${\bf W}$.

16.2 Atmosphere/Ocean Surface Fluxes

The Coupler calculates the fluxes between the atmosphere and ocean for the following reasons: By convention, fluxes between two models with different resolutions are calculated on the grid with the higher resolution. The atmosphere model can not calculate the fluxes since that would require it to know the ocean's grid. The ocean only communicates with the Coupler once per day so to update the fluxes as often sub-diurnally, the fluxes are calculated in the Coupler. The Coupler receives ocean state data and holds it constant while calculating new fluxes with each receive of new atmosphere data. The new fluxes are mapped back to the atmosphere grid and sent to the atmosphere. The Coupler also keeps a running sum of the fluxes and sends the time average to the ocean. The atmosphere-ocean fluxes are calculated using the formulas below.

16.2.1 General Expressions

The fluxes across the interface are calculated from bulk formulae and general expressions are

$$% latex2html id marker 3874 \eqno (\thesection .1)$$

$$\begin{eqnarray*} \vec \tau & = & \rho_A~~u^{*2}~~~\Delta \vec U ~~ \ \vert \De... ...\ -~\epsilon ~\sigma ~~ T^4 ~~~ \ + ~~ \alpha^L ~~ L\downarrow, \end{eqnarray*}$$

where the turbulent velocity scales are given by

$$% latex2html id marker 3876 \eqno (\thesection .2)$$

$$\begin{eqnarray*} u^* & = & ~~CD^{1/2}~~~\vert \Delta \vec U \vert \\ Q^* & = &... ...= & \ ~~CH~~\vert \Delta \vec U \vert~~(\Delta \theta)~~u^{*-1}, \end{eqnarray*}$$

where $\rho_A$ is atmospheric surface density, $Cp_A$ is the specific heat, $\sigma = 5.67 \times 10^{-8}$W/m$^2$/K$^4$ is the Stefan-Boltzmann constant, $\epsilon$ is the emissivity of the interface, and $\alpha^L$ is the surface albedo for incident longwave radiation, $L\downarrow$. In (16.2) the differences $\Delta \vec U$, $\Delta q$ and $\Delta \theta$ are defined at each interface in accord with the convention of fluxes being positive down. The reflected downward incident longwave radiation is simply accounted for by assuming an emissivity, $\epsilon = 1$, and the water surface albedo for incident longwave radiation, $\alpha^L = 0.0$.

The transfer coefficients in (16.2), shifted to a height, $Z$, and considering the appropriate stability parameter, $\zeta$, are :

$$% latex2html id marker 3908 \eqno (\thesection .3)$$

$$\begin{eqnarray*} CD & = & \kappa ^2 \left[ ln \left( {Z \over Z^o} \right) - \p... ...\left[ ln \left( {Z \over Z^h} \right) - \psi _s \right] ^{-1}, \end{eqnarray*}$$

where $\kappa = 0.4$ is von Karman's constant and the integrated flux profiles, $\psi _m$ for momentum and $\psi _s$ for scalars, are functions of the stability parameter, $\zeta$. These functions as used in the coupler are:

$$\begin{eqnarray*} \psi_m(\zeta)~~ &=& ~~\psi_s(\zeta)~~=~~-5\zeta ~~~~~~~~~~~~~~... ...~~~~~~~~~~~~~~ \zeta~<~0 \\ \ X ~~ &=& ~~( 1 - 16 \zeta )^{1/4} \end{eqnarray*}$$

Above the atmospheric interfaces $i=1,~2$ and $3$ the stability parameter

$$\zeta~~=~~{\kappa~g~Z_A \over u^{*2} } \Big( {\theta^* \over \theta_v} \ + {Q^* \over (Z_v^{-1} + q_A) } \Big) ~~~,$$

where virtual potential temperature is computed as $\theta_v = \theta_A (1 + Z_v q_A)$, $q_A$ and $\theta_A$ are the lowest level atmospheric humidity, and potential temperature, respectively, and $Z_v = (\varrho(water) / \varrho(air)) - 1 = 0.606$.

In addition to surface fluxes, the atmospheric model requires effective surface albedos for both direct $\alpha (dir)$, and diffuse, $\alpha (dif)$, radiation at each wavelength. They are used in a single call to the computationally demanding atmospheric radiation routines. This call gives downward atmospheric albedo for diffuse radiation, $\alpha_a$($dif$). If direct and diffuse albedos, $\alpha^m$($dir$) and $\alpha^m$($dif$) for $m = 1, \, M$ different surfaces below an atmospheric grid cell are known, the multiple scattering coefficients, $R^m = \left( 1 - \alpha_a (dif) \, \alpha^m(dif) \right)^{-1}$ are computed, and with $f^m$ the fractional coverage of each surface type, then the albedos are given by

$$\begin{eqnarray*} \alpha (dif) &=& \left( 1 - F_2 \right) \big/ \left( 1 - F_2... ... (dir) &=& F_1 \left( 1 - \alpha_a (dif) \, \alpha (dif) \right) \end{eqnarray*}$$

$$F_1 = \sum\limits^M_{m=1} \, f^m \, R^m \, \alpha^m (dir)$$

$$F_2 = \sum\limits^M_{m=1} \, f^m \, R^m \, \left( 1 - \alpha^m (dif) \right)$$

Use of these albedos ensures that the solar radiation exiting the bottom of the atmosphere is identical to the sum of $M$ radiation calculations, each using an $\alpha^m (dif)$ and $\alpha^m (dir)$. At present, however, the albedo calculations are greatly simplified by assuming $\alpha_a(dif) = 0$, such that $R^m = 1$. This albedo then does not need to be passed from the atmosphere to the coupler, and the coupler simply computes the effective albedos to be passed to the atmosphere as

$$\begin{eqnarray*} \alpha (dif) &=& \sum\limits^M_{m=1} \, f^m \, \alpha^m (dif) \\ \alpha (dir) &=& \sum\limits^M_{m=1} \, f^m \, \alpha^m (dir) \end{eqnarray*}$$

In general the partition of radiation among the $M$ surfaces is a function of $R^m$, $\alpha^m (dif)$, and $\alpha^m (dir)$, and hence of wavelength. Hence, correct partitioning would need to be performed for each wavelength band within the radiation transfer code of the atmospheric model. This procedure is greatly simplified by partitioning the net solar radiation within the Coupler.

16.2.2 Specific Expressions

At the atmosphere-ocean interface the near-surface air is also assumed to be saturated,

$$q = 0.98 ~~ \rho_A^{-1} ~~ C_5 ~ \exp (C_6/T),$$

where the sea surface temperature is $T = SST$, $C_5 = 640380.0$ kg/m$^3$ and $C_6 = -5107.4$ K, and the factor 0.98 accounts for the salinity of the ocean. The differences (16.2) become

$$\begin{eqnarray*} \Delta \vec U &=& \vec U_A - \vec U \\ \Delta q &=& q_A - q \\ \Delta \theta &=& \theta_A - T, \end{eqnarray*}$$

where the surface current, $\vec U$, is presently assumed to be negligible.

The roughness length for momentum, $Z^o$ in meters, is a function of the atmospheric wind at 10 meters height, $U_{10}$:

$$Z^o = 10~ \exp - \left[ \kappa \left( {C_1 \over U_{10}} + C_2 + C_3 U_{10} \right) ^{-1/2} \right],$$

where $C_1 = 0.0027$ m/s, $C_2 = 0.000142,$ and $C_3 = 0.0000764$ m$^{-1}$s. The corresponding drag coefficient at 10m height and neutral stability is

$$C_{10}^N ~~=~~ C_1 ~ U_{10}^{-1} ~~+~~ C_2 ~~+~~ C_3 ~ U_{10} ~~.$$

The roughness length for heat, $Z^h$, is a function of stability, and for evaporation, $Z^e$, is a different constant:

$$\begin{eqnarray*} Z^h &=& 2.2 \times 10^{-9} m~~~~~\zeta > 0 \\ \ &=& 4.9 \times 10^{-5} m~~~~~\zeta \leq 0 \\ Z^e &=& 9.5 \times 10^{-5} m. \end{eqnarray*}$$

Since $\zeta$ is itself a function of the turbulent scales (16.2), and hence the fluxes, an iterative procedure is generally required to solve (16.1). First, $\zeta$ is set incrementally greater than zero when the air-sea temperature difference suggests stable stratification, otherwise it is set to zero. In either case, $\psi_m = \psi_s = 0$, and the initial transfer coefficients are then found from the roughness lengths at this $\zeta$ and $U_{10} = U_A$. As with sea-ice, these coefficients are used to approximate the initial flux scales (16.2) and the first iteration begins with an updated $\zeta$ and calculations of $\psi _m$ and $\psi _s$. The wind speed, $U_A$, is then shifted to its equivalent neutral value at 10m height :

$$U_{10} = U_A ~~\Big( 1 + { \sqrt{C_{10}^N} \over \kappa} ln( {Z_A \over 10} - \psi_m(\zeta)) ~ \Big)^{-1} ~~.$$

This wind speed is used to update the transfer coefficients and hence the flux scales. The second and final iteration begins with another update of $\zeta$. The final flux scales then give the fluxes calculated by (16.1).

16.3 Surface Albedo and Net Absorbed Solar Radiation

For the ice, land, and ocean components there are four surface albedos that are used by the atmosphere component to compute four corresponding components of downward shortwave. Subsequently, the four albedos, together with the four downward shortwave fields, are used to compute net absorbed shortwave flux from the atmosphere to the surface components:

$$SW_{net} = \sum\limits^4_{m=1} \, SW_{down}^m \, \left( 1 - \alpha^m \right) ~~,$$

where $m = 1,..., 4$ corresponds to near-infrared/diffuse, visible/diffuse, near-infrared/direct, and visible/direct shortwave components.

The ice and land components each compute their own surface albedos and, given the downward shortwave fields, compute their own net absorbed shortwave radiation. For the ocean component, it is the Coupler that computes the ocean surface albedo and computes net absorbed shortwave radiation. The Coupler then sends the net absorbed shortwave field to the ocean component.

16.3.1 Land Surface Albedos

The land surface albedos are computed by the land component and passed on to the atmosphere component. These albedos are not altered by the Coupler in any way. Note that the atmosphere component may be an active model computing downward shortwave once per hour, or it may be a data model feeding the Coupler a daily average downward shortwave. It is the user's responsibility to ensure that the albedos the land model sends are appropriate considering what type of downward shortwave fields the atmosphere component is providing.

16.3.2 Ice Surface Albedo

The ice surface albedos are computed by the ice component and passed through the Coupler and on to the atmosphere component. These albedos are ``60 degree reference albedos'' that have no diurnal cycle. Based on an input namelist variable, flx_albav (see the Section 4 of the User's Guide), the Coupler will either pass these albedos on to the atmosphere component unaltered (in which case they are considered daily average albedos), or impose a diurnal cycle on the albedos (in which case they are considered instantaneous albedos). When the Coupler does add a diurnal cycle to the ice albedo, this consists of merely setting the albedos to 1.0 on the dark side of the earth.

16.3.3 Ocean Surface Albedos

Unlike the ice and land components, the Coupler computes the ocean component's surface albedo. There are two ways the Coupler can compute the ocean albedo: with a diurnal cycle (instantaneous) or without a diurnal cycle (daily average). An input namelist variable (the flx_albav namelist variable) selects which option the Coupler implements.

If the albedos are computed as daily average albedos, then all four ocean albedos are set to 0.06 everywhere, regardless of time of day or time of year.

If the albedos are computed as instantaneous albedos, then all four ocean will be set to 1.0 on the dark side of the earth, and where the solar angle is greater than zero, the albedos are set to a value which has both an annual and a diurnal cycle. Ocean albedo distinguishes between direct and diffuse radiation. The direct albedo is solar zenith angle dependent, while the diffuse is not. There is no spectral dependence of the albedo, nor dependence on surface wind speed. The expressions for both direct and diffuse albedo are taken from Briegleb et al. (1986), based on fits to observations of ocean albedo, good to within $\pm 0.3\%$. The albedo expressions are valid for open ocean, and do not include the effects of suspended hydrosols in near-surface waters.

For complete details see Briegleb, B.P., P.Minnis, V.Ramanathan, and E.Harrison, 1986. ``Comparison of regional clear-sky albedos inferred from satellite observations and model comparisons.'' Journal of Climate and Applied Meteorology, Vol. 25, pp.214-226.

16.4 Area normalizing

Area normalizing is done in the Coupler to correct for slight differences between the total area of the sphere assumed by a model and the SCRIP program (Sec. 16.1). SCRIP calculates area weights for each grid as well as the mapping weights between two grids. However each model may have its own method for calculating the area of a grid cell. Thus when the conservative remapping is performed to interpolate a flux from one grid to another, it preserves the total flux over the sphere but the sphere may have a slightly different total area in the Coupler compared to each model. This will in turn effect global conservation of fluxed quantities. To correct for this effect, the Coupler multiplies all received fluxes by the ratio of the two areas immediately after receiving fluxes from a component:

$$Flux_{in\ coupler} = Flux_{from\ model} \frac{Area_{model}}{Area_{SCRIP}}$$

The corrected flux, $Flux_{in\ coupler}$ is then used within the coupler for mapping and any other calculations.

The Coupler receives $Area_{model}$ from each model during the CONTRACT initialization (Sec. 10.2) and stores it in a DOMAIN (Sec. 8.1) under the ``area'' attribute while $Area_{SCRIP}$ is read from the SCRIP mapping weight files during the MAP initialization and stored in the DOMAIN under the ``aream'' attribute. The area fractions are calculated using the areafact_init method from areafact_mod.F90. Before calculated or mapped fluxes are sent to a model from the Coupler, they are multiplied by $\frac{Area_{SCRIP}}{Area_{model}}$.

16.5 Merging and Fractional Weights

When two or more model's supply input to another model, the input field is formed by merging the two outputs. An example is the atmosphere where a atmosphere grid cell may overly both open ocean and sea ice covered ocean. In that case in CCSM3, the atmosphere-ice flux is calculated by the ice model while the atmosphere-ocean flux is calculated by the Coupler. Before sending to the atmosphere, these fluxes must be merged:

$$\ F_a = f_{ia} F_{i2a} + f_{la} F_{l2a} + f_{oa} F_{o2a}$$

where

$$\ f_{ia} + f_{la} + f_{oa} = 1.$$

$F_a$ is the total flux into the atmosphere, $f_{ia}$ is the fraction of ice on the atmosphere grid, an example of a surface fractional weight, and $F_{i2a}$ is an ice-to-atmosphere flux calculated by the ice model. The other terms or for the land-atmosphere and ocean-atmosphere fluxes. States, such as surface temperature, are also merged this way. Note that this merge is performed after all fluxes have been mapped to the atmosphere grid. In that case $f_{la}$ is always 1. Ocean fluxes must also be merged for ocean grid cells that are only partially ice covered. For example, the total momentum flux into the ocean is the weighted sum of the wind stress and the sea-ice drag.

The Coupler calculates all the time-invariant surface fractions using frac_init and updates the ice surface fractions to account for the time-varying extent of sea ice after each receive of data from the sea ice model using frac_set. Both methods are in frac_mod.F90.