diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7a9a3da01ee573485423173ed313c9e6ef7f0438..752418b24cdba9fcbd15df44cf6f6c356a814080 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -57,6 +57,7 @@ set(SOURCES
     obl_run_kato.f90
     obl_run_papa_fluxes.f90
     obl_run_papa_meteo.f90
+    obl_run_cbl.f90
     obl_config.f90
     vermix_inmom.f90
 )
diff --git a/obl_config.f90 b/obl_config.f90
index 4c4adaf7defea79fcc9e8f8f2091cabcfcd6ebae..49adfc178c5dd6a193b4046b2fe08cd10d3f86cc 100644
--- a/obl_config.f90
+++ b/obl_config.f90
@@ -22,10 +22,12 @@ module obl_config
     integer, parameter :: obl_config_kato = 0             !< Kato-Phillips setup
     integer, parameter :: obl_config_papa_fluxes = 1      !< Papa-station (fluxes) setup
     integer, parameter :: obl_config_papa_meteo = 2       !< Papa-station (meteo) setup
+    integer, parameter :: obl_config_cbl = 3              !< CBL setup
 
     character(len = 16), parameter :: obl_config_kato_tag = 'kato'
     character(len = 16), parameter :: obl_config_papa_fluxes_tag = 'papa-fluxes'
     character(len = 16), parameter :: obl_config_papa_meteo_tag = 'papa-meteo'
+    character(len = 16), parameter :: obl_config_cbl_tag = 'cbl'
 
     !> @brief model enum def.
     integer, parameter :: obl_model_pph = 0             !< pacanowski-philander 
@@ -52,6 +54,8 @@ contains
             id = obl_config_papa_fluxes
         else if (trim(tag) == trim(obl_config_papa_meteo_tag)) then
             id = obl_config_papa_meteo
+        else if (trim(tag) == trim(obl_config_cbl_tag)) then
+            id = obl_config_cbl
         end if
 
     end function
@@ -68,6 +72,8 @@ contains
             tag = obl_config_papa_fluxes_tag
         else if (id == obl_config_papa_meteo) then
             tag = obl_config_papa_meteo_tag
+        else if (id == obl_config_cbl) then
+            tag = obl_config_cbl_tag
         end if 
 
     end function
@@ -113,6 +119,7 @@ contains
         use obl_run_kato, only : set_grid_kato => set_grid
         use obl_run_papa_fluxes, only : set_grid_papa_fluxes => set_grid
         use obl_run_papa_meteo, only : set_grid_papa_meteo => set_grid
+        use obl_run_cbl, only : set_grid_cbl => set_grid
 
         type (gridDataType), intent(inout) :: grid
         integer, intent(in) :: config_id
@@ -134,6 +141,10 @@ contains
             call set_grid_papa_meteo(grid)
             return
         endif
+        if (config_id == obl_config_cbl) then
+            call set_grid_cbl(grid)
+            return
+        endif
 
 #ifdef USE_CONFIG_PARSER
         block 
@@ -169,6 +180,7 @@ contains
         use obl_run_kato, only : set_time_kato => set_time
         use obl_run_papa_fluxes, only : set_time_papa_fluxes => set_time
         use obl_run_papa_meteo, only : set_time_papa_meteo => set_time
+        use obl_run_cbl, only : set_time_cbl => set_time
 
         real, intent(out) :: time_begin, time_end, dt
         integer, intent(in) :: config_id
@@ -190,6 +202,10 @@ contains
             call set_time_papa_meteo(time_begin, time_end, dt)
             return
         endif
+        if (config_id == obl_config_cbl) then
+            call set_time_cbl(time_begin, time_end, dt)
+            return
+        endif
 
 #ifdef USE_CONFIG_PARSER
         block 
@@ -228,6 +244,7 @@ contains
         use obl_run_kato, only : set_phys_kato => set_phys
         use obl_run_papa_fluxes, only : set_phys_papa_fluxes => set_phys
         use obl_run_papa_meteo, only : set_phys_papa_meteo => set_phys
+        use obl_run_cbl, only : set_phys_cbl => set_phys
 
         integer, intent(in) :: config_id
         integer, intent(out) :: ierr
@@ -247,6 +264,10 @@ contains
             call set_phys_papa_meteo
             return
         endif
+        if (config_id == obl_config_cbl) then
+            call set_phys_cbl
+            return
+        endif
 
 #ifdef USE_CONFIG_PARSER
         block 
@@ -298,6 +319,7 @@ contains
         use obl_run_kato, only : set_forcing_kato => set_forcing
         use obl_run_papa_fluxes, only : set_forcing_papa_fluxes => set_forcing
         use obl_run_papa_meteo, only : set_forcing_papa_meteo => set_forcing
+        use obl_run_cbl, only : set_forcing_cbl => set_forcing
 
         integer, intent(in) :: config_id
         integer, intent(out) :: ierr
@@ -317,6 +339,10 @@ contains
             call set_forcing_papa_meteo
             return
         endif
+        if (config_id == obl_config_cbl) then
+            call set_forcing_cbl
+            return
+        endif
 
         !< assuming that forcing def. is optional
         block
@@ -400,6 +426,7 @@ contains
         use obl_run_kato, only : set_initial_conditions_kato => set_initial_conditions
         use obl_run_papa_fluxes, only : set_initial_conditions_papa_fluxes => set_initial_conditions
         use obl_run_papa_meteo, only : set_initial_conditions_papa_meteo => set_initial_conditions
+        use obl_run_cbl, only : set_initial_conditions_cbl => set_initial_conditions
 
         type (gridDataType), intent(in) :: grid
 
@@ -421,6 +448,10 @@ contains
             call set_initial_conditions_papa_meteo(grid)
             return
         endif
+        if (config_id == obl_config_cbl) then
+            call set_initial_conditions_cbl(grid)
+            return
+        endif
 
         block 
             !< *: will fail without configuration file
diff --git a/obl_diag.f90 b/obl_diag.f90
index 62dc45d89eef40239409984337210e6f93894da4..d25ff0ae4035d709ccccff7e260e98861f49c259 100644
--- a/obl_diag.f90
+++ b/obl_diag.f90
@@ -98,7 +98,7 @@ module obl_diag
         k_min = cz
 
         do k = cz - 1, 1, -1
-            if (B(k) < B_min) then
+            if (B(k) <= B_min) then
                 B_min = B(k)
                 k_min = k
             end if
diff --git a/obl_main.f90 b/obl_main.f90
index dbc36936264e2d23e5ee19fe7e9c95aed2cb447e..9931cc081cd1931eadf245b9a113094a352e42dc 100644
--- a/obl_main.f90
+++ b/obl_main.f90
@@ -55,6 +55,7 @@ program obl_main
                                 ! = obl_config_kato: Kato-Phillips
                                 ! = obl_config_papa_fluxes: Papa-station (fluxes)
                                 ! = obl_config_papa_meteo: Papa-station (meteo)
+                                ! = obl_config_cbl: convective boundary layer (Willis exp.)
 
     integer :: obl_model_id     ! --- OBL model def.
                                 ! = obl_model_pph: pacanowski-philander
@@ -94,6 +95,7 @@ program obl_main
     real, parameter :: Cd0 = 0.001      ! bottom drag coefficient
 
     real :: mld, mld_ref                    ! mixed layer depth (model & reference), [m]
+    real :: eld, eld_ref                    ! entrainment layer depth (model & reference), [m]
     real, parameter :: N2_ref = 0.000044    ! reference N**2 (using in mld ref. estimate), [1/s**2]
     
     ! command line arguments & configuration file variables
@@ -115,6 +117,7 @@ program obl_main
     !<      = obl_config_kato
     !<      = obl_config_papa_fluxes
     !<      = obl_config_papa_meteo
+    !<      = obl_config_cbl
     obl_config_id = obl_config_kato
 
     !< default closure = 4, k-epsilon
@@ -141,7 +144,7 @@ program obl_main
             write(*, *) ' --help'
             write(*, *) '    print usage options'
             write(*, *) ' --config [key] || [filename]'
-            write(*, *) '    builtin setup [key] = kato (default) || papa-fluxes || papa'
+            write(*, *) '    builtin setup [key] = kato (default) || papa-fluxes || papa || cbl'
             write(*, *) '    use configuration file [filename]'
             write(*, *) ' --model [key]'
             write(*, *) '    [key] = pph || pph-dyn || k-epsilon (default)'
@@ -335,10 +338,11 @@ program obl_main
         ! ----------------------------------------------------------------------------
         if (mod(time_index, nscreen) == 0) then
             call get_mld(mld, N2, grid%dz, grid%cz)
+            call get_eld(eld, TKE_buoyancy, grid%dz, grid%cz)
 
             ! *: add finite check
 
-            write(*, '(a,g0)') ' mld = ', mld
+            write(*, '(a,g0,a,g0)') ' mld = ', mld, ' eld = ', eld
             write(*, '(a,g0)') ' Theta(surface) = ', Theta_dev(grid%cz) + Theta_ref
             write(*, '(a,g0,a,g0,a)') ' current time = ', time_current / 3600.0, ' HRS [ ', & 
                 (time_current / time_end) * 100.0, '% ]' 
diff --git a/obl_run_cbl.f90 b/obl_run_cbl.f90
new file mode 100644
index 0000000000000000000000000000000000000000..74e41087f40f4c59f6b0e68dfb3210c25f1d62b0
--- /dev/null
+++ b/obl_run_cbl.f90
@@ -0,0 +1,127 @@
+module obl_run_cbl
+    !< @brief obl scm convective boundary layer (Willis experiment) setup
+    ! --------------------------------------------------------------------------------
+
+    ! TO DO:
+    !   -- ***
+
+    ! modules used
+    ! --------------------------------------------------------------------------------
+
+    ! directives list
+    ! --------------------------------------------------------------------------------
+    implicit none
+    private
+
+    ! public interface
+    ! --------------------------------------------------------------------------------
+    public :: set_phys
+    public :: set_grid
+    public :: set_time
+    public :: set_forcing
+    public :: set_initial_conditions
+    ! --------------------------------------------------------------------------------
+
+    ! --------------------------------------------------------------------------------
+    ! --------------------------------------------------------------------------------
+    
+
+    contains
+
+    ! --------------------------------------------------------------------------------
+    subroutine set_phys
+        !> @brief phys parameters setup
+        ! ----------------------------------------------------------------------------
+        use obl_scm
+        ! ----------------------------------------------------------------------------
+
+        !< coriolis frequency
+        f = 0.0
+
+    end subroutine set_phys
+
+    ! --------------------------------------------------------------------------------
+    subroutine set_grid(grid)
+        !> @brief grid parameters setup
+        ! ----------------------------------------------------------------------------
+        use obl_grid
+
+        type (gridDataType), intent(inout) :: grid
+        ! ----------------------------------------------------------------------------
+
+        !< in: [zpos, height, cz]
+        call set_uniform_grid(grid, - 32.0, 32.0, 32)
+
+    end subroutine set_grid
+
+    ! --------------------------------------------------------------------------------
+    subroutine set_time(time_begin, time_end, dt)
+        !> @brief time parameters setup
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: time_begin, time_end, dt
+        ! ----------------------------------------------------------------------------
+
+        time_begin = 0.0
+        time_end = 72.0 * 3600.0
+        dt = 1.0
+
+    end subroutine set_time
+
+    ! --------------------------------------------------------------------------------
+    subroutine set_forcing
+        !> @brief forcing setup
+        ! ----------------------------------------------------------------------------
+        use obl_fluxes
+        use obl_tforcing
+        ! ----------------------------------------------------------------------------
+
+        !< setting atmospheric forcing
+        ! ---------------------------------------------------------------------------- 
+        !< using 'flux' mode
+        is_meteo_setup = 0
+
+        call set_const_tforcing(sensible_hflux_surf, -100.0)
+        call set_const_tforcing(latent_hflux_surf, 0.0)
+
+        call set_const_tforcing(salin_flux_surf, 0.0)
+
+        call set_const_tforcing(tau_x_surf, 0.0)
+        call set_const_tforcing(tau_y_surf, 0.0)
+
+        call set_const_tforcing(sw_flux_surf, 0.0)
+
+        call set_const_tforcing(lw_net_surf, 0.0)
+        ! ---------------------------------------------------------------------------- 
+
+        !< setting bottom forcing
+        ! ----------------------------------------------------------------------------
+        call set_const_tforcing(hflux_bot, 0.0)
+
+        call set_const_tforcing(salin_flux_bot, 0.0)
+
+        call set_const_tforcing(tau_x_bot, 0.0)
+        call set_const_tforcing(tau_y_bot, 0.0)
+        ! ----------------------------------------------------------------------------
+
+    end subroutine set_forcing
+
+    ! --------------------------------------------------------------------------------
+    subroutine set_initial_conditions(grid)
+        !> @brief initial_conditions setup
+        ! ----------------------------------------------------------------------------
+        use obl_state
+        use obl_init
+        use obl_grid
+
+        type (gridDataType), intent(in) :: grid
+        ! ----------------------------------------------------------------------------
+
+        call set_linear_profile(Theta, 295.15, 0.1, grid)
+        call set_const_profile(Salin, 35.0, grid)
+
+        call set_const_profile(U, 0.0, grid)
+        call set_const_profile(V, 0.0, grid)
+
+    end subroutine set_initial_conditions
+
+end module