F. Flux__Formulae___________
Flux computations in the Flux Coupler are denoted by a circle in Fig. 1. Many of the
fluxes across the interfaces i = 1; 2; 3 and 4 are calculated from bulk formulae and general
expressions are
"oi= aeA u*2i U"i |U"i |-1
Ei = aeA u*i Q*i
(F:1)
Hi = aeA CpA u*i *i
L "i = - ffli oe Ti4 + ffLi L #;
where the turbulent velocity scales are given by
u*i = CD1=2i |U"i |
Q*i = CEi |U"i | (qi) u*-1i (F:2)
*i = CHi |U"i | (i) u*-1i ;
where aeA is atmospheric surface density, CpA is the specific heat, oe = 5:67x10-8 W/m2 /K4
is the Stefan-Boltzmann constant, ffli is the emissivity of the interface and ffLi is the surface
albedo for incident longwave radiation, L #. In (F.2) the differences U"i , qi and i are
defined at each interface in accord with the convention of fluxes being positive down.
The transfer coefficients in (F.2), shifted to a height, Zi, (Zi = ZA for i = 1, 2, 3),
and to the stability, ii, appropriate to the particular interface, are :
-2
CDi = 2 ln Zi___Zo- m
i
-1 Z -1
CEi = 2 ln Zi___Zo- m ln _i__e - s
i Zi
-1 Z -1
CHi = 2 ln Zi___Zo- m ln _i___h - s ;
i Zi
where = 0:4 is von Karman's constant and the integrated flux profiles, m for momentum
and s for scalars, are functions of the stability parameter, ii. These functions as used in
the coupler are:
m (i) = s (i) = -5i i > 0
m (i) = 2ln[0:5(1 + X )] + ln[0:5(1 + X 2)] - 2tan-1 X + 0:5ss i < 0
s (i) = 2ln[0:5(1 + X 2)] i < 0
X = (1 - 16i)1=4
Above the atmospheric interfaces i = 1; 2 and 3 the stability parameter
i * Q* j
ii = _g_ZA______u*2i__+ ________i________-1 ;
iv (Zv + qA )
F-1
where virtual potential temperature is computed as v = A (1 + Zv qA ), qA and A are
the lowest level atmospheric humidity, and potential temperature, respectively, and Zv =
(%(water)=%(air)) - 1 = 0:606. Since ii is itself a function of the turbulent scales (F.2),
and hence the fluxes, an iterative procedure is generally required to solve (F.1).
In addition to surface fluxes, the atmospheric model requires effective surface albedos
for both direct ff (dir), and diffuse, ff (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, ffa (dif ). If direct and diffuse
albedos, ffm (dir) and ffm (dif ) for m = 1; M different surfaces below an atmospheric grid
cell are known, the multiple scattering coefficients, Rm = (1 - ffa (dif ) ffm (dif )) -1 are
computed, and with f m the fractional coverage of each surface type, then the albedos are
given by ffi
ff(dif ) = (1 - F2 ) ( 1 - F2 =; ffa (dif ))
ff(dir) = F1 ( 1 - ffa (dif ) ff(dif ))
XM
F1 = f m Rm ffm (dir)
m=1
XM
F2 = f m Rm (1 - ffm (dif ))
m=1
Use of these albedos ensures that the solar radiation exiting the bottom of the at-
mosphere is identical to the sum of M radiation calculations, each using an ffm (dif ) and
ffm (dir). At present, however, the albedo calculations are greatly simplified by assuming
ffa (dif ) = 0, such that Rm = 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
XM
ff(dif ) = f m ffm (dif )
m=1
XM
ff(dir) = f m ffm (dir)
m=1
In general the partition of radiation among the M surfaces is a function of Rm ,
ffm (dif ), and ffm (dir), and hence of wavelength. Hence, correct partitioning would need
to be performed for each wavelength band within the radiation transfer code of the atmo-
spheric model. This procedure is greatly simplified by partitioning the net solar radiation
within the Flux Coupler.
F-2
I. Fluxes Across the Atmosphere - Cryosphere Interface (1)
The state variables required to calculate the fluxes (F.1) for i = 1, are the atmospheric
wind velocity, potential temperature, density and specific humidity in (kg=kg) at a height
ZA , and the ice surface temperature T1 = TC . The ice surface humidity is assumed to be
saturated :
q1 = ae-1A C5 exp (C6 =T1 );
where C5 = 640380: kg/m3 and C6 = -5107:4 K. The differences in (F.2) become
U"1 = U"A - "U1
q1 = qA - q1
1 = A - T1 ;
where the ice velocity, U"1 , is presently assumed to be negligible.
Presently the coupler assumes constant and equal roughness lengths for ice :
Zo1 = Zh1 = Ze1 = 0:04m:
To start the iterative solver for the fluxes (F.1), i1 is set to zero : m = s = 0, and the
neutral coefficients are used to approximate the flux scales (F.2). The iteration begins by
computing i1 , finding the non-neutral transfer coefficients and updating the flux scales.
After two such iterations there is sufficient convergence to compute the fluxes (F.1). Also
in (F.1), the reflected downward incident longwave radiation is simply accounted for by
assuming an emissivity, ffl1 = 1, and the ice surface albedo for incident longwave radiation,
ffL1 = 0:0.
II. Fluxes Across the Atmosphere - Biosphere Interface (2)
The fluxes across the Atmosphere - Biosphere Interface (2) are computed in the bio-
sphere model. Their formulae can be written in the forms of (F.1) and are detailed in
Bonan (1996), "A Land Surface Model (LSM version 1.0) For Ecological, Hydrological,
and Atmospheric Studies : Technical Description and User's Guide".
III. Fluxes Across the Atmosphere - Hydrosphere Interface (3)
At the atmosphere-hydrosphere interface the near-surface air is also assumed to be
saturated,
q1 = 0:98 ae-1A C5 exp (C6 =T3 );
where the sea surface temperature is T3 = TH , and the factor 0.98 accounts for the salinity
of the ocean. The differences (F.2) become
U"3 = U"A - "U3
q3 = qA - q3
3 = A - T3 ;
F-3
where the surface current, U"3 , is presently assumed to be negligible.
The roughness length for momentum, Zo3 in meters, is a function of the atmospheric
wind at 10 meters height, U10 :
" -1 #
Zo3 = 10 exp - C1____U+ C2 + C3 U10 ;
10
where C1 = 0:0027 m/s, C2 = 0:000142; and C3 = 0:0000764 m-1 s. The corresponding
drag coefficient at 10m height and neutral stability is
CN10 = C1 U1-10 + C2 + C3 U10 :
The roughness length for heat, Zh3, is a function of stability, and for evaporation, Zo3, is a
different constant:
Zh3 = 2:2 x 10-9 m i3 > 0
= 4:9 x 10-5 m i3 0
Ze3 = 9:5 x 10-5 m :
Since the roughness lengths are neither constant nor equal, the iterative solution for
the fluxes (F.1) differs somewhat than that for sea-ice. First, i3 is set incrementally greater
than zero when the air-sea temperature difference suggests stable stratification, otherwise
it is set to zero. In either case, m = s = 0, and the initial transfer coefficients are
then found from the roughness lengths at this i3 and U10 = UA . As with sea-ice, these
coefficients are used to approximate the initial flux scales (F.2) and the first iteration
begins with an updated i3 and calculations of m and s . The wind speed, UA , is then
shifted to its equivalent neutral value at 10m height :
i p ______CN Z j -1
U10 = UA 1 + _________ln10( _A___10- m (i3 )) :
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 i3 . The final flux scales then give
the fluxes calculated by (F.1). As for the sea-ice interface (1), the reflected downward
incident longwave radiation is simply accounted for by assuming an emissivity, ffl3 = 1, and
the water surface albedo for incident longwave radiation, ffL3 = 0:0.
IV. Fluxes Across the Cryosphere - Hydrosphere Interface (4)
Only the stress between the water and ice can be computed using bulk formulae sim-
ilar to (F.1). The upper water column stability is assumed to be near neutral (i4 =
0:0; m (0) = 0:0), so that CD4 , and hence o"4 depend only on a constant roughness length
and the velocity difference,
Zo4 = 0:05m
U"4 = U"C - U"H ;
F-4
where U"C and UH" are, respectively, the ice velocity and surface current. The implied drag
coefficient for Z4 = 10m is CD4 = 0:0057.
However, with the present sea-ice model, the stress is not computed in this fashion,
but is given by a linear function of U"4 , as described in Section G.
The heat and freshwater fluxes between the cryosphere and hydrosphere due to the
freezing and melting of sea-ice are computed, respectively, in the hydrosphere (see ocean
model documentation) and cryosphere (see ice model documentation) component models.
Each column of the hydrospheric component needs to be scanned after each time step to see
if any ice is formed. The cumulative heat of fusion associated with the freezing is output
as QH . The sense of this heat transfer is downwards, so QH is positive. The heat and
freshwater fluxes associated with freezing are applied directly by the hydrospheric model
and not passed to the coupler other than by the cumulative QH > 0.
If there is no freezing,
QH = aeH___CpH____1___(Tf__-_TH__)________t < 0 ;
m
is used to represent the potential of the upper ocean layer (thickness 1 ) to melt sea-ice,
where tm is a melting time scale. Since the hydrosphere itself does not know how much
ice is available to be melted it does not respond to QH < 0 in any way. Instead the heat
and freshwater fluxes associated with melting, both due to warm surface water QH , and
atmospheric heating are computed within the cryosphere and passed to the Flux Coupler
as H4 and F4 , respectively.
The solar radiation passing through the bottom of the ice, S4 , is found in the cryosphere
component by reducing the atmospheric solar, S1 , to account for heat absorption by the
ice and snow.
V. Fluxes Across the Biosphere - Hydrosphere Interface (5)
The only important direct exchange between the biosphere and hydrosphere (across
interface 5) is river runoff. The runoff from each of the N land model grid points is passed
to the coupler as a vector Rn ; n = 1; N . The freshwater flux associated with runoff is
passed to the ocean model as a vector F5 (m); m = 1; M , where each element of F5 is the
runoff into one of the M ocean surface grid points. This flux is computed by a matrix
multiplication in the coupler
F5 = T R ;
where T is an (M by N ) runoff transport matrix that is specified in the coupler initialization
and remains constant in time. Matrix element Tm;n is the fraction of runoff from land point
n that is transported to ocean grid point m.
For an ocean grid of dimension Im longitudes and Jm latitudes, M = Im x Jm and
m = i + (j - 1)Im at grid point (i, j). Similarly, N = In x Jn and n = i + (j - 1)In at
grid point (i, j) of a land model with In longitudes and Jn latitudes. Note, runoff from a
F-5
land gridpoint can be transported to different ocean basins, including marginal seas where
it can be made to balance the local evaporation and precipitation. River runoff can be
distributed over a number of ocean grid points. Conservation of freshwater is satisfied by
the condition M
X
Tm;n = 1:0 ;
m=1
for all n.
At present runoff is not transported by the coupler, i.e., Tn;m = 0 for all n and m.
The freshwater contained in the active ocean is conserved in the coupler by multiplying
the precipitation by the uniform factor needed to balance the global average evaporation
every coupling period, see section G.
F-6