JAEA-Data/Code 2008-033

Improvement of a Land Surface Model for Accurate Prediction ofSurface Energy and Water Balances
Genki KATATA
Division of Environment and Radiation Sciences
Nuclear Science and Engineering Directorate
Japan Atomic Energy Agency
Tokai-mura, Naka-gun, Ibaraki-ken(Received December 19, 2008)
In order to predict energy and water balances between the biosphere and atmosphereaccurately, sophisticated schemes to calculate evaporation and adsorption processes in thesoil and cloud (fog) water deposition on vegetation were implemented in the
one-dimensional atmosphere-soil-vegetation model including CO2 exchange process(SOLVEG2). Performance tests in arid areas showed that the above schemes have asignificant effect on surface energy and water balances. The framework of the above
schemes incorporated in the SOLVEG2 and instruction for running the model are
documented. With further modifications of the model to implement the carbon exchanges
between the vegetation and soil, deposition processes of materials on the land surface,vegetation stress-growth-dynamics etc., the model is suited to evaluate an effect ofenvironmental loads to ecosystems by atmospheric pollutants and radioactive substances
under climate changes such as global warming and drought.

Keywords: Land Surface Model, Atmosphere, Soil, Vegetation, Cloud Water Deposition,
Evaporation, AdsorptioniJAEA-Data/Code 2008-033

地表面エネルギー・水収支の高精度な予測のための地表面モデルの高度化
日本原子力研究開発機構

原子力基礎工学研究部門
堅田

環境・放射線工学ユニット

元喜（2008 年 12 月 19 日

受理）

大気－陸面間の熱・水交換過程を高精度に評価するために、土壌中の蒸発・吸着過程や植生への雲（霧）水沈着過程などの新しいスキームを、CO2交換過程を含む大気－土壌－植生 1 次
元モデル(SOLVEG2)に導入した。著者は、これまでの研究でこのモデルを乾燥地域に適用し、これらの物理過程が地表面エネルギー・水収支に重要な影響を与えることを示した。本報告では、大気・土壌・植生サブモデルのそれぞれ導入された新しい物理スキームに関する物理方程式、および改良したモデルの利用方法について、詳細に記述した。今後、開発したモデルに土
壌－植生間の炭素交換、大気から陸面への物質沈着、植物のストレス・生長・動態などを組み込むことによって、地球温暖化や干ばつなどの気候変動の影響を含めた放射性物質等の環境負荷物質による生態系への影響を評価・予測することが可能となる。

原子力科学研究所（駐在）：〒319-1195

茨城県那珂郡東海村白方白根 2-4ii

JAEA-Data/Code 2008-033

Contents
1.

Introduction-------------------------------------------------------------------------------------------- 1

2.

Model overview --------------------------------------------------------------------------------------- 3

3.

Modifications of soil sub-model------------------------------------------------------------------ 4
3.1. Basic equations ------------------------------------------------------------------------------- 4
3.2. Soil thermal characteristics -------------------------------------------------------------- 5
3.3. Soil water characteristics ----------------------------------------------------------------- 6
3.4. Phase change processes of soil water -------------------------------------------------- 7
3.5. Surface runoff model ----------------------------------------------------------------------- 13

4.

Modifications of atmosphere and vegetation sub-models--------------------------------- 15
4.1. Basic equations for cloud liquid water ------------------------------------------------ 15
4.2. Capture efficiency of cloud water droplets ------------------------------------------- 15
4.3. Gravitational flux of cloud water ------------------------------------------------------- 16
4.4. Droplet size distribution of cloud water ---------------------------------------------- 17
4.5. Phase change processes for cloud water ---------------------------------------------- 19
4.6. Evaporation and condensation process on the leaf surface --------------------- 19

5.

Model code---------------------------------------------------------------------------------------------- 21
5.1. Structure of model code--------------------------------------------------------------------- 21
5.2. Compiling the model-------------------------------------------------------------------------26
5.3. Running the model--------------------------------------------------------------------------- 27

6.

Summary -----------------------------------------------------------------------------------------------33

References --------------------------------------------------------------------------------------------------- 35
Appendix: Description of model code----------------------------------------------------------------- 39iii

JAEA-Data/Code 2008-033

目

次

1.

はじめに------------------------------------------------------------------------------------------------- 1

2.

モデル概要---------------------------------------------------------------------------------------------- 3

3.

土壌サブモデルの改良------------------------------------------------------------------------------- 4

4.

5.

6.

3.1.

基礎方程式 ------------------------------------------------------------------------------------ 43.2.

土壌中熱移動に関するパラメータ-------------------------------------------------------- 53.3.

土壌中液体水移動に関するパラメータ ------------------------------------------------ 6

3.4.

土壌水の相変化過程------------------------------------------------------------------------- 7

3.5.

地表面流出モデル---------------------------------------------------------------------------- 13

大気・植生サブモデルの改良---------------------------------------------------------------------- 15
4.1.

雲（霧）水に関する基礎方程式 --------------------------------------------------------- 15

4.2.

植生の雲水捕集係数 ------------------------------------------------------------------------ 15

4.3.

雲水の重力沈降フラックス --------------------------------------------------------------- 16

4.4.

雲水の粒径分布 ------------------------------------------------------------------------------ 17

4.5.

雲水の相変化過程---------------------------------------------------------------------------- 19

4.6.

葉面水の蒸発・凝結過程 ------------------------------------------------------------------ 19

モデルコード------------------------------------------------------------------------------------------- 21
5.1.

モデルコードの構成------------------------------------------------------------------------- 21

5.2.

モデルのコンパイル------------------------------------------------------------------------- 26

5.3.

モデルの実行---------------------------------------------------------------------------------- 27

まとめ -------------------------------------------------------------------------------------------------- 33

参考文献 ----------------------------------------------------------------------------------------------------- 35
付録：モデルコードの詳細------------------------------------------------------------------------------- 39iv

JAEA-Data/Code 2008-033

1.

Introduction
The energy and water exchanges between the atmosphere and land surface areimportant for environmental studies such as the prediction of climate change, but are notwell understood. Especially, it is a challenging study to model and clarify these exchanges
among the atmosphere, soil, and vegetation under climate changes (such as global
warming) and human activities. In addition, the processes of heat and water exchanges at
the air-land surface strongly affect the material exchange of air pollutions, plant nutrients,
radionuclide substances, etc. to the land surface. Therefore, it is necessary to develop adetailed land surface model which is applicable to various environmental issues.
A detailed one-dimensional model for atmosphere-soil-vegetation interaction called
SOLVEG2 has been developed (Nagai 2004) 1) . The SOLVEG2 is unique in the way it treatsphase change process between liquid and water vapor in the soil and turbulent transfer ofcloud liquid water in the atmosphere (fog), momentum, heat, and water vapor. The newschemes of phase change processes in the soil and cloud water deposition on vegetation
have been lately developed (Katata et al. 2007 2) and 2008 3) ). These schemes are especiallysuitable for accurate prediction of surface energy and water balances. For example, dew
deposition and vapor adsorption have been shown by Agam (Ninari) and Belriner (2006) 4)
to be the most important water source in desert areas during the dry season. Furthermore,
dew at the soil surface has been shown to be an important source of moisture for plants,
biological soil crusts, insects, and small animals in desert areas (Jacobs et al., 1999) 5).
Moreover, latent heat fluxes due to dew and adsorption significantly affect the energybalance at the soil surface (Agam (Ninari) et al. 2004) 6) . For woody plants, fog deposition
has long been recognized to be as an important factor in determining water balance of
mountainous vegetation, especially in arid and semi-arid regions. Fog occurring in thesurface boundary layer is transported downward by turbulence generated by strong wind
shears at terrestrial surfaces such as plant canopies, and is captured by them. If the cloudwater captured on foliar and woody surfaces of plants does not exceed the storage capacityof the canopy during fog deposition, the water is either lost from the canopy to theatmosphere via surface evaporation, or it is directly absorbed by the leaves (Burgess andDawson 2004) 7). If the intercepted water increases above the storage capacity, the water on
the plant surface drips to the soil via throughfall and stemflow. This phenomenon is known
as fog precipitation (Hutley et al. 1997) 8). Fog precipitation can be considered to be a
crucial water resource and should be quantified for accurate prediction of water budget inarid areas. Considering the difficulty in collecting accurate meteorological or hydrological
data in arid environments on an annual or interannual basis, it is necessary to develop anovel land surface model which is applicable to the arid environment for better predictionsof surface energy and water balances underlying material exchanges between theatmosphere and biosphere.
－ 1 －JAEA-Data/Code 2008-033

The objective of this study is to summarize the framework of the new schemes
incorporated in the original SOLVEG2. Future prospects for an application of the modified
SOLVEG2 to environmental issues are also described in the paper.

－ 2 －JAEA-Data/Code 2008-033

2.

Model overview
The SOLVEG2 model consists of one-dimensional multilayer sub-models foratmosphere, soil, and vegetation with a radiation transfer scheme for calculating thetransmission of solar and long-wave radiation fluxes in the canopy layer. The variablesfrom the bottom of soil layer to the top of air layer were integrated numerically using animplicit finite difference method and Gaussian elimination method. A detailed descriptionof SOLVEG2 can be found in Nagai (2004) 1) . In the present study, the soil, atmosphere andvegetation components of SOLVEG2 were modified to model the processes of evaporation
and adsorption in the soil and cloud water deposition on vegetation more precisely. Basicequations for sub-models and the newly incorporated processes are described here. Detailsof modifications and model performance tests are described in Katata et al. (2007) 2) and(2008) 3) .

－ 3 －JAEA-Data/Code 2008-033

3.

Modifications of soil sub-model

3.1. Basic equations
The temporal change in soil temperature is expressed by the heat conduction equationasCs ρ s∂Ts∂ ⎛ ∂T ⎞= ⎜ λ s ⎟ − lE b ,∂t∂z ⎝ ∂z ⎠(3-1)

where t is the time [sec], z the depth of the soil [m], Ts the soil temperature [K], Cs and ρs the
specific heat [J kg-1 K-1 ] and the density of the bulk soil [kg m-3 ], respectively, λ the thermal
conductivity [Wm-1K-1 ], l the latent heat of vaporization [J kg -1 ], and Eb the phase changesof soil water [kg m-2 s-1 ], respectively.
The mass balance equation for liquid water is given as

ρw

∂η w
∂η w
∂ ⎛
⎞
= ⎜ Dw+ K ⎟ − Eb ,
∂t
∂z ⎝
∂z
⎠

(3-2)

where η w is the volumetric soil water content [m3 m -3 ], D w is the soil water diffusivity [m 2s-1 ], K the unsaturated hydraulic conductivity [m s -1 ], and ρ w the density of liquid water [kgm-3]. The soil water diffusivity D w is expressed by

Dw = K

∂ψ
,
∂η w

(3-3)

where ψ is the water potential [m].
Water vapor diffusion in the soil is considered in SOLVEG2. According to Fick's law,the diffusion equation of water vapor in the soil pores is expressed asρ∂q ⎞∂(η ws − η w )q s∂⎛
= ρε t ⎜ Dv (η ws − η w ) s ⎟ + Eb ,
∂t
∂z
∂z ⎠
⎝

(3-4)

where qs is the specific humidity in the soil pores [kg kg-1 ], Dv the diffusion coefficient ofwater vapor [m2 s -1], ε t the tortuosity, ρ the density of water vapor [kg m -3 ] and η ws thesaturated volumetric water content [m3 m -3 ]. The variable (η ws － η w) represents thevolumetric content of gaseous phase in the soil. The εt was chosen as two-thirds, as
recommended by Jackson et al. (1974) 9) . Convection of water vapor is neglected in
SOLVEG2 because its contribution to water vapor transport near the soil surface is small
under the natural condition that moderate heating and cooling processes are caused by
diurnal change of solar radiation (Grifoll et al. 2005) 10) . Although empirical relations forwater vapor enhancement factor have been proposed to match measurements withpredictions of Philip and DeVries (1957) 11) model, such modifications have not yieldedsatisfactory agreement with field data (Cahill and Parlange 1998) 12) . SOLVEG2 simulates
water vapor movement in the soil without introducing such empirical relations.
－ 4 －JAEA-Data/Code 2008-033

3.2. Soil thermal characteristics
The following formulae of λ and Csρs are used in the original SOLVEG2⎧max[1.0, 419 exp(− pF − 2.7 )]0.172⎩( pF ≤ 5.1),
( pF > 5.1)

λ=⎨C s ρ s = (C s ρ s )soil + η w (C s ρ s )water ,

(3-5)

(3-6)

where pF=log 10 (－ψ)+2. Although the above functions suggested by McCumber and Pielke(1981) 13) has been used in numerous land surface models, it is known that Eq. (3-5) tends
to overestimate (underestimate) λ which affects surface energy fluxes during wet (dry)periods (Peters-Lidard et al. 1998) 14) . Thus, in the modified version, Eq. (3-5) is replaced by
the more precise formulation of λ by McInnes (1981) 15)[]

λ = A + Bη w − ( A − D )exp − (Cη w )E ,(3-7)

where A, B, C, D, and E are constants derived from De Vries (1963) 16) given as

A = 0.65 − 0.78 ρ d + 0.6 ρ d2 ,B = 1.06 ρ d ,C = 1.0 + 2.6(mc / 100 )

− 0.5

(3-8)

,

D = 0.03 + 0.1ρ d2 ,E = 4.0,
where ρd is the dry bulk density of the soil [kg m-3], and mc the clay fraction [%]. Since the ρdfor most soils ranges between 1.1 and 1.6 kg m-3, the value of 1.45 kg m-3 for sandy loam
(Agam (Ninari) and Berliner 2004) 6) is applied to all soils. The formulation of Csρs validated
in Katata et al. (2007) 2) is also replaced as the following formulation by Brutsaert (1982) 17)

C s ρ s = (1.095 + 4.18η w ) × 10 6 .

(3-9)

Sample calculations of λ and Csρs using Eq. (3-5) and (3-6) together with Eq. (3-7) and(3-8) for a volcanic soil in Spain (‘Other’ soil in Table 3-1) is shown in Fig.3-1.

Fig.3-1 Comparisons of (a) soil thermal conductivity and (b) volumetric heat capacity
between original (red lines) and modified formulations (green lines).
－ 5 －JAEA-Data/Code 2008-033

3.3. Soil water characteristics
The soil water retention curve is essential for the simulation of liquid and water vapor
flow in an unsaturated zone. The original SOLVEG2 uses the commonly employed curve
suggested by Brooks and Corey (1964) 18) and expressed as−b⎛η ⎞ψ = ψ s ⎜⎜ w ⎟⎟ ,⎝ η ws ⎠and⎛η ⎞K = K s ⎜⎜ w ⎟⎟⎝ η ws ⎠(3-10)

2b + 3,(3-11)

where ψs and Ks are the saturated values of water potential [m] and K [m s-1], respectively,
and b is the empirical constant. These equations are widely used in various land surface
models because of their simplicity. The function of (3-10) has, however, significant
limitations at low soil water content. Since the Brooks and Corey's curve has been designed
for wet soil, it is generally used to predict ψ when η w is larger than the wilting point. The
curve, hence, cannot be applied to dry soil. Moreover, water vapor adsorption is controlled
by water retention curves in dry soil. The direct use of the curve for dry soil leads to errors
not only in the prediction of the actual soil water content but also in the evaluation ofadsorption processes in dry soil. A more applicable curve for dry soil is therefore needed.
To avoid this problem, the authoer adopted a modified soil water retention curve
suggested by Webb (2000) 19) consisting of two regions: the capillary region (van Genuchten1980) 20) and the adsorption region (Campbell and Shiozawa 1992) 21). The functions of bothcurves are described by Webb (2000) 19) as()⎧⎪α −1 S w−(1 / m ) − 1 (1 / n )ψ = ⎨ [γ ( S − S )+log ψ ]
wm
10 m
⎪⎩ 10

(η w ≥ η wm ) , S = (S − S ) / (1 − S ), ,rr
(η w < η wm ) w(3-12)

where S, Sr, S w, and S wm are the saturation ratio, residual and effective saturation ratio, andthe corresponding value of S w when η w = η wm , respectively, and α, n, and m (= 1－1/n) are
fitting parameters, and ψm water potential when η w = η wm [m]. The K can be described by
combining Mualem's pore-size distribution model (Mualem 1976) 22)

[ (

2

) ].

1/ m m. K = KsSw 1− 1− Sw(3-13)As an example, modified soil water retention curve and unsaturated hydraulicconductivity using Eq. (3-12) and (3-13) together with Eq. (3-10) and (3-11) for a volcanic
soil in Spain (‘Other’ soil in Table 3-1) are shown in Fig.3-2. As seen in the Fig.3-2 (a), thevalues of ψ calculated by Eq. (3-10) approach infinity with decreasing amounts of water inthe dry soil. In contrast, the modified curves can be applied to the regions from water
saturation to extreme drying. The summary of parameters for Eq. (3-12) and (3-13) is listed－ 6 －JAEA-Data/Code 2008-033

in Table 3-1.
(a)

(b)

Fig.3-2 Comparisons of (a) soil water retention curves and (b) unsaturated hydraulicconductivity between original (red lines) and modified formulations (green lines).
Table 3-1 List of parameters for modified soil sub-model (from Carsel and Parrish 1988) 23).
Soil type
1

Sand

ηw

α-1

(m3 m -3 )(m-1 )

0.430

0.069

κ
1.70×10

η wr

mc

Ks

(m3 m -3 )(%)(ms-1)

2.68

0.045

2.90

8.25×10-5

2.28

0.057

6.40

4.05×10-5

1.89

0.065

11.1

1.23×10-5

1.41

0.067

18.5

1.25×10-6

1.37

0.034

9.50

6.94×10-7

1.56

0.078

19.7

2.89×10-6

1.48

0.100

27.4

3.64×10-6

1.23

0.089

33.2

1.94×10-7

1.31

0.095

32.6

7.22×10-7

1.23

0.100

41.0

3.33×10-7

1.19

0.000

46.8

5.56×10-9n32Loamy sand

0.410

0.081

2.15×103

3

Sandy loam

0.410

0.133

1.42×103

4Silt loam

0.450

0.500

1.15×10
2

5

Silt

0.460

0.625

3.00×10
1

6

Loam

0.430

0.278

5.20×102

7

Sandy clay loam

0.390

0.169

1.00×103

8Silty clay loam

0.430

1.000

1.65×101

9Clay loam

0.410

0.526

9.80×101

10Sandy clay

0.380

0.370

6.00×101

11Silty clay

0.549

0.374

2.00×101－ 7 －JAEA-Data/Code 2008-033

12

Clay

0.446

6.667

3.50×10

1.17

0.000

65.0

9.49×10-9

1.56

0.078

19.7

2.89×10-6

0.00

0.000

0.00

0.00×100

1.04

0.150

5.00

1.25×10-5

1.50

0.198

19.7

4.00×10-50

13Organic

0.430

0.278

5.35×10
2

14

Water

1.000

0.000

0.00×10015Bedrock

0.250

1.234

3.68×10
2

16

Other

0.666

0.377

1.90×102

3.4. Phase change processes of soil water
The explicit calculation of phase change processes (evaporation, condensation andadsorption) in the soil (Eb) incorporated in the original SOLVEG2 are described by (3-1),(3-2), and (3-4). The formulation of Eb is expressed as

E b = ρ rb−1 [q sat (Ts ) − q s ](3-14)

where r b is the function of volumetric soil water content empirically determined by Kondo
et al. (1992) 24) as

rb = 0.02 F1 (η ws − η w ) F2 Dv−1

(3-15)

where F1 and F2 are the constants. The above formulations have, however, no theoreticalbackground which can be used to phase change processes in the soil. Thus, a new
formulation for E b is proposed and introduced to the modified SOLVEG2.
In the modified SOLVEG2, the soil is formed by aggregation of "cylindrical pores" eachof which has a different radius as shown in Fig.3-3a. We assume that only two patterns ofpores exist. One is filled with capillary water, and the other is filled with air and adsorbs
water on its wall. The drying mechanism of the soil in the modified SOLVEG2 is illustrated
by Fig.3-3b. When the soil is almost saturated, evaporation occurs at soil pores exposed tothe air at the ground surface (direct evaporation, E dir defined later in (3-30)). In this case,evaporation does not occur in the soil since almost all the pores are filled with capillarywater (Fig.3-3b, 1). When the soil dries out, pores with a large radius are dehydrated andadsorb water films on their walls. As a result, evaporation of water adsorbed by large pores(evaporation in the soil, Eb) contributes to the water vapor flux in addition to evaporationfrom small water-filled pores at the ground surface (Fig.3-3b, 2). When the soil is extremely
dry, all pores except for micro-pores are dehydrated and evaporation mainly occurs in thesoil instead of at the ground surface (Fig.3-3b, 3). In modified SOLVEG2, water vapor flux
at the air-land interface (i.e., total evaporation at the soil surface; z = 0) represents the sum－ 8 －JAEA-Data/Code 2008-033

of direct evaporation (Edir) and water vapor flux from the pores with adsorbed water film tothe atmosphere (E0 defined later in (3-34)), which is resulted in evaporation in the soil
throughout all soil layers. Similar models on the basis of cylindrical capillaries have been
widely used to describe hydraulic and thermodynamic characteristics in unsaturated
porous media (e.g., Mualem 1976) 22. In the present study, we apply this concept to describe
the phase change processes in the soil.

Fig.3-3 Schematic diagrams of (a) the cylindrical pore model and (b) the drying process ofsoil due to evaporation at the ground surface or in the soil
Considering the two patterns of pores (Fig.3-3a), a threshold radius of pores, r k [m],which determines whether a pore is filled with capillary water or not, is expressed asrk =2τ cos φ
ρ w gψ

(3-16)

where τ is the surface tension of the liquid water ( = 0.072) [N m-2], φ the contact angle ofwater ( = 0) [degree], and g the gravity acceleration [m s-2]. The largest pore radius in thefilled pores, r [m], is calculated by r = r k + t, where t is the film thickness of liquid water
adsorbed on walls. t is given by (Derjaguin et al. 1978 25) ; Iwamatsu and Horii 1996 26))

t=3

Asvl
6πρ w gΠ

(3-17)

where Asvl is the Hamaker constant [J] for solid-vapor interactions through an intervening－ 9 －JAEA-Data/Code 2008-033

liquid (for condensation Asvl < 0) and Π the disjoining pressure [m]. Or and Tuller (1999) 27)
determined Asvl = －6.0×10-20 [J] using soil data from Campbell and Shiozawa (1992) 21) . Ifinterfacial interaction between the adsorbed water and the pore's wall surface is only
induced by van der Waals forces, we can assume that the disjoining pressure Π is replacedby the water potential ψ (Tuller et al., 1999) 28) . It is assumed that only liquid film exists inthe pores whose radius is larger than rk. Using the soil water retention curve and the above
formulations, we obtain the pore-size distribution of the soil. The cylindrical pore model
detailed above has often been used to estimate pore size distribution of absorbents, e.g.,
activated carbon and silica gel.
The liquid-vapor interfacial area in the soil (Aia [m2 m -3]) is simultaneously calculatedwith the discrete pore-size distribution, and is described as⎞⎛ r⎜ ∑ Aia (r ) ≤ SA ⎟
⎟
⎜r
⎠
⎝ min

Aia (r ) = 2πκ (r − t )L(r )(3-18)

where κ is the fitting parameter, L(r) the pore length with r per unit volume [m m-3], r min theminimum radius in all of pores [m], and SA the specific surface area [m2 m -3 ]. The κ isdetermined from the fact that cumulative A ia does not exceed measured SA. The values of κfor various soil textures are given in Table 3-1. From data from several soils (Petersen et al.
1996 29) ; Campbell and Shiozawa 1992 21) ; Banin and Amiel 1970 30) ), we obtained thefollowing equation for SA()

SA = 0.06 mc2 + 2.01mc + 5.0 × 10 3 ,(3-19)

where mc the clay content [%]. The following function of log-normal pore radius distribution,f(r) = dη w / dr, is used in this study (Kosugi 1994) 31)⎡ [ln (r / rm )] ⎤η ws
f (r ) =
exp
⎢−
⎥
2ω 2
(2π )1 / 2 ωr ⎢⎣⎥⎦2(3-20)

where rm is the geometric mean radius [m], which is determined by ψrm [m] by (3-16), and ω
standard deviation [m]. ψrm and ω are obtained from van Genuchten's parameters of α andmψ rm = −α −1 (21 / m − 1)

1− m

and

,

(3-21)

ω 2 = (1 − m) ln[(21 / m − 1) / m],(3-22)

Kosugi (1994) 31) has reported that this model performs as well as any existing
empirical model for determining retention curves of various soils. Using (3-20), L(r) isexpressed from the surface area of a cylindrical pore－ 10 －

JAEA-Data/Code 2008-033

L(r ) =f (r ) f (r )Δr
,
≈
πr 2
πr 2

(3-23)

Discrete calculation of f(r) and L(r) in the regions divided into four hundred from water
saturation to extreme drying is carried out at each pore radius.
Assuming that the distribution of the relative humidity of air adjacent to adsorbed
water is similar in all pores, when the relative humidity of air adjacent to adsorbed wateris smaller than 35 %, almost all pores are filled with air and adsorb water according to theBrunauer-Emmett-Teller (BET) theory (Brunauer et al. 1938) 32) . Thus, we use the value ofκ when the adsorbed water content, η wads , is almost equal to η w. η wads is described asrη wads = η w∑ A (r )
ia

r min

SA

.

(3-24)

The capillary water content, η wmat , is also determined as η wmat (=η w－η wads ). The changesof η wads and η wmat with ψ are shown in Fig.3-4a.
(a)

(b)

Fig.3-4 Samples of calculated (a) adsorbed (blue line), capillary (green line), and total soil
water contents (red line), and (b) the cumulative air-water interfacial area in the soil (redline) and the product of mass transfer coefficient (k v) and liquid-vapor interfacial area (Aia )(green line).
To describe evaporation process in the soil, mass transfer model in the cylindrical tube
(Bird et al. 2001) 33) was applied to the modified SOLVEG2. Since we assumed that specifichumidity of air adjacent to adsorbed water is equally distributed in all pores, theevaporation rate of adsorbed water in a pore with r, Eb (r), is determined by

Eb = ρk v (r )Aia (r )(q wall − q s ) ,(3-25)

where kv (r) is the mass transfer coefficient [m s -1 ] for a pore with r, q wall , and qs the specifichumidity at the surface of water adsorbed on the wall and at the center of the cylindricalpore, respectively. kv is calculated as follows

kv =1.83Dv,r −t(3-26)－ 11 －

JAEA-Data/Code 2008-033

where a value of 1.83 represents the Sherwood number of developed laminar flow of acircular tube. It can be assumed that gaseous phase in the soil is quiescent or laminar if
any flow exists under small difference in pressure between the atmosphere and the soil
(Rolsten 1986) 34) , and liquid-vapor interfaces remains stable (Ranshoff and Radke 1988 35) ;Blunt and Scher 1995 36) ). Under such condition, mass transfer coefficient between the porewall and bulk phase is constant (i.e., Sherwood number = 3.66) in the cylindrical tube
regardless of bulk flow velocity if any; in other words, mass transfer coefficient is
proportional to inverse r. Thus, if pore size distribution in the soil is given, the modified
SOLVEG2 can provide the amount of evaporation in the soil, which is an integrated valueof Eb (r) throughout all pore radiuses. This is a new approach to describe detailed
movements of liquid and water vapor including phase change processes in the soil.
In the modified SOLVEG2, evaporation occurs in pores not filled with capillary water
whose radius is larger than rk as calculated by (3-16). Therefore, the total evaporation ratein the soil at a certain η w is represented by
rmax

rmax

rmin

rmin

(k v Aia )tot = ∫ k v (r )Aia (r )dr ≈ ∑ k v (r )Aia (r )Δr ,(3-27)

where r max is the maximum radius in all of the pores [m]. The change of cumulative Aia and(kv Aia ) tot for η w are shown in Fig.3-4b. Result reveals the general concept of decreasing (k v
Aia) tot with an increase in soil water. This indicates that an increase of interfacial areas witha decrease in pores filled with capillary water allows more evaporation conductance of soil
water; in other words, a thermodynamic equilibrium is formed between water vapor andliquid water when the soil dries out.
The q wall is assumed to be equal to the specific humidity of air at the surface of liquidwater adsorbed in the pores. When the soil water is in thermodynamic equilibrium, therelative humidity of the air adjacent to the adsorbed water in the pores, hae , is defined as(Israelachvili 1992) 37)⎛ ψg ⎞⎟⎟ ,hae = exp⎜⎜⎝ R w Ts ⎠(3-28)

where R w is the gas constant of water vapor [J mol -1 K-1 ]. q wall is thus expressed as

q wall = hae q sat (Ts ) ,(3-29)

Equation (3-28), combined with Eq. (3-16) forms the Kelvin equation, which is widelyused to calculate the vapor pressure in thermodynamic equilibrium within a
Representative Elementary Volume (REV). It also provides the relation of capillary radius(r) and saturated vapor pressure (q wall) on spherical meniscus of the same curvature radius.
Here, Kelvin's law is assumed to express vapor pressure of adsorbate on the cylindrical－ 12 －

JAEA-Data/Code 2008-033

pore wall since it can be considered that vapor pressure of the water film adsorbed on thecylindrical pore wall with the same curvature radius, where van der Waals force works, isin the same order as that on the spherical meniscus. In the modified SOLVEG2, the soil
water potential expressed by Eq. (3-12) which includes these effects is used to calculate ψ in
(3-28) instead of (3-16).
Vapor diffusion with coefficient (k v Aia ) tot in (3-27) between the surface of adsorbed
water and the air at the center of the cylindrical pore (evaporation or condensation in thesoil) is caused by (q wall－qs). The q wall calculated by (3-29) rapidly decreases with decreasing
soil water. When q wall finally becomes qs, evaporation (or condensation) in the soil by (3-25)
stop.
We use the following equation for direct evaporation at the ground surface from thepores filled with capillary water, Edir, which is regarded as evaporation

E dir = ρσc E 0 u [q sat (Ts 0 ) − q r ] ,(3-30)

where σ is the fractional area of pores filled with capillary water exposed to the air (≈ η wmat /η ws ), c E0 the bulk coefficient, u the horizontal wind speed, Ts0 the soil surface temperature[℃], and qr the specific humidity of the air [kg kg-1 ]. The σ becomes almost zero when η w issmall because there are almost no pores exposed to the air at the ground surface.
Considering the term Edir, the upper boundary conditions are expressed as− Cs ρ s λ∂Ts= H 0 + lE dir − Rnet ,∂z z =0(3-31)

for soil temperature,∂η⎛⎞= − Pr + E r + E dir ,− ⎜ Dw w + K ⎟∂z⎝⎠ z =0(3-32)

for soil liquid water content, and− ρεDv (η ws − η w )∂q s= − E dir + E 0 ,∂z z =0(3-33)

for water vapor in the soil, respectively, where Rnet and H0 are the net radiation and sensible
heat fluxes [W m -2 ], Pr the precipitation [kg m-2s -1 ], Er the surface runoff [kg m -2s -1],
respectively. Note that only the latent heat by direct evaporation Edir is included inEq.(3-31) since the latent heat caused by phase change processes in the soil, Eb, areexplicitly calculated at an each soil layer in the SOLVEG2. E0 [kg m -2s -1 ] is water vapor flux
from the pores with adsorbed water film to the atmosphere expressed by

E 0 = ρ (1 − σ )c E 0 u (q s 0 − q r )(3-34)－ 13 －

JAEA-Data/Code 2008-033

where qs0 the specific humidity at the soil surface [kg kg -1].
3.5. Surface runoff model
In the original SOLVEG2, total surface runoff can be calculated as an integration valueof surface runoff per unit of time, Er [kg m -2 s-1 ], and is excluded from the object of
SOLVEG2 calculation at each time step. In order to better represent the water budget atthe ground surface, a sophisticated surface runoff model is incorporated to the modified
model. It is assumed that Er only occurs when water storage at the ground surface at thecertain time (S w [kg m -2]) exceeds the water capacity of the ground surface (S wmax [kg m -2]);
that is, Er = S w－S wmax . S wmax depends on surface characteristics and is mathematicallydescribed (Driessen 1986) 38) as

S w max = 0.5 z r

sin 2 (σ − θ ) cot (σ + θ ) + cot (σ − θ )⋅sin σ
2 cos σ cos θ(3-35)

where zr is the surface roughness [mm], σ and θ the clod (or furrow) and slope angle of theland [degree], respectively. z r and σ are set to have typical values for untilled land as 15 mm
and 30°suggested by Driessen (1986) 38) . Using Er, S w at next time step (S w*) is determinedby the following equation of water balance at the ground surface for η w

S w* − S w= E pr 0 − E w − E r
dt(3-36)

where Epr0 is the precipitation rate at the ground surface [kg m-2 s -1], and E w the water
infiltration to the soil [kg m-2 s -1 ] determined by the simple water balance model (Schaake
et al. 1996) 39) used in the Noah Land Surface Model (Chen and Dudhia 2001) 40) ,respectively. Noted that evapo-transpiration is considered in water balance equation at the
ground surface for specific humidity in soil.

－ 14 －

JAEA-Data/Code 2008-033

4.

Modifications of atmosphere and vegetation sub-models

4.1. Basic equations for cloud liquid water
In the atmosphere sub-model, there is a one-dimensional diffusion equation for cloudliquid water content in the atmosphere, wf [kg kg-1]

∂w f
∂t

=

∂w f Fsed ⎞ aE cap + E f + E col − E pr
∂ ⎛⎜⎜ K z
⎟−
+
,
∂z ⎝
∂z
ρ ⎟⎠
ρ

(4-1)

where z is the height in the atmosphere [m], Kz the vertical turbulence diffusivity for w f
calculated by the turbulence closure model (Yamada 1981) 40) [m 2 s -1], Fsed the gravitational
flux of cloud water [kg m-2 s -1 ], a the leaf area density at each canopy layer [m -2 m -3 ], Ecap the
capture rate of cloud water by leaves [kg m -2 s -1 ], Ef evaporation or condensation rate ofcloud water [kg m-3 s -1 ], Ecol and Epr the capture rate of cloud water by rain droplets and theevaporation rate of rain droplets [kg m -3 s -1 ], respectively.
The vegetation sub-model calculates the leaf surface water for each canopy layer, w d[kg m -2 ], vertical liquid water flux in the canopy, Pr [kg m -2 s -1 ], and leaf temperature, Tc [K].
The wd is determined by the following equation∂wd= Eint − E d + E cap − Pd ,
∂t

(4.2)

where Eint is the water exchange resulting from the interception of precipitation by leaves[kg m -2 s-1 ], Ed the evaporation rate of leaf surface water [kg m -2 s-1 ], and Pd the drip fromleaves [kg m-2 s -1 ], respectively. Pr is calculated by the following equation∂Pr= a[E d + E s ] + E pr − E col .
∂z

(4.3)

The Tc is calculated by solving the leaf surface heat budget equation

R c = H c + l (E d + E s ) + H p ,(4.4)

where Rc, Hc , Es , and Hp are the net radiation, sensible heat flux [W m-2], transpiration
rate [kg m -2 s -1], and cooling by precipitation [W m-2], respectively. Each of these terms isdetermined by the leaf temperature and variables from the atmosphere and radiation
sub-models.
4.2. Capture efficiency of cloud water droplets
Cloud water deposition on the canopy, Ecap , is expressed by the following equations

E cap = εF f u ρw f ,
and

Ff =

1 − exp(− k p aΔz )
a Δz

(4.5)
,

(4.6)

where ε is the capture efficiency of leaves for cloud water, Ff the shielding coefficient for－ 15 －

JAEA-Data/Code 2008-033

cloud water in horizontal direction, |u| the wind speed [m s-1 ], k p the averaged projection
coefficient of an individual leaf for cloud water, and Δz the thickness of canopy layer atheight z [m], respectively. The average projection coefficient of radiation flux on an
individual leaf was introduced in the original SOLVEG2 (Nagai 2003) 42). In the present
study, we applied this coefficient of radiation flux to cloud water (kp ). This results in the factthat Ff incorporates the decrease of effective leaf surface area intercepting cloud water
below the total leaf area because of the inclination of leaf surface and the overlap of leaves.
The original SOLVEG2 is based on the assumption that all cloud water droplet trajectories
are perpendicular to the leaves and that the droplets are all captured, i.e., ε = 1. However, εis usually < 1 because the droplets moving towards the leaf surfaces travel along curved
streamlines that lead beyond the leaf, and only by loss of the inertia of the droplets, they
are intercepted by leaves, according to a theory of droplet impaction. In the present study, εis calculated with the following empirical function of the Stokes number Stk (Peters andEiden 1992) 43)a2

and⎛ a3 S tk ⎞⎟⎟ ,ε = ⎜⎜⎝ a3 S tk + a1 ⎠(4-7)

ρ w d p2 u,
9υd leaf(4-8)

S tk =

where a1 , a2 , and a3 are fitting parameters, dp the droplet diameter of cloud water [μm], νthe viscosity coefficient of air [kg m -1s-1 ], and dleaf the characteristic leaf length [m],respectively. A summary of the parameters used in (4-7) for several plant species is given in
Table 4-1.
Table 4-1 List of parameters for capture efficiency by Eq.(4-7).
Vegetation species

a1

a2

A3

dleaf

Reference(mm)
Needle leaf

5.0

1.05

1.0

1.0

Thorne et al. (1982) 44)

Broad leaf

0.5

1.90

5.0

Belot and Gauthier (1975) 45)

Grass

0.6

3.20

1.0

30.0
1.8

Crop

0.06

10.0

1.0

10.0

Aylor (1982) 47)

Others

5.0

1.05

1.0

1.0

Same as Needle leaf

Chamberlain (1967) 46)

4.3. Gravitational flux of cloud water
Total cloud water flux over the canopy is the sum of the turbulent flux, caused byatmospheric turbulent mixing, and the gravitational flux, Fsed in Eq.(4-1), which cannot be
neglected when cloud droplets are large or wind speed is low (Lovett 1984) 48) . According to－ 16 －

JAEA-Data/Code 2008-033

Stokes' law, Fsed for each cloud water droplet is given as

Fsed (d p ) = ρw f (d p )v s (d p ) ,andv s (d p ) =gd p (ρ w − ρ )
18μ

(4-9)

,

(4-10)

where vs is the sedimentation velocity of cloud water droplets [m s-1], g the gravitational
acceleration [m s -2 ], and μ the viscosity of air [m 2 s -1], respectively. Considering Fsed, Ecap foreach cloud water droplet calculated by Eq. (4-5) can be expressed as[]

E cap (d p ) = u ρw f (d p ) F f ε (d p ) + Fv v s (d p ) ,(4-11)

where Fv is the shielding coefficient for cloud water in the vertical direction, and isassumed to be equal to Ff.
4.4. Droplet size distribution of cloud water
Since Ecap varies with the droplet diameter of cloud water d p , the droplet size
distribution (DSD) is required for the calculation of ε for all cloud droplets. Since the DSD
is generally a site-specific parameter and varies greatly with studies, several DSD
parameterizations are introduced to the modified SOLVEG2. Best (1951) 49) suggested thefollowing empirical formulation
n −1

n

B
B
⎛ dp ⎞⎛ n ⎞⎛ d p ⎞⎟ ,g (d p ) = ⎜⎜ B ⎟⎟⎜⎜ ⎟⎟exp⎜⎜ −⎟
⎝ a B ⎠⎝ a B ⎠⎝ aB ⎠with[ (a B = max 1.0, 9.091 × 10 5 ρw f

)

0.559

],

(4-12)(4-13)

where g(d p ) is the probability density for cloud water droplets with d p, n B the empiricalconstant (nB = 3.27). Lovett (1984) 48) used the parameterization in his model to quantify
the amount of fog water deposition on the balsam fir forest.
The modified gamma distribution (Deirmendjian 1969) 50) is applied to describe the
DSD

⎛dp ⎞
g (d p ) = a D ⎜⎜ ⎟⎟⎝ 2 ⎠with bD =pDpD

q D (d mod / 2) DqqD⎡
⎛dp ⎞ ⎤
exp ⎢− bD ⎜⎜ ⎟⎟ ⎥ ,⎢⎣⎝ 2 ⎠ ⎥⎦

(4-14),(4-15)

where aD , pD , and qD are the constants (aD , pD , qD) = (3.041×10-4, 4.0, 1.77), and bD thefunction of dmod , which is the mode diameter of droplet [μm], respectively. The dmod isapproximately equal to mean droplet diameter dmean (= 17.03ρwf×10 3+9.72) [μm] which the－ 17 －

JAEA-Data/Code 2008-033

author derived from the data collected at the Waldstein forest in Germany (Katata et al.
2008) 3) .
Chaumerliac et al. (1987) 51) suggested the following expression using the dmean⎡⎛ d p ⎞⎤1⎟⎥ ,ln 2 ⎜⎜exp⎢−
2⎟2π σ D d p⎢⎣ 2σ D⎝ d mean ⎠⎥⎦g (d p ) =ND(4-16)

The above formulation has been used in the one-dimensional model for radiation fog
and low-level stratiform clouds developed by Bott and Trautmann (2002) 52) .
Klemm et al. (2005) 53) proposed the two log-normal distributions for eight classes ofliquid water content

[

]

[

]

⎡
⎧⎪ log10 (d p / 2) − e K 2 ⎫⎪⎤⎧⎪ log 10 (d p / 2) − bK 2 ⎫⎪g (d p ) = ⎢a K exp⎨−⎬⎥ LWC ,⎬ + d K exp⎨−
fK
cK
⎢⎣
⎪⎭⎥⎦
⎪⎩
⎪⎭
⎪⎩

(4-17)

where aK , b K , c K, dK , eK , fK, and LWC are the constants given as Table 4-2.
Table 4-2 List of parameters for droplet size distribution by Eq.(4-17).
bK

cK

eK

fK

0.001

0.798

0.415

0.176

0.006

0.809

0.304

0.823

0.186

0.003

0.837

0.514

0.050

0.857

0.186

0.003

0.889

0.529

0.5

0.044

0.893

0.167

0.018

0.917

0.312

0.6

0.064

0.926

0.201

0.006

2.911

2.122

0.7

0.064

0.951

0.183

0.010

1.079

0.521

0.8

0.039

0.996

0.339

0.027

1.013

0.154

LWC

aK

(g m -3)(g m -3)

0.1

0.008

0.722

0.167

0.2

0.021

0.769

0.3

0.039

0.4

dK
(g m -3)

Considering the DSD, total capture of cloud water by leaves is represented byd p max

E cap =∫ E (d )dd ,capp(4-18)pd p minwhere dpmax and dpmin are the maximum and minimum diameters among cloud water in allpresent cloud droplets [μm], set to 5 × dmean and 0 μm, respectively. In the modified
SOLVEG2 model, hundred bins from d pmin to d pmax , each with an increment of 1 μm, areintegrated.

－ 18 －

JAEA-Data/Code 2008-033

4.5. Phase change processes for cloud water
To describe evaporation and condensation processes of cloud water, the following masstransfer sub-model was introduced in the modified SOLVEG2. The evaporation orcondensation rate of cloud water droplets with d p, Ef (d p), is determined by

E f (d p ) = ρk f (d p )A f (d p )[q sat (Ta ) − q a ] ,(4-19)

where kf (d p ) is the mass transfer coefficient for cloud droplets [m s -1 ], Af (d p ) theliquid-vapor interfacial area per weight of cloud water droplets [m 2 m -3], and qsat and qa thesaturated and air specific humidity [kg kg-1 ], respectively. Assuming that all cloud water
droplets are of spherical shape, kf (d p) is calculated as follows

k f (d p ) = (2 Dv ) / d p ,(4-20)

where a value of 2 represents the Sherwood number of liquid water droplet at rest (Ranzand Marshall 1952) 54) , and Dv the diffusion coefficient of water vapor [m2 s -1 ], respectively.
Using the volume [4/3π(dp/2)3 ] and surface area of a cloud droplet [4π(dp /2) 2], Af (d p ) can be
described as⎤ 6 ρw f (d p )⎡ρw f (d p )2,
A f (d p ) = 4π (d p / 2) ⎢
⎥=
3
ρwd p
⎣⎢ ρ w (4 / 3)π (d p / 2) ⎦⎥(4-21)

where the second term on the middle of Eq. (4-21) is the number density of cloud water
droplets with dp [m -3]. Finally, total evaporation/condensation of cloud water Ef in Eq. (4.1)
expressed as Ef (dp ) is integrated from dpmin to d pmax , using Eqs. (4.20) and (4.21), as in thecase of the Ecap calculation.
4.6. Evaporation and condensation process on the leaf surface
The water vapor flux from leaf surface to the atmosphere is divided into two
components; transpiration Es and evaporation from leaf surface water Ed in Eq. (4-2),
expressed as

andEs = ρ

rd[q sat (Tc ) − q a ] ,R'(4-22)

Ed = ρ

rs[q sat (Tc ) − q a ] ,R'(4-23)

where R’ = (r a rs + ra r d + rs rd ), and ra , r d, and rs are evaporation resistances of the laminarleaf boundary layer, the leaf surface water, and the stomata of the leaves, respectively. r a
and the total of r d and r s are in series, but r d and rs are parallel because evaporation fromleaf surface water and transpiration only occur from leaf surface water and stomata,respectively. rs used in the model has been given by Nagai (2004) 1). ra is expressed using ageneric transfer coefficient that accounts for the influence of wind speed and object shape－ 19 －

JAEA-Data/Code 2008-033

and size estimated by mass transport theory (Magarey et al. 2005) 55) as

ra = (c u )

−0.5

and

,

(4-24)

c = c w (wd / wdw ) + cd (1 − wd / wdw ) ,(4-25)

where wdw is the leaf surface water amount where there is maximum evaporation(Deardorff 1978) 56) [kg m-3], and c w and cd the vapor exchange coefficients between leaf andcanopy air for wet and dry leaf surfaces, respectively [m 0.5 s-0.5]. The value of c varies fromc w to c d with a change of leaf wetness fraction (wd /wdw). The values of c w and cd were
determined from the experiments of the drying of water droplets on an artificial leaf
surface representing a typical hydrophobic leaf surface; the best estimates were 2.582 ×
10-2 and 1.123 × 10-1 [m0.5 s- 0.5 ], respectively (Magarey et al. 2005) 55) . The modified
SOLVEG2 uses the formulation of r d expressed as⎧(1 − x d )x d−1 rard = ⎨
0⎩(q sat (Tc ) > qa )2/3
with x d = min [1, (wd / wdw )] . (4-26)(q sat (Tc ) < q a )When water partially exists on the leaf, the balance of Es and Ed is controlled by total
resistance of r d and rs connected in parallel. When the leaf is completely dry, r d has thevalue of infinity because xd representing leaf wetness fraction is zero; this results in Ed = 0
by Eq. (4-23). In contrast, when leaf surface is covered with water sufficiently (wd > wdw), r d= 0 due to the fact that x d = 1. As a result, Es = 0 as calculated by Eq. (4-22).

－ 20 －

JAEA-Data/Code 2008-033

5.

Model code
The model code is written in fortran77 and 90, and executable on linux. The model has

a function of parallel calculation by MPI. Test calculations have been done on PC cluster.
Details of the model code and procedure to run the model are described here. Note that
underlined files represent newly incorporated routines in the present study.
5.1. Structure of model code
The SOLVEG2 model mainly consists of the three directories of source code (SRC),
input meteorological data (INPUT), and output directory (OUTPUT). Inside the SRC
directory, the include files and modules, program files, model setting files, parameter files,
and make files are found. Below are descriptions of the files and directories for a one
dimensional calculation. The modules for three-dimensional and coupling calculations aredescribed in the prior paper by Nagai (2004) 1) .
a)b)

Root directory: SOLVEG2/- NQS shell-script

go_1D.sh- Executablezsolveg_1D.exe- Parameter fileparam_1D

Source code directory: SOLVEG2/SRC/- Make shell-scripts

zmake_xxx- Makefiles

zmakefile_xxx••

Include and module files (Incl*, prm_*)- Inclcon1

Common blocks for constants and input meteorological data

- Inclcon2

Common blocks for constants defined in BLOCKDATA- Inclnum

Parameters for grid numbers- Inclnum_1D

Inclnum files for 1-D calculation- InclvariCommon blocks for variables- prm_soil

Module for the dry soil model- prm_fog

Module for fog deposition

Program files (*.f)- efalbedo.f

Subroutine EALBED: soil surface albedo- ehws.f

Subroutine EHWS: saturated soil water content- eli2va.f

Subroutine ELI2VA: specific humidity in soil pore- eppara.f

Subroutine EPPARA: leaf projection coefficient- espara.f

Subroutine ESPARA: soil heat capacity and conductivity－ 21 －

JAEA-Data/Code 2008-033- evpara.f

Subroutine EVPARA: soil vapor diffusivity/evaporation resistance- ewpara.f

Subroutine EWPARA: soil water conductivity and diffusivity- faipsy.f

Functions FAIM, FAIH, PSYM, PSYH, SHMD, and SHMDD:
soil surface exchange functions- fcpair.f

Function FCPAIR: specific heat of air- fcw.fFunction FCW: specific heat of water- feps.fFunction EPS: Capture efficiency of cloud water droplets
Function FDIST: Droplet size distribution (DSD) of cloud water
Function DGL10: Gauss-Legendre quadrature integration
Function DEIR: Deirmendjian parameterization of DSD- fdensa.f

Function FDENSA: air density- fl.fFunction FL: latent heat of vaporization- gtable.f

Subroutine GTABLE: soil parameters- gvprofile.f

Subroutine GVPROFILE: vegetation profile data

- gvtable.f

Subroutine GVTABLE: vegetation parameters- gzsolveg.f

Subroutine GZSOLVEG: soil and vegetation grid- gzsolveg.f_1D

gzsolveg.f files for 1-D calculations- main.f

Main routine SOLVEG2- main.f_1D

main.f for 1-D calculation- pblkd01.f

BLOCK DATA: atmospheric parameters- pdebugw0.f

Subroutine DEBUGW: atmospheric variable output- pfluxcal.f

Subroutines SFPR13, FLXCAL, KMHCAL, and SAVEOD:
turbulence and variable for the next time step

- pgener.f

Subroutines GENER, DIREC1, and DIREC2: diffusion scheme- pinit01.f

Subroutine MSHINT: atmosphere grid

- pinit01.f_1D

pinit01.f files for 1-D calculation- pinitpf.f

Subroutines INITPF and CLSL2A: initial atmospheric variables- pmain03.f

Subroutines UMAIN, TMAIN, EMAIN, and CMAIN: wind,temperature, specific humidity, fog water, turbulence, and CO2- ppcal.f

Subroutine PCAL: air pressure- ppread.f

Subroutines PREAD and DEWTMP: parameters and initial data

- ppread.f_1D

ppread.f files for 1-D calculations- ptint.f

Subroutine TIMEINT: boundary data

- ptint.f_1Dptint.f files for 1-D calculations- shifi1.f

Subroutine HIFI1: advection in soil

- slco2.f

Subroutine SLCO2: soil CO2- sliqu.f

Subroutine SLIQU: soil water- solveg.f

Subroutines SLVGIN and SOLVEG: soil and vegetation control－ 22 －

JAEA-Data/Code 2008-033- solver1.f

Subroutine SOLV1: diffusion scheme- solver2.f

Subroutine SOLV2: diffusion scheme- srad.f

Subroutine SFCRAD: canopy radiation transmission- sradiatn.f

Subroutine RADIATION: solar and long-wave radiation- stemp.f

Subroutine STEMP: soil temperature- svapo.f

Subroutine SVAPO: specific humidity in soil pore- svfogcp.f

Subroutine FOGCAP: cloud water collection rate by leaves- svliqu.f

Subroutine VLIQU: leaf surface water and canopy water flux

- svrsco2.f

Subroutine RSCO2: CO2 assimilation and stomatal resistance- svrsst.f

Subroutine RESISTS: stomatal resistance- svsed.f

Subroutine SED: gravitational flux of cloud water- svtemp.f

Subroutine VTEMP: vegetation temperature- swadsp.f

Subroutine WADSP0: water retention curve and phase change
processes in the soil for the dry condition
Subroutine LINEINT: Linear interpolation of data
Function POTEV: potential surface evaporation•c)

Parameter files- Vertical grid files

zmesh.model_1D- Soil parameter file for Brooks and Corey

BCsoil.table (previously zsoil.table)- Soil parameter file for the dry soil modelvGsoil.table- Vegetation profile filezvege.profile_1D- Vegetation parameter filezvege.table_1D

Input meteorological data directory: SOLVEG2/INPUT/- Meteorology data filemetdata.dat (file format is shown in Table 5-1)

Table 5-1 Meteorology data file for 1-D calculation: metdata.dat.
TIME

P

RS

RL

RR

U

V

T

Q

WL

CO2

1997-04-20_0015 964.36 108.28 317.14 0.00

-1.79

2.86

293.016.93

0.12

370.001997-04-20_0045 964.53

-1.75

2.41

291.966.81

0.32

370.0024.64 313.94 0.00

… repeat until the end time of the calculation period
P: surface pressure [hPa], RS: solar radiation flux [W m-2], RL: long-wave radiation flux [W m-2], RR:
rain intensity [mm h -1 ], U: wind u-component [m s-1], V: wind v-component [m s -1], T: air temperature[K], Q: specific humidity [g kg-1 ], WL: cloud liquid water content [g m-3 ], CO2: CO 2 concentration[ppmv].

－ 23 －

JAEA-Data/Code 2008-033d)

Output data directory: SOLVEG2/OUTPUT/- Input meteorological data file (fu23)

METout- Surface flux file (fu24)

FLXout- Wind speed file (fu25)

WNDout- Atmospheric CO2 file (fu26)

ACO2out- Atmospheric CO2 budget file (fu27)

BACO2out- Soil variable file (fu30)

mnout- Soil temperature file (fu31)

TSout- Soil water content file (fu32)

HWout- Soil humidity file (fu33)

QSout- Soil evaporation file (fu34)

EBout- Soil surface flux file (fu35)

SFout- Canopy variable fale (fu40)

VGout- Canopy water budget file (fu41)

VWout- Canopy heat budget file (fu42)

VTout- Canopy radiation file (fu43)

RADout- Soil CO2 file (fu50)

SCO2out- Canopy CO2 file (fu51)

VCO2out- Accumulated fog and precipitation (fu52)

PREout- Soil CO2 production file (fu52)

PSCO2out- Soil CO2 budget file (fu53)

BSCO2out- Atmospheric variable file (fu20)

dbout- Standard output files (fu06)

outlist- Water retention curve (fu36: only for IFDSL = 1)

wcurve－ 24 －

JAEA-Data/Code 2008-033****************************************************************************************Main| Primary| Secondary| Functions and libraries---------------------------------------------------------------------------------------SOLVEG2*├--- PREAD├--- MSHINT├--- SLVGIN* ┬-------------------------------------- SHMD│├------ GZSOLVEG│├------ GTABLE│├------ GVTABLE│├------ WADSP0│├------ LINEINT│└------ GVPROFILE├--- INITPF ┬------ PCAL│└------ CLSL2A├--- FLXCAL├--- SFPR13 ---------------------------------------- PSYM, PSYH├--- DEBUGW(Start of time loop)┌-┐↑ ├-- TIMEINT ┬│├------ DEWTMP ------------------------- SHMD│└------ PCAL├-- KMHCAL├-- FLXCAL├-- SOLVEG* ┬-------------------------------------- PSYM, PSYH, SHMD, FL, FCPAIR│├------ DEWTMP ------------------------- SHMD│├------ GVPROFILE│├------ SFCRAD ┬------ RADIATION││└------ EPPARA│├------ SVAPO* ┬---------------------- FDENSA, SHMD, FL││├------ EVPARA││├------ EHWS││├------ LINEINT││└------ SOLV1│├------ SLIQU* ┬------ EWPARA*---LINEINT││└------ SOLV2│├------ STEMP* ┬---------------------- FCPAIR, FCW, FL, POTEV││├------ ESPARA*--------- FCW││├------ EALBED││└------ SOLV1│├------ SLCO2┬------ EHWS││├------ HIFI1││└------ SOLVE2│├------ RSCO2------------------------ SHMD│├------ RESISTS ------------------------ SHMD│├------ VLIQU* ┬---------------------- SHMD││├------ EPPARA││└------FOGCAP---------- EPS, FDIST, DGL10, DEIR, VGRV│└------ VTEMP┬---------------------- FCPAIR, FCW, FL, SHMD, SHMDD│└------ EPPARA├-- SFPR13 ---------------------------------------- PSYM, PSYH├-- UMAIN----------------┬---------------------- PSYM│├------ GENER│├------ DIREC2

Continued on the following page.

－ 25 －

JAEA-Data/Code 2008-033│└------ SED------------- FDIST, VGRV├-- TMAIN* ----------------┬---------------------- PSYM, PSYH, SHMD, FCPAIR, FL│├------ GENER│├------ DIREC2│└------ ELI2VA -------- SHMD, SHMDD├-- EMAIN----------------┬---------------------- FAIM│├------ GENER│└------ DIREC1├-- CMAIN----------------┬---------------------- PSYM, PSYH│├------ GENER│└------ DIREC2├-- SAVEOD├-- DEBUGW↑ │└-┘(End of time loop)
STOP
END****************************************************************************************

Fig. 5-1 Calculation flow of the modified SOLVEG2. Underlined files are newly
incorporated routines to the original SOLVEG2. Single asterisks are added to the routines
modified in the present study. Note that the routines for MPI libraries and coupling
modules for the Numerical Environmental System SPEEDI-MP are omitted in the flow.
5.2. Compiling the model
Before the model users go on to compile and run the model, the Inclnum_1D file at the
SOLVEG2/SRC directory should be set to correct values of vertical grids for atmosphere(N1 and M1), vegetation (NC) and soil layers (NS). An example of this file is shown in
Tables 5-2. Note that the number of vegetation layers are smaller than those of
atmosphere; i.e., NC < N1. The numbers of total horizontal grids NX (x-direction) and NY(y-direction), and the grids for output file IX and JY are set to unity because of aone-dimensional calculation. The parameter NA should be set to the layer number where
atmospheric variables are generated at the directory of SOLVEG2/OUTPUT. Under above
settings, a compilation of the model is started by running the script of zmake_1D.sh at the
SOLVEG2/SRC directory. After the execution of this script normally terminates, theexecution file of SOLVEG2 named zsolveg_1D.exe must be created at the root directory
SOLVEG2/.
Table 5-2 Model grid number setting file: Inclnum_1D.
IMPLICIT DOUBLEPRECISION ( A-H , O-Z )PARAMETER ( N1 =

9 , M1 = N1+1 , N2 = 2*M1 )PARAMETER ( NS =7 , NC =9 )PARAMETER ( NX =1 , NY =1 )PARAMETER ( IX =1 , JY =1 , NA =－ 26 －

7 )JAEA-Data/Code 2008-033

5.3. Running the model
Following the setting of grid number in the file of Inclnum_1D (Table 5-2), the verticalmesh sizes are edit in the zmesh.model_1D file (Table 5-3). In addition to this, an input fileof hourly (or half hourly) meteorological data covering throughout the calculation period isnecessary before the model execution (Table 5-1). In the source directory of SOLVEG2/SRC,
the vegetation profile file (zvege.profile_1D: Table 5-4) needs to be modified to specify
variations in the whole calculation period. In this file, vertical distributions of vegetation
type (VTYPE), leaf area density (AZ), and root fraction for each vegetation type need to be
specified.

The

vegetation

type

is

chosenfromthevegetation

parameter

file(zvege.table_1D), which is specified by two integers: the former represents the category ofvegetation, and the latter the spatial or temporal variation of them.
Table 5-3 Model vertical grid setting file: zmesh.model_1D
SOLVEG Z-MESH DATA: 2002/02/12 CREATED BY NAGAI
SOIL LAYER-0.000

ZS(K) K=1,NS+1 : NS=7-0.020-0.050-0.100-0.200-0.500ATMOSPHERE LAYERS (M1) AND CANOPY LAYERS (NC)
0.1

0.3

0.5

0.7

1.0

-1.000-2.000Z(K) K=1,M1
1.5

3.0

5.0

8.0

12.0

Table 5-4 Vegetation profile file: zvege.profile_1D.
***** DATA FORMATT *************************************************************

yymmdd hrmnsc : DATE AND TIME -----------------------------------+ repeat**

N : NUMBER OF LANDUSE CATEGORIES

|

*

*

N1 / LANDUSE (COMMENTS)------------------------+ repeat N times |**VTYPE(K),K=1,NC : NO in zvege.tabele

|

|

*

*

AZ(K),K=1,NC: leaf area density

|

|

*

*

L : number of VTYPE (ROOT FRACTION)

|

|

*

*

V1 RZ(K,V1),K=1,NS -------+ repeat L times

|

|

*

|

|

*

*

|

|

*

*

---END OF THIS LANDUSE CATEGORY ----------------+

|

*

*

N2 / LANDUSE (COMMENTS)

|

*

|

*

*

|

*

*

***END OF THIS TIME ---------------------------------------------+**

yymmdd hrmnsc : DATE AND TIME

*

*

*

.....

.....

*
.....
Continued on the following page.

*－ 27 －

JAEA-Data/Code 2008-033******* DATA FORMATT ************************************************************!!!!! DATA START !!!!!
030101 000000 : DATE AND TIME
1 : NUMBER OF LANDUSE CATEGORIES
1 / Canary Islands
31

31

31

31

31

31

31

0.050 0.050 0.100 1.000 1.600 1.200 0.630
1 : NUMBER OF ROOT (ROOT FRACTION)
31

0.030

0.045

0.075

0.09250.2775

0.330

0.150

0.330

0.150

---END OF THIS LANDUSE CATEGORY***END OF THIS TIME
030201 000000 : DATE AND TIME
1 : NUMBER OF LANDUSE CATEGORIES
1 / Canary Islands
31

31

31

31

31

31

31

0.050 0.050 0.150 1.200 1.800 2.200 0.800
1 : NUMBER OF ROOT (ROOT FRACTION)
31

0.030

0.045

0.075

0.09250.2775---END OF THIS LANDUSE CATEGORY***END OF THIS TIME

…repeat until the end time of the calculation period
In the parameter file (param_1D), the simulation condition such as calculation period,
output interval is specified (Table 5-5). Parameters of DX, X00, and Y00 are meaninglessfor one-dimensional calculations.
Table 5-5 Sample parameter file: param_1D.
IPRINT=

1800<=========( S )PRINT OUT INTERVAL

*

GRTOP

=

3.5000000000D-03(K/M)

V.P.T.G. AT MODEL TOP

*

ROU

=

1.2000000000D+00(KG/M3) AIR DENSITY

*

TREF

=

12.0000000000D+00(C)

REFERENCE TEMPERATURE*

TBOTOM=

12.0000000000D+00(C)

SOIL BOTTOM TEMPERATURE *TWATER=

12.0000000000D+00(C)

WATER TEMPERATUREEMSVTY

=

0.9800000000D-00() EMMISSIVITY OF G. SFC.

*ALBEDO=

0.1400000000D-00() ALBEDO*TURBID=

0.1000000000D-00(

Z0=1.0000000000D-02
Continued on the following page.

)(M)－ 28 －OF G. SFC.

TURBIDITY OF AIR**

SURFACE ROUGHNESS (WIND)*JAEA-Data/Code 2008-033

ZT
DELT

=
=

1.0000000000D-03(M)

SURFACE ROUGHNESS (TEMP)*

30.0000000000D+00(S)

TIME INCREAMENT (S)(Y4M2D2)INITIAL DATE*

(H2M2S2)INITIAL TIME (STD.T.)**

LDATES=

20030207<=========LTIMES=

140000<=========TINTEGD=

358.0000000000D+00(DAY)

INTEGRATION PERIOD (DAY)*TINTEGH=

10.0000000000D+00(H)

INTEGRATION PERIOD (H)TINTINP= 1800.0000000000D+00(S)

INPUT DATA INTERVAL (S) *

FLON

=

17.2330000000D+00(DEG)

LONGITUDE

*

FLAT

=

28.1670000000D+00(DEG)

LATITUDE*

STDLON=

0.0000000000D+00(DEG)

LON. of STANDARD TIME

*

STYPE

=

16.0000000000D+00() SOIL TEXTURE NUMBER

*

UTYPE

=

1.0000000000D+00() LANDUSE TYPE NUMBER*

SFMOIS=

0.6400000000D+00(M3/M3) SURFACE SOIL WATER CONT.*SBMOIS=

0.3000000000D+00(M3/M3) SOIL BOTTOM WATER CONT. *SFTEMP=

280.1600000000D+00(K)

SOIL SURFACE TEMPERATURE*

CO2AP

=

0.0000000000D+00(ppm)

AIR CO2 CONCENTRATION + *CO2SI= 5000.0000000000D+00(ppm)

SOIL SURFACE CO2 CONC.

SS0=

0.0600000000D+00(m-mol~)SOIL CO2 PRODUCTION RATE*

SP0=

0.0100000000D+00(m-mol~)ROOT CO2 PRO. RATE / LAI*

AAZ=

10.5000000000D+00(M-1)

SOIL CO2 PARAMETER

*

DX

=

1.0000000000D+03(M)

GRID INTERVAL*

X00=-2673.0000000000D+03(M)

X COORDINATE of SW-POINT*

Y00=-2673.0000000000D+03(M)

Y COORDINATE of SW-POINT***

The model execution is done by submitting the NQS shell-script (go_1D.sh). Sample ofthe shell-script for three-dimensional MPI parallel calculation is shown in Table 5-6. The
SOLVEG2 offers several options for CO2, fog, and soil models. Here we outline thoseavailable in the modified version.
a)iffog:

fog deposition calculation (0 = no fog, 1 = include fog deposition)b)

npdsd:

Droplet size distribution of cloud water when iffog = 1 (1 = Best 48) , 2 =Deirmendjian 49) , 3 = Caumerlic et al. 50) , 4 = Klemm et al. 52))c)

ifdsl:

Soil water retention curve and phase change processes of soil water (0 =do not include, 1 = include)d)npwrc:

Soil water retention curve when ifdsl = 0 (1 = old scheme by Brooks-Corey18) , 2 = new scheme of complete water retention curve based on van

Genuchten 20) and Campbell and Shiozawa 21))e)nptms:

Soil thermal conductivity when ifdsl = 1 (1 = old scheme by McCumber
and Pielke 13) , 2 = new scheme using clay content by McInnes 15) )－ 29 －

JAEA-Data/Code 2008-033

f)

ifco2:

CO2 exchange process (0 = not include, 1 = include)

Table 5-6 Sample NQS shell-script for 3-D MPI parallel calculation: go_1D.sh.
#!/bin/csh –f######################################################## SOLVEG grid calculation########################################################BSUB -J SOL1D#BSUB -n 1
#BSUB -a intelmpi#BSUB -o sol.out#BSUB -e sol.err#######################################################
set hmdir=$HOME/SOLVEG_1D
set srdir=${hmdir}"/SRC"
set indir=${hmdir}"/INPUT"
set otdir=${hmdir}"/OUTPUT"
set grdir=${hmdir}"/GRIDOUT"# FOG CALCULATION (iffog) / 1:yes , 0:no##

> iffog = 1: CHOOSE DSD PARAMETERIZATION(npdsd=1-4)#> 1:BEST, 2: DEIRMENDJIAN(PRECISE BUT HEAVY),#> 3: CHAUMERLIAC, 4: KLEMM (EFFICTIVE SCHEME)#
set iffog="1"
set npdsd="4"# USE DRY SOIL MODEL (ifdsl) / 1:yes , 0:no##

> ifdsl = 0: CHOOSE WATER RETENTION CURVE (nswrc=1-2)#> 1: OLD (B-C), 2: NEW (van-G + C-S COMPLETE CURVE)##> ifdsl = 1: CHOOSE SOIL THERMAL CONDUCT.(nstms=1-2)#> 1: OLD, 2: NEW (MCINNES 1981)#
Continued on the following page.

－ 30 －

JAEA-Data/Code 2008-033

set ifdsl="1"
set nswrc="2"
set nstms="2"
set ifco2="0"# CO2 CALCULATION/ 1:yes , 0:no

set ifrst="0"# RESTRT CALCULATION / 1:yes , 0:no/bin/rm ${otdir}/*#/bin/rm -r ${grdir}/grads#/bin/rm -r ${grdir}/netcdf#/bin/mkdir ${grdir}/grads#/bin/mkdir ${grdir}/netcdf#input files
setenv FORT10 ${hmdir}/param_1D
setenv FORT12 ${indir}/metdata.dat
setenv FORT14 ${srdir}/zmesh.model_1Dsetenv FORT15 ${srdir}/BCsoil.tablesetenv FORT16 ${srdir}/zvege.table_1Dsetenv FORT17 ${srdir}/zvege.profile_1Dsetenv FORT18 ${srdir}/vGsoil.table#output files#---in pzroot.f
setenv FORT20 ${otdir}/dbout
setenv FORT23 ${otdir}/METout
setenv FORT24 ${otdir}/FLXout
setenv FORT25 ${otdir}/WNDout
setenv FORT26 ${otdir}/ACO2out
setenv FORT27 ${otdir}/BACO2out#---in solveg.f
setenv FORT30 ${otdir}/mnout
setenv FORT31 ${otdir}/TSout
setenv FORT32 ${otdir}/HWout
setenv FORT33 ${otdir}/QSout
setenv FORT34 ${otdir}/EBout
Continued on the following page.

－ 31 －

JAEA-Data/Code 2008-033

setenv FORT35 ${otdir}/SFout
setenv FORT36 ${otdir}/wcurve#---within the canopy
setenv FORT40 ${otdir}/VGout
setenv FORT41 ${otdir}/VWout
setenv FORT42 ${otdir}/VTout
setenv FORT43 ${otdir}/RADout#---fog and precipitation
setenv FORT44 ${otdir}/PREout#---CO2
setenv FORT50 ${otdir}/SCO2out
setenv FORT51 ${otdir}/VCO2out
setenv FORT52 ${otdir}/PSCO2out
setenv FORT53 ${otdir}/BSCO2out#go
cd ${hmdir}echo ${indir} >! EXpara_1D# INPUT DATA PATH

echo ${grdir} >> EXpara_1D# GRIDecho ${iffog} >> EXpara_1D# FOG CALCULATION

echo ${npdsd} >> EXpara_1D# IFFOG -> DSD PARAM.

echo ${ifdsl} >> EXpara_1D# USE ADSORPTION MODEL

echo ${nswrc} >> EXpara_1D# WATER RETENTION CURVE

echo ${nstms} >> EXpara_1D# SOIL THERMAL COND.

echo ${ifco2} >> EXpara_1D# CO2 CALCULATION

echo ${ifrst} >> EXpara_1D# RESTART CALCULATIONDATA PATH

setenv FORT9 ${hmdir}/EXpara_1D./zsolveg_1D.exe >& ${otdir}/outlist/bin/rm EXpara_1D－ 32 －

JAEA-Data/Code 2008-033

6.

Summary
In order to evaluate energy and water balances between atmosphere and land surface

accurately, sophisticated physical schemes such as evaporation and adsorption processes
in the soil and cloud (fog) water deposition on vegetation were implemented in the
one-dimensional atmosphere-soil-vegetation model including CO2 exchange process(SOLVEG2).

The

above

schemesaredesignedtosimulate

processes

intheatmosphere-soil-vegetation system under arid environment in a more realistic treatmentof physical processes than the commonly used land surface models. Performance tests inarid areas showed that the above schemes have a significant effect on surface energy andwater balances. The framework of the novel schemes for arid environments incorporated
to SOLVEG2 was documented in the present paper.
The concept of further modifications of SOVLEG2 is summarized in Fig 6-1. In the
future works, the processes of matter cycle due to atmospheric deposition (dry, wet andcloud water depositions) and plant uptake of nutrients, atmospheric pollutants, andradionuclide substances will be incorporated to the model for studies of material exchanges
between atmosphere and land surface. The carbon exchanges between the vegetation andsoil vegetation (such as litterfall) and stresses-growth-dynamics model which calculatebiomass, LAI and canopy height from available organic substances and nutrients will bealso implemented in the model. Meanwhile, the schemes of the snowfall interception bycanopy and soil water freezing should be developed to adapt the model for cold
environments. With such modifications made in the model, the model is suited to evaluate
environmental loads to ecosystems, and their responses under extreme environmental
conditions due to climate changes and human activities (e.g., global warming, air pollution,
drought).

－ 33 －

JAEA-Data/Code 2008-033

Fig.6-1 The concept of the atmosphere-soil-vegetation system based on the SOLVEG2. The
modules squared with dashed lines will be developed in future works.
Acknowledgement
This research was performed using the computer system of Altix3700Bx2 in the Centerfor Computational Science and e-Systems. This study was also supported by the Research
Revolution Plan 2002 (RR2002) of the Ministry of Education, Culture, Sports, Science, andTechnology.

－ 34 －

JAEA-Data/Code 2008-033

References
1) H. Nagai: "Atmosphere-soil-vegetation model including CO2 exchange processes:
SOLVEG2", JAERI-Data/Code 2004-014 (2004).
2) G. Katata, H, Nagai, H. Ueda, N. Agam, P.R. Berliner: "Development of a Land Surface
Model Including Evaporation and Adsorption Processes in the Soil for the Land-Air
Exchange in Arid Regions", Journal of Hydrometeorology, 8, pp. 1307-1324 (2007).
3) G. Katata, H. Nagai, T. Wrzesinsky, O. Klemm, W. Eugster, R. Burkard: "Development ofa Land Surface Model Including Cloud Water Deposition on Vegetation", Journal of Applied
Meteorology and Climatology, 47, pp. 2129-2146 (2008).
4) N. Agam (Ninari), P.R. Berliner: "Dew formation and water vapor adsorption insemi-arid environments - A review", Journal of Arid Environments, 65, pp. 572-590 (2006).
5) A.F.G. Jacobs, B.G. Heusinkveld, S.M. Berkowicz: "Dew measurements along a
longitudinal sand dune transect, Negev Desert, Israel", International Journal of
Biometeorology, 43, pp. 184-190 (1999).
6) N. Agam (Ninari), P.R. Berliner, A. Zangvil, E. Ben Dor: "Soil water evaporation during
the dry season in an arid zone", Journal of Geophysical Research, 109, D16103 (2004).
7) S.S.O. Burgess, T.E. Dawson: "The contribution of fog to the water relations of Sequoia
sempervirens (D. Don): Foliar uptake and prevention of dehydration", Plant, Cell andEnvironment, 27, pp. 1023-1034 (2004).
8) L.B. Hutley, D. Doley, D.J. Yates, A. Boonsaner: "Water balance of an australiansubtropical rainforest at altitude: The ecological and physiological significance ofintercepted cloud and fog", Australian Journal of Botany, 45, pp. 311-329 (1997).
9) R.D. Jackson, R.J. Reginato, B.A. Kimball, F.S. Nakayama: "Diurnal soil-water
evaporation: comparison of measured and calculated soil-water fluxes", Soil Science Societyof America proceedings, 38, pp. 861-866 (1974).
10) J. Grifoll, J.Ma. Gastó, Y. Cohen: "Non-isothermal soil water transport andevaporation", Advances in Water Resources, 28, pp. 1254-1266 (2005).
11) J. R. Philip, D. A. De Vries: "Moisture movement in porous materials undertemperature gradients", Transactions, American Geophysical Union, 38, pp. 222-232(1957).
12) A. T. Cahill, M. B. Parlange: "On water transport in field soils", Water Resources
Research, 34, pp. 731-739 (1998).
13) M. C. McCumber, R. A. Pielke: "Simulation of the effects of surface fluxes of heat andmoisture in a mesoscale numerical model 1. Soil layer", Journal of Geophysical Research,86, pp. 9929-9938 (1981).
14) C.D. Peters-Lidard, E. Blackburn, X. Liang, E.F. Wood: "The effect of soil thermal
conductivity parameterization on surface energy fluxes and temperatures", Journal of the－ 35 －

JAEA-Data/Code 2008-033

Atmospheric Sciences, 55, pp. 1209-1224 (1998).
15) K. McInnes: "Thermal conductivities of soils from dryland wheat regions in eastern
Washington", M.S. Thesis, Dept. of Agronomy and Soils, Washington State University,
Pullman, WA, 51 p [Available from Owen Science Library, Washington State University,
Pullman, WA 99163-3200.].
16) D.A. de Vries: "Thermal properties of soils", Physics of Plant Environment, 210-235,
W.R. van Wijk North Holland Publ. Co., Amsterdam (1963).
17) W. Brutsaert: "Evaporation Into the Atmosphere", D. Reidel, 299 p (1982).
18) R.H. Brooks, A.T. Corey: "Hydraulic properties of porous media", Hydrology paper 3,
Colorado State University, Fort Collins, 27 p (1964).
19) S. W. Webb: "A simple extension of two-phase characteristic curves to include the dry
region", Water Resources Research, 36, pp. 1425-1430 (2000).
20) M.Th. van Genuchten: "Closed-form equation for predicting the hydraulic conductivityof unsaturated soils", Soil Science Society of America Journal, 44, pp. 892-898 (1980).
21) G.S. Campbell, S. Shiozawa: "Prediction of hydraulic properties of soils using
particle-size distribution and bulk density data", International Workshop on Indirect
Methods for Estimating the Hydraulic Properties of Unsaturated Soils: 317-328, Univ. of
Calif., Riverside, Calif (1992).
22) Y. Mualem: "New model for predicting the hydraulic conductivity of unsaturated porous
media", Water Resources Research, 12, pp. 513-522 (1976).
23) R.F. Carsel, R.S. Parrish: "Developing joint probability distributions of soil water
retention characteristics", Water Resources Research, 24, pp. 755-769 (1988).
24) J. Kondo, T. Watanabe: "Studies on the bulk transfer coefficients over a vegetatedsurface with a multilayer energy budget model", Journal of the Atmospheric Sciences, 49,
pp. 2183-2199 (1992).
25) B.V. Derjaguin, N.V. Churaev: "On the question of determining the concept of disjoiningpressure and its role in the equilibrium and flow of thin films", Journal of Colloid andInterface Science, 66, pp. 389-398 (1978).
26) M. Iwamatsu, K. Horii: "Capillary condensation and adhesion of two wetter surfaces",
Journal of Colloid and Interface Science, 182, pp. 400-406 (1996).
27) D. Or, M. Tuller: "Liquid retention and interfacial area in variably saturated porous
media: Upscaling from single-pore to sample-scale model", Water Resources Research, 35,
pp. 3591-3605 (1999).
28) M. Tuller, O. Dani, L.M. Dudley: "Adsorption and capillary condensation in porous
media: Liquid retention and interfacial configurations in angular pores", Water Resources
Research, 35, pp. 1949-1964 (1999).
29) L.W. Petersen, P. Moldrup, O.H. Jacobsen, D.E. Rolston: "Relations between specificsurface area and soil physical and chemical properties", Soil Science, 161, pp. 9-20 (1996).
－ 36 －

JAEA-Data/Code 2008-033

30) A. Banin, A. Amiel: "A correlative study of the chemical and physical properties of agroup of natural soils of Israel", Geoderma, 3, pp. 185-198 (1970).
31) K. Kosugi: "Three-parameter lognormal distribution model for soil water retention",
Water Resources Research, 30, pp. 891-901 (1994).
32) S. Brunauer, P.H. Emmett, E. Teller: "Adsorption of gases in multimolecular layers",
Journal of American Chemical Society, 60, pp. 309-319 (1938).
33) R.B. Bird, W.E. Stewart, E.N. Lightfoot: "Transport Phenomena", Wiley, 912 p (2001).
34) D.E. Rolston: "Gas flux", In: Klute, A. ed., "Methods of Soil Analysis, Part 1: Physicaland Mineralogical Methods", Agronomy Monograph No.9, American Society of Agronomy,
Madison, pp. 1103-1119 (1986).
35) T.C. Ransohoff, C.J. Radke: "Laminar flow of a wetting liquid along the corners of apredominantly gas-occupied noncircular pore", Journal of Colloid and Interface Science,
121, pp. 392-401 (1988).
36) M.J. Blunt, H. Scher: "Pore-level modeling of wetting", Physical Review E, 52, pp.
6387-6403 (1995).
37) J.N. Israelachvili: "Intermolecular and Surface Forces, 2nd Ed.", Academic Press,
London, 450 p (1992).
38) P.M. Driessen: "The water balance of the soil. In: Van Keulen H, Wolf J (Eds). Modellingof agricultural production: weather, soils and crops", Simulation Monographs, Pudoc,,
Wageningen, pp. 76-116 (1986).
39) J.C. Schaake, V.I. Koren, Q.Y. Duan, K. Mitchell, F. Chen: "Simple water balance modelfor estimating runoff at different spatial and temporal scales", Journal of GeophysicalResearch, 101, pp. 7461-7475 (1996).
40) F. Chen, J. Dudhia: "Coupling an advanced land surface-hydrology model with the Penn
State-NCAR MM5 modeling system. Part I: Model implementation and sensitivity",Monthly Weather Review, 129, pp. 569-585 (2001).
41) T. Yamada: "A numerical simulation of nocturnal drainage flow", Journal of the
Meteorological Society of Japan, 59, pp. 108-122 (1981).
42) H. Nagai: "Validation and sensitivity analysis of a new atmosphere-soil-vegetation
model. Part II: Impacts on in-canopy latent heat flux over a winter wheat field determinedby detailed calculation of canopy radiation transmission and stomatal resistance", Journalof Applied Meteorology, 42, pp. 434-451 (2003).
43) K. Peters, R. Eiden: "Modelling the dry deposition velocity of aerosol particles to a
spruce forest", Atmospheric Environment, 26A, pp. 2555-2564 (1992).
44) P.G. Thorne, G.M. Lovett, W.A. Reiners: "Experimental determination of droplet
impaction on canopy components of balsam fir", Journal of Applied Meteorology, 21, pp.
1413-1416 (1982).
45) Y. Belot, D. Gauthier: "Transport of micronic particles from atmosphere to folier－ 37 －

JAEA-Data/Code 2008-033

surfaces", Semin on Heat and Mass Transfer in the Environ of Veg, Heat and Mass Transfer
in the Biosphere (Pt 1) pp. 583-591 (1975).
46) A.C. Chamberlain: "Transport of Lycopodium spores and other small particles to roughsurfaces", Proceedings of the Royal Society of London Series A-Mathmatical Physical, 296,
pp. 45-70 (1967).
47) D.E. Aylor: "Modeling spore dispersal in a barley crop", Agricultural Meteorology, 26, pp.
215-219 (1982).
48) G.M. Lovett: "Rates and mechanisms of cloud water deposition to a subalpine balsam
fir forest", Atmospheric Environment, 18, pp. 361-371 (1984).
49) A.C. Best: "Drop-size distribution in cloud and fog", Quarterly Journal of the Royal
Meteorological Society, 77, pp. 418-426 (1951).
50) D. Deirmendjian: "Electromagnetic Scattering on Spherical Polydispersions", Elsevier,
New York, 290 p (1969).
51) N. Chaumerliac, E. Richard, J.-P., Pinty, E.C. Nickerson: "Sulfur scavenging in amesoscale model with quasi-spectral microphysics: Two-dimensional results for continental
and maritime clouds", Journal of Geophysical Research, 92, pp. 3114-3126 (1987).
52) A. Bott, T. Trautmann: "PAFOG - a new efficient forecast model of radiation fog andlow-level stratiform clouds", Atmospheric Research, 64, pp. 191-203 (2002).
53) O. Klemm, T. Wrzesinsky, C. Scheer: "Fog water flux at a canopy top: Direct
measurement versus one-dimensional model", Atmospheric Environment, 39, pp.
5375-5386 (2005).
54) W.E. Ranz, W.R. Marshall: "Evaporation from drops, Part II", Chemical Engineering
Progress, 48, pp. 173-180 (1952).
55) R.D. Magarey, J.M. Russo, R.C. Seem, D.M. Gadoury: "Surface wetness duration undercontrolled environmental conditions", Agricultural and Forest Meteorology, 128, pp.
111-122 (2005).
56) J.W. Deardorff: "Efficient prediction of ground surface temperature and moisture, withinclusion of a layer of vegetation", Journal of Geophysical Research, 83, pp. 1889-1903(1978).

－ 38 －

JAEA-Data/Code 2008-033

Appendix: Description of model code
The source files of the SOLVEG2 are described based on the categories of main routine,
primary routines, secondary routines, and functions.
A-1 Main Routine
SOLVEG2 (main.f)
Function:- Reading and setting the input and output directory paths and physics options for thedry soil model, fog deposition, CO2 and restart calculations.
- Controlling calculation flow and time step.
Include files: Inclnum, Inclcon1, Inclcon2, Inclvari, prm_soil, prm_fog.

A-2 Primary Routines
PREAD (ppread.f)
Function:- Reading parameters from param_1D.
- Reading initial data from files in the input directory.
Include files: Inclnum, Inclcon1, Inclcon2, Inclvari, prm_soil, prm_fog.
Variables: None.
MSHINT (NX, NY, N2, UTY, Z, ZZD) (main.f)
Function:- Setting the atmosphere grid from data in zmesh.model_1D.
Include files: None.
Variables:- NX:

Integer Horizontal x-grid size.

- NY:

Integer Horizontal y-grid size.

- N2:

Integer Vertical grid size in the atmosphere.

- UTY(NX,NY):

Real*8 Land-use type [1 to 24 USGS categories].

- Z(NX,NY,N2):

Real*8 3-D grid height in the atmosphere [m].

- ZZD(NX,NY,N2):

Real*8 [z(k+1)-z(k-1)]-1 [m -1 ].

SLVGIN (ZI , DELTI , LDATEYI , TIMESI, LDATEI , LTIMEI , FLATI , TLAGLI , Z0I , ZTI,TBOTMI , TWATRI , TSI , HWI , STY, UTY, CO2SI, IFRST, PATHG, EMSVTYI,ALBEDOI , TURBIDI , PHPS, T0I , QS0I , CO20I , CTIME1 , ZSO) (solveg.f)－ 39 －

JAEA-Data/Code 2008-033

Function:- Initializing soil and vegetation routines.
Include files: Inclnum, prm_soil.
Variables:- ZI(NX,NY,N2):

Real*8 3-D grid height in the atmosphere [m].

- DELTI:

Real*8 Time increment [s].

- LDATEYI:

Integer Julian day [day].

- TIMESI:

Real*8 Calculation start time [h].

- LDATEI:

Integer Date [yyyymmdd].

- LTIMEI:

Integer Time [hhmmss].

- FLATI(NX,NY):

Real*8 Latitude of model grid [deg.].

- TLAGLI(NX,NY):

Real*8 Local time difference from standard time [h].

- Z0I(NX,NY):

Real*8 Roughness height for momentum [m].

- ZTI(NX,NY):

Real*8 Roughness height for heat [m].

- TBOTMI(NX,NY):

Real*8 Soil bottom temperature [K].

- TWATRI(NX,NY):

Real*8 Sea surface temperature [K].

- TSI(NX,NY):

Real*8 Initial soil surface temperature [K].

- HWI(NX,NY):

Real*8 Initial soil water content [m3 m -3].

- STY(NX,NY):

Real*8 Soil type [1 to 16].

- UTY(NX,NY):

Real*8 Land-use type [1 to 24 USGS categories].

- CO2SI:

Real*8 Initial soil CO2 concentration [ppmv].

- IFRST:

Integer An index for restart calculations (0 or 1).

- PATHG:Char.

- EMSVTYI(NX,NY):

Real*8 Surface emissivity.

- ALBEDOI(NX,NY):

Real*8 Surface albedo.

- TURBIDI(NX,NY):

Real*8 Turbidity of air.

- PHPS(NX,NY):

Real*8 Surface air pressure [hPa].

- T0I(NX,NY):

Real*8 Initial soil surface temperature [°C].

- QS0I(NX,NY):

Real*8 Initial soil surface specific humidity [kg kg-1].

- CO20I(NX,NY):

Real*8 Initial soil surface CO2 concentration [ppmv].

- CTIME1

Char*15 Calculation time [yyyy-mm-dd_hhmm].

- ZSO(NS+1)

Real*8 Soil grid depth [m].

The directory path for restart files.

INITPF (pinitpf.f)
Function:- Initializing atmospheric variables.
Include files: Inclnum, Inclcon1, Inclcon2, and Inclvari.
Variables: None.
－ 40 －

JAEA-Data/Code 2008-033

FLXCAL (pfluxcal.f)Function:- Calculating turbulent flux.
Include files: Inclnum, Inclcon1, Inclcon2, and Inclvari.
Variables: None.
SFPR13 (pfluxcal.f)Function:- Calculating stability, friction velocity and temperature.
Include files: Inclnum, Inclcon1, Inclcon2, and Inclvari.
Variables: None.
DEBUGW (IU, IS, IE, N, K1, K2, K3, K4, IPG) (pdebugw0.f)
Function:- Outputting atmospheric variables.
Include files: Inclnum, Inclcon1, Inclcon2, and Inclvari.
Variables:- IU:

Integer File output unit number.

- IS:

Integer Start layer number of output.

- IE:

Integer End layer number of output.

- N:

Integer Interval of output layers.

- K1:

Integer Mean variable output flag [1: on, 0: off].

- K2:

Integer Turbulence variable output flag [1: on, 0: off].

- K3:

Integer Turbulence diffusivity output flag [1: on, 0: off].

- K4:

Integer Turbulence term output flag [1: on, 0: off].

- IPG:

Integer Paging flag.

TIMEINT (ptint.f)
Function:- Reading boundary condition from files in the input directory or receive them from
MM5.
- Interpolating input data into the model time step.
Include files: Inclnum, Inclcon1, Inclcon2, Inclvari, prm_fog.
Variables: None.
KMHCAL (pfluxcal.f)Function:- Calculating eddy diffusivity by level 2.5 turbulence closure model.
－ 41 －

JAEA-Data/Code 2008-033

Include files: Inclnum, Inclcon1, Inclcon2, and Inclvari.
Variables: None.
SOLVEG (IFDBG, IFCO2, LDATEY, TIMEH, TIMET, LDATE, LTIME, ROU, PHPS, U,
V, T, Q, WL, CO2, PRECIP, STBL, CNL, CNM, CNH, RSOL, RINF, RRATE,RSOLX, RINFX, T0, QS0, CO20, ALB, TEAS, FSH, FLH, FLQ, FGH, FLW,FNR, HWO, TSO, SOUV, SOTF, SOTA, SOQF, SOQA, SOWL, SOE2, AZCD,
SOCO2, SS0, SP0, AAZ) (solveg.f)
Function:- Controlling soil and vegetation routines.
- Interfacing atmosphere, soil, and vegetation sub-models.
Include files: Inclnum, prm_soil.
Variables:- IFDBG:

Integer Output flag [1: on, 0: off].

- IFCO2:Integer CO2 calculation flag [1: on, 0: off].

- LDATEY:

Integer Julian day [day].

- TIMEH:

Real*8 Local time [h].

- TIMET

Real*8 Integration time [h].

- LDATE:

Integer Date [yyyymmdd].

- LTIME:

Integer Time [hhmmss].

- ROU(NX,NY):

Real*8 Air density [kg m-3].

- PHPS(NX,NY):

Real*8 Surface air pressure [hPa].

- U(NX,NY,M1):

Real*8 Wind u-component [m s -1].

- V(NX,NY,M1):

Real*8 Wind v-component [m s-1].

- T(NX,NY,M1):

Real*8 Air temperature [°C].

- Q(NX,NY,M1):

Real*8 Specific humidity [kg kg-1 ].

- WL(NX,NY,M1):

Real*8 Fog water content [kg kg-1].

- CO2(NX,NY,M1):

Real*8 CO2 concentration [ppmv].

- PRECIP(NX,NY):

Real*8 Precipitability [cm].

- STBL(NX,NY):

Real*8 Stability length [m].

- CNL(NX,NY):

Real*8 Low-level cloudiness [0 to 1 fraction].

- CNM(NX,NY):

Real*8 Middle-level cloudiness [0 to 1 fraction].

- CNH(NX,NY):

Real*8 High-level cloudiness [0 to 1 fraction].

- RSOL(NX,NY):

Real*8 Solar radiation flux [W m-2].

- RINF(NX,NY):

Real*8 Long-wave radiation flux [W m -2 ].

- RRATE(NX,NY):

Real*8 Rain intensity [mm h-1].

- RSOLX(NX,NY):

Real*8 Model solar radiation flux [W m-2 ].

- RINFX(NX,NY):

Real*8 Model long-wave radiation flux [W m -2].
－ 42 －

JAEA-Data/Code 2008-033- T0(NX,NY):

Real*8 Soil surface temperature [°C].

- QS0(NX,NY):

Real*8 Soil surface specific humidity [kg kg-1 ].

- CO20(NX,NY):

Real*8 Soil surface CO2 concentration [ppmv].

- ALB(NX,NY):

Real*8 Surface albedo.

- TEAS(NX,NY):

Real*8 Energy consumption by photosynthesis [W m-2].

- FSH(NX,NY):

Real*8 Sensible heat flux [W m -2].

- FLH(NX,NY):

Real*8 Latent heat flux [W m -2 ].

- FLQ(NX,NY):

Real*8 Vapor flux [kg m -2 s-1 ].

- FGH(NX,NY):

Real*8 Soil heat flux [W m -2 ].

- FLW(NX,NY):

Real*8 Upward long-wave radiation flux [W m -2].

- FNR(NX,NY):

Real*8 Net radiation [W m -2].

- HWO(NX,NY,NS):

Real*8 Soil water content [m 3 m-3].

- TSO(NX,NY,NS):

Real*8 Soil temperature [K].

- SOUV(NX,NY,M1):

Real*8 Momentum exchange term by vegetation [s -1 ].

- SOTF(NX,NY,M1):

Real*8 Heat source term by vegetation [K s -1 ].

- SOTA(NX,NY,M1):

Real*8 Heat exchange term by vegetation [s-1].

- SOQF(NX,NY,M1):

Real*8 Vapor source term by vegetation [kg kg -1 s-1 ].

- SOQA(NX,NY,M1):

Real*8 Vapor exchange term by vegetation [s-1].

- SOWL(NX,NY,N1):

Real*8 Fog source term by vegetation [kg kg -1 s-1].

- SOE2(NX,NY,M1):

Real*8 Turbulence source term by vegetation [m 0.5 s -1.5 ].

- AZCD(NX,NY,M1):

Real*8 a×cD [m 2 m -3].

- SOCO2(NX,NY,M1):

Real*8 CO2 source term by vegetation [ppmv s -1 ].

- SS0:

Real*8 Optimal soil CO2 production rate [μmol m -2 s -1 ].

- SP0:

Real*8 Optimal root CO2 production rate [μmol m -2 s -1 ].

- AAZ:

Real*8 CO2 production distribution parameter in thesoil [m-1 ].

UMAIN (SOUV) (pmain03.f)
Function:- Calculating wind components u and v.
Include files: Inclnum, Inclcon1, Inclcon2, and Inclvari.
Variables:- SOUV(NX,NY,M1):

Real*8 Momentum exchange term by vegetation [s -1 ].

TMAIN (SOTF, SOTA, SOQF, SOQA, SOWL) (pmain03.f)
Function:- Calculating temperature, specific humidity, and fog water.
Include files: Inclnum, Inclcon1, Inclcon2, and Inclvari.
－ 43 －

JAEA-Data/Code 2008-033

Variables:- SOTF(NX,NY,M1):

Real*8 Heat source term by vegetation [K s -1 ].

- SOTA(NX,NY,M1):

Real*8 Heat exchange term by vegetation [s-1].

- SOQF(NX,NY,M1):

Real*8 Vapor source term by vegetation [kg kg -1 s-1 ].

- SOQA(NX,NY,M1):

Real*8 Vapor exchange term by vegetation [s-1].

- SOWL(NX,NY,N1):

Real*8 Fog source term by vegetation [kg kg -1 s-1].

EMAIN (SOE2) (pmain03.f)
Function:- Calculating turbulence energy and length scale.
Include files: Inclnum, Inclcon1, Inclcon2, and Inclvari.
Variables:- SOE2(NX,NY,M1):

Real*8 Turbulence source term by vegetation [m 0.5 s -1.5 ].

CMAIN (SOCO2) (pmain03.f)
Function:- Calculating atmospheric CO2 concentration.
Include files: Inclnum, Inclcon1, Inclcon2, and Inclvari.
Variables:- SOCO2(NX,NY,M1):

Real*8 CO2 source term by vegetation [ppmv s -1 ].

SAVEOD (AZCD) (pfluxcal.f)Function:- Setting variables for the next time step.
- Calculating the top boundary values of turbulence energy and length scale.
Include files: Inclnum, Inclcon1, Inclcon2, and Inclvari.
Variables:- AZCD(NX,NY,M1):

Real*8 a×cD [m 2 m -3].

A-3 Secondary Routines
GZSOLVEG (NX, NY, NS, NC, UTY, ZS, ZC) (gzsolveg.f)
Function:- Setting soil and vegetation grids from data in zmesh.grid_3D and zmesh.model_1D/3D.
Include files: None.
Variables:- NX:

Integer Horizontal x-grid size.
－ 44 －

JAEA-Data/Code 2008-033- NY:

Integer Horizontal y-grid size.

- NS:

Integer Vertical grid size in soil.

- NC:

Integer Vertical grid size in vegetation canopy.

- UTY(NX,NY):

Real*8 Land-use type [1 to 24 USGS categories].

- ZS(0:NS+1):

Real*8 Vertical grid in soil [m].

- ZC(NX,NY,0:NC+1):

Real*8 3-D grid height in vegetation canopy [m].

GTABLE (gtable.f)Function:- Setting soil parameter from data in zoil.table.
Include files: None.
Variables: None.
GVTABLE(CHANGD,CHANGEH,

CHANGEE,

WDSAT,WDEVP,FANGLE,

REFLECT, TRANSMT, EMISSIV, RSSTMIN, PARCRT, FGRVGE, VMAX25,
SLOPE) (gvtable.f)Function:- Setting vegetation parameter from data in zvege.table_1D/3D.
Include files: None.
Variables:- CHANGD(0:99):

Real*8 Aerodynamic resistance of leaf.

- CHANGEH(0:99):

Real*8 Heat exchange coefficient between leaf and air.

- CHANGEE(0:99):

Real*8 Vapor exchange coefficient between leaf and air.

- WDSAT(0:99):

Real*8 Maximum leaf surface water [kg m -2 ].

- WDEVP(0:99):

Real*8 Leaf surface water for max. evaporation [kg m -2 ].

- FANGLE(0:99):

Real*8 Mean zenith angle of leaf orientation [deg.].

- REFLECT(2,0:99):

Real*8 Leaf reflectivity for visible and near-infrared.

- TRANSMT(2,0:99):

Real*8 Leaf transmissivity for visible and near-infrared.

- EMISSIV(0:99):

Real*8 Leaf emissivity.

- RSSTMIN(0:99):

Real*8 Minimum stomatal resistance [s m -1].

- PARCRT(0:99):

Real*8 Critical value of PAR [W m -2 ].

- FGRVGE(0:99):

Real*8 Green leaf fraction [0 to 1 fraction].

- VMAX25(0:99):

Real*8 Vm at 25 °C [μmol m -2 s-1].

- SLOPE (0:99):

Real*8 Empirical constant m for stomatal resistance.

WADSP0 (swadsp.f)
Function:- Calculating water retention curve and phase change processes in the soil.
－ 45 －

JAEA-Data/Code 2008-033

Include files: Inclnum and prm_soil.
Variables: None.
LINEINT (X, Y, XX, YY, N) (swadsp.f)
Function:- Linearly interpolating the data.
Include files: prm_soil.
Variables:- MAXA:

Integer Maximum number of arrays.

- X:

Real*8 Target of x value.

- Y:

Real*8 Target of y value.

- XX(0:MAXA):

Real*8 X arrays for linear interpolation.

- YY(0:MAXA):

Real*8 Y arrays for linear interpolation.

- N:

Integer Number of interpolation data.

GVPROFILE (NS, NC, IDATEV, ITIMEV, VTYPEI, AZI, RZI) (gvprofile.f)Function:- Setting vegetation profile from data in zvege.profile_1D/3D.
Include files: None.
Variables:- NS:

Integer Vertical grid size in soil.

- NC:

Integer Vertical grid size in vegetation canopy.

- IDATEV:

Integer Date in vegetation profile data [yymmdd].

- ITIMEV:

Integer Time in vegetation profile data [hhmmss].

- VTYPEI(0:NC,30,2):Integer Vegetation type.

- AZI(0:NC,30,2):

Real*8 Leaf area density [m 2 m -3 ].

- RZI(0:NC,0:NS,30,2):

Real*8 Root distribution [0 to 1 fraction].

PCAL (ppcal.f)Function:- Calculating air pressure profile.
Include files: Inclnum, Inclcon1, Inclcon2, and Inclvari.
Variables: None.
CLSL2A (pinitpf.f)
Function:- Calculating turbulence variables by level 2.0 turbulence closure model.
Include files: Inclnum, Inclcon1, Inclcon2, and Inclvari.
－ 46 －

JAEA-Data/Code 2008-033

Variables: None.
DEWTMP (QA, PHPS, TDEW) (ppread.f)
Function:- Calculating dew temperature from specific humidity.
Include files: None.
Variables:- QA

Real*8 Specific humidity [kg kg -1 ].

- PHPS

Real*8 Air pressure [hPa].

- TDEW

Real*8 Dew temperature [K].

SFCRAD (ZC, AZ, FLAT, TLAGL, LDATEY, TIMEH, CNL, CNM, CNH, RSOL, RINF,
PHPS,

EMSVTY,ALBEDOD,ALBEDOS,PRECIP,TURBID,

REFLECT,

TRANSMT, EMISSIV, FANGLE, ROU, T, WL, TTC, WLC, TC, T0, PR, VF, RDDN,
RSDN, RSUP, RLDN, RLUP, PARABS, SUNABS, FLIT, HR, ZENITH, RSMX,RSOLX, RINFX) (srad.f)
Function:- Calculating radiation transmission in vegetation canopy.
Include files: Inclnum.
Variables:- ZC(NX,NY,0:NC+1):

Real*8 3-D grid height in vegetation canopy [m].

- AZ(NX,NY,0:NC):

Real*8 Leaf area density [m2 m -3 ].

- FLAT(NX,NY):

Real*8 Latitude of model grid [deg.].

- TLAGL(NX,NY):

Real*8 Local time difference from standard time [h].

- LDATEY:

Integer Julian day [day].

- TIMEH:

Real*8 Local time [h].

- CNL(NX,NY):

Real*8 Low-level cloudiness [0 to 1 fraction].

- CNM(NX,NY):

Real*8 Middle-level cloudiness [0 to 1 fraction].

- CNH(NX,NY):

Real*8 High-level cloudiness [0 to 1 fraction].

- RSOL(NX,NY):

Real*8 Solar radiation flux [W m-2].

- RINF(NX,NY):

Real*8 Long-wave radiation flux [W m -2 ].

- PHPS(NX,NY):

Real*8 Surface air pressure [hPa].

- EMSVTY(NX,NY):

Real*8 Surface emissivity.

- ALBEDOD(2,NX,NY):

Real*8 Surface albedo for direct solar radiation.

- ALBEDOS(2,NX,NY):

Real*8 Surface albedo for diffuse solar radiation.

- PRECIP(NX,NY):

Real*8 Precipitability [cm].

- TURBID(NX,NY):

Real*8 Turbidity of air.

- REFLECT(2,NX,NY,0:NC): Real*8 Leaf reflectivity for visible and near-infrared.
－ 47 －

JAEA-Data/Code 2008-033- TRANSMT(2,NX,NY,0:NC): Real*8 Leaf transmissivity for visible and near-infrared.
- EMISSIV(NX,NY,0:NC):

Real*8 Leaf emissivity.

- FANGLE(NX,NY,0:NC):

Real*8 Mean zenith angle of leaf orientation [deg.].

- ROU(NX,NY):

Real*8 Air density [kg m-3].

- T(NX,NY,M1):

Real*8 Air temperature [°C].

- WL(NX,NY,M1):

Real*8 Fog water content [kg kg-1].

- TTC(NX,NY,0:NC):

Real*8 Canopy air temperature [°C].

- WLC(NX,NY,0:NC):

Real*8 Canopy air fog water content [kg kg -1].

- TC(2,NX,NY,0:NC):

Real*8 Temperature for sun-lit/shaded leaf [°C].

- T0(NX,NY):

Real*8 Soil surface temperature [°C].

- PR(NX,NY,0:NC):

Real*8 Vertical water flux in canopy [kg m -2 s-1 ].

- VF(NX,NY):

Real*8 Vertical speed of precipitation [m s-1].

- RDDN(2,NX,NY,0:NC):

Real*8 Downward direct solar radiation flux [W m -2].

- RSDN(2,NX,NY,0:NC):

Real*8 Downward diffuse solar radiation flux [W m -2].

- RSUP(2,NX,NY,0:NC):

Real*8 Upward diffuse solar radiation flux [W m-2].

- RLDN(2,NX,NY,0:NC):

Real*8 Downward long-wave radiation flux [W m -2].

- RLUP(2,NX,NY,0:NC):

Real*8 Upward long-wave radiation flux [W m -2].

- PARABS(2,NX,NY,NC):

Real*8 PAR absorbed by leaves [W m-2].

- SUNABS(2,NX,NY,NC):

Real*8 Solar radiation absorbed by leaves [W m-2].

- FLIT(NX,NY,NC):

Real*8 Fraction of sun-lit leaves [0 to 1 fraction].

- HR(NX,NY,0:NC):

Real*8 Heating of canopy air by radiation [J m -3 s -1].

- ZENITH(NX,NY):

Real*8 Solar zenith angle [deg.].

- RSMX:

Real*8 Maximum solar radiation flux [W m -2].

- RSOLX(NX,NY):

Real*8 Model solar radiation flux [W m-2 ].

- RINFX(NX,NY):

Real*8 Model long-wave radiation flux [W m -2].

RADIATION (DFLAT, TLAGL, DSLAT, TIME, PRECIP, TURBID, ALBEDO, TEMP,
CNL, CNM, CNH, ZENITH, RSMX, RSOLX, RINFX) (sradiatn.f)
Function:- Calculating solar and long-wave radiation flux.
Include files: None.
Variables:- DFLAT:Real*8 Latitude [deg.].

- TLAGL:

Real*8 Local time difference from standard time [h].

- DSLAT:Real*8 Solar latitude [deg.].

- TIME:

Real*8 Model time [h].

- PRECIP:

Real*8 Precipitability [cm].

- TURBID:

Real*8 Turbidity of air.
－ 48 －

JAEA-Data/Code 2008-033- ALBEDO:

Real*8 Surface albedo.

- TEMP:

Real*8 Surface layer temperature [K].

- CNL:

Real*8 Low-level cloudiness [0 to 1 fraction].

- CNM:

Real*8 Middle-level cloudiness [0 to 1 fraction].

- CNH:

Real*8 High-level cloudiness [0 to 1 fraction].

- ZENITH:

Real*8 Solar zenith angle [deg.].

- RSMX:

Real*8 Maximum solar radiation flux [W m -2].

- RSOLX:

Real*8 Model solar radiation flux [W m-2 ].

- RINFX:

Real*8 Model long-wave radiation flux [W m -2].

EPPARA (NC, DZC, AZ, FA, FR) (eppara.f)
Function:- Calculating shielding coefficient of canopy layers.
Include files: None.
Variables:- NC:

Integer Vertical grid size in vegetation canopy.

- DZC(0:NC):

Real*8 Depth of vegetation layers [m].

- AZ(0:NC):

Real*8 Leaf area density [m 2 m -3 ].

- FA(0:NC):

Real*8 Mean projection coefficient of leaves.

- FR(0:NC):

Real*8 Shielding coefficient of canopy layers.

SVAPO (Z, TEXTURE, MYU, DELT, PHPS, TR, UR, QR, CE0, AIRDEN, HW, EB, T,
QS, VFS, RES) (svapo.f)Function:- Calculating vapor transport in soil.
Include files: Inclnum, prm_soil.
Variables:- Z(0:NS+1):

Real*8 Vertical grid in soil [m].

- TEXTURE(NX,NY,0:NS):

Integer Soil type [1 to 24 USGS soil type].

- MYU:

Real*8 Parameter for numerical scheme.

- DELT:

Real*8 Time increment [s].

- PHPS(NX,NY):

Real*8 Surface air pressure [hPa].

- TR(NX,NY):

Real*8 Lowest atmosphere layer temperature [K].

- UR(NX,NY):

Real*8 Lowest atmosphere layer wind speed [m s -1 ].

- QR(NX,NY):

Real*8 Lowest atmosphere layer humidity [kg kg -1 ].

- CE0(NX,NY):

Real*8 Vapor exchange coefficient at soil surface.

- AIRDEN(NX,NY):

Real*8 Air density [kg m-3].

- HW(NX,NY,0:NS):

Real*8 Soil water content [m 3 m-3].
－ 49 －

JAEA-Data/Code 2008-033- EB(NX,NY,0:NS):

Real*8 Evaporation rate of soil water [kg m-3 s-1 ].

- T(NX,NY,0:NS):

Real*8 Soil temperature [K].

- QS(NX,NY,0:NS):

Real*8 Specific humidity of soil pore [kg kg-1].

- VFS(NX,NY,0:NS):

Real*8 Vapor flux in soil [kg m -2 s -1].

- RES(NX,NY):

Real*8 Soil surface water [kg m -2 ].

SOLVE1 (NX, NY, NZ, A, B, C, D, Q, WK1, WK2) (solver1.f)
Function:- Numerical scheme for diffusion equation in soil with constant bottom boundary.
Include files: None.
Variables:- NX:

Integer Horizontal x-grid size.

- NY:

Integer Horizontal y-grid size.

- NZ:

Integer Vertical grid size.

- A(NX,NY,0:NZ):

Real*8 Coefficient (7-8) in the numerical scheme.

- B(NX,NY,0:NZ):

Real*8 Coefficient (7-9) in the numerical scheme.

- C(NX,NY,0:NZ):

Real*8 Coefficient (7-10) in the numerical scheme.

- D(NX,NY,0:NZ):

Real*8 Coefficient (7-11) in the numerical scheme.

- Q(NX,NY,0:NZ):

Real*8 Solution.

- WK1(NX,NY,0:NZ):

Real*8 Temporal variable.

- WK2(NX,NY,0:NZ):

Real*8 Temporal variable.

EVPARA (NX, NY, NS, TEXTURE, HW, T, CDIF, RB) (evpara.f)
Function:- Calculating vapor diffusivity in soil and evaporation resistance of soil water.
Include files: None.
Variables:- NX:

Integer Horizontal x-grid size.

- NY:

Integer Horizontal y-grid size.

- NS:

Integer Vertical grid size in soil.

- TEXTURE(NX,NY,0:NS):

Integer Soil type [1 to 24 USGS soil type].

- HW(NX,NY,0:NS):

Real*8 Soil water content [m 3 m-3].

- T(NX,NY,0:NS):

Real*8 Soil temperature [K].

- CDIF(NX,NY,0:NS):

Real*8 Vapor diffusivity in soil [m2 s -1 ].

- RB(NX,NY,0:NS):

Real*8 Evaporation resistance of soil water [s m-1 ].

EHWS (NX, NY, NS, TEXTURE, HWS) (ehws.f)
Function:－ 50 －

JAEA-Data/Code 2008-033- Setting saturated soil water content.
Include files: None.
Variables:- NX:

Integer Horizontal x-grid size.

- NY:

Integer Horizontal y-grid size.

- NS:

Integer Vertical grid size in soil.

- TEXTURE(NX,NY,0:NS):

Integer Soil type [1 to 24 USGS soil type].

- HWS(NX,NY,0:NS):

Real*8 Saturated soil water content [m 3 m -3].

SLIQU (Z, TEXTURE, MYU, DELT, ET, PR0, HW, EB, ER, T, RES, VFS, WF, WFR,IFLGFEED) (sliqu.f)Function:- Calculating water transport in soil.
Include files: Inclnum, prm_soil.
Variables:- Z(0:NS+1):

Real*8 Vertical grid in soil [m].

- TEXTURE(NX,NY,0:NS):

Integer Soil type [1 to 24 USGS soil type].

- MYU:

Real*8 Parameter for numerical scheme.

- DELT:

Real*8 Time increment [s].

- ET(NX,NY,0:NS):

Real*8 Root uptake (transpiration) rate [kg m -3 s -1 ].

- PR0(NX,NY):

Real*8 Precipitation intensity at soil surface [kg m -2 s -1].

- HW(NX,NY,0:NS):

Real*8 Soil water content [m 3 m-3].

- EB(NX,NY,0:NS):

Real*8 Evaporation rate of soil water [kg m-3 s-1 ].

- ER(NX,NY,0:NS):

Real*8 Runoff [kg m-2 s-1 (surface), kg m -3 s-1 (in-soil)].

- T(NX,NY,0:NS):

Real*8 Soil temperature [K].

- RES(NX,NY):

Real*8 Soil surface water [kg m -2 ].

- VFS(NX,NY,0:NS):

Real*8 Vapor flux in soil [kg m -2 s -1].

- WF(NX,NY,0:NS):

Real*8 Water flux in soil [kg m -2 s -1].

- WFR(NX,NY):

Real*8 Soil bottom water up-flow [kg m -2 s-1 ].

- IFLGFEED:

Integer Hydrology model coupling flag [1: on, 0: off].

SOLVE2 (NX, NY, NZ, A, B, C, D, Q, WK1, WK2) (solver2.f)Function:- Numerical scheme for diffusion equation in soil with variable bottom boundary.
Include files: None.
Variables:- NX:

Integer Horizontal x-grid size.

- NY:

Integer Horizontal y-grid size.
－ 51 －

JAEA-Data/Code 2008-033- NZ:

Integer Vertical grid size.

- A(NX,NY,0:NZ):

Real*8 Coefficient (7-8) in the numerical scheme.

- B(NX,NY,0:NZ):

Real*8 Coefficient (7-9) in the numerical scheme.

- C(NX,NY,0:NZ):

Real*8 Coefficient (7-10) in the numerical scheme.

- D(NX,NY,0:NZ):

Real*8 Coefficient (7-11) in the numerical scheme.

- Q(NX,NY,0:NZ):

Real*8 Solution.

- WK1(NX,NY,0:NZ):

Real*8 Temporal variable.

- WK2(NX,NY,0:NZ):

Real*8 Temporal variable.

EWPARA (NX, NY, NS, TEXTURE, HW, DW, DDW, KW, DKW) (ewpara.f)
Function:- Calculating diffusivity and hydraulic conductivity of soil water.
Include files: prm_soil.
Variables:- NX:

Integer Horizontal x-grid size.

- NY:

Integer Horizontal y-grid size.

- NS:

Integer Vertical grid size in soil.

- TEXTURE(NX,NY,0:NS):

Integer Soil type [1 to 24 USGS soil type].

- HW(NX,NY,0:NS):

Real*8 Soil water content [m 3 m-3].

- DW(NX,NY,NS):

Real*8 Soil water diffusivity [m 2 s -1].

- DDW(NX,NY,NS):

Real*8 Derivative of soil water diffusivity [m s -1].

- KW(NX,NY,NS+1):

Real*8 Hydraulic conductivity of soil [m s-1 ].

- DKW(NX,NY,NS+1):

Real*8 Derivative of hydraulic conductivity of soil [s -1].

STEMP (Z, TEXTURE, MYU, DELT, ZENITH, RS, RSD, RSS, RL, TR, UR, QR, CH0,
AIRDEN, PR0, TP, HTOPF, HRADF, HTOPL, HGTOP, HW, EB, T, RES, VFS,ALBDD, ALBDS, EMIS) (stemp.f)
Function:- Calculating heat conduction in soil.
Include files: Inclnum, prm_soil.
Variables:- Z(0:NS+1):

Real*8 Vertical grid in soil [m].

- TEXTURE(NX,NY,0:NS):

Integer Soil type [1 to 24 USGS soil type].

- MYU:

Real*8 Parameter for numerical scheme.

- DELT:

Real*8 Time increment [s].

- ZENITH(NX,NY):

Real*8 Solar zenith angle [deg.].

- RS(NX,NY):

Real*8 Soil surface solar radiation flux [W m -2 ].

- RSD(2,NX,NY):

Real*8 Soil surface direct solar radiation flux [W m -2].
－ 52 －

JAEA-Data/Code 2008-033- RSS(2,NX,NY):

Real*8 Soil surface diffuse solar radiation flux [W m-2 ].

- RL(NX,NY):

Real*8 Soil surface long-wane radiation flux [W m -2].

- TR(NX,NY):

Real*8 Lowest atmosphere layer temperature [K].

- UR(NX,NY):

Real*8 Lowest atmosphere layer wind speed [m s -1 ].

- QR(NX,NY):

Real*8 Lowest atmosphere layer humidity [kg kg -1 ].

- CH0(NX,NY):

Real*8 Heat exchange coefficient at soil surface.

- AIRDEN(NX,NY):

Real*8 Air density [kg m-3].

- PR0(NX,NY):

Real*8 Precipitation intensity at soil surface [kg m -2 s -1].

- TP(NX,NY):

Real*8 Precipitation temperature [K].

- HTOPF(NX,NY):

Real*8 Sensible heat flux at soil surface[W m -2].

- HRADF(NX,NY):

Real*8 Net radiation at soil surface [W m -2].

- HTOPL(NX,NY):

Real*8 Latent heat flux at soil surface [W m -2 ].

- HGTOP(NX,NY):

Real*8 Soil heat flux [W m -2 ].

- HW(NX,NY,0:NS):

Real*8 Soil water content [m 3 m-3].

- EB(NX,NY,0:NS):

Real*8 Evaporation rate of soil water [kg m-3 s-1 ].

- T(NX,NY,0:NS):

Real*8 Soil temperature [K].

- RES(NX,NY):

Real*8 Soil surface water [kg m -2 ].

- VFS(NX,NY,0:NS):

Real*8 Vapor flux in soil [kg m -2 s -1].

- ALBDD(2,NX,NY):

Real*8 Soil surface albedo for direct solar radiation.

- ALBDS(2,NX,NY):

Real*8 Soil surface albedo for diffuse solar radiation.

- EMIS(NX,NY):

Real*8 Soil surface emissivity.

ESPARA (NX, NY, NS, TEXTURE, HW, T, RES, CSRS, KS, WORK) (espara.f)
Function:- Calculating heat capacity and thermal conductivity of soil.
Include files: prm_soil.
Variables:- NX:

Integer Horizontal x-grid size.

- NY:

Integer Horizontal y-grid size.

- NS:

Integer Vertical grid size in soil.

- TEXTURE(NX,NY,0:NS):

Integer Soil type [1 to 24 USGS soil type].

- HW(NX,NY,0:NS):

Real*8 Soil water content [m 3 m-3].

- T(NX,NY,0:NS):

Real*8 Soil temperature [K].

- RES(NX,NY):

Real*8 Soil surface water [kg m -2 ].

- CSRS(NX,NY,0:NS):

Real*8 Heat capacity of bulk soil [J m-3 K-1 ].

- KS(NX,NY,NS):

Real*8 Thermal conductivity of soil [m2 s -1 ].

- WORK(NX,NY,0:NS):

Real*8 Thermal conductivity of soil [m2 s -1 ].

－ 53 －

JAEA-Data/Code 2008-033

EALBED (NX, NY, NS, TEXTURE, HW, ZENITH, ALBDD, ALBDS, EMIS) (efalbedo.f)
Function:- Calculating soil surface albedo and emissivity.
Include files: None.
Variables:- NX:

Integer Horizontal x-grid size.

- NY:

Integer Horizontal y-grid size.

- NS:

Integer Vertical grid size in soil.

- TEXTURE(NX,NY,0:NS):

Integer Soil type [1 to 24 USGS soil type].

- HW(NX,NY,0:NS):

Real*8 Soil water content [m 3 m-3].

- ZENITH(NX,NY):

Real*8 Solar zenith angle [deg.].

- ALBDD(2,NX,NY):

Real*8 Soil surface albedo for direct solar radiation.

- ALBDS(2,NX,NY):

Real*8 Soil surface albedo for diffuse solar radiation.

- EMIS(NX,NY):

Real*8 Soil surface emissivity.

SLCO2 (Z, TEXTURE, MYU, DELT, FRZ, SRZ, IFDBG, PHPS, TR, UR, QR, CE0,
AIRDEN, PR0, HW, TS, QS, ET, EB, ER, RES, WF, CO2BOT, CO2R, SS0, SP0,
AAZ, CO2S, FCO2S, ETCO2, SSCO2, RTCO2, ACO2G, ACO2A, DCO2A,
FSCO2) (slco2.f)Function:- Calculating CO2 transport in soil.
Include files: Inclnum.
Variables:- Z(0:NS+1):

Real*8 Vertical grid in soil [m].

- TEXTURE(NX,NY,0:NS):

Integer Soil type [1 to 24 USGS soil type].

- MYU:

Real*8 Parameter for numerical scheme.

- DELT:

Real*8 Time increment [s].

- FRZ(NX,NY,0:NS):

Real*8 Fractional root distribution [0 to 1 fraction].

- SRZ(NX,NY):

Real*8 Accumulated leaf area density (LAI) [m2 m -2 ].

- IFDBG

Integer Debug output flag [1: on, 0: off].

- PHPS(NX,NY):

Real*8 Surface air pressure [hPa].

- TR(NX,NY):

Real*8 Lowest atmosphere layer temperature [K].

- UR(NX,NY):

Real*8 Lowest atmosphere layer wind speed [m s -1 ].

- QR(NX,NY):

Real*8 Lowest atmosphere layer humidity [kg kg -1 ].

- CE0(NX,NY):

Real*8 Vapor exchange coefficient at soil surface.

- AIRDEN(NX,NY):

Real*8 Air density [kg m-3].

- PR0(NX,NY):

Real*8 Precipitation intensity at soil surface [kg m -2 s -1].

- HW(NX,NY,0:NS):

Real*8 Soil water content [m 3 m-3].
－ 54 －

JAEA-Data/Code 2008-033- TS(NX,NY,0:NS):

Real*8 Soil temperature [K].

- QS(NX,NY,0:NS):

Real*8 Specific humidity of soil pore [kg kg-1].

- ET(NX,NY,0:NS):

Real*8 Root uptake (transpiration) rate [kg m -3 s -1 ].

- EB(NX,NY,0:NS):

Real*8 Evaporation rate of soil water [kg m-3 s-1 ].

- ER(NX,NY,0:NS):

Real*8 Runoff [kg m-2 s-1 (surface), kg m -3 s-1 (in-soil)].

- RES(NX,NY):

Real*8 Soil surface water [kg m -2 ].

- WF(NX,NY,0:NS):

Real*8 Water flux in soil [kg m -2 s -1].

- CO2BOT(NX,NY):

Real*8 Soil bottom CO2 concentration [ppmv].

- CO2R(NX,NY):

Real*8 Surface air CO 2 concentration [ppmv].

- SS0:

Real*8 Optimal soil CO2 production rate [μmol m -2 s -1 ].

- SP0:

Real*8 Optimal root CO2 production rate [μmol m -2 s -1 ].

- AAZ:

Real*8 CO2 production distribution parameter in thesoil [m-1 ].

- CO2S(NX,NY,0:NS):

Real*8 Soil CO2 concentration [ppmv].

- FCO2S(NX,NY,NS):

Real*8 Soil CO2 flux [ppmv m -2 s-1 ].

- ETCO2(NX,NY,0:NS):

Real*8 Root up-take CO2 [ppmv s-1 ].

- SSCO2(NX,NY,0:NS):

Real*8 Soil respiration rate [ppmv s-1 ].

- RTCO2(NX,NY,0:NS):

Real*8 Root respiration rate [ppmv s-1].

- ACO2G(NX,NY,0:NS):

Real*8 Gas phase CO2 concentration [ppmv].

- ACO2A(NX,NY,0:NS):

Real*8 Aqueous phase CO2 concentration [ppmv].

- DCO2A(NX,NY):

Real*8 Drain rate of CO 2 [ppmv m s -1 ].

- FSCO2(NX,NY):

Real*8 Soil surface CO2 flux [ppmv m s -1 ].

HIFI1(QN, Q, QQ0, QQ1, WI, DELT, DZ) (shifi1.f)
Function:- One-dimensional HIFI scheme58) for advection term of CO2 transport in soil.
Include files: Inclnum.
Variables:- QN(NX,NY,0:NS):

Real*8 Change in Q during DELT.

- Q(NX,NY,0:NS):

Real*8 Dependent variable at the time step.

- QQ0(NX,NY):

Real*8 Top boundary value of Q.

- QQ1(NX,NY):

Real*8 Bottom boundary value of Q.

- WI(NX,NY,0:NS):

Real*8 Advection velocity [m s-1 ].

- DELT:

Real*8 Time increment.

- DZ(0:NS):

Real*8 Vertical grid [m].

RSCO2 (NX, NY, NS, NC, ZS, ZC, TEXTURE, RZ, AZ, PARABS, HW, PHPS, ROU,
CHANGEE, VMAX25, SLOPE, UU, VV, TT, QQ, CO2, TC, ES, ET, ETCO2,－ 55 －

JAEA-Data/Code 2008-033

RS, AS, AF, AG, AP, IFDBG, TIMET) (svrsco2.f)Function:- Calculating CO2 assimilation by photosynthesis and stomatal resistance.
Include files: None.
Variables:- NX:

Integer Horizontal x-grid size.

- NY:

Integer Horizontal y-grid size.

- NS:

Integer Vertical grid size in soil.

- NC:

Integer Vertical grid size in vegetation canopy.

- ZS(0:NS+1):

Real*8 Vertical grid in soil [m].

- ZC(NX,NY,0:NC+1):

Real*8 3-D grid height in vegetation canopy [m].

- TEXTURE(NX,NY,0:NS):

Integer Soil type [1 to 24 USGS soil type].

- RZ(NX,NY,0:NC,0:NS):

Real*8 Root distribution [0 to 1 fraction].

- AZ(NX,NY,0:NC):

Real*8 Leaf area density [m2 m -3 ].

- PARABS(2,NX,NY,NC):

Real*8 PAR absorbed by leaves [W m-2].

- HW(NX,NY,0:NS):

Real*8 Soil water content [m 3 m-3].

- PHPS(NX,NY):

Real*8 Surface air pressure [hPa].

- ROU(NX,NY):

Real*8 Air density [kg m-3].

- CHANGEE(NX,NY,0:NC):

Real*8 Vapor exchange coefficient.

- VMAX25(NX,NY,0:NC):

Real*8 Vm at 25 °C [μmol m -2 s-1].

- SLOPE(NX,NY,0:NC):

Real*8 Empirical constant m for stomatal resistance.

- UU(NX,NY,NC):

Real*8 Wind u-component [m s -1].

- VV(NX,NY,NC):

Real*8 Wind v-component [m s-1].

- TT(NX,NY,0:NC):

Real*8 Temperature of canopy air [°C].

- QQ(NX,NY,NC):

Real*8 Specific humidity of canopy air [kg kg-1].

- CO2(NX,NY,NC):

Real*8 CO2 concentration [ppmv].

- TC(2,NX,NY,0:NC):

Real*8 Temperature of sun-lit/shaded leaves [°C].

- ES(NX,NY,NC):

Real*8 Transpiration rate [kg m-2 s -1].

- ET(NX,NY,0:NS):

Real*8 Root uptake (transpiration) rate [kg m -3 s -1 ].

- ETCO2(NX,NY,0:NS):

Real*8 Root up-take CO2 [μmol m-3 s -1 or ppmv s -1 ].

- RS(2,NX,NY,0:NC):

Real*8 Stomatal resistance [s m-1].

- AS(2,NX,NY,0:NC):

Real*8 Net CO2 assimilation rate 1 [μmol m -2 s -1].

- AF(2,NX,NY,0:NC):

Real*8 Total CO2 flux from leaf [μmol m-2 s -1 ].

- AG(2,NX,NY,0:NC):

Real*8 Gross CO2 assimilation rate [μmol m -2 s-1 ].

- AP(2,NX,NY,0:NC):

Real*8 Net CO2 assimilation rate 2 [μmol m -2 s -1].

- IFDBG

Integer Debug output flag [1: on, 0: off].

- TIMET

Real*8 Integration time [h].

－ 56 －

JAEA-Data/Code 2008-033

RESISTS (NX, NY, NS, NC, ZS, TEXTURE, RZ, RSSTMIN, PARCRT, SMX, DDN,

SDN,

HW, QA, TA, PHPS, RRF, PARABS, SUNABS, FLIT, RS, FMX, FTX) (svrsst.f)
Function:- Calculating stomatal resistance by Jarvis-type scheme.
Include files: None.
Variables:- NX:

Integer Horizontal x-grid size.

- NY:

Integer Horizontal y-grid size.

- NS:

Integer Vertical grid size in soil.

- NC:

Integer Vertical grid size in vegetation canopy.

- ZS(0:NS+1):

Real*8 Vertical grid in soil [m].

- TEXTURE(NX,NY,0:NS):

Integer Soil type [1 to 24 USGS soil type].

- RZ(NX,NY,0:NC,0:NS):

Real*8 Root distribution [0 to 1 fraction].

- RSSTMIN(NX,NY,0:NC):

Real*8 Minimum stomatal resistance [s m -1].

- PARCRT(NX,NY,0:NC):

Real*8 Critical value of PAR [W m -2 ].

- SMX:

Real*8 Maximum solar radiation flux [W m -2].

- DDN(2,NX,NY,0:NC):

Real*8 Downward direct solar radiation flux [W m -2].

- SDN(2,NX,NY,0:NC):

Real*8 Downward diffuse solar radiation flux [W m -2].

- HW(NX,NY,0:NS):

Real*8 Soil water content [m 3 m-3].

- QA(NX,NY,NC):

Real*8 Specific humidity of canopy air [kg kg-1].

- TA(NX,NY,0:NC):

Real*8 Temperature of canopy air [°C].

- PHPS(NX,NY):

Real*8 Surface air pressure [hPa].

- RRF(2,NX,NY,NC):

Real*8 Fractionoftranspirationbysun-lit/shaded

leaves.
- PARABS(2,NX,NY,NC):

Real*8 PAR absorbed by leaves [W m-2].

- SUNABS(2,NX,NY,NC):

Real*8 Solar radiation absorbed by leaves [W m-2].

- FLIT(NX,NY,NC):

Real*8 Fraction of sun-lit leaves [0 to 1 fraction].

- RS(2,NX,NY,0:NC):

Real*8 Stomatal resistance [s m-1].

- FMX(2,NX,NY,0:NC):

Real*8 Stress function by humidity.

- FTX(2,NX,NY,0:NC):

Real*8 Stress function by temperature.

VLIQU (ZC, AZ, DELT, CHANGEE, WDSAT, WDEVP, FANGLE, FGRVEGE, ROU, PHPS,
UUC, VVC, TTC, QQC, WLC, PRT, TP, RP, TC, WDL, RS, FLIT, PD, RRD, RRS, RRF, EPR,
ECL, ECP, PR, WDLN, VEGTYPE) (svliqu.f)Function:- Calculating leaf surface water and water flux in canopy.
Include files: Inclnum.
Variables:－ 57 －

JAEA-Data/Code 2008-033- ZC(NX,NY,0:NC+1):

Real*8 3-D grid height in vegetation canopy [m].

- AZ(NX,NY,0:NC):

Real*8 Leaf area density [m 2 m -3 ].

- DELT:

Real*8 Time increment [s].

- CHANGEE(NX,NY,0:NC):

Real*8 Vapor exchange coefficient.

- WDSAT(NX,NY,0:NC):

Real*8 Maximum leaf surface water [kg m -2 ].

- WDEVP(NX,NY,0:NC):

Real*8 Leaf surface water for max. evaporation [kg m -2 ].

- FANGLE(NX,NY,0:NC):

Real*8 Mean zenith angle of leaf orientation [deg.].

- FGRVEGE(NX,NY,0:NC):

Real*8 Green leaf fraction [0 to 1 fraction].

- ROU(NX,NY):

Real*8 Air density [kg m-3].

- PHPS(NX,NY):

Real*8 Surface air pressure [hPa].

- UUC(NX,NY,NC):

Real*8 Wind u-component [m s -1].

- VVC(NX,NY,NC):

Real*8 Wind v-component [m s-1].

- TTC(NX,NY,0:NC):

Real*8 Temperature of canopy air [°C].

- QQC(NX,NY,NC):

Real*8 Specific humidity of canopy air [kg kg-1].

- WLC(NX,NY,NC):

Real*8 Fog water in canopy air [kg kg-1 ].

- PRT(NX,NY):

Real*8 Precipitation intensity [kg m-2 s -1].

- TP(NX,NY):

Real*8 Precipitation temperature [K].

- RP(NX,NY):

Real*8 Radius of rain drop [m].

- TC(2,NX,NY,0:NC):

Real*8 Temperature of sun-lit/shaded leaves [°C].

- WDL(NX,NY,0:NC):

Real*8 Leaf surface water [kg m -2].

- RS(2,NX,NY,0:NC):

Real*8 Stomatal resistance [s m-1].

- FLIT(NX,NY,NC):

Real*8 Fraction of sun-lit leaves [0 to 1 fraction].

- PD(NX,NY,NC):

Real*8 Drip from leaf surface water [kg m -2 s-1].

- RRD(2,NX,NY,NC):

Real*8 Total conductance for transpiration [m s -1 ].

- RRS(2,NX,NY,NC):

Real*8 Total conductance for evaporation [m s-1 ].

- RRF(2,NX,NY,NC):

Real*8 Fractionoftranspirationbysun-lit/shaded

leaves.
- EPR(NX,NY,NC):

Real*8 Evaporation from precipitation [kg m-3 s -1 ].

- ECL(NX,NY,NC):

Real*8 Accretion of fog by precipitation [kg m -3 s-1].

- ECP(NX,NY,NC):

Real*8 Accretion of fog by leaf [kg m -2 s -1].

- PR(NX,NY,0:NC):

Real*8 Vertical water flux in canopy [kg m -2 s-1 ].

- WDLN(NX,NY,0:NC):

Real*8 Leaf surface water for next time step [kg m -2].

- VEGTYPE(NX,NY,0:NC):

Integer Vegetation type in zvege.table_1D.

FOGCAP (NL, DELT, FRW, ROU, WS, WL, ECP, ECL, FRV, VTYPE, AZ) (svfogcp.f)
Function:- Calculating cloud water collection rate by leaves.
Include files: Inclcon2, prm_fog.
－ 58 －

JAEA-Data/Code 2008-033

Variables:- NL:

Integer Canopy layer number.

- DELT:

Real*8 Time increment [s].

- FRW:

Real*8 the product of projection coefficient for verticaldirection and horizontal wind speed [m s-1 ].

- ROU:

Real*8 Air density [kg m-3].

- WS:

Real*8 Wind speed at the layer of NL [m s -1 ].

- WL:

Real*8 Liquid water in canopy air at the layer of NL [kgm-3].

- ECP:

Real*8 Accretion of fog by leaf at the layer of NL [kg m -2
s-1 ].

- ECL:

Real*8 Accretion of fog by precipitation at the layer of
NL [kg m-3 s-1 ].

- FRV:

Real*8 Projection coefficient for vertical direction.

- VTYPE:

Integer Vegetation type in zvege.table_1D at the layer of
NL.

- AZ:

Real*8 Leaf area density at the layer of NL [m 2 m -3 ].

VTEMP (ZC, AZ, DELT, CHANGEH, FANGLE, REFLECT, TRANSMT, EMISSIV,

ROU,

PHPS, UUC, VVC, TTC, QQC, WLC, PR, TP, RDDN, RSDN, RSUP, RLDN, RLUP, RRD,
RRS, ED, ES, WDL, ZENITH, FLIT, TET, HS, HC, TC, ESL, ESS, ASCO2, EAS, RES)
(svtemp.f)
Function:- Calculating leaf surface heat budget.
Include files: Inclnum.
Variables:- ZC(NX,NY,0:NC+1):

Real*8 3-D grid height in vegetation canopy [m].

- AZ(NX,NY,0:NC):

Real*8 Leaf area density [m2 m -3 ].

- DELT:

Real*8 Time increment [s].

- CHANGEH(NX,NY,0:NC):

Real*8 Heat exchange coefficient.

- FANGLE(NX,NY,0:NC):

Real*8 Mean zenith angle of leaf orientation [deg.].

- REFLECT(2,NX,NY,0:NC): Real*8 Leaf reflectivity for visible and near-infrared.
- TRANSMT(2,NX,NY,0:NC): Real*8 Leaf transmissivity for visible and near-infrared.
- EMISSIV(NX,NY,0:NC):

Real*8 Leaf emissivity.

- ROU(NX,NY):

Real*8 Air density [kg m-3].

- PHPS(NX,NY):

Real*8 Surface air pressure [hPa].

- UUC(NX,NY,NC):

Real*8 Wind u-component [m s -1].

- VVC(NX,NY,NC):

Real*8 Wind v-component [m s-1].
－ 59 －

JAEA-Data/Code 2008-033- TTC(NX,NY,0:NC):

Real*8 Temperature of canopy air [°C].

- QQC(NX,NY,NC):

Real*8 Specific humidity of canopy air [kg kg-1].

- WLC(NX,NY,NC):

Real*8 Fog water in canopy air [kg kg-1 ].

- PR(NX,NY,0:NC):

Real*8 Vertical water flux in canopy [kg m -2 s-1 ].

- TP(NX,NY):

Real*8 Precipitation temperature [K].

- RDDN(2,NX,NY,0:NC):

Real*8 Downward direct solar radiation flux [W m -2].

- RSDN(2,NX,NY,0:NC):

Real*8 Downward diffuse solar radiation flux [W m -2].

- RSUP(2,NX,NY,0:NC):

Real*8 Upward diffuse solar radiation flux [W m-2].

- RLDN(2,NX,NY,0:NC):

Real*8 Downward long-wave radiation flux [W m -2].

- RLUP(2,NX,NY,0:NC):

Real*8 Upward long-wave radiation flux [W m -2].

- RRD(2,NX,NY,NC):

Real*8 Total conductance for transpiration [m s -1 ].

- RRS(2,NX,NY,NC):

Real*8 Total conductance for evaporation [m s-1 ].

- ED(NX,NY,NC):

Real*8 Evaporation rate of leaf surface water [kgm-2

s-1 ].
- ES(NX,NY,NC):

Real*8 Transpiration rate [kg m-2 s -1].

- WDL(NX,NY,0:NC):

Real*8 Leaf surface water [kg m -2].

- ZENITH(NX,NY):

Real*8 Solar zenith angle [deg.].

- FLIT(NX,NY,NC):

Real*8 Fraction of sun-lit leaves [0 to 1 fraction].

- TET(NX,NY,0:NC):

Real*8 Temperature of root up-take water [K].

- HS(NX,NY,0:NC):

Real*8 Cooling by precipitation [W m-2 ].

- HC(NX,NY,0:NC):

Real*8 Heat exchange between leaf and air [W m -2 ].

- TC(2,NX,NY,0:NC):

Real*8 Temperature of sun-lit/shaded leaves [°C].

- ESL(NX,NY,NC):

Real*8 Transpiration of sun-lit leaves [kg m-2 s -1 ].

- ESS(NX,NY,NC):

Real*8 Transpiration of sun-shaded leaves [kg m -2 s -1].

- ASCO2(2,NX,NY,0:NC):

Real*8 Net CO2 assimilation rate [μmol m -2 s-1 ].

- EAS(NX,NY,0:NC):

Real*8 Energy consumption by photosynthesis [W m -2].

- RES(NX,NY):

Real*8 Soil surface water [kg m -2 ].

GENER (I1, J1, ID, DK, A) (pgener.f)
Function:- Setting variables of numerical scheme for diffusion equation in atmosphere.
Include files: Inclnum and Inclcon1.
Variables:- I1:

Integer Horizontal x-grid number.

- J1:Integer Horizontal y-grid number.

- ID:

Integer Grid type.

- DK(M1):

Real*8 Diffusivity K in (7-1).

- A(N1):

Real*8 Coefficient A in (7-1).
－ 60 －

JAEA-Data/Code 2008-033

DIREC2 (N, TOP, PE1, PF1, Q) (pgener.f)
Function:- Numerical scheme for diffusion equation in atmosphere with variable top boundary.
Include files: Inclnum, Inclcon1.
Variables:- N:

Integer Vertical grid size.

- TOP:

Real*8 Gradient of variable at the top boundary.

- PE1:

Real*8 Surface boundary condition.

- PF1:

Real*8 Surface boundary condition.

- Q(M1):

Real*8 Solution.

DIREC1 (N, TOP, PE1, PF1, Q) (pgener.f)
Function:- Numerical scheme for diffusion equation in atmosphere with constant top boundary.
Include files: Inclnum, Inclcon1.
Variables:- N:

Integer Vertical grid size.

- TOP:

Real*8 Gradient of variable at the top boundary.

- PE1:

Real*8 Surface boundary condition.

- PF1:

Real*8 Surface boundary condition.

- Q(M1):

Real*8 Solution.

SED (NL, ROU, WL, SEDWL) (svsed.f)
Function:- Calculating gravitational flux of cloud water.
Include files: Inclcon2, prm_fog.
Variables:- NL:

Integer Atmospheric layer number.

- ROU:

Real*8 Air density [kg m-3].

- WL:

Real*8 Liquid water in canopy air at the layer of NL [kgm-3].

- SEDWL:

Real*8 Gravitational flux of cloud water [kg m -2 s-1 ].

ELI2VA(TA, QA, PHPA, RO, CP, CL, QASS, EFDT) (eli2va.f)
Function:- Calculating temperature and specific humidity after evaporation.
Include files: None.
－ 61 －

JAEA-Data/Code 2008-033

Variables:- TA:

Real*8 Air temperature before evaporation [K].

- QA :

Real*8 Specific humidity before evaporation [kg kg-1].

- PHPA:

Real*8 Air pressure [hPa].

- RO:

Real*8 Air density [kg m-3].

- CP:

Real*8 Specific heat of air [J K -1 kg -1 ].

- CL:

Real*8 Latent heat of vaporization [J kg -1].

- QASS:

Real*8 Saturated specific humidity [kg kg-1].

- EFDT:

Real*8 Potential evaporation [kg m -3 ].

A-4 Functions
SHMD (TEMP, PRESSURE, JFLG) (faipsy.f)
Input:

Output:- TEMP:

Real*8 Temperature [K].

- PRESSURE:

Real*8 Air pressure [hPa].

- JFLG:

Integer Flag [1: liquid water, 2: ice].

- SHMD:

Real*8 Saturated specific humidity [kg kg-1].

PSYM (Z, STB, ZZER) (faipsy.f)
Input:

Output:- Z:

Real*8 Height [m].

- STB:

Real*8 Stability length scale [m].

- ZZER:

Real*8 Roughness height of wind [m].

- PSYM: Real*8 Integrated shear function of wind.

PSYH (Z, STB, ZZER) (faipsy.f)
Input:- Z:

Real*8 Height [m].

- STB:

Real*8 Stability length scale [m].

- ZZER:

Real*8 Roughness height of heat [m].

- PSYH:

Real*8 Integrated shear function of temperature.

Input:- T:

Real*8 Temperature [K].

Output:- FL:

Real*8 Latent heat of vaporization [J kg -1].

Output:FL (T) (fl.f)

FCPAIR (Q) (fcpair.f)
Input:- Q:

Real*8 Specific humidity [kg kg-1 ].

Output:- FCPAIR:

Real*8 Specific heat of air [J kg -1 K-1 ].
－ 62 －

JAEA-Data/Code 2008-033

FDENSA (D, Q, T, QQ, TT) (fdensa.f)
Input:- D:

Real*8 Air density [kg m-3].

- Q:

Real*8 Specific humidity of air [kg kg -1 ].

- T:

Real*8 Air temperature [K].

- QQ:

Real*8 Specific humidity of soil pore [kg kg-1].

- TT:

Real*8 Soil temperature [K].

- FDENSA:

Real*8 Density of soil pore air [kg m-3].

Input:- T:

Real*8 Temperature [K].

Output:- FCPAIR:

Real*8 Specific heat of water [J kg -1 K-1 ].

Output:FCW (T) (fcw.f)

POTEV (T, BG, GA, QR, BRN, DLT, QSAT) (swadsp.f)
Input:- T:

Real*8 Temperature [K].

- BG:

Real*8 Soil heat flux at soil surface [W m -2].

- GA:

Real*8 Transfer coefficient for sensible heat flux at soil
surface.

Output:- QR:

Real*8 Specific humidity of air [kg kg -1 ].

- BRN:

Real*8 Net radiation at soil surface [W m -2].

- DLT:

Real*8 Coefficient for potential evaporation calculation.

- QSAT:

Real*8 Saturated specific humidity [kg kg-1].

- POTEV:

Real*8 Potential evaporation at soil surface [kg m -2 s-1 ].

EPS (DP, WS, VTYPE) (feps.f)
Input:
Output:- DP:

Real*8 Droplet diameter of cloud water [m].

- WS:

Real*8 Wind speed [m s-1 ].

- EPS:

Real*8 Capture efficiency of cloud water droplets.

FDIST (DDP, WL) (feps.f)
Input:
Output:- DDP:

Real*8 Droplet diameter of cloud water [m].

- WL:

Real*8 Liquid water in canopy air [kg kg -1].

- FDIST:

Real*8 Probability density of cloud water droplet.

DGL10 (a, b, AD, BD, PD, QD) (feps.f)
Input:- a:

Real*8 Constant for Gauss-Legendre integration.

- b:

Real*8 Constant for Gauss-Legendre integration.

- AD:

Real*8 Constant in Deirmendjian’s function.
－ 63 －

JAEA-Data/Code 2008-033

Output:- BD:

Real*8 Constant in Deirmendjian’s function.

- PD:

Real*8 Constant in Deirmendjian’s function.

- QD:

Real*8 Constant in Deirmendjian’s function.

- DGL10:

Real*8 Integrated value of probability density of cloudwater droplet.

DEIR (DP, AD, BD, PD, QD) (feps.f)
Input:

Output:- DP:

Real*8 Mean droplet diameter of cloud water [m].

- AD:

Real*8 Constant in Deirmendjian’s function.

- BD:

Real*8 Constant in Deirmendjian’s function.

- PD:

Real*8 Constant in Deirmendjian’s function.

- QD:

Real*8 Constant in Deirmendjian’s function.

- DEIR:

Real*8 Probability density of cloud water droplet for
Deirmendjian parameterization.

VGRV (X) (svfopcp.f, svsed.f)
Input:- X:

Real*8 Droplet diameter [m].

Output:- VGRV:

Real*8 Terminal settling velocity of cloud droplet [m s -1].

SHMDD (TEMP, PRESSURE, JFLG) (faipsy.f)
Input:

Output:- TEMP:

Real*8 Temperature [K].

- PRESSURE:

Real*8 Air pressure [hPa].

- JFLG:

Integer Flag [1: liquid water, 2: ice].

- SHMDD:

Real*8 Derivative of saturated humidity [kg kg -1 K-1].

FAIM (Z, STB) (faipsy.f)
Input:
Output:- Z:

Real*8 Height [m].

- STB:

Real*8 Stability length scale [m].

- FAIM:

Real*8 Non-dimensional shear function of wind.

－ 64 －