s_outp_pkgs/text/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 in this section of the manual is a listing of
diagnostic quantities available within the main (dynamics) part of the
GCM. Additional diagnostic quantities, defined within the different
GCM packages, are available and are listed in the diagnostic menu
subsection of the manual section associated with each relevant
package. 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
referred 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.  Quantities
referred 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
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 DIAGNOSTICS\_ADDTOLIST} 
(\filelink{pkg/diagnostics/diagnostics\_addtolist.F}{pkg-diagnostics-diagnostics\_addtolist.F}):
%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 routine for defining a new permanent 
diagnostic in the main model or in a package. The calling sequence is:

\begin{verbatim}
       CALL DIAGNOSTICS_ADDTOLIST (
     O     diagNum,
     I     diagName, diagCode, diagUnits, diagTitle, diagMate,
     I     myThid )

     where:
       diagNum   = diagnostic Id number - Output from routine
       diagName  = name of diagnostic to declare
       diagCode  = parser code for this diagnostic
       diagUnits = field units for this diagnostic
       diagTitle = field description for this diagnostic
       diagMate  = diagnostic mate number
       myThid    = my Thread Id number

\end{verbatim}


\noindent
{\bf DIAGNOSTICS\_FILL} 
(\filelink{pkg/diagnostics/diagnostics\_fill.F}{pkg-diagnostics-diagnostics\_fill.F}):
This is the main user interface routine to the diagnostics package. 
This routine will increment the specified
diagnostic quantity with a field sent through the argument list.

\begin{verbatim}
        CALL DIAGNOSTICS_FILL(
       I             inpFld, chardiag,
       I             kLev, nLevs, bibjFlg, bi, bj, myThid )

     where:
        inpFld   = Field to increment diagnostics array
        diagName = diagnostic identificator name (8 characters long)
        kLev     = Integer flag for vertical levels:
                   > 0 (any integer): WHICH single level to increment in qdiag.
                   0,-1 to increment "nLevs" levels in qdiag,
                   0 : fill-in in the same order as the input array
                   -1: fill-in in reverse order.
        nLevs    = indicates Number of levels of the input field array
                   (whether to fill-in all the levels (kLev<1) or just one (kLev>0))
        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 (with overlap regions)
                      (local array here means that it has no bi-bj dimensions)
                 = 3 indicates that the bi-bj loop is done OUTSIDE
                      AND that we have been sent a local array
                      AND that the array has no overlap region (interior only)
                   NOTE - bibjFlg can be NEGATIVE to indicate not to increment counter
        bi       = X-direction tile number - used for bibjFlg=1-3
        bj       = Y-direction tile number - used for bibjFlg=1-3
        myThid   =  my thread Id number
\end{verbatim}

\noindent
{\bf DIAGNOSTICS\_SCALE\_FILL} 
(\filelink{pkg/diagnostics/diagnostics\_scale\_fill.F}{pkg-diagnostics-diagnostics\_scale\_fill.F}):
This is a possible alternative routine to
DIAGNOSTICS\_FILL which performs the same functions and has an additional option
to scale the field before filling or raise the field to a power before filling.

\begin{verbatim}
        CALL DIAGNOSTICS_SCALE_FILL(
       I             inpFld, scaleFact, power, chardiag,
       I             kLev, nLevs, bibjFlg, bi, bj, myThid )


     where all the arguments are the same as for DIAGNOSTICS_FILL with 
     the addition of:
        scaleFact   = Scaling factor to apply to the input field product
        power       = Integer power to which to raise the input field (before scaling)
\end{verbatim}

\noindent
{\bf DIAGNOSTICS\_IS\_ON}: Function call to inquire whether a
diagnostic is active and can be incremented. Useful when there is a
computation that must be done locally before a call to
DIAGNOSTICS\_FILL. The call sequence:

\begin{verbatim}
        flag = DIAGNOSTICS_IS_ON( diagName, myThid )

     where:
        diagName = diagnostic identificator name (8 characters long)
        myThid   = my thread Id number
\end{verbatim}

\noindent
{\bf DIAGNOSTICS\_GET\_POINTERS}
(\filelink{pkg/diagnostics/diagnostics\_utils.F}{pkg-diagnostics-diagnostics\_utils.F}):
This subroutine retrieves the value of a the diagnostics pointers 
that other routines require as input - can be useful if the diagnostics common 
blocks are not local to a routine.

\begin{verbatim}
        CALL DIAGNOSTICS_GET_POINTERS(
     I                       diagName, listId,
     O                       ndId, ip,
     I                       myThid )

     where:
        diagName = diagnostic identificator name (8 characters long)
        listId   = list number that specifies the output frequency
        ndId     = diagnostics  Id number (in available diagnostics list)
        ip       = diagnostics  pointer to storage array
        myThid   = my Thread Id number
\end{verbatim}

\noindent
{\bf GETDIAG} 
(\filelink{pkg/diagnostics/diagnostics\_utils.F}{pkg-diagnostics-diagnostics\_utils.F}):
This subroutine retrieves the value of a model diagnostic.  
This routine is particularly 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:

\begin{verbatim}
        CALL GETDIAG(
       I             levreal, undef,
       O             qtmp,
       I             ndId, mate, ip, im, bi, bj, myThid )

     where:
        lev     = Model Level at which the diagnostic is desired
        undef   = Fill value to be used when diagnostic is undefined
        qtmp    = Time-Averaged Diagnostic Output
        ndId    = diagnostics  Id number (in available diagnostics list)
        mate    = counter diagnostic number if any ; 0 otherwise
        ip      = pointer to storage array location for diag.
        im      = pointer to storage array location for mate
        bi      = X-direction tile number
        bj      = Y-direction tile number
        myThid  = my thread Id number
\end{verbatim}

\noindent
{\bf DIAGNOSTICS\_CLRDIAG} 
(\filelink{pkg/diagnostics/diagnostics\_clear.F}{pkg-diagnostics-diagnostics\_clear.F}):
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:

\begin{verbatim}
        CALL DIAGNOSTICS_CLRDIAG ( ipt, nLev, myThid )

     where:
        ipt    :: diagnostic pointer to storage array
        nLev   :: number of levels (in storage array) to reset
        myThid :: my Thread Id number
\end{verbatim}

\noindent
The diagnostics are computed at various times and places within the
GCM. Because MITgcm 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 {\tt diagnostics} in file {\tt
  DIAGNOSTICS.h}) is internally defined as a character*10 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$   &  Used for Level Integrated output: cumulate levels \\
           & $\rightarrow$ r &  same but cumulate product by model level thickness \\
           & $\rightarrow$ R &  same but cumulate product by hFac \& level thickness \\ \hline
 parse(4)  & $\rightarrow$ P &  Positive Definite Diagnostic      \\ \hline
 parse(5)  & $\rightarrow$ C &  with Counter array                \\
           & $\rightarrow$ D &  Disabled Diagnostic for output    \\ \hline
 parse(6-8)&                 &  retired, formerly: 3-digit mate number \\ \hline
 parse(9)  & $\rightarrow$ U &  model-level plus 1/2  \\
           & $\rightarrow$ M &  model-level middle  \\
           & $\rightarrow$ L &  model-level minus 1/2  \\ \hline
 parse(10) & $\rightarrow$ 0 &  levels = 0  \\
           & $\rightarrow$ 1 &  levels = 1  \\
           & $\rightarrow$ R &  levels = Nr  \\
           & $\rightarrow$ L &  levels = MAX(Nr,NrPhys)  \\
           & $\rightarrow$ M &  levels = MAX(Nr,NrPhys) - 1  \\
           & $\rightarrow$ G &  levels = Ground\_level Number \\
           & $\rightarrow$ I &  levels = sea-Ice\_level Number \\
           & $\rightarrow$ X &  free levels option (need to be set explicitly)\\ \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 ``UUR\hspace{5mm}MR''.  From GDIAG we can determine that this diagnostic is a 
U-vector component located at the C-grid U-point, model mid-level (M)
with Nr levels (last R).

\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.

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

\noindent
To use the diagnostics package, other than enabling it in {\tt packages.conf}
and turning the {\tt useDiagnostics} flag in {\tt data.pkg} to .TRUE., there are two
further steps the user must take to enable the diagnostics package for
output of quantities that are already defined in the GCM under an experiment's
configuration of packages.  
A parameter file {\tt data.diagnostics} must be supplied in the run directory,
and the file DIAGNOSTICS\_SIZE.h must be included in the 
code directory.  The steps for defining a new (permanent or experiment-specific 
temporary) diagnostic quantity will be outlined later. 

\noindent The namelist in parameter file {\tt data.diagnostics} will activate 
a user-defined list of diagnostics quantities to be computed, 
specify the frequency and type of output, the number of levels, and 
the name of all the separate output files. 
A sample {\tt data.diagnostics} namelist file:

\begin{verbatim}
# Diagnostic Package Choices
#--------------------
#  dumpAtLast (logical): always write output at the end of simulation (default=F)
#  diag_mnc   (logical): write to NetCDF files (default=useMNC)
#--for each output-stream:
#  fileName(n) : prefix of the output file name (max 80c long) for outp.stream n
#  frequency(n):< 0 : write snap-shot output every |frequency| seconds
#               > 0 : write time-average output every frequency seconds
#  timePhase(n)     : write at time = timePhase + multiple of |frequency|
#    averagingFreq  : frequency (in s) for periodic averaging interval
#    averagingPhase : phase     (in s) for periodic averaging interval
#    repeatCycle    : number of averaging intervals in 1 cycle
#  levels(:,n) : list of levels to write to file (Notes: declared as REAL)
#                when this entry is missing, select all common levels of this list
#  fields(:,n) : list of selected diagnostics fields (8.c) in outp.stream n
#                (see "available_diagnostics.log" file for the full list of diags)
#  missing_value(n) : missing value for real-type fields in output file "n"
#  fileFlags(n)     : specific code (8c string) for output file "n"
#--------------------
 &DIAGNOSTICS_LIST
  fields(1:2,1) = 'UVEL    ','VVEL    ', 
   levels(1:5,1) = 1.,2.,3.,4.,5., 
   filename(1) = 'diagout1', 
  frequency(1) = 86400., 
  fields(1:3,2) = 'THETA   ','SALT    ',
   filename(2) = 'diagout2', 
  fileflags(2) = ' P1     ', 
  frequency(2) = 3600., 
 &

#--------------------
# Parameter for Diagnostics of per level statistics:
#--------------------
#  diagSt_mnc (logical): write stat-diags to NetCDF files (default=diag_mnc)
#  diagSt_regMaskFile : file containing the region-mask to read-in
#  nSetRegMskFile   : number of region-mask sets within the region-mask file
#  set_regMask(i)   : region-mask set-index that identifies the region "i"
#  val_regMask(i)   : region "i" identifier value in the region mask
#--for each output-stream:
#  stat_fName(n) : prefix of the output file name (max 80c long) for outp.stream n
#  stat_freq(n):< 0 : write snap-shot output every |stat_freq| seconds
#               > 0 : write time-average output every stat_freq seconds
#  stat_phase(n)    : write at time = stat_phase + multiple of |stat_freq|
#  stat_region(:,n) : list of "regions" (default: 1 region only=global)
#  stat_fields(:,n) : list of selected diagnostics fields (8.c) in outp.stream n
#                (see "available_diagnostics.log" file for the full list of diags)
#--------------------
 &DIAG_STATIS_PARMS
 &
\end{verbatim}

\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 86400. seconds,
(frequency is 86400.), and will write fields which are multiple-level
fields at 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 3600. seconds, includes fields
with all levels, and the names of diagnostics quantities are THETA and SALT.

\noindent
The user must assure that enough computer memory is allocated for the
diagnostics and the output streams selected for a particular
experiment.  This is accomplished by modifying the file
DIAGNOSTICS\_SIZE.h and including it in the experiment code directory.
The parameters that should be checked are called numDiags, numLists,
numperList, and diagSt\_size.

\noindent numDiags (and diagSt\_size): \\
\noindent 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:\\
\begin{verbatim}
      _RL  qdiag(1-Olx,sNx+Olx,1-Olx,sNx+Olx,numDiags,nSx,nSy)
      _RL  qSdiag(0:nStats,0:nRegions,diagSt_size,nSx,nSy)
      COMMON / DIAG_STORE_R / qdiag, qSdiag
\end{verbatim}
\noindent
The first two-dimensions of qdiag correspond to the horizontal
dimension of a given diagnostic, and 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 Nr. In order for the User to
enable more than 1 three-dimensional diagnostic, the size of the
diagnostics common must be expanded to accommodate 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.  For the above example,
there are 4 multiple level fields, which the diagnostics menu (see
below) indicates are defined at the GCM vertical resolution, Nr. The
value of numDiags in DIAGNOSTICS\_SIZE.h would therefore be equal to
4*Nr, or, say 40 if $Nr=10$.

\noindent numLists and numperList: \\
\noindent The parameter numLists must be set greater than or equal to
the number of separate output streams that the user specifies in the
namelist file data.diagnostics.  The parameter numperList corresponds
to the maximum number of diagnostics requested per output streams.

\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
{\tt data.diagnostics}. This is accomplished by adding one of the ``User
Diagnostic'' field names (UDIAG1 through UDIAG10, for multi-level
fields, or SDIAG1 through SDIAG10 for single level fields) to the
data.diagnostics namelist in one of the output streams. These fields
are listed in the diagnostics menu. The second step is to add a call
to DIAGNOSTICS\_FILL from the subroutine in which the quantity desired
for diagnostic output is computed.

\noindent
In order to add a new diagnostic to the permanent set of diagnostics
that the main model or any package contains as part of its diagnostics
menu, the subroutine DIAGNOSTICS\_ADDTOLIST should be called during the
initialization phase of the main model or package. For the main model,
the call should be made from subroutine DIAGNOSTICS\_MAIN\_INIT, and
for a package, the call should probably be made from 
%somewhere in the PACKAGES\_INIT\_FIXED sequence (probably 
from inside the particular package's init\_fixed routine. 
A typical code sequence to set the
input arguments to DIAGNOSTICS\_ADDTOLIST would look like:

\begin{verbatim}
      diagName  = 'RHOAnoma'
      diagTitle = 'Density Anomaly (=Rho-rhoConst)'
      diagUnits = 'kg/m^3          '
      diagCode  = 'SMR     MR      '
      CALL DIAGNOSTICS\_ADDTOLIST( diagNum,
     I          diagName, diagCode, diagUnits, diagTitle, 0, myThid )
\end{verbatim}

\noindent If the new diagnostic quantity is associated with either a
vector pair or a diagnostic counter, the diagMate argument must be
provided with the proper index corresponding to the ``mate''. 
The output argument from DIAGNOSTICS\_ADDTOLIST that is called diagNum here 
contains a running total of the number of diagnostics defined in the code up to
any point during the run. The sequence number for the next two
diagnostics defined (the two components of the vector pair, for
instance) will be diagNum+1 and diagNum+2. The definition of the first
component of the vector pair must fill the ``mate'' segment of the
diagCode as diagnostic number diagNum+2.  Since the subroutine
increments diagNum, the definition of the second component of the
vector fills the ``mate'' part of diagCode with diagNum. A code
sequence for this case would look like:

\begin{verbatim}
      diagName  = 'UVEL    '
      diagTitle = 'Zonal Component of Velocity (m/s)'
      diagUnits = 'm/s             '
      diagCode  = 'UUR     MR      '
      diagMate  = diagNum + 2
      CALL DIAGNOSTICS_ADDTOLIST( diagNum,
     I   diagName, diagCode, diagUnits, diagTitle, diagMate, myThid )

      diagName  = 'VVEL    '
      diagTitle = 'Meridional Component of Velocity (m/s)'
      diagUnits = 'm/s             '
      diagCode  = 'VVR     MR      '
      diagMate  = diagNum
      CALL DIAGNOSTICS_ADDTOLIST( diagNum,
     I   diagName, diagCode, diagUnits, diagTitle, diagMate, myThid )
\end{verbatim}

\input{s_outp_pkgs/text/diagnostics-menu.tex}

\newpage
\noindent For a list of the diagnostic fields available in the
different MITgcm packages, follow the link to the diagnostics menu
in the manual section describing the package:

\begin{itemize}
\item aim: \ref{sec:pkg:aim:diagnostics}
\item exf: \ref{sec:pkg:exf:diagnostics}
\item gchem: \ref{sec:pkg:gchem:diagnostics}
\item generic\_advdiff: \ref{sec:pkg:gad:diagnostics}
\item gridalt: \ref{sec:pkg:gridalt:diagnostics}
\item gmredi: \ref{sec:pkg:gmredi:diagnostics}
\item fizhi: \ref{sec:pkg:fizhi:diagnostics}
\item kpp: \ref{sec:pkg:kpp:diagnostics}
\item land: \ref{sec:pkg:land:diagnostics}
\item mom\_common: \ref{sec:pkg:mom_common:diagnostics}
\item obcs: \ref{sec:pkg:obcs:diagnostics}
\item thsice: \ref{sec:pkg:thsice:diagnostics}
\item shap\_filt: \ref{sec:pkg:shap_filt:diagnostics}
\item ptracers: \ref{sec:pkg:ptracers:diagnostics}
\end{itemize}

\subsection{Dos and Donts}

\subsection{Diagnostics Reference}
