diff --git a/io_metadata.f90 b/io_metadata.f90
index 994dfe42fa41d47848afa359793e0d5a21245a94..3b6e37dd5f5b53c11fb9e017ddc86e241bef3b16 100644
--- a/io_metadata.f90
+++ b/io_metadata.f90
@@ -339,6 +339,24 @@ module io_metadata
     real_name = [character(len=32) :: 'missing_value', '', '', '', '', '', '', '', '', ''], &
     real_value = [ -9999.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ])
 
+    ! Metadata for eld
+    type(variable_metadata), public :: meta_eld = variable_metadata( &
+    name = 'eld', &
+    dim_names = [character(len=32) :: 'Time', '', '', ''], &
+    char_name = [character(len=32) :: 'long_name', 'units', '', '', '', '', '', '', '', ''], &
+    char_value = [character(len=32) :: 'Entrainment Layer Depth', 'm', '', '', '', '', '', '', '', ''], &
+    real_name = [character(len=32) :: 'missing_value', '', '', '', '', '', '', '', '', ''], &
+    real_value = [ -9999.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ])
+     
+    ! Metadata for lab_eld
+    type(variable_metadata), public :: meta_lab_eld = variable_metadata( &
+    name = 'lab_eld', &
+    dim_names = [character(len=32) :: 'Time', '', '', ''], &
+    char_name = [character(len=32) :: 'long_name', 'units', '', '', '', '', '', '', '', ''], &
+    char_value = [character(len=32) :: 'Entrainment Layer Depth', 'm', '', '', '', '', '', '', '', ''], &
+    real_name = [character(len=32) :: 'missing_value', '', '', '', '', '', '', '', '', ''], &
+    real_value = [ -9999.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ])
+
     ! Metadata for tau_u
     type(variable_metadata), public :: meta_tau_u = variable_metadata( &
     name = 'tau_u', &
diff --git a/obl_main.f90 b/obl_main.f90
index 9931cc081cd1931eadf245b9a113094a352e42dc..65e7bdfb42fc577e694931682a6d1c74b3bf5663 100644
--- a/obl_main.f90
+++ b/obl_main.f90
@@ -96,7 +96,11 @@ program obl_main
 
     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]
+    real, parameter :: N2_ref = 0.000044    ! reference N**2 (used in mld ref. estimate), [1/s**2]
+                                            ! = 0.000044 (default Kato)
+                                            ! = 2.0 * 1e-4 (default CBL)
+    real, parameter :: Bsurf_ref = 0.52 * 1e-7      ! reference buoyancy flux (used in eld ref. estimate)
+                                                    ! = 0.52 * 1e-7 (default CBL)
     
     ! command line arguments & configuration file variables
     ! --------------------------------------------------------------------------------
@@ -358,9 +362,15 @@ program obl_main
             call get_mld(mld, N2, grid%dz, grid%cz)
             call get_mld_ref(mld_ref, bc%U_dynH, N2_ref, time_current, grid%height)
 
+            call get_eld(eld, TKE_buoyancy, grid%dz, grid%cz)
+            call get_eld_ref(eld_ref, N2_ref, Bsurf_ref, time_current, grid%height)
+
             call push_value_to_tseries(scm_output%mld, mld)
             call push_value_to_tseries(scm_output%mld_ref, mld_ref)
 
+            call push_value_to_tseries(scm_output%eld, eld)
+            call push_value_to_tseries(scm_output%eld_ref, eld_ref)
+
             call push_value_to_tseries(scm_output%tau_x, bc%u_momentum_fluxH * Rho_ref)
             call push_value_to_tseries(scm_output%tau_y, bc%v_momentum_fluxH * Rho_ref)
 
diff --git a/obl_output.f90 b/obl_output.f90
index 5eb2e24e09fda6eae634f54dfeffa1cd808314d6..11cdbbf0db5e86220f8bff9e995bad36e7992aec 100644
--- a/obl_output.f90
+++ b/obl_output.f90
@@ -43,6 +43,9 @@ module obl_output
         type(oblTimeSeries), public :: mld
         type(oblTimeSeries), public :: mld_ref
 
+        type(oblTimeSeries), public :: eld
+        type(oblTimeSeries), public :: eld_ref
+
         type(oblTImeSeries), public :: tau_x, tau_y
 
         type(oblTimeSeries), public :: Theta_surf
@@ -123,6 +126,9 @@ module obl_output
         call write_1d_real_nc(output%mld%data(1:num), trim(path) // 'output.nc', meta_mld)
         call write_1d_real_nc(output%mld_ref%data(1:num), trim(path) // 'output.nc', meta_lab_mld)
 
+        call write_1d_real_nc(output%eld%data(1:num), trim(path) // 'output.nc', meta_eld)
+        call write_1d_real_nc(output%eld_ref%data(1:num), trim(path) // 'output.nc', meta_lab_eld)
+
         call write_1d_real_nc(output%tau_x%data(1:num), trim(path) // 'output.nc', meta_tau_u)
         call write_1d_real_nc(output%tau_y%data(1:num), trim(path) // 'output.nc', meta_tau_v)
 
@@ -184,6 +190,11 @@ module obl_output
         call write_1d_real_ascii(output%mld_ref%data(1:num), &
             trim(path) // trim(meta_lab_mld%name) // '.txt', meta_lab_mld)
 
+        call write_1d_real_ascii(output%eld%data(1:num), &
+            trim(path) // trim(meta_eld%name) // '.txt', meta_eld)
+        call write_1d_real_ascii(output%eld_ref%data(1:num), &
+            trim(path) // trim(meta_lab_eld%name) // '.txt', meta_lab_eld)
+
         call write_1d_real_ascii(output%tau_x%data(1:num), &
             trim(path) // trim(meta_tau_u%name) // '.txt', meta_tau_u)
         call write_1d_real_ascii(output%tau_y%data(1:num), &
@@ -249,6 +260,11 @@ module obl_output
         call write_1d_real_plt(output%mld_ref%data(1:num), output%time%data(1:num), &
             trim(path) // trim(meta_lab_mld%name) // '.plt', meta_lab_mld)
 
+        call write_1d_real_plt(output%eld%data(1:num), output%time%data(1:num), &
+            trim(path) // trim(meta_eld%name) // '.plt', meta_eld)
+        call write_1d_real_plt(output%eld_ref%data(1:num), output%time%data(1:num), &
+            trim(path) // trim(meta_lab_eld%name) // '.plt', meta_lab_eld)
+
         call write_1d_real_plt(output%tau_x%data(1:num), output%time%data(1:num), &
             trim(path) // trim(meta_tau_u%name) // '.plt', meta_tau_u)
         call write_1d_real_plt(output%tau_y%data(1:num), output%time%data(1:num), &
@@ -290,6 +306,9 @@ module obl_output
         call deallocate_tseries(output%mld)
         call deallocate_tseries(output%mld_ref)
 
+        call deallocate_tseries(output%eld)
+        call deallocate_tseries(output%eld_ref)
+
         call deallocate_tseries(output%tau_x)
         call deallocate_tseries(output%tau_y)
 
diff --git a/obl_run_cbl.f90 b/obl_run_cbl.f90
index 74e41087f40f4c59f6b0e68dfb3210c25f1d62b0..3d44e7c786a18d0e8bf8c9ee4d0fe5ea93e61c3b 100644
--- a/obl_run_cbl.f90
+++ b/obl_run_cbl.f90
@@ -80,7 +80,7 @@ module obl_run_cbl
         !< using 'flux' mode
         is_meteo_setup = 0
 
-        call set_const_tforcing(sensible_hflux_surf, -100.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)