diff --git a/CMakeLists.txt b/CMakeLists.txt
index a6c7d111531d8ed72ac61a73a9da762846a22b88..46f863338fa9d31da264fa93b2f1089096fc478a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -35,7 +35,10 @@ set(lib_files
     src/state_utils.f90
     src/scm_state_data.f90
     src/pbl_turb_data.f90
-        src/pbl_solver.f90
+    src/pbl_solver.f90
+    src/fo_stability_functions.f90
+    src/pbl_fo_turb.f90
+    src/pbl_turb_common.f90
     src/diag_pbl.f90)
 
 add_library(abl_lib ${lib_files})
@@ -117,4 +120,6 @@ if (WITH_TESTS)
                     PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -fbacktrace -ffpe-trap=zero,overflow,underflow >)
         endif()
     endif()
-endif()
\ No newline at end of file
+
+endif()
+
diff --git a/config-examples/config-gabls1-ex.txt b/config-examples/config-gabls1-ex.txt
index 5f885834cbdf2e3f3d943aca4226f69cd21ee8a6..5cf74b0bbfd195cd76d8b6f1e4ecd18379b80efd 100644
--- a/config-examples/config-gabls1-ex.txt
+++ b/config-examples/config-gabls1-ex.txt
@@ -1,12 +1,16 @@
 
 Ug = 8.0;
 Vg = 0.0;
+hbl_0 = 100.0;
+Tgrad = 0.01;
+cooling_rate = -0.25 * 3600.0;
+
 time{
     begin = 0.0;
-    end = 9.0 * 3600.0;
+    end = 20.0; #9.0 * 3600.0;
     dt = 20.0;
     out_dt = 360.0;
-    }
+}
 
 fluid {
     pref= 1013.4; #hPa
@@ -14,16 +18,28 @@ fluid {
     g = 9.81; # m/s2
     beta = g /tref;
     kappa = 0.4; # von Karman constant
-    }
+}
 
-surface_model {
+grid {
+    type = "uniform";
 
-        id = "esm";       # optional: "esm" (default), "sheba", "most", "log"
-        type = "land";
-        z0_m = 0.01;                      # aerodynamic roughness [m]
-        z0_h = -1;                        # no prescribed value                                                #       -> using scheme assigned by surface type
-    }
+    nz = 8;
+    h_bot = 0.0;
+    h_top = 400.0;
+}
+
+surface_model {
+    id = "esm";       # optional: "esm" (default), "sheba", "most", "log"
+    type = "land";
+    z0_m = 0.01;                      # aerodynamic roughness [m]
+    z0_h = -1;                        # no prescribed value                                                #       -> using scheme assigned by surface type
+}
 
 output{
     dt = 0.2 *3600.0;
+}
+
+closure {
+      type = "FOM";
+      name = "Louis"; # "Esau"; # "EFB";
 }
\ No newline at end of file
diff --git a/src/config-utils.f90 b/src/config-utils.f90
index d3026263401b6b2e6d450cb8539a0352b10ee787..0076d75c01d458bf38874c933f2f85356c733323 100644
--- a/src/config-utils.f90
+++ b/src/config-utils.f90
@@ -8,6 +8,7 @@ module config_utils
     implicit none
     integer, public, save:: is_config_initialized = 0
     public init_config, get_fluid_params, get_grid_params
+    public get_closure_params
     !public get_geo_forcing, get_heat_forcing
 
 
@@ -86,6 +87,7 @@ module config_utils
                 grid_inmcm21_tag, grid_inmcm73_tag, &
             grid_streached_tag, grid_uniform_tag, &
             set_pbl_grid_uniform
+    use sfx_common, only: char_array2str
 
     implicit none
 
@@ -105,15 +107,18 @@ module config_utils
         if ((status /= 0)) then
             call c_config_get_string("grid.type"//C_NULL_CHAR, config_field, status)
             if (status == 0) then
+                write(*, *) ' FAILURE! > could not set:  grid.type '
                 ierr = 1        ! signal ERROR
                 return
             end if
         end if
+        tag = char_array2str(config_field)
         if (trim(tag) == trim(grid_uniform_tag)) then
             call c_config_is_varname("grid.nz"//C_NULL_CHAR, status)
             if ((status /= 0)) then
                 call c_config_get_int("grid.nz"//C_NULL_CHAR, nz, status)
                 if (status == 0) then
+                    write(*, *) ' FAILURE! > could not set:  grid.nz '
                     ierr = 1        ! signal ERROR
                     return
                 end if
@@ -123,6 +128,7 @@ module config_utils
             if ((status /= 0)) then
                 call c_config_get_float("grid.h_bot"//C_NULL_CHAR, h_bot, status)
                 if (status == 0) then
+                    write(*, *) ' FAILURE! > could not set:  grid.h_bot '
                     ierr = 1        ! signal ERROR
                     return
                 end if
@@ -131,6 +137,7 @@ module config_utils
             if ((status /= 0)) then
                 call c_config_get_float("grid.h_top"//C_NULL_CHAR, h_top, status)
                 if (status == 0) then
+                    write(*, *) ' FAILURE! > could not set:  grid.h_top '
                     ierr = 1        ! signal ERROR
                     return
                 end if
@@ -144,6 +151,63 @@ module config_utils
     end subroutine get_grid_params
 
 
+        subroutine get_closure_params(turb, status, ierr)
+            use pbl_turb_data, only: turbBLDataType, &
+                    pbl_fom_tag, pbl_kl_tag, &
+                    pbl_fom_louis_tag, pbl_fom_esau_tag, &
+                    pbl_fom_efb_tag
+            use sfx_common, only: char_array2str
+
+            implicit none
+            type(turbBLDataType), intent(inout):: turb
+            integer,intent(out):: status, ierr
+            character(len=50) :: tag, tag_stabf
+            character, allocatable :: config_field(:)
+            ierr = 0
+            if( is_config_initialized /= 0 ) then
+                ! closure type
+                call c_config_is_varname("closure.type"//C_NULL_CHAR, status)
+                if ((status /= 0)) then
+                    call c_config_get_string("closure.type"//C_NULL_CHAR, config_field, status)
+                    if (status == 0) then
+                        write(*, *) ' FAILURE! > could not set:  grid.type '
+                        ierr = 1        ! signal ERROR
+                        return
+                    end if
+                end if
+                tag = char_array2str(config_field)
+                if (trim(tag) == trim(pbl_fom_tag)) then
+
+                    turb%cl_type = 1
+                    call c_config_is_varname("closure.name"//C_NULL_CHAR, status)
+                    if ((status /= 0)) then
+                        call c_config_get_string("closure.name"//C_NULL_CHAR, config_field, status)
+                        if (status == 0) then
+                            write(*, *) ' FAILURE! > could not set:  closure.name '
+                            ierr = 1        ! signal ERROR
+                            return
+                        else
+                            tag_stabf = char_array2str(config_field)
+                            if (trim(tag_stabf) == trim(pbl_fom_louis_tag)) then
+                                turb%sf_type = 1
+                            elseif (trim(tag_stabf) == trim(pbl_fom_esau_tag)) then
+                                turb%sf_type = 2
+                            elseif (trim(tag_stabf) == trim(pbl_fom_efb_tag)) then
+                                turb%sf_type = 3
+                            else
+                                status = 1
+                                ierr  = 1
+                                write(*, *) ' FAILURE! > could not find eligable:  closure.name '
+                                write(*,*) tag_stabf, pbl_fom_louis_tag
+                                return
+                            end if
+                        end if
+                    end if
+                else
+                    turb%cl_type = 2 ! KL - type
+                endif
+            end if
+        end subroutine get_closure_params
 
     !> @brief character array to string conversion
     function char_array2str(char_array) result(str)
diff --git a/src/fo_stability_functions.f90 b/src/fo_stability_functions.f90
new file mode 100644
index 0000000000000000000000000000000000000000..0f321774cc7931e5bcaccec789569d182858c81c
--- /dev/null
+++ b/src/fo_stability_functions.f90
@@ -0,0 +1,91 @@
+! Created by Andrey Debolskiy on 24.01.2025.
+module fo_stability_functions
+    !< Constants to for stability functions
+    !< louis stability functions
+    real, parameter:: yb = 5.e0
+    real, parameter:: yc = 5.e0
+    real, parameter:: yd = 5.e0
+    real, parameter:: ye = yb*yc*sqrt(1.e0/3.e0)
+    real, parameter::ricr = 2.e0/(3.e0*yd)
+    real, parameter:: vlinfm = 120.0e0
+    real, parameter:: vlinfh = vlinfm * sqrt(3.e0*yd/2.e0)
+
+
+    !    security parameters
+    !    du2min is a minimal value to calculate wind shear
+    !   cormin is a minimal value of coriolis parameter to calculate the
+    !   pbl extention
+    real, parameter::  du2min = 1.e-02
+    real, parameter::  cormin = 5.e-05
+    contains
+
+    !>
+    subroutine louis_stab_functions(fnriuv,fnritq,rinum,vdlm, &
+            & vdlh,zkh, rolim)
+        implicit none
+        real fnriuv,fnritq
+        real, intent(in):: rinum
+        real, intent(in):: vdlm,vdlh, zkh
+        real scf, ucfm, ucfh
+        !real yb, yc, yd, ye, ricr
+        real, intent(in):: rolim
+
+        if(rinum.ge.0.e0) then
+            scf = sqrt(1.e0 + yd*rinum)
+            fnriuv = 1.e0/(1.e0 + 2.e0*yb * rinum / scf)
+            fnritq = 1.e0/(1.e0 + 3.e0*yb * rinum * scf)
+        else
+            ucfm = 1.e0/(1.e0 + ye * (vdlm/zkh)**2 * sqrt(abs(rinum)))
+            ucfh = 1.e0/(1.e0 + ye * (vdlh/zkh)**2 * sqrt(abs(rinum)))
+            fnriuv = 1.e0 - 2.e0*yb * rinum * ucfm
+            if(rolim.gt.0.99.and.rolim.lt.1.5) then
+                fnritq = 1.e0 - 3.0*3.e0*yb * rinum * ucfh    !enhanced mixing for unstable stratification over land
+            else
+                fnritq = 1.e0 - 3.e0*yb * rinum * ucfh
+            end if
+        end if
+    end subroutine louis_stab_functions
+
+    !>
+    !>esau byrkjedal 2007  +   viterbo 1999  unstable stratification
+    subroutine esau_vit_stab_functions( fnriuv,fnritq,rinum, &
+            &   vdlm,zkh, dz)
+        implicit none
+        real fnriuv,fnritq,rinum, vdlm, zkh, dz
+        real scf, ucfm, ucfh
+        real yb, yc, yd, ye, ricr
+        real gh
+        real, parameter:: p_m = 21.0
+        real, parameter:: q_m = 0.005
+        real, parameter:: p_h = 10.0
+        real, parameter:: q_h = 0.0012
+        real, parameter:: bb = 5.0
+        real, parameter:: cc = 5.0
+
+        yb = 5.e0
+        yc = 5.e0
+        yd = 5.e0
+        ye = yb*yc*sqrt(1.e0/3.e0)
+        ricr = 2.e0/(3.e0*yd)
+
+
+        if(rinum.ge.0.e0) then
+            scf = sqrt(1.e0 + yd*rinum)
+            fnriuv = 1.0/((1.0 + p_m * rinum)*(1.0 + p_m * rinum)) &
+                    &   + q_m * sqrt(rinum)
+            fnritq = 1.0/((1.0+p_h*rinum)*(1.0+p_h*rinum)*(1.0+p_h*rinum)) &
+                    &  + q_h
+        else
+            gh = ((vdlm * vdlm)/(dz * sqrt(abs(dz * zkh)))) &
+                    &     * ((1 + dz / zkh)**(1/3) - 1)**(3/2)
+
+            fnriuv =  1.0 - (2.0 * bb * rinum) &
+                    &  / (1.0 + 3.0 * bb * cc * gh * sqrt(abs(rinum)))
+
+            fnritq = 1.0 - (2.0 * bb * rinum) &
+                    &  / (1.0 + 3.0 * bb * cc * gh * sqrt(abs(rinum)))
+        end if
+    end subroutine esau_vit_stab_functions
+
+
+end module fo_stability_functions
\ No newline at end of file
diff --git a/src/pbl_fo_turb.f90 b/src/pbl_fo_turb.f90
index 959dbccb60bed4a3e57a342506f1b31bed6a659c..cf617e14f5d90bf4edf58bc11c2966b97631eb01 100644
--- a/src/pbl_fo_turb.f90
+++ b/src/pbl_fo_turb.f90
@@ -2,25 +2,93 @@
 
 module pbl_fo_turb
     implicit none
-
+    !constants to for stability functions
+    real, parameter:: yb = 5.e0
+    real, parameter:: yc = 5.e0
+    real, parameter:: yd = 5.e0
+    real:: ye = yb*yc*sqrt(1.e0/3.e0)
+    real:: ricr = 2.e0/(3.e0*yd)
+    real, parameter:: vlinfm = 160.0e0
+    real:: vlinfh = vlinfm * sqrt(3.e0*yd/2.e0)
+    real, parameter::  aka = 0.4e0
+    !real, parameter:: r = 287.05
+    !real, parameter:: cp = 1.005
+    !    security parameters
+    !    du2min is a minimal value to calculate wind shear
+    !   cormin is a minimal value of coriolis parameter to calculate the
+    !   pbl extention
+    real, parameter::  du2min = 1.e-02
+    real, parameter::  cormin = 5.e-05
     public get_fo_coeffs
+    public get_fo_coeffs_louis
     public get_turb_length
-    public get_fo_diff
-    public get_rig
+!    public :: default_stab_functions
+!    public :: esau_stab_functions
+!    private:: get_ph_efb
+!    private:: q_root
     contains
-    subroutine get_fo_diff(turb, bl, fluid, grid)
+
+!    subroutine get_fo_diff(turb, bl, fluid, grid)
+!        use pbl_turb_data, only: turbBLDataType
+!        use scm_state_data, only: stateBLDataType
+!        use pbl_grid, only: pblgridDataType
+!        use phys_fluid, only: fluidParamsDataType
+!    end subroutine get_fo_diff
+    subroutine get_fo_coeffs(turb, bl, fluid,  grid)
         use pbl_turb_data, only: turbBLDataType
         use scm_state_data, only: stateBLDataType
         use pbl_grid, only: pblgridDataType
         use phys_fluid, only: fluidParamsDataType
-    end subroutine get_fo_diff
+        implicit none
+        type(fluidParamsDataType), intent(in):: fluid
+        type(stateBLDataType), intent(inout):: bl
+        type(pblgridDataType), intent(in)::grid
+        type(turbBLDataType), intent(inout):: turb
 
-    subroutine get_fo_coeffs(turb, bl, fluid,  grid)
+        if (turb%sf_type == 1) then !> Louis stab functions
+            call get_fo_coeffs_louis(turb, bl, fluid,  grid)
+        elseif (turb%sf_type == 2) then !> Esau stab functions
+            ! call get_fo_coeffs_esau(turb, bl, fluid,  grid)
+        elseif  (turb%sf_type == 3)    then   !> EFB stab functions
+            ! call get_fo_coeffs_efb(turb, bl, fluid,  grid)
+        else
+            write(*,*) 'FAILURE: stab functions type in FOM not set!'
+        end if
+    end subroutine get_fo_coeffs
+
+    subroutine get_fo_coeffs_louis(turb, bl, fluid,  grid)
         use pbl_turb_data, only: turbBLDataType
         use scm_state_data, only: stateBLDataType
         use pbl_grid, only: pblgridDataType
         use phys_fluid, only: fluidParamsDataType
-    end subroutine get_fo_coeffs
+        use fo_stability_functions
+        implicit none
+        type(fluidParamsDataType), intent(in):: fluid
+        type(stateBLDataType), intent(inout):: bl
+        type(pblgridDataType), intent(in)::grid
+        type(turbBLDataType), intent(inout):: turb
+
+        real fnriuv, fnritq,rho
+        real ocean_flag
+
+        integer k, kmax
+
+        ocean_flag = real(bl%land_type)
+        do k = kmax, 1, -1
+            call louis_stab_functions(fnriuv, fnritq, turb%rig(k), turb%l_turb_m(k), &
+                    & turb%l_turb_m(k), grid%z_cell(k),ocean_flag)
+            turb%km(k) = fnriuv * turb%l_turb_m(k) * turb%l_turb_m(k) * turb%s2(k)
+            turb%kh(k) = fnritq * turb%l_turb_h(k) * turb%l_turb_h(k) * turb%s2(k)
+        end do
+        ! for compliance with old version
+        bl%vdcuv(grid%kmax) = 0.0
+        bl%vdcuv(grid%kmax) = 0.0
+        do k = kmax-1,1, -1
+            bl%vdcuv(k) = 0.5 * (bl%rho(k) + bl%rho(k+1)) * turb%km(k)
+            bl%vdctq(k) = 0.5 * (bl%rho(k) + bl%rho(k+1)) * turb%kh(k)
+        end do
+
+    end subroutine get_fo_coeffs_louis
 
     subroutine get_turb_length(turb, bl, fluid, grid, hbl_option)
         use pbl_turb_data, only: turbBLDataType
@@ -30,19 +98,40 @@ module pbl_fo_turb
         use diag_pbl
         implicit none
         type(fluidParamsDataType), intent(in):: fluid
-        type(stateBLDataType), intent(in):: bl
+        type(stateBLDataType), intent(inout):: bl
         type(pblgridDataType), intent(in)::grid
-        type(turbBLDataTyoe), intent(inout):: turb
+        type(turbBLDataType), intent(inout):: turb
+        integer hbl_option
+        real l_neutral, lm, lh
 
-        call get_rig(turb, bl, fluid, grid)
+        integer k, kmax
 
+        if (hbl_option ==1) then
+            call get_hpbl(bl%hpbla_diag, bl%kpbl, bl%theta_v, grid%z_cell, &
+                    grid%z_edge(kmax), grid%kmax , bl%cor, turb%ustr)
+        else
+            write(*,*) 'wrong hbl diagnostics option'
+            return
+        end if
+
+        do k = kmax,1, -1
+            l_neutral = fluid%kappa * grid%z_cell(k)
+            if ( k <= bl%kpbl) then
+                lm = 0.0
+                lh = 0.0
+            else
+                if (turb%rig(k) >= 0.0) then
+                    lm = (1.0 - grid%z_cell(k) / bl%hpbla_diag)
+                    lh = (1.0 - grid%z_cell(k) / bl%hpbla_diag)
+                else
+                    lm = 1.0
+                    lh = 1.0
+                endif
+            endif
+            turb%l_turb_m = l_neutral * lm
+            turb%l_turb_h = l_neutral * lh
+        end do
 
     end subroutine get_turb_length
 
-    subroutine get_rig(turb, bl, fluid, grid)
-        use pbl_turb_data, only: turbBLDataType
-        use scm_state_data, only: stateBLDataType
-        use pbl_grid, only: pblgridDataType
-        use phys_fluid, only: fluidParamsDataType
-    end subroutine get_rig
 end module pbl_fo_turb
\ No newline at end of file
diff --git a/src/pbl_grid.f90 b/src/pbl_grid.f90
index 90a0d0a992106fec0a76e8b05e1efad16cc7c153..cb1fe277857ada38feb288e995357f31dc61d1f1 100644
--- a/src/pbl_grid.f90
+++ b/src/pbl_grid.f90
@@ -95,6 +95,7 @@ module pbl_grid
             call allocate_pbl_grid(grid, nk)
             dz = (h_top - h_bottom) / real(nk, kind = kind(dz))
             grid%z_edge(nk) = h_bottom
+            grid%z_cell(nk) = grid%z_edge(nk) + 0.5 * dz
             do k = nk-1, 1, -1
                 grid%z_edge(k) = grid%z_edge(k+1) + dz
                 grid%z_cell(k) = 0.5 * (grid%z_edge(k+1) + grid%z_edge(k))
diff --git a/src/pbl_turb_common.f90 b/src/pbl_turb_common.f90
new file mode 100644
index 0000000000000000000000000000000000000000..948abbccd7ee0808301df7cb2e2c85b1cfc51e26
--- /dev/null
+++ b/src/pbl_turb_common.f90
@@ -0,0 +1,112 @@
+! Created by Andrey Debolskiy on 11.02.2025.
+
+!> Collection of utilities for pbl_turb_data
+module pbl_turb_common
+use pbl_turb_data
+implicit none
+    public get_s2, get_n2, get_rig
+    public get_coeffs_general
+
+contains
+    subroutine get_n2(turb, bl, fluid, grid)
+        use pbl_turb_data, only: turbBLDataType
+        use scm_state_data, only: stateBLDataType
+        use pbl_grid, only: pblgridDataType
+        use phys_fluid, only: fluidParamsDataType
+        implicit none
+        type(fluidParamsDataType), intent(in):: fluid
+        type(stateBLDataType), intent(in):: bl
+        type(pblgridDataType), intent(in)::grid
+        type(turbBLDataType), intent(inout):: turb
+
+        real dtvir, buoyp
+        integer k
+        do  k = 1,grid%kmax-1
+            dtvir = bl%theta_v(k) - bl%theta_v(k+1)
+            buoyp = fluid%g / ( (bl%theta_v(k+1) + bl%theta_v(k))/2.0e0 )
+            turb%n2 = buoyp*(grid%z_cell(k) - grid%z_cell(k+1))*dtvir
+        end do
+    end subroutine get_n2
+
+    subroutine get_s2(turb, bl, fluid, grid, du2min)
+        use pbl_turb_data, only: turbBLDataType
+        use scm_state_data, only: stateBLDataType
+        use pbl_grid, only: pblgridDataType
+        use phys_fluid, only: fluidParamsDataType
+        implicit none
+        type(fluidParamsDataType), intent(in):: fluid
+        type(stateBLDataType), intent(in):: bl
+        type(pblgridDataType), intent(in)::grid
+        type(turbBLDataType), intent(inout):: turb
+        real, intent(in) :: du2min
+
+        real du2
+        integer k
+        do  k = 1,grid%kmax-1
+            write(*,*) 'k=', k
+            du2 = amax1(du2min, &
+                    &    (bl%u(k)-bl%u(k+1))**2 + (bl%v(k)-bl%v(k+1))**2)
+            turb%s2 = du2 / ((grid%z_cell(k) - grid%z_cell(k+1)) &
+                        * (grid%z_cell(k) - grid%z_cell(k+1)))
+        enddo
+    end subroutine get_s2
+
+    subroutine get_rig(turb, bl, fluid, grid, du2min, ricr)
+    use pbl_turb_data, only: turbBLDataType
+    use scm_state_data, only: stateBLDataType
+    use pbl_grid, only: pblgridDataType
+    use phys_fluid, only: fluidParamsDataType
+    implicit none
+    type(fluidParamsDataType), intent(in):: fluid
+    type(stateBLDataType), intent(in):: bl
+    type(pblgridDataType), intent(in)::grid
+    real, intent(in):: du2min, ricr
+    type(turbBLDataType), intent(inout):: turb
+
+    real dtvir, buoyp, du2, rinum
+    integer k
+    do  k = 1,grid%kmax-1
+        write(*,*)'here rig'
+        du2 = amax1(du2min, &
+                &    (bl%u(k)-bl%u(k+1))**2 + (bl%v(k)-bl%v(k+1))**2)
+        dtvir = bl%theta_v(k) - bl%theta_v(k+1)
+        buoyp = fluid%g / ( (bl%theta_v(k+1) + bl%theta_v(k))/2.0e0 )
+        rinum = buoyp*(grid%z_cell(k) - grid%z_cell(k+1))*dtvir/du2
+        rinum = amin1(ricr,rinum)
+        turb%rig(k) = rinum
+    end do
+    end subroutine get_rig
+
+    subroutine get_coeffs_general(turb, bl, fluid, hbl_option,  grid)
+        use pbl_turb_data, only: turbBLDataType
+        use scm_state_data, only: stateBLDataType
+        use pbl_grid, only: pblgridDataType
+        use phys_fluid, only: fluidParamsDataType
+        use pbl_fo_turb, only: get_fo_coeffs, get_turb_length, du2min, ricr
+        implicit none
+        type(fluidParamsDataType), intent(in):: fluid
+        type(stateBLDataType), intent(inout):: bl
+        type(pblgridDataType), intent(in)::grid
+        type(turbBLDataType), intent(inout):: turb
+        integer, intent(in):: hbl_option
+
+
+        write(*,*) 'here1'
+        call  get_n2(turb, bl, fluid, grid)
+        write(*,*) 'here2'
+        call get_s2(turb, bl, fluid, grid, du2min)
+        write(*,*) 'here3'
+        call get_rig(turb, bl, fluid, grid, du2min, ricr)
+
+        call  get_turb_length(turb, bl, fluid, grid, hbl_option)
+
+        if (turb%cl_type == 1) then ! FOM closures
+            write(*,*) 'here4'
+            call get_fo_coeffs(turb, bl, fluid,  grid)
+        end if
+        if (turb%cl_type == 2) then ! KL closures
+            !call get_kl_coeffs(turb, bl, fluid,  grid)
+        end if
+
+    end subroutine get_coeffs_general
+end module pbl_turb_common
\ No newline at end of file
diff --git a/src/pbl_turb_data.f90 b/src/pbl_turb_data.f90
index 9e55a8af34844f7a2d37ca1d16e092a0e8e48202..7f6e88db5af7ea8710571452a7f12cef4424c3bf 100644
--- a/src/pbl_turb_data.f90
+++ b/src/pbl_turb_data.f90
@@ -1,16 +1,43 @@
 ! Created by Andrey Debolskiy on 26.11.2024.
 
 module pbl_turb_data
+    !<@brief turbulent fields
     type turbBLDataType
-        real, allocatable:: uw(:), vw(:), thetaw(:), thetavw(:), qvw(:)
-        real, allocatable:: s2(:),n2(:)
-        real, allocatable :: qcw(:), cloud_frac_w(:)
-        real, allocatable:: tke(:), eps(:), l_turb(:)
-        real, allocatable:: rig(:), prt(:), rif(:)
-        real, allocatable:: km(:), kh(:)
+        real, allocatable:: uw(:) !< u-component of kinematic vertical momentum flux, [m2/s2]
+        real, allocatable:: vw(:) !< v-component of kinematic vertical momentum flux, [m2/s2]
+        real, allocatable:: thetaw(:) !<  kinematic vertical heat flux flux, [K m/s2]
+        real, allocatable::thetavw(:) !<  kinematic vertical virtual temperature  flux, [ K m/s2]
+        real, allocatable:: qvw(:)  !<  kinematic vertical moisture  flux, [ bar m/s]
+        real, allocatable:: s2(:)   !< square of shear frequancy [s-2]
+        real, allocatable:: n2(:)   !< Brunt-Vaisala frequancy [s-2]
+        real, allocatable :: qcw(:) !< turbulent vertical liquid water flux [kg/kg m /s2]
+        real, allocatable:: cloud_frac_w(:) !< turbulent vertical cloud fraction flux [ m /s2]
+        real, allocatable:: tke(:) !< turbulent kinetic energy [m2/s2]
+        real, allocatable:: eps(:) !< dissipation rate of turbulent kinetic energy [m2/s3]
+        real, allocatable:: l_turb_m(:) !<  turbulent length scale for momentum [m]
+        real, allocatable:: l_turb_h(:) !<  turbulent length scale for heat  [m]
+        real, allocatable:: rig(:)  !< gradient Richardson number
+        real, allocatable::prt(:)  !< turbulent Prandtl number
+        real, allocatable::rif(:)  !< flux Richardson number
+        real, allocatable:: km(:) !< kinematic turbulent viscosity coefficient [ m2 / s]
+        real, allocatable:: kh(:) !< kinematic turbulent diffusion coefficient [ m2 / s]
+        real:: ustr !< surface dynamic velocity scale
+        real:: tstar !< surface dynamic temperature scale
         integer :: kmax
+        integer :: cl_type !< closure type - "FOM"/"KL"
+        integer :: sf_type !< stability functions type - "Louis"/"Esau"/"EFB"
+
     end type turbBLDataType
-    public turb_data_allocate, turb_data_deallocate
+    !closure types
+    character(len = 16), parameter :: pbl_fom_tag = 'FOM'
+    character(len = 16), parameter :: pbl_kl_tag = 'KL'
+    !FOM stab functions
+    character(len = 16), parameter :: pbl_fom_louis_tag = 'Louis'
+    character(len = 16), parameter :: pbl_fom_esau_tag= 'Esau'
+    character(len = 16), parameter :: pbl_fom_efb_tag = 'EFB'
+    public turb_data_allocate
+    public turb_data_deallocate
+    public update_fields
     contains
     subroutine turb_data_allocate(turb_data, kmax)
     implicit none
@@ -29,7 +56,8 @@ module pbl_turb_data
         allocate(turb_data%n2(kmax))
         allocate(turb_data%tke(kmax))
         allocate(turb_data%eps(kmax))
-        allocate(turb_data%l_turb(kmax))
+        allocate(turb_data%l_turb_h(kmax))
+        allocate(turb_data%l_turb_m(kmax))
         allocate(turb_data%rig(kmax))
         allocate(turb_data%rif(kmax))
         allocate(turb_data%prt(kmax))
@@ -37,7 +65,7 @@ module pbl_turb_data
         allocate(turb_data%kh(kmax))
     end subroutine turb_data_allocate
 
-    subroutine scm_data_deallocate(turb_data)
+    subroutine pbl_data_deallocate(turb_data)
         implicit none
         type(turbBLDataType), intent(inout):: turb_data
         deallocate(turb_data%uw)
@@ -53,10 +81,10 @@ module pbl_turb_data
         deallocate(turb_data%kh)
         deallocate(turb_data%tke)
         deallocate(turb_data%eps)
-        deallocate(turb_data%l_turb)
+        deallocate(turb_data%l_turb_m)
+        deallocate(turb_data%l_turb_h)
         deallocate(turb_data%rig)
         deallocate(turb_data%rif)
         deallocate(turb_data%prt)
-    end subroutine scm_data_deallocate
-
+    end subroutine pbl_data_deallocate
 end module pbl_turb_data
\ No newline at end of file
diff --git a/src/scm_state_data.f90 b/src/scm_state_data.f90
index ce08dc9ad7ffae5c7c214ffd012b5c67d7b3a6bd..bf65975d3f6e191ff8832058ba80b0325240416e 100644
--- a/src/scm_state_data.f90
+++ b/src/scm_state_data.f90
@@ -1,21 +1,31 @@
 ! Created by Andrey Debolskiy on 20.11.2024.
 
 module scm_state_data
+    !< @brief abl state def.
     implicit none
     type stateBLDataType
-        real, allocatable:: u(:), v(:), temp(:), theta(:) ,qv(:)
-        real, allocatable:: rho(:)
-        real, allocatable :: lwp(:), cloud_frac(:), s_e(:)
-        real, allocatable:: km(:), kh(:)
-        real, allocatable:: vdctq(:)
-        real, allocatable:: vdcuv(:)
-        real :: HPBLA_diag
-        real :: p0
-        integer :: ktvdm=10
-        integer :: ktvd
-        integer :: kpbl
-        integer :: ktop
-        integer :: kmax
+        real, allocatable:: u(:) !< u-componet velocity, [m/s]
+        real, allocatable:: v(:) !< v-componet velocity, [m/s]
+        real, allocatable:: temp(:) !< absolute temperature, [K]
+        real, allocatable:: theta(:) !< potential temperature, [K]
+        real, allocatable:: qv(:) ! specific humidity,[kg/m3]
+        real, allocatable:: theta_v(:) !< virtual potential temperature, [K]
+        real, allocatable:: rho(:) !< air density, [K]
+        real, allocatable:: lwp(:) !< liquid water path [kg/kg]
+        real, allocatable:: cloud_frac(:) !< cloud fraction [0-1]
+        real, allocatable:: s_e(:)  !< Static energy [ m2/ s2]
+        real, allocatable:: km(:) !< kinematic viscosity koefficient  [m2 / s2]
+        real, allocatable:: kh(:) !< kinematic diffusivity koefficient [ K m / s2]
+        real, allocatable:: vdcuv(:) !< viscosity koefficient  [kg /m s2]
+        real, allocatable:: vdctq(:) !< diffusivity koefficient [ kg K / m2 s2]
+        real :: p0 !< surface pressure [bar]
+        real :: cor !< in column corioulis parameter [s-1]
+        integer :: kmax !< number of vertical cells
+        real :: hpbla_diag !< diagnostic boundary layer height [m]
+        integer :: ktvdm !< top layer index for diffusion
+        integer :: kpbl !< boundary layer height index
+        integer :: land_type !< ocean/land flag
+
     end type stateBLDataType
     public scm_data_allocate
     public scm_data_deallocate
@@ -30,14 +40,15 @@ module scm_state_data
             allocate(bl_data%v(kmax))
             allocate(bl_data%temp(kmax))
             allocate(bl_data%theta(kmax))
+            allocate(bl_data%theta_v(kmax))
             allocate(bl_data%qv(kmax))
             allocate(bl_data%lwp(kmax))
             allocate(bl_data%cloud_frac(kmax))
             allocate(bl_data%s_e(kmax))
             allocate(bl_data%km(kmax))
             allocate(bl_data%kh(kmax))
-            allocate(bl_data%vdcuv(kmax))
-            allocate(bl_data%vdctq(kmax))
+            allocate(bl_data%vdcuv(0:kmax))
+            allocate(bl_data%vdctq(0:kmax))
         end subroutine scm_data_allocate
 
         subroutine scm_data_deallocate(bl_data)
@@ -47,6 +58,7 @@ module scm_state_data
             deallocate(bl_data%v)
             deallocate(bl_data%temp)
             deallocate(bl_data%theta)
+            deallocate(bl_data%theta_v)
             deallocate(bl_data%qv)
             deallocate(bl_data%lwp)
             deallocate(bl_data%cloud_frac)
@@ -60,14 +72,15 @@ module scm_state_data
         subroutine scm_data_copy(bl, bl_old)
             implicit none
             type(stateBLDataType),intent(in):: bl_old
-            type(stateBLDataType),intent(out):: bl
+            type(stateBLDataType),intent(inout):: bl
             if (bl%kmax /= bl_old%kmax) then
-                write(*,*) 'error in copy BLData: size is not compatible'
+                write(*,*) 'error in copy BLData: size is not compatible', bl%kmax, bl_old%kmax
             else
                 bl%u  = bl_old%u
                 bl%v = bl_old%v
                 bl%temp = bl_old%temp
                 bl%theta = bl_old%theta
+                bl%theta_v = bl_old%theta_v
                 bl%qv = bl_old%qv
                 bl%lwp = bl_old%lwp
                 bl%cloud_frac= bl_old%cloud_frac
@@ -80,7 +93,7 @@ module scm_state_data
                 bl%p0 = bl_old%p0
                 bl%ktvdm = bl_old%ktvdm
                 bl%kpbl = bl_old%kpbl
-                bl%ktop = bl_old%ktop
+                bl%cor = bl_old%cor
             end if
 
         end subroutine scm_data_copy
diff --git a/src/state_utils.f90 b/src/state_utils.f90
index d4c3af35bbe0efa953110dbe6c1fba8e194b815f..8c19a71baa623e664d151ee44297185cacdac1b9 100644
--- a/src/state_utils.f90
+++ b/src/state_utils.f90
@@ -9,10 +9,12 @@ module state_utils
     public get_exner
     public get_sigma_from_z
     public get_z_from_sigma
+    public get_z_from_sigma_theta
     public get_density_theta
     public get_density
     public dqsdt_sc, qmax_sc, qmax
     public DQSDT2, DQSDT, CLDT
+    public get_theta_v
     contains
         subroutine get_geopotential_moist(F,t, qv, sig, fs, rd, eps, kmax)
             implicit none
@@ -186,6 +188,22 @@ module state_utils
         get_density_theta = (p / (fluid%R_d * t)) * (1.0 + qv) / (1 + qv * fluid%eps_v)
     end function get_density_theta
 
+    subroutine get_theta_v( theta_v, theta, qv,  fluid, kmax)
+        use phys_fluid, only: fluidParamsDataType
+        implicit none
+        integer,intent(in):: kmax
+        type(fluidParamsDataType), intent(in):: fluid
+        real,dimension(kmax), intent(in):: theta
+        real,dimension(kmax), intent(in):: qv
+        real,dimension(kmax), intent(out):: theta_v
+
+        integer k
+
+        do k = kmax,1,-1
+            theta_v(k) = theta(k) * (1.0 + fluid%R_d / fluid%R_v * qv(k))
+        end do
+    end subroutine get_theta_v
+
     subroutine get_exner( exnr, p, fluid, kmax)
         use phys_fluid, only: fluidParamsDataType
         implicit none
@@ -233,10 +251,12 @@ module state_utils
         kmax = grid%kmax
         grid%sigma_edge(kmax) = 1.0
         rho = get_density(ps, bl%temp(kmax), bl%qv(kmax), fluid)
+
         dp = - rho * fluid%g * grid%dzc(kmax)
-        grid%sigma(kmax) = (ps + dp/2) / ps
-        grid%sigma_edge(kmax) = (ps + dp) / ps
+        grid%sigma(kmax) = 1.0 + (dp/2) / ps
+        grid%sigma_edge(kmax) = 1.0
         p = ps + dp
+
         do k = kmax-1,1,-1
             rho = get_density(p, bl%temp(k), bl%qv(k), fluid)
             dp = - rho * fluid%g * grid%dzc(k)
@@ -244,7 +264,7 @@ module state_utils
             rho = get_density(p, 0.5*(bl%temp(k+1) + bl%temp(k)), &
                     0.5*(bl%qv(k+1) + bl%qv(k)), fluid)
             grid%sigma(k) = grid%sigma(k+1)  &
-                    - rho * fluid%g * grid%dze(k)
+                    - rho * fluid%g * grid%dze(k) / ps
             p = grid%sigma_edge(k) * ps
         end do
         grid%dsigmac(2:kmax) = grid%sigma(2:kmax) - grid%sigma(1:kmax-1)
@@ -252,6 +272,47 @@ module state_utils
         grid%dsigmae(1:kmax) = grid%sigma_edge(1:kmax) - grid%sigma_edge(1:kmax-1)
     end subroutine get_sigma_from_z
 
+    subroutine get_sigma_from_z_theta( grid, bl, ps, fluid)
+        use pbl_grid, only : pblgridDataType
+        use phys_fluid, only: fluidParamsDataType
+        use scm_state_data, only: stateBLDataType
+        implicit none
+        type(fluidParamsDataType), intent(in):: fluid
+        type(stateBLDataType), intent(in):: bl
+        real, intent(in):: ps
+        type(pblgridDataType), intent(inout)::grid
+
+
+        real dp, rho, p, temp, appa
+        integer k, kmax
+        kmax = grid%kmax
+        grid%sigma_edge(kmax) = 1.0
+        !convert theta -> temp
+        appa = fluid%R_d / fluid%cp
+        temp = bl%theta(kmax) * (100000.0 / ps )**appa
+        rho = get_density(ps, temp, bl%qv(kmax), fluid)
+
+        dp = - rho * fluid%g * grid%dzc(kmax)
+        grid%sigma(kmax) = 1.0 + (dp/2) / ps
+        grid%sigma_edge(kmax) = 1.0
+        p = ps + dp
+
+        do k = kmax-1,1,-1
+            temp  = bl%theta(k) * (100000.0 / p )**appa
+            rho = get_density(p, bl%temp(k), bl%qv(k), fluid)
+            dp = - rho * fluid%g * grid%dzc(k)
+            grid%sigma_edge(k) = grid%sigma_edge(k+1) + dp / ps
+            rho = get_density(p, 0.5*(bl%temp(k+1) + bl%temp(k)), &
+                    0.5*(bl%qv(k+1) + bl%qv(k)), fluid)
+            grid%sigma(k) = grid%sigma(k+1)  &
+                    - rho * fluid%g * grid%dze(k) / ps
+            p = grid%sigma_edge(k) * ps
+        end do
+        grid%dsigmac(2:kmax) = grid%sigma(2:kmax) - grid%sigma(1:kmax-1)
+        grid%sigma_edge(0) = 0.0
+        grid%dsigmae(1:kmax) = grid%sigma_edge(1:kmax) - grid%sigma_edge(1:kmax-1)
+    end subroutine get_sigma_from_z_theta
+
     subroutine get_z_from_sigma( grid, bl, zs, fluid, kmax)
         use pbl_grid, only : pblgridDataType
         use phys_fluid, only: fluidParamsDataType
diff --git a/src/tests/gabls1.f90 b/src/tests/gabls1.f90
index 0ddb027467bd967ddd9ce3419f6338765abae95e..d4b852ac4c5aa7fc2b25ea291aabda70e0ebfa87 100644
--- a/src/tests/gabls1.f90
+++ b/src/tests/gabls1.f90
@@ -10,8 +10,10 @@ program gabls1
     use scm_sfx_data, only: taux, tauy, umst, hflx, qflx, cu
     use scm_state_data
     use pbl_turb_data
+    use pbl_turb_common
     use pbl_solver
-    use state_utils, only : geo, theta2ta, ta2theta, get_sigma_from_z
+    use state_utils, only : geo, theta2ta, ta2theta, get_sigma_from_z, &
+            get_sigma_from_z_theta, get_theta_v
     use diag_pbl
     use pbl_grid
     use sfx_data, only: meteoDataType, sfxDataType
@@ -35,6 +37,9 @@ program gabls1
     real(kind=rf):: z0, zh, f_cor
     real(kind=rf):: tsurf
     real(kind=rf):: shir
+    real(kind=rf):: hbl_0
+    real(kind=rf):: tgrad
+    real(kind=rf):: cooling_rate
     type(fluidParamsDataType) fluid_params
     type(pblgridDataType) grid
 
@@ -89,8 +94,9 @@ program gabls1
 
     call init_config(fname_config,status, ierr)
     if (status == 0) then
+        write(*, *) ' FAILURE! > could not initialize config file: ', fname_config
         ierr = 1        ! signal ERROR
-        return
+        stop
     end if
 
     call c_config_is_varname("time.begin"//C_NULL_CHAR, status)
@@ -98,8 +104,9 @@ program gabls1
         !< mandatory in user dataset
         call c_config_get_float("time.begin"//C_NULL_CHAR, time_begin, status)
         if (status == 0) then
+            write(*, *) ' FAILURE! > could not set:  time.begin '
             ierr = 1        ! signal ERROR
-            return
+            stop
         end if
     end if
     call c_config_is_varname("time.end"//C_NULL_CHAR, status)
@@ -107,8 +114,9 @@ program gabls1
         !< mandatory in user dataset
         call c_config_get_float("time.end"//C_NULL_CHAR, time_end, status)
         if (status == 0) then
+            write(*, *) ' FAILURE! > could not set:  time.end '
             ierr = 1        ! signal ERROR
-            return
+            stop
         end if
     end if
     call c_config_is_varname("time.dt"//C_NULL_CHAR, status)
@@ -116,18 +124,20 @@ program gabls1
         !< mandatory in user dataset
         call c_config_get_float("time.dt"//C_NULL_CHAR, dt, status)
         if (status == 0) then
+            write(*, *) ' FAILURE! > could not set:  time.dt '
             ierr = 1        ! signal ERROR
-            return
+            stop
         end if
     end if
 
     call c_config_is_varname("time.out_dt"//C_NULL_CHAR, status)
     if ((status /= 0)) then
         !< mandatory in user dataset
-        call c_config_get_float("time.out_dt"//C_NULL_CHAR, output_dt, status)
+        call c_config_get_float("output.dt"//C_NULL_CHAR, output_dt, status)
         if (status == 0) then
+            write(*, *) ' FAILURE! > could not set:  output.dt '
             ierr = 1        ! signal ERROR
-            return
+            stop
         end if
     end if
 
@@ -136,42 +146,88 @@ program gabls1
         !< mandatory in user dataset
         call c_config_get_float("Ug"//C_NULL_CHAR, ug, status)
         if (status == 0) then
+            write(*, *) ' FAILURE! > could not set:  Ug '
             ierr = 1        ! signal ERROR
-            return
+            stop
         end if
     end if
-        call c_config_is_varname("Vg"//C_NULL_CHAR, status)
+    call c_config_is_varname("Vg"//C_NULL_CHAR, status)
     if ((status /= 0)) then
         !< mandatory in user dataset
         call c_config_get_float("Vg"//C_NULL_CHAR, vg, status)
         if (status == 0) then
+            write(*, *) ' FAILURE! > could not set:  Vg '
             ierr = 1        ! signal ERROR
-            return
+            stop
+        end if
+    end if
+    call c_config_is_varname("hbl_0"//C_NULL_CHAR, status)
+    if ((status /= 0)) then
+        !< mandatory in user dataset
+        call c_config_get_float("hbl_0"//C_NULL_CHAR, hbl_0, status)
+        if (status == 0) then
+            write(*, *) ' FAILURE! > could not set:  hbl_0 '
+            ierr = 1        ! signal ERROR
+            stop
+        end if
+    end if
+    call c_config_is_varname("Tgrad"//C_NULL_CHAR, status)
+    if ((status /= 0)) then
+        !< mandatory in user dataset
+        call c_config_get_float("Tgrad"//C_NULL_CHAR, tgrad, status)
+        if (status == 0) then
+            write(*, *) ' FAILURE! > could not set:  Tgrad '
+            ierr = 1        ! signal ERROR
+            stop
         end if
     end if
     call set_fluid_default(fluid_params)
     call get_fluid_params(fluid_params, status, ierr)
     if (status == 0) then
+        write(*, *) ' FAILURE! > could not set:  Fluid params '
         ierr = 1        ! signal ERROR
-        return
+        stop
+    end if
+    call c_config_is_varname("cooling_rate"//C_NULL_CHAR, status)
+    if ((status /= 0)) then
+        !< mandatory in user dataset
+        call c_config_get_float("cooling_rate"//C_NULL_CHAR, cooling_rate, status)
+        if (status == 0) then
+            write(*, *) ' FAILURE! > could not set:  cooling_rate '
+            ierr = 1        ! signal ERROR
+            stop
+        end if
+    end if
+    call set_fluid_default(fluid_params)
+    call get_fluid_params(fluid_params, status, ierr)
+    if (status == 0) then
+        write(*, *) ' FAILURE! > could not set:  Fluid params '
+        ierr = 1        ! signal ERROR
+        stop
     end if
     !z coordinate
     call get_grid_params(grid, status, ierr)
     if (status == 0) then
+        write(*, *) ' FAILURE! > could not set:  grid '
         ierr = 1        ! signal ERROR
-        return
+        stop
     end if
+
     kmax = grid%kmax
     time_current = time_begin
     ! hack to insure shir array is good
     shir = 3.141592653589793238 * 72.3798 / 180.0
     f_cor = 2.0 * fluid_params%omega * sin(shir)
+    bl_old%cor= f_cor
+    bl%cor = f_cor
+
     write(*,*) 'F coriolis: =  ', f_cor
     !fs(1, 1) = 0.0
     !setup surface
     z0 = 0.01
     zh = z0 / 10.0
     ! init blData states
+
     call scm_data_allocate(bl, kmax)
     call scm_data_allocate(bl_old, kmax)
     call init_state(bl, ug, vg, fluid_params%tref)
@@ -179,22 +235,27 @@ program gabls1
 
 
 
-    !TODO why zet is not initialized in beirit properly?????
-    !CALL GEO(, FS(1, 1), grid%z_cell, grid%kmax)  !change for clima
-    DO K = 1, kmax
-        grid%z_cell(k) = grid%z_cell(k) / fluid_params%g
 
+    DO K = 1, kmax
         if (grid%z_cell(k) > 100.0) then
-            bl_old%theta(k) = fluid_params%tref + 0.01* (grid%z_cell(k) - 100.0)
+            bl_old%theta(k) = fluid_params%tref + tgrad * (grid%z_cell(k) - 100.0)
         else
             bl_old%theta(k) = fluid_params%tref
         end if
         !write(*,*) k, sig(k)* fluid_params%pref* 100.0, grid%z_cell(k)
     END DO
+    !finish updating grid
+    call get_sigma_from_z_theta( grid, bl, fluid_params%p0 * 100.0, fluid_params)
     call theta2ta(bl_old%temp, bl_old%theta, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
+    bl_old%theta_v(1:kmax) = bl_old%theta(1:kmax) * (1.0+ fluid_params%eps_v * bl_old%qv(1:kmax))
+    call scm_data_copy(bl,bl_old)
+
+    !Initialize turbulent closure
+    call turb_data_allocate(turb, kmax)
+    ! getclosurefromconfig
+    call get_closure_params(turb, status, ierr)
+
 
-    !finish updating grid
-    call get_sigma_from_z( grid, bl, fluid_params%p0 * 100.0, fluid_params)
 
 
     !io setup
@@ -208,7 +269,8 @@ program gabls1
     nt=0
     write(*,*)'nt=0', nt
     do while (time_current < time_end)
-        tsurf = fluid_params%tref - 0.25 * (time_current - time_begin) / 3600.0;
+
+        tsurf = fluid_params%tref - cooling_rate * (time_current - time_begin) / 3600.0;
         meteo_cell%U = sqrt(bl_old%u(kmax)**2 + bl_old%u(kmax)**2)
         meteo_cell%dT = -tsurf + bl_old%theta(kmax)
         meteo_cell%Tsemi = 0.5 * (tsurf + bl_old%theta(kmax))
@@ -216,16 +278,20 @@ program gabls1
         meteo_cell%h = grid%z_cell(kmax)
         meteo_cell%z0_m = z0
 
+
         call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
         cu = 1000.0 * sfx_cell%Cm**2
         taux = 1.1 * sfx_cell%Cm**2 * meteo_cell%U *  bl_old%u(kmax)
         tauy = 1.1 * sfx_cell%Cm**2 * meteo_cell%U *  bl_old%v(kmax)
         hflx =  sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U * meteo_cell%dT
-        umst = sfx_cell%Cm * meteo_cell%U
+        turb%ustr = sfx_cell%Cm * meteo_cell%U
+
+        call  get_theta_v(bl%theta_v, bl%theta, bl%qv, fluid_params, grid%kmax)
+        call  get_coeffs_general(turb, bl, fluid_params, 1,  grid)
+        !call solv_diffusion(bl, bl_old, turb, fluid_params, grid)
+        !call add_coriolis(bl, ug, vg, f_cor, dt, grid)
+        call scm_data_copy(bl,bl_old)
 
-        call get_fo_diff(turb, bl, grid)
-        call solv_diffusion(bl, bl_old, turb, fluid_params, grid)
-        call add_coriolis(bl, ug, vg, f_cor, dt, grid)
 
         time_current = time_current + dt
         if (time_current >=out_dt) then
@@ -237,11 +303,12 @@ program gabls1
 
             call ta2theta( bl_old%theta, bl_old%temp, &
                     fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
+
             call diag_pblh_inmcm(grid%z_cell,bl_old%theta,shir,0.0,umst,bl_old%kmax,hpbl)
             !create output
             series_val(1) = time_current/3600.0
             series_val(2) = hflx
-            series_val(3) = umst
+            series_val(3) = turb%ustr
             series_val(4) = hpbl
             series_val(5) = 0.0
 
@@ -275,10 +342,9 @@ program gabls1
     call close_file(scan_f)
 
     call deallocate_pbl_grid(grid)
-    !deallocate(ttemp)
-    !deallocate(ttime)
-    !deallocate(utemp)
-    !deallocate(vtemp)
+    call scm_data_deallocate(bl)
+    call scm_data_deallocate(bl_old)
+    call pbl_data_deallocate(turb)
 end program gabls1
 
 
@@ -288,7 +354,10 @@ subroutine init_state(bl, ug_,vg_,tref)
 
     implicit none
     real(kind=rf):: ug_,vg_,tref
-    type(stateBLDataType), intent(out) :: bl
+    type(stateBLDataType), intent(inout) :: bl
+    integer kmax
+    kmax= bl%kmax
+    write(*,*) 'blkmax=', kmax
     bl%u(1:bl%kmax) = ug_
     bl%v(1:bl%kmax) = vg_
     bl%temp(1:bl%kmax) = tref
@@ -372,7 +441,7 @@ subroutine get_surface_from_config(model, surface, z0,zh )
         call c_config_get_string("model.id"//C_NULL_CHAR, config_field, status)
         if (status == 0) then
             ierr = 1        ! signal ERROR
-            return
+            stop
         end if
         model = get_model_id(char_array2str(config_field))
         if (model == -1) then