Repository: parflow/parflow (language: python)

Issue:
print_van_Genuchten
Stores the van Genuchten parameters `alpha`, `n`, `sres`, `ssat` in `problem_data`
and prints them in `solver_richards.c` if `print_subsurf_data` is set. Memory usage is increased by these four vectors which  store one value for every gridpoint. Performance should not be affected because the values are generated anyway.
Print van Genuchten data
Rework existing PR #292 so it does not store vectors that can be temporary. 
Compute the vectors to be stored a single time rather than each invocation of Saturation.



Discussion / hints:
> Stores the van Genuchten parameters `alpha`, `n`, `sres`, `ssat` in `problem_data`
> and prints them in `solver_richards.c` if `print_subsurf_data` is set. Memory usage is increased by these four vectors which store one value for every gridpoint. Performance should not be affected because the values are generated anyway.

this is a nice feature, we should get this PR through.  Is it ready to go?
Branch has conflicts that have to be resolved.

> Branch has conflicts that have to be resolved.

thanks, yes, that should have been obvious : )

@ben-bou this is a great feature add, it would be nice to have the conflicts resolved and get this pushed
Not sure what is going on with the changes, there are a lot of differences in file not related to this patch.
> Not sure what is going on with the changes, there are a lot of differences in file not related to this patch.

@smithsg84 I believe the file changes are there because I branched of @hokkanen 's CUDA11 branch. None of the conflicting files are actually relevant for this PR. How do we proceed with this?
The PR needs to have the non-related changes removed.
Why is the reading of the new vectors necessary?   It looks to me like the pd_alphas/N/sres/ssat arrays are being computed in the problem_saturation.c loop after the VGInit method is called so the read values will be overwritten.

This is going to impact performance more than suggested.   Yes the values are calculated but this is dramatically increasing the memory writes in the loop.   Going from 1 write to 5.    Memory movement is more important than FLOPS these days.


I have reworked this feature to avoid keeping the values in problem data and reassigning the vectors each time saturation is invoked.



Fix patch:
diff --git a/pfsimulator/parflow_lib/parflow_proto.h b/pfsimulator/parflow_lib/parflow_proto.h
index c7fc543d4..239175c61 100644
--- a/pfsimulator/parflow_lib/parflow_proto.h
+++ b/pfsimulator/parflow_lib/parflow_proto.h
@@ -1,4 +1,3 @@
-/* Header.c */
 typedef void (*AdvectionConcentrationInvoke) (ProblemData *problem_data, int phase, int concentration, Vector *old_concentration, Vector *new_concentration, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity, Vector *solid_mass_factor, double time, double deltat, int order);
 typedef PFModule *(*AdvectionConcentrationInitInstanceXtraType) (Problem *problem, Grid *grid, double *temp_data);
 
@@ -342,7 +341,6 @@ int KinsolPCSizeOfTempData(void);
 
 
 typedef void (*L2ErrorNormInvoke) (double time, Vector *pressure, ProblemData *problem_data, double *l2_error_norm);
-
 /* kokkos.cpp */
 void kokkosInit();
 void kokkosFinalize();
@@ -503,6 +501,7 @@ void PerturbSystem(Lattice *lattice, Problem *problem);
 
 /* pf_module.c */
 PFModule *NewPFModule(void *call, void *init_instance_xtra, void *free_instance_xtra, void *new_public_xtra, void *free_public_xtra, void *sizeof_temp_data, void *instance_xtra, void *public_xtra);
+PFModule *NewPFModuleExtended(void *call, void *init_instance_xtra, void *free_instance_xtra, void *new_public_xtra, void *free_public_xtra, void *sizeof_temp_data,  void *output, void *output_static,void *instance_xtra, void *public_xtra);
 PFModule *DupPFModule(PFModule *pf_module);
 void FreePFModule(PFModule *pf_module);
 
@@ -989,6 +988,7 @@ int RichardsBCInternalSizeOfTempData(void);
 
 typedef void (*SaturationInvoke) (Vector *phase_saturation, Vector *phase_pressure, Vector *phase_density, double gravity, ProblemData *problem_data, int fcn);
 typedef PFModule *(*SaturationInitInstanceXtraInvoke) (Grid *grid, double *temp_data);
+typedef PFModule *(*SaturationOutputStaticInvoke) (char *file_prefix, ProblemData *problem_data);
 
 /* problem_saturation.c */
 void Saturation(Vector *phase_saturation, Vector *phase_pressure, Vector *phase_density, double gravity, ProblemData *problem_data, int fcn);
@@ -997,6 +997,8 @@ void SaturationFreeInstanceXtra(void);
 PFModule *SaturationNewPublicXtra(void);
 void SaturationFreePublicXtra(void);
 int SaturationSizeOfTempData(void);
+void  SaturationOutput(void);
+void  SaturationOutputStatic(char *file_prefix, ProblemData *problem_data);
 
 typedef void (*SaturationConstitutiveInvoke) (Vector **phase_saturations);
 typedef PFModule *(*SaturationConstitutiveInitInstanceXtraInvoke) (Grid *grid);
diff --git a/pfsimulator/parflow_lib/pf_module.c b/pfsimulator/parflow_lib/pf_module.c
index 598643b57..66c2a1ced 100644
--- a/pfsimulator/parflow_lib/pf_module.c
+++ b/pfsimulator/parflow_lib/pf_module.c
@@ -65,6 +65,37 @@ PFModule  *NewPFModule(
   return new_module;
 }
 
+PFModule  *NewPFModuleExtended(
+                       void *call,
+                       void *init_instance_xtra,
+                       void *free_instance_xtra,
+                       void *new_public_xtra,
+                       void *free_public_xtra,
+                       void *sizeof_temp_data,
+		       void *output,
+		       void *output_static,
+                       void *instance_xtra,
+                       void *public_xtra)
+{
+  PFModule         *new_module;
+
+  new_module = talloc(PFModule, 1);
+
+  (new_module->call) = (void (*)())call;
+  (new_module->init_instance_xtra) = (void (*)())init_instance_xtra;
+  (new_module->free_instance_xtra) = (void (*)())free_instance_xtra;
+  (new_module->new_public_xtra) = (void (*)())new_public_xtra;
+  (new_module->free_public_xtra) = (void (*)())free_public_xtra;
+  (new_module->sizeof_temp_data) = (int (*)())sizeof_temp_data;
+  (new_module->output) = (void (*)())output;
+  (new_module->output_static) = (void (*)())output_static;
+
+  PFModuleInstanceXtra(new_module) = instance_xtra;
+  PFModulePublicXtra(new_module) = public_xtra;
+
+  return new_module;
+}
+
 
 /*--------------------------------------------------------------------------
  * DupPFModule
@@ -72,14 +103,16 @@ PFModule  *NewPFModule(
 
 PFModule  *DupPFModule(PFModule *pf_module)
 {
-  return NewPFModule((void*)(pf_module->call),
-                     (void*)(pf_module->init_instance_xtra),
-                     (void*)(pf_module->free_instance_xtra),
-                     (void*)(pf_module->new_public_xtra),
-                     (void*)(pf_module->free_public_xtra),
-                     (void*)(pf_module->sizeof_temp_data),
-                     PFModuleInstanceXtra(pf_module),
-                     PFModulePublicXtra(pf_module));
+  return NewPFModuleExtended((void*)(pf_module->call),
+			     (void*)(pf_module->init_instance_xtra),
+			     (void*)(pf_module->free_instance_xtra),
+			     (void*)(pf_module->new_public_xtra),
+			     (void*)(pf_module->free_public_xtra),
+			     (void*)(pf_module->sizeof_temp_data),
+			     (void*)(pf_module->output),
+			     (void*)(pf_module->output_static),
+			     PFModuleInstanceXtra(pf_module),
+			     PFModulePublicXtra(pf_module));
 }
 
 
diff --git a/pfsimulator/parflow_lib/pf_module.h b/pfsimulator/parflow_lib/pf_module.h
index bb11efa04..7f506c026 100644
--- a/pfsimulator/parflow_lib/pf_module.h
+++ b/pfsimulator/parflow_lib/pf_module.h
@@ -33,84 +33,291 @@
 #ifndef _PFMODULE_HEADER
 #define _PFMODULE_HEADER
 
-
 /*--------------------------------------------------------------------------
  * PFModule structure
  *--------------------------------------------------------------------------*/
 
+/**
+ * Object like structure to achieve polymorphism.
+ *
+ * ParFlow existed before C++ was in widespread use.  Polymorphism is
+ * achieved by using function pointers stored in a standard C
+ * structure. This pattern is used across ParFlow to enable the use of
+ * multiple algorithms without need for switch like statements at the
+ * invocation site.
+ *
+ * The Parflow Module pattern has a single main "call" method which
+ * invokes the algorithm.  It is a little like a functor in C++.  This
+ * is limited so two Output methods have been added in Extended
+ * modules.  And yes, moving to C++ would make this cleaner.
+ * 
+ * Construction is broken down into two calls.   PublicExtra is data that is 
+ * shared by all instances of that type of module.   The Extra concept is
+ * much like a static member of a C++ class.   InstanceExtra is the 
+ * data that is unique to each instance of an algorithm.
+ *
+ * The output methods can be called to output state that is 
+ * time variant (called at each timestep as indicated by the user)
+ * or time invariant (called only once).
+ */
 typedef struct {
+  /**
+   * Virtual method invoked to do the function of the algorithm.
+   */
   void (*call)();
+
+  /**
+   * Virtual method invoked to construct instance data.
+   * 
+   * Instance data is per instantiation the module type; it is 
+   * not shared between instances.
+   */
   void (*init_instance_xtra)();
+  
+  /**
+   * Virtual method invoked to free the instance data.
+   */
   void (*free_instance_xtra)();
+
+  /**
+   * Virtual method invoked to construct the class data.
+   *
+   * Public data is allocated once per module type.   It is
+   * similar to class data in C++.
+   */
   void (*new_public_xtra)();
+
+  /**
+   * Virtual method invoked to free the class data.
+   */
   void (*free_public_xtra)();
+  
+  /**
+   * Virtual method to compute the size of the temporary data.
+   *
+   * This returns the size of the temporary data used by 
+   * the instance.
+   */
   int (*sizeof_temp_data)();
 
+  /**
+   * Output time variant state.
+   *
+   * Output data from the module instance, once per time
+   * interval specified by the user.
+   */
+  void (*output)();
+
+  /**
+   * Output static (time invariant) state.
+   *
+   * Output data from the module instance, once per simulation.
+   * Time invariant data should be outputted in this method.
+   */
+  void (*output_static)();
+
+  /**
+   * The instance data for each instance of a module.
+   * This data is unique to each instance; it is not shared.
+   * Like object data members in C++.
+   */
   void  *instance_xtra;
+  /**
+   * The class data for each type of module.
+   * This data is shared by all instances.
+   * Like class data members in C++.
+   */
   void  *public_xtra;
 } PFModule;
 
-/*--------------------------------------------------------------------------
- * Global xtra pointer
- *   used by module macros below to pass module xtra data to module routines
- *--------------------------------------------------------------------------*/
-
+/**
+ * A 'this' pointer for module methods.
+ * 
+ * Used by module macros below to pass module xtra data to module
+ * routines This is similar in purpose to the 'this' pointer in C++.
+ * Rather than passing this in each function it is done through a
+ * thread pointer.  For non-threaded impls this is just a global
+ * value.
+ *
+ * This value is set before the by the PFModule macros before the
+ * function pointer is extracted from the module structure and
+ * invoked.
+ */
 #ifdef PARFLOW_GLOBALS
 amps_ThreadLocalDcl(PFModule  *, global_ptr_this_pf_module);
 #else
 amps_ThreadLocalDcl(extern PFModule *, global_ptr_this_pf_module);
 #endif
 
-/*--------------------------------------------------------------------------
- * Accessor macros
- *--------------------------------------------------------------------------*/
+/**
+ * Device copy of the 'this' module pointer.
+ */ 
+#if (PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA) && defined(__CUDACC__)
+#ifdef PARFLOW_GLOBALS
+__device__ PFModule *dev_global_ptr_this_pf_module;
+#else
+/* This extern requires CUDA separate compilation, otherwise nvcc compiler 
+   treats the pointer as static variable for each compilation unit          */
+extern __device__ PFModule *dev_global_ptr_this_pf_module;
+#endif // PARFLOW_GLOBALS
+#endif // PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA && __CUDACC__
 
+/**
+ * Return the module instance data.
+ *
+ * @param pf_module The module instance.
+ * @return Pointer to the instance data
+ */
 #define PFModuleInstanceXtra(pf_module)      (pf_module->instance_xtra)
+
+/**
+ * Return the module class (public) data.
+ *
+ * @param pf_module The module instance
+ * @return Pointer to the class data.
+ */
 #define PFModulePublicXtra(pf_module)        (pf_module->public_xtra)
 
+/**
+ * Accessor macro for the global/thread variable holding the 'this' pointer
+ *
+ * These accessor macros depends on host/device context.
+ */
+#ifdef __CUDA_ARCH__
+#define dev_global_this_pf_module amps_ThreadLocal(dev_global_ptr_this_pf_module)
+#define ThisPFModule dev_global_this_pf_module
+#else
 #define global_this_pf_module amps_ThreadLocal(global_ptr_this_pf_module)
 #define ThisPFModule global_this_pf_module
+#endif
 
 /*--------------------------------------------------------------------------
  * PFModule interface macros
  *--------------------------------------------------------------------------*/
 
+/**
+ * Invoke the main algorithm method in the module.
+ *
+ * Modules are based on the idea the module has single method that
+ * does the work of the algorithm.  A bit like the functor pattern in
+ * C++.   
+ * 
+ * @param type Function type to invoke.
+ * @param pf_module The module instance
+ * @param args Arguments for the method to invoke.
+ * @return The invoked method return
+ *
+ * @note
+ * Rational: This was done so all PFModules would have the same
+ * calling convention.
+ */
 #define PFModuleInvoke(type, pf_module, args) \
   (                                           \
    ThisPFModule = pf_module,                  \
    (*(type (*)())(ThisPFModule->call))args    \
   )
 
+/**
+ * Invoke the main algorithm method in the module.
+ *
+ * Modules are based on the idea the module has single method that
+ * does the work of the algorithm.  A bit like the functor pattern in
+ * C++.   
+ *
+ * @note
+ * Rational: This was done so all PFModules would have the same
+ * calling convention.
+ *
+ * @param type Function type to invoke.
+ * @param pf_module The module instance
+ * @param args Arguments for the method to invoke.
+ * @return The invoked method return
+ */
 #define PFModuleInvokeType(type, pf_module, args) \
   (                                               \
    ThisPFModule = pf_module,                      \
    (*(type)(ThisPFModule->call))args              \
   )
 
+/**
+ * Create a new instance of a module.
+ *
+ * Note we have done some renaming to better match OO (C++) naming
+ * schemes.  This invokes the "init" method in the method which is a
+ * ctor for the method instance.
+ * 
+ * @param type Function type to invoke.
+ * @param pf_module The module instance
+ * @param args Arguments for the method to invoke.
+ * @return The new module instance pointer.
+ */
 #define PFModuleNewInstance(pf_module, args)                   \
   (                                                            \
    ThisPFModule = DupPFModule(pf_module),                      \
    (*(PFModule * (*)())(ThisPFModule->init_instance_xtra))args \
   )
 
+/**
+ * Create a new instance of a module.
+ *
+ * Note we have done some renaming to better match OO (C++) naming
+ * schemes.  This invokes the "init" method in the method which is a
+ * ctor for the method instance.
+ *
+ * @param type Function type to invoke.
+ * @param pf_module The module instance
+ * @param args Arguments for the method to invoke.
+ * @return The new module instance pointer.
+ */
 #define PFModuleNewInstanceType(type, pf_module, args) \
   (                                                    \
    ThisPFModule = DupPFModule(pf_module),              \
    (*(type)(ThisPFModule->init_instance_xtra))args     \
   )
 
+/**
+ * 'ReNew' the instance of a module.
+ * 
+ * 'ReNew'ing a module is done when the module has already been
+ * constructed but some state may have changed.
+ * 
+ * \TODO this should be better documented.
+ *
+ * @param pf_module The module instance
+ * @param args Arguments for the method to invoke.
+ * @return The new module instance pointer.
+ */
 #define PFModuleReNewInstance(pf_module, args)                 \
   (                                                            \
    ThisPFModule = pf_module,                                   \
    (*(PFModule * (*)())(ThisPFModule->init_instance_xtra))args \
   )
 
+/**
+ * 'ReNew' the instance of a module.
+ * 
+ * 'ReNew'ing a module is done when the module has already been
+ * constructed but some state may have changed.
+ * 
+ * TODO this should be better documented.
+ *
+ * @param pf_module The module instance
+ * @param args Arguments for the method to invoke.
+ * @return The new module instance pointer.
+ */
 #define PFModuleReNewInstanceType(type, pf_module, args) \
   (                                                      \
    ThisPFModule = pf_module,                             \
    (*(type)(ThisPFModule->init_instance_xtra))args       \
   )
 
+/**
+ * Free the module.
+ *
+ * Similiar to a C++ destructor.
+ *
+ * @param pf_module The module instance
+ */
 #define PFModuleFreeInstance(pf_module)                 \
   (                                                     \
    ThisPFModule = pf_module,                            \
@@ -118,6 +325,51 @@ amps_ThreadLocalDcl(extern PFModule *, global_ptr_this_pf_module);
    FreePFModule(pf_module)                              \
   )
 
+/**
+ * Output time variant state associated with this module instance.
+ * 
+ * This method is invoked at the output time intervals specified by
+ * the user input.  The method should write any vectors or other state
+ * out for the current timestep.
+ *
+ * @param type Type of method to invoke for outputing
+ * @param pf_module The module instance
+ * @param args Arguments for the output method
+ */
+#define PFModuleOutputType(type, pf_module, args) \
+  (                                           \
+   ThisPFModule = pf_module,                  \
+   (*(type (*)())(ThisPFModule->output))args  \
+  )
+
+
+/**
+ * Output static (time invariant) state associated with this module instance.
+ * 
+ * This method is invoked at the start of a run and should write any
+ * vectors or other state for the module that does not require a
+ * timestamp.
+ *
+ * @param type Type of method to invoke for outputing
+ * @param pf_module The module instance
+ * @param args Arguments for the output method
+ */
+#define PFModuleOutputStaticType(type, pf_module, args) \
+  (                                           \
+   ThisPFModule = pf_module,                  \
+   (*(type (*)())(ThisPFModule->output_static))args  \
+  )
+
+/**
+ * Create a class of module.
+ *
+ * This allocates the public/class data for the module
+ * class.
+ *
+ * @param name PFModule name
+ * @param args Arguments for the module class constructor
+ * @return The new module instance
+ */
 #define PFModuleNewModule(name, args)                          \
   (                                                            \
    ThisPFModule = NewPFModule((void*)name,                     \
@@ -130,6 +382,16 @@ amps_ThreadLocalDcl(extern PFModule *, global_ptr_this_pf_module);
    (*(PFModule * (*)())(ThisPFModule->new_public_xtra))args    \
   )
 
+/**
+ * Create a class of module.
+ *
+ * This allocates the public/class data for the module
+ * class.
+ *
+ * @param name PFModule name
+ * @param args Arguments for the module class constructor
+ * @return The new module instance
+ */
 #define PFModuleNewModuleType(type, name, args)                \
   (                                                            \
    ThisPFModule = NewPFModule((void*)name,                     \
@@ -142,6 +404,66 @@ amps_ThreadLocalDcl(extern PFModule *, global_ptr_this_pf_module);
    (*(type)(ThisPFModule->new_public_xtra))args                \
   )
 
+/**
+ * Create a class of extended module.
+ * 
+ * For use with modules that implement the Extended module API.
+ * Currently this is the Output methods.
+ *
+ * This allocates the public/class data for the module
+ * class.
+ *
+ * @param name PFModule name
+ * @param args Arguments for the module class constructor
+ * @return The new module instance
+ */
+#define PFModuleNewModuleExtended(name, args)                  \
+  (                                                            \
+   ThisPFModule = NewPFModuleExtended((void*)name,                     \
+                              (void*)name ## InitInstanceXtra, \
+                              (void*)name ## FreeInstanceXtra, \
+                              (void*)name ## NewPublicXtra,    \
+                              (void*)name ## FreePublicXtra,   \
+                              (void*)name ## SizeOfTempData,   \
+                              (void*)name ## Output,	       \
+                              (void*)name ## OutputStatic,     \
+			      NULL, NULL),		       \
+   (*(PFModule * (*)())(ThisPFModule->new_public_xtra))args    \
+  )
+
+/**
+ * Create a class of extended module.
+ * 
+ * For use with modules that implement the Extended module API.
+ * Currently this is the Output methods.
+ *
+ * This allocates the public/class data for the module
+ * class.
+ *
+ * @param name PFModule name
+ * @param args Arguments for the module class constructor
+ * @return The new module instance
+ */
+#define PFModuleNewModuleExtendedType(type, name, args)        \
+  (                                                            \
+   ThisPFModule = NewPFModuleExtended((void*)name,                     \
+                              (void*)name ## InitInstanceXtra, \
+                              (void*)name ## FreeInstanceXtra, \
+                              (void*)name ## NewPublicXtra,    \
+                              (void*)name ## FreePublicXtra,   \
+                              (void*)name ## SizeOfTempData,   \
+                              (void*)name ## Output,	       \
+                              (void*)name ## OutputStatic,     \
+                              NULL, NULL),                     \
+   (*(type)(ThisPFModule->new_public_xtra))args                \
+  )
+
+
+/**
+ * Invoke the destructor of the module.
+ * 
+ * @param pf_module The module instance
+ */
 #define PFModuleFreeModule(pf_module)                 \
   (                                                   \
    ThisPFModule = pf_module,                          \
@@ -149,6 +471,18 @@ amps_ThreadLocalDcl(extern PFModule *, global_ptr_this_pf_module);
    FreePFModule(pf_module)                            \
   )
 
+/**
+ * Return this size of the temporary data needed by this module instance.
+ * 
+ * @note
+ * Rational: PF originally operated on very low memory compute nodes
+ * (4MB-16MB) so had to carefully manage space for vectors.   This 
+ * method was part of that mechanism.
+ *
+ * \TODO Believe this is deprecated.
+ *
+ * @param pf_module The module instance
+ */
 #define PFModuleSizeOfTempData(pf_module)           \
   (                                                 \
    ThisPFModule = pf_module,                        \
diff --git a/pfsimulator/parflow_lib/problem.c b/pfsimulator/parflow_lib/problem.c
index 8998ac8e5..f6a9a4bda 100644
--- a/pfsimulator/parflow_lib/problem.c
+++ b/pfsimulator/parflow_lib/problem.c
@@ -253,7 +253,7 @@ Problem   *NewProblem(
   else  /* Richards case */
   {
     ProblemSaturation(problem) =
-      PFModuleNewModule(Saturation, ());
+      PFModuleNewModuleExtended(Saturation, ());
   }
 
   /*-----------------------------------------------------------------------
diff --git a/pfsimulator/parflow_lib/problem_saturation.c b/pfsimulator/parflow_lib/problem_saturation.c
index 4673589af..30231738f 100644
--- a/pfsimulator/parflow_lib/problem_saturation.c
+++ b/pfsimulator/parflow_lib/problem_saturation.c
@@ -28,6 +28,7 @@
 
 #include "parflow.h"
 
+#include <string.h>
 #include <float.h>
 
 /*--------------------------------------------------------------------------
@@ -115,7 +116,7 @@ void     Saturation(
                                        *                  value
                                        * fcn = CALCDER => calculate the function
                                        *                  derivative */
-{
+{ 
   PFModule      *this_module = ThisPFModule;
   PublicXtra    *public_xtra = (PublicXtra*)PFModulePublicXtra(this_module);
 
@@ -804,7 +805,7 @@ void  SaturationFreeInstanceXtra()
       Type5* dummy5 = (Type5*)(public_xtra->data);
       FreeVector(dummy5->satRF);
     }
-    
+
     tfree(instance_xtra);
   }
 }
@@ -1130,7 +1131,7 @@ void  SaturationFreePublicXtra()
 	  tfree(dummy1->s_ress);
 	  tfree(dummy1->s_difs);
 	}
-	
+
         tfree(dummy1);
 
         break;
@@ -1221,3 +1222,219 @@ int  SaturationSizeOfTempData()
 
   return sz;
 }
+
+void  SaturationOutput()
+{
+  //PFModule    *this_module = ThisPFModule;
+  //PublicXtra  *public_xtra = (PublicXtra*)PFModulePublicXtra(this_module);
+  printf("SaturationOutput does nothing\n");
+}
+
+
+void  SaturationOutputStatic(
+			     char *file_prefix,
+			     ProblemData *problem_data /* Contains geometry info. for the problem */
+			     )
+{
+  PFModule      *this_module = ThisPFModule;
+  PublicXtra    *public_xtra = (PublicXtra*)PFModulePublicXtra(this_module);
+
+  Type1         *dummy1;
+
+  Grid          *grid = VectorGrid(ProblemDataSpecificStorage(problem_data));
+
+  GrGeomSolid   *gr_solid;
+
+  Subvector     *n_values_sub;
+  Subvector     *alpha_values_sub;
+  Subvector     *s_res_values_sub;
+  Subvector     *s_sat_values_sub;
+
+  double        *n_values_dat, *alpha_values_dat;
+  double        *s_res_values_dat, *s_sat_values_dat;
+
+  SubgridArray  *subgrids = GridSubgrids(grid);
+
+  Subgrid       *subgrid;
+
+  int sg;
+
+  int ix, iy, iz, r;
+  int nx, ny, nz;
+
+  int i, j, k;
+
+  int            *region_indices, num_regions, ir;
+
+  Subvector      *pd_alpha_sub;         //BB
+  Subvector      *pd_n_sub;         //BB
+  Subvector      *pd_sres_sub;         //BB
+  Subvector      *pd_ssat_sub;         //BB
+  double *pd_alpha_dat, *pd_n_dat, *pd_sres_dat, *pd_ssat_dat;    //BB
+
+  Vector *pd_alpha = NewVectorType(grid, 1, 1, vector_cell_centered);
+  Vector *pd_n = NewVectorType(grid, 1, 1, vector_cell_centered);
+  Vector *pd_sres = NewVectorType(grid, 1, 1, vector_cell_centered);
+  Vector *pd_ssat = NewVectorType(grid, 1, 1, vector_cell_centered);
+
+  /* Initialize saturations */
+  InitVector(pd_alpha, 0.0); 
+  InitVector(pd_n, 0.0);
+  InitVector(pd_sres, 0.0);
+  InitVector(pd_ssat, 0.0);
+
+  switch ((public_xtra->type))
+  {
+    case 1: /* Van Genuchten saturation curve */
+    {
+      int data_from_file;
+      double *alphas, *ns, *s_ress;
+
+      Vector *n_values, *alpha_values, *s_res_values, *s_sat_values;
+
+      dummy1 = (Type1*)(public_xtra->data);
+
+      num_regions = (dummy1->num_regions);
+      region_indices = (dummy1->region_indices);
+      alphas = (dummy1->alphas);
+      ns = (dummy1->ns);
+      s_ress = (dummy1->s_ress);
+      data_from_file = (dummy1->data_from_file);
+
+      if (data_from_file == 0) /* Soil parameters given by region */
+      {
+        for (ir = 0; ir < num_regions; ir++)
+        {
+          gr_solid = ProblemDataGrSolid(problem_data, region_indices[ir]);
+
+          ForSubgridI(sg, subgrids)
+          {
+            subgrid = SubgridArraySubgrid(subgrids, sg);
+
+            pd_alpha_sub = VectorSubvector(pd_alpha, sg);   //BB
+            pd_n_sub = VectorSubvector(pd_n, sg);           //BB
+            pd_sres_sub = VectorSubvector(pd_sres, sg);     //BB
+            pd_ssat_sub = VectorSubvector(pd_ssat, sg);     //BB
+
+            ix = SubgridIX(subgrid);
+            iy = SubgridIY(subgrid);
+            iz = SubgridIZ(subgrid);
+
+            nx = SubgridNX(subgrid);
+            ny = SubgridNY(subgrid);
+            nz = SubgridNZ(subgrid);
+
+            r = SubgridRX(subgrid);
+
+            pd_alpha_dat = SubvectorData(pd_alpha_sub);   //BB
+            pd_n_dat = SubvectorData(pd_n_sub);           //BB
+            pd_sres_dat = SubvectorData(pd_sres_sub);     //BB
+            pd_ssat_dat = SubvectorData(pd_ssat_sub);     //BB
+
+	    GrGeomInLoop(i, j, k, gr_solid, r, ix, iy, iz, nx, ny, nz,
+	    {
+	      int ips = SubvectorEltIndex(pd_alpha_sub, i, j, k);
+	      
+	      double alpha = alphas[ir];
+	      double n = ns[ir];
+	      double s_res = s_ress[ir];
+	      
+	      pd_alpha_dat[ips] = alpha;  //BB
+	      pd_n_dat[ips] = n;  //BB
+	      pd_sres_dat[ips] = s_res;  //BB  // no ssat???
+	      
+	    });
+          }     /* End subgrid loop */
+        }       /* End loop over regions */
+      }         /* End if data not from file */
+      else
+      {
+        gr_solid = ProblemDataGrDomain(problem_data);
+        n_values = dummy1->n_values;
+        alpha_values = dummy1->alpha_values;
+        s_res_values = dummy1->s_res_values;
+        s_sat_values = dummy1->s_sat_values;
+
+        ForSubgridI(sg, subgrids)
+        {
+          subgrid = SubgridArraySubgrid(subgrids, sg);
+
+          n_values_sub = VectorSubvector(n_values, sg);
+          alpha_values_sub = VectorSubvector(alpha_values, sg);
+          s_res_values_sub = VectorSubvector(s_res_values, sg);
+          s_sat_values_sub = VectorSubvector(s_sat_values, sg);
+
+          pd_alpha_sub = VectorSubvector(pd_alpha, sg);   //BB
+          pd_n_sub = VectorSubvector(pd_n, sg);           //BB
+          pd_sres_sub = VectorSubvector(pd_sres, sg);     //BB
+          pd_ssat_sub = VectorSubvector(pd_ssat, sg);     //BB
+
+          ix = SubgridIX(subgrid);
+          iy = SubgridIY(subgrid);
+          iz = SubgridIZ(subgrid);
+
+          nx = SubgridNX(subgrid);
+          ny = SubgridNY(subgrid);
+          nz = SubgridNZ(subgrid);
+
+          r = SubgridRX(subgrid);
+
+          n_values_dat = SubvectorData(n_values_sub);
+          alpha_values_dat = SubvectorData(alpha_values_sub);
+          s_res_values_dat = SubvectorData(s_res_values_sub);
+          s_sat_values_dat = SubvectorData(s_sat_values_sub);
+
+          pd_alpha_dat = SubvectorData(pd_alpha_sub);   //BB
+          pd_n_dat = SubvectorData(pd_n_sub);           //BB
+          pd_sres_dat = SubvectorData(pd_sres_sub);     //BB
+          pd_ssat_dat = SubvectorData(pd_ssat_sub);     //BB
+
+	  GrGeomInLoop(i, j, k, gr_solid, r, ix, iy, iz, nx, ny, nz,
+          {
+	    int ips = SubvectorEltIndex(pd_alpha_sub, i, j, k);
+	    
+	    int n_index = SubvectorEltIndex(n_values_sub, i, j, k);
+	    int alpha_index = SubvectorEltIndex(alpha_values_sub, i, j, k);
+	    int s_res_index = SubvectorEltIndex(s_res_values_sub, i, j, k);
+	    int s_sat_index = SubvectorEltIndex(s_sat_values_sub, i, j, k);
+	    
+	    double alpha = alpha_values_dat[alpha_index];
+	    double n = n_values_dat[n_index];
+	    double s_res = s_res_values_dat[s_res_index];
+	    double s_sat = s_sat_values_dat[s_sat_index];
+	    
+	    pd_alpha_dat[ips] = alpha;  //BB
+	    pd_n_dat[ips] = n;  //BB
+	    pd_sres_dat[ips] = s_res;  //BB
+	    pd_ssat_dat[ips] = s_sat;  //BB
+	  });
+        }       /* End subgrid loop */
+      }         /* End if data_from_file */
+      break;
+    }        /* End case 1 */
+  }          /* End switch */
+
+  char file_postfix[2048];
+  
+  strcpy(file_postfix, "alpha");
+  WritePFBinary(file_prefix, file_postfix,
+		pd_alpha);
+  
+  strcpy(file_postfix, "n");
+  WritePFBinary(file_prefix, file_postfix,
+		pd_n);
+  
+  strcpy(file_postfix, "sres");
+  WritePFBinary(file_prefix, file_postfix,
+		pd_sres);
+  
+  strcpy(file_postfix, "ssat");
+  WritePFBinary(file_prefix, file_postfix,
+		pd_ssat);
+
+  FreeVector(pd_alpha);
+  FreeVector(pd_n);
+  FreeVector(pd_sres);
+  FreeVector(pd_ssat);
+}
+
diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c
index 66dbdcc14..e5be05971 100644
--- a/pfsimulator/parflow_lib/solver_richards.c
+++ b/pfsimulator/parflow_lib/solver_richards.c
@@ -491,6 +491,8 @@ SetupRichards(PFModule * this_module)
     WritePFBinary(file_prefix, file_postfix,
                   ProblemDataSpecificStorage(problem_data));
 
+    PFModuleOutputStaticType(SaturationOutputStaticInvoke, ProblemSaturation(problem), (file_prefix, problem_data));
+    
     // Now add metadata entries:
     static const char* permeability_filenames[] = {
       "perm_x", "perm_y", "perm_z"


Test patch:
diff --git a/test/correct_output/crater2D_vangtable_linear.out.alpha.pfb b/test/correct_output/crater2D_vangtable_linear.out.alpha.pfb
new file mode 100644
index 000000000..30d39c602
Binary files /dev/null and b/test/correct_output/crater2D_vangtable_linear.out.alpha.pfb differ
diff --git a/test/correct_output/crater2D_vangtable_linear.out.n.pfb b/test/correct_output/crater2D_vangtable_linear.out.n.pfb
new file mode 100644
index 000000000..a3b2de275
Binary files /dev/null and b/test/correct_output/crater2D_vangtable_linear.out.n.pfb differ
diff --git a/test/correct_output/crater2D_vangtable_linear.out.sres.pfb b/test/correct_output/crater2D_vangtable_linear.out.sres.pfb
new file mode 100644
index 000000000..e9a477bb4
Binary files /dev/null and b/test/correct_output/crater2D_vangtable_linear.out.sres.pfb differ
diff --git a/test/correct_output/crater2D_vangtable_linear.out.ssat.pfb b/test/correct_output/crater2D_vangtable_linear.out.ssat.pfb
new file mode 100644
index 000000000..d3a5da4d0
Binary files /dev/null and b/test/correct_output/crater2D_vangtable_linear.out.ssat.pfb differ
diff --git a/test/.gitignore b/test/tcl/.gitignore
similarity index 100%
rename from test/.gitignore
rename to test/tcl/.gitignore
diff --git a/test/tcl/crater2D_vangtable_linear.tcl b/test/tcl/crater2D_vangtable_linear.tcl
index 21f9fcc87..7d7bcff77 100644
--- a/test/tcl/crater2D_vangtable_linear.tcl
+++ b/test/tcl/crater2D_vangtable_linear.tcl
@@ -512,17 +512,10 @@ pfundist $runname
 source pftest.tcl
 set passed 1
 
-if ![pftestFile $runname.out.perm_x.pfb "Max difference in perm_x" $sig_digits] {
-    set passed 0
-}
-if ![pftestFile $runname.out.perm_y.pfb "Max difference in perm_y" $sig_digits] {
-    set passed 0
-}
-if ![pftestFile $runname.out.perm_z.pfb "Max difference in perm_z" $sig_digits] {
-    set passed 0
-}
-if ![pftestFile $runname.out.porosity.pfb "Max difference in porosity" $sig_digits] {
-    set passed 0
+foreach file { perm_x perm_y perm_z alpha n sres ssat } {
+    if ![pftestFile $runname.out.$file.pfb "Max difference in $file" $sig_digits] {
+	set passed 0
+    }
 }
 
 foreach i "00000 00001 00002" {