diff --git a/config-examples/config-gabls1-ex.txt b/config-examples/config-gabls1-ex.txt
index 3c6a9bd3b4d2f90dfa5401637b8fc0faac37d316..b677501ce38bd3c97e7ff11d8b3f78d8a3f91041 100644
--- a/config-examples/config-gabls1-ex.txt
+++ b/config-examples/config-gabls1-ex.txt
@@ -24,7 +24,7 @@ fluid {
 grid {
     type = "uniform";
 
-    nz = 16;
+    nz = 64;
     h_bot = 0.0;
     h_top = 800.0;
 }
diff --git a/src/diag_pbl.f90 b/src/diag_pbl.f90
index dbabc0ed60e8fe47e8fc53b2ea3597f13652fab4..2bef69b389ee13723934feefac5a4cb6ab5841cb 100644
--- a/src/diag_pbl.f90
+++ b/src/diag_pbl.f90
@@ -46,7 +46,7 @@ contains
 
         write(*,*) 'here_hpbl'
         dz_low = z(kl) - z_surf
-        hdyn = max(dz_low , 0.5E0 * ustar/cor)
+        hdyn = max(dz_low , 0.1E0 * ustar/cor)
         write(*,*) 'hdyn', hdyn
         kpblc = kl
         kpbld = kl
@@ -55,7 +55,7 @@ contains
             dz_hdyn =  z(k) - z_surf - HDYN
             dz_conv = theta_v(k) - theta_v(kl)
             if(kpbld.EQ.kl .AND. dz_hdyn.GE.0.E0) kpbld = k
-            if(kpblc.EQ.kl .AND. dz_conv.GT.0.E0) kpblc = k
+            if(kpblc.EQ.kl .AND. dz_conv.GT.0.E0) kpblc = k + 1
         enddo
         kpbl = MIN0(kpblc, kpbld, kl-2)
         hpbla = z(kpbl) - z_surf
diff --git a/src/pbl_dry_contrgradient.f90 b/src/pbl_dry_contrgradient.f90
index 42bc01f79431f66cb824a401708b8ebb4d2d71a5..4afbbccb5013da353bd2446db9d1c1400d6124bc 100644
--- a/src/pbl_dry_contrgradient.f90
+++ b/src/pbl_dry_contrgradient.f90
@@ -18,7 +18,7 @@ module pbl_dry_contrgradient
 
     real, parameter :: c_entr = 0.5 !< dry entrainment parameterization
     real, parameter :: b1 = 1.0; !< vertical velocity in updraft Bernulli equation constant
-    real, parameter :: b2 = 2.0; !< vertical velocity in updraft Bernulli equation constant
+    real, parameter :: b2 = 2.0/3.0; !< vertical velocity in updraft Bernulli equation constant
 
     public cbl_allocate, cbl_deallocate
 
@@ -66,7 +66,7 @@ module pbl_dry_contrgradient
 
             do k = grid%kmax, 2, -1
                 dz =(grid%z_cell(k-1) - grid%z_cell(k))
-                if (grid%z_cell(k)  > bl%hpbla_diag) then
+                if (k  <= bl%kpbl) then
                     cbl%entr(k) = 0.0
                 else
                     cbl%entr(k) = c_entr * ( 1/(grid%z_cell(k) + dz)  &
@@ -96,8 +96,8 @@ module pbl_dry_contrgradient
             cbl%thetav_up0 = (-bl%surf%hs/bl%rho(kmax)) / wstr
             cbl%theta_up0 =  (-bl%surf%hs/bl%rho(kmax)) / wstr
             cbl%qv_up0 =  (-bl%surf%es/bl%rho(kmax)) / wstr
-            cbl%ua_up0 =  (-bl%surf%cm2u * bl%u(kmax)/bl%rho(kmax)) / wstr
-            cbl%va_up0  =  (-bl%surf%cm2u * bl%v(kmax)/bl%rho(kmax)) / wstr
+            cbl%ua_up0 =  0.0
+            cbl%va_up0  =  0.0
 
             cbl%thetav_up(kmax) = cbl%thetav_up0
             cbl%theta_up(kmax) = cbl%theta_up0
@@ -106,7 +106,7 @@ module pbl_dry_contrgradient
             cbl%va_up(kmax) = cbl%va_up0
 
             do k = kmax-1 ,1,-1
-                if (grid%z_cell(k) >= bl%hpbla_diag ) then
+                if (k <= bl%kpbl ) then
                     cbl%theta_up(k) = 0.0
                     cbl%thetav_up(k) = 0.0
                     cbl%qv_up(k) = 0
@@ -205,8 +205,8 @@ module pbl_dry_contrgradient
             rhs_up(:) = 0
             do k = kmax-1, 1, -1
                 if (abs(cbl%wu(k)) > 1.0e-8) then
-                    rhs_up(k) =  0.05 * (cbl%wu(k)) *  &
-                            (cbl%theta_up(k+1) - cbl%theta_up(k) + bl%theta(k+1) - bl%theta(k))  &
+                    rhs_up(k) =  0.1 * (cbl%wu(k)) *  &
+                            (cbl%theta_up(k+1) - cbl%theta_up(k) - bl%theta(k+1) + bl%theta(k))  &
                             / (grid%z_cell(k) - grid%z_cell(k+1))
                 else
                     rhs_up(k)=0
@@ -218,13 +218,37 @@ module pbl_dry_contrgradient
             do k = kmax-1, 1, -1
                 if (abs(cbl%wu(k)) > 1.0e-7) then
                     rhs_up(k) =  0.05 * (cbl%wu(k)) *  &
-                            (cbl%qv_up(k+1) - cbl%qv_up(k)  &
-                                    + bl%qv(k+1) - bl%qv(k)) / (grid%z_cell(K) -grid%z_cell(K+1))
+                            (cbl%qv_up(+1) - cbl%qv_up(k)  &
+                                    - bl%qv(k+1) + bl%qv(k)) / (grid%z_cell(k) -grid%z_cell(k+1))
                 else
                     rhs_up(k)=0
                 end if
             end do
             bl%qv(:) = bl%qv(:) + dt * rhs_up(:)
+            !apply upwind
+            rhs_up(:) = 0
+            do k = kmax-1, 1, -1
+                if (abs(cbl%wu(k)) > 1.0e-8) then
+                    rhs_up(k) =  0.00 * (cbl%wu(k)) *  &
+                            (cbl%ua_up(k+1) - cbl%ua_up(k) - bl%u(k+1) + bl%u(k))  &
+                            / (grid%z_cell(k) - grid%z_cell(k+1))
+                else
+                    rhs_up(k)=0
+                end if
+            end do
+            bl%u(:) = bl%u(:) + dt * rhs_up(:)
+            rhs_up(:) = 0
+            do k = kmax-1, 1, -1
+                if (abs(cbl%wu(k)) > 1.0e-8) then
+                    rhs_up(k) =  0.00 * (cbl%wu(k)) *  &
+                            (cbl%va_up(k+1) - cbl%va_up(k) - bl%v(k+1) + bl%v(k))  &
+                            / (grid%z_cell(k) - grid%z_cell(k+1))
+                else
+                    rhs_up(k)=0
+                end if
+            end do
+            bl%v(:) = bl%v(:) + dt * rhs_up(:)
+
             end subroutine cbl_apply_cntrg
 
 end module pbl_dry_contrgradient
\ No newline at end of file
diff --git a/src/pbl_fo_turb.f90 b/src/pbl_fo_turb.f90
index ab33affd6b9330278576b141a645748b47f4a597..0ed5830daaca73240fcc9ebda3d3525d5120be9f 100644
--- a/src/pbl_fo_turb.f90
+++ b/src/pbl_fo_turb.f90
@@ -8,8 +8,8 @@ module pbl_fo_turb
     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:: vlinfm = 60.0e0
+    real, parameter:: vlinfh = vlinfm * sqrt(3.e0*yd/2.e0)
     real, parameter::  aka = 0.4e0
     !real, parameter:: r = 287.05
     !real, parameter:: cp = 1.005
@@ -46,7 +46,8 @@ module pbl_fo_turb
         type(turbBLDataType), intent(inout):: turb
 
         if (turb%sf_type == 1) then !> Louis stab functions
-            call get_fo_coeffs_louis(turb, bl, fluid,  grid)
+            !call get_fo_coeffs_louis(turb, bl, fluid,  grid)
+            call get_fo_coeffs_esau(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
@@ -91,6 +92,42 @@ module pbl_fo_turb
         end do
 
     end subroutine get_fo_coeffs_louis
+    subroutine get_fo_coeffs_esau(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
+        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 = grid%kmax-1, 1, -1
+            call esau_vit_stab_functions( fnriuv,fnritq,turb%rig(k), turb%l_turb_m(k), &
+                    &   grid%z_cell(k), grid%z_cell(k) - grid%z_cell(k+1))
+
+            turb%km(k) = fnriuv * turb%l_turb_m(k) * turb%l_turb_m(k) * sqrt(turb%s2(k))
+            turb%kh(k) = fnritq * turb%l_turb_h(k) * turb%l_turb_h(k) * sqrt(turb%s2(k))
+            !write(*,*) 'here_coeffs', k, turb%l_turb_m(k),  fnriuv, fnritq
+        end do
+        ! for compliance with old version
+        bl%vdctq(grid%kmax) = 0.0
+        bl%vdcuv(grid%kmax) = 0.0
+        do k = grid%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_esau
+
 
     subroutine get_turb_length(turb, bl, fluid, grid, hbl_option)
         use pbl_turb_data, only: turbBLDataType
@@ -139,8 +176,8 @@ module pbl_fo_turb
                     lh = 1.0
                 endif
             endif
-            turb%l_turb_m(k) = l_neutral * lm
-            turb%l_turb_h (k) = l_neutral * lh
+            turb%l_turb_m(k) = l_neutral * lm / ( 1.0 +  l_neutral / vlinfm )
+            turb%l_turb_h(k) = l_neutral * lh / ( 1.0 +  l_neutral / vlinfh )
         end do
 
     end subroutine get_turb_length
diff --git a/src/pbl_solver.f90 b/src/pbl_solver.f90
index f7b99a2a9aab00433a07534c2790472d8b6e45bc..ae92ec25911560046fd821f5a4605ce4f06021f2 100644
--- a/src/pbl_solver.f90
+++ b/src/pbl_solver.f90
@@ -206,7 +206,6 @@ module pbl_solver
         if (.not.(allocated(subs_qv))) then
             allocate(subs_qv(grid%kmax), source=0.0)
         end if
-
         ! central differences
         do k = kmax-1, 2, -1
             rhop = 0.5*(bl%rho(k) + bl%rho(k+1))
@@ -238,7 +237,7 @@ module pbl_solver
         end do
 
         ! bottom boundary
-        write(*,*) 'ws', grid%z_cell(kmax), grid%z_cell(kmax-1)
+        write(*,*) 'ws', grid%z_cell(kmax), grid%z_cell(kmax-1), w_s(kmax)
         k = kmax
         subs_u(k) = w_s(k-1) * (bl%u(k)-bl%u(k-1)) &
                     /(grid%z_cell(k)-grid%z_cell(k-1))
@@ -261,10 +260,10 @@ module pbl_solver
 
         !apply
         do k=kmax, 1, -1
-            bl%u(k) = bl%u(k) + subs_u(k) * dt
-            bl%v(k) = bl%v(k) + subs_v(k) * dt
-            bl%theta(k) = bl%theta(k) + subs_theta(k) * dt
-            bl%qv(k) = bl%qv(k) + subs_qv(k) * dt
+            bl%u(k) = bl%u(k) - subs_u(k) * dt
+            bl%v(k) = bl%v(k) - subs_v(k) * dt
+            bl%theta(k) = bl%theta(k) - subs_theta(k) * dt
+            bl%qv(k) = bl%qv(k) - subs_qv(k) * dt
         end do
     end subroutine apply_subsidence
 end module pbl_solver
\ No newline at end of file
diff --git a/src/tests/bomex.f90 b/src/tests/bomex.f90
new file mode 100644
index 0000000000000000000000000000000000000000..b06eabe3dd6523323e50c41b387e259beaf0280e
--- /dev/null
+++ b/src/tests/bomex.f90
@@ -0,0 +1,5 @@
+! Created by Andrey Debolskiy on 04.04.2025.
+
+program bomex
+
+end program bomex
\ No newline at end of file
diff --git a/src/tests/cbl_exp.f90 b/src/tests/cbl_exp.f90
index 9ea944c4ef78320df8844affc2c42afe53fc1558..824cc965013874576db26aee6d443b052e645058 100644
--- a/src/tests/cbl_exp.f90
+++ b/src/tests/cbl_exp.f90
@@ -299,6 +299,7 @@ program cbl_exp
         meteo_cell%dQ  = 0.0
         meteo_cell%h = grid%z_cell(kmax)
         meteo_cell%z0_m = z0
+        meteo_cell%z0_t = z0/10
 
 
         call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
@@ -308,14 +309,14 @@ program cbl_exp
         tauy =  bl%rho(kmax)* sfx_cell%Cm**2 * meteo_cell%U *  bl_old%v(kmax)
         hflx =  - bl%rho(kmax) &
                 * tw_flux0
-        bl%surf%hs = bl%rho(grid%kmax) * hflx
+        bl%surf%hs =  hflx
         bl%surf%es = sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U *meteo_cell%dQ
         turb%ustr = sfx_cell%Cm * meteo_cell%U
         umst = turb%ustr
 
         !call  get_density_theta_vector( bl, fluid_params, grid)
         call  get_theta_v(bl%theta_v, bl%theta, bl%qv, fluid_params, grid%kmax)
-        call  get_coeffs_general(turb, bl, fluid_params, 2,  grid)
+        call  get_coeffs_general(turb, bl, fluid_params, 1,  grid)
         write(*,*) 'bl%kpbl', bl%kpbl, grid%z_cell(bl%kpbl), sfx_cell%Cm, sfx_cell%Ct, tsurf
         call to_file_1d_3var('coeffs1.dat',grid%z_cell, turb%km, turb%km, bl_old%kmax)
         !do  k=1,grid%kmax
diff --git a/src/tests/dycoms2.f90 b/src/tests/dycoms2.f90
new file mode 100644
index 0000000000000000000000000000000000000000..6e0ed5bbff433b55c04cfca0b6c104edafd86e8a
--- /dev/null
+++ b/src/tests/dycoms2.f90
@@ -0,0 +1,5 @@
+! Created by Andrey Debolskiy on 04.04.2025.
+
+program dycoms2
+
+end program dycoms2
\ No newline at end of file
diff --git a/src/tests/gabls2.f90 b/src/tests/gabls2.f90
index 90739d3b1255e6431ad398225247c98be8a3744b..4ad2edf92e10ea7d93ccab8f9d13e034d82f723e 100644
--- a/src/tests/gabls2.f90
+++ b/src/tests/gabls2.f90
@@ -248,12 +248,13 @@ program gabls2
 
 
         tsurf = get_ts_gabls2(time_current)
-        meteo_cell%U = sqrt(bl_old%u(kmax)**2 + bl_old%u(kmax)**2)
+        meteo_cell%U = sqrt(bl_old%u(kmax)**2 + bl_old%v(kmax)**2)
         meteo_cell%dT = -tsurf + bl_old%theta(kmax)
         meteo_cell%Tsemi = 0.5 * (tsurf + bl_old%theta(kmax))
         meteo_cell%dQ  = 0.0
         meteo_cell%h = grid%z_cell(kmax)
         meteo_cell%z0_m = z0
+        meteo_cell%z0_t= z0 / 10.0
 
 
         call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
@@ -267,22 +268,22 @@ program gabls2
         bl%surf%es = sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U *meteo_cell%dQ
         turb%ustr = sfx_cell%Cm * meteo_cell%U
         umst = turb%ustr
-        write(*,*) 'surf', bl%surf%hs, turb%ustr, tsurf, bl%theta(kmax)
+        !write(*,*) 'surf', bl%surf%hs, turb%ustr, tsurf, bl%theta(kmax)
         call  get_density_theta_vector( bl, fluid_params, grid)
         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)
-        write(*,*) 'bl%kpbl', bl%kpbl, grid%z_cell(bl%kpbl), sfx_cell%Cm, sfx_cell%Ct, tsurf
+        !write(*,*) 'bl%kpbl', bl%kpbl, grid%z_cell(bl%kpbl), sfx_cell%Cm, sfx_cell%Ct, tsurf
         call to_file_1d_3var('coeffs1.dat',grid%z_cell, turb%km, turb%km, bl_old%kmax)
         !do  k=1,grid%kmax
         !    write(*,*) 'kdiffs', k, bl%vdctq(k), bl%vdcuv(k)
         !end do
         call solve_diffusion(bl, bl_old, turb, fluid_params, grid, dt)
-        call cbl_apply_cntrg(cbl, bl, fluid_params, grid, dt)
-        write(*,*) 'dttheta,', bl%theta- bl_old%theta
+        !call cbl_apply_cntrg(cbl, bl, fluid_params, grid, dt)
+        !write(*,*) 'dttheta,', bl%theta - bl_old%theta
         call add_coriolis(bl, bl_old, ug, vg, f_cor, dt, grid)
-        write(*,*) 'aftercoriolis,'
+        !write(*,*) 'aftercoriolis,'
         call    get_wsub(w_sub, time_current, grid)
-        !call    apply_subsidence(bl, w_sub, fluid_params, grid, dt)
+        call    apply_subsidence(bl, w_sub, fluid_params, grid, dt)
         !call to_file_1d_2var('temperature.dat', grid%z_cell, bl_old%U - bl%U, bl_old%kmax)
         call scm_data_copy(bl_old,bl)
 
@@ -479,10 +480,11 @@ subroutine get_wsub(w, curtime, grid)
     if(curtime >= 24.0* 3600.0) then
         do k= grid%kmax, 1, -1
             if (grid%z_edge(k) <= 1000.0) then
-                w(k) = -0.005 * (grid%z_edge(k) - 1000.0)
+                w(k) = - 0.005 * grid%z_edge(k) / 1000.0
             else
-                w(k) = -0.005
+                w(k) = - 0.005
             end if
+            write(*,*) 'ws', k, w(k)
         end do
     else
         w(1:grid%kmax) = 0.0