diff --git a/Makefile b/Makefile
index 38f354d..da64a96 100644
--- a/Makefile
+++ b/Makefile
@@ -4,7 +4,7 @@ CXXFLAGS = 	$(PRJCXXFLAGS)
 LDLIBS   =  ./src/libac_core.a $(PRJLDFLAGS)
 
 VPATH=./src
-MAIN = apfel_comb src/cfac_scale 
+MAIN = apfel_comb apfel_evol src/cfac_scale 
 DEV = appl_optgrid ftdy_hcx
 
 .PHONY: all dev core clean
diff --git a/src/apfel_comb.cc b/src/apfel_comb.cc
index f9e2473..22e873e 100644
--- a/src/apfel_comb.cc
+++ b/src/apfel_comb.cc
@@ -65,6 +65,7 @@ int main(int argc, char* argv[]) {
 
   // // Initialise QCD
   if (table.GetSource() == FKTarget::DYP) QCD::setFTDYmode(true);
+  if (table.GetSource() == FKTarget::DIS) QCD::setDISmode(true);
   QCD::initQCD(par, table.GetPositivity(), table.GetQ2max());
   QCD::initEvolgrid(table.GetNX(), table.GetXmin());
 
diff --git a/src/apfel_evol.cc b/src/apfel_evol.cc
new file mode 100644
index 0000000..3d55639
--- /dev/null
+++ b/src/apfel_evol.cc
@@ -0,0 +1,110 @@
+// apfel_evol.cc
+// Code to generate evolution-grid FK tables
+
+#include "LHAPDF/LHAPDF.h"
+
+#include <vector>
+#include <iostream>
+#include <cstdlib>
+#include <fstream>
+#include <string>
+#include <time.h>
+#include <sys/time.h>
+#include <cstdio>
+#include <map>
+
+#include "apfelcomb/fk_evol.h"
+
+#include "apfelcomb/fk_utils.h"
+#include "apfelcomb/fk_pdf.h"
+#include "apfelcomb/fk_qcd.h"
+
+using namespace std;
+
+
+int main(int argc, char* argv[]) { if ( argc != 2 )
+    {
+        std::cerr << "Usage: " <<argv[0] <<" <target ThID>" <<std::endl;
+        exit(-1);
+    }
+
+    Splash();
+
+    // target theoryIDs
+    const int th = atoi(argv[1]);
+    QCD::qcd_param qcdPar;
+    QCD::parse_input(th, qcdPar);
+    setupDir(th);
+
+    const double Q2Min = pow(qcdPar.Q0,2);
+    const double Q2Max = 1E5*1E5;
+
+    const int nx_in   = 195;
+    const int nx_out  = 100;
+
+    const double xmin = 1E-9;
+    const double xmed = 0.1;
+    const double xmax = 1.0;
+
+    // NNPDF3.1 standard grid uses a = 30 (fished out from an old email)
+    QCD::initQCD(qcdPar, false, Q2Max);
+    QCD::initEvolgrid(nx_in, xmin, 30);
+
+    // Generate target x-grid
+    const vector<double> xg_out = EVL::generate_xgrid(nx_out, xmin, xmed, xmax);
+    const vector< vector<double> > q2_subgrids = EVL::generate_q2subgrids(qcdPar, 50, Q2Min, Q2Max);
+
+    // Generate and loop over subgrids
+    for (size_t i=0; i< q2_subgrids.size(); i++)
+    {
+        const vector<double> q2_out = q2_subgrids[i];
+        const int nfl     = 14; // All for now
+        const int npoints = q2_out.size()*xg_out.size()*nfl;
+        const dataInfoRaw subgrid_info = {npoints, 0, "EVOL", "PDF"};
+        const EVL::EvolutionData ed(subgrid_info, xg_out, q2_out);
+
+        // Initialise empty mFK table
+        NNPDF::FKHeader FKhead;
+        //
+        //table.SetFKHeader(FKhead);
+        QCD::set_params(qcdPar, FKhead);
+
+        const string setname     = "PDFEVOL"+std::to_string(i);
+        const string description = "PDF Evolution test grid";
+        FKhead.AddTag(NNPDF::FKHeader::BLOB,     "GridDesc", description);
+        FKhead.AddTag(NNPDF::FKHeader::GRIDINFO, "SETNAME", setname);
+        FKhead.AddTag(NNPDF::FKHeader::GRIDINFO, "NDATA",   ed.GetNData());
+        FKhead.AddTag(NNPDF::FKHeader::VERSIONS, "APPLrepo", applCommit() );
+        FKhead.AddTag(NNPDF::FKHeader::GRIDINFO, "HADRONIC", false);
+        FKhead.ResetFlavourMap();
+
+        std::stringstream IO; FKhead.Print(IO);
+        NNPDF::FKGenerator* FK = new NNPDF::FKGenerator( IO );
+
+        for (int i=0; i<ed.GetNData(); i++)
+        {
+            const double xo = ed.GetKinematics(i,0);
+            const double Q  = sqrt(ed.GetKinematics(i,1));
+            const double fl = ed.GetKinematics(i,2);
+            std::array<double, 14> pdf;
+
+            for(int xi=0; xi<FK->GetNx(); xi++)
+                for(int ifl=0; ifl<14; ifl++)
+                {
+                    // Including photon
+                    QCD::EvolutionOperator(false, true, xi, xo, ifl, Q, pdf.data());
+                    FK->Fill(i, xi, ifl, pdf[fl]);
+                }
+            std::cout << i <<"/"<<ed.GetNData() <<"\r";
+            std::cout.flush();
+        }
+
+        // Write to file
+        const std::string outFile = resultsPath()+"theory_" + to_string(th) + "/evolution/FK_"+setname+ ".dat";
+        FK->Print(outFile, true);
+        delete FK;
+    }
+
+    exit(0);
+}
+
diff --git a/src/apfelcomb/fk_evol.h b/src/apfelcomb/fk_evol.h
new file mode 100644
index 0000000..f21e6fe
--- /dev/null
+++ b/src/apfelcomb/fk_evol.h
@@ -0,0 +1,43 @@
+#pragma once
+
+#include <stdlib.h>
+#include <fstream>
+#include <vector>
+
+// NNPDF
+#include "NNPDF/fkgenerator.h"
+#include "NNPDF/commondata.h"
+#include "NNPDF/nnpdfdb.h"
+
+#include "fk_utils.h"
+#include "fk_qcd.h"
+#include "fk_grids.h"
+
+using NNPDF::dataInfoRaw;
+using NNPDF::CommonData;
+
+namespace EVL
+{
+    // Class storing the required kinematic information for PDF evolution
+    class EvolutionData: public CommonData
+    {
+        public:
+        EvolutionData( dataInfoRaw const& di,
+                       vector<double> const& xgrid,
+                       vector<double> const& q2grid );
+
+    };
+
+    // Generate output (evolved scale) x-grid
+    vector<double> generate_xgrid( const int nx,
+            const double xmin,
+            const double xmed,
+            const double xmax);
+
+    // Generate output q2 subgrids
+    // One subgrid is generated for every HQ threshold crossed
+    vector< vector<double> > generate_q2subgrids( QCD::qcd_param par, const int nq2,
+            const double q2min, const double q2max);
+
+}
+
diff --git a/src/apfelcomb/fk_qcd.h b/src/apfelcomb/fk_qcd.h
index 5cb9cf3..aa352ac 100644
--- a/src/apfelcomb/fk_qcd.h
+++ b/src/apfelcomb/fk_qcd.h
@@ -103,7 +103,7 @@ namespace QCD
 
   // Initialise the APFEL interface
   void initQCD(qcd_param& par, const bool& positivity, const double& Q2max);
-  void initEvolgrid(const int& nx, const double& xmin);
+  void initEvolgrid(const int nx, const double xmin, const double a_factor = 6);
   void initPDF(std::string const& setname, int const& i);
 
   // APFEL PDF access functions
@@ -119,6 +119,7 @@ namespace QCD
   double disobs(std::string const& obs, double const& x, double const& Q, double const& y);
 
   // mode setters
+  void setDISmode (bool const& mode);
   void setFTDYmode(bool const& mode);
-  void setSIAmode(bool const& mode);
+  void setSIAmode (bool const& mode);
 }
diff --git a/src/apfelcomb/fk_xgrid.h b/src/apfelcomb/fk_xgrid.h
index d5a1620..f910ae0 100644
--- a/src/apfelcomb/fk_xgrid.h
+++ b/src/apfelcomb/fk_xgrid.h
@@ -1,6 +1,7 @@
 /*
  *  fk_xgrid.h
  *  PDF x grid distributions
+ *  Code adapted from applgrid 1.4.70
  * *  nph 09/14
  */
 
@@ -12,33 +13,45 @@
 
 namespace XGrid
 {
-  static double m_transvar = 6.0;
-
-  static double appl_fy(double x) { return -std::log(x)+m_transvar*(1-x); }
-  static double appl_fx(double y) {
-    // use Newton-Raphson: y = ln(1/x)
-    // solve   y - yp - a(1 - exp(-yp)) = 0
-    // deriv:  - 1 -a exp(-yp)
-    
-    if ( m_transvar==0 )  return std::exp(-y);
-    
-    const double eps  = 1e-12;  // our accuracy goal
-    const int    imax = 100;    // for safety (avoid infinite loops)
-    
-    double yp = y;
-    double x, delta, deriv;
-    for ( int iter=imax ; iter-- ; ) {
-      x = std::exp(-yp);
-      delta = y - yp - m_transvar*(1-x);
-      if ( std::fabs(delta)<eps ) return x; // we have found good solution
-      deriv = -1 - m_transvar*x;
-      yp  -= delta/deriv;
+  // Generates an x-grid according to the APPLgrid prescription.
+  // Points are distributed linearly in y(x) = ln(1/x) + a(1-x).
+  // Parameter `a` controls point density at high-x.
+  class Generator
+  {
+    public:
+      Generator(const double a = 6):
+      m_transvar(a){}
+
+      double appl_fy(double x) const { return -std::log(x)+m_transvar*(1-x); }
+      double appl_fx(double y) const {
+        // use Newton-Raphson: y = ln(1/x)
+        // solve   y - yp - a(1 - exp(-yp)) = 0
+        // deriv:  - 1 -a exp(-yp)
+
+        if ( m_transvar==0 )  return std::exp(-y);
+
+        const double eps  = 1e-12;  // our accuracy goal
+        const int    imax = 100;    // for safety (avoid infinite loops)
+
+        double yp = y;
+        double x, delta, deriv;
+        for ( int iter=imax ; iter-- ; ) {
+          x = std::exp(-yp);
+          delta = y - yp - m_transvar*(1-x);
+          if ( std::fabs(delta)<eps ) return x; // we have found good solution
+          deriv = -1 - m_transvar*x;
+          yp  -= delta/deriv;
+        }
+        // exceeded maximum iterations
+        std::cerr << "_fx2() iteration limit reached y=" << y << std::endl;
+        return std::exp(-yp);
     }
-    // exceeded maximum iterations
-    std::cerr << "_fx2() iteration limit reached y=" << y << std::endl;
-    std::cout << "_fx2() iteration limit reached y=" << y << std::endl;
-    return std::exp(-yp);
-}
+
+    private:
+      const double m_transvar;
+
+  };
+
 
 }
 
diff --git a/src/core/Makefile b/src/core/Makefile
index bd4198c..e89cd73 100644
--- a/src/core/Makefile
+++ b/src/core/Makefile
@@ -6,7 +6,7 @@ CXXFLAGS = 	$(PRJCXXFLAGS) -I ..
 LDFLAGS	 = 	$(PRJLDFLAGS)
 
 CORELIB = ../libac_core.a
-COREOBJ = fk_appl.o fk_dis.o  fk_qcd.o fk_utils.o fk_grids.o fk_ftdy.o
+COREOBJ = fk_appl.o fk_dis.o  fk_qcd.o fk_utils.o fk_grids.o fk_ftdy.o fk_evol.o
 all : $(CORELIB)
 
 $(CORELIB) : $(COREOBJ)
diff --git a/src/core/fk_evol.cc b/src/core/fk_evol.cc
new file mode 100644
index 0000000..7117ac3
--- /dev/null
+++ b/src/core/fk_evol.cc
@@ -0,0 +1,111 @@
+#include "apfelcomb/fk_evol.h"
+#include <array>
+
+
+namespace EVL
+{
+    EvolutionData::EvolutionData( dataInfoRaw const& di,
+            vector<double> const& xgrid,
+            vector<double> const& q2grid ):
+        CommonData(di)
+    {
+        const int nx = xgrid.size();
+        const int nq = q2grid.size();
+        const int nf = 14;
+
+        if (nx*nq*nf != fNData)
+            throw std::runtime_error("Mismatch in nData for Evolution grid kinematics.");
+
+        // Fill data kinematics
+        int iDat = 0;
+        for (auto q2 : q2grid)
+            for (auto x : xgrid)
+                for (int f=0; f < nf; f++)
+                {
+                    fData[iDat] = 0;
+                    fStat[iDat] = 0;
+                    fKin2[iDat] = x;
+                    fKin2[iDat] = q2;
+                    fKin3[iDat] = f;
+                    iDat++;
+                }
+    }
+
+    // allocate grid in x
+    vector<double> generate_xgrid( const int nx,
+            const double xmin,
+            const double xmed,
+            const double xmax)
+    {
+        vector<double> xg;
+        const int nxm = nx/2;
+        for (int ix = 1; ix <= nx; ix++)
+        {
+            if (ix <= nxm)
+                xg.push_back(xmin*pow(xmed/xmin,2.0*(ix-1.0)/(nx-1.0)));
+            else
+                xg.push_back(xmed+(xmax-xmed)*((ix-nxm-1.0)/(nx-nxm-1.0)));
+        }
+        return xg;
+    };
+
+    // Generates a list of nq2 Q2 points between q2min and q2max.
+    // Distributed as linear in tau = log( log( Q2 / lambda2 ) )
+    vector<double> generate_q2grid( const int nq2,
+            const double q2min,
+            const double q2max,
+            const double lambda2
+            )
+    {
+        vector<double> q2grid;
+        const double tau_min   = log( log( q2min / lambda2 ));
+        const double tau_max   = log( log( q2max / lambda2 ));
+        const double delta_tau = (tau_max - tau_min) /( (double) nq2 - 1);
+
+        for (int iq2 = 0; iq2 < nq2; iq2++)
+            q2grid.push_back(lambda2 * exp(exp( tau_min + iq2*delta_tau)));
+        return q2grid;
+    }
+
+    // Compute Q2 grid
+    // This is suitable for PDFs not AlphaS (see maxNF used below)
+    vector< vector<double> > generate_q2subgrids( QCD::qcd_param par, const int nq2,
+            const double q2min, const double q2max)
+    {
+        const double eps     = 1E-4;
+        const double lambda2 = 0.0625e0;
+        const int nfpdf = atoi(par.thMap["MaxNfPdf"].c_str());
+
+        const std::array<double, 3> hqth = {
+            pow(atof(par.thMap["Qmc"].c_str())*atof(par.thMap["kcThr"].c_str()),2),
+            pow(atof(par.thMap["Qmb"].c_str())*atof(par.thMap["kbThr"].c_str()),2),
+            pow(atof(par.thMap["Qmt"].c_str())*atof(par.thMap["ktThr"].c_str()),2),
+        };
+
+        // Determine subgrid edges
+        vector<double> subgrid_edges={q2min};
+        for (int i=0; i<(nfpdf-3); i++)
+            if (hqth[i] > q2min) subgrid_edges.push_back(hqth[i]);
+        subgrid_edges.push_back(q2max);
+
+        // Determine point distribution across subgrids and generate subgrids
+        const std::vector< double > point_dist = generate_q2grid(nq2, q2min, q2max, lambda2);
+        std::vector< vector< double > > q2grid_subgrids;
+        for (int i=0; i<subgrid_edges.size()-1; i++)
+        {
+            const double min_edge = subgrid_edges[i];
+            const double max_edge = subgrid_edges[i+1];
+            auto in_sub = [=](double q2pt){return (q2pt - min_edge) > -eps && (q2pt - max_edge) < eps;};
+            const int npoints_sub = std::count_if(point_dist.begin(), point_dist.end(), in_sub);
+
+            if (npoints_sub < 2)
+                throw std::runtime_error("Too few points to sufficiently populate subgrids. More Q points required");
+
+            const vector< double > subgrid = generate_q2grid(npoints_sub, min_edge, max_edge, lambda2);
+            q2grid_subgrids.push_back(subgrid);
+        }
+
+        return q2grid_subgrids;
+    }
+}
+
diff --git a/src/core/fk_qcd.cc b/src/core/fk_qcd.cc
index ca98905..ca463a5 100644
--- a/src/core/fk_qcd.cc
+++ b/src/core/fk_qcd.cc
@@ -25,21 +25,24 @@ using NNPDF::FKHeader;
 
 namespace QCD
 {
+  static double Q0 = 0; // Initial scale  Q0
+  static double QM = 0; // Maximum scale  QM
+  static double QC = 0; // Cached scale   QC
 
-  // Set static DIS mode
-  static bool FTDY_mode = false;
-  static bool SIA_mode = false;
+  // Required if hadronic scale variations are being performed
   static bool HOPPET_COMPATIBILITY = false;
 
   static std::string FKObs; //!< Cached FK observable
 
-  static double Q0 = 0; // Initial scale  Q0
-  static double QM = 0; // Maximum scale  QM
-  static double QC = 0; // Cached scale   QC
+  // Process specific modes
+  static bool DIS_mode  = false;
+  static bool FTDY_mode = false;
+  static bool SIA_mode  = false;
 
   // Set modes
+  void setDISmode (bool const& mode) {DIS_mode =mode;};
   void setFTDYmode(bool const& mode) {FTDY_mode=mode;};
-  void setSIAmode(bool const& mode) {SIA_mode=mode;};
+  void setSIAmode (bool const& mode) {SIA_mode =mode;};
 
   // *********************** BASIS ROTATION *****************************
 
@@ -139,7 +142,7 @@ namespace QCD
     FK.AddTag(FKHeader::VERSIONS, "GenTime", ctime(&start_time));
   }
 
-  // *********************** EVOLUTON FUNCTIONS *****************************
+  // *********************** EVOLUTION FUNCTIONS *****************************
 
   // Initialise QCD according to parameters
   void initQCD(qcd_param& par, const bool& positivity, const double& Q2max)
@@ -178,7 +181,10 @@ namespace QCD
     }
 
     // Start APFEL
-    APFEL::InitializeAPFEL_DIS();
+    if (DIS_mode)
+        APFEL::InitializeAPFEL_DIS();
+    else
+        APFEL::InitializeAPFEL();
 
     return;
   }
@@ -186,37 +192,41 @@ namespace QCD
   // Initialise APFEL for evolution factors
   // Note here nx is the number of x-points to be output to the Fk table
   // Therefore it doesn't include x=1. This is added manually in this function
-  void initEvolgrid(int const& nx, double const& xmin)
+  // The arguments are the number of desired x-points, the minimum x-value and
+  // the `a` variable used in the grid point distribution (see XGrid::Generator)
+  void initEvolgrid(const int nx, const double xmin, const double a_factor)
   {
     // Reset cache
     QC = 0;
     double* xg = new double[nx+1];
     std::cout << " Initialising  "<< nx << " points starting from x = " << xmin <<std::endl;
 
+    // Initialise x-grid generator
+    const XGrid::Generator xgen(a_factor);
 
     // Special requirements for FTDY in APFEL
     if (FTDY_mode)
     {
       // Requires two x-grid points below x-min
-      xg[0] = 0.9*xmin;
+      xg[0] = 0.90*xmin;
       xg[1] = 0.95*xmin;
 
-      const double ymin = XGrid::appl_fy(xmin);
-      const double ymax = XGrid::appl_fy(1.0);
+      const double ymin = xgen.appl_fy(xmin);
+      const double ymax = xgen.appl_fy(1.0);
 
       // Populate grid
       for (int i=2; i<=nx; i++)
-        xg[i] = XGrid::appl_fx(ymin + ((ymax-ymin)/((double) nx))*(i-2));
+        xg[i] = xgen.appl_fx(ymin + ((ymax-ymin)/((double) nx))*(i-2));
 
     }
     else
     { // Normal mode
-      const double ymin = XGrid::appl_fy(0.99*xmin);
-      const double ymax = XGrid::appl_fy(1.0);
+      const double ymin = xgen.appl_fy(0.99*xmin);
+      const double ymax = xgen.appl_fy(1.0);
 
       // Populate grid
       for (int i=0; i<=nx; i++)
-        xg[i] = XGrid::appl_fx(ymin + ((ymax-ymin)/((double) nx))*i);
+        xg[i] = xgen.appl_fx(ymin + ((ymax-ymin)/((double) nx))*i);
     }
 
     // Just in case of numerical trouble
@@ -238,7 +248,10 @@ namespace QCD
     APFEL::EnableWelcomeMessage(false);
 
     // Start APFEL
-    APFEL::InitializeAPFEL_DIS();
+    if (DIS_mode)
+        APFEL::InitializeAPFEL_DIS();
+    else
+        APFEL::InitializeAPFEL();
 
     // Needed to join the grids
     APFEL::EvolveAPFEL(Q0,Q0);
diff --git a/src/core/fk_utils.cc b/src/core/fk_utils.cc
index a990e6c..5c3812b 100644
--- a/src/core/fk_utils.cc
+++ b/src/core/fk_utils.cc
@@ -121,6 +121,7 @@ std::string setupDir(int const& theoryID)
   mkdir(resultsPath().c_str(),0777);
   mkdir(theoryDir.str().c_str(),0777);
   mkdir((theoryDir.str() + "subgrids/").c_str(),0777);
+  mkdir((theoryDir.str() + "evolution/").c_str(),0777);
   mkdir((theoryDir.str() + "fastkernel/").c_str(),0777);
   return theoryDir.str();
 }