manual/s_phys_pkgs/diagnostics.tex

\section{Diagnostics--A Flexible Infrastructure}
\label{sec:pkg:diagnostics}
\begin{rawhtml}
<!-- CMIREDIR:package_diagnostics: -->
\end{rawhtml}

\subsection{Introduction}

\noindent
This section of the documentation describes the Diagnostics package available within 
the GCM.  A large selection of model diagnostics is available for output.  
In addition to the diagnostic quantities pre-defined in the GCM, there exists
the option, in any experiment, to define a new diagnostic quantity and include it
as part of the diagnostic output with the addition of a single subroutine call in the
routine where the field is computed. As a matter of philosophy, no diagnostic is enabled 
as default, thus each user must specify the exact diagnostic information required for an 
experiment.  This is accomplished by enabling the specific diagnostic of interest cataloged 
in the Diagnostic Menu (see Section \ref{sec:diagnostics:menu}). Instructions for enabling
diagnostic output and defining new diagnostic quantities are found in Section 
\ref{sec:diagnostics:usersguide} of this document.

\noindent
The Diagnostic Menu is a hard-wired enumeration of diagnostic quantities available within 
the GCM. Once a diagnostic is enabled, the GCM will continually increment an array
specifically allocated for that diagnostic whenever the appropriate quantity is computed.  
A counter is defined which records how many times each diagnostic quantity has been 
incremented.  Several special diagnostics are included in the menu. Quantities refered to 
as ``Counter Diagnostics'', are defined for selected diagnostics which record the 
frequency at which a diagnostic is incremented separately for each model grid location.
Quantitied refered to as ``User Diagnostics'' are included in the menu to facilitate
defining new diagnostics for a particular experiment.

\subsection{Equations}
Not relevant.

\subsection{Key Subroutines and Parameters}
\label{sec:diagnostics:diagover}

\noindent
The diagnostics are computed at various times and places within the GCM. Because the
MIT GCM may employ a staggered grid, diagnostics may be computed at grid box centers,
corners, or edges, and at the middle or edge in the vertical. Some diagnostics are scalars, 
while others are components of vectors. An internal array is defined which contains 
information concerning various grid attributes of each diagnostic. The GDIAG
array (in common block \\diagnostics in file diagnostics.h) is internally defined as a 
character*8 variable, and is equivalenced to a character*1 "parse" array in output in 
order to extract the grid-attribute information.  The GDIAG array is described in 
Table \ref{tab:diagnostics:gdiag.tabl}.

\begin{table}
\caption{Diagnostic Parsing Array}
\label{tab:diagnostics:gdiag.tabl}
\begin{center}
\begin{tabular}{ |c|c|l| }
\hline
\multicolumn{3}{|c|}{\bf Diagnostic Parsing Array} \\ 
\hline
\hline
Array & Value & Description \\
\hline
  parse(1)   & $\rightarrow$ S &  Scalar Diagnostic                 \\ 
             & $\rightarrow$ U &  U-vector component Diagnostic     \\ 
             & $\rightarrow$ V &  V-vector component Diagnostic     \\ \hline
  parse(2)   & $\rightarrow$ U &  C-Grid U-Point                    \\ 
             & $\rightarrow$ V &  C-Grid V-Point                    \\ 
             & $\rightarrow$ M &  C-Grid Mass Point                 \\ 
             & $\rightarrow$ Z &  C-Grid Vorticity (Corner) Point   \\ \hline
  parse(3)   & $\rightarrow$ R &  Not Currently in Use              \\ \hline
  parse(4)   & $\rightarrow$ P &  Positive Definite Diagnostic      \\ \hline
  parse(5)   & $\rightarrow$ C &  Counter Diagnostic                \\
             & $\rightarrow$ D &  Disabled Diagnostic for output    \\ \hline
  parse(6-8) & $\rightarrow$ C &  3-digit integer corresponding to  \\
             &                 &  vector or counter component mate  \\ \hline
\end{tabular}
\addcontentsline{lot}{section}{Table 3:  Diagnostic Parsing Array}
\end{center}
\end{table}


\noindent
As an example, consider a diagnostic whose associated GDIAG parameter is equal
to ``UU  002''.  From GDIAG we can determine that this diagnostic is a 
U-vector component located at the C-grid U-point.
Its corresponding V-component diagnostic is located in Diagnostic \# 002.


\noindent
In this way, each Diagnostic in the model has its attributes (ie. vector or scalar,
C-grid location, etc.) defined internally.  The Output routines use this information 
in order to determine what type of transformations need to be performed.  Any 
interpolations are done at the time of output rather than during each model step.
In this way the User has flexibility in determining the type of gridded data which 
is output.


\noindent
There are several utilities within the GCM available to users to enable, disable,
clear, write and retrieve model diagnostics, and may be called from any routine.  
The available utilities and the CALL sequences are listed below.


\noindent
{\bf fill\_diagnostics}:  This routine will increment the specified diagnostic
quantity with a field sent through the argument list.


\noindent
\begin{tabbing}
XXXXXXXXX\=XXXXXX\= \kill
\>        call fill\_diagnostics (myThid, chardiag, levflg, nlevs, \\
                bibjflg, bi, bj, arrayin) \\
\\
where \>  myThid   \>= Current Process(or) \\
      \>  chardiag \>= Character *8 expression for diag to fill \\
      \>  levflg   \>= Integer flag for vertical levels: \\
      \>           \> 0 indicates multiple levels incremented in qdiag \\
      \>           \> non-0 (any integer) - WHICH single level to increment. \\
      \>           \> negative integer - the input data array is single-leveled \\
      \>           \> positive integer - the input data array is multi-leveled \\
      \>  nlevs    \>= indicates Number of levels to be filled (1 if levflg <> 0) \\
      \>           \> positive: fill in "nlevs" levels in the same order as \\
      \>           \> the input array \\
      \>           \> negative: fill in -nlevs levels in reverse order. \\
      \>  bibjflg  \>= Integer flag to indicate instructions for bi bj loop \\
      \>           \> 0 indicates that the bi-bj loop must be done here \\
      \>           \> 1 indicates that the bi-bj loop is done OUTSIDE \\
      \>           \> 2 indicates that the bi-bj loop is done OUTSIDE \\
      \>           \>    AND that we have been sent a local array \\
      \>           \> 3 indicates that the bi-bj loop is done OUTSIDE \\
      \>           \>    AND that we have been sent a local array \\
      \>           \>    AND that the array has the shadow regions \\
      \>  bi       \>= X-direction process(or) number - used for bibjflg=1-3 \\
      \>  bj       \>= Y-direction process(or) number - used for bibjflg=1-3 \\
      \>  arrayin  \>= Field to increment diagnostics array \\
\end{tabbing}


\noindent
{\bf setdiag}:  This subroutine enables a diagnostic from the Diagnostic Menu, meaning 
that space is allocated for the diagnostic and the model routines will increment the 
diagnostic value during execution.  This routine is the underlying interface
between the user and the desired diagnostic.  The diagnostic is referenced by its diagnostic
number from the menu, and its calling sequence is given by:

\noindent
\begin{tabbing}
XXXXXXXXX\=XXXXXX\= \kill
\>        call setdiag (num) \\
\\
where \>  num   \>= Diagnostic number from menu \\
\end{tabbing}

\noindent
{\bf getdiag}:  This subroutine retrieves the value of a model diagnostic.  This routine 
is particulary useful when called from a user output routine, although it can be called 
from any routine.  This routine returns the time-averaged value of the diagnostic by
dividing the current accumulated diagnostic value by its corresponding counter.  This 
routine does not change the value of the diagnostic itself, that is, it does not replace 
the diagnostic with its time-average.  The calling sequence for this routine is givin by:

\noindent
\begin{tabbing}
XXXXXXXXX\=XXXXXX\= \kill
\>        call getdiag (lev,num,qtmp,undef) \\
\\
where \>  lev   \>= Model Level at which the diagnostic is desired \\
      \>  num   \>= Diagnostic number from menu \\
      \>  qtmp  \>= Time-Averaged Diagnostic Output \\
      \>  undef \>= Fill value to be used when diagnostic is undefined \\
\end{tabbing}

\noindent
{\bf clrdiag}:  This subroutine initializes the values of model diagnostics to zero, and is
particularly useful when called from user output routines to re-initialize diagnostics 
during the run.  The calling sequence is:

\noindent
\begin{tabbing}
XXXXXXXXX\=XXXXXX\= \kill
\>        call clrdiag (num) \\
\\
where \>  num   \>= Diagnostic number from menu \\
\end{tabbing}

\noindent
{\bf zapdiag}:  This entry into subroutine SETDIAG disables model diagnostics, meaning 
that the diagnostic is no longer available to the user.  The memory previously allocated 
to the diagnostic is released when ZAPDIAG is invoked.  The calling sequence is given by:

\noindent
\begin{tabbing}
XXXXXXXXX\=XXXXXX\= \kill
\>        call zapdiag (NUM) \\
\\
where \>  num   \>= Diagnostic number from menu \\
\end{tabbing}


\subsection{Usage Notes}
\label{sec:diagnostics:usersguide}

\noindent
We begin this section with a discussion on the manner in which computer 
memory is allocated for diagnostics. All GCM diagnostic quantities are stored in the 
single diagnostic array QDIAG which is located in the file \\
\filelink{pkg/diagnostics/diagnostics.h}{pkg-diagnostics-diagnostics.h}.
and has the form:

common /diagnostics/ qdiag(1-Olx,sNx+Olx,1-Olx,sNx+Olx,numdiags,Nsx,Nsy)

\noindent
where numdiags is an Integer variable which should be set equal to the number of 
enabled diagnostics, and qdiag is a three-dimensional array.  The first two-dimensions 
of qdiag correspond to the horizontal dimension of a given diagnostic, while the third 
dimension of qdiag is used to identify diagnostic fields and levels combined. In order 
to minimize the memory requirement of the model for diagnostics, the default GCM 
executable is compiled with room for only one horizontal diagnostic array, or with
numdiags set to 1. In order for the User to enable more than 1 two-dimensional diagnostic,
the size of the diagnostics common must be expanded to accomodate the desired diagnostics.
This can be accomplished by manually changing the parameter numdiags in the
file \filelink{pkg/diagnostics/diagnostics\_SIZE.h}{pkg-diagnostics-diagnostics_SIZE.h}.
numdiags should be set greater than or equal to the sum of all the diagnostics activated
for output each multiplied by the number of levels defined for that diagnostic quantity.
This is illustrated in the example below:

\noindent
To use the diagnostics package, other than enabling it in packages.conf
and turning the usediagnostics flag in data.pkg to .TRUE., a namelist
must be supplied in the run directory called data.diagnostics. The namelist
will activate a user-defined list of diagnostics quantities to be computed,
specify the frequency of output, the number of levels, and the name of
up to 10 separate output files. A sample data.diagnostics namelist file:

\noindent
$\#$ Diagnostic Package Choices \\
 $\&$diagnostics\_list \\
  frequency(1) = 10, \ \\
   levels(1,1) = 1.,2.,3.,4.,5., \ \\
   fields(1,1) = 'UVEL    ','VVEL    ', \ \\
   filename(1) = 'diagout1', \ \\
  frequency(2) = 100, \ \\
   levels(1,2) = 1.,2.,3.,4.,5., \ \\
   fields(1,2) = 'THETA   ','SALT    ', \ \\
   filename(2) = 'diagout2', \ \\
 $\&$end \ \\

\noindent
In this example, there are two output files that will be generated
for each tile and for each output time. The first set of output files
has the prefix diagout1, does time averaging every 10 time steps
(frequency is 10), they will write fields which are multiple-level 
fields and output levels 1-5. The names of diagnostics quantities are 
UVEL and VVEL.  The second set of output files
has the prefix diagout2, does time averaging every 100 time steps,
they include fields which are multiple-level fields, levels output are 1-5,
and the names of diagnostics quantities are THETA and SALT.

\noindent
In order to define and include as part of the diagnostic output any field
that is desired for a particular experiment, two steps must be taken. The
first is to enable the ``User Diagnostic'' in data.diagnostics. This is
accomplished by setting one of the fields slots to either UDIAG1 through 
UDIAG10, for multi-level fields, or SDIAG1 through SDIAG10 for single level
fields. These are listed in the diagnostics menu. The second step is to
add a call to fill\_diagnostics from the subroutine in which the quantity
desired for diagnostic output is computed. 

\newpage

\subsubsection{GCM Diagnostic Menu}
\label{sec:diagnostics:menu}

\begin{tabular}{lllll}
\hline\hline
N & NAME & UNITS & LEVELS & DESCRIPTION \\
\hline

&\\
1 & UFLUX    &   $Newton/m^2$  &    1  
         &\begin{minipage}[t]{3in}
          {Surface U-Wind Stress on the atmosphere}
         \end{minipage}\\
2 & VFLUX    &   $Newton/m^2$  &    1  
         &\begin{minipage}[t]{3in}
          {Surface V-Wind Stress on the atmosphere}
         \end{minipage}\\
3 & HFLUX    &   $Watts/m^2$  &    1  
         &\begin{minipage}[t]{3in}
          {Surface Flux of Sensible Heat}
         \end{minipage}\\
4 & EFLUX    &   $Watts/m^2$  &    1  
         &\begin{minipage}[t]{3in}
          {Surface Flux of Latent Heat}
         \end{minipage}\\
5 & QICE     &   $Watts/m^2$  &    1  
         &\begin{minipage}[t]{3in}
          {Heat Conduction through Sea-Ice}
         \end{minipage}\\
6 & RADLWG   &   $Watts/m^2$ &    1  
         &\begin{minipage}[t]{3in}
          {Net upward LW flux at the ground}
         \end{minipage}\\
7 & RADSWG   &   $Watts/m^2$  &    1 
         &\begin{minipage}[t]{3in}
          {Net downward SW flux at the ground} 
         \end{minipage}\\
8 & RI       &  $dimensionless$ &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Richardson Number}
         \end{minipage}\\
9 & CT       &  $dimensionless$ &  1 
         &\begin{minipage}[t]{3in}
          {Surface Drag coefficient for T and Q}
         \end{minipage}\\
10 & CU       & $dimensionless$ &  1 
         &\begin{minipage}[t]{3in}
          {Surface Drag coefficient for U and V}
         \end{minipage}\\
11 & ET       &  $m^2/sec$ &  Nrphys
         &\begin{minipage}[t]{3in}
          {Diffusivity coefficient for T and Q}
         \end{minipage}\\
12 & EU       &  $m^2/sec$ &  Nrphys
         &\begin{minipage}[t]{3in}
          {Diffusivity coefficient for U and V}
         \end{minipage}\\
13 & TURBU    &  $m/sec/day$ &  Nrphys 
         &\begin{minipage}[t]{3in}
          {U-Momentum Changes due to Turbulence}
         \end{minipage}\\
14 & TURBV    &  $m/sec/day$ &  Nrphys 
         &\begin{minipage}[t]{3in}
          {V-Momentum Changes due to Turbulence}
         \end{minipage}\\
15 & TURBT    &  $deg/day$ &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Temperature Changes due to Turbulence}
         \end{minipage}\\
16 & TURBQ    &  $g/kg/day$ &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Specific Humidity Changes due to Turbulence}
         \end{minipage}\\
17 & MOISTT   &   $deg/day$ &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Temperature Changes due to Moist Processes}
         \end{minipage}\\
18 & MOISTQ   &  $g/kg/day$ &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Specific Humidity Changes due to Moist Processes}
         \end{minipage}\\
19 & RADLW    &  $deg/day$ &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Net Longwave heating rate for each level}
         \end{minipage}\\
20 & RADSW    &  $deg/day$ &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Net Shortwave heating rate for each level}
         \end{minipage}\\
21 & PREACC   &  $mm/day$ &  1
         &\begin{minipage}[t]{3in}
          {Total Precipitation}
         \end{minipage}\\
22 & PRECON   &  $mm/day$ &  1
         &\begin{minipage}[t]{3in}
          {Convective Precipitation}
         \end{minipage}\\
23 & TUFLUX   &  $Newton/m^2$ &  Nrphys
         &\begin{minipage}[t]{3in}
          {Turbulent Flux of U-Momentum}
         \end{minipage}\\
24 & TVFLUX   &  $Newton/m^2$ &  Nrphys
         &\begin{minipage}[t]{3in}
          {Turbulent Flux of V-Momentum}
         \end{minipage}\\
25 & TTFLUX   &  $Watts/m^2$ &  Nrphys
         &\begin{minipage}[t]{3in}
          {Turbulent Flux of Sensible Heat}
         \end{minipage}\\
26 & TQFLUX   &  $Watts/m^2$ &  Nrphys
         &\begin{minipage}[t]{3in}
          {Turbulent Flux of Latent Heat}
         \end{minipage}\\
27 & CN       &  $dimensionless$ &  1
         &\begin{minipage}[t]{3in}
          {Neutral Drag Coefficient}
         \end{minipage}\\
28 & WINDS     &  $m/sec$ &  1
         &\begin{minipage}[t]{3in}
          {Surface Wind Speed}
         \end{minipage}\\
29 & DTSRF     &  $deg$ &  1
         &\begin{minipage}[t]{3in}
          {Air/Surface virtual temperature difference}
         \end{minipage}\\
30 & TG        &  $deg$ &  1
         &\begin{minipage}[t]{3in}
          {Ground temperature}
         \end{minipage}\\
31 & TS        &  $deg$ &  1
         &\begin{minipage}[t]{3in}
          {Surface air temperature (Adiabatic from lowest model layer)}
         \end{minipage}\\
32 & DTG       &  $deg$ &  1
         &\begin{minipage}[t]{3in}
          {Ground temperature adjustment}
         \end{minipage}\\

\end{tabular}

\newpage
\vspace*{\fill}
\begin{tabular}{lllll}
\hline\hline
N & NAME & UNITS & LEVELS & DESCRIPTION \\
\hline

&\\
33 & QG        &  $g/kg$ &  1
         &\begin{minipage}[t]{3in}
          {Ground specific humidity}
         \end{minipage}\\
34 & QS        &  $g/kg$ &  1
         &\begin{minipage}[t]{3in}
          {Saturation surface specific humidity}
         \end{minipage}\\

&\\
35 & TGRLW    &    $deg$   &    1  
         &\begin{minipage}[t]{3in}
          {Instantaneous ground temperature used as input to the
           Longwave radiation subroutine} 
         \end{minipage}\\
36 & ST4      &   $Watts/m^2$  &    1  
         &\begin{minipage}[t]{3in}
          {Upward Longwave flux at the ground ($\sigma T^4$)}
         \end{minipage}\\
37 & OLR      &   $Watts/m^2$  &    1  
         &\begin{minipage}[t]{3in}
          {Net upward Longwave flux at the top of the model}
         \end{minipage}\\
38 & OLRCLR   &   $Watts/m^2$  &    1  
         &\begin{minipage}[t]{3in}
          {Net upward clearsky Longwave flux at the top of the model}
         \end{minipage}\\
39 & LWGCLR   &   $Watts/m^2$  &    1  
         &\begin{minipage}[t]{3in}
          {Net upward clearsky Longwave flux at the ground}
         \end{minipage}\\
40 & LWCLR    &  $deg/day$ &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Net clearsky Longwave heating rate for each level}
         \end{minipage}\\
41 & TLW      &    $deg$   &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Instantaneous temperature used as input to the Longwave radiation
          subroutine} 
         \end{minipage}\\
42 & SHLW     &    $g/g$   &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Instantaneous specific humidity used as input to the Longwave radiation
          subroutine} 
         \end{minipage}\\
43 & OZLW     &    $g/g$   &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Instantaneous ozone used as input to the Longwave radiation
          subroutine} 
         \end{minipage}\\
44 & CLMOLW   &    $0-1$   &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Maximum overlap cloud fraction used in the Longwave radiation
          subroutine} 
         \end{minipage}\\
45 & CLDTOT   &    $0-1$   &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Total cloud fraction used in the Longwave and Shortwave radiation
          subroutines} 
         \end{minipage}\\
46 & RADSWT   &    $Watts/m^2$   &  1 
         &\begin{minipage}[t]{3in}
          {Incident Shortwave radiation at the top of the atmosphere}
         \end{minipage}\\
47 & CLROSW   &    $0-1$   &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Random overlap cloud fraction used in the shortwave radiation
          subroutine} 
         \end{minipage}\\
48 & CLMOSW   &    $0-1$   &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Maximum overlap cloud fraction used in the shortwave radiation
          subroutine} 
         \end{minipage}\\
49 & EVAP     &    $mm/day$   &  1 
         &\begin{minipage}[t]{3in}
          {Surface evaporation}
         \end{minipage}\\
\end{tabular}
\vfill

\newpage
\vspace*{\fill}
\begin{tabular}{lllll}
\hline\hline
N & NAME & UNITS & LEVELS & DESCRIPTION \\
\hline

&\\
50 & DUDT     &    $m/sec/day$ &  Nrphys
         &\begin{minipage}[t]{3in}
          {Total U-Wind tendency}
         \end{minipage}\\
51 & DVDT     &    $m/sec/day$ &  Nrphys
         &\begin{minipage}[t]{3in}
          {Total V-Wind tendency}
         \end{minipage}\\
52 & DTDT     &    $deg/day$ &  Nrphys
         &\begin{minipage}[t]{3in}
          {Total Temperature tendency}
         \end{minipage}\\
53 & DQDT     &    $g/kg/day$ &  Nrphys
         &\begin{minipage}[t]{3in}
          {Total Specific Humidity tendency}
         \end{minipage}\\
54 & USTAR    &    $m/sec$ &  1
         &\begin{minipage}[t]{3in}
          {Surface USTAR wind}
         \end{minipage}\\
55 & Z0       &    $m$ &  1
         &\begin{minipage}[t]{3in}
          {Surface roughness}
         \end{minipage}\\
56 & FRQTRB   &    $0-1$ &  Nrphys-1
         &\begin{minipage}[t]{3in}
          {Frequency of Turbulence}
         \end{minipage}\\
57 & PBL      &    $mb$ &  1
         &\begin{minipage}[t]{3in}
          {Planetary Boundary Layer depth}
         \end{minipage}\\
58 & SWCLR    &  $deg/day$ &  Nrphys 
         &\begin{minipage}[t]{3in}
          {Net clearsky Shortwave heating rate for each level}
         \end{minipage}\\
59 & OSR      &   $Watts/m^2$  &    1 
         &\begin{minipage}[t]{3in}
          {Net downward Shortwave flux at the top of the model}
         \end{minipage}\\
60 & OSRCLR   &   $Watts/m^2$  &    1  
         &\begin{minipage}[t]{3in}
          {Net downward clearsky Shortwave flux at the top of the model}
         \end{minipage}\\
61 & CLDMAS   &   $kg / m^2$  &    Nrphys
         &\begin{minipage}[t]{3in}
          {Convective cloud mass flux}
         \end{minipage}\\
62 & UAVE     &   $m/sec$  &    Nrphys
         &\begin{minipage}[t]{3in}
          {Time-averaged $u-Wind$}
         \end{minipage}\\
63 & VAVE     &   $m/sec$  &    Nrphys
         &\begin{minipage}[t]{3in}
          {Time-averaged $v-Wind$}
         \end{minipage}\\
64 & TAVE     &   $deg$  &    Nrphys
         &\begin{minipage}[t]{3in}
          {Time-averaged $Temperature$}
         \end{minipage}\\
65 & QAVE     &   $g/g$  &    Nrphys
         &\begin{minipage}[t]{3in}
          {Time-averaged $Specific \, \, Humidity$}
         \end{minipage}\\
66 & PAVE     &   $mb$  &    1
         &\begin{minipage}[t]{3in}
          {Time-averaged $p_{surf} - p_{top}$}
         \end{minipage}\\
67 & QQAVE    &   $(m/sec)^2$  &    Nrphys
         &\begin{minipage}[t]{3in}
          {Time-averaged $Turbulent Kinetic Energy$}
         \end{minipage}\\
68 & SWGCLR   &   $Watts/m^2$  &    1  
         &\begin{minipage}[t]{3in}
          {Net downward clearsky Shortwave flux at the ground} 
         \end{minipage}\\
69 & SDIAG1   &             &    1  
         &\begin{minipage}[t]{3in}
          {User-Defined Surface Diagnostic-1} 
         \end{minipage}\\
70 & SDIAG2   &             &    1  
         &\begin{minipage}[t]{3in}
          {User-Defined Surface Diagnostic-2} 
         \end{minipage}\\
71 & UDIAG1   &             &    Nrphys
         &\begin{minipage}[t]{3in}
          {User-Defined Upper-Air Diagnostic-1} 
         \end{minipage}\\
72 & UDIAG2   &             &    Nrphys
         &\begin{minipage}[t]{3in}
          {User-Defined Upper-Air Diagnostic-2} 
         \end{minipage}\\
73 & DIABU    & $m/sec/day$ &    Nrphys
         &\begin{minipage}[t]{3in}
          {Total Diabatic forcing on $u-Wind$} 
         \end{minipage}\\
74 & DIABV    & $m/sec/day$ &    Nrphys
         &\begin{minipage}[t]{3in}
          {Total Diabatic forcing on $v-Wind$} 
         \end{minipage}\\
75 & DIABT    & $deg/day$ &    Nrphys
         &\begin{minipage}[t]{3in}
          {Total Diabatic forcing on $Temperature$} 
         \end{minipage}\\
76 & DIABQ    & $g/kg/day$ &    Nrphys
         &\begin{minipage}[t]{3in}
          {Total Diabatic forcing on $Specific \, \, Humidity$} 
         \end{minipage}\\

\end{tabular}
\vfill

\newpage
\vspace*{\fill}
\begin{tabular}{lllll}
\hline\hline
N & NAME & UNITS & LEVELS & DESCRIPTION \\
\hline

77 & VINTUQ  & $m/sec \cdot g/kg$ &    1
         &\begin{minipage}[t]{3in}
          {Vertically integrated $u \, q$} 
         \end{minipage}\\
78 & VINTVQ  & $m/sec \cdot g/kg$ &    1
         &\begin{minipage}[t]{3in}
          {Vertically integrated $v \, q$} 
         \end{minipage}\\
79 & VINTUT  & $m/sec \cdot deg$ &    1
         &\begin{minipage}[t]{3in}
          {Vertically integrated $u \, T$} 
         \end{minipage}\\
80 & VINTVT  & $m/sec \cdot deg$ &    1
         &\begin{minipage}[t]{3in}
          {Vertically integrated $v \, T$} 
         \end{minipage}\\
81 & CLDFRC  & $0-1$ &    1
         &\begin{minipage}[t]{3in}
          {Total Cloud Fraction} 
         \end{minipage}\\
82 & QINT    & $gm/cm^2$ &    1
         &\begin{minipage}[t]{3in}
          {Precipitable water} 
         \end{minipage}\\
83 & U2M     & $m/sec$ &    1
         &\begin{minipage}[t]{3in}
          {U-Wind at 2 meters}
         \end{minipage}\\
84 & V2M     & $m/sec$ &    1
         &\begin{minipage}[t]{3in}
          {V-Wind at 2 meters}
         \end{minipage}\\
85 & T2M     & $deg$ &    1
         &\begin{minipage}[t]{3in}
          {Temperature at 2 meters}
         \end{minipage}\\
86 & Q2M     & $g/kg$ &    1
         &\begin{minipage}[t]{3in}
          {Specific Humidity at 2 meters}
         \end{minipage}\\
87 & U10M    & $m/sec$ &    1
         &\begin{minipage}[t]{3in}
          {U-Wind at 10 meters}
         \end{minipage}\\
88 & V10M    & $m/sec$ &    1
         &\begin{minipage}[t]{3in}
          {V-Wind at 10 meters}
         \end{minipage}\\
89 & T10M    & $deg$ &    1
         &\begin{minipage}[t]{3in}
          {Temperature at 10 meters}
         \end{minipage}\\
90 & Q10M    & $g/kg$ &    1
         &\begin{minipage}[t]{3in}
          {Specific Humidity at 10 meters}
         \end{minipage}\\
91 & DTRAIN  & $kg/m^2$ &    Nrphys
         &\begin{minipage}[t]{3in}
          {Detrainment Cloud Mass Flux}
         \end{minipage}\\
92 & QFILL   & $g/kg/day$ &    Nrphys
         &\begin{minipage}[t]{3in}
          {Filling of negative specific humidity}
         \end{minipage}\\

\end{tabular}
\vspace{1.5in}
\vfill

\newpage

\subsubsection{Diagnostic Description}

In this section we list and describe the diagnostic quantities available within the 
GCM.  The diagnostics are listed in the order that they appear in the 
Diagnostic Menu, Section \ref{sec:diagnostics:menu}.
In all cases, each diagnostic as currently archived on the output datasets
is time-averaged over its diagnostic output frequency:

\[
{\bf DIAGNOSTIC} = {1 \over TTOT} \sum_{t=1}^{t=TTOT} diag(t)
\]
where $TTOT = {{\bf NQDIAG} \over \Delta t}$, {\bf NQDIAG} is the 
output frequency of the diagnostic, and $\Delta t$ is
the timestep over which the diagnostic is updated.  

{\bf 1)  \underline {UFLUX} Surface Zonal Wind Stress on the Atmosphere ($Newton/m^2$) } 

The zonal wind stress is the turbulent flux of zonal momentum from 
the surface. See section 3.3 for a description of the surface layer parameterization.
\[
{\bf UFLUX} =  - \rho C_D W_s u \hspace{1cm}where: \hspace{.2cm}C_D = C^2_u
\]
where $\rho$ = the atmospheric density at the surface, $C_{D}$ is the surface
drag coefficient, $C_u$ is the dimensionless surface exchange coefficient for momentum 
(see diagnostic number 10), $W_s$ is the magnitude of the surface layer wind, and $u$ is 
the zonal wind in the lowest model layer.
\\


{\bf 2)  \underline {VFLUX} Surface Meridional Wind Stress on the Atmosphere ($Newton/m^2$) } 

The meridional wind stress is the turbulent flux of meridional momentum from 
the surface. See section 3.3 for a description of the surface layer parameterization.
\[
{\bf VFLUX} =  - \rho C_D W_s v \hspace{1cm}where: \hspace{.2cm}C_D = C^2_u
\]
where $\rho$ = the atmospheric density at the surface, $C_{D}$ is the surface
drag coefficient, $C_u$ is the dimensionless surface exchange coefficient for momentum 
(see diagnostic number 10), $W_s$ is the magnitude of the surface layer wind, and $v$ is 
the meridional wind in the lowest model layer.
\\

{\bf 3)  \underline {HFLUX} Surface Flux of Sensible Heat ($Watts/m^2$) } 

The turbulent flux of sensible heat from the surface to the atmosphere is a function of the
gradient of virtual potential temperature and the eddy exchange coefficient:
\[
{\bf HFLUX} =  P^{\kappa}\rho c_{p} C_{H} W_s (\theta_{surface} - \theta_{Nrphys})
\hspace{1cm}where: \hspace{.2cm}C_H = C_u C_t
\]
where $\rho$ = the atmospheric density at the surface, $c_{p}$ is the specific
heat of air, $C_{H}$ is the dimensionless surface heat transfer coefficient, $W_s$ is the 
magnitude of the surface layer wind, $C_u$ is the dimensionless surface exchange coefficient 
for momentum (see diagnostic number 10), $C_t$ is the dimensionless surface exchange coefficient 
for heat and moisture (see diagnostic number 9), and $\theta$ is the potential temperature 
at the surface and at the bottom model level.
\\


{\bf 4)  \underline {EFLUX} Surface Flux of Latent Heat ($Watts/m^2$) } 

The turbulent flux of latent heat from the surface to the atmosphere is a function of the
gradient of moisture, the potential evapotranspiration fraction and the eddy exchange coefficient:
\[
{\bf EFLUX} =  \rho \beta L C_{H} W_s (q_{surface} - q_{Nrphys})
\hspace{1cm}where: \hspace{.2cm}C_H = C_u C_t
\]
where $\rho$ = the atmospheric density at the surface, $\beta$ is the fraction of
the potential evapotranspiration actually evaporated, L is the latent
heat of evaporation, $C_{H}$ is the dimensionless surface heat transfer coefficient, $W_s$ is the 
magnitude of the surface layer wind, $C_u$ is the dimensionless surface exchange coefficient 
for momentum (see diagnostic number 10), $C_t$ is the dimensionless surface exchange coefficient 
for heat and moisture (see diagnostic number 9), and $q_{surface}$ and $q_{Nrphys}$ are the specific
humidity at the surface and at the bottom model level, respectively.
\\

{\bf 5)  \underline {QICE} Heat Conduction Through Sea Ice ($Watts/m^2$) } 

Over sea ice there is an additional source of energy at the surface due to the heat
conduction from the relatively warm ocean through the sea ice. The heat conduction
through sea ice represents an additional energy source term for the ground temperature equation.

\[
{\bf QICE} = {C_{ti} \over {H_i}} (T_i-T_g)
\]

where $C_{ti}$ is the thermal conductivity of ice, $H_i$ is the ice thickness, assumed to
be $3 \hspace{.1cm} m$ where sea ice is present, $T_i$ is 273 degrees Kelvin, and
$T_g$ is the temperature of the sea ice.

NOTE: QICE is not available through model version 5.3, but is available in subsequent versions.
\\
 

{\bf 6) \underline {RADLWG} Net upward Longwave Flux at the surface ($Watts/m^2$)}

\begin{eqnarray*}
{\bf RADLWG} & =  & F_{LW,Nrphys+1}^{Net} \\
             & =  & F_{LW,Nrphys+1}^\uparrow - F_{LW,Nrphys+1}^\downarrow
\end{eqnarray*}
\\
where Nrphys+1 indicates the lowest model edge-level, or $p = p_{surf}$.
$F_{LW}^\uparrow$ is
the upward Longwave flux and $F_{LW}^\downarrow$ is the downward Longwave flux.
\\

{\bf 7) \underline {RADSWG} Net downard shortwave Flux at the surface ($Watts/m^2$)}

\begin{eqnarray*}
{\bf RADSWG} & =  & F_{SW,Nrphys+1}^{Net} \\
             & =  & F_{SW,Nrphys+1}^\downarrow - F_{SW,Nrphys+1}^\uparrow
\end{eqnarray*}
\\
where Nrphys+1 indicates the lowest model edge-level, or $p = p_{surf}$.
$F_{SW}^\downarrow$ is
the downward Shortwave flux and $F_{SW}^\uparrow$ is the upward Shortwave flux.
\\


\noindent
{\bf 8)  \underline {RI} Richardson Number} ($dimensionless$)

\noindent
The non-dimensional stability indicator is the ratio of the buoyancy to the shear:
\[
{\bf RI} = { { {g \over \theta_v} \pp {\theta_v}{z} } \over { (\pp{u}{z})^2 + (\pp{v}{z})^2 } }
 =  {  {c_p \pp{\theta_v}{z} \pp{P^ \kappa}{z} } \over { (\pp{u}{z})^2 + (\pp{v}{z})^2 } }
\]
\\
where we used the hydrostatic equation: 
\[
{\pp{\Phi}{P^ \kappa}} = c_p \theta_v
\]
Negative values indicate unstable buoyancy {\bf{AND}} shear, small positive values ($<0.4$)
indicate dominantly unstable shear, and large positive values indicate dominantly stable
stratification.
\\

\noindent
{\bf 9)  \underline {CT}  Surface Exchange Coefficient for Temperature and Moisture ($dimensionless$) }

\noindent
The surface exchange coefficient is obtained from the similarity functions for the stability
 dependant flux profile relationships:
\[
{\bf CT} = -{( {\overline{w^{\prime}\theta^{\prime}}}) \over {u_* \Delta \theta }} = 
-{( {\overline{w^{\prime}q^{\prime}}}) \over {u_* \Delta q }} = 
{ k \over { (\psi_{h} + \psi_{g}) } } 
\]
where $\psi_h$ is the surface layer non-dimensional temperature change and $\psi_g$ is the
viscous sublayer non-dimensional temperature or moisture change:
\[
\psi_{h} = {\int_{\zeta_{0}}^{\zeta} {\phi_{h} \over \zeta} d \zeta} \hspace{1cm} and 
\hspace{1cm} \psi_{g} = { 0.55 (Pr^{2/3} - 0.2) \over \nu^{1/2} } 
(h_{0}u_{*} - h_{0_{ref}}u_{*_{ref}})^{1/2}
\]
and:
$h_{0} = 30z_{0}$ with a maximum value over land of 0.01

\noindent
$\phi_h$ is the similarity function of $\zeta$, which expresses the stability dependance of
the temperature and moisture gradients, specified differently for stable and unstable 
layers according to Helfand and Schubert, 1993. k is the Von Karman constant, $\zeta$ is the 
non-dimensional stability parameter, Pr is the Prandtl number for air, $\nu$ is the molecular 
viscosity, $z_{0}$ is the surface roughness length, $u_*$ is the surface stress velocity 
(see diagnostic number 67), and the subscript ref refers to a reference value.
\\

\noindent
{\bf 10)  \underline {CU}  Surface Exchange Coefficient for Momentum ($dimensionless$) }

\noindent
The surface exchange coefficient is obtained from the similarity functions for the stability
 dependant flux profile relationships:
\[
{\bf CU} = {u_* \over W_s} = { k \over \psi_{m} } 
\]
where $\psi_m$ is the surface layer non-dimensional wind shear: 
\[
\psi_{m} = {\int_{\zeta_{0}}^{\zeta} {\phi_{m} \over \zeta} d \zeta}
\]
\noindent
$\phi_m$ is the similarity function of $\zeta$, which expresses the stability dependance of
the temperature and moisture gradients, specified differently for stable and unstable layers
according to Helfand and Schubert, 1993. k is the Von Karman constant, $\zeta$ is the 
non-dimensional stability parameter, $u_*$ is the surface stress velocity 
(see diagnostic number 67), and $W_s$ is the magnitude of the surface layer wind.
\\

\noindent
{\bf 11)  \underline {ET}  Diffusivity Coefficient for Temperature and Moisture ($m^2/sec$) }

\noindent
In the level 2.5 version of the Mellor-Yamada (1974) hierarchy, the turbulent heat or
moisture flux for the atmosphere above the surface layer can be expressed as a turbulent 
diffusion coefficient $K_h$ times the negative of the gradient of potential temperature 
or moisture. In the Helfand and Labraga (1988) adaptation of this closure, $K_h$ 
takes the form:
\[
{\bf ET} = K_h = -{( {\overline{w^{\prime}\theta_v^{\prime}}}) \over {\pp{\theta_v}{z}} }
 = \left\{ \begin{array}{l@{\quad\mbox{for}\quad}l} q \, \ell \, S_H(G_M,G_H) & \mbox{decaying turbulence}
\\ { q^2 \over {q_e} } \, \ell \, S_{H}(G_{M_e},G_{H_e}) & \mbox{growing turbulence} \end{array} \right.
\]
where $q$ is the turbulent velocity, or $\sqrt{2*turbulent \hspace{.2cm} kinetic \hspace{.2cm} 
energy}$, $q_e$ is the turbulence velocity derived from the more simple level 2.0 model, 
which describes equilibrium turbulence, $\ell$ is the master length scale related to the layer 
depth, 
$S_H$ is a function of $G_H$ and $G_M$, the dimensionless buoyancy and
wind shear parameters, respectively, or a function of $G_{H_e}$ and $G_{M_e}$, the equilibrium 
dimensionless buoyancy and wind shear
parameters.   Both $G_H$ and $G_M$, and their equilibrium values $G_{H_e}$ and $G_{M_e}$, 
are functions of the Richardson number.

\noindent
For the detailed equations and derivations of the modified level 2.5 closure scheme,
see Helfand and Labraga, 1988.

\noindent
In the surface layer, ${\bf {ET}}$ is the exchange coefficient for heat and moisture,
in units of $m/sec$, given by:
\[
{\bf ET_{Nrphys}} =  C_t * u_* = C_H W_s
\]
\noindent
where $C_t$ is the dimensionless exchange coefficient for heat and moisture from the 
surface layer similarity functions (see diagnostic number 9), $u_*$ is the surface 
friction velocity (see diagnostic number 67), $C_H$ is the heat transfer coefficient,
and $W_s$ is the magnitude of the surface layer wind.
\\
 
\noindent
{\bf 12)  \underline {EU}  Diffusivity Coefficient for Momentum ($m^2/sec$) }
 
\noindent  
In the level 2.5 version of the Mellor-Yamada (1974) hierarchy, the turbulent heat
momentum flux for the atmosphere above the surface layer can be expressed as a turbulent
diffusion coefficient $K_m$ times the negative of the gradient of the u-wind.
In the Helfand and Labraga (1988) adaptation of this closure, $K_m$
takes the form:
\[
{\bf EU} = K_m = -{( {\overline{u^{\prime}w^{\prime}}}) \over {\pp{U}{z}} }
 = \left\{ \begin{array}{l@{\quad\mbox{for}\quad}l} q \, \ell \, S_M(G_M,G_H) & \mbox{decaying turbulence}
\\ { q^2 \over {q_e} } \, \ell \, S_{M}(G_{M_e},G_{H_e}) & \mbox{growing turbulence} \end{array} \right.
\]
\noindent
where $q$ is the turbulent velocity, or $\sqrt{2*turbulent \hspace{.2cm} kinetic \hspace{.2cm}
energy}$, $q_e$ is the turbulence velocity derived from the more simple level 2.0 model,
which describes equilibrium turbulence, $\ell$ is the master length scale related to the layer
depth, 
$S_M$ is a function of $G_H$ and $G_M$, the dimensionless buoyancy and
wind shear parameters, respectively, or a function of $G_{H_e}$ and $G_{M_e}$, the equilibrium 
dimensionless buoyancy and wind shear
parameters.   Both $G_H$ and $G_M$, and their equilibrium values $G_{H_e}$ and $G_{M_e}$, 
are functions of the Richardson number.

\noindent
For the detailed equations and derivations of the modified level 2.5 closure scheme,
see Helfand and Labraga, 1988.
 
\noindent
In the surface layer, ${\bf {EU}}$ is the exchange coefficient for momentum,
in units of $m/sec$, given by:
\[
{\bf EU_{Nrphys}} = C_u * u_* = C_D W_s
\]
\noindent
where $C_u$ is the dimensionless exchange coefficient for momentum from the surface layer 
similarity functions (see diagnostic number 10), $u_*$ is the surface friction velocity 
(see diagnostic number 67), $C_D$ is the surface drag coefficient, and $W_s$ is the 
magnitude of the surface layer wind.
\\
 
\noindent
{\bf 13)  \underline {TURBU}  Zonal U-Momentum changes due to Turbulence ($m/sec/day$) }
 
\noindent
The tendency of U-Momentum due to turbulence is written:
\[
{\bf TURBU} = {\pp{u}{t}}_{turb} = {\pp{}{z} }{(- \overline{u^{\prime}w^{\prime}})}
 = {\pp{}{z} }{(K_m \pp{u}{z})}
\]

\noindent
The Helfand and Labraga level 2.5 scheme models the turbulent
flux of u-momentum in terms of $K_m$, and the equation has the form of a diffusion
equation.
 
\noindent
{\bf 14)  \underline {TURBV}  Meridional V-Momentum changes due to Turbulence ($m/sec/day$) }
 
\noindent
The tendency of V-Momentum due to turbulence is written:
\[
{\bf TURBV} = {\pp{v}{t}}_{turb} = {\pp{}{z} }{(- \overline{v^{\prime}w^{\prime}})}
 = {\pp{}{z} }{(K_m \pp{v}{z})}
\]

\noindent
The Helfand and Labraga level 2.5 scheme models the turbulent
flux of v-momentum in terms of $K_m$, and the equation has the form of a diffusion
equation.
\\
 
\noindent
{\bf 15)  \underline {TURBT}  Temperature changes due to Turbulence ($deg/day$) }
 
\noindent
The tendency of temperature due to turbulence is written:
\[
{\bf TURBT} = {\pp{T}{t}} = P^{\kappa}{\pp{\theta}{t}}_{turb} = 
P^{\kappa}{\pp{}{z} }{(- \overline{w^{\prime}\theta^{\prime}})}
 = P^{\kappa}{\pp{}{z} }{(K_h \pp{\theta_v}{z})}
\]

\noindent
The Helfand and Labraga level 2.5 scheme models the turbulent
flux of temperature in terms of $K_h$, and the equation has the form of a diffusion
equation.
\\
 
\noindent
{\bf 16)  \underline {TURBQ}  Specific Humidity changes due to Turbulence ($g/kg/day$) }
 
\noindent
The tendency of specific humidity due to turbulence is written:
\[
{\bf TURBQ} = {\pp{q}{t}}_{turb} = {\pp{}{z} }{(- \overline{w^{\prime}q^{\prime}})}
 = {\pp{}{z} }{(K_h \pp{q}{z})}
\]

\noindent
The Helfand and Labraga level 2.5 scheme models the turbulent
flux of temperature in terms of $K_h$, and the equation has the form of a diffusion
equation.
\\
 
\noindent
{\bf 17)  \underline {MOISTT} Temperature Changes Due to Moist Processes ($deg/day$) } 

\noindent
\[
{\bf MOISTT} = \left. {\pp{T}{t}}\right|_{c} + \left. {\pp{T}{t}} \right|_{ls}
\]
where:
\[
\left.{\pp{T}{t}}\right|_{c} = R \sum_i \left( \alpha { m_B \over c_p} \Gamma_s \right)_i 
\hspace{.4cm} and 
\hspace{.4cm} \left.{\pp{T}{t}}\right|_{ls} = {L \over c_p } (q^*-q)
\]
and
\[
\Gamma_s = g \eta \pp{s}{p}
\]

\noindent
The subscript $c$ refers to convective processes, while the subscript $ls$ refers to large scale
precipitation processes, or supersaturation rain. 
The summation refers to contributions from each cloud type called by RAS.  
The dry static energy is given 
as $s$, the convective cloud base mass flux is given as $m_B$, and the cloud entrainment is
given as $\eta$, which are explicitly defined in Section \ref{sec:fizhi:mc}, 
the description of the convective parameterization.  The fractional adjustment, or relaxation
parameter, for each cloud type is given as $\alpha$, while
$R$ is the rain re-evaporation adjustment.
\\

\noindent
{\bf 18)  \underline {MOISTQ} Specific Humidity Changes Due to Moist Processes ($g/kg/day$) } 

\noindent
\[
{\bf MOISTQ} = \left. {\pp{q}{t}}\right|_{c} + \left. {\pp{q}{t}} \right|_{ls}
\]
where:
\[
\left.{\pp{q}{t}}\right|_{c} = R \sum_i \left( \alpha { m_B \over {L}}(\Gamma_h-\Gamma_s) \right)_i 
\hspace{.4cm} and 
\hspace{.4cm} \left.{\pp{q}{t}}\right|_{ls} = (q^*-q)
\]
and
\[
\Gamma_s = g \eta \pp{s}{p}\hspace{.4cm} and \hspace{.4cm}\Gamma_h = g \eta \pp{h}{p}
\]
\noindent
The subscript $c$ refers to convective processes, while the subscript $ls$ refers to large scale
precipitation processes, or supersaturation rain. 
The summation refers to contributions from each cloud type called by RAS.  
The dry static energy is given as $s$, 
the moist static energy is given as $h$, 
the convective cloud base mass flux is given as $m_B$, and the cloud entrainment is
given as $\eta$, which are explicitly defined in Section \ref{sec:fizhi:mc}, 
the description of the convective parameterization.  The fractional adjustment, or relaxation
parameter, for each cloud type is given as $\alpha$, while
$R$ is the rain re-evaporation adjustment.
\\

\noindent
{\bf 19)  \underline {RADLW} Heating Rate due to Longwave Radiation ($deg/day$) }

\noindent
The net longwave heating rate is calculated as the vertical divergence of the
net terrestrial radiative fluxes.
Both the clear-sky and cloudy-sky longwave fluxes are computed within the
longwave routine.
The subroutine calculates the clear-sky flux, $F^{clearsky}_{LW}$, first.
For a given cloud fraction,
the clear line-of-sight probability $C(p,p^{\prime})$ is computed from the current level pressure $p$ 
to the model top pressure, $p^{\prime} = p_{top}$, and the model surface pressure, $p^{\prime} = p_{surf}$,
for the upward and downward radiative fluxes.
(see Section \ref{sec:fizhi:radcloud}).
The cloudy-sky flux is then obtained as:
   
\noindent
\[
F_{LW} = C(p,p') \cdot F^{clearsky}_{LW},
\]

\noindent
Finally, the net longwave heating rate is calculated as the vertical divergence of the
net terrestrial radiative fluxes:
\[
\pp{\rho c_p T}{t} = - {\partial \over \partial z} F_{LW}^{NET} ,
\]
or
\[
{\bf RADLW} = \frac{g}{c_p \pi} {\partial \over \partial \sigma} F_{LW}^{NET} .
\]

\noindent
where $g$ is the accelation due to gravity,
$c_p$ is the heat capacity of air at constant pressure,
and
\[
F_{LW}^{NET} = F_{LW}^\uparrow - F_{LW}^\downarrow
\]
\\


\noindent
{\bf 20)  \underline {RADSW} Heating Rate due to Shortwave Radiation ($deg/day$) }

\noindent
The net Shortwave heating rate is calculated as the vertical divergence of the
net solar radiative fluxes.
The clear-sky and cloudy-sky shortwave fluxes are calculated separately.
For the clear-sky case, the shortwave fluxes and heating rates are computed with
both CLMO (maximum overlap cloud fraction) and
CLRO (random overlap cloud fraction) set to zero (see Section \ref{sec:fizhi:radcloud}).
The shortwave routine is then called a second time, for the cloudy-sky case, with the
true time-averaged cloud fractions CLMO
and CLRO being used.  In all cases, a normalized incident shortwave flux is used as
input at the top of the atmosphere.

\noindent
The heating rate due to Shortwave Radiation under cloudy skies is defined as:
\[
\pp{\rho c_p T}{t} = - {\partial \over \partial z} F(cloudy)_{SW}^{NET} \cdot {\rm RADSWT},
\]
or
\[
{\bf RADSW} = \frac{g}{c_p \pi} {\partial \over \partial \sigma} F(cloudy)_{SW}^{NET}\cdot {\rm RADSWT} .
\]

\noindent
where $g$ is the accelation due to gravity,
$c_p$ is the heat capacity of air at constant pressure, RADSWT is the true incident
shortwave radiation at the top of the atmosphere (See Diagnostic \#48), and
\[
F(cloudy)_{SW}^{Net} = F(cloudy)_{SW}^\uparrow - F(cloudy)_{SW}^\downarrow
\]
\\

\noindent
{\bf 21)  \underline {PREACC} Total (Large-scale + Convective) Accumulated Precipition ($mm/day$) } 

\noindent
For a change in specific humidity due to moist processes, $\Delta q_{moist}$, 
the vertical integral or total precipitable amount is given by:   
\[
{\bf PREACC} = \int_{surf}^{top} \rho \Delta q_{moist} dz = - \int_{surf}^{top} \Delta  q_{moist}
{dp \over g} = {1 \over g} \int_0^1 \Delta q_{moist} dp
\]
\\

\noindent
A precipitation rate is defined as the vertically integrated moisture adjustment per Moist Processes
time step, scaled to $mm/day$.
\\

\noindent
{\bf 22)  \underline {PRECON} Convective Precipition ($mm/day$) } 

\noindent
For a change in specific humidity due to sub-grid scale cumulus convective processes, $\Delta q_{cum}$, 
the vertical integral or total precipitable amount is given by:   
\[
{\bf PRECON} = \int_{surf}^{top} \rho \Delta q_{cum} dz = - \int_{surf}^{top} \Delta  q_{cum}
{dp \over g} = {1 \over g} \int_0^1 \Delta q_{cum} dp
\]
\\

\noindent
A precipitation rate is defined as the vertically integrated moisture adjustment per Moist Processes
time step, scaled to $mm/day$.
\\

\noindent
{\bf 23)  \underline {TUFLUX}  Turbulent Flux of U-Momentum ($Newton/m^2$) }

\noindent
The turbulent flux of u-momentum is calculated for $diagnostic \hspace{.2cm} purposes
 \hspace{.2cm} only$ from the eddy coefficient for momentum:

\[
{\bf TUFLUX} =  {\rho } {(\overline{u^{\prime}w^{\prime}})} =  
{\rho } {(- K_m \pp{U}{z})}
\]
 
\noindent
where $\rho$ is the air density, and $K_m$ is the eddy coefficient.
\\

\noindent
{\bf 24)  \underline {TVFLUX}  Turbulent Flux of V-Momentum ($Newton/m^2$) }

\noindent
The turbulent flux of v-momentum is calculated for $diagnostic \hspace{.2cm} purposes 
\hspace{.2cm} only$ from the eddy coefficient for momentum:

\[
{\bf TVFLUX} =  {\rho } {(\overline{v^{\prime}w^{\prime}})} = 
 {\rho } {(- K_m \pp{V}{z})}
\]
 
\noindent
where $\rho$ is the air density, and $K_m$ is the eddy coefficient.
\\


\noindent
{\bf 25)  \underline {TTFLUX}  Turbulent Flux of Sensible Heat ($Watts/m^2$) }

\noindent
The turbulent flux of sensible heat is calculated for $diagnostic \hspace{.2cm} purposes 
\hspace{.2cm} only$ from the eddy coefficient for heat and moisture:

\noindent
\[
{\bf TTFLUX} = c_p {\rho }  
P^{\kappa}{(\overline{w^{\prime}\theta^{\prime}})}
 = c_p  {\rho } P^{\kappa}{(- K_h \pp{\theta_v}{z})}
\]
 
\noindent
where $\rho$ is the air density, and $K_h$ is the eddy coefficient.
\\


\noindent
{\bf 26)  \underline {TQFLUX}  Turbulent Flux of Latent Heat ($Watts/m^2$) }

\noindent
The turbulent flux of latent heat is calculated for $diagnostic \hspace{.2cm} purposes 
\hspace{.2cm} only$ from the eddy coefficient for heat and moisture:

\noindent
\[
{\bf TQFLUX} = {L {\rho } (\overline{w^{\prime}q^{\prime}})} = 
{L {\rho }(- K_h \pp{q}{z})}
\]
 
\noindent
where $\rho$ is the air density, and $K_h$ is the eddy coefficient.
\\

 
\noindent
{\bf 27)  \underline {CN}  Neutral Drag Coefficient ($dimensionless$) }

\noindent
The drag coefficient for momentum obtained by assuming a neutrally stable surface layer:
\[
{\bf CN} = { k \over { \ln({h \over {z_0}})} }
\]

\noindent
where $k$ is the Von Karman constant, $h$ is the height of the surface layer, and
$z_0$ is the surface roughness. 

\noindent
NOTE: CN is not available through model version 5.3, but is available in subsequent
versions.
\\

\noindent
{\bf 28)  \underline {WINDS}  Surface Wind Speed ($meter/sec$) }

\noindent
The surface wind speed is calculated for the last internal turbulence time step:
\[
{\bf WINDS} = \sqrt{u_{Nrphys}^2 + v_{Nrphys}^2}
\]

\noindent
where the subscript $Nrphys$ refers to the lowest model level.
\\
 
\noindent
{\bf 29)  \underline {DTSRF}  Air/Surface Virtual Temperature Difference ($deg \hspace{.1cm} K$) }

\noindent
The air/surface virtual temperature difference measures the stability of the surface layer:
\[
{\bf DTSRF} = (\theta_{v{Nrphys+1}} - \theta{v_{Nrphys}}) P^{\kappa}_{surf}
\]
\noindent
where
\[
\theta_{v{Nrphys+1}} = { T_g \over {P^{\kappa}_{surf}} } (1 + .609 q_{Nrphys+1}) \hspace{1cm}
and \hspace{1cm} q_{Nrphys+1} = q_{Nrphys} + \beta(q^*(T_g,P_s) - q_{Nrphys})
\]

\noindent
$\beta$ is the surface potential evapotranspiration coefficient ($\beta=1$ over oceans),
$q^*(T_g,P_s)$ is the saturation specific humidity at the ground temperature 
and surface pressure, level $Nrphys$ refers to the lowest model level and level $Nrphys+1$ 
refers to the surface.
\\

 
\noindent
{\bf 30)  \underline {TG}  Ground Temperature ($deg \hspace{.1cm} K$) }

\noindent
The ground temperature equation is solved as part of the turbulence package
using a backward implicit time differencing scheme:
\[
{\bf TG} \hspace{.1cm} is \hspace{.1cm} obtained \hspace{.1cm} from: \hspace{.1cm}
C_g\pp{T_g}{t} = R_{sw} - R_{lw} + Q_{ice} - H - LE
\]

\noindent
where $R_{sw}$ is the net surface downward shortwave radiative flux, $R_{lw}$ is the
net surface upward longwave radiative flux, $Q_{ice}$ is the heat conduction through
sea ice, $H$ is the upward sensible heat flux, $LE$ is the upward latent heat
flux, and $C_g$ is the total heat capacity of the ground. 
$C_g$ is obtained by solving a heat diffusion equation 
for the penetration of the diurnal cycle into the ground (Blackadar, 1977), and is given by:
\[
C_g = \sqrt{ {\lambda C_s \over {2 \omega} } } = \sqrt{(0.386 + 0.536W + 0.15W^2)2x10^{-3}
{ 86400. \over {2 \pi} } } \, \, .
\]
\noindent
Here, the thermal conductivity, $\lambda$, is equal to $2x10^{-3}$ ${ly\over{ sec}} 
{cm \over {^oK}}$, 
the angular velocity of the earth, $\omega$, is written as $86400$ $sec/day$ divided 
by $2 \pi$ $radians/
day$, and the expression for $C_s$, the heat capacity per unit volume at the surface, 
is a function of the ground wetness, $W$. 
\\

\noindent
{\bf 31)  \underline {TS}  Surface Temperature ($deg \hspace{.1cm} K$) }

\noindent
The surface temperature estimate is made by assuming that the model's lowest
layer is well-mixed, and therefore that $\theta$ is constant in that layer.
The surface temperature is therefore:
\[
{\bf TS} = \theta_{Nrphys} P^{\kappa}_{surf}
\]
\\
 
\noindent
{\bf 32)  \underline {DTG}  Surface Temperature Adjustment ($deg \hspace{.1cm} K$) }

\noindent
The change in surface temperature from one turbulence time step to the next, solved
using the Ground Temperature Equation (see diagnostic number 30) is calculated:
\[
{\bf DTG} = {T_g}^{n} - {T_g}^{n-1}
\]

\noindent
where superscript $n$ refers to the new, updated time level, and the superscript $n-1$
refers to the value at the previous turbulence time level.
\\
 
\noindent
{\bf 33)  \underline {QG}  Ground Specific Humidity ($g/kg$) }

\noindent
The ground specific humidity is obtained by interpolating between the specific
humidity at the lowest model level and the specific humidity of a saturated ground.
The interpolation is performed using the potential evapotranspiration function:
\[
{\bf QG} = q_{Nrphys+1} = q_{Nrphys} + \beta(q^*(T_g,P_s) - q_{Nrphys})
\]

\noindent
where $\beta$ is the surface potential evapotranspiration coefficient ($\beta=1$ over oceans), 
and $q^*(T_g,P_s)$ is the saturation specific humidity at the ground temperature and surface
pressure.
\\
 
\noindent
{\bf 34)  \underline {QS}  Saturation Surface Specific Humidity ($g/kg$) }

\noindent
The surface saturation specific humidity is the saturation specific humidity at
the ground temprature and surface pressure:
\[
{\bf QS} = q^*(T_g,P_s)
\]
\\
 
\noindent
{\bf 35)  \underline {TGRLW} Instantaneous ground temperature used as input to the Longwave
 radiation subroutine (deg)}
\[
{\bf TGRLW}  = T_g(\lambda , \phi ,n)
\]
\noindent
where $T_g$ is the model ground temperature at the current time step $n$.
\\
 
 
\noindent
{\bf 36)  \underline {ST4} Upward Longwave flux at the surface ($Watts/m^2$) }
\[
{\bf ST4} = \sigma T^4
\]
\noindent
where $\sigma$ is the Stefan-Boltzmann constant and T is the temperature.
\\
 
\noindent
{\bf 37)  \underline {OLR} Net upward Longwave flux at $p=p_{top}$ ($Watts/m^2$) }
\[
{\bf OLR}  =  F_{LW,top}^{NET}
\]
\noindent
where top indicates the top of the first model layer.
In the GCM, $p_{top}$ = 0.0 mb.
\\


\noindent
{\bf 38)  \underline {OLRCLR} Net upward clearsky Longwave flux at $p=p_{top}$ ($Watts/m^2$) }
\[
{\bf OLRCLR}  =  F(clearsky)_{LW,top}^{NET}
\]
\noindent
where top indicates the top of the first model layer.
In the GCM, $p_{top}$ = 0.0 mb.
\\

\noindent
{\bf 39)  \underline {LWGCLR} Net upward clearsky Longwave flux at the surface ($Watts/m^2$) }

\noindent
\begin{eqnarray*}
{\bf LWGCLR} & =  & F(clearsky)_{LW,Nrphys+1}^{Net} \\
             & =  & F(clearsky)_{LW,Nrphys+1}^\uparrow - F(clearsky)_{LW,Nrphys+1}^\downarrow
\end{eqnarray*}
where Nrphys+1 indicates the lowest model edge-level, or $p = p_{surf}$.
$F(clearsky)_{LW}^\uparrow$ is
the upward clearsky Longwave flux and the $F(clearsky)_{LW}^\downarrow$ is the downward clearsky Longwave flux.
\\

\noindent
{\bf 40)  \underline {LWCLR} Heating Rate due to Clearsky Longwave Radiation ($deg/day$) }

\noindent
The net longwave heating rate is calculated as the vertical divergence of the
net terrestrial radiative fluxes.
Both the clear-sky and cloudy-sky longwave fluxes are computed within the
longwave routine.
The subroutine calculates the clear-sky flux, $F^{clearsky}_{LW}$, first.
For a given cloud fraction,
the clear line-of-sight probability $C(p,p^{\prime})$ is computed from the current level pressure $p$ 
to the model top pressure, $p^{\prime} = p_{top}$, and the model surface pressure, $p^{\prime} = p_{surf}$,
for the upward and downward radiative fluxes.
(see Section \ref{sec:fizhi:radcloud}).
The cloudy-sky flux is then obtained as:
   
\noindent
\[
F_{LW} = C(p,p') \cdot F^{clearsky}_{LW},
\]

\noindent
Thus, {\bf LWCLR} is defined as the net longwave heating rate due to the 
vertical divergence of the
clear-sky longwave radiative flux:
\[
\pp{\rho c_p T}{t}_{clearsky} = - {\partial \over \partial z} F(clearsky)_{LW}^{NET} ,
\]
or
\[
{\bf LWCLR} = \frac{g}{c_p \pi} {\partial \over \partial \sigma} F(clearsky)_{LW}^{NET} .
\]

\noindent
where $g$ is the accelation due to gravity,
$c_p$ is the heat capacity of air at constant pressure,
and
\[
F(clearsky)_{LW}^{Net} = F(clearsky)_{LW}^\uparrow - F(clearsky)_{LW}^\downarrow
\]
\\

 
\noindent
{\bf 41)  \underline {TLW} Instantaneous temperature used as input to the Longwave
 radiation subroutine (deg)}
\[
{\bf TLW}  = T(\lambda , \phi ,level, n)
\]
\noindent
where $T$ is the model temperature at the current time step $n$.
\\
 
 
\noindent
{\bf 42)  \underline {SHLW} Instantaneous specific humidity used as input to
 the Longwave radiation subroutine (kg/kg)}
\[
{\bf SHLW}  = q(\lambda , \phi , level , n)
\]
\noindent
where $q$ is the model specific humidity at the current time step $n$.
\\
 
 
\noindent
{\bf 43)  \underline {OZLW} Instantaneous ozone used as input to
 the Longwave radiation subroutine (kg/kg)}
\[
{\bf OZLW}  = {\rm OZ}(\lambda , \phi , level , n)
\]
\noindent
where $\rm OZ$ is the interpolated ozone data set from the climatological monthly
mean zonally averaged ozone data set.
\\
 

\noindent
{\bf 44) \underline {CLMOLW} Maximum Overlap cloud fraction used in LW Radiation ($0-1$) }

\noindent
{\bf CLMOLW} is the time-averaged maximum overlap cloud fraction that has been filled by the Relaxed
Arakawa/Schubert Convection scheme and will be used in the Longwave Radiation algorithm.  These are
convective clouds whose radiative characteristics are assumed to be correlated in the vertical.
For a complete description of cloud/radiative interactions, see Section \ref{sec:fizhi:radcloud}.
\[
{\bf CLMOLW} = CLMO_{RAS,LW}(\lambda, \phi,  level )
\]
\\
 

{\bf 45) \underline {CLDTOT} Total cloud fraction used in LW and SW Radiation ($0-1$) }

{\bf CLDTOT} is the time-averaged total cloud fraction that has been filled by the Relaxed
Arakawa/Schubert and Large-scale Convection schemes and will be used in the Longwave and Shortwave
Radiation packages.
For a complete description of cloud/radiative interactions, see Section \ref{sec:fizhi:radcloud}.
\[
{\bf CLDTOT} = F_{RAS} + F_{LS}
\]
\\
where $F_{RAS}$ is the time-averaged cloud fraction due to sub-grid scale convection, and $F_{LS}$ is the
time-averaged cloud fraction due to precipitating and non-precipitating large-scale moist processes.
\\


\noindent
{\bf 46) \underline {CLMOSW} Maximum Overlap cloud fraction used in SW Radiation ($0-1$) }

\noindent
{\bf CLMOSW} is the time-averaged maximum overlap cloud fraction that has been filled by the Relaxed
Arakawa/Schubert Convection scheme and will be used in the Shortwave Radiation algorithm.  These are
convective clouds whose radiative characteristics are assumed to be correlated in the vertical.
For a complete description of cloud/radiative interactions, see Section \ref{sec:fizhi:radcloud}.
\[
{\bf CLMOSW} = CLMO_{RAS,SW}(\lambda, \phi,  level )
\]
\\

\noindent
{\bf 47) \underline {CLROSW} Random Overlap cloud fraction used in SW Radiation ($0-1$) }

\noindent
{\bf CLROSW} is the time-averaged random overlap cloud fraction that has been filled by the Relaxed
Arakawa/Schubert and Large-scale Convection schemes and will be used in the Shortwave 
Radiation algorithm.  These are
convective and large-scale clouds whose radiative characteristics are not 
assumed to be correlated in the vertical.
For a complete description of cloud/radiative interactions, see Section \ref{sec:fizhi:radcloud}.
\[
{\bf CLROSW} = CLRO_{RAS,Large Scale,SW}(\lambda, \phi,  level )
\]
\\

\noindent
{\bf 48)  \underline {RADSWT} Incident Shortwave radiation at the top of the atmosphere ($Watts/m^2$) }
\[
{\bf RADSWT} = {\frac{S_0}{R_a^2}} \cdot cos \phi_z
\]
\noindent
where $S_0$, is the extra-terrestial solar contant,
$R_a$ is the earth-sun distance in Astronomical Units,
and $cos \phi_z$ is the cosine of the zenith angle.
It should be noted that {\bf RADSWT}, as well as
{\bf OSR} and {\bf OSRCLR}, 
are calculated at the top of the atmosphere (p=0 mb).  However, the
{\bf OLR} and {\bf OLRCLR} diagnostics are currently
calculated at $p= p_{top}$ (0.0 mb for the GCM).
\\
   
\noindent
{\bf 49)  \underline {EVAP}  Surface Evaporation ($mm/day$) }

\noindent
The surface evaporation is a function of the gradient of moisture, the potential 
evapotranspiration fraction and the eddy exchange coefficient:
\[
{\bf EVAP} =  \rho \beta K_{h} (q_{surface} - q_{Nrphys})
\]
where $\rho$ = the atmospheric density at the surface, $\beta$ is the fraction of
the potential evapotranspiration actually evaporated ($\beta=1$ over oceans), $K_{h}$ is the 
turbulent eddy exchange coefficient for heat and moisture at the surface in $m/sec$ and 
$q{surface}$ and $q_{Nrphys}$ are the specific humidity at the surface (see diagnostic
number 34) and at the bottom model level, respectively.
\\

\noindent
{\bf 50)  \underline {DUDT} Total Zonal U-Wind Tendency  ($m/sec/day$) }

\noindent
{\bf DUDT} is the total time-tendency of the Zonal U-Wind due to Hydrodynamic, Diabatic,
and Analysis forcing.
\[
{\bf DUDT} = \pp{u}{t}_{Dynamics} + \pp{u}{t}_{Moist} + \pp{u}{t}_{Turbulence} + \pp{u}{t}_{Analysis} 
\]
\\

\noindent
{\bf 51)  \underline {DVDT} Total Zonal V-Wind Tendency  ($m/sec/day$) }

\noindent
{\bf DVDT} is the total time-tendency of the Meridional V-Wind due to Hydrodynamic, Diabatic,
and Analysis forcing.
\[
{\bf DVDT} = \pp{v}{t}_{Dynamics} + \pp{v}{t}_{Moist} + \pp{v}{t}_{Turbulence} + \pp{v}{t}_{Analysis} 
\]
\\

\noindent
{\bf 52)  \underline {DTDT} Total Temperature Tendency  ($deg/day$) }

\noindent
{\bf DTDT} is the total time-tendency of Temperature due to Hydrodynamic, Diabatic,
and Analysis forcing.
\begin{eqnarray*}
{\bf DTDT} & = & \pp{T}{t}_{Dynamics} + \pp{T}{t}_{Moist Processes} + \pp{T}{t}_{Shortwave Radiation} \\
           & + & \pp{T}{t}_{Longwave Radiation} + \pp{T}{t}_{Turbulence} + \pp{T}{t}_{Analysis} 
\end{eqnarray*}
\\

\noindent
{\bf 53)  \underline {DQDT} Total Specific Humidity Tendency  ($g/kg/day$) }

\noindent
{\bf DQDT} is the total time-tendency of Specific Humidity due to Hydrodynamic, Diabatic,
and Analysis forcing.
\[
{\bf DQDT} = \pp{q}{t}_{Dynamics} + \pp{q}{t}_{Moist Processes} 
+ \pp{q}{t}_{Turbulence} + \pp{q}{t}_{Analysis} 
\]
\\
   
\noindent
{\bf 54)  \underline {USTAR}  Surface-Stress Velocity ($m/sec$) }

\noindent
The surface stress velocity, or the friction velocity, is the wind speed at 
the surface layer top impeded by the surface drag:
\[
{\bf USTAR} = C_uW_s \hspace{1cm}where: \hspace{.2cm} 
C_u = {k \over {\psi_m} }
\]

\noindent
$C_u$ is the non-dimensional surface drag coefficient (see diagnostic
number 10), and $W_s$ is the surface wind speed (see diagnostic number 28).
 
\noindent
{\bf 55)  \underline {Z0}  Surface Roughness Length ($m$) }

\noindent
Over the land surface, the surface roughness length is interpolated to the local
time from the monthly mean data of Dorman and Sellers (1989). Over the ocean,
the roughness length is a function of the surface-stress velocity, $u_*$.
\[
{\bf Z0} = c_1u^3_* + c_2u^2_* + c_3u_* + c_4 + {c_5 \over {u_*}}
\]

\noindent
where the constants are chosen to interpolate between the reciprocal relation of
Kondo(1975) for weak winds, and the piecewise linear relation of Large and Pond(1981) 
for moderate to large winds.
\\
 
\noindent
{\bf 56)  \underline {FRQTRB}  Frequency of Turbulence ($0-1$) }

\noindent
The fraction of time when turbulence is present is defined as the fraction of
time when the turbulent kinetic energy exceeds some minimum value, defined here
to be $0.005 \hspace{.1cm}m^2/sec^2$. When this criterion is met, a counter is
incremented. The fraction over the averaging interval is reported.
\\
 
\noindent
{\bf 57)  \underline {PBL}  Planetary Boundary Layer Depth ($mb$) }

\noindent
The depth of the PBL is defined by the turbulence parameterization to be the
depth at which the turbulent kinetic energy reduces to ten percent of its surface
value.

\[
{\bf PBL} = P_{PBL} - P_{surface}
\]

\noindent
where $P_{PBL}$ is the pressure in $mb$ at which the turbulent kinetic energy
reaches one tenth of its surface value, and $P_s$ is the surface pressure.
\\
 
\noindent
{\bf 58)  \underline {SWCLR} Clear sky Heating Rate due to Shortwave Radiation ($deg/day$) }

\noindent
The net Shortwave heating rate is calculated as the vertical divergence of the
net solar radiative fluxes.
The clear-sky and cloudy-sky shortwave fluxes are calculated separately.
For the clear-sky case, the shortwave fluxes and heating rates are computed with
both CLMO (maximum overlap cloud fraction) and
CLRO (random overlap cloud fraction) set to zero (see Section \ref{sec:fizhi:radcloud}).
The shortwave routine is then called a second time, for the cloudy-sky case, with the
true time-averaged cloud fractions CLMO
and CLRO being used.  In all cases, a normalized incident shortwave flux is used as
input at the top of the atmosphere.

\noindent
The heating rate due to Shortwave Radiation under clear skies is defined as:
\[
\pp{\rho c_p T}{t} = - {\partial \over \partial z} F(clear)_{SW}^{NET} \cdot {\rm RADSWT},
\]
or
\[
{\bf SWCLR} = \frac{g}{c_p } {\partial \over \partial p} F(clear)_{SW}^{NET}\cdot {\rm RADSWT} .
\]

\noindent
where $g$ is the accelation due to gravity,
$c_p$ is the heat capacity of air at constant pressure, RADSWT is the true incident
shortwave radiation at the top of the atmosphere (See Diagnostic \#48), and
\[
F(clear)_{SW}^{Net} = F(clear)_{SW}^\uparrow - F(clear)_{SW}^\downarrow
\]
\\

\noindent
{\bf 59)  \underline {OSR} Net upward Shortwave flux at the top of the model ($Watts/m^2$) }
\[
{\bf OSR}  =  F_{SW,top}^{NET}
\]                                                                                       
\noindent
where top indicates the top of the first model layer used in the shortwave radiation
routine.
In the GCM, $p_{SW_{top}}$ = 0 mb.
\\

\noindent
{\bf 60)  \underline {OSRCLR} Net upward clearsky Shortwave flux at the top of the model ($Watts/m^2$) }
\[
{\bf OSRCLR}  =  F(clearsky)_{SW,top}^{NET}
\]
\noindent
where top indicates the top of the first model layer used in the shortwave radiation
routine.
In the GCM, $p_{SW_{top}}$ = 0 mb.
\\


\noindent
{\bf 61)  \underline {CLDMAS} Convective Cloud Mass Flux ($kg/m^2$) } 

\noindent
The amount of cloud mass moved per RAS timestep from all convective clouds is written:
\[
{\bf CLDMAS} = \eta m_B
\]
where $\eta$ is the entrainment, normalized by the cloud base mass flux, and $m_B$ is
the cloud base mass flux. $m_B$ and $\eta$ are defined explicitly in Section \ref{sec:fizhi:mc}, the 
description of the convective parameterization.
\\


\noindent
{\bf 62)  \underline {UAVE} Time-Averaged Zonal U-Wind ($m/sec$) }

\noindent
The diagnostic {\bf UAVE} is simply the time-averaged Zonal U-Wind over
the {\bf NUAVE} output frequency.  This is contrasted to the instantaneous
Zonal U-Wind which is archived on the Prognostic Output data stream.
\[
{\bf UAVE} = u(\lambda, \phi, level , t)
\]
\\
Note, {\bf UAVE} is computed and stored on the staggered C-grid.
\\

\noindent
{\bf 63)  \underline {VAVE} Time-Averaged Meridional V-Wind ($m/sec$) }

\noindent
The diagnostic {\bf VAVE} is simply the time-averaged Meridional V-Wind over
the {\bf NVAVE} output frequency.  This is contrasted to the instantaneous
Meridional V-Wind which is archived on the Prognostic Output data stream.
\[
{\bf VAVE} = v(\lambda, \phi, level , t)
\]
\\
Note, {\bf VAVE} is computed and stored on the staggered C-grid.
\\

\noindent
{\bf 64)  \underline {TAVE} Time-Averaged Temperature ($Kelvin$) }

\noindent
The diagnostic {\bf TAVE} is simply the time-averaged Temperature over
the {\bf NTAVE} output frequency.  This is contrasted to the instantaneous
Temperature which is archived on the Prognostic Output data stream.
\[
{\bf TAVE} = T(\lambda, \phi, level , t)
\]
\\

\noindent
{\bf 65)  \underline {QAVE} Time-Averaged Specific Humidity ($g/kg$) }

\noindent
The diagnostic {\bf QAVE} is simply the time-averaged Specific Humidity over
the {\bf NQAVE} output frequency.  This is contrasted to the instantaneous
Specific Humidity which is archived on the Prognostic Output data stream.
\[
{\bf QAVE} = q(\lambda, \phi, level , t)
\]
\\

\noindent
{\bf 66)  \underline {PAVE} Time-Averaged Surface Pressure - PTOP ($mb$) }

\noindent
The diagnostic {\bf PAVE} is simply the time-averaged Surface Pressure - PTOP over
the {\bf NPAVE} output frequency.  This is contrasted to the instantaneous
Surface Pressure - PTOP which is archived on the Prognostic Output data stream.
\begin{eqnarray*}
{\bf PAVE} & =  & \pi(\lambda, \phi, level , t) \\
           & =  & p_s(\lambda, \phi, level , t) - p_T
\end{eqnarray*}
\\

 
\noindent
{\bf 67)  \underline {QQAVE} Time-Averaged Turbulent Kinetic Energy $(m/sec)^2$ }
 
\noindent
The diagnostic {\bf QQAVE} is simply the time-averaged prognostic Turbulent Kinetic Energy 
produced by the GCM Turbulence parameterization over
the {\bf NQQAVE} output frequency.  This is contrasted to the instantaneous
Turbulent Kinetic Energy which is archived on the Prognostic Output data stream.
\[
{\bf QQAVE} = qq(\lambda, \phi, level , t)
\]
\\
Note, {\bf QQAVE} is computed and stored at the ``mass-point'' locations on the staggered C-grid.
\\
 
\noindent
{\bf 68)  \underline {SWGCLR} Net downward clearsky Shortwave flux at the surface ($Watts/m^2$) }

\noindent
\begin{eqnarray*}
{\bf SWGCLR} & =  & F(clearsky)_{SW,Nrphys+1}^{Net} \\
             & =  & F(clearsky)_{SW,Nrphys+1}^\downarrow - F(clearsky)_{SW,Nrphys+1}^\uparrow
\end{eqnarray*}
\noindent
\\
where Nrphys+1 indicates the lowest model edge-level, or $p = p_{surf}$.
$F(clearsky){SW}^\downarrow$ is
the downward clearsky Shortwave flux and $F(clearsky)_{SW}^\uparrow$ is 
the upward clearsky Shortwave flux.
\\

\noindent
{\bf 69)  \underline {SDIAG1} User-Defined Surface Diagnostic-1 }

\noindent
The GCM provides Users with a built-in mechanism for archiving user-defined
diagnostics.  The generic diagnostic array QDIAG located in COMMON /DIAG/, and the associated 
diagnostic counters and pointers located in COMMON /DIAGP/,
must be accessable in order to use the user-defined diagnostics (see Section \ref{sec:diagnostics:diagover}).  
A convenient method for incorporating all necessary COMMON files is to
include the GCM {\em vstate.com} file in the routine which employs the
user-defined diagnostics.

\noindent
In addition to enabling the user-defined diagnostic (ie., CALL SETDIAG(84)), the User must fill 
the QDIAG array with the desired quantity within the User's
application program or within modified GCM subroutines, as well as increment
the diagnostic counter at the time when the diagnostic is updated.  
The QDIAG location index for {\bf SDIAG1} and its corresponding counter is 
automatically defined as {\bf ISDIAG1} and {\bf NSDIAG1}, respectively, after the 
diagnostic has been enabled.  
The syntax for its use is given by
\begin{verbatim}
      do j=1,jm
      do i=1,im
      qdiag(i,j,ISDIAG1) = qdiag(i,j,ISDIAG1) + ...
      enddo
      enddo

      NSDIAG1 = NSDIAG1 + 1
\end{verbatim}
The diagnostics defined in this manner will automatically be archived by the output routines.
\\

\noindent
{\bf 70)  \underline {SDIAG2} User-Defined Surface Diagnostic-2 }

\noindent
The GCM provides Users with a built-in mechanism for archiving user-defined
diagnostics.  For a complete description refer to Diagnostic \#84.
The syntax for using the surface SDIAG2 diagnostic is given by
\begin{verbatim}
      do j=1,jm
      do i=1,im
      qdiag(i,j,ISDIAG2) = qdiag(i,j,ISDIAG2) + ...
      enddo
      enddo

      NSDIAG2 = NSDIAG2 + 1
\end{verbatim}
The diagnostics defined in this manner will automatically be archived by the output routines.
\\

\noindent
{\bf 71)  \underline {UDIAG1} User-Defined Upper-Air Diagnostic-1 }

\noindent
The GCM provides Users with a built-in mechanism for archiving user-defined
diagnostics.  For a complete description refer to Diagnostic \#84.
The syntax for using the upper-air UDIAG1 diagnostic is given by
\begin{verbatim}
      do L=1,Nrphys
      do j=1,jm
      do i=1,im
      qdiag(i,j,IUDIAG1+L-1) = qdiag(i,j,IUDIAG1+L-1) + ...
      enddo
      enddo
      enddo

      NUDIAG1 = NUDIAG1 + 1
\end{verbatim}
The diagnostics defined in this manner will automatically be archived by the 
output programs.
\\

\noindent
{\bf 72)  \underline {UDIAG2} User-Defined Upper-Air Diagnostic-2 }

\noindent
The GCM provides Users with a built-in mechanism for archiving user-defined
diagnostics.  For a complete description refer to Diagnostic \#84.
The syntax for using the upper-air UDIAG2 diagnostic is given by
\begin{verbatim}
      do L=1,Nrphys
      do j=1,jm
      do i=1,im
      qdiag(i,j,IUDIAG2+L-1) = qdiag(i,j,IUDIAG2+L-1) + ...
      enddo
      enddo
      enddo

      NUDIAG2 = NUDIAG2 + 1
\end{verbatim}
The diagnostics defined in this manner will automatically be archived by the 
output programs.
\\


\noindent
{\bf 73)  \underline {DIABU} Total Diabatic Zonal U-Wind Tendency  ($m/sec/day$) }

\noindent
{\bf DIABU} is the total time-tendency of the Zonal U-Wind due to Diabatic processes
and the Analysis forcing.
\[
{\bf DIABU} = \pp{u}{t}_{Moist} + \pp{u}{t}_{Turbulence} + \pp{u}{t}_{Analysis} 
\]
\\

\noindent
{\bf 74)  \underline {DIABV} Total Diabatic Meridional V-Wind Tendency  ($m/sec/day$) }

\noindent
{\bf DIABV} is the total time-tendency of the Meridional V-Wind due to Diabatic processes
and the Analysis forcing.
\[
{\bf DIABV} = \pp{v}{t}_{Moist} + \pp{v}{t}_{Turbulence} + \pp{v}{t}_{Analysis} 
\]
\\

\noindent
{\bf 75)  \underline {DIABT} Total Diabatic Temperature Tendency  ($deg/day$) }

\noindent
{\bf DIABT} is the total time-tendency of Temperature due to Diabatic processes
and the Analysis forcing.
\begin{eqnarray*}
{\bf DIABT} & = & \pp{T}{t}_{Moist Processes} + \pp{T}{t}_{Shortwave Radiation} \\
           & + & \pp{T}{t}_{Longwave Radiation} + \pp{T}{t}_{Turbulence} + \pp{T}{t}_{Analysis} 
\end{eqnarray*}
\\
If we define the time-tendency of Temperature due to Diabatic processes as
\begin{eqnarray*}
\pp{T}{t}_{Diabatic} & = & \pp{T}{t}_{Moist Processes} + \pp{T}{t}_{Shortwave Radiation} \\
                     & + & \pp{T}{t}_{Longwave Radiation} + \pp{T}{t}_{Turbulence}
\end{eqnarray*}
then, since there are no surface pressure changes due to Diabatic processes, we may write
\[
\pp{T}{t}_{Diabatic} = {p^\kappa \over \pi }\pp{\pi \theta}{t}_{Diabatic}
\]
where $\theta = T/p^\kappa$.  Thus, {\bf DIABT} may be written as
\[
{\bf DIABT} = {p^\kappa \over \pi } \left( \pp{\pi \theta}{t}_{Diabatic} + \pp{\pi \theta}{t}_{Analysis} \right)
\]
\\

\noindent
{\bf 76)  \underline {DIABQ} Total Diabatic Specific Humidity Tendency  ($g/kg/day$) }

\noindent
{\bf DIABQ} is the total time-tendency of Specific Humidity due to Diabatic processes
and the Analysis forcing.
\[
{\bf DIABQ} = \pp{q}{t}_{Moist Processes} + \pp{q}{t}_{Turbulence} + \pp{q}{t}_{Analysis} 
\]
If we define the time-tendency of Specific Humidity due to Diabatic processes as
\[
\pp{q}{t}_{Diabatic} = \pp{q}{t}_{Moist Processes} + \pp{q}{t}_{Turbulence}
\]
then, since there are no surface pressure changes due to Diabatic processes, we may write
\[
\pp{q}{t}_{Diabatic} = {1 \over \pi }\pp{\pi q}{t}_{Diabatic}
\]
Thus, {\bf DIABQ} may be written as
\[
{\bf DIABQ} = {1 \over \pi } \left( \pp{\pi q}{t}_{Diabatic} + \pp{\pi q}{t}_{Analysis} \right)
\]
\\

\noindent
{\bf 77)  \underline {VINTUQ} Vertically Integrated Moisture Flux ($m/sec \cdot g/kg$) }

\noindent
The vertically integrated moisture flux due to the zonal u-wind is obtained by integrating
$u q$ over the depth of the atmosphere at each model timestep, 
and dividing by the total mass of the column.
\[
{\bf VINTUQ} = \frac{ \int_{surf}^{top} u q \rho dz  } { \int_{surf}^{top} \rho dz  }
\]
Using $\rho \delta z = -{\delta p \over g} = - {1 \over g} \delta p$, we have 
\[
{\bf VINTUQ} = { \int_0^1 u q dp  }
\]
\\


\noindent
{\bf 78)  \underline {VINTVQ} Vertically Integrated Moisture Flux ($m/sec \cdot g/kg$) }

\noindent
The vertically integrated moisture flux due to the meridional v-wind is obtained by integrating
$v q$ over the depth of the atmosphere at each model timestep, 
and dividing by the total mass of the column.
\[
{\bf VINTVQ} = \frac{ \int_{surf}^{top} v q \rho dz  } { \int_{surf}^{top} \rho dz  }
\]
Using $\rho \delta z = -{\delta p \over g} = - {1 \over g} \delta p$, we have 
\[
{\bf VINTVQ} = { \int_0^1 v q dp  }
\]
\\


\noindent
{\bf 79)  \underline {VINTUT} Vertically Integrated Heat Flux ($m/sec \cdot deg$) }

\noindent
The vertically integrated heat flux due to the zonal u-wind is obtained by integrating
$u T$ over the depth of the atmosphere at each model timestep, 
and dividing by the total mass of the column.
\[
{\bf VINTUT} = \frac{ \int_{surf}^{top} u T \rho dz  } { \int_{surf}^{top} \rho dz  }
\]
Or,
\[
{\bf VINTUT} = { \int_0^1 u T dp  }
\]
\\

\noindent
{\bf 80)  \underline {VINTVT} Vertically Integrated Heat Flux ($m/sec \cdot deg$) }

\noindent
The vertically integrated heat flux due to the meridional v-wind is obtained by integrating
$v T$ over the depth of the atmosphere at each model timestep, 
and dividing by the total mass of the column.
\[
{\bf VINTVT} = \frac{ \int_{surf}^{top} v T \rho dz  } { \int_{surf}^{top} \rho dz  }
\]
Using $\rho \delta z = -{\delta p \over g} $, we have 
\[
{\bf VINTVT} = { \int_0^1 v T dp  }
\]
\\

\noindent
{\bf 81 \underline {CLDFRC} Total 2-Dimensional Cloud Fracton ($0-1$) }

If we define the
time-averaged random and maximum overlapped cloudiness as CLRO and
CLMO respectively, then the probability of clear sky associated 
with random overlapped clouds at any level is (1-CLRO) while the probability of
clear sky associated with maximum overlapped clouds at any level is (1-CLMO). 
The total clear sky probability is given by (1-CLRO)*(1-CLMO), thus
the total cloud fraction at each  level may be obtained by 
1-(1-CLRO)*(1-CLMO).

At any given level, we may define the clear line-of-site probability by
appropriately accounting for the maximum and random overlap
cloudiness.  The clear line-of-site probability is defined to be
equal to the product of the clear line-of-site probabilities
associated with random and maximum overlap cloudiness.  The clear
line-of-site probability $C(p,p^{\prime})$ associated with maximum overlap clouds, 
from the current pressure $p$ 
to the model top pressure, $p^{\prime} = p_{top}$, or the model surface pressure, $p^{\prime} = p_{surf}$,
is simply 1.0 minus the largest maximum overlap cloud value along  the
line-of-site, ie.

$$1-MAX_p^{p^{\prime}} \left( CLMO_p \right)$$

Thus, even in the time-averaged sense it is assumed that the
maximum overlap clouds are correlated in the vertical.  The clear
line-of-site probability associated with random overlap clouds is
defined to be the product of the clear sky probabilities at each
level along the line-of-site, ie. 

$$\prod_{p}^{p^{\prime}} \left( 1-CLRO_p \right)$$

The total cloud fraction at a given level associated with a line-
of-site calculation is given by

$$1-\left( 1-MAX_p^{p^{\prime}} \left[ CLMO_p \right] \right)
    \prod_p^{p^{\prime}} \left( 1-CLRO_p \right)$$


\noindent
The 2-dimensional net cloud fraction as seen from the top of the
atmosphere is given by
\[
{\bf CLDFRC} = 1-\left( 1-MAX_{l=l_1}^{Nrphys} \left[ CLMO_l \right] \right)
    \prod_{l=l_1}^{Nrphys} \left( 1-CLRO_l \right)
\]
\\
For a complete description of cloud/radiative interactions, see Section \ref{sec:fizhi:radcloud}.


\noindent
{\bf 82)  \underline {QINT} Total Precipitable Water ($gm/cm^2$) }

\noindent
The Total Precipitable Water is defined as the vertical integral of the specific humidity,
given by:
\begin{eqnarray*}
{\bf QINT} & = & \int_{surf}^{top} \rho q dz \\
           & = & {\pi \over g} \int_0^1 q dp
\end{eqnarray*}
where we have used the hydrostatic relation 
$\rho \delta z = -{\delta p \over g} $.
\\


\noindent
{\bf 83)  \underline {U2M}  Zonal U-Wind at 2 Meter Depth ($m/sec$) }

\noindent
The u-wind at the 2-meter depth is determined from the similarity theory:
\[
{\bf U2M} = {u_* \over k} \psi_{m_{2m}} {u_{sl} \over {W_s}} =
{ \psi_{m_{2m}} \over {\psi_{m_{sl}} }}u_{sl}
\]

\noindent
where $\psi_m(2m)$ is the non-dimensional wind shear at two meters, and the subscript
$sl$ refers to the height of the top of the surface layer. If the roughness height
is above two meters, ${\bf U2M}$ is undefined.
\\
 
\noindent
{\bf 84)  \underline {V2M}  Meridional V-Wind at 2 Meter Depth ($m/sec$) }

\noindent
The v-wind at the 2-meter depth is a determined from the similarity theory:
\[
{\bf V2M} = {u_* \over k} \psi_{m_{2m}} {v_{sl} \over {W_s}} =
{ \psi_{m_{2m}} \over {\psi_{m_{sl}} }}v_{sl}
\]

\noindent
where $\psi_m(2m)$ is the non-dimensional wind shear at two meters, and the subscript
$sl$ refers to the height of the top of the surface layer. If the roughness height
is above two meters, ${\bf V2M}$ is undefined.
\\
 
\noindent
{\bf 85)  \underline {T2M}  Temperature at 2 Meter Depth ($deg \hspace{.1cm} K$) }

\noindent
The temperature at the 2-meter depth is a determined from the similarity theory:
\[
{\bf T2M} = P^{\kappa} ({\theta* \over k} ({\psi_{h_{2m}}+\psi_g}) + \theta_{surf} ) = 
P^{\kappa}(\theta_{surf} + { {\psi_{h_{2m}}+\psi_g} \over {{\psi_{h_{sl}}+\psi_g}} }
(\theta_{sl} - \theta_{surf})) 
\]
where:
\[
\theta_* = - { (\overline{w^{\prime}\theta^{\prime}}) \over {u_*} }
\]

\noindent
where $\psi_h(2m)$ is the non-dimensional temperature gradient at two meters, $\psi_g$ is
the non-dimensional temperature gradient in the viscous sublayer, and the subscript
$sl$ refers to the height of the top of the surface layer. If the roughness height
is above two meters, ${\bf T2M}$ is undefined.
\\
 
\noindent
{\bf 86)  \underline {Q2M}  Specific Humidity at 2 Meter Depth ($g/kg$) }

\noindent
The specific humidity at the 2-meter depth is determined from the similarity theory:
\[
{\bf Q2M} = P^{\kappa} ({q_* \over k} ({\psi_{h_{2m}}+\psi_g}) + q_{surf} ) = 
P^{\kappa}(q_{surf} + { {\psi_{h_{2m}}+\psi_g} \over {{\psi_{h_{sl}}+\psi_g}} }
(q_{sl} - q_{surf})) 
\]
where:
\[
q_* = - { (\overline{w^{\prime}q^{\prime}}) \over {u_*} }
\]

\noindent
where $\psi_h(2m)$ is the non-dimensional temperature gradient at two meters, $\psi_g$ is
the non-dimensional temperature gradient in the viscous sublayer, and the subscript
$sl$ refers to the height of the top of the surface layer. If the roughness height
is above two meters, ${\bf Q2M}$ is undefined.
\\
 
\noindent
{\bf 87)  \underline {U10M}  Zonal U-Wind at 10 Meter Depth ($m/sec$) }

\noindent
The u-wind at the 10-meter depth is an interpolation between the surface wind
and the model lowest level wind using the ratio of the non-dimensional wind shear
at the two levels:
\[
{\bf U10M} = {u_* \over k} \psi_{m_{10m}} {u_{sl} \over {W_s}} =
{ \psi_{m_{10m}} \over {\psi_{m_{sl}} }}u_{sl}
\]

\noindent
where $\psi_m(10m)$ is the non-dimensional wind shear at ten meters, and the subscript
$sl$ refers to the height of the top of the surface layer.
\\
 
\noindent
{\bf 88)  \underline {V10M}  Meridional V-Wind at 10 Meter Depth ($m/sec$) }

\noindent
The v-wind at the 10-meter depth is an interpolation between the surface wind
and the model lowest level wind using the ratio of the non-dimensional wind shear
at the two levels:
\[
{\bf V10M} = {u_* \over k} \psi_{m_{10m}} {v_{sl} \over {W_s}} =
{ \psi_{m_{10m}} \over {\psi_{m_{sl}} }}v_{sl}
\]

\noindent
where $\psi_m(10m)$ is the non-dimensional wind shear at ten meters, and the subscript
$sl$ refers to the height of the top of the surface layer.
\\
 
\noindent
{\bf 89)  \underline {T10M}  Temperature at 10 Meter Depth ($deg \hspace{.1cm} K$) }

\noindent
The temperature at the 10-meter depth is an interpolation between the surface potential 
temperature and the model lowest level potential temperature using the ratio of the 
non-dimensional temperature gradient at the two levels:
\[
{\bf T10M} = P^{\kappa} ({\theta* \over k} ({\psi_{h_{10m}}+\psi_g}) + \theta_{surf} ) = 
P^{\kappa}(\theta_{surf} + { {\psi_{h_{10m}}+\psi_g} \over {{\psi_{h_{sl}}+\psi_g}} }
(\theta_{sl} - \theta_{surf})) 
\]
where:
\[
\theta_* = - { (\overline{w^{\prime}\theta^{\prime}}) \over {u_*} }
\]

\noindent
where $\psi_h(10m)$ is the non-dimensional temperature gradient at two meters, $\psi_g$ is
the non-dimensional temperature gradient in the viscous sublayer, and the subscript
$sl$ refers to the height of the top of the surface layer.
\\
 
\noindent
{\bf 90)  \underline {Q10M}  Specific Humidity at 10 Meter Depth ($g/kg$) }

\noindent
The specific humidity at the 10-meter depth is an interpolation between the surface specific 
humidity and the model lowest level specific humidity using the ratio of the 
non-dimensional temperature gradient at the two levels:
\[
{\bf Q10M} = P^{\kappa} ({q_* \over k} ({\psi_{h_{10m}}+\psi_g}) + q_{surf} ) = 
P^{\kappa}(q_{surf} + { {\psi_{h_{10m}}+\psi_g} \over {{\psi_{h_{sl}}+\psi_g}} }
(q_{sl} - q_{surf})) 
\]
where:
\[
q_* =  - { (\overline{w^{\prime}q^{\prime}}) \over {u_*} }
\]

\noindent
where $\psi_h(10m)$ is the non-dimensional temperature gradient at two meters, $\psi_g$ is
the non-dimensional temperature gradient in the viscous sublayer, and the subscript
$sl$ refers to the height of the top of the surface layer.
\\
 
\noindent
{\bf 91)  \underline {DTRAIN} Cloud Detrainment Mass Flux ($kg/m^2$) } 

The amount of cloud mass moved per RAS timestep at the cloud detrainment level is written:
\[
{\bf DTRAIN} = \eta_{r_D}m_B
\]
\noindent
where $r_D$ is the detrainment level, 
$m_B$ is the cloud base mass flux, and $\eta$
is the entrainment, defined in Section \ref{sec:fizhi:mc}.
\\

\noindent
{\bf 92)  \underline {QFILL}  Filling of negative Specific Humidity ($g/kg/day$) }

\noindent
Due to computational errors associated with the numerical scheme used for
the advection of moisture, negative values of specific humidity may be generated.  The
specific humidity is checked for negative values after every dynamics timestep.  If negative
values have been produced, a filling algorithm is invoked which redistributes moisture from
below.  Diagnostic {\bf QFILL} is equal to the net filling needed
to eliminate negative specific humidity, scaled to a per-day rate:
\[
{\bf QFILL} = q^{n+1}_{final} - q^{n+1}_{initial}
\]
where
\[
q^{n+1} = (\pi q)^{n+1} / \pi^{n+1}
\]

\subsection{Dos and Donts}

\subsection{Diagnostics Reference}
