diff --git a/src/3d/epic3d.f90 b/src/3d/epic3d.f90
index d8f23e1a2..170c3ac91 100644
--- a/src/3d/epic3d.f90
+++ b/src/3d/epic3d.f90
@@ -21,7 +21,7 @@ program epic3d
     use field_netcdf, only : field_io_timer
     use field_diagnostics, only : field_stats_timer
     use field_diagnostics_netcdf, only : field_stats_io_timer
-    use inversion_mod, only : vor2vel_timer, db_timer
+    use inversion_mod, only : vor2vel_timer, vtend_timer
     use inversion_utils, only : init_fft
     use parcel_interpl, only : grid2par_timer, par2grid_timer
     use parcel_init, only : init_parcels, init_timer
@@ -72,7 +72,7 @@ subroutine pre_run
             call register_timer('field diagnostics', field_stats_timer)
             call register_timer('field diagnostics I/O', field_stats_io_timer)
             call register_timer('vor2vel', vor2vel_timer)
-            call register_timer('buoyancy derivatives', db_timer)
+            call register_timer('vorticity tendency', vtend_timer)
             call register_timer('parcel push', rk4_timer)
             call register_timer('merge nearest', merge_nearest_timer)
             call register_timer('merge tree resolve', merge_tree_resolve_timer)
diff --git a/src/3d/fields/fields.f90 b/src/3d/fields/fields.f90
index 121560183..b22319b7d 100644
--- a/src/3d/fields/fields.f90
+++ b/src/3d/fields/fields.f90
@@ -17,6 +17,7 @@ module fields
     double precision, allocatable, dimension(:, :, :, :) :: &
         velog,     &   ! velocity vector field (u, v, w)
         vortg,     &   ! vorticity vector field (\omegax, \omegay, \omegaz)
+        vtend,     &   ! vorticity tendency
         velgradg       ! velocity gradient tensor
                        ! ordering: du/dx, du/dy,
                        !                  dv/dy,
@@ -34,8 +35,6 @@ module fields
         dbuoyg,    &   ! dry buoyancy (or liquid-water buoyancy)
 #endif
         tbuoyg,    &   ! buoyancy
-        dbdx,      &   ! buoyancy derivative in x
-        dbdy,      &   ! buoyancy derivative in y
 #ifndef NDEBUG
         sym_volg,  &   ! symmetry volume (debug mode only)
 #endif
@@ -64,9 +63,7 @@ subroutine field_alloc
 
             allocate(vortg(-1:nz+1, 0:ny-1, 0:nx-1, 3))
 
-            allocate(dbdx(-1:nz+1, 0:ny-1, 0:nx-1))
-
-            allocate(dbdy(-1:nz+1, 0:ny-1, 0:nx-1))
+            allocate(vtend(-1:nz+1, 0:ny-1, 0:nx-1, 3))
 
             allocate(tbuoyg(-1:nz+1, 0:ny-1, 0:nx-1))
 
@@ -87,8 +84,7 @@ subroutine field_default
             velgradg = zero
             volg     = zero
             vortg    = zero
-            dbdx     = zero
-            dbdy     = zero
+            vtend    = zero
             tbuoyg   = zero
 #ifndef ENABLE_DRY_MODE
             dbuoyg   = zero
diff --git a/src/3d/inversion/inversion.f90 b/src/3d/inversion/inversion.f90
index 666539a6a..441daa816 100644
--- a/src/3d/inversion/inversion.f90
+++ b/src/3d/inversion/inversion.f90
@@ -1,12 +1,13 @@
 module inversion_mod
     use inversion_utils
     use parameters, only : nx, ny, nz
+    use physics, only : ft_cor, f_cor
     use constants, only : zero, two, f12
     use timer, only : start_timer, stop_timer
     implicit none
 
     integer :: vor2vel_timer,   &
-               db_timer
+               vtend_timer
 
     contains
 
@@ -14,9 +15,10 @@ module inversion_mod
 
         ! Given the vorticity vector field (vortg) in physical space, this
         ! returns the associated velocity field (velog) and the velocity
-        ! gradient tensor (velgradg).
+        ! gradient tensor (velgradg).  Note: the
+        ! vorticity is modified to be solenoidal and spectrally filtered.
         subroutine vor2vel(vortg,  velog,  velgradg)
-            double precision, intent(in)    :: vortg(-1:nz+1, 0:ny-1, 0:nx-1, 3)
+            double precision, intent(inout) :: vortg(-1:nz+1, 0:ny-1, 0:nx-1, 3)
             double precision, intent(out)   :: velog(-1:nz+1, 0:ny-1, 0:nx-1, 3)
             double precision, intent(out)   :: velgradg(-1:nz+1, 0:ny-1, 0:nx-1, 5)
             double precision                :: svelog(0:nz, 0:nx-1, 0:ny-1, 3)
@@ -24,30 +26,93 @@ subroutine vor2vel(vortg,  velog,  velgradg)
                                              , bs(0:nz, 0:nx-1, 0:ny-1) &
                                              , cs(0:nz, 0:nx-1, 0:ny-1)
             double precision                :: ds(0:nz, 0:nx-1, 0:ny-1) &
-                                             , es(0:nz, 0:nx-1, 0:ny-1)
-            double precision                :: dstop(0:nx-1, 0:ny-1), dsbot(0:nx-1, 0:ny-1)
+                                             , es(0:nz, 0:nx-1, 0:ny-1) &
+                                             , fs(0:nz, 0:nx-1, 0:ny-1)
             double precision                :: ubar(0:nz), vbar(0:nz)
             double precision                :: uavg, vavg
             integer                         :: iz
 
             call start_timer(vor2vel_timer)
 
-            ! Copy vorticity to velocity field to perform FFT transforms
-            ! (FFT transforms overwrite the input array)
-            velog = vortg
-
             !------------------------------------------------------------------
-            !Convert vorticity to spectral space as (as, bs, cs): (velog is overwritten in this operation)
-            call fftxyp2s(velog(0:nz, :, :, 1), as)
-            call fftxyp2s(velog(0:nz, :, :, 2), bs)
-            call fftxyp2s(velog(0:nz, :, :, 3), cs)
+            !Convert vorticity to spectral space as (as, bs, cs):
+            call fftxyp2s(vortg(0:nz, :, :, 1), as)
+            call fftxyp2s(vortg(0:nz, :, :, 2), bs)
+            call fftxyp2s(vortg(0:nz, :, :, 3), cs)
+
+            !Add -grad(lambda) where Laplace(lambda) = div(vortg) to
+            !enforce the solenoidal condition on the vorticity field:
+            call diffx(as, ds)
+            call diffy(bs, es)
+
+            !For the vertical parcel vorticity, use 2nd-order
+            !differencing:
+            call diffz(cs, fs)
+
+            !Form div(vortg):
+            !$omp parallel
+            !$omp workshare
+            fs = ds + es + fs
+            !$omp end workshare
+            !$omp end parallel
+
+            !Remove horizontally-averaged part (plays no role):
+            fs(:, 0, 0) = zero
+
+            !Invert Lap(lambda) = div(vortg) assuming dlambda/dz = 0 at the
+            !boundaries (store solution lambda in fs):
+            call lapinv1(fs)
+
+            !Filter lambda:
+            !$omp parallel shared(fs, filt, nz) private(iz)  default(none)
+            !$omp do
+            do iz = 0, nz
+                fs(iz, :, :) = filt * fs(iz, :, :)
+            enddo
+            !$omp end do
+            !$omp end parallel
+
+            !Subtract grad(lambda) to enforce div(vortg) = 0:
+            call diffx(fs, ds)
+            !$omp parallel
+            !$omp workshare
+            as = as - ds
+            !$omp end workshare
+            !$omp end parallel
+
+            call diffy(fs, ds)
+            !$omp parallel
+            !$omp workshare
+            bs = bs - ds
+            !$omp end workshare
+            !$omp end parallel
+
+            call diffz(fs, ds)
+            !$omp parallel
+            !$omp workshare
+            cs = cs - ds
+            !$omp end workshare
+            !$omp end parallel
+            !Ensure horizontal average of vertical vorticity is zero:
+            cs(:, 0, 0) = zero
+
+            !Compute spectrally filtered vorticity in physical space:
+            !$omp parallel shared(ds, es, fs, as, bs, cs, filt, nz) private(iz) default(none)
+            !$omp do
+            do iz = 0, nz
+                ds(iz, :, :) = filt * as(iz, :, :)
+                es(iz, :, :) = filt * bs(iz, :, :)
+                fs(iz, :, :) = filt * cs(iz, :, :)
+            enddo
+            !$omp end do
+            !$omp end parallel
 
             !Define horizontally-averaged flow by integrating horizontal vorticity:
             ubar(0) = zero
             vbar(0) = zero
             do iz = 0, nz-1
-                ubar(iz+1) = ubar(iz) + dz2 * (bs(iz, 0, 0) + bs(iz+1, 0, 0))
-                vbar(iz+1) = vbar(iz) - dz2 * (as(iz, 0, 0) + as(iz+1, 0, 0))
+                ubar(iz+1) = ubar(iz) + dz2 * (es(iz, 0, 0) + es(iz+1, 0, 0))
+                vbar(iz+1) = vbar(iz) - dz2 * (ds(iz, 0, 0) + ds(iz+1, 0, 0))
             enddo
 
             ! remove the mean value to have zero net momentum
@@ -58,71 +123,86 @@ subroutine vor2vel(vortg,  velog,  velgradg)
                 vbar(iz) = vbar(iz) - vavg
             enddo
 
-            !Form source term for inversion of vertical velocity:
-            call diffy(as, ds)
-            call diffx(bs, es)
+            !Return corrected vorticity to physical space:
+            call fftxys2p(ds, vortg(0:nz, :, :, 1))
+            call fftxys2p(es, vortg(0:nz, :, :, 2))
+            call fftxys2p(fs, vortg(0:nz, :, :, 3))
+
+            !-----------------------------------------------------------------
+            !Invert vorticity to find vector potential (A, B, C) -> (as, bs, cs):
+            call lapinv0(as)
+            call lapinv0(bs)
+            call lapinv1(cs)
+
+            !------------------------------------------------------------
+            !Compute x velocity component, u = B_z - C_y:
+            call diffy(cs, ds)
+            call diffz(bs, es)
 
             !$omp parallel
             !$omp workshare
-            ds = ds - es
+            fs = es - ds
             !$omp end workshare
             !$omp end parallel
-
-            !as & bs are now free to re-use
-
-            !Save boundary values for z derivative of w below:
-            dsbot = ds(0,  :, :)
-            dstop = ds(nz, :, :)
-
-            !Invert to find vertical velocity \hat{w} (store in ds, spectrally):
-            call lapinv0(ds)
-
-            !Find \hat{w}' (store in es, spectrally):
-            call diffz0(ds, es, dsbot, dstop)
-
-            !Find x velocity component \hat{u}:
-            call diffx(es, as)
-            call diffy(cs, bs)
-
-            !$omp parallel do
+            !Add horizontally-averaged flow:
+            fs(:, 0, 0) = ubar
+            !$omp parallel shared(svelog, fs, filt, nz) private(iz) default(none)
+            !$omp do
             do iz = 0, nz
-                as(iz, :, :) = k2l2i * (as(iz, :, :) + bs(iz, :, :))
+                svelog(iz, :, :, 1) = filt * fs(iz, :, :)
             enddo
-            !$omp end parallel do
-
-            !Add horizontally-averaged flow:
-            as(:, 0, 0) = ubar
+            !$omp end do
+            !$omp end parallel
 
-            svelog(:, :, :, 1) = as
+            !------------------------------------------------------------
+            !Compute y velocity component, v = C_x - A_z:
+            call diffx(cs, ds)
+            call diffz(as, es)
 
-            !Get u in physical space:
-            call fftxys2p(as, velog(0:nz, :, :, 1))
+            !$omp parallel
+            !$omp workshare
+            fs = ds - es
+            !$omp end workshare
+            !$omp end parallel
+            !Add horizontally-averaged flow:
+            fs(:, 0, 0) = vbar
+            !$omp parallel shared(svelog, fs, filt, nz) private(iz) default(none)
+            !$omp do
+            do iz = 0, nz
+                svelog(iz, :, :, 2) = filt * fs(iz, :, :)
+            enddo
+            !$omp end do
+            !$omp end parallel
 
-            !Find y velocity component \hat{v}:
-            call diffy(es, as)
-            call diffx(cs, bs)
+            !------------------------------------------------------------
+            !Compute z velocity component, w = A_y - B_x:
+            call diffx(bs, ds)
+            call diffy(as, es)
+            !$omp parallel
+            !$omp workshare
+            fs = es - ds
+            !$omp end workshare
+            !$omp end parallel
 
-            !$omp parallel do
+            !$omp parallel shared(svelog, fs, filt, nz) private(iz) default(none)
+            !$omp do
             do iz = 0, nz
-                as(iz, :, :) = k2l2i * (as(iz, :, :) - bs(iz, :, :))
+                svelog(iz, :, :, 3) = filt * fs(iz, :, :)
             enddo
-            !$omp end parallel do
+            !$omp end do
+            !$omp end parallel
 
-            !Add horizontally-averaged flow:
-            as(:, 0, 0) = vbar
+            ! compute the velocity gradient tensor
+            call vel2vgrad(svelog, velgradg)
 
-            svelog(:, :, :, 2) = as
+            !Get u in physical space:
+            call fftxys2p(svelog(0:nz, :, :, 1), velog(0:nz, :, :, 1))
 
             !Get v in physical space:
-            call fftxys2p(as, velog(0:nz, :, :, 2))
-
-            svelog(:, :, :, 3) = ds
+            call fftxys2p(svelog(0:nz, :, :, 2), velog(0:nz, :, :, 2))
 
             !Get w in physical space:
-            call fftxys2p(ds, velog(0:nz, :, :, 3))
-
-            ! compute the velocity gradient tensor
-            call vel2vgrad(svelog, velgradg)
+            call fftxys2p(svelog(0:nz, :, :, 3), velog(0:nz, :, :, 3))
 
             ! use extrapolation in u and v and anti-symmetry in w to fill z grid points outside domain:
             velog(-1, :, :, 1) =  two * velog(0, :, :, 1) - velog(1, :, :, 1) ! u
@@ -184,16 +264,18 @@ end subroutine vel2vgrad
 
         !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 
-        ! Compute the gridded buoyancy derivatives db/dx and db/dy
-        subroutine buoyancy_derivatives(tbuoyg, dbdx, dbdy)
+        ! Compute the gridded vorticity tendency:
+        subroutine vorticity_tendency(vortg, tbuoyg, velgradg, vtend)
+            double precision, intent(in)  :: vortg(-1:nz+1, 0:ny-1, 0:nx-1, 3)
             double precision, intent(in)  :: tbuoyg(-1:nz+1, 0:ny-1, 0:nx-1)
-            double precision, intent(out) :: dbdx(-1:nz+1, 0:ny-1, 0:nx-1)
-            double precision, intent(out) :: dbdy(-1:nz+1, 0:ny-1, 0:nx-1)
-            double precision              :: b(0:nz, 0:ny-1, 0:nx-1)    ! buoyancy in physical space
-            double precision              :: bs(0:nz, 0:nx-1, 0:ny-1)   ! buoyancy in spectral space
-            double precision              :: ds(0:nz, 0:nx-1, 0:ny-1)   ! buoyancy derivative in spectral space
+            double precision              :: b(0:nz, 0:ny-1, 0:nx-1)
+            double precision, intent(out) :: velgradg(-1:nz+1, 0:ny-1, 0:nx-1, 5)
+            double precision, intent(out) :: vtend(-1:nz+1, 0:ny-1, 0:nx-1, 3)
+            double precision              :: bs(0:nz, 0:nx-1, 0:ny-1) ! spectral buoyancy
+            double precision              :: ds(0:nz, 0:nx-1, 0:ny-1) ! spectral derivatives
+            double precision              :: db(0:nz, 0:ny-1, 0:nx-1) ! buoyancy derivatives
 
-            call start_timer(db_timer)
+            call start_timer(vtend_timer)
 
             ! copy buoyancy
             b = tbuoyg(0:nz, :, :)
@@ -202,21 +284,45 @@ subroutine buoyancy_derivatives(tbuoyg, dbdx, dbdy)
             call fftxyp2s(b, bs)
 
             call diffy(bs, ds)                      ! b_y = db/dy in spectral space
-            call fftxys2p(ds, dbdy(0:nz, :, :))     ! db = b_y in physical space
+            call fftxys2p(ds, db)                   ! db = b_y in physical space
+
+            !$omp parallel
+            !$omp workshare
+            vtend(0:nz, :, :, 1) =  vortg(0:nz, :, :, 1)           * velgradg(0:nz, :, :, 1) & ! \omegax * du/dx
+                                 + (vortg(0:nz, :, :, 2) + ft_cor) * velgradg(0:nz, :, :, 2) & ! \omegay * du/dy
+                                 + (vortg(0:nz, :, :, 3) +  f_cor) *                         &
+                                            (vortg(0:nz, :, :, 2) + velgradg(0:nz, :, :, 4)) & ! \omegaz * du/dz
+                                 + db                                                          ! db/dy
+            !$omp end workshare
+            !$omp end parallel
 
             call diffx(bs, ds)                      ! b_x = db/dx in spectral space
-            call fftxys2p(ds, dbdx(0:nz, :, :))     ! db = b_x in physical space
+            call fftxys2p(ds, db)                   ! db = b_x in physical space
 
-            ! Extrapolate to halo grid points
-            dbdy(-1,   :, :) = two * dbdy(0,  :, :) - dbdy(1,    :, :)
-            dbdy(nz+1, :, :) = two * dbdy(nz, :, :) - dbdy(nz-1, :, :)
+            !$omp parallel
+            !$omp workshare
+            vtend(0:nz, :, :, 2) =  vortg(0:nz, :, :, 1)           *                         &
+                                            (vortg(0:nz, :, :, 3) + velgradg(0:nz, :, :, 2)) & ! \omegax * dv/dx
+                                 + (vortg(0:nz, :, :, 2) + ft_cor) * velgradg(0:nz, :, :, 3) & ! \omegay * dv/dy
+                                 + (vortg(0:nz, :, :, 3) + f_cor)  *                         &
+                                            (velgradg(0:nz, :, :, 5) - vortg(0:nz, :, :, 1)) & ! \omegaz * dv/dz
+                                 - db                                                          ! dbdx
+
+            vtend(0:nz, :, :, 3) =  vortg(0:nz, :, :, 1)           * velgradg(0:nz, :, :, 4) & ! \omegax * dw/dx
+                                 + (vortg(0:nz, :, :, 2) + ft_cor) * velgradg(0:nz, :, :, 5) & ! \omegay * dw/dy
+                                 - (vortg(0:nz, :, :, 3) + f_cor)  *                         &
+                                         (velgradg(0:nz, :, :, 1) + velgradg(0:nz, :, :, 3))   ! \omegaz * dw/dz
 
-            dbdx(-1,   :, :) = two * dbdx(0,  :, :) - dbdx(1,    :, :)
-            dbdx(nz+1, :, :) = two * dbdx(nz, :, :) - dbdx(nz-1, :, :)
+            !$omp end workshare
+            !$omp end parallel
 
-            call stop_timer(db_timer)
+            ! Extrapolate to halo grid points
+            vtend(-1,   :, :, :) = two * vtend(0,  :, :, :) - vtend(1,    :, :, :)
+            vtend(nz+1, :, :, :) = two * vtend(nz, :, :, :) - vtend(nz-1, :, :, :)
 
-        end subroutine buoyancy_derivatives
+            call stop_timer(vtend_timer)
+
+        end subroutine vorticity_tendency
 
         !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 
@@ -252,7 +358,11 @@ subroutine diverge(div,  ud, vd, wd)
             call fftxys2p(vs, vd)
 
             ! Compute z derivative by compact differences:
-            call diffz1(ds, ws)
+            call diffz(ds, ws)
+
+            ! Set vertical boundary values to zero
+            ws(0,  :, :) = zero
+            ws(nz, :, :) = zero
 
             ! Add on the x and y-independent part of wd:
             ws(:, 1, 1) = ws(:, 1, 1) + wbar
diff --git a/src/3d/inversion/inversion_utils.f90 b/src/3d/inversion/inversion_utils.f90
index 89c17c7dc..7c0736378 100644
--- a/src/3d/inversion/inversion_utils.f90
+++ b/src/3d/inversion/inversion_utils.f90
@@ -11,31 +11,27 @@ module inversion_utils
     ! Ordering in spectral space: z, x, y
 
     ! Tridiagonal arrays for the horizontal vorticity components:
-    double precision, allocatable :: etdh(:, :, :), htdh(:, :, :), ap(:, :), apb(:, :)
+    double precision, allocatable :: etdh(:, :, :), htdh(:, :, :)
 
     ! Tridiagonal arrays for the vertical vorticity component:
     double precision, allocatable :: etdv(:, :, :), htdv(:, :, :)
 
-    ! Tridiagonal arrays for the compact difference calculation of d/dz:
-    double precision, allocatable :: etd0(:), htd0(:)
-    double precision, allocatable :: etd1(:), htd1(:)
-
-    ! Tridiagonal arrays for integrating in z:
-    double precision, allocatable :: etda(:), htda(:)
-
     !Horizontal wavenumbers:
     double precision, allocatable :: rkx(:), hrkx(:), rky(:), hrky(:)
 
-    ! Note k2l2i = 1/(k^2+l^2) (except k = l = 0, then k2l2i(0, 0) = 0)
-    double precision, allocatable :: k2l2i(:, :)
-
     !Quantities needed in FFTs:
     double precision, allocatable :: xtrig(:), ytrig(:)
     integer :: xfactors(5), yfactors(5)
     integer, parameter :: nsubs_tri = 8 !number of blocks for openmp
     integer :: nxsub
 
-    double precision :: dz, dzi, dz2, dz6, dz24, hdzi, dzisq
+    !De-aliasing filter:
+    double precision, allocatable :: filt(:, :)
+    double precision, allocatable :: skx(:), sky(:)
+
+
+
+    double precision :: dz, dzi, dz2, dz6, dz24, hdzi, dzisq, ap
     integer :: nwx, nwy, nxp2, nyp2
 
     logical :: is_initialised = .false.
@@ -43,20 +39,19 @@ module inversion_utils
     public :: init_fft  &
             , diffx     &
             , diffy     &
-            , diffz0    &
-            , diffz1    &
+            , diffz     &
             , lapinv0   &
             , lapinv1   &
             , vertint   &
             , fftxyp2s  &
             , fftxys2p  &
             , dz2       &
+            , filt      &
             , hdzi      &
             , xfactors  &
             , yfactors  &
             , xtrig     &
-            , ytrig     &
-            , k2l2i
+            , ytrig
 
     contains
 
@@ -65,9 +60,10 @@ module inversion_utils
         !Initialises this module (FFTs, x & y wavenumbers, tri-diagonal
         !coefficients, etc).
         subroutine init_fft
-            double precision, allocatable  :: a0(:, :), a0b(:, :), ksq(:, :)
+            double precision, allocatable  :: a0(:, :), ksq(:, :)
             double precision               :: rkxmax, rkymax
             double precision               :: rksqmax
+            double precision               :: kxmaxi, kymaxi
             integer                        :: kx, ky, iz, isub, ib_sub, ie_sub
 
             if (is_initialised) then
@@ -89,28 +85,21 @@ subroutine init_fft
             nxp2 = nx + 2
 
             allocate(a0(nx, ny))
-            allocate(a0b(nx, ny))
             allocate(ksq(nx, ny))
-            allocate(k2l2i(nx, ny))
+            allocate(skx(nx))
+            allocate(sky(ny))
 
             allocate(etdh(nz-1, nx, ny))
             allocate(htdh(nz-1, nx, ny))
-            allocate(ap(nx, ny))
-            allocate(apb(nx, ny))
             allocate(etdv(0:nz, nx, ny))
             allocate(htdv(0:nz, nx, ny))
-            allocate(etd0(0:nz))
-            allocate(htd0(0:nz))
-            allocate(etd1(nz-1))
-            allocate(htd1(nz-1))
-            allocate(etda(nz))
-            allocate(htda(nz))
             allocate(rkx(nx))
             allocate(hrkx(nx))
             allocate(rky(ny))
             allocate(hrky(ny))
             allocate(xtrig(2 * nx))
             allocate(ytrig(2 * ny))
+            allocate(filt(nx, ny))
 
             nxsub = nx / nsubs_tri
 
@@ -146,16 +135,22 @@ subroutine init_fft
                 enddo
             enddo
 
-            ksq(1, 1) = one
-            k2l2i = one / ksq
-            ksq(1, 1) = zero
+            !--------------------------------------------------------------------
+            ! Define Hou and Li filter:
+            kxmaxi = one / maxval(rkx)
+            skx = -36.d0 * (kxmaxi * rkx) ** 36
+            kymaxi = one/maxval(rky)
+            sky = -36.d0 * (kymaxi * rky) ** 36
+            do ky = 1, ny
+                do kx = 1, nx
+                    filt(kx, ky) = dexp(skx(kx) + sky(ky))
+                enddo
+            enddo
 
             !-----------------------------------------------------------------------
             ! Fixed coefficients used in the tridiagonal problems:
-            a0 = -two * dzisq - f56 * ksq
-            a0b = -dzisq - f13 * ksq
-            ap = dzisq - f112 * ksq
-            apb = dzisq - f16 * ksq
+            a0 = -two * dzisq - ksq
+            ap = dzisq
 
             !-----------------------------------------------------------------------
             ! Tridiagonal arrays for the horizontal vorticity components:
@@ -168,8 +163,8 @@ subroutine init_fft
                 ie_sub = (isub + 1) * nxsub
                 do iz = 2, nz-2
                     htdh(iz, ib_sub:ie_sub, :) = one / (a0(ib_sub:ie_sub, :) &
-                                               + ap(ib_sub:ie_sub, :) * etdh(iz-1, ib_sub:ie_sub, :))
-                    etdh(iz, ib_sub:ie_sub, :) = -ap(ib_sub:ie_sub, :) * htdh(iz, ib_sub:ie_sub, :)
+                                               + ap * etdh(iz-1, ib_sub:ie_sub, :))
+                    etdh(iz, ib_sub:ie_sub, :) = -ap * htdh(iz, ib_sub:ie_sub, :)
                 enddo
             enddo
             !$omp end do
@@ -180,8 +175,8 @@ subroutine init_fft
             etdh(:, 1, 1) = zero
 
             ! Tridiagonal arrays for the vertical vorticity component:
-            htdv(0, :, :) = one / a0b
-            etdv(0, :, :) = -apb * htdv(0, :, :)
+            htdv(0, :, :) = one / a0
+            etdv(0, :, :) = -two * ap * htdv(0, :, :)
             !$omp parallel shared(a0, ap, etdv, htdv, nz, nxsub) private(isub, ib_sub, ie_sub, iz) default(none)
             !$omp do
             do isub = 0, nsubs_tri-1
@@ -189,8 +184,8 @@ subroutine init_fft
                 ie_sub = (isub + 1) * nxsub
                 do iz = 1, nz-1
                     htdv(iz, ib_sub:ie_sub, :) = one / (a0(ib_sub:ie_sub, :) &
-                                               + ap(ib_sub:ie_sub, :) * etdv(iz-1, ib_sub:ie_sub, :))
-                    etdv(iz, ib_sub:ie_sub, :) = -ap(ib_sub:ie_sub, :) * htdv(iz, ib_sub:ie_sub, :)
+                                               + ap * etdv(iz-1, ib_sub:ie_sub, :))
+                    etdv(iz, ib_sub:ie_sub, :) = -ap * htdv(iz, ib_sub:ie_sub, :)
                 enddo
             enddo
             !$omp end do
@@ -198,43 +193,12 @@ subroutine init_fft
 
             etdv(nz-1, 1, 1) = zero
 
-            htdv(nz, :, :) = one / (a0b + apb * etdv(nz-1, :, :))
+            htdv(nz, :, :) = one / (a0 + two * ap * etdv(nz-1, :, :))
             ! Remove horizontally-averaged part (done separately):
             htdv(:, 1, 1) = zero
             etdv(:, 1, 1) = zero
 
-            ! Tridiagonal arrays for the compact difference calculation of d/dz
-            ! for fields f for which f  =  0 at the boundaries:
-            htd0(0) = one / f23
-            etd0(0) = -f13 * htd0(0)
-            do iz = 1, nz-1
-                htd0(iz) = one / (f23 + f16 * etd0(iz-1))
-                etd0(iz) = -f16 * htd0(iz)
-            enddo
-            htd0(nz) = one / (f23 + f13 * etd0(nz-1))
-
-            ! Tridiagonal arrays for the compact difference calculation of d /dz
-            ! for fields f for which df/dz = 0 at the boundaries:
-            htd1(1) = one / f23
-            etd1(1) = -f16 * htd1(1)
-            do iz = 2, nz-2
-                htd1(iz) = one / (f23 + f16 * etd1(iz-1))
-                etd1(iz) = -f16 * htd1(iz)
-            enddo
-            htd1(nz-1) = one / (f23 + f16 * etd1(nz-2))
-
-            ! Tridiagonal arrays used for integrating in z (see vertint):
-            htda(1) = one / (one + f16)
-            etda(1) = -f16 * htda(1)
-            do iz = 2, nz-1
-                htda(iz) = one / (one + f16 * etda(iz-1))
-                etda(iz) = -f16 * htda(iz)
-            enddo
-            htda(nz) = one / (one + f16 + f16 * etda(nz-2))
-
-
             deallocate(a0)
-            deallocate(a0b)
             deallocate(ksq)
         end subroutine
 
@@ -293,141 +257,51 @@ subroutine diffy(fs,ds)
 
         !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 
-        !Calculates df/dz for a field f which has f = 0 at the boundaries
-        !using 4th-order compact differencing.  Here fs = f, ds = df/dz,
-        !lapfsbot = Laplace(f) at z_min and lapfstop = Laplace(f) at z_max,
-        !both given.  *** All quantities must be in semi-spectral space ***
-        subroutine diffz0(fs, ds, lapfsbot, lapfstop)
-            double precision, intent(in)  :: fs(0:nz, nx, ny), &
-                                             lapfsbot(nx, ny), &
-                                             lapfstop(nx, ny)
+        !Calculates df/dz for a field f using 2nd-order differencing.
+        !Here fs = f, ds = df/dz.
+        subroutine diffz(fs, ds)
+            double precision, intent(in)  :: fs(0:nz, nx, ny)
             double precision, intent(out) :: ds(0:nz, nx, ny)
-            integer                       :: iz, isub, ib_sub, ie_sub
+            integer                       :: iz
 
-            ds(0, :, :) = fs(1, :, :) * dzi - dz6 * lapfsbot
-            ds(1, :, :) = fs(2, :, :) * hdzi
+            ! forward and backward differencing for boundary cells
+            ! iz = 0:  (fs(1) - fs(0)) / dz
+            ! iz = nz: (fs(nz) - fs(nz-1)) / dz
+            ds(0,  :, :) = dzi * (fs(1,    :, :) - fs(0,    :, :))
+            ds(nz, :, :) = dzi * (fs(nz,   :, :) - fs(nz-1, :, :))
 
+            ! central differencing for interior cells
             !$omp parallel shared(ds, fs, hdzi, nz) private(iz) default(none)
             !$omp do
-            do iz = 2, nz-2
-                ds(iz, :, :) = (fs(iz+1, :, :) - fs(iz-1, :, :)) * hdzi
-            enddo
-            !$omp end do
-            !$omp end parallel
-
-            ds(nz-1, :, :) =  - fs(nz-2, :, :) * hdzi
-            ds(nz, :, :) = f16 * dz * lapfstop - fs(nz-1, :, :) * dzi
-
-            ds(0, :, :) = ds(0, :, :) * htd0(0)
-
-            !$omp parallel shared(ds, htd0, nz, nxsub) private(isub, ib_sub, ie_sub,iz) default(none)
-            !$omp do
-            do isub = 0, nsubs_tri-1
-                ib_sub = isub * nxsub+1
-                ie_sub = (isub + 1) * nxsub
-                do iz = 1, nz-1
-                    ds(iz, ib_sub:ie_sub, :) = (ds(iz, ib_sub:ie_sub, :) &
-                                             - f16 * ds(iz-1, ib_sub:ie_sub, :)) * htd0(iz)
-                enddo
-            enddo
-            !$omp end do
-            !$omp end parallel
-
-            ds(nz, :, :) = (ds(nz, :, :) - f13 * ds(nz-1, :, :)) * htd0(nz)
-
-            !$omp parallel shared(ds, etd0, nz, nxsub) private(isub, ib_sub, ie_sub, iz) default(none)
-            !$omp do
-            do isub = 0, nsubs_tri-1
-                ib_sub = isub * nxsub + 1
-                ie_sub = (isub + 1) * nxsub
-                do iz = nz-1, 0, -1
-                    ds(iz, ib_sub:ie_sub, :) = etd0(iz) * ds(iz+1, ib_sub:ie_sub, :) &
-                                             + ds(iz, ib_sub:ie_sub, :)
-                enddo
-            enddo
-            !$omp end do
-            !$omp end parallel
-        end subroutine
-
-        !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
-
-        !Calculates df/dz for a field f which has df/dz = 0 at the boundaries
-        !using 4th-order compact differencing.  Here fs = f and ds = df/dz.
-        subroutine diffz1(fs, ds)
-            double precision, intent(in)  :: fs(0:nz, nx, ny)
-            double precision, intent(out) :: ds(0:nz, nx, ny)
-            integer                       :: iz, isub, ib_sub, ie_sub
-
-            !$omp parallel shared(ds, fs, nz, hdzi) private(iz) default(none)
-            !$omp do
             do iz = 1, nz-1
                 ds(iz, :, :) = (fs(iz+1, :, :) - fs(iz-1, :, :)) * hdzi
             enddo
             !$omp end do
             !$omp end parallel
 
-            ds(0, :, :) = zero
-            ds(1, :, :) = ds(1, :, :) * htd1(1)
-
-            !$omp parallel shared(ds, htd1, nz, nxsub) private(isub, ib_sub, ie_sub, iz) default(none)
-            !$omp do
-            do isub = 0, nsubs_tri-1
-                ib_sub = isub * nxsub + 1
-                ie_sub = (isub + 1) * nxsub
-                do iz = 2, nz-1
-                    ds(iz, ib_sub:ie_sub, :) = (ds(iz, ib_sub:ie_sub, :) &
-                                             - f16 * ds(iz-1, ib_sub:ie_sub, :)) * htd1(iz)
-                enddo
-            enddo
-            !$omp end do
-            !$omp end parallel
-
-            ds(nz, :, :) = zero
-
-            !$omp parallel shared(ds, etd1, nz, nxsub) private(isub, ib_sub, ie_sub, iz) default(none)
-            !$omp do
-            do isub = 0, nsubs_tri-1
-                ib_sub = isub * nxsub + 1
-                ie_sub = (isub + 1) * nxsub
-                do iz = nz-2, 1, -1
-                    ds(iz, ib_sub:ie_sub, :) = etd1(iz) * ds(iz+1, ib_sub:ie_sub, :) &
-                                             + ds(iz, ib_sub:ie_sub, :)
-                enddo
-            end do
-            !$omp end do
-            !$omp end parallel
         end subroutine
 
         !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 
         !Inverts Laplace's operator on fs in semi-spectral space.
         !Here fs = 0 on the z boundaries.
-        !Uses 4th-order compact differencing
+        !Uses 2nd-order differencing
         !*** Overwrites fs ***
         subroutine lapinv0(fs)
             double precision, intent(inout) :: fs(0:nz, nx, ny)
-            double precision                :: rs(nz-1, nx, ny)
             integer                         :: iz, isub, ib_sub, ie_sub
 
-            !$omp parallel shared(rs, fs, nz) private(iz) default(none)
-            !$omp do
-            do iz = 1, nz-1
-                rs(iz, :, :) = f112 * (fs(iz-1, :, :) + fs(iz+1, :, :)) + f56 * fs(iz, :, :)
-            enddo
-            !$omp end do
-            !$omp end parallel
-
             fs(0, :, :) = zero
-            fs(1, :, :) = rs(1, :, :) * htdh(1, :, :)
+            fs(1, :, :) = fs(1, :, :) * htdh(1, :, :)
 
-            !$omp parallel shared(rs, fs, ap, htdh, nz, nxsub) private(isub, ib_sub, ie_sub, iz) default(none)
+            !$omp parallel shared(fs, ap, htdh, nz, nxsub) private(isub, ib_sub, ie_sub, iz) default(none)
             !$omp do
             do isub = 0, nsubs_tri-1
                 ib_sub = isub * nxsub + 1
                 ie_sub = (isub + 1) * nxsub
                 do iz = 2, nz-1
-                    fs(iz, ib_sub:ie_sub, :) = (rs(iz, ib_sub:ie_sub, :) &
-                                             - ap(ib_sub:ie_sub, :) * fs(iz-1, ib_sub:ie_sub, :)) &
+                    fs(iz, ib_sub:ie_sub, :) = (fs(iz, ib_sub:ie_sub, :) &
+                                             - ap * fs(iz-1, ib_sub:ie_sub, :)) &
                                              * htdh(iz, ib_sub:ie_sub, :)
                 enddo
             enddo
@@ -454,23 +328,14 @@ subroutine lapinv0(fs)
 
         !Inverts Laplace's operator on fs in semi-spectral space.
         !Here dfs/dz = 0 on the z boundaries.
-        !Uses 4th-order compact differencing
+        !Uses 2nd-order differencing
         !*** Overwrites fs ***
         subroutine lapinv1(fs)
             double precision, intent(inout) :: fs(0:nz, nx, ny)
             double precision                :: rs(0:nz, nx, ny)
             integer                         :: iz, isub, ib_sub, ie_sub
 
-            rs(0, :, :) = f13 * fs(0, :, :) + f16 * fs(1, :, :)
-            !$omp parallel shared(rs, fs) private(iz)
-            !$omp do
-            do iz = 1, nz-1
-                rs(iz, :, :) = f112 * (fs(iz-1, :, :) + fs(iz+1, :, :)) + f56 * fs(iz, :, :)
-            enddo
-            !$omp end do
-            !$omp end parallel
-            rs(nz, :, :) = f13 * fs(nz, :, :) + f16 * fs(nz-1, :, :)
-
+            rs = fs
             fs(0, :, :) = rs(0, :, :) * htdv(0, :, :)
 
             !$omp parallel shared(rs, fs, ap, htdv, nz, nxsub) private(isub, ib_sub, ie_sub, iz) default(none)
@@ -480,14 +345,14 @@ subroutine lapinv1(fs)
                 ie_sub = (isub + 1) * nxsub
                 do iz = 1, nz-1
                     fs(iz, ib_sub:ie_sub, :) = (rs(iz, ib_sub:ie_sub, :) &
-                                             - ap(ib_sub:ie_sub, :) * fs(iz-1, ib_sub:ie_sub, :)) &
+                                             - ap * fs(iz-1, ib_sub:ie_sub, :)) &
                                              * htdv(iz, ib_sub:ie_sub, :)
                 enddo
             enddo
             !$omp end do
             !$omp end parallel
 
-            fs(nz, :, :) = (rs(nz, :, :) - apb * fs(nz-1, :, :)) * htdv(nz, :, :)
+            fs(nz, :, :) = (rs(nz, :, :) - two * ap * fs(nz-1, :, :)) * htdv(nz, :, :)
 
             !$omp parallel shared(fs, etdv, nz, nxsub) private(isub, ib_sub, ie_sub, iz) default(none)
             !$omp do
@@ -509,42 +374,37 @@ subroutine lapinv1(fs)
         !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 
         !Finds f by integrating df/dz = d, ensuring f = 0 at the boundaries
-        !using 4th-order compact differencing.  Here ds = df/dz and fs = f.
+        !using trapezoidal rule.  Here ds = df/dz and fs = f.
         subroutine vertint(ds, fs)
             double precision, intent(in)  :: ds(0:nz)
             double precision, intent(out) :: fs(0:nz)
-            double precision              :: es(nz), esum
+            double precision              :: c
             integer                       :: iz
 
-            !-------------------------------------------
-            !First interpolate ds to a half grid as es:
+            ! set lower boundary value
+            fs(0)  = zero
+
             do iz = 1, nz
-                es(iz) = f23 * (ds(iz-1) + ds(iz))
+                fs(iz) = fs(iz-1) + dz2 * (ds(iz) + ds(iz-1))
             enddo
 
-            es(1) = es(1) * htda(1)
-            do iz = 2, nz
-                es(iz) = (es(iz) - f16 * es(iz-1)) * htda(iz)
-            enddo
+            ! shift to adjust f(nz) to be zero
+            c = fs(nz) / dble(nz)
 
-            do iz = nz-1, 1, -1
-                es(iz) = etda(iz) * es(iz+1) + es(iz)
+            !$omp parallel private(iz)
+            !$omp do
+            do iz = 1, nz-1
+                fs(iz) = fs(iz) - c * dble(iz)
             enddo
+            !$omp end do
+            !$omp end parallel
 
-            !-------------------------------------------
-            !Next adjust es to ensure f(nz) = 0:
-            esum = (f1112 * (es(1) + es(nz)) + sum(es(2:nz-1))) / (dble(nz) - f16)
-            es = es - esum
+            ! set upper boundary value
+            fs(nz)  = zero
 
-            !Integrate:
-            fs(0) = zero
-            fs(1) = dz24 * (23.0d0 * es(1) + es(2))
-            do iz = 2, nz-1
-                fs(iz) = fs(iz-1) + dz24 * (es(iz-1) + 22.d0 * es(iz) + es(iz+1))
-            enddo
-            fs(nz) = zero
         end subroutine
 
+
         !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 
         ! Computes a 2D FFT (in x & y) of a 3D array fp in physical space
diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90
index 700707331..d597f1e95 100644
--- a/src/3d/parcels/parcel_interpl.f90
+++ b/src/3d/parcels/parcel_interpl.f90
@@ -292,7 +292,6 @@ subroutine grid2par(vel, vortend, vgrad, add)
             double precision,     intent(inout) :: vel(:, :), vortend(:, :), vgrad(:, :)
             logical, optional, intent(in)       :: add
             integer                             :: n, l
-            double precision                    :: dudz, dvdz, dwdz, dvdx
 
             call start_timer(grid2par_timer)
 
@@ -320,14 +319,10 @@ subroutine grid2par(vel, vortend, vgrad, add)
             endif
 
             !$omp parallel default(shared)
-            !$omp do private(n, l, is, js, ks, weights, dudz, dvdz, dwdz, dvdx)
+            !$omp do private(n, l, is, js, ks, weights)
             do n = 1, n_parcels
 
                 vgrad(:, n) = zero
-                dudz = zero
-                dvdz = zero
-                dwdz = zero
-                dvdx = zero
 
                 ! ensure point is within the domain
                 call apply_periodic_bc(parcels%position(:, n))
@@ -341,47 +336,8 @@ subroutine grid2par(vel, vortend, vgrad, add)
 
                     vgrad(:, n) = vgrad(:, n) + weights(l) * velgradg(ks(l), js(l), is(l), :)
 
-                    ! du/dz = \omegay + dw/dx
-                    dudz = dudz &
-                         + weights(l) * (vortg(ks(l), js(l), is(l), 2) + velgradg(ks(l), js(l), is(l), 4))
-
-                    ! dv/dz = dw/dy - \omegax
-                    dvdz = dvdz &
-                         + weights(l) * (velgradg(ks(l), js(l), is(l), 5) - vortg(ks(l), js(l), is(l), 1))
-
-                    ! dw/dz = - (du/dx + dv/dy)
-                    dwdz = dwdz &
-                         - weights(l) * (velgradg(ks(l), js(l), is(l), 1) + velgradg(ks(l), js(l), is(l), 3))
-
-                    ! dv/dx = \omegaz + du/dy
-                    dvdx = dvdx &
-                         + weights(l) * (vortg(ks(l), js(l), is(l), 3) + velgradg(ks(l), js(l), is(l), 2))
-
-                    ! add buoyancy part of vorticity tendency to x-vorticity
-                    vortend(1, n) = vortend(1, n) + weights(l) * dbdy(ks(l), js(l), is(l))
-
-                    ! subtract buoyancy part of vorticity tendency to y-vorticity
-                    vortend(2, n) = vortend(2, n) - weights(l) * dbdx(ks(l), js(l), is(l))
+                    vortend(:, n) = vortend(:, n) + weights(l) * vtend(ks(l), js(l), is(l), :)
                 enddo
-
-                ! add strain part of vorticity tendency to x-vorticity
-                vortend(1, n) =  vortend(1, n)                                   &
-                              +  parcels%vorticity(1, n)           * vgrad(1, n) & ! \omegax * du/dx
-                              + (parcels%vorticity(2, n) + ft_cor) * vgrad(2, n) & ! \omegay * du/dy
-                              + (parcels%vorticity(3, n) +  f_cor) * dudz          ! \omegaz * du/dz
-
-                ! add strain part of vorticity tendency to y-vorticity
-                vortend(2, n) =  vortend(2, n)                                   &
-                              +  parcels%vorticity(1, n)           * dvdx        & ! \omegax * dv/dx
-                              + (parcels%vorticity(2, n) + ft_cor) * vgrad(3, n) & ! \omegay * dv/dy
-                              + (parcels%vorticity(3, n) + f_cor)  * dvdz          ! \omegaz * dv/dz
-
-                ! add strain part of vorticity tendency to z-vorticity
-                vortend(3, n) =  vortend(3, n)                                   &
-                              +  parcels%vorticity(1, n)           * vgrad(4, n) & ! \omegax * dw/dx
-                              + (parcels%vorticity(2, n) + ft_cor) * vgrad(5, n) & ! \omegay * dw/dy
-                              + (parcels%vorticity(3, n) + f_cor)  * dwdz          ! \omegaz * dw/dz
-
             enddo
             !$omp end do
             !$omp end parallel
diff --git a/src/3d/stepper/ls_rk4.f90 b/src/3d/stepper/ls_rk4.f90
index 11b70ff3f..307cb623f 100644
--- a/src/3d/stepper/ls_rk4.f90
+++ b/src/3d/stepper/ls_rk4.f90
@@ -9,8 +9,8 @@ module ls_rk4
     use rk4_utils, only: get_dBdt, get_time_step
     use utils, only : write_step
     use parcel_interpl, only : par2grid, grid2par, grid2par_add
-    use fields, only : velgradg, velog, vortg, dbdx, dbdy, tbuoyg
-    use inversion_mod, only : vor2vel, buoyancy_derivatives
+    use fields, only : velgradg, velog, vortg, vtend, tbuoyg
+    use inversion_mod, only : vor2vel, vorticity_tendency
     use parcel_diagnostics, only : calculate_parcel_diagnostics
     use field_diagnostics, only : calculate_field_diagnostics
     use parameters, only : nx, nz
@@ -80,7 +80,7 @@ subroutine ls_rk4_step(t)
             ! this is also needed for the first ls-rk4 substep
             call vor2vel(vortg, velog, velgradg)
 
-            call buoyancy_derivatives(tbuoyg, dbdx, dbdy)
+            call vorticity_tendency(vortg, tbuoyg, velgradg, vtend)
 
             ! update the time step
             dt = get_time_step(t)
@@ -146,7 +146,7 @@ subroutine ls_rk4_substep(ca, cb, dt, step)
             else
                 call vor2vel(vortg, velog, velgradg)
 
-                call buoyancy_derivatives(tbuoyg, dbdx, dbdy)
+                call vorticity_tendency(vortg, tbuoyg, velgradg, vtend)
 
                 call grid2par_add(delta_pos, delta_vor, strain)
 
diff --git a/src/3d/stepper/rk4_utils.f90 b/src/3d/stepper/rk4_utils.f90
index 7074e44f6..01d865aa7 100644
--- a/src/3d/stepper/rk4_utils.f90
+++ b/src/3d/stepper/rk4_utils.f90
@@ -2,7 +2,7 @@ module rk4_utils
     use parcel_ellipsoid, only : get_B33
     use fields, only : velgradg, tbuoyg, vortg
     use constants, only : zero, one, two, f12
-    use parameters, only : nx, ny, nz, dxi
+    use parameters, only : nx, ny, nz, dxi, vcell
     use jacobi, only : jacobi_eigenvalues
 #ifdef ENABLE_VERBOSE
     use options, only : output
diff --git a/src/3d/utils/utils.f90 b/src/3d/utils/utils.f90
index 0007e6de0..bba92aa39 100644
--- a/src/3d/utils/utils.f90
+++ b/src/3d/utils/utils.f90
@@ -8,7 +8,7 @@ module utils
     use parcel_diagnostics_netcdf
     use parcel_diagnostics
     use parcel_container, only : n_parcels
-    use inversion_mod, only : vor2vel, buoyancy_derivatives
+    use inversion_mod, only : vor2vel, vorticity_tendency
     use parcel_interpl, only : par2grid, grid2par
     use netcdf_reader, only : get_file_type, get_num_steps, get_time, get_netcdf_box
     use parameters, only : lower, extent, update_parameters
@@ -69,7 +69,7 @@ subroutine write_last_step(t)
             ! this is also needed for the first ls-rk4 substep
             call vor2vel(vortg, velog, velgradg)
 
-            call buoyancy_derivatives(tbuoyg, dbdx, dbdy)
+            call vorticity_tendency(vortg, tbuoyg, velgradg, vtend)
 
             call grid2par(velocity, vorticity, strain)
 
diff --git a/unit-tests/3d/Makefile.am b/unit-tests/3d/Makefile.am
index 657e6e724..bf9dfec75 100644
--- a/unit-tests/3d/Makefile.am
+++ b/unit-tests/3d/Makefile.am
@@ -45,6 +45,7 @@ unittests_PROGRAMS = 				\
 	test_velgradg				\
 	test_lapinv0_1				\
 	test_lapinv0_2				\
+	test_lapinv1				\
 	test_diffz0				\
 	test_diffz1
 
@@ -117,6 +118,9 @@ test_lapinv0_1_LDADD = libcombi.la
 test_lapinv0_2_SOURCES = test_lapinv0_2.f90
 test_lapinv0_2_LDADD = libcombi.la
 
+test_lapinv1_SOURCES = test_lapinv1.f90
+test_lapinv1_LDADD = libcombi.la
+
 test_diffz0_SOURCES = test_diffz0.f90
 test_diffz0_LDADD = libcombi.la
 
diff --git a/unit-tests/3d/test_diffz0.f90 b/unit-tests/3d/test_diffz0.f90
index 8b83e5609..fe9299286 100644
--- a/unit-tests/3d/test_diffz0.f90
+++ b/unit-tests/3d/test_diffz0.f90
@@ -1,25 +1,24 @@
 ! =============================================================================
-!                       Test subroutine diffz0
+!                       Test subroutine diffz
 !
-!  This unit test checks the subroutine diffz0 using the
+!  This unit test checks the subroutine diffz using the
 !  function:
 !               cos(k * x) * sin(l * y) * sin(m * z)
 !  where k = 2pi/L_x, l = 2pi/L_y and m = pi/L_z and where x, y and z all start
 !  at 0 (one could start at -pi for x and y just as well).
-!  The subroutine diffz0 should return
+!  The subroutine diffz should return
 !               m * cos(k * x) * sin(l * y) * cos(m * z)
 ! =============================================================================
-program test_diffz0
+program test_diffz
     use unit_test
     use constants, only : zero, one, two, pi, twopi
     use parameters, only : lower, update_parameters, dx, nx, ny, nz, extent
     use inversion_utils, only : init_fft
-    use inversion_mod, only : diffz0
+    use inversion_mod, only : diffz
     implicit none
 
     double precision              :: error
     double precision, allocatable :: fs(:, :, :), ds(:, :, :), &
-                                     fsbot(:, :), fstop(:, :), &
                                      ref_sol(:, :, :)
     integer                       :: ix, iy, iz
     double precision              :: x, y, z, k, l, m, prefactor
@@ -32,8 +31,6 @@ program test_diffz0
 
     allocate(fs(0:nz, nx, ny))
     allocate(ds(0:nz, nx, ny))
-    allocate(fsbot(nx, ny))
-    allocate(fstop(nx, ny))
     allocate(ref_sol(0:nz, nx, ny))
 
     call update_parameters
@@ -58,21 +55,16 @@ program test_diffz0
         enddo
     enddo
 
-    fsbot = fs(0, :, :)
-    fstop = fs(nz, :, :)
-
     call init_fft
 
-    call diffz0(fs, ds, fsbot, fstop)
+    call diffz(fs, ds)
 
     error = maxval(dabs(ds - ref_sol))
 
-    call print_result_dp('Test inversion (diffz0)', error, atol=6.0e-10)
+    call print_result_dp('Test inversion (diffz)', error, atol=2.6e-5)
 
     deallocate(fs)
     deallocate(ds)
-    deallocate(fsbot)
-    deallocate(fstop)
     deallocate(ref_sol)
 
-end program test_diffz0
+end program test_diffz
diff --git a/unit-tests/3d/test_diffz1.f90 b/unit-tests/3d/test_diffz1.f90
index e14440edd..5c83c0d3c 100644
--- a/unit-tests/3d/test_diffz1.f90
+++ b/unit-tests/3d/test_diffz1.f90
@@ -1,20 +1,20 @@
 ! =============================================================================
-!                       Test subroutine diffz1
+!                       Test subroutine diffz
 !
-!  This unit test checks the subroutine diffz1 using the
+!  This unit test checks the subroutine diffz using the
 !  function:
 !               cos(k * x) * sin(l * y) * cos(m * z)
 !  where k = 2pi/L_x, l = 2pi/L_y and m = pi/L_z and where x, y and z all start
 !  at 0 (one could start at -pi for x and y just as well).
-!  The subroutine diffz1 should return
+!  The subroutine diffz should return
 !               - m * cos(k * x) * sin(l * y) * sin(m * z)
 ! =============================================================================
-program test_diffz1
+program test_diffz
     use unit_test
     use constants, only : zero, one, two, pi, twopi
     use parameters, only : lower, update_parameters, dx, nx, ny, nz, extent
     use inversion_utils, only : init_fft
-    use inversion_mod, only : diffz1
+    use inversion_mod, only : diffz
     implicit none
 
     double precision              :: error
@@ -57,14 +57,18 @@ program test_diffz1
 
     call init_fft
 
-    call diffz1(fs, ds)
+    call diffz(fs, ds)
+
+    ! we need to explicitly set to zero
+    ds(0,  :, :) = zero
+    ds(nz, :, :) = zero
 
     error = maxval(dabs(ds - ref_sol))
 
-    call print_result_dp('Test inversion (diffz1)', error, atol=1.0e-7)
+    call print_result_dp('Test inversion (diffz)', error, atol=5.0e-4)
 
     deallocate(fs)
     deallocate(ds)
     deallocate(ref_sol)
 
-end program test_diffz1
+end program test_diffz
diff --git a/unit-tests/3d/test_grid2par.f90 b/unit-tests/3d/test_grid2par.f90
index 3acea835a..33074ee13 100644
--- a/unit-tests/3d/test_grid2par.f90
+++ b/unit-tests/3d/test_grid2par.f90
@@ -6,12 +6,11 @@
 program test_trilinear
     use unit_test
     use constants, only : pi, zero, one, two, three, four, five, f12, f23
-    use physics, only : ft_cor, f_cor
     use parcel_container
     use parcel_interpl, only : grid2par, grid2par_timer
     use parcel_ellipsoid, only : get_abc
     use parameters, only : lower, update_parameters, vcell, dx, nx, ny, nz
-    use fields, only : velog, vortg, velgradg, dbdx, dbdy, field_alloc
+    use fields, only : velog, vortg, velgradg, field_alloc
     use timer
     implicit none
 
@@ -81,9 +80,6 @@ program test_trilinear
     vortg(:, :, :, 2) = two
     vortg(:, :, :, 3) = three
 
-    dbdx = one
-    dbdy = two
-
     do l = 1, 5
         velgradg(:, :, :, l) = dble(l)
     enddo
diff --git a/unit-tests/3d/test_lapinv0_1.f90 b/unit-tests/3d/test_lapinv0_1.f90
index 483cb4ade..2f0e6ebc4 100644
--- a/unit-tests/3d/test_lapinv0_1.f90
+++ b/unit-tests/3d/test_lapinv0_1.f90
@@ -64,7 +64,7 @@ program test_lapinv0_1
 
     error = maxval(dabs(fp - ref_sol))
 
-    call print_result_dp('Test inversion (lapinv0)', error, atol=4.0e-12)
+    call print_result_dp('Test inversion (lapinv0)', error, atol=2.0e-7)
 
     deallocate(fs)
     deallocate(fp)
diff --git a/unit-tests/3d/test_lapinv0_2.f90 b/unit-tests/3d/test_lapinv0_2.f90
index babc5f287..f0a88e21e 100644
--- a/unit-tests/3d/test_lapinv0_2.f90
+++ b/unit-tests/3d/test_lapinv0_2.f90
@@ -1,13 +1,7 @@
 ! =============================================================================
 !                       Test subroutine lapinv0
 !
-!  This unit test checks the subroutine lapinv0 using the the
-!  function:
-!              cos(k * x) * sin(l * y) * sin(m * z)
-!  where k = 2pi/L_x, l = 2pi/L_y and m = pi/L_z and where x, y and z all start
-!  at 0 (one could start at -pi for x and y just as well).
-!  The subroutine lapinv0 should return the same function multiplied by
-!  -1/(k^2 + l^2 + m^2).
+!  This unit test checks the subroutine lapinv0.
 ! =============================================================================
 program test_lapinv0_2
     use unit_test
@@ -90,7 +84,7 @@ program test_lapinv0_2
     endif
 
 
-    call print_result_dp('Test inversion (lapinv0)', emax, atol=2.0e-11)
+    call print_result_dp('Test inversion (lapinv0)', emax, atol=7.0e-7)
 
     deallocate(vorx)
     deallocate(potx)
diff --git a/unit-tests/3d/test_lapinv1.f90 b/unit-tests/3d/test_lapinv1.f90
new file mode 100644
index 000000000..0380c0784
--- /dev/null
+++ b/unit-tests/3d/test_lapinv1.f90
@@ -0,0 +1,73 @@
+! =============================================================================
+!                       Test subroutine lapinv1
+!
+!  This unit test checks the subroutine lapinv1 using the the
+!  function:
+!              cos(k * x) * sin(l * y) * cos(m * z)
+!  where k = 2pi/L_x, l = 2pi/L_y and m = pi/L_z and where x, y and z all start
+!  at 0 (one could start at -pi for x and y just as well).
+!  The subroutine lapinv1 should return the same function multiplied by
+!  -1/(k^2 + l^2 + m^2).
+! =============================================================================
+program test_lapinv1
+    use unit_test
+    use constants, only : zero, one, two, pi, twopi
+    use parameters, only : lower, update_parameters, dx, nx, ny, nz, extent
+    use inversion_utils, only : init_fft, fftxyp2s, fftxys2p
+    use inversion_mod, only : lapinv1
+    implicit none
+
+    double precision              :: error
+    double precision, allocatable :: fp(:, :, :), fs(:, :, :), ref_sol(:, :, :)
+    integer                       :: ix, iy, iz
+    double precision              :: x, y, z, k, l, m, prefactor
+
+    nx = 32
+    ny = 64
+    nz = 128
+    lower  = (/zero, zero, zero/)
+    extent =  (/pi, twopi, two * twopi/)
+
+    allocate(fs(0:nz, nx, ny))
+    allocate(fp(0:nz, ny, nx))
+    allocate(ref_sol(0:nz, ny, nx))
+
+    call update_parameters
+
+    k = twopi / extent(1)
+    l = twopi / extent(2)
+    m =    pi / extent(3)
+
+    prefactor = - one / (k ** 2 + l ** 2 + m ** 2)
+
+    do ix = 1, nx
+        x = lower(1) + (ix - 1) * dx(1)
+        do iy = 1, ny
+            y = lower(2) + (iy - 1) * dx(2)
+            do iz = 0, nz
+                z = lower(3) + iz * dx(3)
+
+                fp(iz, iy, ix) = dcos(k * x) * dsin(l * y) * dcos(m * z)
+
+                ref_sol(iz, iy, ix) = prefactor * fp(iz, iy, ix)
+            enddo
+        enddo
+    enddo
+
+    call init_fft
+
+    call fftxyp2s(fp, fs)
+
+    call lapinv1(fs)
+
+    call fftxys2p(fs, fp)
+
+    error = maxval(dabs(fp - ref_sol))
+
+    call print_result_dp('Test inversion (lapinv1)', error, atol=2.0e-7)
+
+    deallocate(fs)
+    deallocate(fp)
+    deallocate(ref_sol)
+
+end program test_lapinv1