diff --git a/applications/solvers/combustion/coldEngineFoam/Make/options b/applications/solvers/combustion/coldEngineFoam/Make/options
index a59f3b8053149415dcc2d7972afaf005f0e18971..25217ab23f79c72d9f81c2aba8eb51c0ae138c30 100644
--- a/applications/solvers/combustion/coldEngineFoam/Make/options
+++ b/applications/solvers/combustion/coldEngineFoam/Make/options
@@ -7,7 +7,10 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lengine \
@@ -16,4 +19,7 @@ EXE_LIBS = \
     -lcompressibleLESModels \
     -lfluidThermophysicalModels \
     -lspecie \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lsampling \
+    -lmeshTools \
+    -lfieldSources
diff --git a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
index 87474e778ffe25124a56264d8a95bcaae382d6aa..401690c3e49a263443da462ebca49c6a11863329 100644
--- a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
+++ b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
@@ -35,6 +35,7 @@ Description
 #include "psiThermo.H"
 #include "turbulenceModel.H"
 #include "OFstream.H"
+#include "IObasicSourceList.H"
 #include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/combustion/coldEngineFoam/createFields.H b/applications/solvers/combustion/coldEngineFoam/createFields.H
index f78e5bde2969cd6935e3bccec61887da6f9a2905..72e88435aceacbf2d6fa546b068157dc7caf0de3 100644
--- a/applications/solvers/combustion/coldEngineFoam/createFields.H
+++ b/applications/solvers/combustion/coldEngineFoam/createFields.H
@@ -69,3 +69,6 @@
 
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));
+
+    Info<< "Creating sources\n" << endl;
+    IObasicSourceList sources(mesh);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/EEqn.H b/applications/solvers/compressible/rhoPimpleFoam/EEqn.H
index 6288b8a0af0842842ab061aa1148c804bfb9b7a4..1cc6ed584c6c3806864cc0a0f461182796b9799a 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/EEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/EEqn.H
@@ -16,9 +16,12 @@
           : -dpdt
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        sources(rho, he)
     );
 
     EEqn.relax();
+    sources.constrain(EEqn);
     EEqn.solve();
 
     thermo.correct();
diff --git a/applications/solvers/compressible/rhoPimpleFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/Make/options
index f21e7394e4db3c8d5cf3770ecebbe46e30e4ec50..c96a146eca26d9bfdb69c0526632707b4cea13e4 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/Make/options
@@ -2,7 +2,10 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -11,4 +14,6 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
index 0f5ec2487b30854d2982e39196ea194bd9d32c4e..397e8930357dfd548a96f62b9f63acfc943f2423 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
@@ -9,8 +9,10 @@ tmp<fvVectorMatrix> UEqn
 
 UEqn().relax();
 
+sources.constrain(UEqn());
+
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn() == -fvc::grad(p) + sources(rho, U));
     K = 0.5*magSqr(U);
 }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
index 67cc0c3e456573e71fe9f4937a5fd5d3e9253692..8b3cfd10b9a1180d574c1610ccec487b3ab5ae68 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
@@ -69,3 +69,6 @@
 
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));
+
+    Info<< "Creating sources\n" << endl;
+    IObasicSourceList sources(mesh);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index d40c66d40252239b172f7d493da57b51c92e7d93..5f28a3456d82c8d6c1b8f05d90b1519e9a5f5931 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -5,7 +5,7 @@ rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
+HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -33,8 +33,12 @@ if (pimple.transonic())
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
           - fvm::laplacian(Dp, p)
+          ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -65,8 +69,12 @@ else
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
           - fvm::laplacian(Dp, p)
+          ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -92,6 +100,7 @@ Info<< "rho max/min : " << max(rho).value()
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index adb6882268a7315be173b4b8b5c175937fc7ada4..f9b8a901f721265fe10cd8c478651f8334ba8497 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -38,6 +38,7 @@ Description
 #include "turbulenceModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options
index 6d3c41f0dfa6bde2b2e91e36a9458c172ca6cf13..c78d474b4fb3d4649d7afad5957132ee2da37253 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options
@@ -3,7 +3,10 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -12,4 +15,6 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
index 75ed20e271981cb2d350171947c2bc680593b389..628eb71f2961392d5c1173fbbf30aeb07cd5f81b 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
@@ -7,7 +7,7 @@ volScalarField rAU(1.0/UEqn().A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
+HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
 if (pimple.nCorrPIMPLE() <= 1)
 {
@@ -44,11 +44,15 @@ if (pimple.transonic())
           + fvm::div(phid, p)
           + fvc::div(phic)
           - fvm::laplacian(Dp, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
         // Relax the pressure equation to maintain diagonal dominance
         pEqn.relax();
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve();
 
         if (pimple.finalNonOrthogonalIter())
@@ -81,8 +85,12 @@ else
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
           - fvm::laplacian(Dp, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve();
 
         if (pimple.finalNonOrthogonalIter())
@@ -100,6 +108,7 @@ p.relax();
 
 U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C
index cf8aed0bcad3870312599f76237f1ae6cd6f8be2..7e1664dbeac745cf8a24e634cfd72b87e4f41213 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C
@@ -38,6 +38,7 @@ Description
 #include "turbulenceModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/Make/options
index 5c83c910d95a4aa6b5de30075e4f764309a4bec3..2976366c08ac366a8c0abe6c3291714ed86f1908 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/Make/options
@@ -5,7 +5,9 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -14,4 +16,6 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
index 3a8eb884ba8ca7eaed6d7e8a6ea28ca834393447..82ba5cbac17371b9147b4057d5f657426cff1afc 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
@@ -38,7 +38,8 @@ Description
 #include "psiThermo.H"
 #include "turbulenceModel.H"
 #include "MRFZones.H"
-#include "porousZones.H"
+#include "IOporosityModelList.H"
+#include "IObasicSourceList.H"
 #include "fvcSmooth.H"
 #include "pimpleControl.H"
 #include "bound.H"
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/options
index f01ebda533b3a30f359683e4b72fde503798f830..c78d474b4fb3d4649d7afad5957132ee2da37253 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/options
@@ -4,7 +4,9 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -13,4 +15,6 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H
index c24b1f587a59bac6cd7b305b5dc6ccdf55986613..9651610157c483596253639cb1380be1cc8e2f9b 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H
@@ -11,9 +11,12 @@ tmp<fvVectorMatrix> UEqn
 UEqn().relax();
 
 mrfZones.addCoriolis(rho, UEqn());
+
 pZones.addResistance(UEqn());
 
+sources.constrain(UEqn());
+
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn() == -fvc::grad(p) + sources(rho, U));
 }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/createZones.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/createZones.H
index 487a55d8eddc44a6fe987ec15b06165ed4b4008e..d2522b8c752ada18d083b3ebf450d0ae553aac64 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/createZones.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/createZones.H
@@ -1,5 +1,5 @@
     MRFZones mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U);
 
-    porousZones pZones(mesh);
+    IOporosityModelList pZones(mesh);
     Switch pressureImplicitPorosity(false);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
index b68823c6e3d1beaa68743ba7a973e5a17b2d44eb..304798b9662d6c7bb8dd3bfc0bded47e552db19e 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
@@ -5,7 +5,7 @@ rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
+HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -34,8 +34,12 @@ if (pimple.transonic())
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
           - fvm::laplacian(Dp, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -68,8 +72,12 @@ else
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
           - fvm::laplacian(Dp, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -94,6 +102,7 @@ Info<< "rho max/min : " << max(rho).value()
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
index 16e6c0ad746313bcd41a51ab34d48868dfaebc71..ac5dc1eadcb3d6c9108160928e93fa67d6e4569d 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
@@ -38,7 +38,8 @@ Description
 #include "turbulenceModel.H"
 #include "bound.H"
 #include "MRFZones.H"
-#include "porousZones.H"
+#include "IOporosityModelList.H"
+#include "IObasicSourceList.H"
 #include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/compressible/rhoSimpleFoam/EEqn.H b/applications/solvers/compressible/rhoSimpleFoam/EEqn.H
index e4c79b40cc5685ddbeeab18c32e8611d745656b8..e496906d082b96a7e2dc8a9aa731dba5372608aa 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/EEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/EEqn.H
@@ -10,9 +10,12 @@
           : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        sources(rho, he)
     );
 
     EEqn.relax();
+    sources.constrain(EEqn);
     EEqn.solve();
 
     thermo.correct();
diff --git a/applications/solvers/compressible/rhoSimpleFoam/Make/options b/applications/solvers/compressible/rhoSimpleFoam/Make/options
index 339cc53bd9e9e2e045b859cec7a26b8076168565..741dc4f822bc8343080f8afa237be831ba110e00 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/Make/options
@@ -3,7 +3,10 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -11,4 +14,6 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lsampling \
+    -lmeshTools \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
index dea35b76575477e953c521d77f6bb12cf28c3dc0..3f4b2fd0217cb800ae9c25793630f8266ca4d8d4 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@@ -60,3 +60,6 @@
     );
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+    
+    Info<< "Creating sources\n" << endl;
+    IObasicSourceList sources(mesh);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 4d7334c5f8640b44f31a06802f845a7d7c730407..2ea228f172d80bc4d4a330662723e9754840c0a2 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -6,7 +6,7 @@
 
     volScalarField rAU(1.0/UEqn().A());
     volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn().H();
+    HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
     UEqn.clear();
 
@@ -25,12 +25,16 @@
             fvScalarMatrix pEqn
             (
                 fvm::div(phid, p)
-                - fvm::laplacian(rho*rAU, p)
+              - fvm::laplacian(rho*rAU, p)
+              ==
+                sources(psi, p, rho.name())
             );
 
             // Relax the pressure equation to ensure diagonal-dominance
             pEqn.relax();
 
+            sources.constrain(pEqn, rho.name());
+
             pEqn.setReference(pRefCell, pRefValue);
 
             pEqn.solve();
@@ -56,11 +60,15 @@
             fvScalarMatrix pEqn
             (
                 fvc::div(phiHbyA)
-                - fvm::laplacian(rho*rAU, p)
+              - fvm::laplacian(rho*rAU, p)
+              ==
+                sources(psi, p, rho.name())
             );
 
             pEqn.setReference(pRefCell, pRefValue);
 
+            sources.constrain(pEqn, rho.name());
+
             pEqn.solve();
 
             if (simple.finalNonOrthogonalIter())
@@ -78,6 +86,7 @@
 
     U = HbyA - rAU*fvc::grad(p);
     U.correctBoundaryConditions();
+    sources.correct(U);
 
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H
index ff467c0382cfcb444f66d64cbbcb549c0dac1736..e496906d082b96a7e2dc8a9aa731dba5372608aa 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H
@@ -10,11 +10,12 @@
           : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        sources(rho, he)
     );
 
-    pZones.addEnergySource(thermo, rho, EEqn);
-
     EEqn.relax();
+    sources.constrain(EEqn);
     EEqn.solve();
 
     thermo.correct();
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/Make/options b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/Make/options
index 2d328b3ed1f86ecd981c6c3fa3bc0225bd67aa4a..4ba7ed1f7d10bc843dcc0c9d55cd310b08bf98b0 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/Make/options
@@ -1,18 +1,20 @@
 EXE_INC = \
     -I.. \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermalPorousZone/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
-    -lthermalPorousZone \
     -lspecie \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lsampling \
+    -lmeshTools \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/UEqn.H
index c0e087b536c3ffdffb01214e7204b13abfba8a46..52e15e07b5143f0e078cec89c9d1e57ebf04108a 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/UEqn.H
@@ -24,11 +24,13 @@
         trTU = inv(tTU());
         trTU().rename("rAU");
 
+        sources.constrain(UEqn());
+
         volVectorField gradp(fvc::grad(p));
 
         for (int UCorr=0; UCorr<nUCorr; UCorr++)
         {
-            U = trTU() & (UEqn().H() - gradp);
+            U = trTU() & ((UEqn() == sources(rho, U))().H() - gradp);
         }
         U.correctBoundaryConditions();
     }
@@ -36,7 +38,9 @@
     {
         pZones.addResistance(UEqn());
 
-        solve(UEqn() == -fvc::grad(p));
+        sources.constrain(UEqn());
+
+        solve(UEqn() == -fvc::grad(p) + sources(rho, U));
 
         trAU = 1.0/UEqn().A();
         trAU().rename("rAU");
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H
index 4fff74d224cadb920dec4f6c42b936dff55f0452..ad924283eef67e97e6fcb193e0a12a6ce834ce9d 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H
@@ -59,3 +59,6 @@
     );
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+    Info<< "Creating sources\n" << endl;
+    IObasicSourceList sources(mesh);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createZones.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createZones.H
index 0ed5afa2745897eea93992ba7cad5c8e6a5c6fc2..4eb2cb193cc43ea551ec1a0e03548bd8f498bdf2 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createZones.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createZones.H
@@ -1,12 +1,12 @@
     MRFZones mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U);
 
-    thermalPorousZones pZones(mesh);
+    IOporosityModelList pZones(mesh);
     Switch pressureImplicitPorosity(false);
 
     // nUCorrectors used for pressureImplicitPorosity
     int nUCorr = 0;
-    if (pZones.size())
+    if (pZones.active())
     {
         // nUCorrectors for pressureImplicitPorosity
         simple.dict().readIfPresent("nUCorrectors", nUCorr);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
index 81b194823cb3c22012db5117ce7b7f7f5382b604..708449eb3d872a522f3a6c234937cef21aafe782 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
@@ -4,15 +4,17 @@
     rho = min(rho, rhoMax);
     rho.relax();
 
+    const volScalarField& psi = thermo.psi();
+
     volVectorField HbyA("HbyA", U);
 
     if (pressureImplicitPorosity)
     {
-        HbyA = trTU() & UEqn().H();
+        HbyA = trTU() & (UEqn() == sources(rho, U))().H();
     }
     else
     {
-        HbyA = trAU()*UEqn().H();
+        HbyA = trAU()*(UEqn() == sources(rho, U))().H();
     }
 
     UEqn.clear();
@@ -35,15 +37,29 @@
 
         if (pressureImplicitPorosity)
         {
-            tpEqn = (fvm::laplacian(rho*trTU(), p) == fvc::div(phiHbyA));
+            tpEqn =
+            (
+                fvm::laplacian(rho*trTU(), p)
+              + sources(psi, p, rho.name())
+             ==
+                fvc::div(phiHbyA)
+            );
         }
         else
         {
-            tpEqn = (fvm::laplacian(rho*trAU(), p) == fvc::div(phiHbyA));
+            tpEqn =
+            (
+                fvm::laplacian(rho*trAU(), p)
+              + sources(psi, p, rho.name())
+             ==
+                fvc::div(phiHbyA)
+            );
         }
 
         tpEqn().setReference(pRefCell, pRefValue);
 
+        sources.constrain(tpEqn(), rho.name());
+
         tpEqn().solve();
 
         if (simple.finalNonOrthogonalIter())
@@ -67,12 +83,12 @@
     }
 
     U.correctBoundaryConditions();
+    sources.correct(U);
 
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
     {
-        const volScalarField& psi = thermo.psi();
         p += (initialMass - fvc::domainIntegrate(psi*p))
             /fvc::domainIntegrate(psi);
     }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
index 3898295b915061669daa94fa773218a5877ef357..b552cadfcef253000d0e5352d941ae912d65348b 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
@@ -35,7 +35,8 @@ Description
 #include "rhoThermo.H"
 #include "RASModel.H"
 #include "MRFZones.H"
-#include "thermalPorousZones.H"
+#include "IObasicSourceList.H"
+#include "IOporosityModelList.H"
 #include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
index 312196583eaae9a90ce43b22350f2ecb5d4ac383..efc1067db233ca92b1a286fad6f799a6c77ea6c7 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
@@ -34,6 +34,7 @@ Description
 #include "psiThermo.H"
 #include "RASModel.H"
 #include "simpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options
index 8a8553525f39f3e90f00669cec5ee8eb4633f149..facff4ac649182a75207b28def3408506e51bc90 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options
@@ -4,11 +4,16 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
     -lspecie \
     -lcompressibleRASModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lsampling \
+    -lmeshTools \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
index 8c7405c346ffdebce86d052fe60be505cb52d584..1ffacc1afaf9c5a2495e8e321c84f65a7ea592ac 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
@@ -7,7 +7,7 @@ volScalarField rAU(1.0/UEqn().A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
+HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
 UEqn.clear();
 
@@ -38,11 +38,15 @@ if (simple.transonic())
             fvm::div(phid, p)
           + fvc::div(phic)
           - fvm::laplacian(Dp, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
         // Relax the pressure equation to maintain diagonal dominance
         pEqn.relax();
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.setReference(pRefCell, pRefValue);
 
         pEqn.solve();
@@ -74,8 +78,12 @@ else
         (
             fvc::div(phiHbyA)
           - fvm::laplacian(Dp, p)
+          ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.setReference(pRefCell, pRefValue);
 
         pEqn.solve();
@@ -96,6 +104,7 @@ p.relax();
 
 U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C
index f707570f60cd0e0a73162bd7ebabb4857b2b325e..41e15dbfa842e5c8b75228e95c37edf5cc3d2234 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C
@@ -36,6 +36,7 @@ Description
 #include "mixedFvPatchFields.H"
 #include "bound.H"
 #include "simpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options b/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
index f6e12a3b7ce6dcd564b0eacedd671e7f2d1b75c1..5d86b094de32587a79d2bade587ef6f47ccbeb90 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
@@ -2,7 +2,10 @@ EXE_INC = \
     -I../../compressible/rhoPimpleFoam \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
 
 EXE_LIBS = \
     -lmeshTools \
@@ -11,4 +14,8 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lsampling \
+    -lmeshTools \
+    -lfieldSources
+
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
index 9597f6ac08ddcac42bba70642fe61c013f1c87ed..d8f3bb56ce13a0455c52f7ccbe598f353ea21b2e 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
@@ -22,6 +22,7 @@
                   - fvc::snGrad(p_rgh)
                 )*mesh.magSf()
             )
+          + sources(rho, U)
         );
         K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index dfa88b59a9aebb32486807f80705b1381925946c..79d3fb7f2f7668b4dbf1696c25ac7f38a6a0c25b 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -37,6 +37,7 @@ Description
 #include "rhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
+#include "IObasicSourceList.H"
 #include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
index f6b183e3fcbf95c113c1a81d7da789c70478ecaa..9838978028edd00a3f9fbce7563b59bce0dcc9cb 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
@@ -89,3 +89,6 @@
 
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));
+
+    Info<< "Creating sources\n" << endl;
+    IObasicSourceList sources(mesh);
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index d67011dba107dfaa529a214691903cfb81db6dc7..86c314c9c4bab41f858844d961d85eeaac511ccf 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -9,7 +9,7 @@
     surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn.H();
+    HbyA = rAU*(UEqn == sources(rho, U))().H();
 
     surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
@@ -28,6 +28,8 @@
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
       + fvc::div(phiHbyA)
+      ==
+        sources(psi, p_rgh, rho.name())
     );
 
     while (pimple.correctNonOrthogonal())
@@ -38,6 +40,8 @@
           - fvm::laplacian(rhorAUf, p_rgh)
         );
 
+        sources.constrain(p_rghEqn, rho.name());
+
         p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -52,6 +56,7 @@
             // calculated from the relaxed pressure
             U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
+            sources.correct(U);
             K = 0.5*magSqr(U);
         }
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
index 25e661dfc3260667acb167e9781651dc33e6b1e1..557017ccad4c7254298b1cc13622623b11d846d0 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@@ -45,7 +45,7 @@ Description
 #include "solidRegionDiffNo.H"
 #include "solidThermo.H"
 #include "radiationModel.H"
-#include "porousZones.H"
+#include "IOporosityModelList.H"
 #include "IObasicSourceList.H"
 
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
index c875d850cd8c16f2c4344c7648a104fdfdde703f..2f4b801cb398d5002de1ca0ec748199efc80eccf 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
@@ -36,7 +36,7 @@ Description
 #include "regionProperties.H"
 #include "solidThermo.H"
 #include "radiationModel.H"
-#include "porousZones.H"
+#include "IOporosityModelList.H"
 #include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
index c34b547350610a684bbf77d2de057c44782d9bd4..3d66fa798b6bd653efc93ea38b1e35347b72bd9c 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
@@ -18,7 +18,7 @@
     PtrList<dimensionedScalar> rhoMin(fluidRegions.size());
 
     PtrList<IObasicSourceList> heatSources(fluidRegions.size());
-    PtrList<porousZones> porousZonesFluid(fluidRegions.size());
+    PtrList<IOporosityModelList> porousZonesFluid(fluidRegions.size());
 
     // Populate fluid field pointer lists
     forAll(fluidRegions, i)
@@ -205,7 +205,7 @@
         porousZonesFluid.set
         (
             i,
-            new porousZones(fluidRegions[i])
+            new IOporosityModelList(fluidRegions[i])
         );
     }
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
index 3abad394170710cf0748f64e42665bd773574e0c..4f15d8c6192c87acd9401c2ce0184d5c19faa1cb 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
@@ -14,7 +14,7 @@
 
     IObasicSourceList& sources = heatSources[i];
 
-    const porousZones& pZones = porousZonesFluid[i];
+    const IOporosityModelList& pZones = porousZonesFluid[i];
 
     const dimensionedScalar initialMass
     (
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index 81cb541ea4b646f33b2432334043cf81515d4606..a2021a5edc66f83255886fc86c18535c4c2e689d 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -15,7 +15,7 @@
     List<scalar> initialMassFluid(fluidRegions.size());
 
     PtrList<IObasicSourceList> heatSources(fluidRegions.size());
-    PtrList<porousZones> porousZonesFluid(fluidRegions.size());
+    PtrList<IOporosityModelList> porousZonesFluid(fluidRegions.size());
 
     // Populate fluid field pointer lists
     forAll(fluidRegions, i)
@@ -202,6 +202,6 @@
         porousZonesFluid.set
         (
             i,
-            new porousZones(fluidRegions[i])
+            new IOporosityModelList(fluidRegions[i])
         );
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index 04986d1289847fc88d09c615e8a24702ce879d2c..28f0ea96834746cf0dfafaea0951bbd8a5b5503c 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -22,7 +22,7 @@
 
     IObasicSourceList& sources = heatSources[i];
 
-    const porousZones& pZones = porousZonesFluid[i];
+    const IOporosityModelList& pZones = porousZonesFluid[i];
 
     const dimensionedScalar initialMass
     (
diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/createPorousZones.H b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/createPorousZones.H
index e4614b7063469dc05af8911a540300939c1c3a96..5c6dbea19c4122682311089f3a5f3ddb16d870f6 100644
--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/createPorousZones.H
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/createPorousZones.H
@@ -1,9 +1,9 @@
-    porousZones pZones(mesh);
+    IOporosityModelList pZones(mesh);
     Switch pressureImplicitPorosity(false);
 
     // nUCorrectors used for pressureImplicitPorosity
     int nUCorr = 0;
-    if (pZones.size())
+    if (pZones.active())
     {
         // nUCorrectors for pressureImplicitPorosity
         nUCorr = simple.dict().lookupOrDefault<int>("nUCorrectors", 0);
diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C
index ef10119a9e372f0f8c6023d75c5f14d5bbf3bc41..bde78bcbe155eb56a0a213c688e1ee317a1f0d4a 100644
--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C
@@ -33,7 +33,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
-#include "porousZones.H"
+#include "IOporosityModelList.H"
 #include "simpleControl.H"
 #include "IObasicSourceList.H"
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C
index fabb0bba6090f482877f029a4e97cab13fad88a2..4ea17727d9e6cc9088a8952cf6c2175841fd6cf5 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C
@@ -40,7 +40,7 @@ Description
 #include "basicReactingMultiphaseCloud.H"
 #include "rhoCombustionModel.H"
 #include "radiationModel.H"
-#include "porousZones.H"
+#include "IOporosityModelList.H"
 #include "IObasicSourceList.H"
 #include "SLGThermo.H"
 #include "fvcSmooth.H"
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createPorousZones.H b/applications/solvers/lagrangian/reactingParcelFoam/createPorousZones.H
index 90506856d2a072cad5ee3b4be584ee0aa42fda81..05f0f044ce4314fc37f7c64f0a8808099199a213 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/createPorousZones.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createPorousZones.H
@@ -1,3 +1,3 @@
     Info<< "Creating porous zones" << nl << endl;
 
-    porousZones pZones(mesh);
+    IOporosityModelList pZones(mesh);
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index 5b902273234aa924389682cd19bccb285bd1448d..1bcdf7c92f917c1a1514b50ab0ccbe95f1d38f15 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -10,7 +10,7 @@
     HbyA = rAU*(UEqn == sources(rho, U))().H();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
-    if (pZones.size() == 0)
+    if (!pZones.active())
     {
         // ddtPhiCorr only used without porosity
         phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi);
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
index 8ef6838c597de6b67c85ee7175f0a64643f44319..a8f43285433a4494c275f66cd556c65fddbc7f4d 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
@@ -44,7 +44,7 @@ Description
 #include "basicReactingMultiphaseCloud.H"
 #include "rhoCombustionModel.H"
 #include "radiationModel.H"
-#include "porousZones.H"
+#include "IOporosityModelList.H"
 #include "IObasicSourceList.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"
diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/createPorousZones.H b/applications/solvers/multiphase/interFoam/porousInterFoam/createPorousZones.H
index 430b466aac6e6cd3c634a3e75257d9efa477561c..983906074aa69f21dcf5a16c5f4ce8622909ba34 100644
--- a/applications/solvers/multiphase/interFoam/porousInterFoam/createPorousZones.H
+++ b/applications/solvers/multiphase/interFoam/porousInterFoam/createPorousZones.H
@@ -1 +1 @@
-    porousZones pZones(mesh);
+    IOporosityModelList pZones(mesh);
diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
index 891736c77fcc0f4872f92f76558b17b56f7bfe72..244ea75faef0e11ab52aaa6f4c40a42ef6db4004 100644
--- a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -44,7 +44,7 @@ Description
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
-#include "porousZones.H"
+#include "IOporosityModelList.H"
 #include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/fieldSources/Make/files b/src/fieldSources/Make/files
index eced302bfded50b2988169b216e37a6ac27d068a..41fecea8008b9daed1f4be8e3654e50121066bff 100644
--- a/src/fieldSources/Make/files
+++ b/src/fieldSources/Make/files
@@ -1,34 +1,35 @@
-basicSource/basicSource/basicSource.C
-basicSource/basicSource/basicSourceIO.C
-basicSource/basicSource/basicSourceList.C
-basicSource/basicSource/IObasicSourceList.C
+basicSource/basicSource.C
+basicSource/basicSourceIO.C
+basicSource/basicSourceList.C
+basicSource/IObasicSourceList.C
 
-basicSource/pressureGradientExplicitSource/pressureGradientExplicitSource.C
-basicSource/pressureGradientExplicitSource/pressureGradientExplicitSourceIO.C
+general/explicitSource/explicitSource.C
+general/explicitSetValue/explicitSetValue.C
+general/codedSource/codedSource.C
 
-basicSource/explicitSource/explicitSource.C
-basicSource/explicitSetValue/explicitSetValue.C
+derived/pressureGradientExplicitSource/pressureGradientExplicitSource.C
+derived/pressureGradientExplicitSource/pressureGradientExplicitSourceIO.C
 
-basicSource/rotorDiskSource/rotorDiskSource.C
-basicSource/rotorDiskSource/bladeModel/bladeModel.C
-basicSource/rotorDiskSource/profileModel/profileModel.C
-basicSource/rotorDiskSource/profileModel/profileModelList.C
-basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C
-basicSource/rotorDiskSource/profileModel/series/seriesProfile.C
-basicSource/rotorDiskSource/trimModel/trimModel/trimModel.C
-basicSource/rotorDiskSource/trimModel/trimModel/trimModelNew.C
-basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C
-basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
+derived/rotorDiskSource/rotorDiskSource.C
+derived/rotorDiskSource/bladeModel/bladeModel.C
+derived/rotorDiskSource/profileModel/profileModel.C
+derived/rotorDiskSource/profileModel/profileModelList.C
+derived/rotorDiskSource/profileModel/lookup/lookupProfile.C
+derived/rotorDiskSource/profileModel/series/seriesProfile.C
+derived/rotorDiskSource/trimModel/trimModel/trimModel.C
+derived/rotorDiskSource/trimModel/trimModel/trimModelNew.C
+derived/rotorDiskSource/trimModel/fixed/fixedTrim.C
+derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
 
-basicSource/actuationDiskSource/actuationDiskSource.C
-basicSource/radialActuationDiskSource/radialActuationDiskSource.C
+derived/actuationDiskSource/actuationDiskSource.C
+derived/radialActuationDiskSource/radialActuationDiskSource.C
 
-interRegion = basicSource/interRegionHeatTransferModel
+interRegion = derived/interRegionHeatTransferModel
 $(interRegion)/interRegionHeatTransferModel/interRegionHeatTransferModel.C
 $(interRegion)/constantHeatTransfer/constantHeatTransfer.C
 $(interRegion)/tabulatedHeatTransfer/tabulatedHeatTransfer.C
 $(interRegion)/variableHeatTransfer/variableHeatTransfer.C
 
-basicSource/codedSource/codedSource.C
+derived/fixedTemperatureSource/fixedTemperatureSource.C
 
 LIB = $(FOAM_LIBBIN)/libfieldSources
diff --git a/src/fieldSources/basicSource/basicSource/IObasicSourceList.C b/src/fieldSources/basicSource/IObasicSourceList.C
similarity index 100%
rename from src/fieldSources/basicSource/basicSource/IObasicSourceList.C
rename to src/fieldSources/basicSource/IObasicSourceList.C
diff --git a/src/fieldSources/basicSource/basicSource/IObasicSourceList.H b/src/fieldSources/basicSource/IObasicSourceList.H
similarity index 100%
rename from src/fieldSources/basicSource/basicSource/IObasicSourceList.H
rename to src/fieldSources/basicSource/IObasicSourceList.H
diff --git a/src/fieldSources/basicSource/basicSource/basicSource.C b/src/fieldSources/basicSource/basicSource.C
similarity index 98%
rename from src/fieldSources/basicSource/basicSource/basicSource.C
rename to src/fieldSources/basicSource/basicSource.C
index 2d218de752afdf67cd63e169d01a2e714924edd3..9ebcdc5834f30cf42cab164395a404db33206693 100644
--- a/src/fieldSources/basicSource/basicSource/basicSource.C
+++ b/src/fieldSources/basicSource/basicSource.C
@@ -56,6 +56,12 @@ namespace Foam
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+bool Foam::basicSource::alwaysApply() const
+{
+    return false;
+}
+
+
 void Foam::basicSource::setSelection(const dictionary& dict)
 {
     switch (selectionMode_)
@@ -314,6 +320,7 @@ Foam::basicSource::~basicSource()
     }
 }
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::basicSource::isActive()
@@ -341,6 +348,11 @@ bool Foam::basicSource::isActive()
 
 Foam::label Foam::basicSource::applyToField(const word& fieldName) const
 {
+    if (alwaysApply())
+    {
+        return 0;
+    }
+
     forAll(fieldNames_, i)
     {
         if (fieldNames_[i] == fieldName)
diff --git a/src/fieldSources/basicSource/basicSource/basicSource.H b/src/fieldSources/basicSource/basicSource.H
similarity index 99%
rename from src/fieldSources/basicSource/basicSource/basicSource.H
rename to src/fieldSources/basicSource/basicSource.H
index ca105f90b1883bf444bc9588556e123d3d306101..69e9dcf408b9588b8ef24eea9007afca51bb7e42 100644
--- a/src/fieldSources/basicSource/basicSource/basicSource.H
+++ b/src/fieldSources/basicSource/basicSource.H
@@ -150,6 +150,9 @@ protected:
 
     // Protected functions
 
+        //- Flag to bypass the apply flag list checking
+        virtual bool alwaysApply() const;
+
         //- Set the cellSet or points selection
         void setSelection(const dictionary& dict);
 
diff --git a/src/fieldSources/basicSource/basicSource/basicSourceI.H b/src/fieldSources/basicSource/basicSourceI.H
similarity index 100%
rename from src/fieldSources/basicSource/basicSource/basicSourceI.H
rename to src/fieldSources/basicSource/basicSourceI.H
diff --git a/src/fieldSources/basicSource/basicSource/basicSourceIO.C b/src/fieldSources/basicSource/basicSourceIO.C
similarity index 100%
rename from src/fieldSources/basicSource/basicSource/basicSourceIO.C
rename to src/fieldSources/basicSource/basicSourceIO.C
diff --git a/src/fieldSources/basicSource/basicSource/basicSourceList.C b/src/fieldSources/basicSource/basicSourceList.C
similarity index 100%
rename from src/fieldSources/basicSource/basicSource/basicSourceList.C
rename to src/fieldSources/basicSource/basicSourceList.C
diff --git a/src/fieldSources/basicSource/basicSource/basicSourceList.H b/src/fieldSources/basicSource/basicSourceList.H
similarity index 99%
rename from src/fieldSources/basicSource/basicSource/basicSourceList.H
rename to src/fieldSources/basicSource/basicSourceList.H
index 6a7a0126ed49d126775285fd9ecd0b6a6717ad23..5f5e053a942ece384be0b04617e7ccec75a042e0 100644
--- a/src/fieldSources/basicSource/basicSource/basicSourceList.H
+++ b/src/fieldSources/basicSource/basicSourceList.H
@@ -157,7 +157,7 @@ public:
             //- Read dictionary
             virtual bool read(const dictionary& dict);
 
-            //- Write data to Istream
+            //- Write data to Ostream
             virtual bool writeData(Ostream& os) const;
 
             //- Ostream operator
diff --git a/src/fieldSources/basicSource/basicSource/basicSourceListTemplates.C b/src/fieldSources/basicSource/basicSourceListTemplates.C
similarity index 100%
rename from src/fieldSources/basicSource/basicSource/basicSourceListTemplates.C
rename to src/fieldSources/basicSource/basicSourceListTemplates.C
diff --git a/src/fieldSources/basicSource/basicSource/makeBasicSource.H b/src/fieldSources/basicSource/makeBasicSource.H
similarity index 100%
rename from src/fieldSources/basicSource/basicSource/makeBasicSource.H
rename to src/fieldSources/basicSource/makeBasicSource.H
diff --git a/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.C b/src/fieldSources/derived/actuationDiskSource/actuationDiskSource.C
similarity index 100%
rename from src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.C
rename to src/fieldSources/derived/actuationDiskSource/actuationDiskSource.C
diff --git a/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.H b/src/fieldSources/derived/actuationDiskSource/actuationDiskSource.H
similarity index 100%
rename from src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.H
rename to src/fieldSources/derived/actuationDiskSource/actuationDiskSource.H
diff --git a/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C b/src/fieldSources/derived/actuationDiskSource/actuationDiskSourceTemplates.C
similarity index 100%
rename from src/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C
rename to src/fieldSources/derived/actuationDiskSource/actuationDiskSourceTemplates.C
diff --git a/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.C b/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.C
new file mode 100644
index 0000000000000000000000000000000000000000..bef68f3c4e927ddba13685181a863bd26f66c3bd
--- /dev/null
+++ b/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.C
@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "fixedTemperatureSource.H"
+#include "fvMesh.H"
+#include "fvMatrices.H"
+#include "basicThermo.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(fixedTemperatureSource, 0);
+    addToRunTimeSelectionTable
+    (
+        basicSource,
+        fixedTemperatureSource,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fixedTemperatureSource::fixedTemperatureSource
+(
+    const word& name,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    basicSource(name, modelType, dict, mesh),
+    T_(readScalar(coeffs_.lookup("temperature")))
+{
+    fieldNames_.setSize(1, "energy");
+    applied_.setSize(1, false);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::fixedTemperatureSource::alwaysApply() const
+{
+    return true;
+}
+
+
+void Foam::fixedTemperatureSource::setValue
+(
+    fvMatrix<scalar>& eqn,
+    const label
+)
+{
+    const basicThermo& thermo =
+        mesh_.lookupObject<basicThermo>("thermophysicalProperties");
+
+    if (eqn.psi().name() == thermo.he().name())
+    {
+        const scalarField Tfield(cells_.size(), T_);
+        eqn.setValues(cells_, thermo.he(thermo.p(), Tfield, cells_));
+    }
+}
+
+
+void Foam::fixedTemperatureSource::writeData(Ostream& os) const
+{
+    os  << indent << name_ << endl;
+    dict_.write(os);
+}
+
+
+bool Foam::fixedTemperatureSource::read(const dictionary& dict)
+{
+    if (basicSource::read(dict))
+    {
+        coeffs_.readIfPresent("T", T_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.H b/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.H
new file mode 100644
index 0000000000000000000000000000000000000000..4281a38e2ccc5990db25be28b7b447c0e120889f
--- /dev/null
+++ b/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.H
@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::fixedTemperatureSource
+
+Description
+    Fixed temperature equation constraint
+
+    Sources described by:
+
+        fixedTemperatureSourceCoeffs
+        {
+            fieldNames      (h e hs);   // names of fields to apply source
+            temperature     500;        // fixed temperature [K]
+        }
+
+
+SourceFiles
+    basicSource.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedTemperatureSource_H
+#define fixedTemperatureSource_H
+
+#include "basicSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class fixedTemperatureSource Declaration
+\*---------------------------------------------------------------------------*/
+
+class fixedTemperatureSource
+:
+    public basicSource
+{
+
+protected:
+
+    // Protected data
+
+        //- Fixed temperature [K]
+        scalar T_;
+
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        fixedTemperatureSource(const fixedTemperatureSource&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const fixedTemperatureSource&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fixedTemperatureSource");
+
+
+    // Constructors
+
+        //- Construct from components
+        fixedTemperatureSource
+        (
+            const word& name,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~fixedTemperatureSource()
+    {}
+
+
+    // Member Functions
+
+        virtual bool alwaysApply() const;
+
+
+        // Set values directly
+
+            //- Scalar
+            virtual void setValue(fvMatrix<scalar>& eqn, const label fieldI);
+
+
+        // I-O
+
+            //- Write data
+            virtual void writeData(Ostream&) const;
+
+            //- Read dictionary
+            virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/fieldSources/basicSource/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C b/src/fieldSources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
similarity index 100%
rename from src/fieldSources/basicSource/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
rename to src/fieldSources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
diff --git a/src/fieldSources/basicSource/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H b/src/fieldSources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H
similarity index 100%
rename from src/fieldSources/basicSource/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H
rename to src/fieldSources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H
diff --git a/src/fieldSources/basicSource/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C b/src/fieldSources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C
similarity index 100%
rename from src/fieldSources/basicSource/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C
rename to src/fieldSources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C
diff --git a/src/fieldSources/basicSource/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H b/src/fieldSources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H
similarity index 100%
rename from src/fieldSources/basicSource/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H
rename to src/fieldSources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H
diff --git a/src/fieldSources/basicSource/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C b/src/fieldSources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C
similarity index 100%
rename from src/fieldSources/basicSource/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C
rename to src/fieldSources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C
diff --git a/src/fieldSources/basicSource/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H b/src/fieldSources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H
similarity index 100%
rename from src/fieldSources/basicSource/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H
rename to src/fieldSources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H
diff --git a/src/fieldSources/basicSource/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C b/src/fieldSources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C
similarity index 100%
rename from src/fieldSources/basicSource/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C
rename to src/fieldSources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C
diff --git a/src/fieldSources/basicSource/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H b/src/fieldSources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H
similarity index 100%
rename from src/fieldSources/basicSource/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H
rename to src/fieldSources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H
diff --git a/src/fieldSources/basicSource/pressureGradientExplicitSource/pressureGradientExplicitSource.C b/src/fieldSources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.C
similarity index 100%
rename from src/fieldSources/basicSource/pressureGradientExplicitSource/pressureGradientExplicitSource.C
rename to src/fieldSources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.C
diff --git a/src/fieldSources/basicSource/pressureGradientExplicitSource/pressureGradientExplicitSource.H b/src/fieldSources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.H
similarity index 100%
rename from src/fieldSources/basicSource/pressureGradientExplicitSource/pressureGradientExplicitSource.H
rename to src/fieldSources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.H
diff --git a/src/fieldSources/basicSource/pressureGradientExplicitSource/pressureGradientExplicitSourceIO.C b/src/fieldSources/derived/pressureGradientExplicitSource/pressureGradientExplicitSourceIO.C
similarity index 100%
rename from src/fieldSources/basicSource/pressureGradientExplicitSource/pressureGradientExplicitSourceIO.C
rename to src/fieldSources/derived/pressureGradientExplicitSource/pressureGradientExplicitSourceIO.C
diff --git a/src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.C b/src/fieldSources/derived/radialActuationDiskSource/radialActuationDiskSource.C
similarity index 100%
rename from src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.C
rename to src/fieldSources/derived/radialActuationDiskSource/radialActuationDiskSource.C
diff --git a/src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.H b/src/fieldSources/derived/radialActuationDiskSource/radialActuationDiskSource.H
similarity index 100%
rename from src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.H
rename to src/fieldSources/derived/radialActuationDiskSource/radialActuationDiskSource.H
diff --git a/src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C b/src/fieldSources/derived/radialActuationDiskSource/radialActuationDiskSourceTemplates.C
similarity index 100%
rename from src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C
rename to src/fieldSources/derived/radialActuationDiskSource/radialActuationDiskSourceTemplates.C
diff --git a/src/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C b/src/fieldSources/derived/rotorDiskSource/bladeModel/bladeModel.C
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C
rename to src/fieldSources/derived/rotorDiskSource/bladeModel/bladeModel.C
diff --git a/src/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.H b/src/fieldSources/derived/rotorDiskSource/bladeModel/bladeModel.H
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.H
rename to src/fieldSources/derived/rotorDiskSource/bladeModel/bladeModel.H
diff --git a/src/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C b/src/fieldSources/derived/rotorDiskSource/profileModel/lookup/lookupProfile.C
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C
rename to src/fieldSources/derived/rotorDiskSource/profileModel/lookup/lookupProfile.C
diff --git a/src/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.H b/src/fieldSources/derived/rotorDiskSource/profileModel/lookup/lookupProfile.H
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.H
rename to src/fieldSources/derived/rotorDiskSource/profileModel/lookup/lookupProfile.H
diff --git a/src/fieldSources/basicSource/rotorDiskSource/profileModel/profileModel.C b/src/fieldSources/derived/rotorDiskSource/profileModel/profileModel.C
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/profileModel/profileModel.C
rename to src/fieldSources/derived/rotorDiskSource/profileModel/profileModel.C
diff --git a/src/fieldSources/basicSource/rotorDiskSource/profileModel/profileModel.H b/src/fieldSources/derived/rotorDiskSource/profileModel/profileModel.H
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/profileModel/profileModel.H
rename to src/fieldSources/derived/rotorDiskSource/profileModel/profileModel.H
diff --git a/src/fieldSources/basicSource/rotorDiskSource/profileModel/profileModelList.C b/src/fieldSources/derived/rotorDiskSource/profileModel/profileModelList.C
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/profileModel/profileModelList.C
rename to src/fieldSources/derived/rotorDiskSource/profileModel/profileModelList.C
diff --git a/src/fieldSources/basicSource/rotorDiskSource/profileModel/profileModelList.H b/src/fieldSources/derived/rotorDiskSource/profileModel/profileModelList.H
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/profileModel/profileModelList.H
rename to src/fieldSources/derived/rotorDiskSource/profileModel/profileModelList.H
diff --git a/src/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.C b/src/fieldSources/derived/rotorDiskSource/profileModel/series/seriesProfile.C
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.C
rename to src/fieldSources/derived/rotorDiskSource/profileModel/series/seriesProfile.C
diff --git a/src/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.H b/src/fieldSources/derived/rotorDiskSource/profileModel/series/seriesProfile.H
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.H
rename to src/fieldSources/derived/rotorDiskSource/profileModel/series/seriesProfile.H
diff --git a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C b/src/fieldSources/derived/rotorDiskSource/rotorDiskSource.C
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C
rename to src/fieldSources/derived/rotorDiskSource/rotorDiskSource.C
diff --git a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H b/src/fieldSources/derived/rotorDiskSource/rotorDiskSource.H
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H
rename to src/fieldSources/derived/rotorDiskSource/rotorDiskSource.H
diff --git a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceI.H b/src/fieldSources/derived/rotorDiskSource/rotorDiskSourceI.H
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceI.H
rename to src/fieldSources/derived/rotorDiskSource/rotorDiskSourceI.H
diff --git a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceTemplates.C b/src/fieldSources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceTemplates.C
rename to src/fieldSources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
diff --git a/src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C b/src/fieldSources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.C
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C
rename to src/fieldSources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.C
diff --git a/src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.H b/src/fieldSources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.H
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.H
rename to src/fieldSources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.H
diff --git a/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C b/src/fieldSources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
rename to src/fieldSources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
diff --git a/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H b/src/fieldSources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H
rename to src/fieldSources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H
diff --git a/src/fieldSources/basicSource/rotorDiskSource/trimModel/trimModel/trimModel.C b/src/fieldSources/derived/rotorDiskSource/trimModel/trimModel/trimModel.C
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/trimModel/trimModel/trimModel.C
rename to src/fieldSources/derived/rotorDiskSource/trimModel/trimModel/trimModel.C
diff --git a/src/fieldSources/basicSource/rotorDiskSource/trimModel/trimModel/trimModel.H b/src/fieldSources/derived/rotorDiskSource/trimModel/trimModel/trimModel.H
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/trimModel/trimModel/trimModel.H
rename to src/fieldSources/derived/rotorDiskSource/trimModel/trimModel/trimModel.H
diff --git a/src/fieldSources/basicSource/rotorDiskSource/trimModel/trimModel/trimModelNew.C b/src/fieldSources/derived/rotorDiskSource/trimModel/trimModel/trimModelNew.C
similarity index 100%
rename from src/fieldSources/basicSource/rotorDiskSource/trimModel/trimModel/trimModelNew.C
rename to src/fieldSources/derived/rotorDiskSource/trimModel/trimModel/trimModelNew.C
diff --git a/src/fieldSources/basicSource/codedSource/CodedSource.C b/src/fieldSources/general/codedSource/CodedSource.C
similarity index 100%
rename from src/fieldSources/basicSource/codedSource/CodedSource.C
rename to src/fieldSources/general/codedSource/CodedSource.C
diff --git a/src/fieldSources/basicSource/codedSource/CodedSource.H b/src/fieldSources/general/codedSource/CodedSource.H
similarity index 100%
rename from src/fieldSources/basicSource/codedSource/CodedSource.H
rename to src/fieldSources/general/codedSource/CodedSource.H
diff --git a/src/fieldSources/basicSource/codedSource/CodedSourceIO.C b/src/fieldSources/general/codedSource/CodedSourceIO.C
similarity index 100%
rename from src/fieldSources/basicSource/codedSource/CodedSourceIO.C
rename to src/fieldSources/general/codedSource/CodedSourceIO.C
diff --git a/src/fieldSources/basicSource/codedSource/codedSource.C b/src/fieldSources/general/codedSource/codedSource.C
similarity index 100%
rename from src/fieldSources/basicSource/codedSource/codedSource.C
rename to src/fieldSources/general/codedSource/codedSource.C
diff --git a/src/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.C b/src/fieldSources/general/explicitSetValue/ExplicitSetValue.C
similarity index 100%
rename from src/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.C
rename to src/fieldSources/general/explicitSetValue/ExplicitSetValue.C
diff --git a/src/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.H b/src/fieldSources/general/explicitSetValue/ExplicitSetValue.H
similarity index 98%
rename from src/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.H
rename to src/fieldSources/general/explicitSetValue/ExplicitSetValue.H
index 1c3adcb279fab9deb8f133aa8e8d3b4a36cc516e..653029340ba088d69bac55f63f4301571adde02f 100644
--- a/src/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.H
+++ b/src/fieldSources/general/explicitSetValue/ExplicitSetValue.H
@@ -105,7 +105,7 @@ public:
 
         // Evaluation
 
-            //- Set value on vector field
+            //- Set value on field
             virtual void setValue(fvMatrix<Type>& eqn, const label fieldI);
 
 
diff --git a/src/fieldSources/basicSource/explicitSetValue/ExplicitSetValueIO.C b/src/fieldSources/general/explicitSetValue/ExplicitSetValueIO.C
similarity index 100%
rename from src/fieldSources/basicSource/explicitSetValue/ExplicitSetValueIO.C
rename to src/fieldSources/general/explicitSetValue/ExplicitSetValueIO.C
diff --git a/src/fieldSources/basicSource/explicitSetValue/explicitSetValue.C b/src/fieldSources/general/explicitSetValue/explicitSetValue.C
similarity index 100%
rename from src/fieldSources/basicSource/explicitSetValue/explicitSetValue.C
rename to src/fieldSources/general/explicitSetValue/explicitSetValue.C
diff --git a/src/fieldSources/basicSource/explicitSource/ExplicitSource.C b/src/fieldSources/general/explicitSource/ExplicitSource.C
similarity index 100%
rename from src/fieldSources/basicSource/explicitSource/ExplicitSource.C
rename to src/fieldSources/general/explicitSource/ExplicitSource.C
diff --git a/src/fieldSources/basicSource/explicitSource/ExplicitSource.H b/src/fieldSources/general/explicitSource/ExplicitSource.H
similarity index 100%
rename from src/fieldSources/basicSource/explicitSource/ExplicitSource.H
rename to src/fieldSources/general/explicitSource/ExplicitSource.H
diff --git a/src/fieldSources/basicSource/explicitSource/ExplicitSourceI.H b/src/fieldSources/general/explicitSource/ExplicitSourceI.H
similarity index 100%
rename from src/fieldSources/basicSource/explicitSource/ExplicitSourceI.H
rename to src/fieldSources/general/explicitSource/ExplicitSourceI.H
diff --git a/src/fieldSources/basicSource/explicitSource/ExplicitSourceIO.C b/src/fieldSources/general/explicitSource/ExplicitSourceIO.C
similarity index 100%
rename from src/fieldSources/basicSource/explicitSource/ExplicitSourceIO.C
rename to src/fieldSources/general/explicitSource/ExplicitSourceIO.C
diff --git a/src/fieldSources/basicSource/explicitSource/explicitSource.C b/src/fieldSources/general/explicitSource/explicitSource.C
similarity index 100%
rename from src/fieldSources/basicSource/explicitSource/explicitSource.C
rename to src/fieldSources/general/explicitSource/explicitSource.C
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 58c0b576737cd27884fa4063967ebeab8ab871d8..8f83ea99b95550838f6561e848c7319f932b31bd 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -364,9 +364,13 @@ $(solutionControl)/solutionControl/solutionControl.C
 $(solutionControl)/simpleControl/simpleControl.C
 $(solutionControl)/pimpleControl/pimpleControl.C
 
-porousMedia = $(general)/porousMedia
-$(porousMedia)/porousZone.C
-$(porousMedia)/porousZones.C
+porosity = $(general)/porosityModel
+$(porosity)/porosityModel/porosityModel.C
+$(porosity)/porosityModel/porosityModelNew.C
+$(porosity)/porosityModel/porosityModelList.C
+$(porosity)/porosityModel/IOporosityModelList.C
+$(porosity)/DarcyForchheimer/DarcyForchheimer.C
+$(porosity)/powerLaw/powerLaw.C
 
 MRF = $(general)/MRF
 $(MRF)/MRFZone.C
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C
new file mode 100644
index 0000000000000000000000000000000000000000..f292cb63c85741f2f12496502b9baff933ffa945
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C
@@ -0,0 +1,202 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "addToRunTimeSelectionTable.H"
+#include "DarcyForchheimer.H"
+#include "geometricOneField.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace porosityModels
+    {
+        defineTypeNameAndDebug(DarcyForchheimer, 0);
+        addToRunTimeSelectionTable(porosityModel, DarcyForchheimer, mesh);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::porosityModels::DarcyForchheimer::DarcyForchheimer
+(
+    const word& name,
+    const word& modelType,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    porosityModel(name, modelType, mesh, dict),
+    coordSys_(coeffs_, mesh),
+    D_("D", dimless/sqr(dimLength), tensor::zero),
+    F_("F", dimless/dimLength, tensor::zero),
+    rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho")),
+    muName_(coeffs_.lookupOrDefault<word>("mu", "mu")),
+    nuName_(coeffs_.lookupOrDefault<word>("nu", "nu"))
+{
+    // local-to-global transformation tensor
+    const tensor& E = coordSys_.R();
+
+    dimensionedVector d(coeffs_.lookup("d"));
+    if (D_.dimensions() != d.dimensions())
+    {
+        FatalIOErrorIn
+        (
+            "Foam::porosityModels::DarcyForchheimer::DarcyForchheimer"
+            "("
+                "const dictionary&, "
+                "const coordinateSystem&, "
+                "const keyType&"
+            ")",
+            coeffs_
+        )   << "incorrect dimensions for d: " << d.dimensions()
+            << " should be " << D_.dimensions()
+            << exit(FatalIOError);
+    }
+
+    adjustNegativeResistance(d);
+
+    D_.value().xx() = d.value().x();
+    D_.value().yy() = d.value().y();
+    D_.value().zz() = d.value().z();
+    D_.value() = (E & D_ & E.T()).value();
+
+    dimensionedVector f(coeffs_.lookup("f"));
+    if (F_.dimensions() != f.dimensions())
+    {
+        FatalIOErrorIn
+        (
+            "Foam::porosityModels::DarcyForchheimer::DarcyForchheimer"
+            "("
+                "const dictionary&, "
+                "const coordinateSystem&, "
+                "const keyType&"
+            ")",
+            coeffs_
+        )   << "incorrect dimensions for f: " << f.dimensions()
+            << " should be " << F_.dimensions()
+            << exit(FatalIOError);
+    }
+
+    adjustNegativeResistance(f);
+
+    // leading 0.5 is from 1/2*rho
+    F_.value().xx() = 0.5*f.value().x();
+    F_.value().yy() = 0.5*f.value().y();
+    F_.value().zz() = 0.5*f.value().z();
+    F_.value() = (E & F_ & E.T()).value();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::porosityModels::DarcyForchheimer::~DarcyForchheimer()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::porosityModels::DarcyForchheimer::correct
+(
+    fvVectorMatrix& UEqn
+) const
+{
+    const vectorField& U = UEqn.psi();
+    const scalarField& V = mesh_.V();
+    scalarField& Udiag = UEqn.diag();
+    vectorField& Usource = UEqn.source();
+ 
+    if (UEqn.dimensions() == dimForce)
+    {
+        const volScalarField& rho =
+            mesh_.lookupObject<volScalarField>(rhoName_);
+        const volScalarField& mu =
+            mesh_.lookupObject<volScalarField>(muName_);
+
+        apply(Udiag, Usource, V, rho, mu, U);
+    }
+    else
+    {
+        const volScalarField& nu =
+            mesh_.lookupObject<volScalarField>(nuName_);
+
+        apply(Udiag, Usource, V, geometricOneField(), nu, U);
+    }
+}
+
+
+void Foam::porosityModels::DarcyForchheimer::correct
+(
+    fvVectorMatrix& UEqn,
+    const volScalarField& rho,
+    const volScalarField& mu
+) const
+{
+    const vectorField& U = UEqn.psi();
+    const scalarField& V = mesh_.V();
+    scalarField& Udiag = UEqn.diag();
+    vectorField& Usource = UEqn.source();
+ 
+    apply(Udiag, Usource, V, rho, mu, U);
+}
+
+
+void Foam::porosityModels::DarcyForchheimer::correct
+(
+    const fvVectorMatrix& UEqn,
+    volTensorField& AU
+) const
+{
+    const vectorField& U = UEqn.psi();
+
+    if (UEqn.dimensions() == dimForce)
+    {
+        const volScalarField& rho =
+            mesh_.lookupObject<volScalarField>(rhoName_);
+        const volScalarField& mu =
+            mesh_.lookupObject<volScalarField>(muName_);
+
+        apply(AU, rho, mu, U);
+    }
+    else
+    {
+        const volScalarField& nu =
+            mesh_.lookupObject<volScalarField>(nuName_);
+
+        apply(AU, geometricOneField(), nu, U);
+    }
+}
+
+
+void Foam::porosityModels::DarcyForchheimer::writeData(Ostream& os) const
+{
+    os  << indent << name_ << endl;
+    dict_.write(os);
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H
new file mode 100644
index 0000000000000000000000000000000000000000..928422938a1da0357a1563335e93f20b618409f4
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H
@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::DarcyForchheimer
+
+Description
+    Darcy-Forchheimer law porosity model, given by:
+
+        \f[
+            S = - (\mu d + \frac{\rho |U|}{2} f) U
+        \f]
+
+    where
+    \vartable
+        d        | Darcy coefficient [1/m2]
+        f        | Forchheimer coefficient [1/m]
+    \endvartable
+
+    Since negative Darcy/Forchheimer parameters are invalid, they can be used
+    to specify a multiplier (of the max component).
+
+    The orientation of the porous region is defined with the same notation as
+    a co-ordinate system, but only a Cartesian co-ordinate system is valid.
+
+SourceFiles
+    DarcyForchheimer.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DarcyForchheimer_H
+#define DarcyForchheimer_H
+
+#include "porosityModel.H"
+#include "coordinateSystem.H"
+#include "dimensionedTensor.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace porosityModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class DarcyForchheimer Declaration
+\*---------------------------------------------------------------------------*/
+
+class DarcyForchheimer
+:
+    public porosityModel
+{
+private:
+
+    // Private data
+
+        //- Local co-ordinate system
+        coordinateSystem coordSys_;
+
+        //- Darcy coefficient [1/m2]
+        dimensionedTensor D_;
+
+        //- Forchheimer coefficient [1/m]
+        dimensionedTensor F_;
+
+        //- Name of density field
+        word rhoName_;
+
+        //- Name of dynamic viscosity field
+        word muName_;
+
+        //- Name of kinematic viscosity field
+        word nuName_;
+
+
+    // Private Member Functions
+
+        //- Apply
+        template<class RhoFieldType>
+        void apply
+        (
+            scalarField& Udiag,
+            vectorField& Usource,
+            const scalarField& V,
+            const RhoFieldType& rho,
+            const scalarField& mu,
+            const vectorField& U
+        ) const;
+
+        //- Apply
+        template<class RhoFieldType>
+        void apply
+        (
+            tensorField& AU,
+            const RhoFieldType& rho,
+            const scalarField& mu,
+            const vectorField& U
+        ) const;
+
+        //- Disallow default bitwise copy construct
+        DarcyForchheimer(const DarcyForchheimer&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const DarcyForchheimer&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("DarcyForchheimer");
+
+    //- Constructor
+    DarcyForchheimer
+    (
+        const word& name,
+        const word& modelType,
+        const fvMesh& mesh,
+        const dictionary& dict
+    );
+
+    //- Destructor
+    virtual ~DarcyForchheimer();
+
+
+    // Member Functions
+
+        //- Add resistance
+        virtual void correct(fvVectorMatrix& UEqn) const;
+
+        //- Add resistance
+        virtual void correct
+        (
+            fvVectorMatrix& UEqn,
+            const volScalarField& rho,
+            const volScalarField& mu
+        ) const;
+
+        //- Add resistance
+        virtual void correct
+        (
+            const fvVectorMatrix& UEqn,
+            volTensorField& AU            
+        ) const;
+
+
+    // I-O
+
+        //- Write
+        void writeData(Ostream& os) const;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace porosityModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "DarcyForchheimerTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..184206a010dccd166878507f0b002e1e8d35f9b9
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C
@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class RhoFieldType>
+void Foam::porosityModels::DarcyForchheimer::apply
+(
+    scalarField& Udiag,
+    vectorField& Usource,
+    const scalarField& V,
+    const RhoFieldType& rho,
+    const scalarField& mu,
+    const vectorField& U
+) const
+{
+    const tensor& D = D_.value();
+    const tensor& F = F_.value();
+
+    forAll(cellZoneIds_, zoneI)
+    {
+        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+
+        forAll(cells, i)
+        {
+            const label cellI = cells[i];
+
+            const tensor Cd = mu[cellI]*D + (rho[cellI]*mag(U[cellI]))*F;
+
+            const scalar isoCd = tr(Cd);
+
+            Udiag[cellI] += V[cellI]*isoCd;
+            Usource[cellI] -= V[cellI]*((Cd - I*isoCd) & U[cellI]);
+        }
+    }
+}
+
+
+template<class RhoFieldType>
+void Foam::porosityModels::DarcyForchheimer::apply
+(
+    tensorField& AU,
+    const RhoFieldType& rho,
+    const scalarField& mu,
+    const vectorField& U
+) const
+{
+    const tensor& D = D_.value();
+    const tensor& F = F_.value();
+
+    forAll(cellZoneIds_, zoneI)
+    {
+        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+
+        forAll(cells, i)
+        {
+            const label cellI = cells[i];
+            AU[cellI] += mu[cellI]*D + (rho[cellI]*mag(U[cellI]))*F;
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/IOporosityModelList.C
similarity index 55%
rename from src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.C
rename to src/finiteVolume/cfdTools/general/porosityModel/porosityModel/IOporosityModelList.C
index 8900a7350a2e1417b0eea6f9b7efffdfa1604709..2a655eefddec9d7fedb8efc01ac4ba1f16c986d9 100644
--- a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.C
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/IOporosityModelList.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,41 +23,68 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "thermalPorousZones.H"
-#include "volFields.H"
+#include "IOporosityModelList.H"
+#include "fvMesh.H"
+#include "Time.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-namespace Foam
+Foam::IOobject Foam::IOporosityModelList::createIOobject
+(
+    const fvMesh& mesh
+) const
 {
-    defineTemplateTypeNameAndDebug(IOPtrList<thermalPorousZone>, 0);
+    IOobject io
+    (
+        "porosityProperties",
+        mesh.time().constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (io.headerOk())
+    {
+        Info<< "Creating porosity model list from " << io.name() << nl << endl;
+
+        io.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
+        return io;
+    }
+    else
+    {
+        Info<< "No porosity models present" << nl << endl;
+
+        io.readOpt() = IOobject::NO_READ;
+        return io;
+    }
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::thermalPorousZones::thermalPorousZones
+Foam::IOporosityModelList::IOporosityModelList
 (
     const fvMesh& mesh
 )
 :
-    PorousZones<thermalPorousZone>(mesh)
+    IOdictionary(createIOobject(mesh)),
+    porosityModelList(mesh, *this)
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::thermalPorousZones::addEnergySource
-(
-    const fluidThermo& thermo,
-    const volScalarField& rho,
-    fvScalarMatrix& hEqn
-) const
+bool Foam::IOporosityModelList::read()
 {
-    forAll(*this, i)
+    if (regIOobject::read())
+    {
+        porosityModelList::read(*this);
+        return true;
+    }
+    else
     {
-        operator[](i).addEnergySource(thermo, rho, hEqn);
+        return false;
     }
 }
 
 
 // ************************************************************************* //
+
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/noThermalModel/noThermalModel.H b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/IOporosityModelList.H
similarity index 62%
rename from src/thermophysicalModels/thermalPorousZone/thermalModel/noThermalModel/noThermalModel.H
rename to src/finiteVolume/cfdTools/general/porosityModel/porosityModel/IOporosityModelList.H
index 1b9c1d8b83242a182b99fac406642d38c24f43fe..e9e534926ab338adf3107a0a13bfab93f3ec368c 100644
--- a/src/thermophysicalModels/thermalPorousZone/thermalModel/noThermalModel/noThermalModel.H
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/IOporosityModelList.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,66 +22,68 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::porousMedia::noThermalModel
+    Foam::IOporosityModelList
 
 Description
-    A dummy thermal model porous media, corresponding to the 'none' option
+    List of porosity models with IO functionality
+
+SourceFiles
+    IOporosityModelList.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef noThermalModel_H
-#define noThermalModel_H
+#ifndef IOporosityModelList_H
+#define IOporosityModelList_H
 
-#include "thermalModel.H"
-#include "runTimeSelectionTables.H"
+#include "IOdictionary.H"
+#include "porosityModelList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-namespace porousMedia
-{
 
 /*---------------------------------------------------------------------------*\
-                       Class noThermalModel Declaration
+                     Class IOporosityModelList Declaration
 \*---------------------------------------------------------------------------*/
 
-class noThermalModel
+class IOporosityModelList
 :
-    public thermalModel
+    public IOdictionary,
+    public porosityModelList
 {
+private:
 
-public:
+    // Private Member Functions
 
-    //- Runtime type information
-    TypeName("none");
+        //- Create IO object if dictionary is present
+        IOobject createIOobject(const fvMesh& mesh) const;
 
+        //- Disallow default bitwise copy construct
+        IOporosityModelList(const IOporosityModelList&);
 
-    // Constructors
+        //- Disallow default bitwise assignment
+        void operator=(const IOporosityModelList&);
 
-        //- Construct from porous zone
-        noThermalModel(const porousZone&);
 
+public:
 
-    //- Destructor
-    virtual ~noThermalModel();
+    // Constructors
 
+        //- Construct from mesh
+        IOporosityModelList(const fvMesh& mesh);
 
-    // Member Functions
 
-        //- Add the thermal source to the enthalpy equation
-        virtual void addEnergySource
-        (
-            const fluidThermo&,
-            const volScalarField& rho,
-            fvScalarMatrix& hEqn
-        ) const;
-};
+        //- Destructor
+        virtual ~IOporosityModelList()
+        {}
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // Member Functions
 
-} // End namespace porousMedia
+        //- Read dictionary
+        virtual bool read();
+};
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
new file mode 100644
index 0000000000000000000000000000000000000000..f0997c822211a116787e1990364ad82955b952f0
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
@@ -0,0 +1,182 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "porosityModel.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(porosityModel, 0);
+    defineRunTimeSelectionTable(porosityModel, mesh);
+}
+
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+void Foam::porosityModel::adjustNegativeResistance(dimensionedVector& resist)
+{
+    scalar maxCmpt = max(0, cmptMax(resist.value()));
+
+    if (maxCmpt < 0)
+    {
+        FatalErrorIn
+        (
+            "void Foam::porosityModel::adjustNegativeResistance"
+            "("
+                "dimensionedVector&"
+            ")"
+        )   << "Negative resistances are invalid, resistance = " << resist
+            << exit(FatalError);
+    }
+    else
+    {
+        vector& val = resist.value();
+        for (label cmpt = 0; cmpt < vector::nComponents; cmpt++)
+        {
+            if (val[cmpt] < 0)
+            {
+                val[cmpt] *= -maxCmpt;
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::porosityModel::porosityModel
+(
+    const word& name,
+    const word& modelType,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    name_(name),
+    mesh_(mesh),
+    dict_(dict),
+    coeffs_(dict.subDict(modelType + "Coeffs")),
+    active_(readBool(dict_.lookup("active"))),
+    zoneName_(dict_.lookup("cellZone")),
+    cellZoneIds_(mesh_.cellZones().findIndices(zoneName_))
+{
+    Info<< "    creating porous zone: " << zoneName_ << endl;
+
+    bool foundZone = !cellZoneIds_.empty();
+    reduce(foundZone, orOp<bool>());
+
+    if (!foundZone && Pstream::master())
+    {
+        FatalErrorIn
+        (
+            "Foam::porosityModel::porosityModel"
+            "("
+                "const word&, "
+                "const word&, "
+                "const fvMesh&, "
+                "const dictionary&"
+            ")"
+        )   << "cannot find porous cellZone " << zoneName_
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::porosityModel::~porosityModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::porosityModel::porosityModel::addResistance
+(
+    fvVectorMatrix& UEqn
+) const
+{
+    if (cellZoneIds_.empty())
+    {
+        return;
+    }
+
+    this->correct(UEqn);
+}
+
+
+void Foam::porosityModel::porosityModel::addResistance
+(
+    fvVectorMatrix& UEqn,
+    const volScalarField& rho,
+    const volScalarField& mu
+) const
+{
+    if (cellZoneIds_.empty())
+    {
+        return;
+    }
+
+    this->correct(UEqn, rho, mu);
+}
+
+
+void Foam::porosityModel::porosityModel::addResistance
+(
+    const fvVectorMatrix& UEqn,
+    volTensorField& AU,
+    bool correctAUprocBC         
+) const
+{
+    if (cellZoneIds_.empty())
+    {
+        return;
+    }
+
+    this->correct(UEqn, AU);
+
+    if (correctAUprocBC)
+    {
+        // Correct the boundary conditions of the tensorial diagonal to ensure
+        // processor boundaries are correctly handled when AU^-1 is interpolated
+        // for the pressure equation.
+        AU.correctBoundaryConditions();
+    }
+}
+
+
+bool Foam::porosityModel::read(const dictionary& dict)
+{
+    active_ = readBool(dict.lookup("active"));
+    coeffs_ = dict.subDict(type() + "Coeffs");
+    dict.lookup("cellZone") >> zoneName_;
+    cellZoneIds_ = mesh_.cellZones().findIndices(zoneName_);
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..56131ba3dd97cd5c287f699acc9eb3f9fcce153c
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H
@@ -0,0 +1,241 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::porosityModel
+
+Description
+    Top level model for porosity models
+
+SourceFiles
+    porosityModel.C
+    porosityModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef porosityModel_H
+#define porosityModel_H
+
+#include "fvMesh.H"
+#include "dictionary.H"
+#include "fvMatricesFwd.H"
+#include "runTimeSelectionTables.H"
+#include "dimensionedVector.H"
+#include "keyType.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class porosityModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class porosityModel
+{
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        porosityModel(const porosityModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const porosityModel&);
+
+
+protected:
+
+    // Protected data
+
+        //- Porosity name
+        word name_;
+
+        //- Reference to the mesh database
+        const fvMesh& mesh_;
+
+        //- Dictionary used for model construction
+        const dictionary dict_;
+
+        //- Model coefficients dictionary
+        dictionary coeffs_;
+
+        //- Porosity active flag
+        bool active_;
+
+        //- Name(s) of cell-zone
+        keyType zoneName_;
+
+        //- Cell zone Ids
+        labelList cellZoneIds_;
+
+
+    // Protected Member Functions
+
+        //- Adjust negative resistance values to be multiplier of max value
+        void adjustNegativeResistance(dimensionedVector& resist);
+
+        virtual void correct(fvVectorMatrix& UEqn) const = 0;
+
+        virtual void correct
+        (
+            fvVectorMatrix& UEqn,
+            const volScalarField& rho,
+            const volScalarField& mu
+        ) const = 0;
+
+        virtual void correct
+        (
+            const fvVectorMatrix& UEqn,
+            volTensorField& AU
+        ) const = 0;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("porosityModel");
+
+    //- Selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        porosityModel,
+        mesh,
+        (
+            const word& modelName,
+            const word& name,
+            const fvMesh& mesh,
+            const dictionary& dict
+        ),
+        (modelName, name, mesh, dict)
+    );
+
+    //- Constructor
+    porosityModel
+    (
+        const word& name,
+        const word& modelType,
+        const fvMesh& mesh,
+        const dictionary& dict
+    );
+
+    //- Return pointer to new porosityModel object created on the freestore
+    //  from an Istream
+    class iNew
+    {
+        //- Reference to the mesh database
+        const fvMesh& mesh_;
+        const word& name_;
+
+    public:
+
+        iNew
+        (
+            const fvMesh& mesh,
+            const word& name
+        )
+        :
+            mesh_(mesh),
+            name_(name)
+        {}
+
+        autoPtr<porosityModel> operator()(Istream& is) const
+        {
+            const dictionary dict(is);
+
+            return autoPtr<porosityModel>
+            (
+                porosityModel::New
+                (
+                    name_,
+                    mesh_,
+                    dict
+                )
+            );
+        }
+    };
+
+    //- Selector
+    static autoPtr<porosityModel> New
+    (
+        const word& name,
+        const fvMesh& mesh,
+        const dictionary& dict
+    );
+
+    //- Destructor
+    virtual ~porosityModel();
+
+
+    // Member Functions
+
+        //- Return const access to the porosity model name
+        inline const word& name() const;
+
+        //- Return const access to the porosity active flag
+        inline bool active() const;
+
+        //- Add resistance
+        virtual void addResistance(fvVectorMatrix& UEqn) const;
+
+        //- Add resistance
+        virtual void addResistance
+        (
+            fvVectorMatrix& UEqn,
+            const volScalarField& rho,
+            const volScalarField& mu
+        ) const;
+
+        //- Add resistance
+        virtual void addResistance
+        (
+            const fvVectorMatrix& UEqn,
+            volTensorField& AU,
+            bool correctAUprocBC
+        ) const;
+
+
+    // I-O
+
+        //- Write
+        virtual void writeData(Ostream& os) const = 0;
+
+        //- Read porosity dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "porosityModelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelI.H
similarity index 82%
rename from src/finiteVolume/cfdTools/general/porousMedia/porousZones.C
rename to src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelI.H
index 5130ca90e98e68758bb048d805557b5dcd9284c9..1df2d1197298852ceb8523be29a096f87085b6c5 100644
--- a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.C
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,13 +23,16 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "porousZones.H"
+inline const Foam::word& Foam::porosityModel::name() const
+{
+    return name_;
+}
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-namespace Foam
+inline bool Foam::porosityModel::active() const
 {
-    defineTemplateTypeNameAndDebug(IOPtrList<porousZone>, 0);
+    return active_;
 }
 
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.C
new file mode 100644
index 0000000000000000000000000000000000000000..14ab39d96436521ad20c9db8cdeef36972d4fc86
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.C
@@ -0,0 +1,183 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "porosityModelList.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+/*
+void Foam::porosityModelList::XXX()
+{}
+*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::porosityModelList::porosityModelList
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    PtrList<porosityModel>(),
+    mesh_(mesh)
+{
+    reset(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::porosityModelList::~porosityModelList()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::porosityModelList::active() const
+{
+    bool a = false;
+    forAll(*this, i)
+    {
+        a = a || this->operator[](i).active();
+    }
+
+    if (!a)
+    {
+        Info<< "No porosity models active" << endl;
+    }
+
+    return a;
+}
+
+
+void Foam::porosityModelList::reset(const dictionary& dict)
+{
+    label count = 0;
+    forAllConstIter(dictionary, dict, iter)
+    {
+        if (iter().isDict())
+        {
+            count++;
+        }
+    }
+
+    this->setSize(count);
+    label i = 0;
+    forAllConstIter(dictionary, dict, iter)
+    {
+        if (iter().isDict())
+        {
+            const word& name = iter().keyword();
+            const dictionary& modelDict = iter().dict();
+
+            this->set
+            (
+                i++,
+                porosityModel::New(name, mesh_, modelDict)
+            );
+        }
+    }
+}
+
+
+bool Foam::porosityModelList::read(const dictionary& dict)
+{
+    bool allOk = true;
+    forAll(*this, i)
+    {
+        porosityModel& pm = this->operator[](i);
+        bool ok = pm.read(dict.subDict(pm.name()));
+        allOk = (allOk && ok);
+    }
+    return allOk;
+}
+
+
+bool Foam::porosityModelList::writeData(Ostream& os) const
+{
+    forAll(*this, i)
+    {
+        os  << nl;
+        this->operator[](i).writeData(os);
+    }
+
+    return os.good();
+}
+
+
+void Foam::porosityModelList::addResistance
+(
+    fvVectorMatrix& UEqn
+) const
+{
+    forAll(*this, i)
+    {
+        this->operator[](i).addResistance(UEqn);
+    }
+}
+
+
+void Foam::porosityModelList::addResistance
+(
+    fvVectorMatrix& UEqn,
+    const volScalarField& rho,
+    const volScalarField& mu
+) const
+{
+    forAll(*this, i)
+    {
+        this->operator[](i).addResistance(UEqn, rho, mu);
+    }
+}
+
+
+void Foam::porosityModelList::addResistance
+(
+    const fvVectorMatrix& UEqn,
+    volTensorField& AU,
+    bool correctAUprocBC         
+) const
+{
+    forAll(*this, i)
+    {
+        this->operator[](i).addResistance(UEqn, AU, correctAUprocBC);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const porosityModelList& models
+)
+{
+    models.writeData(os);
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.H b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.H
new file mode 100644
index 0000000000000000000000000000000000000000..e580a2d0fadd7c46a0df99d1acc83f4c00e203a0
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.H
@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::porosityModelList
+
+Description
+    List container for porosity models
+
+SourceFiles
+    porosityModelList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef porosityModelList_H
+#define porosityModelList_H
+
+#include "fvMesh.H"
+#include "dictionary.H"
+#include "fvMatricesFwd.H"
+#include "porosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+class porosityModelList;
+Ostream& operator<<(Ostream& os, const porosityModelList& models);
+
+/*---------------------------------------------------------------------------*\
+                      Class porosityModelList Declaration
+\*---------------------------------------------------------------------------*/
+
+class porosityModelList
+:
+    PtrList<porosityModel>
+{
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        porosityModelList(const porosityModelList&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const porosityModelList&);
+
+
+protected:
+
+    // Protected data
+
+        //- Reference to the mesh database
+        const fvMesh& mesh_;
+
+
+public:
+
+    //- Constructor
+    porosityModelList(const fvMesh& mesh, const dictionary& dict);
+
+    //- Destructor
+    ~porosityModelList();
+
+
+    // Member Functions
+
+        //- Return active status
+        bool active() const;
+
+        //- Reset the source list
+        void reset(const dictionary& dict);
+
+        //- Add resistance
+        void addResistance(fvVectorMatrix& UEqn) const;
+
+        //- Add resistance
+        void addResistance
+        (
+            fvVectorMatrix& UEqn,
+            const volScalarField& rho,
+            const volScalarField& mu
+        ) const;
+
+        //- Add resistance
+        void addResistance
+        (
+            const fvVectorMatrix& UEqn,
+            volTensorField& AU,
+            bool correctAUprocBC = true
+        ) const;
+
+
+        // I-O
+
+            //- Read dictionary
+            bool read(const dictionary& dict);
+
+            //- Write data to Ostream
+            bool writeData(Ostream& os) const;
+
+            //- Ostream operator
+            friend Ostream& operator<<
+            (
+                Ostream& os,
+                const porosityModelList& models
+            );
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModelNew.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelNew.C
similarity index 56%
rename from src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModelNew.C
rename to src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelNew.C
index 9231ebe4edafb9b4ebca7055ddbdf013c7a6eb4e..1da8c91c4dfb583e022ce27059a0abaee880e63e 100644
--- a/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModelNew.C
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelNew.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,42 +23,43 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "thermalModel.H"
+#include "porosityModel.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::autoPtr<Foam::porousMedia::thermalModel>
-Foam::porousMedia::thermalModel::New
+Foam::autoPtr<Foam::porosityModel> Foam::porosityModel::New
 (
-    const porousZone& pZone
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
 )
 {
-    // a missing thermalModel is the same as type "none"
-    word modelType("none");
+    const word modelType(dict.lookup("type"));
 
-    if (const dictionary* dictPtr = pZone.dict().subDictPtr("thermalModel"))
-    {
-        dictPtr->lookup("type") >> modelType;
-    }
-
-    Info<< "Selecting thermalModel " << modelType << endl;
+    Info<< "Porosity region " << name << ":" << nl
+        << "    selecting model: " << modelType << endl;
 
-    pZoneConstructorTable::iterator cstrIter =
-        pZoneConstructorTablePtr_->find(modelType);
+    meshConstructorTable::iterator cstrIter =
+        meshConstructorTablePtr_->find(modelType);
 
-    if (cstrIter == pZoneConstructorTablePtr_->end())
+    if (cstrIter == meshConstructorTablePtr_->end())
     {
         FatalErrorIn
         (
-            "porousMedia::thermalModel::New(const porousZone&)"
-        )   << "Unknown thermalModel type "
-            << modelType << nl << nl
-            << "Valid thermalModel types are :" << endl
-            << pZoneConstructorTablePtr_->sortedToc()
-            << abort(FatalError);
+            "porosityModel::New"
+            "("
+                "const word& name,"
+                "const fvMesh&, "
+                "const dictionary&"
+            ")"
+        )
+            << "Unknown " << typeName << " type " << modelType << nl << nl
+            << "Valid " << typeName << " types are:" << nl
+            << meshConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
     }
 
-    return autoPtr<porousMedia::thermalModel>(cstrIter()(pZone));
+    return autoPtr<porosityModel>(cstrIter()(name, modelType, mesh, dict));
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.C b/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.C
new file mode 100644
index 0000000000000000000000000000000000000000..15fbae39dfb78486fe133002bb94ef1990956ec1
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.C
@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "addToRunTimeSelectionTable.H"
+#include "powerLaw.H"
+#include "geometricOneField.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace porosityModels
+    {
+        defineTypeNameAndDebug(powerLaw, 0);
+        addToRunTimeSelectionTable(porosityModel, powerLaw, mesh);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::porosityModels::powerLaw::powerLaw
+(
+    const word& name,
+    const word& modelType,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    porosityModel(name, modelType, mesh, dict),
+    C0_(readScalar(coeffs_.lookup("C0"))),
+    C1_(readScalar(coeffs_.lookup("C1"))),
+    rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::porosityModels::powerLaw::~powerLaw()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::porosityModels::powerLaw::correct
+(
+    fvVectorMatrix& UEqn
+) const
+{
+    const vectorField& U = UEqn.psi();
+    const scalarField& V = mesh_.V();
+    scalarField& Udiag = UEqn.diag();
+ 
+    if (UEqn.dimensions() == dimForce)
+    {
+        const volScalarField& rho =
+            mesh_.lookupObject<volScalarField>(rhoName_);
+
+        apply(Udiag, V, rho, U);
+    }
+    else
+    {
+        apply(Udiag, V, geometricOneField(), U);
+    }
+}
+
+
+void Foam::porosityModels::powerLaw::correct
+(
+    fvVectorMatrix& UEqn,
+    const volScalarField& rho,
+    const volScalarField& mu
+) const
+{
+    const vectorField& U = UEqn.psi();
+    const scalarField& V = mesh_.V();
+    scalarField& Udiag = UEqn.diag();
+ 
+    apply(Udiag, V, rho, U);
+}
+
+
+void Foam::porosityModels::powerLaw::correct
+(
+    const fvVectorMatrix& UEqn,
+    volTensorField& AU
+) const
+{
+    const vectorField& U = UEqn.psi();
+
+    if (UEqn.dimensions() == dimForce)
+    {
+        const volScalarField& rho =
+            mesh_.lookupObject<volScalarField>(rhoName_);
+
+        apply(AU, rho, U);
+    }
+    else
+    {
+        apply(AU, geometricOneField(), U);
+    }
+}
+
+
+void Foam::porosityModels::powerLaw::writeData(Ostream& os) const
+{
+    os  << indent << name_ << endl;
+    dict_.write(os);
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.H b/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.H
new file mode 100644
index 0000000000000000000000000000000000000000..070a00e9a021d4dc49d743f240f9e08b5edf03e5
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.H
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::powerLaw
+
+Description
+    Power law porosity model, given by:
+
+        \f[
+            S = - \rho C_0 |U|^{(C_1 - 1)} U
+        \f]
+
+    where
+    \vartable
+        C_0      | model coefficient
+        C_1      | model coefficient
+    \endvartable
+
+
+SourceFiles
+    powerLaw.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef powerLaw_H
+#define powerLaw_H
+
+#include "porosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace porosityModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class powerLaw Declaration
+\*---------------------------------------------------------------------------*/
+
+class powerLaw
+:
+    public porosityModel
+{
+private:
+
+    // Private data
+
+        //- C0 coefficient
+        scalar C0_;
+
+        //- C1 coefficient
+        scalar C1_;
+
+        //- Name of density field
+        word rhoName_;
+
+
+    // Private Member Functions
+
+        //- Apply resistance
+        template<class RhoFieldType>
+        void apply
+        (
+            scalarField& Udiag,
+            const scalarField& V,
+            const RhoFieldType& rho,
+            const vectorField& U
+        ) const;
+
+        //- Apply resistance
+        template<class RhoFieldType>
+        void apply
+        (
+            tensorField& AU,
+            const RhoFieldType& rho,
+            const vectorField& U
+        ) const;
+
+        //- Disallow default bitwise copy construct
+        powerLaw(const powerLaw&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const powerLaw&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("powerLaw");
+
+    //- Constructor
+    powerLaw
+    (
+        const word& name,
+        const word& modelType,
+        const fvMesh& mesh,
+        const dictionary& dict
+    );
+
+    //- Destructor
+    virtual ~powerLaw();
+
+
+    // Member Functions
+
+        //- Add resistance
+        virtual void correct(fvVectorMatrix& UEqn) const;
+
+        //- Add resistance
+        virtual void correct
+        (
+            fvVectorMatrix& UEqn,
+            const volScalarField& rho,
+            const volScalarField& mu
+        ) const;
+
+        //- Add resistance
+        virtual void correct
+        (
+            const fvVectorMatrix& UEqn,
+            volTensorField& AU            
+        ) const;
+
+
+    // I-O
+
+        //- Write
+        void writeData(Ostream& os) const;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace porosityModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "powerLawTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLawTemplates.C b/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLawTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..0e1d7d2343f9e558013e4fab619203e862ec0876
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLawTemplates.C
@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class RhoFieldType>
+void Foam::porosityModels::powerLaw::apply
+(
+    scalarField& Udiag,
+    const scalarField& V,
+    const RhoFieldType& rho,
+    const vectorField& U
+) const
+{
+    const scalar C0 = C0_;
+    const scalar C1m1b2 = (C1_ - 1.0)/2.0;
+
+    forAll(cellZoneIds_, zoneI)
+    {
+        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+
+        forAll(cells, i)
+        {
+            const label cellI = cells[i];
+
+            Udiag[cellI] +=
+                V[cellI]*rho[cellI]*C0*pow(magSqr(U[cellI]), C1m1b2);
+        }
+    }
+}
+
+
+template<class RhoFieldType>
+void Foam::porosityModels::powerLaw::apply
+(
+    tensorField& AU,
+    const RhoFieldType& rho,
+    const vectorField& U
+) const
+{
+    const scalar C0 = C0_;
+    const scalar C1m1b2 = (C1_ - 1.0)/2.0;
+
+    forAll(cellZoneIds_, zoneI)
+    {
+        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+
+        forAll(cells, i)
+        {
+            const label cellI = cells[i];
+
+            AU[cellI] =
+                AU[cellI] + I*(rho[cellI]*C0*pow(magSqr(U[cellI]), C1m1b2));
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.C b/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.C
deleted file mode 100644
index b219da6c00dc92188d7dbaec7b56e50531ef0cc1..0000000000000000000000000000000000000000
--- a/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.C
+++ /dev/null
@@ -1,223 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "PorousZones.H"
-#include "Time.H"
-#include "volFields.H"
-#include "fvm.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class ZoneType>
-template<class Type>
-void Foam::PorousZones<ZoneType>::modifyDdt(fvMatrix<Type>& m) const
-{
-    forAll(*this, i)
-    {
-        this->operator[](i).modifyDdt(m);
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class ZoneType>
-Foam::PorousZones<ZoneType>::PorousZones
-(
-    const fvMesh& mesh
-)
-:
-    IOPtrList<ZoneType>
-    (
-        IOobject
-        (
-            "porousZones",
-            mesh.time().constant(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
-        typename ZoneType::iNew(mesh)
-    ),
-    mesh_(mesh)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class ZoneType>
-template<class Type>
-Foam::tmp<Foam::fvMatrix<Type> >
-Foam::PorousZones<ZoneType>::ddt
-(
-    GeometricField<Type, fvPatchField, volMesh>& vf
-)
-{
-    tmp<fvMatrix<Type> > tres = fvm::ddt(vf);
-    modifyDdt(tres());
-    return tres;
-}
-
-
-template<class ZoneType>
-template<class Type>
-Foam::tmp<Foam::fvMatrix<Type> >
-Foam::PorousZones<ZoneType>::ddt
-(
-    const geometricOneField&,
-    GeometricField<Type, fvPatchField, volMesh>& vf
-)
-{
-    tmp<fvMatrix<Type> > tres = fvm::ddt(vf);
-    modifyDdt(tres());
-    return tres;
-}
-
-
-template<class ZoneType>
-template<class Type>
-Foam::tmp<Foam::fvMatrix<Type> >
-Foam::PorousZones<ZoneType>::ddt
-(
-    const dimensionedScalar& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
-)
-{
-    tmp<fvMatrix<Type> > tres = fvm::ddt(rho,vf);
-    modifyDdt(tres());
-    return tres;
-}
-
-
-template<class ZoneType>
-template<class Type>
-Foam::tmp<Foam::fvMatrix<Type> >
-Foam::PorousZones<ZoneType>::ddt
-(
-    const volScalarField& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
-)
-{
-    tmp<fvMatrix<Type> > tres = fvm::ddt(rho,vf);
-    modifyDdt(tres());
-    return tres;
-}
-
-template<class ZoneType>
-void Foam::PorousZones<ZoneType>::addResistance(fvVectorMatrix& UEqn) const
-{
-    forAll(*this, i)
-    {
-        this->operator[](i).addResistance(UEqn);
-    }
-}
-
-
-template<class ZoneType>
-void Foam::PorousZones<ZoneType>::addResistance
-(
-    fvVectorMatrix& UEqn,
-    const volScalarField& rho,
-    const volScalarField& mu
-) const
-{
-    forAll(*this, i)
-    {
-        this->operator[](i).addResistance(UEqn, rho, mu);
-    }
-}
-
-
-template<class ZoneType>
-void Foam::PorousZones<ZoneType>::addResistance
-(
-    const fvVectorMatrix& UEqn,
-    volTensorField& AU
-) const
-{
-    // addResistance for each zone, delaying the correction of the
-    // processor BCs of AU
-    forAll(*this, i)
-    {
-        this->operator[](i).addResistance(UEqn, AU, false);
-    }
-
-    // Correct the boundary conditions of the tensorial diagonal to ensure
-    // processor bounaries are correctly handled when AU^-1 is interpolated
-    // for the pressure equation.
-    AU.correctBoundaryConditions();
-}
-
-
-template<class ZoneType>
-bool Foam::PorousZones<ZoneType>::readData(Istream& is)
-{
-    this->clear();
-
-    IOPtrList<ZoneType> newLst
-    (
-        IOobject
-        (
-            "porousZones",
-            mesh_.time().constant(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false     // Don't re-register new zones with objectRegistry
-        ),
-        typename ZoneType::iNew(mesh_)
-    );
-
-    this->transfer(newLst);
-
-    return is.good();
-}
-
-
-template<class ZoneType>
-bool Foam::PorousZones<ZoneType>::writeData(Ostream& os, bool subDict) const
-{
-    // Write size of list
-    os << nl << this->size();
-
-    // Write beginning of contents
-    os << nl << token::BEGIN_LIST;
-
-    // Write list contents
-    forAll(*this, i)
-    {
-        os << nl;
-        this->operator[](i).writeDict(os, subDict);
-    }
-
-    // Write end of contents
-    os << token::END_LIST << nl;
-
-    // Check state of IOstream
-    return os.good();
-}
-
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.H b/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.H
deleted file mode 100644
index 2e1624c4e6f212b9ddaa91c372eb86ccc42dafc6..0000000000000000000000000000000000000000
--- a/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.H
+++ /dev/null
@@ -1,170 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::PorousZones
-
-Description
-    A centralized ZoneType collection.
-
-    Container class for a set of ZoneType with the ZoneType member
-    functions implemented to loop over the functions for each ZoneType.
-
-SourceFiles
-    PorousZones.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef PorousZones_H
-#define PorousZones_H
-
-#include "IOPtrList.H"
-
-#include "volFieldsFwd.H"
-#include "fvMatricesFwd.H"
-#include "dimensionedScalarFwd.H"
-#include "geometricOneField.H"
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of friend functions and operators
-class fvMesh;
-
-
-/*---------------------------------------------------------------------------*\
-                           Class PorousZones Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class ZoneType>
-class PorousZones
-:
-    public IOPtrList<ZoneType>
-{
-    // Private data
-
-        //- Reference to the finite volume mesh this zone is part of
-        const fvMesh& mesh_;
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        PorousZones(const PorousZones<ZoneType>&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const PorousZones<ZoneType>&);
-
-
-        //- modify time derivative elements
-        template<class Type>
-        void modifyDdt(fvMatrix<Type>&) const;
-
-public:
-
-    // Constructors
-
-        //- Construct from fvMesh
-        //  with automatically constructed coordinate systems list
-        PorousZones(const fvMesh&);
-
-
-    // Member Functions
-
-        //- mirror fvm::ddt with porosity
-        template<class Type>
-        tmp<fvMatrix<Type> > ddt
-        (
-            GeometricField<Type, fvPatchField, volMesh>&
-        );
-
-        //- mirror fvm::ddt with porosity
-        template<class Type>
-        tmp<fvMatrix<Type> > ddt
-        (
-            const geometricOneField&,
-            GeometricField<Type, fvPatchField, volMesh>&
-        );
-
-        //- mirror fvm::ddt with porosity
-        template<class Type>
-        tmp<fvMatrix<Type> > ddt
-        (
-            const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
-        );
-
-        //- mirror fvm::ddt with porosity
-        template<class Type>
-        tmp<fvMatrix<Type> > ddt
-        (
-            const volScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
-        );
-
-        //- Add the viscous and inertial resistance force contribution
-        //  to the momentum equation
-        void addResistance(fvVectorMatrix& UEqn) const;
-
-        //- Add the viscous and inertial resistance force contribution
-        //  to the momentum equation using rho and mu provided
-        void addResistance
-        (
-            fvVectorMatrix& UEqn,
-            const volScalarField& rho,
-            const volScalarField& mu
-        ) const;
-
-        //- Add the viscous and inertial resistance force contribution
-        //  to the tensorial diagonal
-        void addResistance
-        (
-            const fvVectorMatrix& UEqn,
-            volTensorField& AU
-        ) const;
-
-        //- read modified data
-        virtual bool readData(Istream&);
-
-        //- write data
-        bool writeData(Ostream&, bool subDict = true) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "PorousZones.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C
deleted file mode 100644
index edb0940cc7b419a1a82461a22a55c505c416f195..0000000000000000000000000000000000000000
--- a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C
+++ /dev/null
@@ -1,512 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*----------------------------------------------------------------------------*/
-
-#include "porousZone.H"
-#include "fvMesh.H"
-#include "fvMatrices.H"
-#include "geometricOneField.H"
-
-// * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
-
-// adjust negative resistance values to be multiplier of max value
-void Foam::porousZone::adjustNegativeResistance(dimensionedVector& resist)
-{
-    scalar maxCmpt = max(0, cmptMax(resist.value()));
-
-    if (maxCmpt < 0)
-    {
-        FatalErrorIn
-        (
-            "Foam::porousZone::porousZone::adjustNegativeResistance"
-            "(dimensionedVector&)"
-        )   << "negative resistances! " << resist
-            << exit(FatalError);
-    }
-    else
-    {
-        vector& val = resist.value();
-        for (label cmpt=0; cmpt < vector::nComponents; ++cmpt)
-        {
-            if (val[cmpt] < 0)
-            {
-                val[cmpt] *= -maxCmpt;
-            }
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::porousZone::porousZone
-(
-    const keyType& key,
-    const fvMesh& mesh,
-    const dictionary& dict
-)
-:
-    key_(key),
-    mesh_(mesh),
-    dict_(dict),
-    cellZoneIds_(mesh_.cellZones().findIndices(key_)),
-    coordSys_(dict, mesh),
-    porosity_(1),
-    intensity_(0),
-    mixingLength_(0),
-    C0_(0),
-    C1_(0),
-    D_("D", dimensionSet(0, -2, 0, 0, 0), tensor::zero),
-    F_("F", dimensionSet(0, -1, 0, 0, 0), tensor::zero)
-{
-    Info<< "Creating porous zone: " << key_ << endl;
-
-    bool foundZone = !cellZoneIds_.empty();
-    reduce(foundZone, orOp<bool>());
-
-    if (!foundZone && Pstream::master())
-    {
-        FatalErrorIn
-        (
-            "Foam::porousZone::porousZone"
-            "(const keyType&, const fvMesh&, const dictionary&)"
-        )   << "cannot find porous cellZone " << key_
-            << exit(FatalError);
-    }
-
-
-    // porosity
-    if
-    (
-        dict_.readIfPresent("porosity", porosity_)
-     && (porosity_ <= 0.0 || porosity_ > 1.0)
-    )
-    {
-        FatalIOErrorIn
-        (
-            "Foam::porousZone::porousZone"
-            "(const keyType&, const fvMesh&, const dictionary&)",
-            dict_
-        )
-            << "out-of-range porosity value " << porosity_
-            << exit(FatalIOError);
-    }
-
-    // turbulent intensity
-    if
-    (
-        dict_.readIfPresent("intensity", intensity_)
-     && (intensity_ <= 0.0 || intensity_ > 1.0)
-    )
-    {
-        FatalIOErrorIn
-        (
-            "Foam::porousZone::porousZone"
-            "(const keyType&, const fvMesh&, const dictionary&)",
-            dict_
-        )
-            << "out-of-range turbulent intensity value " << intensity_
-            << exit(FatalIOError);
-    }
-
-    // turbulent length scale
-    if
-    (
-        dict_.readIfPresent("mixingLength", mixingLength_)
-     && (mixingLength_ <= 0.0)
-    )
-    {
-        FatalIOErrorIn
-        (
-            "Foam::porousZone::porousZone"
-            "(const keyType&, const fvMesh&, const dictionary&)",
-            dict_
-        )
-            << "out-of-range turbulent length scale " << mixingLength_
-            << exit(FatalIOError);
-    }
-
-
-    // powerLaw coefficients
-    if (const dictionary* dictPtr = dict_.subDictPtr("powerLaw"))
-    {
-        dictPtr->readIfPresent("C0", C0_);
-        dictPtr->readIfPresent("C1", C1_);
-    }
-
-    // Darcy-Forchheimer coefficients
-    if (const dictionary* dictPtr = dict_.subDictPtr("Darcy"))
-    {
-        // local-to-global transformation tensor
-        const tensor& E = coordSys_.R();
-
-        dimensionedVector d(vector::zero);
-        if (dictPtr->readIfPresent("d", d))
-        {
-            if (D_.dimensions() != d.dimensions())
-            {
-                FatalIOErrorIn
-                (
-                    "Foam::porousZone::porousZone"
-                    "(const keyType&, const fvMesh&, const dictionary&)",
-                    dict_
-                )   << "incorrect dimensions for d: " << d.dimensions()
-                    << " should be " << D_.dimensions()
-                    << exit(FatalIOError);
-            }
-
-            adjustNegativeResistance(d);
-
-            D_.value().xx() = d.value().x();
-            D_.value().yy() = d.value().y();
-            D_.value().zz() = d.value().z();
-            D_.value() = (E & D_ & E.T()).value();
-        }
-
-        dimensionedVector f(vector::zero);
-        if (dictPtr->readIfPresent("f", f))
-        {
-            if (F_.dimensions() != f.dimensions())
-            {
-                FatalIOErrorIn
-                (
-                    "Foam::porousZone::porousZone"
-                    "(const keyType&, const fvMesh&, const dictionary&)",
-                    dict_
-                )   << "incorrect dimensions for f: " << f.dimensions()
-                    << " should be " << F_.dimensions()
-                    << exit(FatalIOError);
-            }
-
-            adjustNegativeResistance(f);
-
-            // leading 0.5 is from 1/2 * rho
-            F_.value().xx() = 0.5*f.value().x();
-            F_.value().yy() = 0.5*f.value().y();
-            F_.value().zz() = 0.5*f.value().z();
-            F_.value() = (E & F_ & E.T()).value();
-        }
-    }
-
-    // it is an error not to define anything
-    if
-    (
-        C0_ <= VSMALL
-     && magSqr(D_.value()) <= VSMALL
-     && magSqr(F_.value()) <= VSMALL
-    )
-    {
-        FatalIOErrorIn
-        (
-            "Foam::porousZone::porousZone"
-            "(const keyType&, const fvMesh&, const dictionary&)",
-            dict_
-        )   << "neither powerLaw (C0/C1) "
-               "nor Darcy-Forchheimer law (d/f) specified"
-            << exit(FatalIOError);
-    }
-
-    // feedback for the user
-    if (dict.lookupOrDefault("printCoeffs", false))
-    {
-        writeDict(Info, false);
-    }
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::porousZone::addResistance(fvVectorMatrix& UEqn) const
-{
-    if (cellZoneIds_.empty())
-    {
-        return;
-    }
-
-    bool compressible = false;
-    if (UEqn.dimensions() == dimensionSet(1, 1, -2, 0, 0))
-    {
-        compressible = true;
-    }
-
-    const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
-    vectorField& Usource = UEqn.source();
-    const vectorField& U = UEqn.psi();
-
-    if (C0_ > VSMALL)
-    {
-        if (compressible)
-        {
-            addPowerLawResistance
-            (
-                Udiag,
-                V,
-                mesh_.lookupObject<volScalarField>("rho"),
-                U
-            );
-        }
-        else
-        {
-            addPowerLawResistance
-            (
-                Udiag,
-                V,
-                geometricOneField(),
-                U
-            );
-        }
-    }
-
-    const tensor& D = D_.value();
-    const tensor& F = F_.value();
-
-    if (magSqr(D) > VSMALL || magSqr(F) > VSMALL)
-    {
-        if (compressible)
-        {
-            addViscousInertialResistance
-            (
-                Udiag,
-                Usource,
-                V,
-                mesh_.lookupObject<volScalarField>("rho"),
-                mesh_.lookupObject<volScalarField>("mu"),
-                U
-            );
-        }
-        else
-        {
-            addViscousInertialResistance
-            (
-                Udiag,
-                Usource,
-                V,
-                geometricOneField(),
-                mesh_.lookupObject<volScalarField>("nu"),
-                U
-            );
-        }
-    }
-}
-
-
-void Foam::porousZone::addResistance
-(
-    fvVectorMatrix& UEqn,
-    const volScalarField& rho,
-    const volScalarField& mu
-) const
-{
-    if (cellZoneIds_.empty())
-    {
-        return;
-    }
-
-    const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
-    vectorField& Usource = UEqn.source();
-    const vectorField& U = UEqn.psi();
-
-    if (C0_ > VSMALL)
-    {
-        addPowerLawResistance
-        (
-            Udiag,
-            V,
-            rho,
-            U
-        );
-    }
-
-    const tensor& D = D_.value();
-    const tensor& F = F_.value();
-
-    if (magSqr(D) > VSMALL || magSqr(F) > VSMALL)
-    {
-        addViscousInertialResistance
-        (
-            Udiag,
-            Usource,
-            V,
-            rho,
-            mu,
-            U
-        );
-    }
-}
-
-void Foam::porousZone::addResistance
-(
-    const fvVectorMatrix& UEqn,
-    volTensorField& AU,
-    bool correctAUprocBC
-) const
-{
-    if (cellZoneIds_.empty())
-    {
-        return;
-    }
-
-    bool compressible = false;
-    if (UEqn.dimensions() == dimensionSet(1, 1, -2, 0, 0))
-    {
-        compressible = true;
-    }
-
-    const vectorField& U = UEqn.psi();
-
-    if (C0_ > VSMALL)
-    {
-        if (compressible)
-        {
-            addPowerLawResistance
-            (
-                AU,
-                mesh_.lookupObject<volScalarField>("rho"),
-                U
-            );
-        }
-        else
-        {
-            addPowerLawResistance
-            (
-                AU,
-                geometricOneField(),
-                U
-            );
-        }
-    }
-
-    const tensor& D = D_.value();
-    const tensor& F = F_.value();
-
-    if (magSqr(D) > VSMALL || magSqr(F) > VSMALL)
-    {
-        if (compressible)
-        {
-            addViscousInertialResistance
-            (
-                AU,
-                mesh_.lookupObject<volScalarField>("rho"),
-                mesh_.lookupObject<volScalarField>("mu"),
-                U
-            );
-        }
-        else
-        {
-            addViscousInertialResistance
-            (
-                AU,
-                geometricOneField(),
-                mesh_.lookupObject<volScalarField>("nu"),
-                U
-            );
-        }
-    }
-
-    if (correctAUprocBC)
-    {
-        // Correct the boundary conditions of the tensorial diagonal to ensure
-        // processor boundaries are correctly handled when AU^-1 is interpolated
-        // for the pressure equation.
-        AU.correctBoundaryConditions();
-    }
-}
-
-
-void Foam::porousZone::writeDict(Ostream& os, bool subDict) const
-{
-    if (subDict)
-    {
-        os  << indent << token::BEGIN_BLOCK << incrIndent << nl;
-        os.writeKeyword("name")
-            << zoneName() << token::END_STATEMENT << nl;
-    }
-    else
-    {
-        os  << indent << zoneName() << nl
-            << indent << token::BEGIN_BLOCK << incrIndent << nl;
-    }
-
-    if (dict_.found("note"))
-    {
-        os.writeKeyword("note")
-            << string(dict_.lookup("note")) << token::END_STATEMENT << nl;
-    }
-
-    coordSys_.writeDict(os, true);
-
-    if (dict_.found("porosity"))
-    {
-        os.writeKeyword("porosity")
-            << porosity() << token::END_STATEMENT << nl;
-    }
-
-    if (dict_.found("intensity"))
-    {
-        os.writeKeyword("intensity")
-            << intensity() << token::END_STATEMENT << nl;
-    }
-
-    if (dict_.found("mixingLength"))
-    {
-        os.writeKeyword("mixingLength")
-            << mixingLength() << token::END_STATEMENT << nl;
-    }
-
-    // powerLaw coefficients
-    if (const dictionary* dictPtr = dict_.subDictPtr("powerLaw"))
-    {
-        os  << indent << "powerLaw";
-        dictPtr->write(os);
-    }
-
-    // Darcy-Forchheimer coefficients
-    if (const dictionary* dictPtr = dict_.subDictPtr("Darcy"))
-    {
-        os  << indent << "Darcy";
-        dictPtr->write(os);
-    }
-
-    // thermalModel
-    if (const dictionary* dictPtr = dict_.subDictPtr("thermalModel"))
-    {
-        os  << indent << "thermalModel";
-        dictPtr->write(os);
-    }
-
-    os  << decrIndent << indent << token::END_BLOCK << endl;
-}
-
-
-// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
-
-Foam::Ostream& Foam::operator<<(Ostream& os, const porousZone& pz)
-{
-    pz.writeDict(os);
-    return os;
-}
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H
deleted file mode 100644
index 431b36fe25033eb2ab5acb048d1fa447e830f4c6..0000000000000000000000000000000000000000
--- a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H
+++ /dev/null
@@ -1,389 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Namespace
-    Foam::porousMedia
-
-Description
-    Namespace for models related to porous media
-
-Class
-    Foam::porousZone
-
-Description
-    Porous zone definition based on cell zones.
-
-    Porous zone definition based on cell zones and parameters obtained from a
-    control dictionary constructed from the given stream. The orientation of
-    the porous region is defined with the same notation as a coordinateSystem,
-    but only a Cartesian coordinate system is valid.
-
-    Implemented porosity models:
-
-    powerLaw (\e C0 and \e C1 parameters)
-    \f[
-        S = - \rho C_0 |U|^{(C_1 - 1)} U
-    \f]
-
-    Darcy-Forchheimer (@e d and \e f parameters)
-    \f[
-        S = - (\mu \, d + \frac{\rho |U|}{2} \, f) U
-    \f]
-
-
-    Since negative Darcy/Forchheimer parameters are invalid, they can be used
-    to specify a multiplier (of the max component).
-
-    The porousZones method porousZones::ddt() mirrors the normal fvm::ddt()
-    method, but accounts for the effective volume of the cells.
-
-    An example dictionary entry:
-    \verbatim
-    cat1
-    {
-        note            "some catalyst";
-        coordinateSystem    system_10;
-        porosity        0.809;
-        intensity       0.001;      // optional
-        mixingLength    0.0001;     // optional
-        printCoeffs     yes;        // optional: feedback for the user
-
-        Darcy
-        {
-            d   d [0 -2 0 0 0]  (-1000 -1000 5.3756e+07);
-            f   f [0 -1 0 0 0]  (-1000 -1000 15.83);
-        }
-    }
-    \endverbatim
-
-See Also
-    porousZones and coordinateSystems
-
-SourceFiles
-    porousZone.C
-    porousZoneTemplates.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef porousZone_H
-#define porousZone_H
-
-#include "dictionary.H"
-#include "coordinateSystem.H"
-#include "coordinateSystems.H"
-#include "wordList.H"
-#include "labelList.H"
-#include "dimensionedScalar.H"
-#include "dimensionedTensor.H"
-#include "primitiveFieldsFwd.H"
-#include "volFieldsFwd.H"
-#include "fvMatricesFwd.H"
-
-#include "fvMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-class fvMesh;
-
-/*---------------------------------------------------------------------------*\
-                        Class porousZone Declaration
-\*---------------------------------------------------------------------------*/
-
-class porousZone
-{
-    // Private data
-
-        //- Name of this zone, or a regular expression
-        keyType key_;
-
-        //- Reference to the finite volume mesh this zone is part of
-        const fvMesh& mesh_;
-
-        //- Dictionary containing the parameters
-        dictionary dict_;
-
-        //- Cell zone Ids
-        labelList cellZoneIds_;
-
-        //- Coordinate system used for the zone (Cartesian)
-        coordinateSystem coordSys_;
-
-        //- porosity of the zone (0 < porosity <= 1)
-        //  Placeholder for treatment of temporal terms.
-        //  Currently unused.
-        scalar porosity_;
-
-        //- Turbulent intensity as fraction of the velocity
-        scalar intensity_;
-
-        //- Turbulent length scale
-        scalar mixingLength_;
-
-        //- powerLaw coefficient C0
-        scalar C0_;
-
-        //- powerLaw coefficient C1
-        scalar C1_;
-
-        //- Darcy coefficient
-        dimensionedTensor D_;
-
-        //- Forchheimer coefficient
-        dimensionedTensor F_;
-
-
-    // Private Member Functions
-
-        //- adjust negative resistance values to be multiplier of max value
-        static void adjustNegativeResistance(dimensionedVector& resist);
-
-        //- Power-law resistance
-        template<class RhoFieldType>
-        void addPowerLawResistance
-        (
-            scalarField& Udiag,
-            const scalarField& V,
-            const RhoFieldType& rho,
-            const vectorField& U
-        ) const;
-
-        //- Viscous and inertial resistance
-        template<class RhoFieldType>
-        void addViscousInertialResistance
-        (
-            scalarField& Udiag,
-            vectorField& Usource,
-            const scalarField& V,
-            const RhoFieldType& rho,
-            const scalarField& mu,
-            const vectorField& U
-        ) const;
-
-
-        //- Power-law resistance
-        template<class RhoFieldType>
-        void addPowerLawResistance
-        (
-            tensorField& AU,
-            const RhoFieldType& rho,
-            const vectorField& U
-        ) const;
-
-        //- Viscous and inertial resistance
-        template<class RhoFieldType>
-        void addViscousInertialResistance
-        (
-            tensorField& AU,
-            const RhoFieldType& rho,
-            const scalarField& mu,
-            const vectorField& U
-        ) const;
-
-
-        //- Disallow default bitwise copy construct
-        porousZone(const porousZone&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const porousZone&);
-
-
-public:
-
-    // Constructors
-
-        //- Construct from components
-        porousZone(const keyType& key, const fvMesh&, const dictionary&);
-
-        //- Return clone
-        autoPtr<porousZone> clone() const
-        {
-            notImplemented("autoPtr<porousZone> clone() const");
-            return autoPtr<porousZone>(NULL);
-        }
-
-        //- Return pointer to new porousZone created on freestore from Istream
-        class iNew
-        {
-            //- Reference to the finite volume mesh this zone is part of
-            const fvMesh& mesh_;
-
-        public:
-
-            iNew(const fvMesh& mesh)
-            :
-                mesh_(mesh)
-            {}
-
-            autoPtr<porousZone> operator()(Istream& is) const
-            {
-                keyType key(is);
-                dictionary dict(is);
-
-                return autoPtr<porousZone>(new porousZone(key, mesh_, dict));
-            }
-        };
-
-
-    //- Destructor
-    virtual ~porousZone()
-    {}
-
-
-    // Member Functions
-
-        // Access
-
-            //- cellZone name
-            const keyType& zoneName() const
-            {
-                return key_;
-            }
-
-            //- Return mesh
-            const fvMesh& mesh() const
-            {
-                return mesh_;
-            }
-
-            //- cellZone numbers
-            const labelList& zoneIds() const
-            {
-                return cellZoneIds_;
-            }
-
-            //- dictionary values used for the porousZone
-            const dictionary& dict() const
-            {
-                return dict_;
-            }
-
-            //- Return coordinate system
-            const coordinateSystem& coordSys() const
-            {
-                return coordSys_;
-            }
-
-            //- Return origin
-            const point& origin() const
-            {
-                return coordSys_.origin();
-            }
-
-            //- Return axis
-            vector axis() const
-            {
-                return coordSys_.axis();
-            }
-
-            //- Return porosity
-            scalar porosity() const
-            {
-                return porosity_;
-            }
-
-            //- Edit access to porosity
-            scalar& porosity()
-            {
-                return porosity_;
-            }
-
-            //- Return turbulent intensity
-            scalar intensity() const
-            {
-                return intensity_;
-            }
-
-            //- Edit access to turbulent intensity
-            scalar& intensity()
-            {
-                return intensity_;
-            }
-
-            //- Return turbulent length scale
-            scalar mixingLength() const
-            {
-                return mixingLength_;
-            }
-
-            //- Edit access to turbulent length scale
-            scalar& mixingLength()
-            {
-                return mixingLength_;
-            }
-
-
-        //- Modify time derivative elements according to porosity
-        template<class Type>
-        void modifyDdt(fvMatrix<Type>&) const;
-
-        //- Add the viscous and inertial resistance force contribution
-        //  to the momentum equation
-        void addResistance(fvVectorMatrix& UEqn) const;
-
-        //- Add the viscous and inertial resistance force contribution
-        //  to the momentum equation giving rho and mu.
-        void addResistance
-        (
-            fvVectorMatrix& UEqn,
-            const volScalarField& rho,
-            const volScalarField& mu
-        ) const;
-
-        //- Add the viscous and inertial resistance force contribution
-        //  to the tensorial diagonal.
-        //  Optionally correct the processor BCs of AU.
-        void addResistance
-        (
-            const fvVectorMatrix& UEqn,
-            volTensorField& AU,
-            bool correctAUprocBC = true
-        ) const;
-
-        //- Write the porousZone dictionary
-        virtual void writeDict(Ostream&, bool subDict = true) const;
-
-
-    // Ostream Operator
-
-        friend Ostream& operator<<(Ostream&, const porousZone&);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "porousZoneTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZoneTemplates.C b/src/finiteVolume/cfdTools/general/porousMedia/porousZoneTemplates.C
deleted file mode 100644
index 0e737c00bd645c701a35cae3270e7f48e467d0c4..0000000000000000000000000000000000000000
--- a/src/finiteVolume/cfdTools/general/porousMedia/porousZoneTemplates.C
+++ /dev/null
@@ -1,156 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*----------------------------------------------------------------------------*/
-
-#include "porousZone.H"
-#include "fvMesh.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class Type>
-void Foam::porousZone::modifyDdt(fvMatrix<Type>& m) const
-{
-    if (porosity_ < 1)
-    {
-        forAll(cellZoneIds_, zoneI)
-        {
-            const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
-
-            forAll(cells, i)
-            {
-                m.diag()[cells[i]]   *= porosity_;
-                m.source()[cells[i]] *= porosity_;
-            }
-        }
-    }
-}
-
-
-template<class RhoFieldType>
-void Foam::porousZone::addPowerLawResistance
-(
-    scalarField& Udiag,
-    const scalarField& V,
-    const RhoFieldType& rho,
-    const vectorField& U
-) const
-{
-    const scalar C0 = C0_;
-    const scalar C1m1b2 = (C1_ - 1.0)/2.0;
-
-    forAll(cellZoneIds_, zoneI)
-    {
-        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
-
-        forAll(cells, i)
-        {
-            Udiag[cells[i]] +=
-                V[cells[i]]*rho[cells[i]]*C0*pow(magSqr(U[cells[i]]), C1m1b2);
-        }
-    }
-}
-
-
-template<class RhoFieldType>
-void Foam::porousZone::addViscousInertialResistance
-(
-    scalarField& Udiag,
-    vectorField& Usource,
-    const scalarField& V,
-    const RhoFieldType& rho,
-    const scalarField& mu,
-    const vectorField& U
-) const
-{
-    const tensor& D = D_.value();
-    const tensor& F = F_.value();
-
-    forAll(cellZoneIds_, zoneI)
-    {
-        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
-
-        forAll(cells, i)
-        {
-            const tensor dragCoeff = mu[cells[i]]*D
-                + (rho[cells[i]]*mag(U[cells[i]]))*F;
-
-            const scalar isoDragCoeff = tr(dragCoeff);
-
-            Udiag[cells[i]] += V[cells[i]]*isoDragCoeff;
-            Usource[cells[i]] -=
-                V[cells[i]]*((dragCoeff - I*isoDragCoeff) & U[cells[i]]);
-        }
-    }
-}
-
-
-template<class RhoFieldType>
-void Foam::porousZone::addPowerLawResistance
-(
-    tensorField& AU,
-    const RhoFieldType& rho,
-    const vectorField& U
-) const
-{
-    const scalar C0 = C0_;
-    const scalar C1m1b2 = (C1_ - 1.0)/2.0;
-
-    forAll(cellZoneIds_, zoneI)
-    {
-        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
-
-        forAll(cells, i)
-        {
-            AU[cells[i]] = AU[cells[i]]
-                + I*(rho[cells[i]]*C0*pow(magSqr(U[cells[i]]), C1m1b2));
-        }
-    }
-}
-
-
-template<class RhoFieldType>
-void Foam::porousZone::addViscousInertialResistance
-(
-    tensorField& AU,
-    const RhoFieldType& rho,
-    const scalarField& mu,
-    const vectorField& U
-) const
-{
-    const tensor& D = D_.value();
-    const tensor& F = F_.value();
-
-    forAll(cellZoneIds_, zoneI)
-    {
-        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
-
-        forAll(cells, i)
-        {
-            AU[cells[i]] += mu[cells[i]]*D + (rho[cells[i]]*mag(U[cells[i]]))*F;
-        }
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H b/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H
deleted file mode 100644
index 17c439d1745ddba3fd6a7f1a166f64210bf0c3d7..0000000000000000000000000000000000000000
--- a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H
+++ /dev/null
@@ -1,46 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Typedef
-    Foam::porousZones
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef porousZones_H
-#define porousZones_H
-
-#include "PorousZones.H"
-#include "porousZone.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    typedef PorousZones<porousZone> porousZones;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/SubModelBase.H b/src/lagrangian/intermediate/submodels/SubModelBase.H
index 0427c162f7e305d317fd21fdd8e4f04ee80617d6..608cee614ebd59ec6950bf6ee3e83a950b0f0741 100644
--- a/src/lagrangian/intermediate/submodels/SubModelBase.H
+++ b/src/lagrangian/intermediate/submodels/SubModelBase.H
@@ -61,7 +61,7 @@ protected:
         CloudType& owner_;
 
         //- Reference to the cloud dictionary
-        const dictionary& dict_;
+        const dictionary dict_;
 
         //- Name of the sub-model base class
         const word baseName_;
@@ -73,7 +73,7 @@ protected:
         const word modelName_;
 
         //- Coefficients dictionary
-        const dictionary& coeffDict_;
+        const dictionary coeffDict_;
 
 
     // Protected Member Functions
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
index 7e85e6723e3db77181ac9ac223de233a26373196..1f9f4c66ec68de7c06a16114e741bd9d5449a2f0 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
@@ -266,7 +266,6 @@ Foam::label Foam::AMIInterpolation<SourcePatch, TargetPatch>::findTargetFace
     const point srcPt = srcFace.centre(srcPts);
     const scalar srcFaceArea = srcMagSf_[srcFaceI];
 
-//    pointIndexHit sample = treePtr_->findNearest(srcPt, sqr(0.1*bb.mag()));
     pointIndexHit sample = treePtr_->findNearest(srcPt, 10.0*srcFaceArea);
 
 
@@ -295,7 +294,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::appendNbrFaces
     DynamicList<label>& faceIDs
 ) const
 {
-//    const labelList& nbrFaces = patch.pointFaces()[faceI];
     const labelList& nbrFaces = patch.faceFaces()[faceI];
 
     // filter out faces already visited from src face neighbours
@@ -401,7 +399,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
     const DynamicList<label>& visitedFaces
 ) const
 {
-//    const labelList& srcNbrFaces = srcPatch0.pointFaces()[srcFaceI];
     const labelList& srcNbrFaces = srcPatch0.faceFaces()[srcFaceI];
 
     // set possible seeds for later use
@@ -610,9 +607,12 @@ restartUncoveredSourceFace
         }
     }
 
-    Info<< "AMIInterpolation : restarting search on "
-        << returnReduce(lowWeightFaces.size(), sumOp<label>())
-        << " faces since sum of weights < 0.5" << endl;
+    if (debug)
+    {
+        Pout<< "AMIInterpolation: restarting search on "
+            << lowWeightFaces.size() << " faces since sum of weights < 0.5"
+            << endl;
+    }
 
     if (lowWeightFaces.size() > 0)
     {
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
index fb15396f0f40b2031f4ab1d30beba7909a7b0226..59ff974315322ef3ac6bcad5e39631ee047aa9e0 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
@@ -211,6 +211,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
 
     dictionary forcesDict;
 
+    forcesDict.add("type", forces::typeName);
     forcesDict.add("patches", wordList(1, ptPatch.name()));
     forcesDict.add("rhoInf", rhoInf_);
     forcesDict.add("rhoName", rhoName_);
diff --git a/src/regionModels/thermoBaffleModels/derivedFvPatchFields/temperatureThermoBaffle/temperatureThermoBaffleFvPatchScalarField.H b/src/regionModels/thermoBaffleModels/derivedFvPatchFields/temperatureThermoBaffle/temperatureThermoBaffleFvPatchScalarField.H
index 8d58f73aa5a6d9eb012eb2d2dc279719c847ed00..b32a0127c4aea39ee7eb3d22dffcee9a3d325ffc 100644
--- a/src/regionModels/thermoBaffleModels/derivedFvPatchFields/temperatureThermoBaffle/temperatureThermoBaffleFvPatchScalarField.H
+++ b/src/regionModels/thermoBaffleModels/derivedFvPatchFields/temperatureThermoBaffle/temperatureThermoBaffleFvPatchScalarField.H
@@ -28,78 +28,91 @@ Group
     grpThermoBoundaryConditions
 
 Description
-    This boundary condition provides a coupled condition between the primary
-    baffle mesh regions.
+    This boundary condition provides a coupled temperature condition between
+    multiple mesh regions.  The regions are generally referred to as the:
+    - primary region, and
+    - baffle region
 
-    The primary region creates it and evolves the thermal baffle heat transfer
-    equation.
+    The primary region creates the baffle region and evolves its energy
+    equation either:
+    - 1-D, normal to each patch face
+    - 2-D, normal and tangential components
 
-    The solid thermo and baffle dictionaries are located on the
-    primary region.
+    The thermodynamic properties of the baffle material are specified via
+    dictionary entries.
 
-    type                compressible::temperatureThermoBaffle;
+    \heading Patch usage
 
-    // Coupled BC.
-    neighbourFieldName  T;
-    kappa               fluidThermo;
-    KName               none;
-
-
-    // Thermo baffle model
-    thermoBaffleModel   thermoBaffle2D;
-    regionName          baffleRegion;
-    infoOutput          yes;
-    active              yes;
-    thermoBaffle2DCoeffs
+    Example of the boundary condition specification:
+    \verbatim
+    myPatch
     {
-    }
+        type                compressible::temperatureThermoBaffle;
 
+        // Coupled boundary condition
+        neighbourFieldName  T;
+        kappa               fluidThermo;
+        KName               none;
 
-    // Solid thermo
-    thermoType
-    heSolidThermo
-        <pureSolidMixture
-            <constIsoSolidTransport
-                <constSolidRad
-                    <thermo
-                        <hConstThermo<incompressible>,sensibleEnthalpy>
-                    >
-                 >
-            >
-         >;
-
-    mixture
-    {
-        specie
-        {
-            nMoles          1;
-            molWeight       20;
-        }
-        transport
+
+        // Thermo baffle model
+        thermoBaffleModel   thermoBaffle2D;
+        regionName          baffleRegion;
+        infoOutput          yes;
+        active              yes;
+        thermoBaffle2DCoeffs
         {
-            kappa           0.01;
         }
-        radiation
+
+
+        // Solid thermo
+        thermoType
         {
-            sigmaS          0;
-            kappaRad        0;
-            emissivity      1;
+            type            heSolidThermo;
+            mixture         pureSolidMixture;
+            transport       constIso;
+            thermo          hConst;
+            equationOfState rhoConst;
+            specie          specie;
+            energy          sensibleEnthalpy;
         }
-        thermodynamics
+
+        mixture
         {
-            Hf              0;
-            Cp              15;
+            specie
+            {
+                nMoles          1;
+                molWeight       20;
+            }
+            transport
+            {
+                kappa           0.01;
+            }
+            thermodynamics
+            {
+                Hf              0;
+                Cp              15;
+            }
+            density
+            {
+                rho             80;
+            }
         }
-        density
+
+        radiation
         {
-            rho             80;
+            radiationModel  opaqueSolid;
+            absorptionEmissionModel none;
+            scatterModel    none;
         }
-    }
 
-    value               uniform 300;
+        value               uniform 300;
+    }
+    \endverbatim
 
 SeeAlso
     Foam::turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
+    Foam::regionModels::thermoBaffleModels::thermoBaffleModel
 
 SourceFiles
     temperatureThermoBaffleFvPatchScalarField.C
diff --git a/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.C b/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.C
index 3a67a8dd23bfc4786e6abff7845ce974ee233451..436457c4694bb2c717822f2ef418c9c9188d8799 100644
--- a/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.C
+++ b/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.C
@@ -87,14 +87,13 @@ void Foam::nastranSurfaceWriter::formatOS(OFstream& os) const
 void Foam::nastranSurfaceWriter::writeCoord
 (
     const point& p,
-    label& nPoint,
-    label& continuation,
-    Ostream& os
+    const label pointI,
+    OFstream& os
 ) const
 {
     // Fixed short/long formats:
     // 1 GRID
-    // 2 ID   : point ID
+    // 2 ID   : point ID - requires starting index of 1
     // 3 CP   : co-ordinate system ID                (blank)
     // 4 X1   : point x cp-ordinate
     // 5 X2   : point x cp-ordinate
@@ -107,35 +106,46 @@ void Foam::nastranSurfaceWriter::writeCoord
     {
         case wfShort:
         {
-            os  << setw(8) << "GRID"
-                << setw(8) << ++nPoint
+            os.setf(ios_base::left);
+            os  << setw(8) << "GRID";
+            os.unsetf(ios_base::left);
+            os.setf(ios_base::right);
+            os  << setw(8) << pointI + 1
                 << "        " 
                 << setw(8) << p.x()
                 << setw(8) << p.y()
                 << setw(8) << p.z()
                 << nl;
+            os.unsetf(ios_base::right);
 
             break;
         }
         case wfLong:
         {
-            os  << setw(8) << "GRID*"
-                << setw(16) << ++nPoint
+            os.setf(ios_base::left);
+            os  << setw(8) << "GRID*";
+            os.unsetf(ios_base::left);
+            os.setf(ios_base::right);
+            os  << setw(16) << pointI + 1
                 << "                "
                 << setw(16) << p.x()
                 << setw(16) << p.y()
-                << setw(8) << ++continuation
-                << nl
-                << setw(8) << continuation
-                << setw(16) << p.z()
                 << nl;
+            os.unsetf(ios_base::right);
+            os.setf(ios_base::left);
+            os  << setw(8) << "*";
+            os.unsetf(ios_base::left);
+            os.setf(ios_base::right);
+            os  << setw(16) << p.z()
+                << nl;
+            os.unsetf(ios_base::right);
 
             break;
         }
         case wfFree:
         {
             os  << "GRID"
-                << ',' << ++nPoint
+                << ',' << pointI + 1
                 << ','
                 << ',' << p.x()
                 << ',' << p.y()
@@ -163,7 +173,7 @@ void Foam::nastranSurfaceWriter::writeFace
     const word& faceType,
     const labelList& facePts,
     label& nFace,
-    Ostream& os
+    OFstream& os
 ) const
 {
     // Only valid surface elements are CTRIA3 and CQUAD4
@@ -172,7 +182,7 @@ void Foam::nastranSurfaceWriter::writeFace
     // 1 CQUAD4
     // 2 EID  : element ID
     // 3 PID  : property element ID; default = EID   (blank)
-    // 4 G1   : grid point index
+    // 4 G1   : grid point index - requires starting index of 1
     // 5 G2   : grid point index
     // 6 G3   : grid point index
     // 7 G4   : grid point index
@@ -183,18 +193,49 @@ void Foam::nastranSurfaceWriter::writeFace
     switch (writeFormat_)
     {
         case wfShort:
-        case wfLong:
         {
-            os  << setw(8) << faceType
-                << setw(8) << ++nFace
+            os.setf(ios_base::left);
+            os  << setw(8) << faceType;
+            os.unsetf(ios_base::left);
+            os.setf(ios_base::right);
+            os  << setw(8) << nFace++
                 << "        ";
 
             forAll(facePts, i)
             {
-                os  << setw(8) << facePts[i];
+                os  << setw(8) << facePts[i] + 1;
+            }
+
+            os  << nl;
+            os.unsetf(ios_base::right);
+
+            break;
+        }
+        case wfLong:
+        {
+            os.setf(ios_base::left);
+            os  << setw(8) << word(faceType + "*");
+            os.unsetf(ios_base::left);
+            os.setf(ios_base::right);
+            os  << setw(16) << nFace++
+                << "                ";
+
+            forAll(facePts, i)
+            {
+                os  << setw(16) << facePts[i] + 1;
+                if (i == 1)
+                {
+                    os  << nl;
+                    os.unsetf(ios_base::right);
+                    os.setf(ios_base::left);
+                    os  << setw(8) << "*";
+                    os.unsetf(ios_base::left);
+                    os.setf(ios_base::right);
+                }
             }
 
             os  << nl;
+            os.unsetf(ios_base::right);
 
             break;
         }
@@ -205,7 +246,7 @@ void Foam::nastranSurfaceWriter::writeFace
 
             forAll(facePts, i)
             {
-                os  << ',' << facePts[i];
+                os  << ',' << facePts[i] + 1;
             }
 
             os  << nl;
@@ -235,7 +276,7 @@ void Foam::nastranSurfaceWriter::writeGeometry
     const pointField& points,
     const faceList& faces,
     List<DynamicList<face> >& decomposedFaces,
-    Ostream& os
+    OFstream& os
 ) const
 {
     // write points
@@ -244,12 +285,9 @@ void Foam::nastranSurfaceWriter::writeGeometry
         << "$ Points" << nl
         << "$" << nl;
 
-    label nPoint = 0;
-    label continuation = 0;
-
     forAll(points, pointI)
     {
-        writeCoord(points[pointI], nPoint, continuation, os);
+        writeCoord(points[pointI], pointI, os);
     }
 
 
@@ -259,7 +297,7 @@ void Foam::nastranSurfaceWriter::writeGeometry
         << "$ Faces" << nl
         << "$" << nl;
 
-    label nFace = 0;
+    label nFace = 1;
 
     forAll(faces, faceI)
     {
@@ -352,7 +390,7 @@ void Foam::nastranSurfaceWriter::write
         Info<< "Writing nastran file to " << os.name() << endl;
     }
 
-    os  << "TITLE=OpeNFOAM " << surfaceName.c_str() << " mesh" << nl
+    os  << "TITLE=OpenFOAM " << surfaceName.c_str() << " mesh" << nl
         << "$" << nl
         << "BEGIN BULK" << nl;
 
diff --git a/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.H b/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.H
index 3667173d65b4da3104e7a98fd02d85a1e3c866d3..4e968ee55f4c53c3c06377590e8a3df73f5c1c3f 100644
--- a/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.H
+++ b/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.H
@@ -85,9 +85,8 @@ private:
         void writeCoord
         (
             const point& p,
-            label& nPoint,
-            label& continuation,
-            Ostream& os
+            const label pointI,
+            OFstream& os
         ) const;
 
         //- Write a face element (CTRIA3 or CQUAD4)
@@ -96,7 +95,7 @@ private:
             const word& faceType,
             const labelList& facePts,
             label& nFace,
-            Ostream& os
+            OFstream& os
         ) const;
 
         //- Main driver to write the surface mesh geometry
@@ -105,7 +104,7 @@ private:
             const pointField& points,
             const faceList& faces,
             List<DynamicList<face> >& decomposedFaces,
-            Ostream& os
+            OFstream& os
         ) const;
 
         //- Write a face-based value
@@ -115,7 +114,7 @@ private:
             const word& nasFieldName,
             const Type& value,
             const label EID,
-            Ostream& os
+            OFstream& os
         ) const;
 
         //- Templated write operation
diff --git a/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriterTemplates.C b/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriterTemplates.C
index a393a5968b37cc0f3befa8e70898ec0528ea6030..e01429fc88cbea5018d7f5bcb128927df49341ea 100644
--- a/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriterTemplates.C
+++ b/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriterTemplates.C
@@ -34,7 +34,7 @@ void Foam::nastranSurfaceWriter::writeFaceValue
     const word& nasFieldName,
     const Type& value,
     const label EID,
-    Ostream& os
+    OFstream& os
 ) const
 {
     // Fixed short/long formats:
@@ -43,26 +43,50 @@ void Foam::nastranSurfaceWriter::writeFaceValue
     // 3 EID  : element ID
     // 4 onwards: load values
 
-    label SID = 0;
+    label SID = 1;
 
-    label w = 16;
     switch (writeFormat_)
     {
         case wfShort:
         {
-            w = 8;
+            os.setf(ios_base::left);
+            os  << setw(8) << nasFieldName;
+            os.unsetf(ios_base::left);
+            os.setf(ios_base::right);
+            os  << setw(8) << SID
+                << setw(8) << EID;
+
+            for (direction dirI = 0; dirI < pTraits<Type>::nComponents; dirI++)
+            {
+                os  << setw(8) << component(value, dirI);
+            }
+
+            os.unsetf(ios_base::right);
+
+            break;
         }
         case wfLong:
         {
-            os  << setw(8) << nasFieldName
-                << setw(8) << SID
-                << setw(8) << EID;
+            os.setf(ios_base::left);
+            os  << setw(8) << word(nasFieldName + "*");
+            os.unsetf(ios_base::left);
+            os.setf(ios_base::right);
+            os  << setw(16) << SID
+                << setw(16) << EID;
 
             for (direction dirI = 0; dirI < pTraits<Type>::nComponents; dirI++)
             {
-                os  << setw(w) << component(value, dirI);
+                os  << setw(16) << component(value, dirI);
             }
 
+            os.unsetf(ios_base::right);
+
+            os  << nl;
+
+            os.setf(ios_base::left);
+            os  << '*';
+            os.unsetf(ios_base::left);
+
             break;
         }
         case wfFree:
diff --git a/src/thermophysicalModels/Allwmake b/src/thermophysicalModels/Allwmake
index 72255b94397fc4131b0f5d90339259fa08a2af6d..3a8da7148a3b1603702abe29aaeffb4c403ffced 100755
--- a/src/thermophysicalModels/Allwmake
+++ b/src/thermophysicalModels/Allwmake
@@ -12,7 +12,6 @@ wmake $makeType reactionThermo
 wmake $makeType laminarFlameSpeed
 wmake $makeType chemistryModel
 wmake $makeType barotropicCompressibilityModel
-wmake $makeType thermalPorousZone
 wmake $makeType SLGThermo
 
 wmake $makeType solidSpecie
diff --git a/src/thermophysicalModels/thermalPorousZone/Make/files b/src/thermophysicalModels/thermalPorousZone/Make/files
deleted file mode 100644
index 32e0e201f556092a8c0a2516f0a943b913ed74de..0000000000000000000000000000000000000000
--- a/src/thermophysicalModels/thermalPorousZone/Make/files
+++ /dev/null
@@ -1,9 +0,0 @@
-thermalPorousZone/thermalPorousZone.C
-thermalPorousZone/thermalPorousZones.C
-
-thermalModel/thermalModel/thermalModel.C
-thermalModel/thermalModel/thermalModelNew.C
-thermalModel/fixedTemperature/fixedTemperature.C
-thermalModel/noThermalModel/noThermalModel.C
-
-LIB = $(FOAM_LIBBIN)/libthermalPorousZone
diff --git a/src/thermophysicalModels/thermalPorousZone/Make/options b/src/thermophysicalModels/thermalPorousZone/Make/options
deleted file mode 100644
index d407126e2aaadef234647e652dbf55b369a7aab4..0000000000000000000000000000000000000000
--- a/src/thermophysicalModels/thermalPorousZone/Make/options
+++ /dev/null
@@ -1,9 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
-
-LIB_LIBS = \
-    -lfluidThermophysicalModels \
-    -lmeshTools \
-    -lfiniteVolume
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.C b/src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.C
deleted file mode 100644
index b7fee8652025cc4cac109f79df3d00eab8b44311..0000000000000000000000000000000000000000
--- a/src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.C
+++ /dev/null
@@ -1,91 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fixedTemperature.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fluidThermo.H"
-#include "volFields.H"
-#include "fvMatrices.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace porousMedia
-{
-    defineTypeNameAndDebug(fixedTemperature, 0);
-
-    addToRunTimeSelectionTable
-    (
-        thermalModel,
-        fixedTemperature,
-        pZone
-    );
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::porousMedia::fixedTemperature::fixedTemperature(const porousZone& pZone)
-:
-    thermalModel(pZone),
-    T_(readScalar(thermalCoeffs_.lookup("T")))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor    * * * * * * * * * * * * * * //
-
-Foam::porousMedia::fixedTemperature::~fixedTemperature()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::porousMedia::fixedTemperature::addEnergySource
-(
-    const fluidThermo& thermo,
-    const volScalarField& rho,
-    fvScalarMatrix& hEqn
-) const
-{
-    const labelList& zones = pZone_.zoneIds();
-    if (zones.empty() || T_ < 0.0)
-    {
-        return;
-    }
-
-    const fvMesh& mesh = pZone_.mesh();
-    const scalarField T(hEqn.diag().size(), T_);
-
-    forAll(zones, zoneI)
-    {
-        const labelList& cells = mesh.cellZones()[zones[zoneI]];
-        hEqn.setValues(cells, thermo.he(thermo.p(), T, cells));
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.H b/src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.H
deleted file mode 100644
index f28317a9d452e45180d97d06b643c427ee2f268d..0000000000000000000000000000000000000000
--- a/src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.H
+++ /dev/null
@@ -1,102 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::porousMedia::fixedTemperature
-
-Description
-    Fixed temperature model for porous media
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef fixedTemperature_H
-#define fixedTemperature_H
-
-#include "thermalModel.H"
-#include "runTimeSelectionTables.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace porousMedia
-{
-
-/*---------------------------------------------------------------------------*\
-                      Class fixedTemperature Declaration
-\*---------------------------------------------------------------------------*/
-
-class fixedTemperature
-:
-    public thermalModel
-{
-
-protected:
-
-    // Protected data
-
-        //- Fixed temperature
-        const scalar T_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("fixedTemperature");
-
-
-    // Constructors
-
-        //- Construct from porous zone
-        fixedTemperature(const porousZone& pZone);
-
-
-    //- Destructor
-    virtual ~fixedTemperature();
-
-
-    // Member Functions
-
-        //- Add the thermal source to the enthalpy equation
-        virtual void addEnergySource
-        (
-            const fluidThermo&,
-            const volScalarField& rho,
-            fvScalarMatrix& hEqn
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace porousMedia
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/noThermalModel/noThermalModel.C b/src/thermophysicalModels/thermalPorousZone/thermalModel/noThermalModel/noThermalModel.C
deleted file mode 100644
index 4beef65f4fa410238ef9df44db0f27eeed6a4753..0000000000000000000000000000000000000000
--- a/src/thermophysicalModels/thermalPorousZone/thermalModel/noThermalModel/noThermalModel.C
+++ /dev/null
@@ -1,77 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "noThermalModel.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fluidThermo.H"
-#include "volFields.H"
-#include "fvMatrices.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace porousMedia
-{
-   defineTypeNameAndDebug(noThermalModel, 0);
-
-   addToRunTimeSelectionTable
-   (
-       thermalModel,
-       noThermalModel,
-       pZone
-   );
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::porousMedia::noThermalModel::noThermalModel(const porousZone& pZone)
-:
-    thermalModel(pZone)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor    * * * * * * * * * * * * * * //
-
-Foam::porousMedia::noThermalModel::~noThermalModel()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::porousMedia::noThermalModel::addEnergySource
-(
-    const fluidThermo&,
-    const volScalarField&,
-    fvScalarMatrix&
-) const
-{
-    // do nothing
-}
-
-
-// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.C b/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.C
deleted file mode 100644
index c8df5b4c24ab7af56dd8c3b3cd09f31714352d19..0000000000000000000000000000000000000000
--- a/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.C
+++ /dev/null
@@ -1,66 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "thermalModel.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace porousMedia
-{
-   defineTypeNameAndDebug(thermalModel, 0);
-   defineRunTimeSelectionTable(thermalModel, pZone);
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::porousMedia::thermalModel::thermalModel(const porousZone& pZone)
-:
-    pZone_(pZone),
-    thermalCoeffs_(pZone.dict().subDictPtr("thermalModel"))
-{}
-
-
-Foam::porousMedia::thermalModel::thermalModel
-(
-    const porousZone& pZone,
-    const dictionary& thermalCoeffs
-)
-:
-    pZone_(pZone),
-    thermalCoeffs_(thermalCoeffs)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor    * * * * * * * * * * * * * * //
-
-Foam::porousMedia::thermalModel::~thermalModel()
-{}
-
-
-// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.H b/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.H
deleted file mode 100644
index 26991dc40ad1d09be827c0da80e6c5df295f7b6b..0000000000000000000000000000000000000000
--- a/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.H
+++ /dev/null
@@ -1,129 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::porousMedia::thermalModel
-
-Description
-    Base class to select the temperature specification models for porousMedia
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef thermalModel_H
-#define thermalModel_H
-
-#include "porousZone.H"
-#include "autoPtr.H"
-#include "runTimeSelectionTables.H"
-#include "volFieldsFwd.H"
-#include "fvMatricesFwd.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class fluidThermo;
-
-namespace porousMedia
-{
-
-/*---------------------------------------------------------------------------*\
-                        Class thermalModel Declaration
-\*---------------------------------------------------------------------------*/
-
-class thermalModel
-{
-
-protected:
-
-    // Protected data
-
-        //- Reference to the porous zone
-        const porousZone& pZone_;
-
-        //- Thermal model coefficients dictionary
-        const dictionary thermalCoeffs_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("thermalModel");
-
-    //- Declare runtime constructor selection table
-    declareRunTimeSelectionTable
-    (
-        autoPtr,
-        thermalModel,
-        pZone,
-        (
-            const porousZone& pZone
-        ),
-        (pZone)
-    );
-
-
-    // Constructors
-
-        //- Construct from porous zone, coefficients from "thermalModel" entry
-        thermalModel(const porousZone&);
-
-        //- Construct from porous zone and thermal model coefficients
-        thermalModel(const porousZone&, const dictionary& thermalCoeffs);
-
-
-    //- Destructor
-    virtual ~thermalModel();
-
-
-    //- Selector
-    static autoPtr<thermalModel> New(const porousZone&);
-
-
-    // Member Functions
-
-        //- Add the thermal source to the enthalpy equation
-        virtual void addEnergySource
-        (
-            const fluidThermo&,
-            const volScalarField& rho,
-            fvScalarMatrix& hEqn
-        ) const = 0;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace porousMedia
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C
deleted file mode 100644
index c636fc9e8760aa93f41d9bed946d59f4a7a090cd..0000000000000000000000000000000000000000
--- a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C
+++ /dev/null
@@ -1,61 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*----------------------------------------------------------------------------*/
-
-#include "thermalPorousZone.H"
-#include "fluidThermo.H"
-#include "volFields.H"
-#include "fvMatrices.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::thermalPorousZone::thermalPorousZone
-(
-    const keyType& key,
-    const fvMesh& mesh,
-    const dictionary& dict
-)
-:
-    porousZone(key, mesh, dict),
-    model_(porousMedia::thermalModel::New(*this))
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::thermalPorousZone::addEnergySource
-(
-    const fluidThermo& thermo,
-    const volScalarField& rho,
-    fvScalarMatrix& hEqn
-) const
-{
-    if (model_.valid())
-    {
-        model_->addEnergySource(thermo, rho, hEqn);
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H
deleted file mode 100644
index 127d7fdf91cb4056b522fd2798bca836d7b1e472..0000000000000000000000000000000000000000
--- a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H
+++ /dev/null
@@ -1,145 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::thermalPorousZone
-
-Description
-    Porous zone definition based on cell zones including terms for energy
-    equations.
-
-See Also
-    porousZone, thermalPorousZones and coordinateSystems with run-time
-    selectable thermal model
-
-SourceFiles
-    thermalPorousZone.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef thermalPorousZone_H
-#define thermalPorousZone_H
-
-#include "porousZone.H"
-#include "thermalModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-class fvMesh;
-class fluidThermo;
-
-/*---------------------------------------------------------------------------*\
-                       Class thermalPorousZone Declaration
-\*---------------------------------------------------------------------------*/
-
-class thermalPorousZone
-:
-    public porousZone
-{
-    // Private data
-
-        //- Disallow default bitwise copy construct
-        thermalPorousZone(const thermalPorousZone&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const thermalPorousZone&);
-
-        //- Thermal model
-        autoPtr<porousMedia::thermalModel> model_;
-
-
-public:
-
-    // Constructors
-
-        //- Construct from components
-        thermalPorousZone
-        (
-            const keyType& key,
-            const fvMesh& mesh,
-            const dictionary& dict
-        );
-
-        //- Return clone
-        autoPtr<thermalPorousZone> clone() const
-        {
-            notImplemented("autoPtr<thermalPorousZone> clone() const");
-            return autoPtr<thermalPorousZone>(NULL);
-        }
-
-        //- Return pointer to new thermalPorousZone
-        //  created on freestore from Istream
-        class iNew
-        {
-            //- Reference to the finite volume mesh this zone is part of
-            const fvMesh& mesh_;
-
-        public:
-
-            iNew(const fvMesh& mesh)
-            :
-                mesh_(mesh)
-            {}
-
-            autoPtr<thermalPorousZone> operator()(Istream& is) const
-            {
-                keyType key(is);
-                dictionary dict(is);
-
-                return autoPtr<thermalPorousZone>
-                (
-                    new thermalPorousZone(key, mesh_, dict)
-                );
-            }
-        };
-
-
-    //- Destructor
-    virtual ~thermalPorousZone()
-    {}
-
-
-    // Member Functions
-
-        //- Add the thermal source to the enthalpy equation
-        void addEnergySource
-        (
-            const fluidThermo&,
-            const volScalarField& rho,
-            fvScalarMatrix& hEqn
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.H b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.H
deleted file mode 100644
index d239e0daf07c96d5fafe3ce420812d5269ee2aeb..0000000000000000000000000000000000000000
--- a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.H
+++ /dev/null
@@ -1,113 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::thermalPorousZones
-
-Description
-    A centralized thermalPorousZone collection.
-
-    Container class for a set of thermalPorousZones with the thermalPorousZone
-    member functions implemented to loop over the functions for each
-    thermalPorousZone.
-
-    The input file \c constant/thermalPorousZone is implemented as
-    IOPtrList\<thermalPorousZone\> (but written as a dictionary)
-    and contains the following type of data:
-
-    \verbatim
-    1
-    (
-    cat1
-    {
-        coordinateSystem    system_10;
-        porosity    0.781;
-        Darcy
-        {
-            d   d [0 -2 0 0 0]  (-1000 -1000 0.50753e+08);
-            f   f [0 -1 0 0 0]  (-1000 -1000 12.83);
-        }
-        thermalModel
-        {
-            type    fixedTemperature;
-            T       T [0 0 1 0 0] 600;
-        }
-    }
-    )
-    \endverbatim
-
-SourceFiles
-    thermalPorousZones.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef thermalPorousZones_H
-#define thermalPorousZones_H
-
-#include "PorousZones.H"
-#include "thermalPorousZone.H"
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                     Class thermalPorousZones Declaration
-\*---------------------------------------------------------------------------*/
-
-class thermalPorousZones
-:
-    public PorousZones<thermalPorousZone>
-{
-
-public:
-
-    // Constructors
-
-        //- Construct from fvMesh
-        thermalPorousZones(const fvMesh&);
-
-
-    // Member Functions
-
-        //- Add the thermal source to the enthalpy equation
-        void addEnergySource
-        (
-            const fluidThermo&,
-            const volScalarField& rho,
-            fvScalarMatrix& hEqn
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/porousZones b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/constant/porosityProperties
similarity index 78%
rename from tutorials/compressible/rhoPimplecFoam/angledDuct/constant/porousZones
rename to tutorials/compressible/rhoPimpleFoam/ras/angledDuct/constant/porosityProperties
index afeb6461e13249ef04a0d449173de79ecdd2111f..650f01268ca7bb4dfe14c25d5d8f2a860cda2544 100644
--- a/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/porousZones
+++ b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/constant/porosityProperties
@@ -11,26 +11,28 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      porousZones;
+    object      porosityProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-1
-(
-    porosity
+porosity1
+{
+    type            DarcyForchheimer;
+    active          yes;
+    cellZone        porosity;
+
+    DarcyForchheimerCoeffs
     {
+        d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
+        f   f [0 -1 0 0 0 0 0] (0 0 0);
+
         coordinateSystem
         {
             e1  (0.70710678 0.70710678 0);
             e2  (0 0 1);
         }
-
-        Darcy
-        {
-            d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
-            f   f [0 -1 0 0 0 0 0] (0 0 0);
-        }
     }
-)
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/constant/porousZones b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/porosityProperties
similarity index 78%
rename from tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/constant/porousZones
rename to tutorials/compressible/rhoPimplecFoam/angledDuct/constant/porosityProperties
index afeb6461e13249ef04a0d449173de79ecdd2111f..22add8ff299919833276bb7b4770c81785d1ace6 100644
--- a/tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/constant/porousZones
+++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/porosityProperties
@@ -11,26 +11,27 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      porousZones;
+    object      porosityProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-1
-(
-    porosity
+porosity1
+{
+    type            DarcyForchheimer;
+    active          yes;
+    cellZone        porosity;
+
+    DarcyForchheimerCoeffs
     {
+        d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
+        f   f [0 -1 0 0 0 0 0] (0 0 0);
+
         coordinateSystem
         {
             e1  (0.70710678 0.70710678 0);
             e2  (0 0 1);
         }
-
-        Darcy
-        {
-            d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
-            f   f [0 -1 0 0 0 0 0] (0 0 0);
-        }
     }
-)
+}
 
 // ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/porousZones b/tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/constant/porosityProperties
similarity index 78%
rename from tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/porousZones
rename to tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/constant/porosityProperties
index afeb6461e13249ef04a0d449173de79ecdd2111f..22add8ff299919833276bb7b4770c81785d1ace6 100644
--- a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/porousZones
+++ b/tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/constant/porosityProperties
@@ -11,26 +11,27 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      porousZones;
+    object      porosityProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-1
-(
-    porosity
+porosity1
+{
+    type            DarcyForchheimer;
+    active          yes;
+    cellZone        porosity;
+
+    DarcyForchheimerCoeffs
     {
+        d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
+        f   f [0 -1 0 0 0 0 0] (0 0 0);
+
         coordinateSystem
         {
             e1  (0.70710678 0.70710678 0);
             e2  (0 0 1);
         }
-
-        Darcy
-        {
-            d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
-            f   f [0 -1 0 0 0 0 0] (0 0 0);
-        }
     }
-)
+}
 
 // ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/porousZones b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/porosityProperties
similarity index 78%
rename from tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/porousZones
rename to tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/porosityProperties
index d3de954c172af5e6d736a533be90608037870272..1a76ef1d80b1b22f34ad537480509ea465b31411 100644
--- a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/porousZones
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/porosityProperties
@@ -11,28 +11,27 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      porousZones;
+    object      porosityProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-(
-/*
-    porousRegion // name of cell zone
+porosity1
+{
+    type            DarcyForchheimer;
+    active          yes;
+    cellZone        stator;
+
+    DarcyForchheimerCoeffs
     {
+        d   d [0 -2 0 0 0 0 0] (1e5 -1000 -1000);
+        f   f [0 -1 0 0 0 0 0] (0 0 0);
+
         coordinateSystem
         {
             e1  (1 0 0);
             e2  (0 1 0);
         }
-
-        Darcy
-        {
-            d   d [0 -2 0 0 0 0 0] (500000 -1000 -1000);
-            f   f [0 -1 0 0 0 0 0] (0 0 0);
-        }
     }
-*/
-);
-
+}
 
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/constant/porousZones b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/porosityProperties
similarity index 78%
rename from tutorials/compressible/rhoPimpleFoam/ras/angledDuct/constant/porousZones
rename to tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/porosityProperties
index afeb6461e13249ef04a0d449173de79ecdd2111f..650f01268ca7bb4dfe14c25d5d8f2a860cda2544 100644
--- a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/constant/porousZones
+++ b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/porosityProperties
@@ -11,26 +11,28 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      porousZones;
+    object      porosityProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-1
-(
-    porosity
+porosity1
+{
+    type            DarcyForchheimer;
+    active          yes;
+    cellZone        porosity;
+
+    DarcyForchheimerCoeffs
     {
+        d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
+        f   f [0 -1 0 0 0 0 0] (0 0 0);
+
         coordinateSystem
         {
             e1  (0.70710678 0.70710678 0);
             e2  (0 0 1);
         }
-
-        Darcy
-        {
-            d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
-            f   f [0 -1 0 0 0 0 0] (0 0 0);
-        }
     }
-)
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/porousZones b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/porousZones
deleted file mode 100644
index fe11354202595c2b3783d11d947f2e1ba775041b..0000000000000000000000000000000000000000
--- a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/porousZones
+++ /dev/null
@@ -1,44 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      porousZones;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-1
-(
-    porosity
-    {
-        coordinateSystem
-        {
-            e1  (0.70710678 0.70710678 0);
-            e2  (0 0 1);
-        }
-
-        Darcy
-        {
-            d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
-            f   f [0 -1 0 0 0 0 0] (0 0 0);
-        }
-
-        thermalModel
-        {
-            type        none; // fixedTemperature;
-
-            // fixedTemperature coefficients
-            T           350;
-        }
-    }
-)
-
-// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/porousZones b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/sourcesProperties
similarity index 74%
rename from tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/porousZones
rename to tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/sourcesProperties
index 11a47989c7f60537aecda39d1f437279dcde2766..4386bc52506b016f5aa089d5535b7d7c81ed7694 100644
--- a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/porousZones
+++ b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/sourcesProperties
@@ -11,26 +11,24 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      porousZones;
+    object      sourcesProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-1
-(
-    stator
-    {
-        coordinateSystem
-        {
-            e1  (1 0 0);
-            e2  (0 1 0);
-        }
+source1
+{
+    type            fixedTemperatureSource;
+    active          true;
+    timeStart       0;
+    duration        1000000;
+    selectionMode   cellZone;
+    cellZone        porosity;
 
-        Darcy
-        {
-            d   d [0 -2 0 0 0 0 0] (1e5 -1000 -1000);
-            f   f [0 -1 0 0 0 0 0] (0 0 0);
-        }
+    fixedTemperatureSourceCoeffs
+    {
+        temperature     350;
     }
-)
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSolution b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSolution
index 7cb8b518286b7a24bd86e560d8bdde399c3dd2db..27f43cb0bd350bc4263453a42298b4f68ae9942d 100644
--- a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSolution
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSolution
@@ -29,7 +29,7 @@ solvers
         mergeLevels     1;
     }
 
-    "(k|epsilon)"
+    "(U|k|epsilon)"
     {
         solver          smoothSolver;
         smoother        GaussSeidel;
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/sourcesProperties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/sourcesProperties
index bcbbcceba00b10a24aab9c4cf18ab083b6dc3b74..715f0ebdd9fc926737bb0f26ff740a4bada1bd4c 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/sourcesProperties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/sourcesProperties
@@ -17,20 +17,16 @@ FoamFile
 
 source1
 {
-    type            scalarExplicitSource;
+    type            fixedTemperatureSource;
     active          true;
     timeStart       0.1;
     duration        0.4;
     selectionMode   cellSet;
     cellSet         ignitionCells;
 
-    scalarExplicitSourceCoeffs
+    fixedTemperatureSourceCoeffs
     {
-        volumeMode      absolute;
-        injectionRate
-        {
-            h               12000;
-        }
+        temperature     2000;
     }
 }
 
diff --git a/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/constant/porousZones b/tutorials/lagrangian/reactingParcelFoam/filter/constant/porosityProperties
similarity index 86%
rename from tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/constant/porousZones
rename to tutorials/lagrangian/reactingParcelFoam/filter/constant/porosityProperties
index d3de954c172af5e6d736a533be90608037870272..fcd5c5d3f1724961ee4de6a065ec6dbdcebad16f 100644
--- a/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/constant/porousZones
+++ b/tutorials/lagrangian/reactingParcelFoam/filter/constant/porosityProperties
@@ -11,28 +11,31 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      porousZones;
+    object      porosityProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-(
-/*
-    porousRegion // name of cell zone
-    {
-        coordinateSystem
-        {
-            e1  (1 0 0);
-            e2  (0 1 0);
-        }
+filter1
+{
+    cellZone        filter;
+    active          true;
+    type            DarcyForchheimer;
 
+    DarcyForchheimerCoeffs
+    {
         Darcy
         {
             d   d [0 -2 0 0 0 0 0] (500000 -1000 -1000);
             f   f [0 -1 0 0 0 0 0] (0 0 0);
         }
+
+        coordinateSystem
+        {
+            e1  (1 0 0);
+            e2  (0 1 0);
+        }
     }
-*/
-);
+}
 
 
 // ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/constant/porousZones b/tutorials/lagrangian/reactingParcelFoam/filter/constant/porousZones
deleted file mode 100644
index 32d3aca74f0525259ac4197dbb2a0386eb5e5140..0000000000000000000000000000000000000000
--- a/tutorials/lagrangian/reactingParcelFoam/filter/constant/porousZones
+++ /dev/null
@@ -1,36 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      porousZones;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-(
-    filter // name of cell zone
-    {
-        coordinateSystem
-        {
-            e1  (1 0 0);
-            e2  (0 1 0);
-        }
-
-        Darcy
-        {
-            d   d [0 -2 0 0 0 0 0] (500000 -1000 -1000);
-            f   f [0 -1 0 0 0 0 0] (0 0 0);
-        }
-    }
-)
-
-
-// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/porousZones b/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/porousZones
deleted file mode 100644
index d3de954c172af5e6d736a533be90608037870272..0000000000000000000000000000000000000000
--- a/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/porousZones
+++ /dev/null
@@ -1,38 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      porousZones;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-(
-/*
-    porousRegion // name of cell zone
-    {
-        coordinateSystem
-        {
-            e1  (1 0 0);
-            e2  (0 1 0);
-        }
-
-        Darcy
-        {
-            d   d [0 -2 0 0 0 0 0] (500000 -1000 -1000);
-            f   f [0 -1 0 0 0 0 0] (0 0 0);
-        }
-    }
-*/
-);
-
-
-// ************************************************************************* //