diff --git a/bin/submit.py b/bin/submit.py
index 5e524d52..b11f8c4a 100755
--- a/bin/submit.py
+++ b/bin/submit.py
@@ -5,36 +5,41 @@
 from snek5000.clusters import Cluster
 
 cluster = Cluster()
-sub_dir = "november"
+sub_dir = "test-abl-cli"
 base_name_run = "hi"
-sub_command = "launch"
-#  sub_command = "launch compile"
-# sub_command = "launch release"; cluster.cmd_run = "echo"
+# sub_command = "launch"
+# sub_command = "launch compile"
+# sub_command = "launch release"
+sub_command = "debug"
+cluster.cmd_run = "echo"
 # sub_command = "debug"
 # sub_command = "show box"
-dry_run = False
+dry_run = not False
 
 for (
-    mesh_nb_nodes_walltime,
+    mesh_nb_nodes,
+    walltime,
     filter_weight,
     filter_cutoff,
     filter_temporal,
     sgs_model,
     sgs_boundary,
     z_wall,
+    z_rough,
 ) in itertools.product(
-    zip([111], [1] * 1, [f"{days}-00:00:00" for days in (7,) * 1]),
-    [0.05],
+    itertools.zip_longest([11, 21], [1], fillvalue=1),
+    ["7-00:00:00"],
+    [0.05, 12],
     [0.75],
     # [0.], [1.],
     [False],
-    ["vreman", "constant", "shear_imp"],
+    ["constant"],  # "vreman"], "shear_imp"],
     # "dynamic"],
     [False],
     [0.1],
+    [0.1],
 ):
-    mesh, nb_nodes, walltime = mesh_nb_nodes_walltime
-
+    mesh, nb_nodes = mesh_nb_nodes
 
     name_run = (
         f"{base_name_run}-{sgs_model}-ft{int(filter_temporal)}-sb{int(sgs_boundary)}"
@@ -44,11 +49,11 @@
         continue
 
     cmd = (
-        f"\n{sys.executable} ./simul.py "
+        f"\n{sys.executable} -m abl.cli "
         f"-d {sub_dir} -m {mesh} -n {name_run} -o {nb_nodes} -w {walltime} "
         f"-fw {filter_weight} -fc {filter_cutoff} -ft {filter_temporal} "
         f"-s {sgs_model} -sb {sgs_boundary} "
-        f"-zw {z_wall} "
+        f"-zw {z_wall} -z0 {z_rough} "
         f"{sub_command}"
     )
     if dry_run:
diff --git a/pyproject.toml b/pyproject.toml
index b01c5ed0..0706ef44 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,7 +1,9 @@
 [build-system]
 requires = ["setuptools>=42.0.0", "wheel", "setuptools_scm[toml]>=3.4.2"]
 
-
 [tool.setuptools_scm]
 write_to = "src/abl/_version.py"
 write_to_template = "__version__ = \"{version}\"\n"
+
+[tool.pytest.ini_options]
+addopts = "--instafail --cov=abl --cov=tests"
diff --git a/setup.cfg b/setup.cfg
index c8b9118b..8db6ac8f 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -25,6 +25,13 @@ install_requires =
 setup_requires =
     setuptools_scm
 
+[options.entry_points]
+console_scripts =
+    abl = abl.cli:main
+
+snek5000.solvers =
+    abl = abl.solver
+
 [options.extras_require]
 docs =
     sphinx
@@ -35,6 +42,7 @@ docs =
 
 tests =
     pytest
+    pytest-instafail
     pytest-xdist
     nox
 
@@ -59,7 +67,7 @@ dev =
 where=src
 
 [options.package_data]
-abl = templates/*.j2, bc/*, sgs/*, toolbox/*, *.usr, *.par, *.box, SIZE, Snakefile, etc/*.yml
+abl = templates/*.j2, bc/*, sgs/*, toolbox/*, *.usr, SIZE, Snakefile, etc/*.yml
 
 [flake8]
 ignore = E501,W503,W505
diff --git a/src/abl/SIZE b/src/abl/SIZE
deleted file mode 100644
index 1cb2115b..00000000
--- a/src/abl/SIZE
+++ /dev/null
@@ -1,42 +0,0 @@
-c
-c     Include file to dimension static arrays
-c     and to set some hardwired run-time parameters
-c
-      integer ldim,lx1,lxd,lx2,lx1m,lelg,lelt,lpmin,lpmax,ldimt
-      integer lpelt,lbelt,toteq,lcvelt
-      integer lelx,lely,lelz,mxprev,lgmres,lorder,lhis
-      integer maxobj,lpert,nsessmax,lxo
-      integer lfdm, ldimt_proj
-
-      ! BASIC
-      parameter (ldim=3)               ! domain dimension (2 or 3)
-      parameter (lx1=8)                ! GLL points per element along each direction
-      parameter (lxd=12)               ! GL  points for over-integration (dealiasing)
-      parameter (lx2=lx1-2)            ! GLL points for pressure (lx1 or lx1-2)
-
-      parameter (lelg=12*24*12+10)     ! max total number of elements
-      parameter (lpmin=12)             ! min MPI ranks
-      parameter (lpmax=512)            ! max MPI ranks
-      parameter (ldimt=1)              ! max auxiliary fields (temperature + scalars)
-
-      ! OPTIONAL
-      parameter (ldimt_proj=1)         ! max auxiliary fields residual projection
-      parameter (lhis=1000)            ! max history/monitoring points
-      parameter (maxobj=4)             ! max number of objects
-      parameter (lpert=1)              ! max number of perturbations
-      parameter (toteq=1)              ! max number of conserved scalars in CMT
-      parameter (nsessmax=1)           ! max sessions to NEKNEK
-      parameter (lxo=lx1)              ! max GLL points on output (lxo>=lx1)
-      parameter (mxprev=20,lgmres=30)  ! max dim of projection & Krylov space
-      parameter (lorder=3)             ! max order in time
-      parameter (lx1m=lx1)             ! GLL points mesh solver
-      parameter (lfdm=0)               ! set to 1 for fast diagonalization method
-      parameter (lelx=1,lely=24,lelz=1) ! global tensor mesh dimensions
-
-      parameter (lelt=500)             ! max number of local elements per MPI rank
-      parameter (lbelt=1)              ! set to lelt for mhd
-      parameter (lpelt=1)              ! set to lelt for linear stability
-      parameter (lcvelt=1)             ! set to lelt for cvode
-
-      ! INTERNALS
-      include 'SIZE.inc'
diff --git a/src/abl/Snakefile b/src/abl/Snakefile
index f6b3aa9c..895e5375 100644
--- a/src/abl/Snakefile
+++ b/src/abl/Snakefile
@@ -179,7 +179,7 @@ rule archive:
         ),
         rest=[
             "SESSION.NAME",
-            "params.xml",
+            "params_simul.xml",
             "SIZE",
             *iglob(f"rs6{CASE}0.f*"),
             f"{CASE}.re2",
diff --git a/src/abl/abl.box b/src/abl/abl.box
deleted file mode 100644
index 5bde7670..00000000
--- a/src/abl/abl.box
+++ /dev/null
@@ -1,13 +0,0 @@
--3                     spatial dimension (will create box.re2)
-1                      number of fields
-#
-#    comments: 
-#
-#========================================================
-#
-Box                                       hillp
--12 -24 -12                                    Nelx  Nely
-0.0 1280 1                                x0 x1 ratio
-0.0 1500 1
-0.0 1280 1
-P  ,P  ,sh ,SYM,P  ,P                    BC's:  (cbx0, cbx1, cby0, cby1)
diff --git a/src/abl/abl.par b/src/abl/abl.par
deleted file mode 100644
index ba17f372..00000000
--- a/src/abl/abl.par
+++ /dev/null
@@ -1,59 +0,0 @@
-#
-# nek parameter file
-#
-[GENERAL]
-# startFrom = abl0.f00000
-stopAt = numSteps
-numSteps  = 10
-# endTime = 10
-
-variableDt=yes
-
-targetCfl=0.4
-
-writeControl = timeStep
-writeInterval = 200000
-
-filtering=HPFRT
-filterCutoffRatio=0.75
-filterWeight=0.05
-
-userParam03 = 5.0  # Geostrophic velocity
-userParam04 = -0.00014  # Coriolis frequency
-userParam05 = 1280
-userParam06 = 1500
-userParam07 = 1280
-
-[PROBLEMTYPE]
-equation = incompNS
-stressFormulation=yes
-variableProperties = yes
-
-[PRESSURE]
-residualTol = 1e-5
-residualProj = yes
-
-[VELOCITY]
-residualTol = 1e-8
-density = 1
-viscosity = -1e10
-
-# KTH toolbox
-#
-[_RUNPAR]               # Runtime parameter section for rprm module
-parfWrite            = no                     # Do we write runtime parameter file
-parfName             = outparfile             # Runtime parameter file name for output (without .par)
-#
-[_MONITOR]              # Runtime parameter section for monitor module
-logLevel             = 4                      # Logging threshold for toolboxes
-wallTime             = 23:45                  # Simulation wall time
-#
-[_CHKPOINT]             # Runtime parameter section for checkpoint module
-readChkpt            = no                     # Restart from checkpoint
-chkpFnumber          = 1                      # Restart file number
-chkpInterval         = 250                    # Checkpiont saving frequency (number of time steps)
-#
-[_STAT]                # Runtime parameter section for statistics module
-avStep               = 5
-ioStep               = 10
-# vim: set ft=cfg
diff --git a/src/abl/abl.usr b/src/abl/abl.usr
index 016e70a6..063f6851 100644
--- a/src/abl/abl.usr
+++ b/src/abl/abl.usr
@@ -62,24 +62,25 @@ c
       implicit none
 
       include 'SIZE'
-      include 'INPUT'           ! uparam
       include 'PARALLEL'        ! GLLEL
       include 'NEKUSE'            ! ffx, ffy, ffz, ux, uz
+      include 'FLOWPHYS'        ! fp_corio_on, fp_corio_freq, fp_u_geo
 
       ! argument list
       integer ix,iy,iz,eg
 
-      ! local variables
-      real U_GEO, f_corio
 
 c     e = gllel(eg)
 
-      U_GEO = uparam(3)
-      f_corio = uparam(4)
-
-      ffx = f_corio * uz
-      ffy = 0.0
-      ffz = -f_corio * (ux - U_GEO)
+      if (fp_corio_on) then
+         ffx = fp_corio_freq * uz
+         ffy = 0.0
+         ffz = -fp_corio_freq * (ux - fp_u_geo)
+      else
+         ffx = 0.0
+         ffy = 0.0
+         ffz = 0.0
+      endif
 
       return
       end
@@ -113,6 +114,7 @@ c> @callgraph
       include 'GEOM'  ! ym1
       include 'INPUT'  ! uparam
       include 'WMLES'  ! KAPPA, PI_, y0
+      include 'FLOWPHYS'  ! fp_u_geo
 
       integer ix,iy,iz,eg
       real Ly
@@ -126,7 +128,7 @@ c      uz   = 0.0
       temp = 0.0
 
 c     Geostrophic velocity
-      U_GEO = uparam(3)
+      U_GEO = fp_u_geo
 c     Get mesh lengths from userParam0{5,6,7}
       XLEN = uparam(5)
       YLEN = uparam(6)
@@ -320,6 +322,14 @@ c> @callgraph
       implicit none
 
       include 'SIZE'
+      include 'INPUT'  ! param
+      include 'FLOWPHYS'  ! fp_corio_on, fp_u_geo
+
+      if (.not. fp_corio_on) then
+         param(54) = -1  ! use >0 for const flowrate or <0 bulk vel
+                         ! flow direction is given by (1=x, 2=y, 3=z)
+         param(55) = fp_u_geo ! flowrate/bulk-velocity
+      endif
 
       return
       end
@@ -341,6 +351,7 @@ c> @{
       call chkpt_register
       call stat_register
       call wmles_register
+      call fp_register
 
       return
       end subroutine
@@ -357,6 +368,7 @@ c> @{
       call chkpt_init
       call stat_init
       call wmles_init
+      call fp_init
 
       return
       end subroutine
diff --git a/bin/simul.py b/src/abl/cli.py
similarity index 84%
rename from bin/simul.py
rename to src/abl/cli.py
index deece7e8..63b3d9e2 100755
--- a/bin/simul.py
+++ b/src/abl/cli.py
@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 """Make a simulation of with solver abl."""
-import sys
+from math import pi
 
 import click
 from abl.output import avail_boundary_conds, avail_sgs_models
@@ -17,11 +17,7 @@
 @click.option("-zw", "--z-wall", default=0.0, type=float, help="wall position")
 @click.option("-z0", "--z-rough", default=0.1, type=float, help="roughness parameter")
 @click.option(
-    "-fw",
-    "--filter-weight",
-    default=0.05,
-    type=float,
-    help="filter weight parameter",
+    "-fw", "--filter-weight", default=0.05, type=float, help="filter weight parameter",
 )
 @click.option(
     "-fc", "--filter-cutoff", default=0.75, type=float, help="filter cutoff ratio"
@@ -100,7 +96,7 @@ def cli(
     elif M in (11, 111):
         oper.ny = 32
         oper.coords_y = (
-            "0.1 5.4 11.6 18.7 26.8 36 46.5 58.5 72 87.5 105 125 147 173 201 "
+            f"{z_wall} 5.4 11.6 18.7 26.8 36 46.5 58.5 72 87.5 105 125 147 173 201 "
             "233 269 310 354 404 458 518 583 654 731 812 900 991 1088 1187 "
             "1290 1394 1500"
         )
@@ -113,20 +109,49 @@ def cli(
         oper.ny = 20
         oper.nz = 4
         oper.coords_y = (
-            "0.1 34.6 72.8 115. 161. 211. 266. 326. 390. 459. 534. 613. 697. "
+            f"{z_wall} 34.6 72.8 115. 161. 211. 266. 326. 390. 459. 534. 613. 697. "
             "786. 879. 976. 1076. 1180. 1285. 1392. 1500."
         )
+    elif M == 21:
+        # Similar to M 11 - scaled by 1500
+        oper.nx = 4
+        oper.ny = 32
+        oper.nz = 4
+        oper.coords_y = (
+            f"{z_wall} 0.0036 0.0077 0.0125 0.0179 0.0240 0.0310 0.0390 0.0480 "
+            "0.0583 0.0700 0.0833 0.0980 0.1153 0.1340 0.1553 0.1793 0.2067 "
+            "0.2360 0.2693 0.3053 0.3453 0.3887 0.4360 0.4873 0.5413 0.6000 "
+            "0.6607 0.7253 0.7913 0.8600 0.9293 1.0000"
+        )
+    elif M in (22, 122):
+        # Similar to M 12 - scaled by 1500
+        oper.ny = 20
+        oper.coords_y = (
+            f"{z_wall} 0.0231 0.0485 0.0767 0.1073 0.1407 0.1773 0.2173 0.2600 "
+            "0.3060 0.3560 0.4087 0.4647 0.5240 0.5860 0.6507 0.7173 0.7867 "
+            "0.8567 0.9280 1.0000"
+        )
+        if M == 22:
+            oper.nx = oper.nz = 4
+        elif M == 122:
+            oper.nx = oper.nz = 8
 
-    if M < 10:
+    if M <= 9:
         oper.origin_y = z_wall
         oper.Lx = 1280
         oper.Ly = 1500
         oper.Lz = 1280
-    elif M < 20:
+    elif M % 100 <= 19:
         oper.origin_y = float(oper.coords_y.split()[0])
         oper.Lx = 640
         oper.Ly = 1500
         oper.Lz = 640
+    elif M % 100 <= 29:
+        # Chatterjee & Peet:
+        oper.origin_y = float(oper.coords_y.split()[0])
+        oper.Lx = round(2 * pi, 4)
+        oper.Ly = 1.0
+        oper.Lz = round(pi, 4)
 
     oper.boundary = "P P sh SYM P P".split()
 
@@ -175,6 +200,7 @@ def cli(
     # Original value:
     # general.target_cfl = 0.8
     #  general.target_cfl = 0.3
+    # TODO: try BDF3
     #  general.time_stepper = "BDF2"
     general.write_control = (
         "timeStep"
@@ -190,8 +216,6 @@ def cli(
     general.filter_weight = filter_weight
     general.filter_cutoff_ratio = filter_cutoff
     general.user_params = {
-        3: 5.0,  # Geostrophic velocity
-        4: -1.4e-4,  # Coriolis frequency at 73 S
         5: oper.Lx,
         6: oper.Ly,
         7: oper.Lz,
@@ -217,7 +241,7 @@ def cli(
     # NOTE: reducing pressure residual tolerance affects velocity divergence
     # TODO: check if w -> O(pressure.residual_tol)
     pressure.residual_tol = 1e-5
-    #  velocity.residual_tol = 1e-8
+    velocity.residual_tol = 1e-8
     #
     if params.output.boundary_cond == "noslip":
         reynolds_number = 1e4
@@ -242,9 +266,15 @@ def cli(
     wmles.bc_z0 = z_rough
     wmles.sgs_delta_max = True
     # wmles.sgs_npow = 3.0
-    wmles.sgs_c0 = 0.18
+    wmles.sgs_c0 = 0.19
     wmles.sgs_bc = sgs_boundary
 
+    # Flow phys parameters
+    # ====================
+    fp = params.nek.flow_phys
+    fp.corio_on = False
+    fp.u_geo = 1.0
+
     # Fluidsim parameters
     # ===================
     params.short_name_type_run = name_run
@@ -265,7 +295,6 @@ def launch(ctx, rule):
     logger.info("Initializing simulation launch...")
 
     sim = Simul(ctx.obj["params"])
-    sim.sanity_check()
     assert sim.make.exec([rule], scheduler="greedy")
     if rule == "release":
         import shutil
@@ -302,7 +331,6 @@ def debug(ctx, rule):
     logger.info("Initializing simulation debug...")
 
     sim = Simul(params)
-    sim.sanity_check()
     logger.info("Executing simulation...")
     sim.make.exec([rule])
     logger.info("Finished simulation...")
@@ -351,6 +379,9 @@ def show(ctx, file):
         )
 
 
-if __name__ == "__main__":
-    print(sys.executable)
+def main():
     cli(obj={})
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/abl/etc/archmage.yml b/src/abl/etc/archmage.yml
index ff8acf47..28cc1fd3 100644
--- a/src/abl/etc/archmage.yml
+++ b/src/abl/etc/archmage.yml
@@ -1,5 +1,5 @@
-CC: gcc-8
-FC: gfortran-8
+CC: gcc
+FC: gfortran
 MPICC: mpicc
 MPIFC: mpif77
 MPIEXEC: "mpiexec --oversubscribe"
diff --git a/src/abl/etc/ponyville.yml b/src/abl/etc/ponyville.yml
new file mode 100644
index 00000000..97fc8356
--- /dev/null
+++ b/src/abl/etc/ponyville.yml
@@ -0,0 +1,8 @@
+CC: gcc-7
+FC: gfortran-7
+MPICC: mpicc
+MPIFC: mpif77
+MPIEXEC: "mpiexec"
+CFLAGS: "-march=native"
+FFLAGS: "-march=native -mcmodel=medium -std=legacy"
+
diff --git a/src/abl/flow_phys/FLOWPHYS b/src/abl/flow_phys/FLOWPHYS
new file mode 100644
index 00000000..1af07a13
--- /dev/null
+++ b/src/abl/flow_phys/FLOWPHYS
@@ -0,0 +1,21 @@
+c
+c     Storage for flow / physics parameters
+c
+      ! initialisation flag
+      logical fp_ifinit
+
+      ! parameter section
+      logical fp_corio_on      !< @var turn coriolis forcing on / off
+      real fp_corio_freq       !< @var coriolis frequency
+      real fp_u_geo            !< @var geostrophic velocity
+
+      ! parameter internals
+      character*(*) fp_sec_name
+      parameter (fp_sec_name="FLOW_PHYS")
+      integer fp_id, fp_sec_id, fp_par_id(3)
+
+      ! common block for storage
+      common /fp_log/ fp_corio_on
+
+      common /fp_real/ fp_corio_freq, fp_u_geo
+c vim: set ft=fortran
diff --git a/src/abl/flow_phys/__init__.py b/src/abl/flow_phys/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/src/abl/flow_phys/flow_phys_init.f b/src/abl/flow_phys/flow_phys_init.f
new file mode 100644
index 00000000..eb6ddd8e
--- /dev/null
+++ b/src/abl/flow_phys/flow_phys_init.f
@@ -0,0 +1,121 @@
+!> @file init.f
+!! @ingroup fp
+!! @brief Flow and physics (FLOWPHYS) module
+!! @details Register and initialize parameters which control flow physics
+!! @{
+!=======================================================================
+!> Register FLOWPHYS module
+!! @note This routine should be called in frame_usr_register
+      subroutine fp_register()
+      implicit none
+
+      include 'SIZE'
+      include 'INPUT'
+      include 'FRAMELP'
+      include 'FLOWPHYS'
+
+      ! local variables
+      integer lpmid
+!-----------------------------------------------------------------------
+      ! check if the current module was already registered
+      call mntr_mod_is_name_reg(lpmid, fp_sec_name)
+      if (lpmid.gt.0) then
+         call mntr_warn(
+     $      lpmid,
+     $     'module ['//trim(fp_sec_name)//'] already registered')
+         return
+      endif
+
+      ! find parent module
+      call mntr_mod_is_name_reg(lpmid,'FRAME')
+      if (lpmid.le.0) then
+         lpmid = 1
+         call mntr_abort(lpmid,
+     $        'parent module ['//'FRAME'//'] not registered')
+      endif
+
+      ! register module
+      call mntr_mod_reg(fp_id, lpmid, fp_sec_name,
+     $      'Flow Physics parameters')
+
+      ! register and set active section
+      call rprm_sec_reg(
+     $    fp_sec_id,fp_id,'_'//adjustl(fp_sec_name),
+     $    'Runtime parameter section for FLOWPHYS module')
+      call rprm_sec_set_act(.true.,fp_sec_id)
+
+      ! register parameters
+      ! subroutine rprm_rp_reg(rpid, mid, pname, pdscr, ptype, ipval, rpval, lpval, cpval)
+      call rprm_rp_reg(
+     $    fp_par_id(1), fp_sec_id,
+     $    'CORIOON', 'Turn coriolis forcing on', rpar_log,
+     $    1, 0.0, .true., ' '
+     $)
+      call rprm_rp_reg(
+     $    fp_par_id(2), fp_sec_id,
+     $    'CORIOFREQ', 'Coriolis frequency', rpar_real,
+     $    1, 1.4e-4, .false., ' '
+     $)
+      call rprm_rp_reg(
+     $    fp_par_id(3), fp_sec_id,
+     $    'UGEO', 'Geostrophic velocity', rpar_real,
+     $    1, 5.0, .false., ' '
+     $)
+
+      ! set initialisation flag
+      fp_ifinit=.false.
+
+      end subroutine
+!=======================================================================
+!> Initiliase FLOWPHYS Module
+!! @note This routine should be called in frame_usr_init
+      subroutine fp_init()
+      implicit none
+
+      include 'SIZE'
+      include 'FRAMELP'
+      include 'FLOWPHYS'
+
+      ! local variables for reading parameters with appropriate types
+      integer itmp
+      real rtmp
+      logical ltmp
+      character*20 ctmp
+!-----------------------------------------------------------------------
+      ! check if the module was already initialised
+      if (fp_ifinit) then
+         call mntr_warn(fp_id,
+     $        'module ['//trim(fp_sec_name)//'] already initiaised.')
+         return
+      endif
+
+      ! get runtime parameters
+      call rprm_rp_get(itmp, rtmp, ltmp, ctmp,
+     $    fp_par_id(1), rpar_log)
+      fp_corio_on = ltmp
+
+      call rprm_rp_get(itmp, rtmp, ltmp, ctmp,
+     $    fp_par_id(2), rpar_real)
+      fp_corio_freq = rtmp
+
+      call rprm_rp_get(itmp, rtmp, ltmp, ctmp,
+     $    fp_par_id(3), rpar_real)
+      fp_u_geo = rtmp
+
+      ! everything is initialised
+      fp_ifinit=.true.
+      end subroutine
+!=======================================================================
+!> Check if module was initialised
+!! @return fp_is_initialised
+      logical function fp_is_initialised()
+      implicit none
+
+      include 'SIZE'
+      include 'FLOWPHYS'
+      fp_is_initialised = fp_ifinit
+
+      return
+      end function
+!=======================================================================
+!! @}
diff --git a/src/abl/io.py b/src/abl/io.py
index a5d2b6cd..998d842a 100644
--- a/src/abl/io.py
+++ b/src/abl/io.py
@@ -2,8 +2,8 @@
 from pathlib import Path
 from tarfile import TarFile
 
-from pymech.dataset import open_mfdataset
 import zstandard as zstd
+from pymech.dataset import open_mfdataset
 
 from .solver import Simul
 
@@ -14,7 +14,7 @@ def open_tar_zst(path_tar_zst, mode="rb"):
     with open(path_tar_zst, mode) as fh:
         if mode == "rb":
             dctx = zstd.ZstdDecompressor()
-            with dctx.stream_reader(fh) as stream: 
+            with dctx.stream_reader(fh) as stream:
                 yield TarFile(fileobj=stream)
         elif mode == "r+b":
             #  FIXME: Writing files into
@@ -27,10 +27,10 @@ def open_tar_zst(path_tar_zst, mode="rb"):
 
 def load_simul(path_dir, prefix="abl", extract_tarballs=True, **kwargs):
     """Load simulation data and parameters.
-    
+
     Parameters
     ----------
-    
+
     path_dir: str
         Path to simulation
 
@@ -39,17 +39,16 @@ def load_simul(path_dir, prefix="abl", extract_tarballs=True, **kwargs):
 
     extract_tarballs: bool
         Extract simulation tarballs
-    
+
     kwargs: dict
         Keyword arguments passed to ``pymech.open_mfdataset``
 
     """
     path_dir = Path(path_dir)
-    params = Simul.load_params_from_file(path_dir / "params.xml")
+    params = Simul.load_params_from_file(path_dir / "params_simul.xml")
     ds = open_mfdataset(path_dir.glob(f"{prefix}0.f*"), **kwargs)
-    
+
     if extract_tarballs:
         raise NotImplementedError("Need to use open_tar_zst")
-                        
+
     return params, ds
-    
\ No newline at end of file
diff --git a/src/abl/makefile_usr.inc b/src/abl/makefile_usr.inc
deleted file mode 100644
index c25cfc24..00000000
--- a/src/abl/makefile_usr.inc
+++ /dev/null
@@ -1,72 +0,0 @@
-
-# vim: set ft=make
-# demo
-$(OBJDIR)/frame.o             :toolbox/frame.f toolbox/FRAMELP
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/mntrlog_block.o     :toolbox/mntrlog_block.f toolbox/MNTRLOGD
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/mntrlog.o           :toolbox/mntrlog.f toolbox/MNTRLOGD
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/mntrtmr_block.o     :toolbox/mntrtmr_block.f toolbox/MNTRLOGD toolbox/MNTRTMRD
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/mntrtmr.o           :toolbox/mntrtmr.f toolbox/MNTRLOGD toolbox/MNTRTMRD toolbox/FRAMELP
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/rprm_block.o        :toolbox/rprm_block.f toolbox/RPRMD
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/rprm.o              :toolbox/rprm.f toolbox/RPRMD toolbox/FRAMELP
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/io_tools_block.o    :toolbox/io_tools_block.f toolbox/IOTOOLD
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/io_tools.o          :toolbox/io_tools.f toolbox/IOTOOLD
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/chkpoint.o          :toolbox/chkpoint.f toolbox/CHKPOINTD
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/chkpt_mstp.o        :toolbox/chkpt_mstp.f toolbox/CHKPTMSTPD toolbox/CHKPOINTD
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/map2D.o             :toolbox/map2D.f toolbox/MAP2D toolbox/FRAMELP
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/stat.o              :toolbox/stat.f toolbox/STATD toolbox/MAP2D toolbox/FRAMELP toolbox/../sgs/SGS
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/stat_IO.o           :toolbox/stat_IO.f toolbox/STATD toolbox/MAP2D toolbox/FRAMELP
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/math_tools.o        :toolbox/math_tools.f
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/utils.o             :sgs/utils.f sgs/SGS
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/wmles_init.o        :sgs/wmles_init.f sgs/WMLES sgs/../toolbox/FRAMELP
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/stat_extras_dummy.o :toolbox/stat_extras_dummy.f
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/smagorinsky.o       :sgs/smagorinsky.f sgs/SGS sgs/WMLES
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/shear_imp_smag.o    :sgs/shear_imp_smag.f sgs/SGS sgs/WMLES
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/dyn_smag.o          :sgs/dyn_smag.f sgs/DYN sgs/SGS sgs/WMLES
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/noslip.o            :bc/noslip.f
-	$(F77) -c $(FL2) -I./ $< -o $@
-
-$(OBJDIR)/moeng.o             :bc/moeng.f bc/../sgs/SGS bc/../sgs/WMLES
-	$(F77) -c $(FL2) -I./ $< -o $@
-
diff --git a/src/abl/output.py b/src/abl/output.py
index 0ea090eb..9bcab212 100644
--- a/src/abl/output.py
+++ b/src/abl/output.py
@@ -1,5 +1,7 @@
 from collections import namedtuple
 
+from abl.templates import box, makefile_usr, size
+from snek5000 import mpi
 from snek5000.output.base import Output as OutputBase
 
 SGS = namedtuple("SGS", ["name", "sources"])
@@ -54,6 +56,9 @@ def makefile_usr_sources(self):
                 ("utils.f", "SGS"),
                 ("wmles_init.f", "WMLES", "../toolbox/FRAMELP"),
             ],
+            "flow_phys": [
+                ("flow_phys_init.f", "FLOWPHYS", "../toolbox/FRAMELP"),
+            ],
             "bc": [],
         }
 
@@ -61,7 +66,7 @@ def makefile_usr_sources(self):
             # Hack to load params from params.xml in current directory
             from abl.solver import Simul
 
-            params = Simul.load_params_from_file(path_xml="params.xml")
+            params = Simul.load_params_from_file(path_xml="params_simul.xml")
         else:
             params = self.sim.params
 
@@ -94,3 +99,17 @@ def fortran_inc_flags(self):
     def _complete_params_with_default(params, info_solver):
         OutputBase._complete_params_with_default(params, info_solver)
         params.output._set_attribs({"sgs_model": "constant", "boundary_cond": "moeng"})
+
+    def post_init(self):
+        params = self.sim.params
+        params.nek.general.user_params[5] = params.oper.Lx
+        params.nek.general.user_params[6] = params.oper.Ly
+        params.nek.general.user_params[7] = params.oper.Lz
+
+        super().post_init()
+
+        # Write additional source files to compile the simulation
+        if mpi.rank == 0 and self._has_to_save and self.sim.params.NEW_DIR_RESULTS:
+            self.write_box(box)
+            self.write_size(size)
+            self.write_makefile_usr(makefile_usr)
diff --git a/src/abl/postproc/__init__.py b/src/abl/postproc/__init__.py
new file mode 100644
index 00000000..8a76d1a7
--- /dev/null
+++ b/src/abl/postproc/__init__.py
@@ -0,0 +1,46 @@
+"""Postprocessing output files
+
+
+**Contents**
+
+.. autosummary::
+   :toctree:
+
+   spatial_means
+   stats
+
+"""
+from inflection import underscore
+
+
+class PostprocABL:
+    @staticmethod
+    def _complete_info_solver(info_solver):
+        """Complete the ParamContainer info_solver."""
+        classes = info_solver.classes.Postproc._set_child("classes")
+
+        classes._set_child(
+            "SpatialMeans",
+            attribs={
+                "module_name": "abl.postproc.spatial_means",
+                "class_name": "SpatialMeansABL",
+            },
+        )
+
+    def __init__(self, sim):
+        self.sim = sim
+
+    def post_init(self, sim):
+        dict_classes = self.sim.info_solver.classes.Postproc.import_classes()
+        for cls_name, Class in dict_classes.items():
+            # only initialize if Class is not the Simul class
+            if not isinstance(self, Class):
+                setattr(self, underscore(cls_name), Class(sim))
+
+    def load(self):
+        """Load results for postprocessing."""
+        dict_classes = self.sim.info_solver.classes.Postproc.import_classes()
+
+        for cls_name in dict_classes:
+            cls = getattr(self, underscore(cls_name))
+            cls.load()
diff --git a/src/abl/postproc/spatial_means.py b/src/abl/postproc/spatial_means.py
new file mode 100644
index 00000000..17e0a218
--- /dev/null
+++ b/src/abl/postproc/spatial_means.py
@@ -0,0 +1,15 @@
+from pathlib import Path
+
+import pandas as pd
+
+
+class SpatialMeansABL:
+    def __init__(self, sim):
+        self.file = Path(sim.path_run / "spatial_means.txt")
+
+    def load(self):
+        df = pd.read_csv(self.file, sep=r"\s+", index_col="it")
+        df["hours"] = df.t / 3600
+        self.df = df
+
+        self.plot = self.df.plot
diff --git a/src/abl/solver.py b/src/abl/solver.py
index 01de7854..7db43890 100644
--- a/src/abl/solver.py
+++ b/src/abl/solver.py
@@ -1,20 +1,11 @@
-from snek5000 import logger, mpi
 from snek5000.info import InfoSolverMake
 from snek5000.solvers.kth import SimulKTH
 
-from .output import OutputABL
-from .templates import box, makefile_usr, size
-
 
 class InfoSolverABL(InfoSolverMake):
-    """Contain the information on a :class:`eturb.solvers.abl.Simul`
+    """Contain the information on a :class:`abl.solver.SimulABL`
     instance.
 
-    .. todo::
-
-        Move Output info to :class:`InfoSolverNek` and only override it in
-        :class:`InfoSolverABL`.
-
     """
 
     def _init_root(self):
@@ -29,11 +20,14 @@ def _init_root(self):
         self.classes.Output.module_name = "abl.output"
         self.classes.Output.class_name = "OutputABL"
 
+        self.classes._set_child(
+            "Postproc",
+            attribs={"module_name": "abl.postproc", "class_name": "PostprocABL"},
+        )
+
 
 class SimulABL(SimulKTH):
-    """A solver which compiles and runs using a Snakefile.
-
-    """
+    """A solver which compiles and runs using a Snakefile."""
 
     InfoSolver = InfoSolverABL
 
@@ -41,7 +35,11 @@ class SimulABL(SimulKTH):
     def _complete_params_with_default(params):
         """Add missing default parameters."""
         params = SimulKTH._complete_params_with_default(params)
-        params.nek.velocity._set_attrib("advection", True)
+
+        params.nek.problemtype.variable_properties = True
+
+        params.nek.velocity._set_attribs({"density": 1.0})
+        params.nek.pressure.residual_proj = True
 
         params.nek._set_child(
             "wmles",
@@ -56,6 +54,11 @@ def _complete_params_with_default(params):
             ),
         )
         params.nek.wmles._set_internal_attr("_enabled", True)
+
+        params.nek._set_child(
+            "flow_phys", dict(corio_on=True, corio_freq=1.4e-4, u_geo=5.0,),
+        )
+        params.nek.flow_phys._set_internal_attr("_enabled", True)
         return params
 
     @classmethod
@@ -66,28 +69,12 @@ def create_default_params(cls):
         """
         params = super().create_default_params()
 
-        # Synchronize baseline parameters as follows:
-        # -----------------------------------------------------------------
-        primary_par_file = OutputABL.get_root() / "abl.par"
-        if mpi.rank == 0:
-            logger.info(f"Reading baseline parameters from {primary_par_file}")
-
-        params.nek._read_par(primary_par_file)
-
         return params
 
     def __init__(self, params):
         super().__init__(params)
-        self.output.write_box(box)
-        self.output.write_size(size)
-        self.output.write_makefile_usr(makefile_usr)
-
-    def sanity_check(self):
-        """Check params for errors"""
-        params = self.params
-        assert params.oper.Lx == params.nek.general.user_params[5]
-        assert params.oper.Ly == params.nek.general.user_params[6]
-        assert params.oper.Lz == params.nek.general.user_params[7]
+
+        self.postproc.post_init(self)
 
 
 Simul = SimulABL
diff --git a/tests/conftest.py b/tests/conftest.py
index f0249382..d15ff82c 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -38,6 +38,9 @@ def sim(request):
 
     params.output.sgs_model = request.param
 
+    reynolds_number = 1e10
+
+    params.nek.velocity.viscosity = -reynolds_number
     params.nek.stat.av_step = 3
     params.nek.stat.io_step = 9
 
diff --git a/tests/test_output.py b/tests/test_output.py
index 9a9c8e03..0b1e65d3 100644
--- a/tests/test_output.py
+++ b/tests/test_output.py
@@ -6,6 +6,10 @@
 def test_copy():
     new_dir = tempfile.mkdtemp(__name__)
     output = Output()
+
+    # Inject name_solver
+    output.name_solver = "abl"
+
     # Fresh copy
     output.copy(new_dir)
     # Redo copy, skip files
diff --git a/tests/test_solver.py b/tests/test_solver.py
index bc089239..2ad9ea40 100644
--- a/tests/test_solver.py
+++ b/tests/test_solver.py
@@ -8,6 +8,11 @@ def test_init():
     params = Simul.create_default_params()
     params.output.sub_directory = "test"
     sim = Simul(params)
+    # Check params for errors
+    params = sim.params
+    assert params.oper.Lx == params.nek.general.user_params[5]
+    assert params.oper.Ly == params.nek.general.user_params[6]
+    assert params.oper.Lz == params.nek.general.user_params[7]
     print(sim.info_solver)