<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">%表題  Model description of AGCM5 of GFD-Dennou-Club edition
%
%履歴  1999-12-01 石渡正樹 環境研モノグラフ  
%                 熱帯に局在する暖水域が降水分布に与える影響
%                 のモデル記述部分として作成
%      2000-01-24 林祥介
%      2006-05-01 石渡正樹 AGCM5 モデルドキュメントとして再構成.

%\section{Model description}
%\Dseclab{モデル}

   AGCM5 of the GFD-DENNOU CLUB edition
        is a three-dimensional primitive system on a 
        sphere (Swamp Project, 1998). 
   The model is basically the same as that used in Numaguti (1993), 
        and its details are described in Numaguti(1992). 
   In the following subsections, we briefly describe the
        model system and its implementation. 

\section{Basic equations}

   The dynamical part of the model consists of 
        the vorticity, divergence,
        hydrostatic, 
        humidity and thermodynamic equations as follows. 
   \begin{eqnarray}
   \label{zeta}
        &amp; &amp;          \DD{\zeta}{t} 
        =  \frac{1}{a \cos \varphi} 
        \DP{}{\lambda} 
        \{ - (\zeta + f) u - \dot{\sigma} \DP{u}{\sigma} 
        - \frac{RT'}{a \cos \varphi} \DP{\pi}{\lambda} \} \nonumber \\
        &amp; &amp; \qquad    - \frac{1}{a \cos \varphi} 
        \DP{}{\varphi} 
        [ \{ (\zeta + f) v - \dot{\sigma} \DP{u}{\sigma} 
        -  \frac{RT'}{a} \DP{\pi}{\varphi} \} \cos \varphi]
        - F^{diff}_{\zeta}, \\
   \label{divergence}
        &amp;  &amp;        \DD{D}{t} 
        =   \frac{1}{a \cos \varphi}  \DP{}{\lambda}
        \{ (\zeta + f) v - \dot{\sigma} \DP{u}{\sigma} 
        - x\frac{RT'}{a} \DP{\pi}{\varphi} \}   \nonumber \\
        &amp; &amp; \qquad +  \frac{1}{a \cos \varphi}  \DP{}{\varphi}
        [  \{ (\zeta + f) v - \dot{\sigma} \DP{u}{\sigma} 
        - \frac{RT'}{a} \DP{\pi}{\varphi} \} \cos \varphi]  \nonumber \\
        &amp; &amp; \qquad - \nabla^2 (\Phi  + R \overline{T} \pi + E )
        - F^{diff}_{D}, \\
   \label{質量}
        &amp;   &amp;    \DD{\pi}{t} 
        =  - \frac{1}{a \cos \varphi} \DP{u}{\lambda}
        - \frac{1}{a \cos \varphi} \DP{}{\varphi} ( v \cos \varphi)
        - \frac{\partial \dot{\sigma}}{\partial \sigma}, \\
   \label{静水圧}
        &amp; &amp;  \frac{\partial \Phi}{\partial \sigma} = - \frac{RT}{\sigma}, \\
   \label{水蒸気}
        &amp; &amp; \DD{q}{t}
        =  \frac{g}{p_s} \DP{F^{vdf}_q}{\sigma}
        + F^{vdf}_q
        + S^{cond}_q, \\
   \label{熱力}
        &amp;   &amp; \DD{T}{t} 
        =   \frac{R T}{c_p}
        \left\{   \DP{\pi}{t} 
          + \frac{u}{a \cos \varphi} \DP{\pi}{\lambda}
          + \frac{v}{a} \DP{\pi}{\varphi}
          + \frac{\dot{\sigma}}{\sigma} \right\}
        \nonumber  \\ 
        &amp; &amp; \qquad \qquad \mbox{}   + \frac{1}{c_p} \left( 
          \frac{g}{p_s} \DP{F^{vdf}_T}{\sigma}
          + \frac{g}{p_s} \DP{F^{vdf}_{rad}}{\sigma} \right)
        +  F^{diff}_{T} + L S^{cond}_q.
   \end{eqnarray}
   In these equations,
   \begin{eqnarray}
        \DD{}{t} &amp; \equiv &amp;  \DP{}{t} 
        + \frac{u}{a \cos \varphi} \DP{}{\lambda}
        + \frac{v}{a} \DP{}{\varphi}
        + \dot{\sigma} \DP{}{\sigma}, \\
        \zeta &amp; \equiv &amp;   \frac{1}{a \cos \varphi} \DP{v}{\lambda}
        - \frac{1}{a \cos \varphi} 
        \DP{}{\varphi} (u \cos \varphi),  \\
        D    &amp; \equiv &amp; \frac{1}{a \cos \varphi} \DP{u}{\lambda}
        + \frac{1}{a \cos \varphi} 
        \DP{}{\varphi} (v \cos \varphi), \\
        \Phi &amp; \equiv &amp; gz, \\
        \pi  &amp; \equiv &amp; \ln p_s, 
   \end{eqnarray}
   where,
   $(\lambda, \varphi)$ are the latitudinal and longitudinal coordinates;
   $\sigma$  is the vertical coordinate; 
   $u$, $v$ are the horizontal wind components;
   $\dot{\sigma}$ is the vertical wind in $\sigma$ coordinate;
   $T$ is temperature;
   $q$ is specific humidity;
   $p_s$ is surface pressure;
   $T_g$ is surface temperature;
   $\Phi$ is geopotential; 
   $f$ is the Coriolis parameter;
   $a$ is the radius of the earth; 
   $g$ is the acceleration of gravity;
   $R$ is the gas constant of dry air;
   $c_p$ is the atmospheric specific heat at constant pressure;
   $L$ is the latent heat of water;
   $F^{diff}_{\zeta}$, $F^{diff}_{D}$, $F^{diff}_{T}$
   $F^{diff}_{q}$ are the horizontal diffusion terms;
   $F^{vdf}_u$, $F^{vdf}_v$, $F^{vdf}_T$, 
   $F^{vdf}_q$ are the vertical diffusion terms; and
   $S^{cond}_q$ is the condensation source (loss) of the specific
   humidity.
    
   The horizontal diffusion terms are represented by a
        hyper-viscosity type formula,
   \begin{equation}
        F^{diff} = - (-1)^n K \nabla^{2n}
   \end{equation}
   for which $n=8$ is adopted. 
   The value of $K$ is determined such that the e-folding time of
   the component with wavenumber 42 is four hours. 

\section{Discretization}

   The equations are discretized horizontally 
        by a non-alias latitude-longitude grid 
        and are integrated by the pseudo-spectral method 
        with a triangular truncation.
%        at wavenumber 42 (T42),
%        i.e., . 
   The scheme of Arakawa and Suarez (1983) is utilized for 
        discretizing in the vertical direction.
%   The number of vertical levels is 16.
%   The positions of the vertical grid points are listed 
%        in \Dtabref{鉛直グリッド}. 
   The time integration is performed by the use of 
        a leap-frog scheme for the terms of dynamical processes
        other than dissipative processes, 
        a backward scheme for the terms of physical processes, 
        and an Euler scheme for the horizontal diffusion terms.
%   The timestep utilized is 20 minutes.
   The time filter of Asselin (1972) is used to remove 
        the computational mode caused by the leap-frog scheme. 

%\vspace{1em}
%\begin{table}[h]
%      \caption{The vertical levels}
%      \Dtablab{鉛直グリッド}
%      \begin{center}
%        \begin{tabular}{ccccccccccccccc} \hline
%          $k$ &amp;  &amp; 1 &amp;  &amp; 2 &amp;  &amp; 3 &amp;  &amp; 4 &amp;  &amp; 5 &amp;  &amp; 6 \\
%          $\sigma_k$ &amp; 1.0 &amp; .995 &amp; .99 &amp; .98 &amp; .97 &amp; .95 &amp; .93 &amp;
%                .90 &amp; .87 &amp; .83 &amp; .79 &amp; .75 \\
%          \hline
%          $k$ &amp;  &amp; 7 &amp;  &amp; 8 &amp;  &amp; 9 &amp;  &amp;10 &amp;  &amp;11 &amp;  &amp;12 \\
%          $\sigma_k$ &amp; .7  &amp; .65  &amp; .6  &amp; .55 &amp;  .5 &amp; .45 &amp; .41 &amp; 
%                .37 &amp; .33 &amp; .29 &amp; .26 &amp; .23 \\
%          \hline
%          $k$ &amp;  &amp;13 &amp;  &amp;14 &amp;  &amp;15 &amp;  &amp;16 &amp; \\
%          $\sigma_k$ &amp; .2 &amp; .17 &amp; .15 &amp; .12 &amp; .1 &amp; .074 &amp; .05 &amp; .021 &amp; .0\\
%          \hline
%        \end{tabular}
%      \end{center}
%\end{table}

\section{Physical processes}

   The model contains the following simplified physical processes.

\subsection{Cumulus parameterization}

   The cumulus parameterization schemes implemented are
      the adjustment scheme (Manabe {\em et al.}, 1965) and Kuo scheme
      (Kuo, 1974). 
   In the followings, the description of adjustment schmes is given.
   Condensed water is removed from the system 
      immediately after condensation and hence there are 
      neither rain drops nor clouds.

   Moist convective adjustment occurs when the following conditions
      are satisfied: 
   \begin{itemize}
       \item In two neighboring levels, the atmosphere becomes moist
             unstable, that is , the temperature lapse rate
             exceeds the moist adiabatic lapse rate. 
       \item Both of the neighboring levels are saturated or
             oversaturated. 
   \end{itemize}
   In the process of adjustment, the temperature and specific humidity of 
      the two neighboring levels are instantaneously adjusted to 
      values so that the levels are saturated and the 
      temperature lapse rate between them coincides with the moist
      adiabatic lapse rate. 
   In due course, moist static energy has to be conserved:
   \begin{equation}
        \Deqlab{ds 保存}
        \sum^{k}_{l=k-1} (C_{p}       T_l + L       q_l) \Delta p_l = 
        \sum^{k}_{l=k-1} (C_{p} \hat{T}_l + L \hat{q}_l) \Delta p_l,
   \end{equation}
   where $k$ is the index of the vertical level. 
   Here, $\hat{T}, \hat{q}$ denote the values before adjustment,
   while $T, q$ denote the values after the adjustment. 

   The moist unstable condition for layer $(k,k-1)$ is evaluated
   by 
   \begin{eqnarray}
          ST_{k-1/2}  &amp; \equiv &amp; 
          \hat{T}_k - \hat{T}_{k-1}
       +  \frac{L}{C_p} 
          \left( q^*(\hat{T}_{k},p_{k}) - q^*(\hat{T}_{k-1},p_{k-1}) \right)
         - \frac{R}{C_p} 
         \frac{\Delta p_{k-1/2}}{p_{k-1/2}} 
         \frac{\hat{T}_k + \hat{T}_{k-1}}{2}
                             \nonumber \\
        &amp; &gt; &amp; 0.
   \end{eqnarray}
   The saturation condition of two neighboring levels is simply
   expressed by 
   \begin{eqnarray}
        \hat{q}_{k-1} &gt; q^*(\hat{T}_{k-1},p_{k-1}), \ \ 
        \hat{q}_{k}   &gt; q^*(\hat{T}_{k}  ,p_{k}). 
   \end{eqnarray}

   The adjusted values $(T_k, q_k)$ are determined as 
   \begin{eqnarray}
         T_{k-1} &amp; = &amp; \hat{T}_{k-1} + \Delta T_{k-1}, \\
         T_{k}   &amp; = &amp; \hat{T}_{k}   + \Delta T_{k},  \\
         q_{k-1} &amp; = &amp; q^* (\hat{T}_{k-1},p_{k-1}) 
                       + \left. \DP{q^*}{T} \right|_{k-1} 
                         \Delta T_{k-1}, \\
         q_{k} &amp; = &amp; q^* (\hat{T}_{k},p_{k}) 
                       + \left. \DP{q^*}{T} \right|_{k} 
                         \Delta T_{k},
   \end{eqnarray}
   \begin{equation}
            \Delta T_{k} 
          =  \frac{
                   (1 + \gamma_{k-1}) \Delta p_{k-1} S_{k-1/2}
                 + [   1 + \gamma_{k-1} 
                     - \kappa \frac{\Delta p_{k-1/2}}{2 p_{k-1/2}}]
                   \frac{L}{C_p} \Delta \hat{Q}
             }{     (1 + \gamma_{k-1}) (1 + \gamma_{k}) 
                  (\Delta p_{k-1} + \Delta p_{k})
                + \kappa \frac{\Delta p_{k-1/2}}{2 p_{k-1/2}}
                  [   (1 + \gamma_{k-1}) \Delta p_{k-1} 
                    - (1 + \gamma_{k}) \Delta p_{k}]
            } \raisebox{-5mm}{,}
   \end{equation}
   \begin{equation}
            \Delta T_{k-1} 
         =  - \frac{ (1 + \gamma_{k}) \Delta p_{k} 
             }{
                  (1 + \gamma_{k-1}) \Delta p_{k-1}
             } \Delta T_{k} 
            + \frac{ 1 }{ (1 + \gamma_{k-1}) \Delta p_{k-1} }
              \frac{L}{C_p}
              \Delta T_{k},
   \end{equation}
   where
   \begin{equation}
          \gamma_{k}   = \frac{L}{C_p} \left. \DP{q^*}{T} \right|_{k}, 
          \ \ 
          \gamma_{k-1} = \frac{L}{C_p} \left. \DP{q^*}{T} \right|_{k-1}, 
       \end{equation}  
       \begin{equation}
            \Delta \hat{Q} 
          =    (\hat{q}_{k-1} - q^*(\hat{T}_{k-1},p_{k-1}) ) \Delta p_{k-1}
            +  (\hat{q}_{k} - q^*(\hat{T}_{k},p_{k}) ) \Delta p_{k}.
       \end{equation}

   Precipitation of one column is calculated as
       \begin{equation}
          P = \sum^{K}_{k=1} (\hat{q}_k - q_k) \frac{\Delta p_k}{g}.
       \end{equation}          
   
   These adjustment procedures are applied to all levels
        until the vertical temperature and moisture profiles 
        converge.


\subsection{Large scale condensation}

   The large scale condensation scheme is applied at each grid point
        when the value of the specific humidity exceeds 
        its saturation value. 
   The large scale condensation scheme demands that
        the value of the specific humidity be adjusted to 
        the saturation value. 
   Condensed water is assumed to fall out as rain and is removed
        from the system immediately.  
   There are no rain drops or clouds, as in the case of 
        the convective adjustment scheme. 
   In the time integration procedure, 
        the large scale condensation scheme is applied just after
        the convective adjustment scheme. 

   The condition for large scale condensation can be denoted by
   \begin{equation}
       \hat{q} &gt; q^*(\hat{T},p),
   \end{equation}
   where $\hat{T}$, $\hat{q}$ are the values of temperature and specific
   humidity before the adjustment. 
   After the large scale condensation adjustment, 
   the values of the temperature and specific humidity become
     \begin{equation}
        q = q^*(T,p)
     \end{equation}
   and
     \begin{equation}
      C_p T + Lq = C_p \hat{T} + L \hat{q}.
     \end{equation}
   
   The solutions to these are obtained by Newton method;
     \begin{eqnarray}
        T &amp; = &amp; 
                \hat{T}
              + \frac{L}{C_p} 
           \frac{\hat{q} - q^* (\hat{T},p)}
                {1 + \frac{L}{C_p} \DP{q^*}{T}}, \\
        q &amp; = &amp; 
                \hat{q}
              - \frac{L}{C_p} \DP{q^*}{T}
           \frac{\hat{q} - q^* (\hat{T},p)}
                {1 + \frac{L}{C_p} \DP{q^*}{T}}.
     \end{eqnarray}
   The amount of precipitation is calculated by
      \begin{equation}
         P = \frac{\Delta p}{g} (\hat{q} - q).
      \end{equation}

\subsection{Radiation}

   Radiative flux is calclated by band model.
   The effects of scattering of solar radiation is not considered.
   It is assumed that there are no cloud nor rain drops. 

   As for longwave radiation, 
        we include only the absorption and emission 
        of water and dry air; 
        the effect of scattering is not included. 
   The absorption coefficients of the longwave radiation are 
      chosen 
      so that the cooling profile of the atmosphere 
      roughly resembles the observed one.

   Radiation flux $F_{R,i}$ is
       \begin{equation}
          F_R(z) =   \left( \pi B(T_g) - \pi B(T_s) \right) {\cal T}^f(z,0)
                   +        \pi B(T(z_T))                  {\cal T}^f(z,z_T)
                   - \int_0^{z_T} \DD{\pi B}{\xi} {\cal T}^f(z,\xi) d \xi, 
       \end{equation}
   where ${\cal T}^f(z_1,z_2)$ is the flux transfer function
        between $z=z_1,z_2$ and 
        $\pi B \equiv \sigma_{SB} T^4$ is the source function.

   The flux transfer function ${\cal T}^f(z_1,z_2)$ is assumed to be
   given as 
       \begin{equation}
          {\cal T}^f(z_1,z_2) 
         = \sum_{i=1}^{N_R} b_i 
          \exp \left( - \delta_R | \tau_{R,i}(z_1) - \tau_{R,i}(z_2) | 
            \right).
       \end{equation}
   Here, $\tau_i(z)$ is the optical depth whose value at the top of
       the atmosphere is 0 and 
       \begin{equation}
          \tau_{R,i}(z) = \int_z^\infty k_{R,i} \rho q dz 
              + \int_z^\infty \bar{k}_{R,i} \rho dz,
       \end{equation}
       where $k_{R,i}$ and $\bar{k}_{R,i}$ are
       the absorption coefficients of water vapor
       and dry air for wavenumber band  $i$, respectively. 
   The values of these absorption coefficients are assumed to be
       constant. 
   $b_i$ is the ratio of the energy for wavenumber band $i$ 
        to the total radiative energy and is 
        also assumed to be constant.
   The default value for $\delta_R$ is 1.5.
   In most cases, $N_R$ is set to be 4, three representing water
       vapor bands and one representing dry air band. 
   The default values of $k_{R,i}$, $\bar{k}_{R,i}$ and $b_i$ 
        are summarized in \Dtabref{放射定数}. 


\vspace*{1em}
\begin{table}[h]
      \caption{The values of constants used
                       in the longwave radiation scheme.}
      \Dtablab{放射定数}
      \begin{center}
        \begin{tabular}{lllll} \hline
          Band Number &amp; 1 &amp; 2 &amp; 3 &amp; 4\\
          $k_{R,i}$   &amp; 8.0   &amp; 1.0   &amp; 0.1   &amp; 0.0 \\
          $\overline{k}_{R,i}$   &amp; 0.0   &amp; 0.0   &amp; 0.0   &amp; $5.0 \times 10^{-5}$ \\
          $b_i$       &amp; 0.2  &amp; 0.1 &amp;  0.1 &amp; 0.6 \\
          \hline
        \end{tabular}
      \end{center}
\end{table}


\subsection{Vertical mixing}

   Vertical turbulent mixing is represented by 
      the level II scheme of Mellor and Yamada (1974).  

   In this scheme of Mellor and Yamada (1974), 
       diffusion coefficients $K_M$, $K_T$, and $K_q$ depends on 
       stability and are given as functions of $u$, $v$, and $T$:
       \begin{eqnarray}
         K_M       &amp; = &amp; l^2 \frac{\Delta |\Dvect{v}|}{\Delta z} S_M,  \\
         K_T = K_q &amp; = &amp; l^2 \frac{\Delta |\Dvect{v}|}{\Delta z} S_H.
       \end{eqnarray}
   Following Blakadar(1962), the mixing length $l$ is evaluated by
       \begin{equation}
         l = \frac{\kappa z}{1+\kappa z/l_0},
       \end{equation}
       where $l_0=300$ m and $\kappa=0.4$ is the K\'{a}rman constant.

   In these expressions, 
        $S_M$ and $S_H$ are given as functions of the Richardson
        number. 
       \begin{eqnarray}
            S_M &amp; = &amp; B_1^{1/2} ( 1- R_f )^{1/2} 
                    \widetilde{S_M}^{3/2},  \\
            S_H &amp; = &amp; B_1^{1/2} ( 1- R_f )^{1/2} 
                    \widetilde{S_M}^{1/2} \widetilde{S_H}, 
       \end{eqnarray} 
   where 
       \begin{equation}
          \widetilde{S_H} = \frac{ \alpha_1-\alpha_2 R_f }{ 1-R_f },
       \end{equation}
       \begin{equation}
          \widetilde{S_M} = \frac{ \beta_1-\beta_2 R_f }
                                 { \beta_3-\beta_4 R_f } 
                  \widetilde{S_H}. 
       \end{equation}

   $R_f$ is the flux Richardson number, given as
       \begin{equation}
           R_f = \frac{1}{2 \beta_2}
                 \left[ \beta_1 + \beta_4 R_{iB}
              - \sqrt{ ( \beta_1 + \beta_4 R_{iB} )^2 
                       - 4 \beta_2 \beta_3 R_{iB} }
              \right] ,
       \end{equation}
       where $R_{iB}$ is the bulk Richardson number and
       \begin{eqnarray}
          \alpha_1 &amp; = &amp; 3 A_2 \gamma_1,  \\
          \alpha_2 &amp; = &amp; 3 A_2 (\gamma_1+\gamma_2), \\
          \beta_1  &amp; = &amp; A_1 B_1 ( \gamma_1 - C_1 ), \\
          \beta_2  &amp; = &amp; A_1 [ B_1 ( \gamma_1 - C_1 ) + 6 A_1 + 3 A_2 ], \\
          \beta_3  &amp; = &amp; A_2 B_1, \\
          \beta_4  &amp; = &amp; A_2 [ B_1 ( \gamma_1 + \gamma_2 ) - 3 A_1 ],
       \end{eqnarray}
       \begin{equation}
          (A_1, B_1, A_2, B_2, C_1 ) = ( 0.92, 16.6, 0.74, 10.1, 0.08 ) ,
       \end{equation}
       \begin{equation}
          \gamma_1 = \frac{1}{3} - \frac{2 A_1}{B_1}\, , \, \, \, 
          \gamma_2 = \frac{B_2}{B_1} + 6\frac{A_1}{B_1}.
       \end{equation}

   $R_{iB}$ is defined as
       \begin{equation}
          R_{iB} = \frac{\displaystyle 
               \frac{q}{\theta_s} \frac{\Delta \theta}{\Delta z} }
              {\displaystyle
                  \left( \frac{\Delta u}{\Delta z} \right)^2 
                + \left( \frac{\Delta v}{\Delta z} \right)^2
               }. 
       \end{equation}
   Here, the value of $(\Delta \bullet)$ at the vertical
        level ${k-1/2}$ is evaluated as $\bullet_{k}
        - \bullet_{k-1}$.


\subsection{Surface flux}

   The surface fluxes of momentum, heat and water vapor are
      evaluated by the bulk formulae
      \begin{eqnarray}
        (F_{\lambda})_{1/2} &amp; = &amp; - \rho_s C_{M} | \Dvect{v}_1 | u_1, \\
        (F_{\varphi})_{1/2} &amp; = &amp; - \rho_s C_{M} | \Dvect{v}_1 | u_1, \\
        (F_{T})_{1/2} &amp; = &amp;  \rho_s C_{p} C_{H} | \Dvect{v}_1 | 
                             (T_g - (\frac{T}{\Pi}_1)), \\
        (F_{q})_{1/2} &amp; = &amp;  \rho_s \beta C_{E} | \Dvect{v}_1 | 
                             (q^*(T_g, p_s) - q_1), 
      \end{eqnarray}
      where the suffix $1/2$ indicates that these are the value at 
      the vertical level $1/2$. 
   
   Following Louis {\it et al.} (1982), the bulk coefficients 
        $C_M$, $C_H$, and $C_E$ are determined by 
       \begin{eqnarray}
        C_M = \left\{ 
               \begin{array}{ll}
               C_{M0} \left( 1 + 10 R_{iB}/\sqrt{1 + 5 R_{iB} } \right)^{-1} 
                          &amp; R_{iB} \geq 0, \\
               C_{M0} \left[1 - 10 R_{iB} 
                          \left( 1 + 75 C_{M0}
                                     \sqrt{\frac{z_1}{z_{M0}} |R_{iB}| }
                          \right)^{-1}
                      \right]
                          &amp; R_{iB} &lt; 0,
              \end{array} \right.
       \end{eqnarray}
       \begin{eqnarray}
        C_H = C_E = \left\{ 
               \begin{array}{ll}
               C_{H0} \left( 1 + 15 R_{iB}/\sqrt{1 + 5 R_{iB} } \right)^{-1} 
                          &amp; R_{iB} \geq 0, \\
               C_{H0} \left[1 - 15 R_{iB} 
                          \left( 1 + 75 C_{H0}
                                     \sqrt{\frac{z_1}{z_{M0}} |R_{iB}| }
                          \right)^{-1}
                      \right]
                          &amp; R_{iB} &lt; 0,
              \end{array} \right.
       \end{eqnarray}
       where $R_{iB}$ is the bulk Richardson number evaluated by
       \begin{equation}
          R_{iB} = \frac{\frac{g}{\theta_s} (\theta_1 - \theta_g)/g}
                        {(u_1/z_1)^2}
                 = \frac{g}{\theta_s}
                   \frac{T_1 (p_s/p_1)^{\kappa} - T_g}
                        {u^2_1/z_1}.
       \end{equation}
       $C_{M0}$ and $C_{H0}$ are the bulk coefficients for neutral case:
       \begin{eqnarray}
          C_{M0} &amp; = &amp; \left( 
                           \frac{\kappa}
                                {\ln \left( \frac{z_1}{z_{M0}} \right)}
                       \right),        \\ 
          C_{H0} &amp; = &amp; \left( 
                           \frac{\kappa}
                                {\ln \left( \frac{z_1}{z_{H0}} \right)}
                       \right),
       \end{eqnarray} 
       where $z_{M0}$ and $z_{H0}$ are the roughness parameters of
       the surface. 
       The values are $10^{-4}$m. 

</pre></body></html>