Commits

Anonymous committed 4f6202f Merge

Merge with pflotran-dev.

Comments (0)

Files changed (23)

 *.stdout
 *.regression
 regression_tests/default/anisothermal/thc_1d.h5
+*.nav
+*.snm
+*.vrb
+shortcourse/exercises/CLM-CN/doc/CLM-CN.pdf

docs/user_manual/governing_eqns.tex

 \end{tabular}
 \end{table}
 
-\section{Method of Solution}
+\section{Appendix: Method of Solution}
 
 \setcounter{equation}{0}
 

docs/user_manual/installation.tex

 \toprule
 Option & \multicolumn{1}{c}{Description}\\
 \midrule
-{\tt -pflotranin} <input file name> & specify input file [Default: {\tt pflotran.in}]\\
-{\tt -file\_prefix} <output file name> & specify output file [Default: {\tt pflotran.out}]\\
-{\tt -input\_prefix} & \\
-{\tt -output\_prefix} & \\
+{\tt -pflotranin <string>} & specify input file [Default: {\tt pflotran.in}]\\
+{\tt -input\_prefix <string>} & specify input (and output) file prefix [default: {\tt pflotran}]\\
+{\tt -output\_prefix <string>} & specify output file prefix [default: {\tt pflotran}]\\
 {\tt -screen\_output off} & turn off screen output\\
 %{\tt -v} &\\
 {\tt -realization\_id <integer>} & run specified realization ID\\

docs/user_manual/local.tex

 \usepackage[table]{xcolor}
 %\usepackage{toc_entr}
 
-%\usepackage{verbatim}
+\usepackage{verbatim}
 %\usepackage{moreverb}
 %\usepackage{sverb}
 \usepackage{fancyvrb}
+\usepackage{enumitem}
 
 \usepackage[framemethod=tikz,outerlinewidth=2,outerlinecolor=teal,backgroundcolor=lightteal]{mdframed}
 

docs/user_manual/reaction_sandbox.tex

+\section{Appendix: Reaction Sandbox}
+
+\setcounter{equation}{0}
+
+\subsection{Background}
+
+Researchers often have a suite of reactions tailored to a unique problem scenario, but these reaction networks only exist in their respective research codes. The ``reaction sandbox'' provides these researchers with a venue for implementing user-defined reactions within PFLOTRAN. Reaction networks developed within the reaction sandbox can leverage existing biogeochemical capability within in PFLOTRAN (e.g. equilibrium aqueous complexation, mineral precipitation--dissolution, etc.) or function independently. Please note that although the reaction sandbox facilitates the integration of user-defined reactions, the process still requires a basic understanding of PFLOTRAN and its approach to solving reaction through the Newton-Raphson method. For instance,  One must understand the purpose and function of the rt\_auxvar and global\_auxvar objects.
+
+\subsection{Implementation}
+The core framework of reaction sandbox leverages Fortran 2003 object--oriented extendable derived types and methods and consists of two modules:
+
+\begin{itemize}
+  \item[] Reaction\_Sandbox\_module (reaction\_sandbox.F90)
+  \item[] Reaction\_Sandbox\_Base\_class (reaction\_sandbox\_base.F90).
+\end{itemize}
+
+To implement a new reaction within the reaction sandbox, one creates a new class by extending the Reaction\_Sandbox\_Base\_class and adds the new class to the Reaction\_Sandbox\_module.  The following steps illustrate this process through the creation of the class Reaction\_Sandbox\_Example\_class that implements a first order decay reaction.
+
+\begin{enumerate}
+  \item Copy reaction\_sandbox\_template.F90 to a new filename (e.g. reaction\_sandbox\_example.F90).
+  \item Replace all references to Template/template with the new reaction name.
+    \begin{itemize}
+      \item[] \verb+Template+ $\rightarrow$ \verb+Example+
+      \item[] \verb+template+ $\rightarrow$ \verb+example+
+    \end{itemize}
+  \item Add necessary variables to the module and/or the extended derived type.
+\begin{verbatim}
+character(len=MAXWORDLENGTH) :: species_name
+PetscInt :: species_id
+PetscReal :: rate_constant
+\end{verbatim}
+  \item Add the necessary functionality within the following subroutines:
+    \begin{enumerate}[label=\alph*]
+      \item ExampleCreate: Allocate the reaction object, initializing all variables to zero and nullifying arrays.  {\bf Be sure to nullify ExampleCreate\%next which comes from the base class.}  E.g.,
+  \begin{verbatim}
+  allocate(ExampleCreate)
+  ExampleCreate%species_name = ''
+  ExampleCreate%species_id = 0
+  ExampleCreate%rate_constant = 0.d0
+  nullify(ExampleCreate%next)
+  \end{verbatim}
+      \item ExampleRead: Read parameters in from the input file block \verb+EXAMPLE+.  E.g.,
+  \begin{verbatim}
+  ...
+  case('SPECIES_NAME')
+    call InputReadWord(input,option,this%species_name, &
+                       PETSC_TRUE)
+    call InputErrorMsg(input,option,'species_name', &
+                   'CHEMISTRY,REACTION_SANDBOX,EXAMPLE')
+  ...
+  \end{verbatim}
+      \item ExampleSetup: Construct the reaction network (e.g. array allocation, establishing linkages, etc.). E.g.,
+  \begin{verbatim}
+  ...
+  this%species_id = &
+    GetPrimarySpeciesIDFromName(this%species_name, &
+                                reaction,option)
+  ...
+  \end{verbatim}
+      \item ExampleReact: Calculate contribution of reaction to the residual (units = moles/sec) and Jacobian (units = kg water/sec). E.g.,
+  \begin{verbatim}
+  ...
+  Residual(this%species_id) = &
+    Residual(this%species_id) - &
+    this%rate_constant*porosity* &
+    global_auxvar%sat(iphase)*volume*1.d3* &
+    rt_auxvar%total(this%species_id,iphase)
+  ...
+  Jacobian(this%species_id,this%species_id) = &
+  Jacobian(this%species_id,this%species_id) + &
+    this%rate_constant*porosity* &
+    global_auxvar%sat(iphase)*volume*1.d3*
+    rt_auxvar%aqueous%dtotal(this%species_id, &
+                             this%species_id,iphase)
+  ...
+  \end{verbatim}
+      \item ExampleDestroy: Deallocate any dynamic memory within the class (without deallocating the object itself).
+    \end{enumerate}
+  \item Ensure that the methods within the extended derived type point to the proper procedures in the module
+  \begin{verbatim}
+  procedure, public :: ReadInput => ExampleRead
+  procedure, public :: Setup => ExampleSetup
+  procedure, public :: Evaluate => ExampleReact
+  procedure, public :: Destroy => ExampleDestroy
+  \end{verbatim}
+  \item Within reaction\_sandbox.F90:
+    \begin{enumerate}[label=\alph*]
+      \item Add Reaction\_Sandbox\_Example\_class+ to the list of modules to be ``used'' at the top of the file.
+      \item Add a case statement in RSandboxRead2 for the keyword defining the new reaction and create the reaction within.  I.e.
+  \begin{verbatim}
+    case('EXAMPLE')
+      new_sandbox => ExampleCreate()
+  \end{verbatim}
+    \end{enumerate}
+\end{enumerate}
+

docs/user_manual/user_manual.tex

 
 \include{governing_eqns}
 
+\include{reaction_sandbox}
+
 \end{document}
 

regression_tests/default/infiltrometer/infiltrometer.cfg

 saturation = 1.0e-12 absolute
 
 [16m]
+pressure = 1.0e-9 absolute

shortcourse/exercises/CLM-CN/CLM-CN.in

+:Description: Batch CLM-CN simulation
+: Rate constants, C/N ratios, respiration fractions from Bonan et al., 2012
+
+:=========================== useful tranport parameters =======================
+UNIFORM_VELOCITY 0.d0 0.d0 0.d0 
+
+REFERENCE_DENSITY 1.d3
+
+:=========================== chemistry ========================================
+CHEMISTRY
+  PRIMARY_SPECIES
+    A(aq)
+  /
+  IMMOBILE_SPECIES
+    N
+    C
+    SOM1
+    SOM2
+    SOM3
+    SOM4
+    LabileC
+    CelluloseC
+    LigninC
+    LabileN
+    CelluloseN
+    LigninN
+  /
+  REACTION_SANDBOX
+    CLM-CN
+      POOLS   ! CN ratio
+        SOM1  12.d0 
+        SOM2  12.d0
+        SOM3  10.d0
+        SOM4  10.d0
+        Labile
+        Cellulose
+        Lignin
+      /
+      REACTION
+        UPSTREAM_POOL Labile
+        DOWNSTREAM_POOL SOM1
+        TURNOVER_TIME 20. h
+        RESPIRATION_FRACTION 0.39d0
+        N_INHIBITION 1.d-10
+      /
+      REACTION
+        UPSTREAM_POOL Cellulose
+        DOWNSTREAM_POOL SOM2
+        TURNOVER_TIME 14. d
+        RESPIRATION_FRACTION 0.55
+        N_INHIBITION 1.d-10
+      /
+      REACTION
+        UPSTREAM_POOL Lignin
+        DOWNSTREAM_POOL SOM3
+        TURNOVER_TIME 71. d
+        RESPIRATION_FRACTION 0.29d0
+        N_INHIBITION 1.d-10
+      /
+      REACTION
+        UPSTREAM_POOL SOM1
+        DOWNSTREAM_POOL SOM2
+        TURNOVER_TIME 14. d
+        RESPIRATION_FRACTION 0.28d0
+        N_INHIBITION 1.d-10
+      /
+      REACTION
+        UPSTREAM_POOL SOM2
+        DOWNSTREAM_POOL SOM3
+        TURNOVER_TIME 71. d
+        RESPIRATION_FRACTION 0.46d0
+        N_INHIBITION 1.d-10
+      /
+      REACTION
+        UPSTREAM_POOL SOM3
+        DOWNSTREAM_POOL SOM4
+        TURNOVER_TIME 2. y
+        RESPIRATION_FRACTION 0.55d0
+        N_INHIBITION 1.d-10
+      /
+      REACTION
+        UPSTREAM_POOL SOM4
+        TURNOVER_TIME 27.4 y
+        RESPIRATION_FRACTION 1.d0
+        N_INHIBITION 1.d-10
+      /
+    /
+  /
+:  LOG_FORMULATION
+  DATABASE ./CLM-CN_database.dat
+  OUTPUT
+    all
+  /
+END
+
+:=========================== solver options ===================================
+LINEAR_SOLVER TRANSPORT
+  SOLVER DIRECT
+END
+
+NEWTON_SOLVER TRANSPORT
+  ATOL 1.d-50
+:  RTOL 1.d-12
+END
+
+:=========================== discretization ===================================
+GRID
+  TYPE structured
+  NXYZ 1 1 1
+  BOUNDS
+    0.d0 1.d0
+    0.d0 1.d0
+    0.d0 1.d0
+  /
+END
+
+:=========================== fluid properties =================================
+FLUID_PROPERTY 
+  DIFFUSION_COEFFICIENT 1.d-9
+END
+
+:=========================== material properties ==============================
+MATERIAL_PROPERTY soil1
+  ID 1
+  POROSITY 0.25d0
+  TORTUOSITY 1.d0
+END
+
+:=========================== output options ===================================
+OUTPUT
+  PERIODIC_OBSERVATION TIMESTEP 1
+END
+
+:=========================== times ============================================
+TIME
+  FINAL_TIME 400.d0 d
+  INITIAL_TIMESTEP_SIZE 1.d0 h
+  MAXIMUM_TIMESTEP_SIZE 10.d0 d
+END
+
+:=========================== regions ==========================================
+REGION all
+  COORDINATES
+    0.d0 0.d0 0.d0
+    1.d0 1.d0 1.d0
+  /
+END
+
+REGION obs_pt
+  COORDINATE 0.5 0.5 0.5
+END
+
+:=========================== observation points ===============================
+OBSERVATION
+  REGION obs_pt
+END
+
+:=========================== transport conditions =============================
+TRANSPORT_CONDITION initial
+  TYPE zero_gradient
+  CONSTRAINT_LIST
+    0.d0 initial
+  /
+END
+
+:=========================== transport constraints ============================
+CONSTRAINT initial
+  CONCENTRATIONS
+    A(aq)  1.d-40     T  ! moles/L
+  /
+  IMMOBILE
+    N     1.d-6  ! moles/m^3
+    C     1.d-6
+    SOM1  1.d-10 ! moles C/m^3
+    SOM2  1.d-10
+    SOM3  1.d-10
+    SOM4  1.d-10
+    LabileC     0.1852d-3
+    CelluloseC  0.4578d-3
+    LigninC     0.2662d-3
+    LabileN     0.00508954d-3
+    CelluloseN  0.01258096d-3
+    LigninN     0.00731553d-3
+  /
+END
+
+:=========================== condition couplers ===============================
+: initial condition
+INITIAL_CONDITION
+  TRANSPORT_CONDITION initial
+  REGION all
+END
+
+:=========================== stratigraphy couplers ============================
+STRATA
+  REGION all
+  MATERIAL soil1
+END
+
+

shortcourse/exercises/CLM-CN/CLM-CN.py

+import sys
+import os
+try:
+  pflotran_dir = os.environ['PFLOTRAN_DIR']
+except KeyError:
+  print('PFLOTRAN_DIR must point to PFLOTRAN installation directory and be defined in system environment variables.')
+  sys.exit(1)
+sys.path.append(pflotran_dir + '/src/python')
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+import math
+import pflotran as pft
+
+filename = 'CLM-CN-obs-0.tec'
+
+f = plt.figure(figsize=(14,10))
+f.suptitle("CLM-CN",fontsize=16)
+plt.subplot(2,2,1)
+plt.xlabel('Time [y]')
+plt.ylabel('Concentration [mol/m^3]')
+
+#plt.xlim(0.,1.)
+#plt.ylim(4.8,8.2)
+#plt.yscale('log')
+#plt.grid(True)
+
+data = pft.Dataset(filename,1,4)
+
+# 12 g C / 1 g N
+CN_ratio_12 = 12./12. / (1. / 14.)
+# 10 g C / 1 g N
+CN_ratio_10 = 10./12. / (1. / 14.)
+
+# if is a dummy array 
+new_array = pft.Dataset(filename,1,2).get_array('y')
+
+arrayN = pft.Dataset(filename,1,3).get_array('y')
+arrayC = pft.Dataset(filename,1,4).get_array('y')
+arraySOM1 = pft.Dataset(filename,1,5).get_array('y')
+arraySOM2 = pft.Dataset(filename,1,6).get_array('y')
+arraySOM3 = pft.Dataset(filename,1,7).get_array('y')
+arraySOM4 = pft.Dataset(filename,1,8).get_array('y')
+arrayLabileC = pft.Dataset(filename,1,9).get_array('y')
+arrayCelluloseC = pft.Dataset(filename,1,10).get_array('y')
+arrayLigninC = pft.Dataset(filename,1,11).get_array('y')
+arrayLabileN = pft.Dataset(filename,1,12).get_array('y')
+arrayCelluloseN = pft.Dataset(filename,1,13).get_array('y')
+arrayLigninN = pft.Dataset(filename,1,14).get_array('y')
+
+for i in range(len(new_array)):
+  new_array[i] = arrayLabileC[i] + arrayLabileN[i]
+plt.plot(data.get_array('x'),new_array,label='Labile')
+for i in range(len(new_array)):
+  new_array[i] = arrayCelluloseC[i] + arrayCelluloseN[i]
+plt.plot(data.get_array('x'),new_array,label='Cellulose')
+for i in range(len(new_array)):
+  new_array[i] = arrayLigninC[i] + arrayLigninN[i]
+plt.plot(data.get_array('x'),new_array,label='Lignin')
+
+#'best'         : 0, (only implemented for axis legends)
+#'upper right'  : 1,
+#'upper left'   : 2,
+#'lower left'   : 3,
+#'lower right'  : 4,
+#'right'        : 5,
+#'center left'  : 6,
+#'center right' : 7,
+#'lower center' : 8,
+#'upper center' : 9,
+#'center'       : 10,
+plt.legend(loc=1)
+# xx-small, x-small, small, medium, large, x-large, xx-large, 12, 14
+plt.setp(plt.gca().get_legend().get_texts(),fontsize='small')
+#plt.setp(plt.gca().get_legend().get_texts(),linespacing=0.)
+plt.setp(plt.gca().get_legend().get_frame().set_fill(False))
+plt.setp(plt.gca().get_legend().draw_frame(False))
+#plt.gca().yaxis.get_major_formatter().set_powerlimits((-1,1))
+
+# SOM
+plt.subplot(2,2,2)
+plt.xlabel('Time [y]')
+plt.ylabel('Concentration [mol/m^3]')
+plt.plot(data.get_array('x'),arraySOM1,label='SOM1')
+plt.plot(data.get_array('x'),arraySOM2,label='SOM2')
+plt.plot(data.get_array('x'),arraySOM3,label='SOM3')
+plt.plot(data.get_array('x'),arraySOM4,label='SOM4')
+plt.legend(loc=1)
+plt.setp(plt.gca().get_legend().get_texts(),fontsize='small')
+plt.setp(plt.gca().get_legend().get_frame().set_fill(False))
+plt.setp(plt.gca().get_legend().draw_frame(False))
+
+# C
+plt.subplot(2,2,3)
+plt.xlabel('Time [y]')
+plt.ylabel('Concentration [mol/m^3]')
+for i in range(len(arrayC)):
+  new_array[i] = arrayC[i] + arrayLabileC[i]  + \
+                 arrayCelluloseC[i] + arrayLigninC[i] + \
+                 (arraySOM1[i] + arraySOM2[i] + \
+                  arraySOM3[i] + arraySOM4[i])
+plt.plot(data.get_array('x'),new_array,label='Total C')
+for i in range(len(arrayC)):
+  new_array[i] = arrayLabileC[i] + \
+                 arrayCelluloseC[i] + arrayLigninC[i] + \
+                 (arraySOM1[i] + arraySOM2[i] + \
+                  arraySOM3[i] + arraySOM4[i])
+plt.plot(data.get_array('x'),new_array,label='Immobile C')
+plt.plot(data.get_array('x'),arrayC,label='Mineral C')
+plt.legend(loc=1)
+plt.setp(plt.gca().get_legend().get_texts(),fontsize='small')
+plt.setp(plt.gca().get_legend().get_frame().set_fill(False))
+plt.setp(plt.gca().get_legend().draw_frame(False))
+
+# N
+plt.subplot(2,2,4)
+plt.xlabel('Time [y]')
+plt.ylabel('Concentration [mol/m^3]')
+for i in range(len(arrayN)):
+  new_array[i] = arrayN[i] + arrayLabileN[i] + \
+                 arrayCelluloseN[i] + arrayLigninN[i] + \
+                 (arraySOM1[i] + arraySOM2[i]) / CN_ratio_12 + \
+                 (arraySOM3[i] + arraySOM4[i]) / CN_ratio_10
+plt.plot(data.get_array('x'),new_array,label='Total N')
+for i in range(len(arrayN)):
+  new_array[i] = arrayLabileN[i] + \
+                 arrayCelluloseN[i] + arrayLigninN[i] + \
+                 (arraySOM1[i] + arraySOM2[i]) / CN_ratio_12 + \
+                 (arraySOM3[i] + arraySOM4[i]) / CN_ratio_10
+plt.plot(data.get_array('x'),new_array,label='Immobile N')
+#plt.plot(data.get_array('x'),arrayLabileN,label='Labile N')
+#plt.plot(data.get_array('x'),arraySOM1*CN_fraction_12_N,label='SOM1 N')
+plt.plot(data.get_array('x'),arrayN,label='Mineral N')
+plt.legend(loc=1)
+plt.setp(plt.gca().get_legend().get_texts(),fontsize='small')
+plt.setp(plt.gca().get_legend().get_frame().set_fill(False))
+plt.setp(plt.gca().get_legend().draw_frame(False))
+
+f.subplots_adjust(hspace=0.2,wspace=0.2,
+                  bottom=.12,top=.9,
+                  left=.12,right=.9)
+
+plt.show()

shortcourse/exercises/CLM-CN/CLM-CN_database.dat

+'temperature points' 8 0. 25. 60. 100. 150. 200. 250. 300.
+'A(aq)' 0. 0. 0.
+'B(aq)' 0. 0. 0.
+'C(aq)' 0. 0. 0.
+'D(aq)' 0. 0. 0.
+'E(aq)' 0. 0. 0.
+'null' 0 0 0
+'null' 1 0. '0' 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
+'null' 0. 1 1. '0' 0. 0. 0. 0. 0. 0. 0. 0. 0.
+'null' 0. 1 0. '0' 0. 0. 0. 0. 0. 0. 0. 0. 0.
+'null' 1 0. '0' 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.

shortcourse/exercises/CLM-CN/doc/CLM-CN.tex

+\documentclass{beamer}
+
+\usepackage{comment}
+\usepackage{color}
+\usepackage{listings}
+\usepackage{verbatim}
+\usepackage{multicol}
+\usepackage{booktabs}
+\definecolor{green}{RGB}{0,128,0}
+
+\def\EQ#1\EN{\begin{equation*}#1\end{equation*}}
+\def\BA#1\EA{\begin{align*}#1\end{align*}}
+\def\BS#1\ES{\begin{split*}#1\end{split*}}
+\newcommand{\bc}{\begin{center}}
+\newcommand{\ec}{\end{center}}
+\newcommand{\eq}{\ =\ }
+\newcommand{\degc}{$^\circ$C}
+
+\def\p{\partial}
+\def\qbs{\boldsymbol{q}}
+\def\Dbs{\boldsymbol{D}}
+\def\A{\mathcal A}
+\def\gC{\mathcal C}
+\def\gD{\mathcal D}
+\def\gL{\mathcal L}
+\def\M{\mathcal M}
+\def\P{\mathcal P}
+\def\Q{\mathcal Q}
+\def\gR{\mathcal R}
+\def\gS{\mathcal S}
+\def\X{\mathcal X}
+\def\bnabla{\boldsymbol{\nabla}}
+\def\bnu{\boldsymbol{\nu}}
+\renewcommand{\a}{{\alpha}}
+%\renewcommand{\a}{{}}
+\newcommand{\s}{{\sigma}}
+\newcommand{\bq}{\boldsymbol{q}}
+\newcommand{\bz}{\boldsymbol{z}}
+\def\bPsi{\boldsymbol{\Psi}}
+
+\def\Li{\textit{L}}
+\def\Fb{\textbf{f}}
+\def\Jb{\textbf{J}}
+\def\cb{\textbf{c}}
+
+\def\Dt{\Delta t}
+\def\tpdt{{t + \Delta t}}
+\def\bpsi{\boldsymbol{\psi}}
+\def\dbpsi{\delta \boldsymbol{\psi}}
+\def\bc{\textbf{c}}
+\def\bx{\textbf{x}}
+\def\dbc{\delta \textbf{c}}
+\def\dbx{\delta \textbf{x}}
+\def\arrows{\rightleftharpoons}
+
+\newcommand{\bGamma}{\boldsymbol{\Gamma}}
+\newcommand{\bOmega}{\boldsymbol{\Omega}}
+%\newcommand{\bPsi}{\boldsymbol{\Psi}}
+%\newcommand{\bpsi}{\boldsymbol{\psi}}
+\newcommand{\bO}{\boldsymbol{O}}
+%\newcommand{\bnu}{\boldsymbol{\nu}}
+\newcommand{\bdS}{\boldsymbol{dS}}
+\newcommand{\bg}{\boldsymbol{g}}
+\newcommand{\bk}{\boldsymbol{k}}
+%\newcommand{\bq}{\boldsymbol{q}}
+\newcommand{\br}{\boldsymbol{r}}
+\newcommand{\bR}{\boldsymbol{R}}
+\newcommand{\bS}{\boldsymbol{S}}
+\newcommand{\bu}{\boldsymbol{u}}
+\newcommand{\bv}{\boldsymbol{v}}
+%\newcommand{\bz}{\boldsymbol{z}}
+\newcommand{\pressure}{P}
+
+\newcommand\gehcomment[1]{{{\color{orange} #1}}}
+\newcommand\add[1]{{{\color{blue} #1}}}
+\newcommand\remove[1]{\sout{{\color{red} #1}}}
+\newcommand\codecomment[1]{{{\color{green} #1}}}
+\newcommand\redcomment[1]{{{\color{red} #1}}}
+\newcommand\bluecomment[1]{{{\color{blue} #1}}}
+\newcommand\greencomment[1]{{{\color{green} #1}}}
+\newcommand\magentacomment[1]{{{\color{magenta} #1}}}
+
+\begin{comment}
+\tiny
+\scriptsize
+\footnotesize
+\small
+\normalsize
+\large
+\Large
+\LARGE
+\huge
+\Huge
+\end{comment}
+
+\begin{document}
+\title{CLM-CN Reaction\ldots in a Nutshell}
+\author{Glenn Hammond}
+\date{\today}
+
+%\frame{\titlepage}
+
+%-----------------------------------------------------------------------------
+\section{Description of CLM-CN Reaction}
+
+\subsection{Batch CLM-CN Reaction Conceptual Model}
+
+\frame{\frametitle{Schematic of CLM-CN Reaction Network}
+\includegraphics[width=\linewidth]{./CLM-CN_cycle}
+}
+
+%-----------------------------------------------------------------------------
+\subsection{Governing Equations}
+
+\frame{\frametitle{Reaction Expression}
+
+\Large
+
+\EQ\label{CN_rxn}
+\text{CN}_u \eq \left(1-f\right) \text{CN}_d + f \text{CO}_2 + n \text{N}_\text{mineral}
+\EN
+
+\EQ\label{n_calc}
+n \eq u - \left(1-f\right) d
+\EN
+
+\footnotesize
+\BA
+\text{CN}_u &\eq \text{upstream carbon pool } [\text{mol C}/m^3] \\
+\text{CN}_d &\eq \text{downstream carbon pool } [\text{mol C}/m^3] \\
+\text{CO}_2 &\eq \text{nitrogen concentration } [\text{mol CO}_2/m^3] \\
+\text{N}_\text{mineral} &\eq \text{nitrogen concentration } [\text{mol N}/m^3] \\
+u &\eq \text{\text{C}/\text{N} atomic weight ratio divided by upstream \text{C}/\text{N} ratio} \\
+d &\eq \text{\text{C}/\text{N} atomic weight ratio divided by downstream \text{C}/\text{N} ratio} \\
+f &\eq \text{respiration fraction} \\
+\EA
+}
+
+\frame{\frametitle{Kinetic Rate Expression}
+
+\Large
+
+\EQ\label{CN_kinetic_rxn}
+rate \eq f_T f_\theta f_\text{pi} k \text{CN}_u
+\EN
+
+%\bigskip
+%\normalsize
+%\footnotesize
+\scriptsize
+\BA
+f_T &\eq \exp\left[308.56 \left(\frac{1}{71.02}-\frac{1}{T - 227.13}\right)\right]\\
+f_\theta &\eq \frac{\log\left(\theta_\text{min}/\theta\right)}{\log\left(\theta_\text{min}/\theta_\text{max}\right)}\\
+f_\text{pi} &\eq \frac{\text{N}_\text{mineral}}{\text{N}_\text{mineral} + K_{\text{N}_\text{mineral}}} \\
+k &\eq \text{kinetic rate constant} [s^{-1}]\\
+T &\eq \text{temperature } [K] \\
+\theta &\eq \text{moisture content or saturation } [-] \\
+\text{CN}_u &\eq \text{upstream carbon pool } [\text{mol C}/m^3] \\
+\text{N}_\text{mineral} &\eq \text{nitrogen concentration } [\text{mol N}/m^3] \\
+K_\text{N} &\eq \text{nitrogen half saturation constant } [\text{mol N}/m^3]\\
+\EA
+
+}
+
+\frame{\frametitle{Kinetic Reaction Equations}
+
+
+\EQ
+\text{CN}_u \eq \left(1-f\right) \text{CN}_d + f \text{CO}_2 + n \text{N}_\text{mineral}
+\EN
+
+\EQ
+rate \eq f_T f_\theta f_\text{pi} k \text{CN}_u
+\EN
+
+\Large
+\BA
+\frac{\p}{\p t} \left(\text{CN}_u\right) &\eq -rate \\
+\frac{\p}{\p t} \left(\text{CN}_d\right) &\eq \left(1-f\right) rate \\
+\frac{\p}{\p t} \left(\text{CO}_2\right) &\eq f \cdot rate \\
+\frac{\p}{\p t} \left(\text{N}_\text{mineral}\right) &\eq n \cdot rate \\
+\EA
+
+}
+
+
+\frame{\frametitle{Newton-Raphson Method}
+\LARGE
+\BA
+\Fb\left(\bx^{k+1,i}\right) &\eq \frac{\bx^{k+1,i}-\bx^k}{\Dt} - \sum_{irxn} \bnu_{irxn} \cdot rate_{irxn} \\
+\\
+\Jb \dbx &\eq -\Fb\left(\bx^{k+1,i}\right) \\
+\\
+\bx^{k+1,i+1} &\eq \bx^{k+1,i} + \dbx
+\EA
+\small
+\BA
+k &\eq \text{time step} \\
+i &\eq \text{iteration}
+\EA
+}
+
+%-----------------------------------------------------------------------------
+\subsection{PFLOTRAN Input Specification}
+
+\begin{frame}[fragile,allowframebreaks]\frametitle{CHEMISTRY}
+\small
+\begin{semiverbatim}
+CHEMISTRY
+  ...
+  IMMOBILE_SPECIES
+    N
+    C
+    SOM1
+    SOM2
+    SOM3
+    SOM4
+    LabileC
+    CelluloseC
+    LigninC
+    LabileN
+    CelluloseN
+    LigninN
+  /
+  ...
+\end{semiverbatim}
+\newpage
+\begin{semiverbatim}
+
+  ...
+  REACTION_SANDBOX
+    CLM-CN
+      POOLS   \bluecomment{! CN ratio}
+        SOM1   12.d0   \bluecomment{! -> u or d}
+        SOM2   12.d0
+        SOM3   10.d0
+        SOM4   10.d0
+        Labile
+        Cellulose
+        Lignin
+      /
+      ...
+\end{semiverbatim}
+\newpage
+\begin{semiverbatim}
+      ...
+      REACTION
+        UPSTREAM_POOL Labile         \bluecomment{! CN_u}
+        DOWNSTREAM_POOL SOM1         \bluecomment{! CN_d}
+        TURNOVER_TIME 20. h          \bluecomment{! -> k}
+        RESPIRATION_FRACTION 0.39d0  \bluecomment{! f}
+        N_INHIBITION 1.d-10          \bluecomment{! K_N}
+      /
+      REACTION
+        UPSTREAM_POOL SOM1
+        DOWNSTREAM_POOL SOM2
+        TURNOVER_TIME 14. d
+        RESPIRATION_FRACTION 0.28d0
+        N_INHIBITION 1.d-10
+      /
+      ...
+    / \bluecomment{! END CLM-CN}
+  / \bluecomment{! END REACTION_SANDBOX}
+/ \bluecomment{! END CHEMISTRY}
+\end{semiverbatim}
+\end{frame}
+
+\begin{frame}[fragile]\frametitle{CONSTRAINT}
+%\small
+\footnotesize
+\begin{semiverbatim}
+CONSTRAINT initial
+  CONCENTRATIONS   \bluecomment{! moles/L}
+    A(aq)  1.d-40  T
+  /
+  IMMOBILE       \bluecomment{! moles/m^3}
+    N     1.d-6
+    C     1.d-6
+    SOM1  1.d-10 \bluecomment{! moles C/m^3}
+    SOM2  1.d-10
+    SOM3  1.d-10
+    SOM4  1.d-10
+    LabileC     0.1852d-3
+    CelluloseC  0.4578d-3
+    LigninC     0.2662d-3
+    LabileN     0.00508954d-3
+    CelluloseN  0.01258096d-3
+    LigninN     0.00731553d-3
+  /
+END
+\end{semiverbatim}
+\end{frame}
+
+%-----------------------------------------------------------------------------
+\subsection{PFLOTRAN Numerical Implementation}
+\frame{\frametitle{Kinetic Reaction Equations}
+\BA
+\text{Lit1C} + \frac{1}{\text{CN}_{\text{Lit1C}}} {\text {Lit1N}} &\eq \left(1-f_1\right) \text{SOM1} + f_1 \text{CO}_2 + n_1 \text{N}_\text{mineral} \\
+\text{Lit2C} + \frac{1}{\text{CN}_{\text{Lit2C}}} {\text {Lit2N}} &\eq \left(1-f_2\right) \text{SOM2} + f_2 \text{CO}_2 + n_2 \text{N}_\text{mineral} \\
+\text{Lit3C} + \frac{1}{\text{CN}_{\text{Lit3C}}} {\text {Lit3N}} &\eq \left(1-f_3\right) \text{SOM3} + f_3 \text{CO}_2 + n_3 \text{N}_\text{mineral} \\
+\text{SOM1} &\eq \left(1-f_4\right) \text{SOM2} + f_4 \text{CO}_2 + n_4 \text{N}_\text{mineral} \\
+\text{SOM2} &\eq \left(1-f_5\right) \text{SOM3} + f_5 \text{CO}_2 + n_5 \text{N}_\text{mineral} \\
+\text{SOM3} &\eq \left(1-f_6\right) \text{SOM4} + f_6 \text{CO}_2 + n_6 \text{N}_\text{mineral} \\
+\text{SOM4} &\eq f_7 \text{CO}_2 + n_7 \text{N}_\text{mineral} \\
+\EA
+}
+
+\frame{\frametitle{Kinetic Reaction Equations}
+
+\BA
+\text{CN}_{\text{Lit1C}} & \eq  \frac{\text{Lit1C}}{\text{Lit1N}} \\
+\text{CN}_{\text{Lit2C}} & \eq \frac{\text{Lit2C}}{\text{Lit2N}} \\
+\text{CN}_{\text{Lit3C}} & \eq \frac{\text{Lit3C}}{\text{Lit3N}} \\
+\EA
+
+{\begin{center}
+\begin{tabular}{ c || c }
+  $\text{CN}_{\text{SOM1}} \eq 12$ & $\text{CN}_{\text{SOM2}} \eq 12$  \\
+  \hline
+  $\text{CN}_{\text{SOM3}} \eq 10$ & $\text{CN}_{\text{SOM4}} \eq 10$  \\
+  \hline
+  $f_1 \eq 0.39 $ & $f_2 \eq 0.55 $ \\
+  \hline
+  $f_3 \eq 0.29 $ & $f_4 \eq 0.28 $ \\
+  \hline
+  $f_5 \eq 0.46 $ & $f_5 \eq 0.55 $ \\
+  \hline
+  $f_7 \eq 1.0 $ & \\
+\end{tabular}
+\end{center}
+}
+}
+
+\frame{\frametitle{Kinetic Reaction Equations}
+\BA
+n_1 &\eq \frac{1}{\text{CN}_{\text{Lit1C}}} - \left(1-f_1\right) \frac{1}{\text{CN}_{\text{SOM1}}} \\
+n_2 &\eq \frac{1}{\text{CN}_{\text{Lit2C}}} - \left(1-f_2\right) \frac{1}{\text{CN}_{\text{SOM2}}} \\
+n_3 &\eq \frac{1}{\text{CN}_{\text{Lit3C}}} - \left(1-f_3\right) \frac{1}{\text{CN}_{\text{SOM3}}} \\
+n_4 &\eq \frac{1}{\text{CN}_{\text{SOM1}}} - \left(1-f_4\right) \frac{1}{\text{CN}_{\text{SOM2}}} \\
+n_5 &\eq \frac{1}{\text{CN}_{\text{SOM2}}} - \left(1-f_5\right) \frac{1}{\text{CN}_{\text{SOM3}}} \\
+n_6 &\eq \frac{1}{\text{CN}_{\text{SOM3}}} - \left(1-f_6\right) \frac{1}{\text{CN}_{\text{SOM4}}} \\
+n_7 &\eq \frac{1}{\text{CN}_{\text{SOM4}}} 
+\EA
+}
+
+\frame{\frametitle{Kinetic Reaction Equations}
+\BA
+R_1 \eq f_T f_\theta f_\text{pi} k \text{Lit1C} \\
+R_2 \eq f_T f_\theta f_\text{pi} k \text{Lit2C} \\
+R_3 \eq f_T f_\theta f_\text{pi} k \text{Lit3C} \\
+R_4 \eq f_T f_\theta f_\text{pi} k \text{SOM1} \\
+R_5 \eq f_T f_\theta f_\text{pi} k \text{SOM2} \\
+R_6 \eq f_T f_\theta f_\text{pi} k \text{SOM3} \\
+R_7 \eq f_T f_\theta f_\text{pi} k \text{SOM4} 
+\EA
+}
+
+
+\frame{\frametitle{Mass Conservation Equations}
+The CLM-CN reaction network results in 12 basis species for which the mass conservation equation are given by: 
+\BA
+\frac{\p}{\p t} \left(\text{N}_\text{mineral}\right) &\eq n_1R_1 + n_2R_2 + n_3R_3 +n_4R_4 +n_5R_5 +n_6R_6 + n_7R_7 \\
+\frac{\p}{\p t} \left(\text{CO}_2\right) &\eq  f_1R_1 + f_2R_2 + f_3R_3 + f_4R_4 + f_5R_5 + f_6R_6 + f_7R_7 \\
+\frac{\p}{\p t} \left(\text{SOM1}\right) &\eq  (1-f_1)R_1 - R_4\\
+\frac{\p}{\p t} \left(\text{SOM2}\right) &\eq  (1-f_2)R_2 + (1-f_4)R_4 - R_5\\
+\frac{\p}{\p t} \left(\text{SOM3}\right) &\eq  (1-f_3)R_3 + (1-f_5)R_5 - R_6\\
+\frac{\p}{\p t} \left(\text{SOM4}\right) &\eq  (1-f_6)R_6 - R_7
+\EA
+}
+
+\frame{\frametitle{Mass Conservation Equations}
+\BA
+\frac{\p}{\p t} \left(\text{Lit1C}\right) &\eq  -R_1 \\
+\frac{\p}{\p t} \left(\text{Lit2C}\right) &\eq  -R_2 \\
+\frac{\p}{\p t} \left(\text{Lit3C}\right) &\eq  -R_3 \\
+\frac{\p}{\p t} \left(\text{Lit1N}\right) &\eq   \frac{-1}{\text{CN}_{\text{Lit1C}}}R_1\\
+\frac{\p}{\p t} \left(\text{Lit3N}\right) &\eq  \frac{-1}{\text{CN}_{\text{Lit2C}}}R_2 \\
+\frac{\p}{\p t} \left(\text{Lit3N}\right) &\eq  \frac{-1}{\text{CN}_{\text{Lit3C}}}R_3 
+\EA
+}
+
+\frame{\frametitle{Numerical Implementation in PFLOTRAN}
+Applying finite-volume spatial discretization on mass conservation equation for Lit1C:
+\EQ
+\int \frac{\p}{\p t} \left(\text{Lit1C}\right) dV \eq  \int -R_1 dV
+\EN
+
+\EQ
+\frac{\p}{\p t} \left(\text{Lit1C }\Delta V\right)  \eq   -R_1 \Delta V
+\EN
+
+Implicit time discretization:
+\EQ
+\frac{\Delta V}{\Delta t} \left(\text{Lit1C}^{t+1} - \text{Lit1C}^t\right)  \eq   -R_1^{t+1} \Delta V
+\EN
+
+Residual equation:
+\EQ
+\mathcal{R}_7 \eq  \frac{\Delta V}{\Delta t} \left(\text{Lit1C}^{t+1} - \text{Lit1C}^t\right)  + R_1^{t+1} \Delta V
+\EN
+Note: The residual equation for Lit1C is $\mathcal{R}_7$ since Lit1C is the 7-th basis specie
+}
+
+\frame{\frametitle{Numerical Implementation in PFLOTRAN}
+Jacobian for Lit1C:
+\BA
+\mathcal{J}_{7,7} & \eq \frac{\p \mathcal{R}_7 }{\p \left(\text{Lit1C}^{t+1} \right)} \\
+\mathcal{J}_{7,1} & \eq \frac{\p \mathcal{R}_7 }{\p \left(\text{N}_\text{mineral}^{t+1} \right)} \\
+\mathcal{J}_{7,i} & \eq 0 \text{\hspace{2cm}i $\neq 1, 7$ }
+\EA
+}
+
+\end{document}

shortcourse/exercises/CLM-CN/doc/CLM-CN_cycle.png

Added
New image

src/pflotran/hdf5.F90

   call h5pclose_f(prop_id,hdf5_err)
   
   ! Open dataset
-  string = 'Region/'//trim(region%name)
+  string = 'Regions/'//trim(region%name)
   call h5dopen_f(file_id,string,data_set_id,hdf5_err)
 
   ! Get dataset's dataspace

src/pflotran/makefile

 	${common_src}reaction.o \
 	${common_src}reaction_aux.o \
 	${common_src}reaction_sandbox.o \
+	${common_src}reaction_sandbox_base.o \
 	${common_src}reaction_sandbox_clm_cn.o \
-	${common_src}reaction_sandbox_base.o \
+	${common_src}reaction_sandbox_example.o \
 	${common_src}solid_solution.o \
 	${common_src}solid_solution_aux.o \
 	${common_src}surface_complexation.o \
                  microbial_aux.o immobile_aux.o
 reaction_sandbox.o : option.o reaction_aux.o reactive_transport_aux.o string.o \
                      input.o utility.o global_aux.o reaction_sandbox_base.o \
-                     reaction_sandbox_clm_cn.o
+                     reaction_sandbox_clm_cn.o reaction_sandbox_example.o
 reaction_sandbox_base.o : reaction_aux.o reactive_transport_aux.o \
                           global_aux.o
 reaction_sandbox_clm_cn.o : option.o reaction_aux.o reactive_transport_aux.o \
                             string.o input.o utility.o global_aux.o \
                             reaction_sandbox_base.o database_aux.o \
                             immobile_aux.o
+reaction_sandbox_example.o : reaction_aux.o option.o reactive_transport_aux.o \
+                             global_aux.o units.o string.o input.o
 #reaction_chunk.o : reactive_transport_aux.o reaction_aux.o option.o string.o
 reactive_transport.o : realization.o reaction.o transport.o matrix_block_aux.o \
                        reactive_transport_aux.o fluid.o logging.o global.o \

src/pflotran/reaction.F90

   PetscInt :: srfcplx_count
   PetscInt :: temp_srfcplx_count
   PetscBool :: found
+  PetscBool :: reaction_sandbox_read
 
   nullify(prev_species)
   nullify(prev_gas)
   nullify(prev_kd_rxn)
   nullify(prev_ionx_rxn)
   
+  reaction_sandbox_read = PETSC_FALSE
+  
   srfcplx_count = 0
   input%ierr = 0
   do
 
       case('REACTION_SANDBOX')
         call RSandboxRead(input,option)
+        reaction_sandbox_read = PETSC_TRUE
       case('MICROBIAL_REACTION')
         call MicrobialRead(reaction%microbial,input,option)
       case('MINERALS')
   endif
   if (reaction%neqcplx + reaction%nsorb + reaction%mineral%nmnrl + &
       reaction%ngeneral_rxn + reaction%microbial%nrxn + &
-      reaction%immobile%nimmobile > 0) then
+      reaction%immobile%nimmobile > 0 .or. &
+      reaction_sandbox_read) then
     reaction%use_full_geochemistry = PETSC_TRUE
   endif
       

src/pflotran/reaction_sandbox.F90

 
   use Reaction_Sandbox_Base_class
   use Reaction_Sandbox_CLM_CN_class
+  use Reaction_Sandbox_Example_class
+  
+  ! Add new reacton sandbox classes here.
   
   implicit none
   
     select case(trim(word))
       case('CLM-CN')
         new_sandbox => CLM_CN_Create()
+      ! Add new cases statements for new reacton sandbox classes here.
+      case('EXAMPLE')
+        new_sandbox => EXAMPLECreate()
       case default
         option%io_buffer = 'CHEMISTRY,REACTION_SANDBOX keyword: ' // &
           trim(word) // ' not recognized.'

src/pflotran/reaction_sandbox_base.F90

     class(reaction_sandbox_base_type), pointer :: next
   contains
 #if 0  
+    procedure(Base_Read), public, deferred :: ReadInput
     procedure(Base_Setup), public, deferred :: Setup 
-    procedure(Base_Read), public, deferred :: ReadInput
     procedure(Base_React), public, deferred :: Evaluate
     procedure(Base_Destroy), public, deferred :: Destroy
 #else

src/pflotran/reaction_sandbox_clm_cn.F90

     type(pool_type), pointer :: pools
     type(clm_cn_reaction_type), pointer :: reactions
   contains
+    procedure, public :: ReadInput => CLM_CN_Read
     procedure, public :: Setup => CLM_CN_Setup
-    procedure, public :: ReadInput => CLM_CN_Read
     procedure, public :: Evaluate => CLM_CN_React
     procedure, public :: Destroy => CLM_CN_Destroy
   end type reaction_sandbox_clm_cn_type
 
 ! ************************************************************************** !
 !
-! CLM_CN_Setup: Sets up CLM-CN reaction after it has been read from input
-! author: Glenn Hammond
-! date: 02/04/13
-!
-! ************************************************************************** !
-subroutine CLM_CN_Setup(this,reaction,option)
-
-  use Reaction_Aux_module, only : reaction_type
-  use Option_module
-  
-  implicit none
-  
-  class(reaction_sandbox_clm_cn_type) :: this
-  type(reaction_type) :: reaction  
-  type(option_type) :: option
-  
-  call CLM_CN_Map(this,reaction,option)
-
-end subroutine CLM_CN_Setup
-
-! ************************************************************************** !
-!
 ! CLM_CN_Read: Reads input deck for reaction sandbox parameters
 ! author: Glenn Hammond
 ! date: 02/04/13
 
 ! ************************************************************************** !
 !
+! CLM_CN_Setup: Sets up CLM-CN reaction after it has been read from input
+! author: Glenn Hammond
+! date: 02/04/13
+!
+! ************************************************************************** !
+subroutine CLM_CN_Setup(this,reaction,option)
+
+  use Reaction_Aux_module, only : reaction_type
+  use Option_module
+  
+  implicit none
+  
+  class(reaction_sandbox_clm_cn_type) :: this
+  type(reaction_type) :: reaction  
+  type(option_type) :: option
+  
+  call CLM_CN_Map(this,reaction,option)
+
+end subroutine CLM_CN_Setup
+
+! ************************************************************************** !
+!
 ! CLM_CN_Map: Maps coefficients to primary dependent variables
 ! author: Glenn Hammond
 ! date: 02/04/13
   PetscReal :: resp_frac
   PetscReal :: stoich_N
   PetscReal :: stoich_C
-  PetscReal :: stoich_downstream_pool
+  PetscReal :: stoich_downstreamC_pool
   PetscReal :: stoich_upstreamC_pool, stoich_upstreamN_pool
   
   PetscReal :: N_inhibition, d_N_inhibition
       ispec_pool_down = this%pool_id_to_species_id(SOM_INDEX,ipool_down)
       CN_ratio_down = this%CN_ratio(ipool_down)
       ! c = (1-resp_frac) * a
-      stoich_downstream_pool = (1.d0-resp_frac) * stoich_upstreamC_pool
+      stoich_downstreamC_pool = (1.d0-resp_frac) * stoich_upstreamC_pool
     else    
       ispec_pool_down = 0
-      stoich_downstream_pool = 0.d0
+      stoich_downstreamC_pool = 0.d0
       CN_ratio_down = 1.d0 ! to prevent divide by zero below.
     endif
       
     ! d = resp_frac * a
     stoich_C = resp_frac * stoich_upstreamC_pool
     ! e = b - c / CN_ratio_dn
-    stoich_N = stoich_upstreamN_pool - stoich_downstream_pool / CN_ratio_down
+    stoich_N = stoich_upstreamN_pool - stoich_downstreamC_pool / CN_ratio_down
  
     ! Inhibition by nitrogen (inhibition concentration > 0 and N is a reactant)
     ! must be calculated here as the sign on the stoichiometry for N is 
       ! downstream pool
       ires_pool_down = reaction%offset_immobile + ispec_pool_down
       Residual(ires_pool_down) = Residual(ires_pool_down) - &
-        stoich_downstream_pool * rate
-      sumC = sumC + stoich_downstream_pool * rate
-      sumN = sumN + stoich_downstream_pool / CN_ratio_down * rate
+        stoich_downstreamC_pool * rate
+      sumC = sumC + stoich_downstreamC_pool * rate
+      sumN = sumN + stoich_downstreamC_pool / CN_ratio_down * rate
     endif
     
     !for debugging
       if (ispec_pool_down > 0) then
         Jacobian(ires_pool_down,iresC_pool_up) = &
           Jacobian(ires_pool_down,iresC_pool_up) - &
-          stoich_downstream_pool * drate
+          stoich_downstreamC_pool * drate
         if (use_N_inhibition) then
           Jacobian(ires_pool_down,ires_N) = &
             Jacobian(ires_pool_down,ires_N) - &
-            stoich_downstream_pool * drate_dN_inhibition
+            stoich_downstreamC_pool * drate_dN_inhibition
         endif
       endif
       
         ! dstoichC_dN_pool_up = 0
 
         ! nitrogen (stoichiometry a function of upstream C/N)
-        ! stoich_N = stoich_upstreamN_pool - stoich_downstream_pool / &
+        ! stoich_N = stoich_upstreamN_pool - stoich_downstreamC_pool / &
         !                                    CN_ratio_down
         ! latter half is constant
         dstoichN_dC_pool_up = dstoich_upstreamN_pool_dC_pool_up

src/pflotran/reaction_sandbox_example.F90

+module Reaction_Sandbox_Example_class
+
+! 1. Change all references to "Example" as desired to rename the module and
+!    and subroutines within the module. 
+
+  use Reaction_Sandbox_Base_class
+  
+  use Global_Aux_module
+  use Reactive_Transport_Aux_module
+  
+  implicit none
+  
+  private
+  
+#include "definitions.h"
+
+! 2. Add module variables here.  Note that one must use the PETSc data types 
+!    PetscInt, PetscReal, PetscBool to declare variables of type integer
+!    float/real*8, and logical respectively.  E.g.
+!
+! PetscReal, parameter :: formula_weight_of_water = 18.01534d0
+  
+  type, public, &
+    extends(reaction_sandbox_base_type) :: reaction_sandbox_example_type
+! 3. Add variables/arrays associated with new reaction.  
+    character(len=MAXWORDLENGTH) :: species_name
+    PetscInt :: species_id
+    PetscReal :: rate_constant
+  contains
+    procedure, public :: ReadInput => ExampleRead
+    procedure, public :: Setup => ExampleSetup
+    procedure, public :: Evaluate => ExampleReact
+    procedure, public :: Destroy => ExampleDestroy
+  end type reaction_sandbox_example_type
+
+  public :: ExampleCreate
+
+contains
+
+! ************************************************************************** !
+!
+! ExampleCreate: Allocates example reaction object.
+! author: John Doe (replace in all subroutine headers with name of developer) 
+! date: 00/00/00 (replace in all subroutine headers with current date)
+!
+! ************************************************************************** !
+function ExampleCreate()
+
+  implicit none
+  
+  class(reaction_sandbox_example_type), pointer :: ExampleCreate
+
+! 4. Add code to allocate object and initialized all variables to zero and
+!    nullify all pointers. E.g.
+  allocate(ExampleCreate)
+  ExampleCreate%species_name = ''
+  ExampleCreate%species_id = 0
+  ExampleCreate%rate_constant = 0.d0
+  nullify(ExampleCreate%next)  
+      
+end function ExampleCreate
+
+! ************************************************************************** !
+!
+! ExampleRead: Reads input deck for example reaction parameters (if any)
+! author: John Doe
+! date: 00/00/00
+!
+! ************************************************************************** !
+subroutine ExampleRead(this,input,option)
+
+  use Option_module
+  use String_module
+  use Input_module
+  use Units_module, only : UnitsConvertToInternal
+  
+  implicit none
+  
+  class(reaction_sandbox_example_type) :: this
+  type(input_type) :: input
+  type(option_type) :: option
+
+  PetscInt :: i
+  character(len=MAXWORDLENGTH) :: word
+  
+  do 
+    call InputReadFlotranString(input,option)
+    if (InputError(input)) exit
+    if (InputCheckExit(input,option)) exit
+
+    call InputReadWord(input,option,word,PETSC_TRUE)
+    call InputErrorMsg(input,option,'keyword', &
+                       'CHEMISTRY,REACTION_SANDBOX,TEMPLATE')
+    call StringToUpper(word)   
+
+    select case(trim(word))
+
+      ! Example Input:
+
+      ! CHEMISTRY
+      !   ...
+      !   REACTION_SANDBOX
+      !   : begin user-defined input
+      !     TEMPLATE
+      !       EXAMPLE_INTEGER 1
+      !       EXAMPLE_INTEGER_ARRAY 2 3 4
+      !     END
+      !   : end user defined input
+      !   END
+      !   ...
+      ! END
+
+! 5. Add case statement for reading variables.
+      case('SPECIES_NAME')
+! 6. Read the variable
+        ! Read the character string indicating which of the primary species
+        ! is being decayed.
+        call InputReadWord(input,option,this%species_name,PETSC_TRUE)  
+! 7. Inform the user of any errors if not read correctly.
+        call InputErrorMsg(input,option,'species_name', &
+                           'CHEMISTRY,REACTION_SANDBOX,EXAMPLE')    
+! 8. Repeat for other variables
+      case('RATE_CONSTANT')
+        ! Read the double precision rate constant
+        call InputReadDouble(input,option,this%rate_constant)
+        call InputErrorMsg(input,option,'rate_constant', &
+                           'CHEMISTRY,REACTION_SANDBOX,EXAMPLE')
+        ! Read the units
+        call InputReadWord(input,option,word,PETSC_TRUE)
+        if (InputError(input)) then
+          ! If units do not exist, assume default units of 1/s which are the
+          ! standard internal PFLOTRAN units for this rate constant.
+          input%err_buf = 'REACTION_SANDBOX,EXAMPLE,RATE CONSTANT UNITS'
+          call InputDefaultMsg(input,option)
+        else              
+          ! If units exist, convert to internal units of 1/s
+          this%rate_constant = this%rate_constant * &
+            UnitsConvertToInternal(word,option)
+        endif
+      case default
+        option%io_buffer = 'CHEMISTRY,REACTION_SANDBOX,TEMPLATE keyword: ' // &
+          trim(word) // ' not recognized.'
+        call printErrMsg(option)
+    end select
+  enddo
+  
+end subroutine ExampleRead
+
+! ************************************************************************** !
+!
+! ExampleSetup: Sets up the example reaction either with parameters either
+!                read from the input deck or hardwired.
+! author: John Doe
+! date: 00/00/00
+!
+! ************************************************************************** !
+subroutine ExampleSetup(this,reaction,option)
+
+  use Reaction_Aux_module, only : reaction_type, GetPrimarySpeciesIDFromName
+  use Option_module
+
+  implicit none
+  
+  class(reaction_sandbox_example_type) :: this
+  type(reaction_type) :: reaction
+  type(option_type) :: option
+
+! 9. Add code to initialize 
+  this%species_id = &
+    GetPrimarySpeciesIDFromName(this%species_name,reaction,option)
+      
+end subroutine ExampleSetup
+
+! ************************************************************************** !
+!
+! ExampleReact: Evaluates reaction storing residual and/or Jacobian
+! author: John Doe
+! date: 00/00/00
+!
+! ************************************************************************** !
+subroutine ExampleReact(this,Residual,Jacobian,compute_derivative, &
+                         rt_auxvar,global_auxvar,porosity,volume,reaction, &
+                         option)
+
+  use Option_module
+  use Reaction_Aux_module
+  
+  implicit none
+  
+  class(reaction_sandbox_example_type) :: this  
+  type(option_type) :: option
+  type(reaction_type) :: reaction
+  PetscBool :: compute_derivative
+  ! the following arrays must be declared after reaction
+  PetscReal :: Residual(reaction%ncomp)
+  PetscReal :: Jacobian(reaction%ncomp,reaction%ncomp)
+  PetscReal :: porosity
+  PetscReal :: volume
+  type(reactive_transport_auxvar_type) :: rt_auxvar
+  type(global_auxvar_type) :: global_auxvar
+
+  PetscInt, parameter :: iphase = 1
+  PetscReal :: L_water
+  
+  ! Description of subroutine arguments:
+
+  ! Residual - 1D array storing residual entries in units mol/sec
+  ! Jacobian - 2D array storing Jacobian entires in units kg water/sec
+  !
+  !  Jacobian [kg water/sec] * dc [mol/kg water] = -Res [mol/sec]
+  !
+  ! compute_derivative - Flag indicating whether analtical derivative should
+  !   be calculated.  The user must provide either the analytical derivatives 
+  !   or a numerical approximation unless always running with 
+  !   NUMERICAL_JACOBIAN_RXN defined in input deck.  If the use of 
+  !   NUMERICAL_JACOBIAN_RXN is assumed, the user should provide an error 
+  !   message when compute_derivative is true.  E.g.
+  !
+  !   option%io_buffer = 'NUMERICAL_JACOBIAN_RXN must always be used ' // &
+  !                      'due to assumptions in Example'
+  !   call printErrMsg(option)
+  !
+  ! rt_auxvar - Object holding chemistry information (e.g. concentrations,
+  !   activity coefficients, mineral volume fractions, etc.).  See
+  !   reactive_transport_aux.F90.  
+  !
+  !   Useful variables:
+  !     rt_auxvar%total(:,iphase) - total component concentrations 
+  !                                 [mol/L water] for phase
+  !     rt_auxvar%pri_molal(:) - free ion concentrations [mol/kg water]
+  !     rt_auxvar%pri_act_coef(:) - activity coefficients for primary species
+  !     rt_auxvar%aqueous%dtotal(:,iphase) - derivative of total component
+  !                 concentration with respect to free ion [kg water/L water]
+  !
+  ! global_auxvar - Object holding information on flow (e.g. saturation,
+  !   density, viscosity, temperature, etc)
+  !
+  !   Useful variables:
+  !     global_auxvar%den(iphase) - fluid density [mol/m^3] 
+  !     global_auxvar%den_kg(iphase) - fluid density [kg/m^3] 
+  !     global_auxvar%sat(iphase) - saturation 
+  !     global_auxvar%temp - temperature [C]
+  !
+  ! porosity - effective porosity of grid cell [m^3 pore/m^3 bulk]                     
+  ! volume - volume of grid cell [m^3]
+  ! reaction - Provides access to variable describing chemistry.  E.g.
+  !   reaction%ncomp - # chemical degrees of freedom (mobile and immobile)
+  !   reaction%naqcomp - # chemical degrees of freedom on water
+  !   reaction%primary_species_names(:) - names of primary species
+  !
+  ! option - Provides handle for controlling simulation, catching and
+  !          reporting errors.
+  
+  ! Unit of the residual must be in moles/second
+  ! global_auxvar%sat(iphase) = saturation of cell
+  ! 1.d3 converts m^3 water -> L water
+  L_water = porosity*global_auxvar%sat(iphase)*volume*1.d3
+  ! alway subtract contribution from residual
+  Residual(this%species_id) = Residual(this%species_id) - &
+    this%rate_constant * &  ! 1/sec
+    L_water * & ! L water
+    ! rt_auxvar%total(this%species_id,iphase) = species total component 
+    !   concentration
+    rt_auxvar%total(this%species_id,iphase) ! mol/L water
+    
+  
+  
+  if (compute_derivative) then
+
+! 11. If using an analytical Jacobian, add code for Jacobian evaluation
+
+    ! always add contribution to Jacobian
+    ! units = (mol/sec)*(kg water/mol) = kg water/sec
+    Jacobian(this%species_id,this%species_id) = &
+    Jacobian(this%species_id,this%species_id) + &
+      this%rate_constant * & ! 1/sec
+      L_water * & ! L water
+                  ! kg water/L water
+      ! rt_auxvar%aqueous%dtotal(this%species_id,this%species_id,iphase) = 
+      !   derivative of total component concentration with respect to the
+      !   free ion concentration of the same species.
+      rt_auxvar%aqueous%dtotal(this%species_id,this%species_id,iphase) 
+
+  endif
+  
+end subroutine ExampleReact
+
+! ************************************************************************** !
+!
+! ExampleDestroy: Destroys allocatable or pointer objects created in this 
+!                  module
+! author: John Doe
+! date: 00/00/00
+!
+! ************************************************************************** !
+subroutine ExampleDestroy(this)
+
+  implicit none
+  
+  class(reaction_sandbox_example_type) :: this  
+
+! 12. Add code to deallocate contents of the example object
+
+end subroutine ExampleDestroy
+
+end module Reaction_Sandbox_Example_class

src/pflotran/reaction_sandbox_example_A.F90

     PetscInt :: species_id
     PetscReal :: rate_constant
   contains
+    procedure, public :: ReadInput => ExampleRead
     procedure, public :: Setup => ExampleSetup
-    procedure, public :: ReadInput => ExampleRead
     procedure, public :: Evaluate => ExampleReact
     procedure, public :: Destroy => ExampleDestroy
   end type reaction_sandbox_example_type
   ExampleCreate%species_name = ''
   ExampleCreate%species_id = 0
   ExampleCreate%rate_constant = 0.d0
-      
+  nullify(ExampleCreate%next)
+  
 end function ExampleCreate
 
 ! ************************************************************************** !
 !
-! ExampleSetup: Sets up the example reaction either with parameters either
-!                read from the input deck or hardwired.
-! author: John Doe
-! date: 00/00/00
-!
-! ************************************************************************** !
-subroutine ExampleSetup(this,reaction,option)
-
-  use Reaction_Aux_module, only : reaction_type, GetPrimarySpeciesIDFromName
-  use Option_module
-
-  implicit none
-  
-  class(reaction_sandbox_example_type) :: this
-  type(reaction_type) :: reaction
-  type(option_type) :: option
-
-! 5. Add code to initialize   
-  this%species_name = 'A(aq)'
-  this%species_id = &
-    GetPrimarySpeciesIDFromName(this%species_name,reaction,option)
-  this%rate_constant = 4.3959d-9 ! 5 yr half-life
-  
-end subroutine ExampleSetup
-
-! ************************************************************************** !
-!
 ! ExampleRead: Reads input deck for example reaction parameters (if any)
 ! author: John Doe
 ! date: 00/00/00
 !   ...
 ! END
  
-! 6. Add case statement for reading variables.  E.g.
+! 5. Add case statement for reading variables.  E.g.
 !     case('EXAMPLE_INTEGER')
-! 7. Read the variable
+! 6. Read the variable
 !       call InputReadInt(input,option,this%example_integer)  
-! 8. Inform the user of any errors if not read correctly.
+! 7. Inform the user of any errors if not read correctly.
 !       call InputErrorMsg(input,option,'example_integer', & 
 !                          'CHEMISTRY,REACTION_SANDBOX,EXAMPLE') 
-! 9. Repeat for other variables
+! 8. Repeat for other variables
 !     case('EXAMPLE_INTEGER_Array')
 !       allocate(this%example_integer_array(3))
 !       this%example_integer_array = 0
 
 ! ************************************************************************** !
 !
+! ExampleSetup: Sets up the example reaction either with parameters either
+!                read from the input deck or hardwired.
+! author: John Doe
+! date: 00/00/00
+!
+! ************************************************************************** !
+subroutine ExampleSetup(this,reaction,option)
+
+  use Reaction_Aux_module, only : reaction_type, GetPrimarySpeciesIDFromName
+  use Option_module
+
+  implicit none
+  
+  class(reaction_sandbox_example_type) :: this
+  type(reaction_type) :: reaction
+  type(option_type) :: option
+
+! 9. Add code to initialize   
+  this%species_name = 'A(aq)'
+  this%species_id = &
+    GetPrimarySpeciesIDFromName(this%species_name,reaction,option)
+  this%rate_constant = 4.3959d-9 ! 5 yr half-life
+  
+end subroutine ExampleSetup
+
+! ************************************************************************** !
+!
 ! ExampleReact: Evaluates reaction storing residual and/or Jacobian
 ! author: John Doe
 ! date: 00/00/00

src/pflotran/reaction_sandbox_example_B.F90

     PetscReal :: rate_constant                    ! Double precision rate 
                                                   !  constant
   contains
+    procedure, public :: ReadInput => ExampleRead
     procedure, public :: Setup => ExampleSetup
-    procedure, public :: ReadInput => ExampleRead
     procedure, public :: Evaluate => ExampleReact
     procedure, public :: Destroy => ExampleDestroy
   end type reaction_sandbox_example_type
   ExampleCreate%species_name = ''
   ExampleCreate%species_id = 0
   ExampleCreate%rate_constant = 0.d0
+  nullify(ExampleCreate%next)
       
 end function ExampleCreate
 
 ! ************************************************************************** !
 !
-! ExampleSetup: Sets up the example reaction either with parameters either
-!                read from the input deck or hardwired.
-! author: John Doe
-! date: 00/00/00
-!
-! ************************************************************************** !
-subroutine ExampleSetup(this,reaction,option)
-
-  use Reaction_Aux_module, only : reaction_type, GetPrimarySpeciesIDFromName
-  use Option_module
-
-  implicit none
-  
-  class(reaction_sandbox_example_type) :: this
-  type(reaction_type) :: reaction
-  type(option_type) :: option
-
-! 5. Add code to initialize   
-  this%species_name = 'A(aq)'
-  this%species_id = &
-    GetPrimarySpeciesIDFromName(this%species_name,reaction,option)
-  this%rate_constant = 4.3959d-9 ! 5 yr half-life
-  
-end subroutine ExampleSetup
-
-! ************************************************************************** !
-!
 ! ExampleRead: Reads input deck for example reaction parameters (if any)
 ! author: John Doe
 ! date: 00/00/00
 !   ...
 ! END
  
-! 6. Add case statement for reading variables.
+! 5. Add case statement for reading variables.
       case('SPECIES_NAME')
-! 7. Read the variable
+! 6. Read the variable
         ! Read the character string indicating which of the primary species
         ! is being decayed.
         call InputReadWord(input,option,this%species_name,PETSC_TRUE)  
-! 8. Inform the user of any errors if not read correctly.
+! 7. Inform the user of any errors if not read correctly.
         call InputErrorMsg(input,option,'species_name', &
                            'CHEMISTRY,REACTION_SANDBOX,EXAMPLE')    
-! 9. Repeat for other variables
+! 8. Repeat for other variables
       case('RATE_CONSTANT')
         ! Read the double precision rate constant
         call InputReadDouble(input,option,this%rate_constant)
 
 ! ************************************************************************** !
 !
+! ExampleSetup: Sets up the example reaction either with parameters either
+!                read from the input deck or hardwired.
+! author: John Doe
+! date: 00/00/00
+!
+! ************************************************************************** !
+subroutine ExampleSetup(this,reaction,option)
+
+  use Reaction_Aux_module, only : reaction_type, GetPrimarySpeciesIDFromName
+  use Option_module
+
+  implicit none
+  
+  class(reaction_sandbox_example_type) :: this
+  type(reaction_type) :: reaction
+  type(option_type) :: option
+
+! 9. Add code to initialize   
+  this%species_id = &
+    GetPrimarySpeciesIDFromName(this%species_name,reaction,option)
+  
+end subroutine ExampleSetup
+
+! ************************************************************************** !
+!
 ! ExampleReact: Evaluates reaction storing residual and/or Jacobian
 ! author: John Doe
 ! date: 00/00/00

src/pflotran/reaction_sandbox_template.F90

 !    Character strings must be sized to either MAXWORDLENGTH or 
 !    MAXSTRINGLENGTH where max string length is a very long string.
   contains
+    procedure, public :: ReadInput => TemplateRead
     procedure, public :: Setup => TemplateSetup
-    procedure, public :: ReadInput => TemplateRead
     procedure, public :: Evaluate => TemplateReact
     procedure, public :: Destroy => TemplateDestroy
   end type reaction_sandbox_template_type
 
 ! 4. Add code to allocate object and initialized all variables to zero and
 !    nullify all pointers. E.g.
-! allocate(TemplateCreate)
+  allocate(TemplateCreate)
 ! TemplateCreate%example_integer = 0
 ! nullify(TemplateCreate%example_integer_array)
-      
+  nullify(TemplateCreate%next)
+  
 end function TemplateCreate
 
 ! ************************************************************************** !
 !
-! TemplateSetup: Sets up the template reaction either with parameters either
-!                read from the input deck or hardwired.
-! author: John Doe
-! date: 00/00/00
-!
-! ************************************************************************** !
-subroutine TemplateSetup(this,reaction,option)
-
-  use Reaction_Aux_module, only : reaction_type
-  use Option_module
-
-  implicit none
-  
-  class(reaction_sandbox_template_type) :: this
-  type(reaction_type) :: reaction
-  type(option_type) :: option
-
-! 5. Add code to initialize 
-      
-end subroutine TemplateSetup
-
-! ************************************************************************** !
-!
 ! TemplateRead: Reads input deck for template reaction parameters (if any)
 ! author: John Doe
 ! date: 00/00/00
       !   ...
       ! END
 
-! 6. Add case statement for reading variables.  E.g.
+! 5. Add case statement for reading variables.  E.g.
 !     case('EXAMPLE_INTEGER')
-! 7. Read the variable
+! 6. Read the variable
 !       call InputReadInt(input,option,this%example_integer)  
-! 8. Inform the user of any errors if not read correctly.
+! 7. Inform the user of any errors if not read correctly.
 !       call InputErrorMsg(input,option,'example_integer', & 
 !                          'CHEMISTRY,REACTION_SANDBOX,TEMPLATE') 
-! 9. Repeat for other variables
+! 8. Repeat for other variables
 !     case('EXAMPLE_INTEGER_Array')
 !       allocate(this%example_integer_array(3))
 !       this%example_integer_array = 0
 
 ! ************************************************************************** !
 !
+! TemplateSetup: Sets up the template reaction either with parameters either
+!                read from the input deck or hardwired.
+! author: John Doe
+! date: 00/00/00
+!
+! ************************************************************************** !
+subroutine TemplateSetup(this,reaction,option)
+
+  use Reaction_Aux_module, only : reaction_type
+  use Option_module
+
+  implicit none
+  
+  class(reaction_sandbox_template_type) :: this
+  type(reaction_type) :: reaction
+  type(option_type) :: option
+
+! 9. Add code to initialize 
+      
+end subroutine TemplateSetup
+
+! ************************************************************************** !
+!
 ! TemplateReact: Evaluates reaction storing residual and/or Jacobian
 ! author: John Doe
 ! date: 00/00/00
   class(reaction_sandbox_template_type) :: this  
 
 ! 12. Add code to deallocate contents of the template object
-  deallocate(this%example_integer_array)
-  nullify(this%example_integer_array) 
+! deallocate(this%example_integer_array)
+! nullify(this%example_integer_array) 
 
 end subroutine TemplateDestroy
 

src/pflotran/region.F90

           case('TOP')
             region%iface = TOP_FACE
         end select
-    case('GRID','SURF_GRID')
+      case('GRID','SURF_GRID')
         call InputReadWord(input,option,word,PETSC_TRUE)
         call InputErrorMsg(input,option,'GRID','REGION')
         call StringToUpper(word)
           case default
             option%io_buffer = 'REGION keyword: GRID = '//trim(word)//'not supported yet'
           call printErrMsg(option)
-      end select
+        end select
       case default
         option%io_buffer = 'REGION keyword: '//trim(keyword)//' not recognized'
         call printErrMsg(option)