diff --git a/demo/transient/transient_snes_mpfaof90.F90 b/demo/transient/transient_snes_mpfaof90.F90
index 127dcea0..8ae7fc31 100644
--- a/demo/transient/transient_snes_mpfaof90.F90
+++ b/demo/transient/transient_snes_mpfaof90.F90
@@ -9,8 +9,8 @@ module snes_mpfaof90mod
    subroutine PorosityFunction(tdy,x,theta,dummy,ierr)
      implicit none
      TDy                    :: tdy
-     PetscReal, intent(in) :: x
-     PetscReal, intent(out):: theta
+     PetscReal, intent(in)  :: x(3)
+     PetscReal, intent(out) :: theta
      integer                :: dummy(*)
      PetscErrorCode :: ierr
 
@@ -49,8 +49,8 @@ end subroutine ResidualSaturation
   subroutine PorosityFunctionPFLOTRAN(tdy,x,theta,dummy,ierr)
     implicit none
     TDy                    :: tdy
-    PetscReal, intent(in) :: x
-    PetscReal, intent(out):: theta
+    PetscReal, intent(in)  :: x(3)
+    PetscReal, intent(out) :: theta
     integer                :: dummy(*)
     PetscErrorCode :: ierr
 
@@ -84,13 +84,13 @@ subroutine MaterialPropAlpha_PFLOTRAN(alpha)
     PetscReal, intent(out) :: alpha
     alpha = 1.d-4
   end subroutine MaterialPropAlpha_PFLOTRAN
-    
+
   subroutine MaterialPropM_PFLOTRAN(m)
     implicit none
     PetscReal, intent(out) :: m
     m = 0.3d0
   end subroutine MaterialPropM_PFLOTRAN
-    
+
   subroutine ResidualSat_PFLOTRAN(resSat)
     implicit none
     PetscReal resSat
@@ -100,8 +100,8 @@ end subroutine ResidualSat_PFLOTRAN
   subroutine PressureFunction(tdy,x,pressure,dummy,ierr)
     implicit none
     TDy                    :: tdy
-    PetscReal, intent(in) :: x(3)
-    PetscReal, intent(out):: pressure
+    PetscReal, intent(in)  :: x(3)
+    PetscReal, intent(out) :: pressure
     integer                :: dummy(*)
     PetscErrorCode :: ierr
 
@@ -160,6 +160,13 @@ program main
 
   call TDyInit(ierr);
   CHKERRA(ierr);
+
+  ! Register some functions.
+  call TDyRegisterFunction("p0", PressureFunction, ierr)
+  call TDyRegisterFunction("porosity", PorosityFunction, ierr)
+  call TDyRegisterFunction("porosity_pflotran", PorosityFunctionPFLOTRAN, ierr)
+  CHKERRA(ierr);
+
   call TDyCreate(tdy, ierr);
   CHKERRA(ierr);
   call TDySetDiscretizationMethod(tdy,MPFA_O,ierr);
@@ -305,7 +312,8 @@ program main
   CHKERRA(ierr);
 
   if (pflotran_consistent) then
-     call TDySetPorosityFunction(tdy,PorosityFunctionPFLOTRAN,0,ierr);
+!     call TDySetPorosityFunction(tdy,PorosityFunctionPFLOTRAN,0,ierr);
+     call TDySelectPorosityFunction(tdy,"porosity_pflotran",ierr);
      CHKERRA(ierr);
 
      do c = 1,ncell
@@ -328,7 +336,8 @@ program main
 
   else
 
-     call TDySetPorosityFunction(tdy,PorosityFunction,0,ierr);
+!     call TDySetPorosityFunction(tdy,PorosityFunction,0,ierr);
+     call TDySelectPorosityFunction(tdy,"porosity",ierr);
      CHKERRA(ierr);
 
      call TDySetBlockPermeabilityValuesLocal(tdy,ncell,index,blockPerm,ierr);
@@ -340,7 +349,8 @@ program main
   end if
 
   if (bc_type == MPFAO_DIRICHLET_BC .OR. bc_type == MPFAO_SEEPAGE_BC ) then
-     call TDySetBoundaryPressureFn(tdy,PressureFunction,0,ierr);
+!     call TDySetBoundaryPressureFn(tdy,PressureFunction,0,ierr);
+     call TDySelectBoundaryPressureFn(tdy,"p0",ierr);
      CHKERRQ(ierr)
   endif
 
diff --git a/include/tdycore.h b/include/tdycore.h
index fcfbca07..e5d0a8f2 100644
--- a/include/tdycore.h
+++ b/include/tdycore.h
@@ -74,8 +74,9 @@ PETSC_EXTERN PetscClassId TDY_CLASSID;
 PETSC_EXTERN PetscLogEvent TDy_ComputeSystem;
 
 /* ---------------------------------------------------------------- */
-PETSC_EXTERN PetscErrorCode TDyInit(int argc, char* argv[]);
+PETSC_EXTERN PetscErrorCode TDyInit(int, char*[]);
 PETSC_EXTERN PetscErrorCode TDyInitNoArguments(void);
+PETSC_EXTERN PetscErrorCode TDyOnFinalize(void (*)(void));
 PETSC_EXTERN PetscErrorCode TDyFinalize(void);
 
 PETSC_EXTERN PetscErrorCode TDyCreate(TDy*);
@@ -102,8 +103,12 @@ PETSC_EXTERN PetscErrorCode TDySetSoilDensityFunction(TDy,PetscErrorCode(*)(TDy,
 PETSC_EXTERN PetscErrorCode TDySetSoilSpecificHeatFunction(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*);
 
 // Set boundary and source-sink: via PETSc operations
+PetscErrorCode TDyRegisterFunction(const char*, PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*));
 PETSC_EXTERN PetscErrorCode TDySetForcingFunction(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*);
 PETSC_EXTERN PetscErrorCode TDySetEnergyForcingFunction(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*);
+PETSC_EXTERN PetscErrorCode TDySelectBoundaryPressureFn(TDy,const char*,void*);
+PETSC_EXTERN PetscErrorCode TDySelectBoundaryTemperatureFn(TDy,const char*,void*);
+PETSC_EXTERN PetscErrorCode TDySelectBoundaryVelocityFn(TDy,const char*,void*);
 PETSC_EXTERN PetscErrorCode TDySetBoundaryPressureFn(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*);
 PETSC_EXTERN PetscErrorCode TDySetBoundaryTemperatureFn(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*);
 PETSC_EXTERN PetscErrorCode TDySetBoundaryVelocityFn(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*);
diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90
index a5b060c1..48bd6d3c 100644
--- a/src/f90-mod/tdycoremod.F90
+++ b/src/f90-mod/tdycoremod.F90
@@ -437,6 +437,28 @@ subroutine TDyPostSolveSNESSolver(a,b,z)
      end subroutine TDyPostSolveSNESSolver
   end interface
 
+  abstract interface
+    subroutine TDyFunction(tdy, x, f, dummy, ierr)
+      use tdycoredef
+      TDy :: tdy
+      PetscReal, intent(in)  :: x(3)
+      PetscReal, intent(out) :: f
+      integer                :: dummy(*)
+      PetscErrorCode         :: ierr
+    end subroutine
+  end interface
+
+  ! We use GetRegFn to retrieve function pointers from the C registry.
+  interface
+    function GetRegFn(name, c_func) bind (c, name="TDyGetFunction") result(ierr)
+      use, intrinsic :: iso_c_binding
+      implicit none
+      type(c_ptr), value :: name
+      type(c_funptr) :: c_func
+      integer(c_int) :: ierr
+    end function
+  end interface
+
   contains
 
   subroutine TDyInit(ierr)
@@ -450,4 +472,145 @@ subroutine TDyInit(ierr)
      call TDyInitNoArguments(ierr)
   end subroutine TDyInit
 
+  subroutine TDyRegisterFunction(name, func, ierr)
+    use, intrinsic :: iso_c_binding
+    implicit none
+    character(len=*), intent(in)   :: name
+    procedure(TDyFunction)         :: func
+    PetscErrorCode                 :: ierr
+
+    interface
+      function RegisterFn(name, func) bind (c, name="TDyRegisterFunction") result(ierr)
+        use, intrinsic :: iso_c_binding
+        implicit none
+        type(c_ptr), value :: name
+        type(c_funptr), value :: func
+        integer(c_int) :: ierr
+      end function
+    end interface
+
+    ierr = RegisterFn(FtoCString(name), c_funloc(func))
+  end subroutine
+
+  subroutine TDySelectPorosityFunction(tdy, name, ierr)
+    use, intrinsic :: iso_c_binding
+    use tdycoredef
+    implicit none
+    TDy :: tdy
+    character(len=*), intent(in)   :: name
+    PetscErrorCode                 :: ierr
+
+    type(c_funptr)                  :: c_func
+    procedure(TDyFunction), pointer :: f_func
+
+    ierr = GetRegFn(FtoCString(name), c_func)
+    call c_f_procpointer(c_func, f_func)
+    call TDySetPorosityFunction(tdy, f_func, 0, ierr)
+  end subroutine
+
+  subroutine TDySelectForcingFunction(tdy, name, ierr)
+    use, intrinsic :: iso_c_binding
+    use tdycoredef
+    implicit none
+    TDy :: tdy
+    character(len=*), intent(in)   :: name
+    PetscErrorCode                 :: ierr
+
+    type(c_funptr)                  :: c_func
+    procedure(TDyFunction), pointer :: f_func
+
+    ierr = GetRegFn(FtoCString(name), c_func)
+    call c_f_procpointer(c_func, f_func)
+    call TDySetForcingFunction(tdy, f_func, 0, ierr)
+  end subroutine
+
+! Uncomment this when we are ready to set energy forcing fns in Fortran.
+!  subroutine TDySelectEnergyForcingFunction(tdy, name, ierr)
+!    use, intrinsic :: iso_c_binding
+!    use tdycoredef
+!    implicit none
+!    TDy :: tdy
+!    character(len=*), intent(in)   :: name
+!    PetscErrorCode                 :: ierr
+!
+!    type(c_funptr)                  :: c_func
+!    procedure(TDyFunction), pointer :: f_func
+!
+!    ierr = GetRegFn(FtoCString(name), c_func)
+!    call c_f_procpointer(c_func, f_func)
+!    call TDySetEnergyForcingFunction(tdy, f_func, 0, ierr)
+!  end subroutine
+
+! Uncomment this when we are ready to set permeability functions programmatically.
+!  subroutine TDySelectPermeabilityFunction(tdy, name, ierr)
+!    use, intrinsic :: iso_c_binding
+!    use tdycoredef
+!    implicit none
+!    TDy :: tdy
+!    character(len=*), intent(in)   :: name
+!    PetscErrorCode                 :: ierr
+!
+!    type(c_funptr)                  :: c_func
+!    procedure(TDyFunction), pointer :: f_func
+!
+!    ierr = GetRegFn(FtoCString(name), c_func)
+!    call c_f_procpointer(c_func, f_func)
+!    call TDySetPermeabilityFunction(tdy, f_func, 0, ierr)
+!  end subroutine
+
+  subroutine TDySelectBoundaryPressureFn(tdy, name, ierr)
+    use, intrinsic :: iso_c_binding
+    use tdycoredef
+    implicit none
+    TDy :: tdy
+    character(len=*), intent(in)   :: name
+    PetscErrorCode                 :: ierr
+
+    type(c_funptr)                  :: c_func
+    procedure(TDyFunction), pointer :: f_func
+
+    ierr = GetRegFn(FtoCString(name), c_func)
+    call c_f_procpointer(c_func, f_func)
+    call TDySetBoundaryPressureFn(tdy, f_func, 0, ierr)
+  end subroutine
+
+  subroutine TDySelectBoundaryVelocityFn(tdy, name, ierr)
+    use, intrinsic :: iso_c_binding
+    use tdycoredef
+    implicit none
+    TDy :: tdy
+    character(len=*), intent(in)   :: name
+    PetscErrorCode                 :: ierr
+
+    type(c_funptr)                  :: c_func
+    procedure(TDyFunction), pointer :: f_func
+
+    ierr = GetRegFn(FtoCString(name), c_func)
+    call c_f_procpointer(c_func, f_func)
+    call TDySetBoundaryVelocityFn(tdy, f_func, 0, ierr)
+  end subroutine
+
+  ! Here's a function that converts a Fortran string to a C string and
+  ! stashes it in TDycore's Fortran string registry. This allows us to
+  ! create more expressive (and standard) Fortran interfaces.
+  function FtoCString(f_string) result(c_string)
+     use, intrinsic :: iso_c_binding
+     implicit none
+     character(len=*), target :: f_string
+     character(len=:), pointer :: f_ptr
+     type(c_ptr) :: c_string
+
+     interface
+       function NewCString(f_str_ptr, f_str_len) bind (c, name="NewCString") result(c_string)
+         use, intrinsic :: iso_c_binding
+         type(c_ptr), value :: f_str_ptr
+         integer(c_int), value :: f_str_len
+         type(c_ptr) :: c_string
+       end function NewCString
+     end interface
+
+     f_ptr => f_string
+     c_string = NewCString(c_loc(f_ptr), len(f_string))
+   end function FtoCString
+
 end module tdycore
diff --git a/src/interface/ftn/tdycoref.c b/src/interface/ftn/tdycoref.c
index 97c591cd..79254367 100644
--- a/src/interface/ftn/tdycoref.c
+++ b/src/interface/ftn/tdycoref.c
@@ -130,6 +130,7 @@
 #define tdydestroy_                                    tdydestroy
 #endif
 
+// Legacy Fortran function pointer registry
 static struct {
   PetscFortranCallbackId porosity;
   PetscFortranCallbackId permeability;
@@ -142,7 +143,6 @@ static struct {
 #endif
 } _cb;
 
-
 #if defined(__cplusplus)
 extern "C" {
 #endif
diff --git a/src/tdyconditions.c b/src/tdyconditions.c
index b10cf179..95c43009 100644
--- a/src/tdyconditions.c
+++ b/src/tdyconditions.c
@@ -1,9 +1,56 @@
 #include <private/tdycoreimpl.h>
+#include <petsc/private/khash/khash.h>
 
 /*
   Boundary and source-sink conditions are set by PETSc operations
 */
 
+// Here's a registry of functions that can be used for boundary conditions and
+// forcing terms.
+typedef PetscErrorCode(*Function)(TDy, PetscReal*, PetscReal*, void*);
+KHASH_MAP_INIT_STR(TDY_FUNC_MAP, Function)
+static khash_t(TDY_FUNC_MAP)* funcs_ = NULL;
+
+// This function is called on finalization to destroy the function registry.
+static void DestroyFunctionRegistry() {
+  kh_destroy(TDY_FUNC_MAP, funcs_);
+}
+
+PetscErrorCode TDyRegisterFunction(const char* name, Function f) {
+  PetscFunctionBegin;
+  if (funcs_ == NULL) {
+    funcs_ = kh_init(TDY_FUNC_MAP);
+    TDyOnFinalize(DestroyFunctionRegistry);
+  }
+
+  int retval;
+  khiter_t iter = kh_put(TDY_FUNC_MAP, funcs_, name, &retval);
+  kh_val(funcs_, iter) = f;
+  PetscFunctionReturn(0);
+}
+
+PetscErrorCode TDyGetFunction(const char* name, Function* f) {
+  PetscFunctionBegin;
+  int ierr;
+
+  if (funcs_ != NULL) {
+    khiter_t iter = kh_get(TDY_FUNC_MAP, funcs_, name);
+    if (iter != kh_end(funcs_)) { // found it!
+      *f = kh_val(funcs_, iter);
+    } else {
+      ierr = -1;
+      SETERRQ(MPI_COMM_WORLD, ierr, "Function not found!");
+      return ierr;
+    }
+  } else {
+    ierr = -1;
+    SETERRQ(MPI_COMM_WORLD, ierr, "No functions have been registered!");
+    return ierr;
+  }
+  PetscFunctionReturn(0);
+}
+
+
 PetscErrorCode TDySetForcingFunction(TDy tdy, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*),void *ctx) {
   PetscFunctionBegin;
   if (f) tdy->ops->computeforcing = f;
@@ -18,6 +65,33 @@ PetscErrorCode TDySetEnergyForcingFunction(TDy tdy, PetscErrorCode(*f)(TDy,Petsc
   PetscFunctionReturn(0);
 }
 
+PetscErrorCode TDySelectBoundaryPressureFn(TDy tdy, const char* name, void* ctx) {
+  PetscFunctionBegin;
+  int ierr;
+  Function f;
+  ierr = TDyGetFunction(name, &f); CHKERRQ(ierr);
+  ierr = TDySetBoundaryPressureFn(tdy, f, ctx); CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+PetscErrorCode TDySelectBoundaryTemperatureFn(TDy tdy, const char* name, void* ctx) {
+  PetscFunctionBegin;
+  int ierr;
+  Function f;
+  ierr = TDyGetFunction(name, &f); CHKERRQ(ierr);
+  ierr = TDySetBoundaryTemperatureFn(tdy, f, ctx); CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+PetscErrorCode TDySelectBoundaryVelocityFn(TDy tdy, const char* name, void* ctx) {
+  PetscFunctionBegin;
+  int ierr;
+  Function f;
+  ierr = TDyGetFunction(name, &f); CHKERRQ(ierr);
+  ierr = TDySetBoundaryVelocityFn(tdy, f, ctx); CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
 PetscErrorCode TDySetBoundaryTemperatureFn(TDy tdy, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*),void *ctx) {
   PetscFunctionBegin;
   if (f) tdy->ops->compute_boundary_temperature = f;
diff --git a/src/tdycore.c b/src/tdycore.c
index ba66fd12..f8926a3f 100644
--- a/src/tdycore.c
+++ b/src/tdycore.c
@@ -140,13 +140,38 @@ PetscBool TDyInitialized(void) {
   return TDyPackageInitialized;
 }
 
+// A registry of functions called at shutdown. This can be used by all
+// subsystems via TDyOnFinalize().
+typedef void (*ShutdownFunc)(void);
+static ShutdownFunc *shutdown_funcs_ = NULL;
+static int num_shutdown_funcs_ = 0;
+static int shutdown_funcs_cap_ = 0;
+
+/// Call this to register a shutdown function that is called during TDyFinalize.
+PetscErrorCode TDyOnFinalize(void (*shutdown_func)(void)) {
+  if (shutdown_funcs_ == NULL) {
+    shutdown_funcs_cap_ = 32;
+    shutdown_funcs_ = malloc(sizeof(ShutdownFunc) * shutdown_funcs_cap_);
+  } else if (num_shutdown_funcs_ == shutdown_funcs_cap_) { // need more space!
+    shutdown_funcs_cap_ *= 2;
+    shutdown_funcs_ = realloc(shutdown_funcs_, sizeof(ShutdownFunc) * shutdown_funcs_cap_);
+  }
+  shutdown_funcs_[num_shutdown_funcs_] = shutdown_func;
+  ++num_shutdown_funcs_;
+
+  PetscFunctionReturn(0);
+}
+
 PetscErrorCode TDyFinalize() {
   PetscFunctionBegin;
 
-  // Dump timing information before we leave.
-  TDyWriteTimingProfile("tdycore_profile.csv");
-
-  TDyDestroyTimers();
+  // Call shutdown functions in reverse order, and destroy the list.
+  if (shutdown_funcs_ != NULL) {
+    for (int i = num_shutdown_funcs_-1; i >= 0; --i) {
+      shutdown_funcs_[i]();
+    }
+    free(shutdown_funcs_);
+  }
 
   // Finalize PETSc.
   PetscFinalize();
@@ -1460,3 +1485,39 @@ PetscErrorCode TDyCreateJacobian(TDy tdy) {
   PetscFunctionReturn(0);
 }
 
+// Here's a registry of C-backed strings created from Fortran.
+KHASH_SET_INIT_STR(TDY_STRING_SET)
+static khash_t(TDY_STRING_SET)* fortran_strings_ = NULL;
+
+// This function is called on finalization to destroy the Fortran ѕtring
+// registry.
+static void DestroyFortranStrings() {
+  kh_destroy(TDY_STRING_SET, fortran_strings_);
+}
+
+// Returns a newly-allocated C string for the given Fortran string pointer with
+// the given length. Resources for this string are managed by the running model.
+const char* NewCString(char* f_str_ptr, int f_str_len) {
+  if (fortran_strings_ == NULL) {
+    fortran_strings_ = kh_init(TDY_STRING_SET);
+    TDyOnFinalize(DestroyFortranStrings);
+  }
+
+  // Copy the Fortran character array to a C string.
+  char c_str[f_str_len+1];
+  memcpy(c_str, f_str_ptr, sizeof(char) * f_str_len);
+  c_str[f_str_len] = '\0';
+
+  // Does this string already exist?
+  khiter_t iter = kh_get(TDY_STRING_SET, fortran_strings_, c_str);
+  if (iter != kh_end(fortran_strings_)) { // yep
+    return kh_key(fortran_strings_, iter);
+  } else { // nope
+    int retval;
+    const char* str = malloc(sizeof(char) * (f_str_len+1));
+    strcpy((char*)str, c_str);
+    iter = kh_put(TDY_STRING_SET, fortran_strings_, str, &retval);
+    return str;
+  }
+}
+
diff --git a/src/tdytimers.c b/src/tdytimers.c
index 3a2649b5..cc62451a 100644
--- a/src/tdytimers.c
+++ b/src/tdytimers.c
@@ -24,6 +24,14 @@ khash_t(TDY_PROFILING_MD_MAP)* TDY_PROFILING_METADATA = NULL;
 // Are timers enabled?
 static PetscBool timersEnabled_ = PETSC_FALSE;
 
+// This function is called at finalization time.
+static void ShutDownTimers() {
+  // Dump timing information before we leave.
+  TDyWriteTimingProfile("tdycore_profile.csv");
+
+  TDyDestroyTimers();
+}
+
 PetscErrorCode TDyInitTimers() {
   // Register timers table.
   if (TDY_TIMERS == NULL) {
@@ -46,6 +54,9 @@ PetscErrorCode TDyInitTimers() {
   ierr = TDyAddProfilingStage("TDycore Stepping"); CHKERRQ(ierr);
   ierr = TDyAddProfilingStage("TDycore I/O"); CHKERRQ(ierr);
 
+  // Register our shutdown function.
+  TDyOnFinalize(ShutDownTimers);
+
   return ierr;
 }