diff --git a/Makefile b/Makefile
index b63a53f20e44f3a2b4000053d8bd543c16a0d5bd..e33d4aa587940bb7e81682b50079cc55d03fa3b9 100644
--- a/Makefile
+++ b/Makefile
@@ -55,7 +55,7 @@ SRCDECOMP = $(DECOMPDIR)/decomp_2d.f90 $(DECOMPDIR)/glassman.f90 $(DECOMPDIR)/ff
 OBJDECOMP = $(SRCDECOMP:%.f90=%.o)
 SRC = $(SRCDIR)/module_param.f90 $(SRCDIR)/variables.f90 $(SRCDIR)/poisson.f90 $(SRCDIR)/derive.f90 $(SRCDIR)/implicit.f90 $(SRCDIR)/schemes.f90 $(SRCDIR)/parameters.f90 $(SRCDIR)/*.f90
 OBJ = $(SRC:%.f90=%.o)
-SRC = $(SRCDIR)/module_param.f90 $(SRCDIR)/variables.f90 $(SRCDIR)/BC-dbg-schemes.f90 $(SRCDIR)/poisson.f90 $(SRCDIR)/ibm.f90 $(SRCDIR)/derive.f90 $(SRCDIR)/implicit.f90 $(SRCDIR)/schemes.f90 $(SRCDIR)/forces.f90 $(SRCDIR)/probes.f90 $(SRCDIR)/navier.f90 $(SRCDIR)/tools.f90 $(SRCDIR)/visu.f90 $(SRCDIR)/BC-TBL.f90 $(SRCDIR)/BC-ABL.f90 $(SRCDIR)/les_models.f90 $(SRCDIR)/BC-Lock-exchange.f90 $(SRCDIR)/time_integrators.f90 $(SRCDIR)/filters.f90 $(SRCDIR)/parameters.f90 $(SRCDIR)/BC-User.f90 $(SRCDIR)/BC-TGV.f90 $(SRCDIR)/BC-Channel-flow.f90 $(SRCDIR)/BC-Periodic-hill.f90 $(SRCDIR)/BC-Cylinder.f90 $(SRCDIR)/BC-Mixing-layer.f90 $(SRCDIR)/BC-Jet.f90 $(SRCDIR)/BC-Sandbox.f90 $(SRCDIR)/BC-Uniform.f90 $(TURBDIR)/constants.f90 $(TURBDIR)/acl_utils.f90 $(TURBDIR)/airfoils.f90 $(TURBDIR)/dynstall.f90 $(TURBDIR)/dynstall_legacy.f90 $(TURBDIR)/acl_elem.f90 $(TURBDIR)/acl_controller.f90 $(TURBDIR)/acl_turb.f90 $(TURBDIR)/acl_out.f90 $(TURBDIR)/acl_farm_controller.f90 $(TURBDIR)/acl_model.f90 $(TURBDIR)/acl_source.f90 $(TURBDIR)/adm.f90 $(TURBDIR)/turbine.f90 $(SRCDIR)/statistics.f90 $(SRCDIR)/case.f90 $(SRCDIR)/transeq.f90 $(SRCDIR)/genepsi3d.f90 $(SRCDIR)/xcompact3d.f90
+SRC = $(SRCDIR)/module_param.f90 $(SRCDIR)/variables.f90 $(SRCDIR)/BC-dbg-schemes.f90 $(SRCDIR)/poisson.f90 $(SRCDIR)/ibm.f90 $(SRCDIR)/derive.f90 $(SRCDIR)/implicit.f90 $(SRCDIR)/schemes.f90 $(SRCDIR)/forces.f90 $(SRCDIR)/probes.f90 $(SRCDIR)/navier.f90 $(SRCDIR)/tools.f90 $(SRCDIR)/visu.f90 $(SRCDIR)/BC-TBL.f90 $(SRCDIR)/BC-ABL.f90 $(SRCDIR)/les_models.f90 $(SRCDIR)/BC-Lock-exchange.f90 $(SRCDIR)/time_integrators.f90 $(SRCDIR)/filters.f90 $(SRCDIR)/parameters.f90 $(SRCDIR)/BC-User.f90 $(SRCDIR)/BC-TGV.f90 $(SRCDIR)/BC-Channel-flow.f90 $(SRCDIR)/BC-Periodic-hill.f90 $(SRCDIR)/BC-Cylinder.f90 $(SRCDIR)/BC-Mixing-layer.f90 $(SRCDIR)/BC-Jet.f90 $(SRCDIR)/BC-Sandbox.f90 $(SRCDIR)/BC-Uniform.f90 $(SRCDIR)/BC-Cavity.f90 $(TURBDIR)/constants.f90 $(TURBDIR)/acl_utils.f90 $(TURBDIR)/airfoils.f90 $(TURBDIR)/dynstall.f90 $(TURBDIR)/dynstall_legacy.f90 $(TURBDIR)/acl_elem.f90 $(TURBDIR)/acl_controller.f90 $(TURBDIR)/acl_turb.f90 $(TURBDIR)/acl_out.f90 $(TURBDIR)/acl_farm_controller.f90 $(TURBDIR)/acl_model.f90 $(TURBDIR)/acl_source.f90 $(TURBDIR)/adm.f90 $(TURBDIR)/turbine.f90 $(SRCDIR)/statistics.f90 $(SRCDIR)/case.f90 $(SRCDIR)/transeq.f90 $(SRCDIR)/genepsi3d.f90 $(SRCDIR)/xcompact3d.f90
 
 
 #######FFT settings##########
diff --git a/examples/Cavity/input.i3d b/examples/Cavity/input.i3d
new file mode 100644
index 0000000000000000000000000000000000000000..c8669b7301dd227ede9efbdcccb3bf1663d8e49c
--- /dev/null
+++ b/examples/Cavity/input.i3d
@@ -0,0 +1,126 @@
+! -*- mode: f90 -*-
+
+!===================
+&BasicParam
+!===================
+
+! Flow type (1=Lock-exchange, 2=TGV, 3=Channel, 4=Periodic hill, 5=Cylinder, 6=dbg-schemes)
+itype = 13
+
+! Domain decomposition
+p_row=0               ! Row partition
+p_col=0               ! Column partition
+
+! Mesh
+nx=257                ! X-direction nodes
+ny=129                ! Y-direction nodes
+nz=1                  ! Z-direction nodes
+istret = 0            ! y mesh refinement (0:no, 1:center, 2:both sides, 3:bottom)
+
+! Domain
+xlx = 1.              ! Lx (Size of the box in x-direction)
+yly = 1.              ! Ly (Size of the box in y-direction)
+zlz = 1.              ! Lz (Size of the box in z-direction)
+
+! Flow parameters
+re  = 14084.5         ! nu=1/re (Kinematic Viscosity)
+
+! Time stepping
+dt = 0.002            ! Time step
+ifirst = 1            ! First iteration
+ilast = 5000          ! Last iteration
+
+! Enable modelling tools
+numscalar = 1         ! How many scalars? (Set to zero to disable scalars)
+
+! Velocity boundary conditions
+nclx1 = 2
+nclxn = 2
+ncly1 = 2
+nclyn = 2
+nclz1 = 0
+nclzn = 0
+
+! Gravity vector
+gravx = 0.
+gravy = 1.
+gravz = 0.
+
+/End
+
+!====================
+&NumOptions
+!====================
+
+! Spatial derivatives
+ifirstder = 4         ! (1->2nd central, 2->4th central, 3->4th compact, 4-> 6th compact)
+isecondder = 4        ! (1->2nd central, 2->4th central, 3->4th compact, 4-> 6th compact, 5->hyperviscous 6th)
+ipinter = 1           ! interpolation scheme (1: classic, 2: optimized, 3: optimized agressive)
+
+! Time scheme
+itimescheme = 2       ! Time integration scheme (1->Euler,2->AB2, 3->AB3, 4->AB4,5->RK3,6->RK4, 7-->CN2+AB3)
+
+! Dissipation control
+nu0nu = 4.0             ! Ratio between hyperviscosity/viscosity at nu
+cnu = 0.44              ! Ratio between hypervisvosity at k_m=2/3pi and k_c= pi
+
+/End
+
+!=================
+&InOutParam
+!=================
+
+! Basic I/O
+irestart = 0          ! Read initial flow field ?
+icheckpoint = 400000  ! Frequency for writing backup file
+ioutput = 500         ! Frequency for visualization
+ilist = 25            ! Frequency for writing to screen        
+nvisu = 1             ! Size for visualisation collection
+nprobes = 1           ! Number of probes
+
+/End
+
+!=================
+&Statistics
+!=================
+
+nstat = 1             ! Size arrays for statistic collection
+initstat = 300000     ! Time steps after which statistics are collected
+
+/End
+
+&CASE
+/End
+
+!=================
+&ProbesParam
+!=================
+
+flag_all_digits = F     ! By default, only 6 digits, 16 when True
+flag_extra_probes = F   ! By default, gradients are not monitored
+xyzprobes(1,1) = 0.1    ! X position of probe 1
+xyzprobes(2,1) = 0.1    ! Y position of probe 1
+xyzprobes(3,1) = 0.1    ! Z position of probe 1
+
+/End
+
+!########################
+! OPTIONAL PARAMETERS
+!#######################
+
+!================
+&ScalarParam
+!================
+
+Sc(1) = 0.71          ! Schmidt number
+sc_skew(1) = F        ! Skew-symetric convection
+ri(1) = 0.71          ! Richardson number
+
+nclxS1 = 2
+nclxSn = 2
+nclyS1 = 2
+nclySn = 2
+nclzS1 = 0
+nclzSn = 0
+
+/End
diff --git a/examples/Taylor-Green-Vortex/input.i3d b/examples/Taylor-Green-Vortex/input.i3d
index 4c352db87cbe9a7edd5a30deacf9382d50d1f821..a9aeba05aa0f6a15695526271f7e9a76d9842c79 100644
--- a/examples/Taylor-Green-Vortex/input.i3d
+++ b/examples/Taylor-Green-Vortex/input.i3d
@@ -105,5 +105,4 @@ iwall=0
 !=================
 &CASE
 !=================
-tgv_twod = .FALSE.
 /End
diff --git a/examples/Taylor-Green-Vortex/input_2D.i3d b/examples/Taylor-Green-Vortex/input_2D.i3d
new file mode 100644
index 0000000000000000000000000000000000000000..b15f20ab7489c79af71bb3b19f7841e432bdc048
--- /dev/null
+++ b/examples/Taylor-Green-Vortex/input_2D.i3d
@@ -0,0 +1,108 @@
+! -*- mode: f90 -*-
+
+!===================
+&BasicParam
+!===================
+
+! Flow type (1=Lock-exchange, 2=TGV, 3=Channel, 4=Periodic hill, 5=Cylinder, 6=dbg-schemes)
+itype = 2
+
+! Domain decomposition
+p_row=0               ! Row partition
+p_col=0               ! Column partition
+
+! Mesh
+nx=65                 ! X-direction nodes
+ny=65                 ! Y-direction nodes
+nz=1                  ! Z-direction nodes
+istret = 0            ! y mesh refinement (0:no, 1:center, 2:both sides, 3:bottom)
+beta = 0.259065151    ! Refinement parameter (beta)
+
+! Domain
+xlx = 3.14159265358979 ! Lx (Size of the box in x-direction)
+yly = 3.14159265358979 ! Ly (Size of the boy in y-direction)
+zlz = 3.14159265358979 ! Lz (Size of the boz in z-direction)
+
+! Boundary conditions
+nclx1 = 1 
+nclxn = 1 
+ncly1 = 1 
+nclyn = 1 
+nclz1 = 1 
+nclzn = 1 
+
+! Flow parameters
+iin = 1               ! Inflow conditions (1: classic, 2: turbinit)
+re = 1600.            ! nu=1/re (Kinematic Viscosity)
+u1 = 8.               ! u1 (max velocity) (for inflow condition)
+u2 = 8.               ! u2 (min velocity) (for inflow condition)
+init_noise  = 0.0     ! Turbulence intensity (1=100%) !! Initial condition
+inflow_noise = 0.0    ! Turbulence intensity (1=100%) !! Inflow condition
+
+! Time stepping
+dt = 0.005            ! Time step
+ifirst = 1            ! First iteration
+ilast = 4000         ! Last iteration
+
+! Enable modelling tools
+ilesmod=0             ! if 0 then DNS
+iscalar=0             ! If iscalar=0 (no scalar), if iscalar=1 (scalar)
+iibm=0                ! Flag for immersed boundary method
+
+! Enable io
+ivisu=1               ! Store snapshots
+ipost=1               ! Do online postprocessing
+
+/End
+
+!====================
+&NumOptions
+!====================
+
+! Spatial derivatives
+ifirstder = 4         ! (1->2nd central, 2->4th central, 3->4th compact, 4-> 6th compact)
+isecondder = 4        ! (1->2nd central, 2->4th central, 3->4th compact, 4-> 6th compact, 5->hyperviscous 6th)
+
+! Time scheme
+itimescheme = 5       ! Time integration scheme (1->Euler,2->AB2, 3->AB3, 4->AB4,5->RK3,6->RK4)
+
+/End
+
+!=================
+&InOutParam
+!=================
+
+! Basic I/O
+irestart = 0          ! Read initial flow field ?
+icheckpoint = 1000    ! Frequency for writing backup file
+ioutput = 1000        ! Frequency for visualization
+ilist = 5             ! Frequency for the output
+nvisu = 1             ! Size for visualisation collection
+
+/End
+
+!=================
+&Statistics
+!=================
+
+nstat = 1             ! Size arrays for statistic collection
+
+/End
+
+!#######################
+! OPTIONAL PARAMETERS
+!#######################
+
+!=================
+&LESModel
+!=================
+jles=1
+smagcst=0.1
+walecst=0.5 
+iwall=0
+/End
+
+!=================
+&CASE
+!=================
+/End
diff --git a/examples/Taylor-Green-Vortex/input_DNS_Re1600.i3d b/examples/Taylor-Green-Vortex/input_DNS_Re1600.i3d
index e60d1b829cc0a0f52ecd015c4363d8652f6f6642..b09d0befbf0ae259553fd07d4a052c8f041ba591 100644
--- a/examples/Taylor-Green-Vortex/input_DNS_Re1600.i3d
+++ b/examples/Taylor-Green-Vortex/input_DNS_Re1600.i3d
@@ -105,5 +105,4 @@ iwall=0
 !=================
 &CASE
 !=================
-tgv_twod = .FALSE.
 /End
diff --git a/examples/Taylor-Green-Vortex/input_ILES_Re5000.i3d b/examples/Taylor-Green-Vortex/input_ILES_Re5000.i3d
index 5feb78efe099728f3d6de092593bd7391208950f..eb0d8e4d8179b0d8692246739f24c5ed20df7568 100644
--- a/examples/Taylor-Green-Vortex/input_ILES_Re5000.i3d
+++ b/examples/Taylor-Green-Vortex/input_ILES_Re5000.i3d
@@ -109,5 +109,4 @@ iwall=0
 !=================
 &CASE
 !=================
-tgv_twod = .FALSE.
 /End
diff --git a/examples/Taylor-Green-Vortex/input_test.i3d b/examples/Taylor-Green-Vortex/input_test.i3d
index 1c784912358b129fc3dfbbabb5d70e1003920e49..03eec3eb4d0990de637bcb11f18d607705ab4de5 100644
--- a/examples/Taylor-Green-Vortex/input_test.i3d
+++ b/examples/Taylor-Green-Vortex/input_test.i3d
@@ -105,5 +105,4 @@ iwall=0
 !=================
 &CASE
 !=================
-tgv_twod = .FALSE.
 /End
diff --git a/src/BC-Cavity.f90 b/src/BC-Cavity.f90
new file mode 100644
index 0000000000000000000000000000000000000000..9cf189fb3d516d5a33334c099696b45a6b791219
--- /dev/null
+++ b/src/BC-Cavity.f90
@@ -0,0 +1,139 @@
+!Copyright (c) 2012-2022, Xcompact3d
+!This file is part of Xcompact3d (xcompact3d.com)
+!SPDX-License-Identifier: BSD 3-Clause
+
+module cavity
+
+   use decomp_2d
+   use variables
+   use param
+
+   implicit none
+
+   ! Temperature difference between both walls
+   real(mytype), parameter :: deltaT = 1._mytype
+   real(mytype), parameter :: temp_l = deltaT/2._mytype
+
+   private ! All functions/subroutines private by default
+   public :: init_cavity, boundary_conditions_cavity, postprocess_cavity
+
+contains
+
+   !############################################################################
+   !!
+   !! Init the cavity case
+   !!
+   !############################################################################
+   subroutine init_cavity(ux1, uy1, uz1, ep1, phi1)
+
+      USE decomp_2d_io
+      USE MPI
+
+      implicit none
+
+      real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: ux1, uy1, uz1, ep1
+      real(mytype), dimension(xsize(1), xsize(2), xsize(3), numscalar) :: phi1
+
+      real(mytype) :: x
+      integer :: i
+
+      ! This does not apply in case of restart
+      if (irestart == 0) then
+
+         ! Velocity is zero
+         ux1 = zero
+         uy1 = zero
+         uz1 = zero
+
+         ! Linear temperature profile
+         if (numscalar >= 1) then
+            do i = 1, xsize(1)
+               phi1(i, :, :, :) = temp_l - deltaT*(i - 1)/real(xsize(1) - 1, kind=mytype)
+            end do
+         end if
+
+      end if
+
+   end subroutine init_cavity
+
+   !############################################################################
+   !!
+   !! Boundary conditions for the cavity case
+   !!
+   !############################################################################
+   subroutine boundary_conditions_cavity(ux, uy, uz, phi)
+
+      implicit none
+
+      real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: ux, uy, uz
+      real(mytype), dimension(xsize(1), xsize(2), xsize(3), numscalar) :: phi
+
+      integer :: i, j, k
+
+      ! Velocity
+      IF (nclx1 == 2) THEN
+      END IF
+      IF (nclxn == 2) THEN
+      END IF
+      IF (ncly1 == 2) THEN
+      END IF
+      IF (nclyn == 2) THEN
+      END IF
+      IF (nclz1 == 2) THEN
+      END IF
+      IF (nclzn == 2) THEN
+      END IF
+
+      ! Scalar
+      if (numscalar >= 1) then
+         ! x=0
+         if (nclxS1 == 2) then
+            phi(1, :, :, :) = temp_l
+         end if
+         ! x=Lx
+         if (nclxSn == 2) then
+            phi(xsize(1), :, :, :) = temp_l - deltaT
+         end if
+         ! y=0
+         if (nclyS1 == 2 .and. xstart(2) == 1) then
+            do i = 1, xsize(1)
+               phi(i, 1, :, :) = temp_l - deltaT*(i - 1)/real(xsize(1) - 1, kind=mytype)
+            end do
+         end if
+         ! y=Ly
+         if (nclySn == 2 .and. xend(2) == ny) then
+            do i = 1, xsize(1)
+               phi(i, xsize(2), :, :) = temp_l - deltaT*(i - 1)/real(xsize(1) - 1, kind=mytype)
+            end do
+         end if
+      end if
+
+      ! Clip
+      if (numscalar >= 1) then
+         do k = 1, xsize(3)
+            do j = 1, xsize(2)
+               do i = 1, xsize(1)
+                  if (phi(i,j,k,1) > temp_l) phi(i,j,k,1) = temp_l
+                  if (phi(i,j,k,1) < -temp_l) phi(i,j,k,1) = -temp_l
+               enddo
+            enddo
+         enddo
+      endif
+
+   end subroutine boundary_conditions_cavity
+
+   !############################################################################
+   !!
+   !! Post-processing for the cavity case
+   !!
+   !############################################################################
+   subroutine postprocess_cavity(ux1, uy1, uz1, phi1)
+
+      implicit none
+
+      real(mytype), intent(in), dimension(xsize(1), xsize(2), xsize(3)) :: ux1, uy1, uz1
+      real(mytype), intent(in), dimension(xsize(1), xsize(2), xsize(3), numscalar) :: phi1
+
+   end subroutine postprocess_cavity
+
+end module cavity
diff --git a/src/BC-TGV.f90 b/src/BC-TGV.f90
index fb81d02d9d2af700f1b8c2ef4372a75b58d9089d..f0963365cada1ce0feb6bd14bb76eb28fcfca66b 100644
--- a/src/BC-TGV.f90
+++ b/src/BC-TGV.f90
@@ -66,18 +66,10 @@ contains
              do i=1,xsize(1)
                 x=real(i-1,mytype)*dx
 
-                if (.not.tgv_twod) then
-                   ux1(i,j,k)=+sin_prec(x)*cos_prec(y)*cos_prec(z)
-                   uy1(i,j,k)=-cos_prec(x)*sin_prec(y)*cos_prec(z)
-                   if (iscalar == 1) then
-                      phi1(i,j,k,1:numscalar)=sin_prec(x)*sin_prec(y)*cos_prec(z)
-                   endif
-                else
-                   ux1(i,j,k)=+sin_prec(x)*cos_prec(y)
-                   uy1(i,j,k)=-cos_prec(x)*sin_prec(y)
-                   if (iscalar == 1) then
-                      phi1(i,j,k,1:numscalar)=sin_prec(x)*sin_prec(y)
-                   endif
+                ux1(i,j,k)=+sin_prec(x)*cos_prec(y)*cos_prec(z)
+                uy1(i,j,k)=-cos_prec(x)*sin_prec(y)*cos_prec(z)
+                if (iscalar == 1) then
+                   phi1(i,j,k,1:numscalar)=sin_prec(x)*sin_prec(y)*cos_prec(z)
                 endif
                 uz1(i,j,k)=zero
              enddo
diff --git a/src/case.f90 b/src/case.f90
index 1656f798c75529698996f7065632cec23677bd56..544c314fe9ede685aeb4ebc4763b1996d59ded85 100644
--- a/src/case.f90
+++ b/src/case.f90
@@ -21,6 +21,7 @@ module case
   use abl
   use uniform
   use sandbox
+  use cavity
 
   use var, only : nzmsize
 
@@ -110,6 +111,10 @@ contains
    
        call init_sandbox (ux1, uy1, uz1, ep1, phi1, 0)
 
+    elseif (itype.eq.itype_cavity) then
+
+       call init_cavity(ux1, uy1, uz1, ep1, phi1)
+
     else
   
          if (nrank.eq.0) then
@@ -194,6 +199,10 @@ contains
    
        call boundary_conditions_sandbox (ux, uy, uz, phi)
 
+    elseif (itype.eq.itype_cavity) then
+
+       call boundary_conditions_cavity(ux, uy, uz, phi)
+
     endif
 
   end subroutine boundary_conditions
@@ -339,6 +348,10 @@ contains
    
        call postprocess_sandbox (ux, uy, uz, phi, ep)
 
+    elseif (itype.eq.itype_cavity) then
+
+       call postprocess_cavity(ux, uy, uz, phi)
+
     endif
 
     if (iforces.eq.1) then
@@ -431,7 +444,7 @@ contains
        call visu_tbl(ux1, uy1, uz1, pp3, phi1, ep1, num)
        called_visu = .true.
        
-   elseif (itype.eq.itype_uniform) then
+    elseif (itype.eq.itype_uniform) then
 
        call visu_uniform(ux1, uy1, uz1, pp3, phi1, ep1, num)
        called_visu = .true.
diff --git a/src/derive.f90 b/src/derive.f90
index 5226509eb845a2f3fecf78ff675eb3d8c5c2838b..5ba54cddc9ac9ced2e48330c9621c7cdca20a140 100644
--- a/src/derive.f90
+++ b/src/derive.f90
@@ -871,6 +871,11 @@ subroutine derz_00(tz,uz,rz,sz,ffz,fsz,fwz,nx,ny,nz,npaire,lind)
   real(mytype), dimension(nz) :: ffz,fsz,fwz
   real(mytype)                      :: lind
 
+  if (nz==1) then
+     tz = 0.
+     return
+  endif
+
   if (iibm.eq.2) call lagpolz(uz)
   if (iibm.eq.3) call cubsplz(uz,lind)
 
@@ -2931,6 +2936,11 @@ subroutine derzz_00(tz,uz,rz,sz,sfz,ssz,swz,nx,ny,nz,npaire,lind)
   real(mytype), dimension(nz) :: sfz,ssz,swz
   real(mytype)                      :: lind
 
+  if (nz==1) then
+     tz = 0.
+     return
+  endif
+
   if (iibm.eq.2) call lagpolz(uz)
   if (iibm.eq.3) call cubsplz(uz,lind)
 
@@ -4924,6 +4934,11 @@ subroutine derzvp(tz,uz,rz,sz,cfz6,csz6,cwz6,nx,ny,nz,nzm,npaire)
   real(mytype), dimension(nzm) :: cfz6,csz6,cwz6
   integer :: i,j,k
 
+  if (nz==1) then
+     tz = 0.
+     return
+  endif
+
   if (nclz) then
      do j=1,ny
         do i=1,nx
@@ -5103,6 +5118,11 @@ subroutine interzvp(tz,uz,rz,sz,cifz6,cisz6,ciwz6,nx,ny,nz,nzm,npaire)
   real(mytype), dimension(nzm) :: cifz6,cisz6,ciwz6
   integer :: i,j,k
 
+  if (nz==1 .and. nzm==1) then
+     tz = uz
+     return
+  endif
+
   if (nclz) then
      do j=1,ny
         do i=1,nx
@@ -5282,6 +5302,11 @@ subroutine derzpv(tz,uz,rz,sz,cfiz6,csiz6,cwiz6,cfz6,csz6,cwz6,&
   real(mytype), dimension(nz) :: cfz6,csz6,cwz6
   integer :: i,j,k
 
+  if (nz==1) then
+     tz = 0.
+     return
+  endif
+
   if (nclz) then
      do j=1,ny
         do i=1,nx
@@ -5417,6 +5442,11 @@ subroutine interzpv(tz,uz,rz,sz,cifiz6,cisiz6,ciwiz6,cifz6,cisz6,ciwz6,&
   real(mytype), dimension(nz) :: cifz6,cisz6,ciwz6
   integer :: i,j,k
 
+  if (nz==1 .and. nzm==1) then
+     tz = uz
+     return
+  endif
+
   if (nclz) then
      do j=1,ny
         do i=1,nx
diff --git a/src/filters.f90 b/src/filters.f90
index 865d7ab8701a57ff8d8812aa5cdc4130319f6845..d71e96fe03bb4410a48d1b07fa900b068d2d86df 100644
--- a/src/filters.f90
+++ b/src/filters.f90
@@ -1004,6 +1004,11 @@ subroutine filz_00(tz,uz,rz,fisz,fiffz,fifsz,fifwz,nx,ny,nz,npaire,lind)
   real(mytype), dimension(nz) :: fiffz,fifsz,fifwz
   real(mytype)                      :: lind
 
+  if (nz==1) then
+    tz = uz
+    return
+  endif
+
   if (iibm.eq.2) call lagpolz(uz)
   if (iibm.eq.3) call cubsplz(uz,lind)
 
diff --git a/src/module_param.f90 b/src/module_param.f90
index a07c6da5c7a23632752ee7d41a9f1a090d1847b3..49b0303cbdcb2775dbc39bee594246524e67924a 100644
--- a/src/module_param.f90
+++ b/src/module_param.f90
@@ -283,7 +283,8 @@ module param
        itype_tbl = 9, &
        itype_abl = 10, &
        itype_uniform = 11, &
-       itype_sandbox = 12
+       itype_sandbox = 12, &
+       itype_cavity = 13
 
   integer :: cont_phi,itr,itime,itest,iprocessing
   integer :: ifft,istret,iforc_entree,iturb
@@ -374,9 +375,6 @@ module param
   character, dimension(100) :: TurbinesPath*80, ActuatorlinesPath*80
   real(mytype) :: eps_factor ! Smoothing factor
   
-  !! Case-specific variables
-  logical :: tgv_twod
-
   character :: filesauve*80, filenoise*80, &
        nchamp*80,filepath*80, fileturb*80, filevisu*80, datapath*80
   real(mytype), dimension(5) :: adt,bdt,cdt,ddt,gdt
diff --git a/src/parameters.f90 b/src/parameters.f90
index 5ad683600841314c408ecf8202f81365af8acce0..ab06f82f9f5f374c8ecc23904c3143e67e40a7a1 100644
--- a/src/parameters.f90
+++ b/src/parameters.f90
@@ -68,7 +68,7 @@ subroutine parameter(input_i3d)
   NAMELIST /ABL/ z_zero, iwallmodel, k_roughness, ustar, dBL, &
        imassconserve, ibuoyancy, iPressureGradient, iCoriolis, CoriolisFreq, &
        istrat, idamping, iheight, TempRate, TempFlux, itherm, gravv, UG, T_wall, T_top, ishiftedper, iconcprec, pdl 
-  NAMELIST /CASE/ tgv_twod, pfront
+  NAMELIST /CASE/ pfront
   NAMELIST/ALMParam/iturboutput,NTurbines,TurbinesPath,NActuatorlines,ActuatorlinesPath,eps_factor,rho_air
   NAMELIST/ADMParam/Ndiscs,ADMcoords,C_T,aind,iturboutput,rho_air
 
@@ -98,6 +98,13 @@ subroutine parameter(input_i3d)
 
   !! These are the 'essential' parameters
   read(10, nml=BasicParam); rewind(10)
+  if (nz==1) then
+     if (nrank==0) write(*,*) "Warning : support for 2D simulations is experimental"
+     nclz1 = 0
+     nclzn = 0
+     p_row = nproc
+     p_col = 1
+  endif
   read(10, nml=NumOptions); rewind(10)
   read(10, nml=InOutParam); rewind(10)
   read(10, nml=Statistics); rewind(10)
@@ -184,6 +191,10 @@ subroutine parameter(input_i3d)
   endif
   if (numscalar.ne.0) then
      read(10, nml=ScalarParam); rewind(10)
+     if (nz==1) then
+        nclzS1 = 0
+        nclzSn = 0
+     endif
   endif
   ! !! These are the 'optional'/model parameters
   ! read(10, nml=ScalarParam)
@@ -333,7 +344,9 @@ subroutine parameter(input_i3d)
      elseif (itype.eq.itype_uniform) then
         print *,'Uniform flow'
      elseif (itype.eq.itype_sandbox) then
-           print *,'Sandbox'
+        print *,'Sandbox'
+     elseif (itype.eq.itype_cavity) then
+        print *,'Cavity'  
      else
         print *,'Unknown itype: ', itype
         stop
@@ -523,8 +536,6 @@ subroutine parameter(input_i3d)
      !! Print case-specific information
      if (itype==itype_lockexch) then
         write(*,*)  "Initial front location: ", pfront
-     elseif (itype==itype_tgv) then
-        write(*,*)  "TGV 2D: ", tgv_twod
      endif
      write(*,*) '==========================================================='
   endif
@@ -673,9 +684,6 @@ subroutine parameter_defaults()
 
   imodulo2 = 1
 
-  !! CASE specific variables
-  tgv_twod = .FALSE.
-
   !! TRIPPING
   A_tr=zero
   xs_tr_tbl=1.402033_mytype
diff --git a/src/schemes.f90 b/src/schemes.f90
index 60d411e20e2d442051ed72b1955d8812d23fe768..1891baf8269c541f2e0b9dd5d44232b7a4d8fe10 100644
--- a/src/schemes.f90
+++ b/src/schemes.f90
@@ -313,6 +313,8 @@ subroutine first_derivative(alfa1,af1,bf1,cf1,df1,alfa2,af2,alfan,afn,bfn,&
      afm  = (three/four)/d
   endif
 
+  if (n==1) return
+
   if     (ncl1.eq.0) then !Periodic
      ff(1)   =alfai
      ff(2)   =alfai
@@ -530,6 +532,8 @@ subroutine second_derivative(alsa1,as1,bs1,&
   bstt = (three/fortyfour)/d2
   cstt = zero
 
+  if (n==1) return
+
   if     (ncl1.eq.0) then !Periodic
      sf(1)   =alsai
      sf(2)   =alsai
@@ -687,6 +691,40 @@ subroutine interpolation(dx,nxm,nx,nclx1,nclxn,&
      bcix6=(17._mytype/62._mytype)/three/dx
   endif
 
+  if (ifirstder == 1) then
+     ailcaix6 = zero
+     aicix6   = half
+     bicix6   = zero
+     cicix6   = zero
+     dicix6   = zero
+  else if (ipinter.eq.1) then
+     ailcaix6=three/ten
+     aicix6=three/four
+     bicix6=one/(two*ten)
+     cicix6=zero
+     dicix6=zero
+  else if (ipinter.eq.2) then
+     ailcaix6=0.461658_mytype
+
+     dicix6=0.00293016_mytype
+     aicix6=one/64._mytype *(75._mytype +70._mytype *ailcaix6-320._mytype *dicix6)
+     bicix6=one/128._mytype *(126._mytype *ailcaix6-25._mytype +1152._mytype *dicix6)
+     cicix6=one/128._mytype *(-ten*ailcaix6+three-640._mytype *dicix6)
+
+     aicix6=aicix6/two
+     bicix6=bicix6/two
+     cicix6=cicix6/two
+     dicix6=dicix6/two
+  else if (ipinter.eq.3) then
+     ailcaix6=0.49_mytype
+     aicix6=one/128._mytype *(75._mytype +70._mytype*ailcaix6)
+     bicix6=one/256._mytype *(126._mytype*ailcaix6-25._mytype)
+     cicix6=one/256._mytype *(-ten*ailcaix6+three)
+     dicix6=zero
+  endif
+
+  if (nx==1) return
+
   cfx6(1)=alcaix6
   cfx6(2)=alcaix6
   cfx6(nxm-2)=alcaix6
@@ -733,38 +771,6 @@ subroutine interpolation(dx,nxm,nx,nclx1,nclxn,&
      cbi6(i)=alcaix6
   enddo
 
-  if (ifirstder == 1) then
-     ailcaix6 = zero
-     aicix6   = half
-     bicix6   = zero
-     cicix6   = zero
-     dicix6   = zero
-  else if (ipinter.eq.1) then
-     ailcaix6=three/ten
-     aicix6=three/four
-     bicix6=one/(two*ten)
-     cicix6=zero
-     dicix6=zero
-  else if (ipinter.eq.2) then
-     ailcaix6=0.461658_mytype
-
-     dicix6=0.00293016_mytype
-     aicix6=one/64._mytype *(75._mytype +70._mytype *ailcaix6-320._mytype *dicix6)
-     bicix6=one/128._mytype *(126._mytype *ailcaix6-25._mytype +1152._mytype *dicix6)
-     cicix6=one/128._mytype *(-ten*ailcaix6+three-640._mytype *dicix6)
-
-     aicix6=aicix6/two
-     bicix6=bicix6/two
-     cicix6=cicix6/two
-     dicix6=dicix6/two
-  else if (ipinter.eq.3) then
-     ailcaix6=0.49_mytype
-     aicix6=one/128._mytype *(75._mytype +70._mytype*ailcaix6)
-     bicix6=one/256._mytype *(126._mytype*ailcaix6-25._mytype)
-     cicix6=one/256._mytype *(-ten*ailcaix6+three)
-     dicix6=zero
-  endif
-
   cifx6(1)=ailcaix6
   cifx6(2)=ailcaix6
   cifx6(nxm-2)=ailcaix6