diff --git a/config-examples/config-cbl-ex.txt b/config-examples/config-cbl-ex.txt
index 9529c82e4cede6e93cffed45c35c81806ec02471..2c7b461d1426e6dff1a5da880007b4d6f129e4ea 100644
--- a/config-examples/config-cbl-ex.txt
+++ b/config-examples/config-cbl-ex.txt
@@ -4,13 +4,14 @@ Vg = 0.0;
 hbl_0 = 100.0;
 Tgrad = 0.03;
 fcoriolis = 1.39 * 0.0001;
-tw_flux_bot = -0.35;
+tw_flux_bot = 0.35;
 
 time{
     begin = 0.0;
     end = 1.0 * 3600.0;
-    dt = 5.0;
-    out_dt = 360.0;
+    dt = 180.0;
+    out_dt = 1.0;
+
 }
 
 fluid {
@@ -37,7 +38,7 @@ surface_model {
 }
 
 output{
-    dt = 0.2 *3600.0;
+    dt = 30.0; #0.2 *3600.0;
 }
 
 closure {
diff --git a/config-examples/config-gabls1-ex.txt b/config-examples/config-gabls1-ex.txt
index b677501ce38bd3c97e7ff11d8b3f78d8a3f91041..998901d751f34dcceca95a0f38f6af5db9504492 100644
--- a/config-examples/config-gabls1-ex.txt
+++ b/config-examples/config-gabls1-ex.txt
@@ -24,7 +24,7 @@ fluid {
 grid {
     type = "uniform";
 
-    nz = 64;
+    nz = 8;
     h_bot = 0.0;
     h_top = 800.0;
 }
diff --git a/config-examples/config-gabls2-ex.txt b/config-examples/config-gabls2-ex.txt
index a4f0a0d0a99b819f2a58e3abe5d5598135ad7747..7340e7a24b3617230ad85314ea1e7183d7eb3cfd 100644
--- a/config-examples/config-gabls2-ex.txt
+++ b/config-examples/config-gabls2-ex.txt
@@ -6,7 +6,7 @@ fcoriolis = 1.39 * 0.0001;
 time{
     begin = 0.0;
     end = 60.0 * 3600.0;
-    dt = 20.0;
+    dt = 60.0;
     out_dt = 360.0;
 }
 
@@ -21,7 +21,7 @@ fluid {
 grid {
     type = "uniform";
 
-    nz = 32;
+    nz = 24;
     h_bot = 0.0;
     h_top = 2000.0;
 }
@@ -34,7 +34,7 @@ surface_model {
 }
 
 output{
-    dt = 0.2 *3600.0;
+    dt = 60.0;
 }
 
 closure {
diff --git a/src/diag_pbl.f90 b/src/diag_pbl.f90
index 2bef69b389ee13723934feefac5a4cb6ab5841cb..33acdaebe2f831a1f632267f1d5256c704ee78f6 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.1E0 * ustar/cor)
+        hdyn = min(dz_low , 0.5E0 * ustar/cor)
         write(*,*) 'hdyn', hdyn
         kpblc = kl
         kpbld = kl
@@ -134,7 +134,7 @@ contains
         real,intent(out)::hpbl
         integer,intent(out)::nkpbl
         real,parameter:: g = 9.81
-        real,parameter:: Ric = 0.25
+        real,parameter:: Ric = 0.5
         real,parameter:: upper_bound = 16000.0 !upper boundary of PBLH
         real::Rib
         real du, dv, dtvirt
@@ -146,7 +146,7 @@ contains
         do k=kl-1,2,-1
             if(z(k).lt.upper_bound)then
 
-                du = (u(k)-u(kl)) * (u(k)-u(kl)) + (v(k)-v(kl)) * (v(k)-v(kl))
+                du = (u(k)) * (u(k)) + (v(k)) * (v(k))
                 dtvirt = theta(k) - theta(kl)
                 if(du <= du2min) du = 0.01
 
@@ -156,16 +156,18 @@ contains
                 kpbl=k
                 if(Rib.ge.Ric) then
                     !write(*,*) 'Rib', k, rib, du, dtvirt
+                        kpbl=k
                         hpbl = z(k)
                     exit
                 end if
             endif
         enddo
         do k=kl-1,2,-1
+            dtvirt = theta(k) - theta(kl)
             if (dtvirt > 0 .and. kpbld == kl ) kpbld = k
         end do
 
-        kpbl  = amin0(kpbl, kpbld, kl-2)
+        kpbl  = min(kpbl, kpbld+1, kl-2)
         hpbl = z(kpbl)
         hpbl = hpbl - zs
         nkpbl = kpbl
diff --git a/src/fo_stability_functions.f90 b/src/fo_stability_functions.f90
index 0f321774cc7931e5bcaccec789569d182858c81c..753b52d75b16d4b801bdf56c40e6d7700a01af2c 100644
--- a/src/fo_stability_functions.f90
+++ b/src/fo_stability_functions.f90
@@ -21,14 +21,14 @@ module fo_stability_functions
 
     !>
     subroutine louis_stab_functions(fnriuv,fnritq,rinum,vdlm, &
-            & vdlh,zkh, rolim)
+            & vdlh,zkh, oc_flag)
         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
+        integer, intent(in):: oc_flag
 
         if(rinum.ge.0.e0) then
             scf = sqrt(1.e0 + yd*rinum)
@@ -38,7 +38,7 @@ module fo_stability_functions
             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
+            if(oc_flag > 0) 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
@@ -87,5 +87,45 @@ module fo_stability_functions
         end if
     end subroutine esau_vit_stab_functions
 
+    !>esau byrkjedal 2007  +   viterbo 1999  unstable stratification
+    subroutine esau_stab_functions( fnriuv,fnritq,rinum, &
+            &   vdlm, vdlh, zkh, dz, oc_flag)
+        implicit none
+        real fnriuv,fnritq,rinum, vdlm, vdlh, zkh, dz
+        integer oc_flag
+        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
+            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(oc_flag >= 1) 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 esau_stab_functions
 
 end module fo_stability_functions
\ No newline at end of file
diff --git a/src/pbl_dry_contrgradient.f90 b/src/pbl_dry_contrgradient.f90
index 4afbbccb5013da353bd2446db9d1c1400d6124bc..e16a94bc74d49daab82cd740276025717893152d 100644
--- a/src/pbl_dry_contrgradient.f90
+++ b/src/pbl_dry_contrgradient.f90
@@ -16,9 +16,14 @@ module pbl_dry_contrgradient
         real :: theta_upH, thetav_upH, qv_upH
     end type pblContrGradDataType
 
-    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/3.0; !< vertical velocity in updraft Bernulli equation constant
+    real, parameter :: c_entr = 0.38 !< dry entrainment parameterization
+    real, parameter :: b1 = 1.8; !< vertical velocity in updraft Bernulli equation constant
+    real, parameter :: b2 = 3.5; !< vertical velocity in updraft Bernulli equation constant
+    real, parameter :: wstrmin = 1.0e-6 !< minimimal deardorff velocity
+    real, parameter :: lmin_entr = 1.0e6 !< maximum entrainement scale above CBL
+    real, parameter :: tau_dry_scale = 500.0 !< Time scale for entrainment
+    real, parameter :: a0 = 0.08 !< Updraft fraction area
+    real, parameter :: alph0 = 0.25 !< amount of surface dispersion to transfer to updraft
 
     public cbl_allocate, cbl_deallocate
 
@@ -64,14 +69,19 @@ module pbl_dry_contrgradient
             integer k
             real dz
 
+            cbl%entr(:) = 0.0
             do k = grid%kmax, 2, -1
-                dz =(grid%z_cell(k-1) - grid%z_cell(k))
-                if (k  <= bl%kpbl) then
+                dz =(grid%z_cell(k-1) - grid%z_cell(k)) + 1.0
+                if (k  < bl%kpbl+2) then
                     cbl%entr(k) = 0.0
+
                 else
-                    cbl%entr(k) = c_entr * ( 1/(grid%z_cell(k) + dz)  &
-                            + 1/(bl%hpbla_diag - grid%z_cell(k) + dz));
+                    cbl%entr(k) = c_entr * ( 1.0/(grid%z_cell(k) + dz)  &
+                            + 1.0/(bl%hpbla_diag - grid%z_cell(k) + dz) &
+                            + 1.0/ lmin_entr)
+
                 end if
+                write(*,*) 'entr, k = ' , cbl%entr(k), k, grid%z_cell(k)
             end do
             cbl%entr(1) = 0
         end subroutine cbl_get_entrainment
@@ -86,51 +96,68 @@ module pbl_dry_contrgradient
             type(stateBLDataType), intent(inout):: bl
             type(fluidParamsDataType), intent(in) :: fluid
             type(pblgridDataType), intent(in) :: grid
-            real phys_beta, wstr, dz
-
+            real phys_beta, wstr, dz, dzentrp, dzentrm, temp
             integer k, kmax
+
+            cbl%theta_up(:) = 0.0
+            cbl%thetav_up(:) = 0
+            cbl%qv_up = 0.0
+            cbl%ua_up = 0.0
+            cbl%va_up = 0.0
+
             kmax = grid%kmax
             phys_beta = 0.5 * fluid%g /(bl%theta_v(kmax) + bl%theta_v(kmax-1))
-            wstr = (phys_beta * (-bl%surf%hs/bl%rho(kmax)) &
-                    * bl%hpbla_diag)**(0.33333)
-            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
+            write(*,*) 'bl%hs, ', (bl%surf%hs), phys_beta
+            if (bl%surf%hs > 0.0) then
+                wstr = (phys_beta * (bl%surf%hs/bl%rho(kmax)) &
+                        * bl%hpbla_diag)**(0.33333)
+                wstr = max(wstrmin, wstr)
+
+                cbl%thetav_up0 = (bl%surf%hs/(bl%rho(kmax)*fluid%cp)) / wstr
+                cbl%theta_up0 =  (bl%surf%hs/(bl%rho(kmax)*fluid%cp)) / wstr
+                cbl%qv_up0 =  (bl%surf%es/bl%rho(kmax)) / wstr
+            else
+                cbl%thetav_up0 = 0.0
+                cbl%theta_up0 = 0.0
+                cbl%qv_up0 = 0.0
+            end if
             cbl%ua_up0 =  0.0
             cbl%va_up0  =  0.0
 
-            cbl%thetav_up(kmax) = cbl%thetav_up0
-            cbl%theta_up(kmax) = cbl%theta_up0
-            cbl%qv_up(kmax) = cbl%qv_up0
-            cbl%ua_up(kmax) = cbl%ua_up0
-            cbl%va_up(kmax) = cbl%va_up0
+            cbl%thetav_up(kmax) = cbl%thetav_up0 + bl%theta_v(kmax)
+            cbl%theta_up(kmax) = cbl%theta_up0  + bl%theta(kmax)
+            cbl%qv_up(kmax) = cbl%qv_up0  + bl%qv(kmax)
+            cbl%ua_up(kmax) = cbl%ua_up0  + bl%u(kmax)
+            cbl%va_up(kmax) = cbl%va_up0  + bl%v(kmax)
 
             do k = kmax-1 ,1,-1
                 if (k <= bl%kpbl ) then
-                    cbl%theta_up(k) = 0.0
-                    cbl%thetav_up(k) = 0.0
-                    cbl%qv_up(k) = 0
+                    cbl%theta_up(k) = bl%theta(k)
+                    cbl%thetav_up(k) = bl%theta_v(k)
+                    cbl%qv_up(k) = bl%qv(k)
+                    cbl%ua_up(k) = bl%u(k)
+                    cbl%ua_up(k) = bl%u(k)
+                    cbl%va_up(k) = bl%v(k)
                 else
                     dz = (grid%z_cell(k) - grid%z_cell(k+1))
-                    cbl%theta_up(k) = cbl%theta_up(k + 1) &
-                            - dz * cbl%entr(k+1) * cbl%theta_up(k+1)
-                    cbl%thetav_up(k) = cbl%thetav_up(k + 1)  &
-                            - dz * cbl%entr(k+1) * (cbl%thetav_up(k+1))
-                    cbl%qv_up(k) = cbl%qv_up(k + 1)  &
-                            - dz * cbl%entr(k+1) * (cbl%qv_up(k+1))
-                    cbl%ua_up(k) = cbl%ua_up(k + 1)  &
-                            - dz * cbl%entr(k+1) * (cbl%ua_up(k+1))
-                    cbl%ua_up(k) = cbl%va_up(k + 1)  &
-                            - dz * cbl%entr(k+1) * (cbl%va_up(k+1))
+                    cbl%thetav_up(k) = cbl%thetav_up(k+1) - &
+                            dz * cbl%entr(k+1) &
+                                    * (cbl%thetav_up(k+1) - bl%theta_v(k+1))
+                    cbl%theta_up(k) = cbl%theta_up(k+1) - &
+                            dz * cbl%entr(k+1) &
+                                    * (cbl%theta_up(k+1) - bl%theta(k+1))
+                    write(*,*) 'theta_up, k = ' , cbl%theta_up(k) - bl%theta(k), k
+                    cbl%qv_up(k) = cbl%qv_up(k+1) - &
+                            dz * cbl%entr(k+1) &
+                                    * (cbl%qv_up(k+1) - bl%qv(k+1))
+                    cbl%ua_up(k) = cbl%ua_up(k+1) + &
+                            dz * cbl%entr(k+1) &
+                                    * (cbl%ua_up(k+1) - bl%u(k+1))
+                    cbl%va_up(k) = cbl%va_up(k+1) - &
+                            dz * cbl%entr(k+1) &
+                                    * (cbl%va_up(k+1) - bl%v(k+1))
                 end if
             end do
-            do k = kmax ,1,-1
-                cbl%thetav_up(k) = cbl%thetav_up(k) + bl%theta_v(k)
-                cbl%theta_up(k) = cbl%theta_up(k) + bl%theta(k)
-                cbl%qv_up(k) = cbl%qv_up(k) + bl%qv(k)
-                cbl%ua_up(k) = cbl%ua_up(k) + bl%u(k)
-                cbl%va_up(k) = cbl%va_up(k) + bl%v(k)
-            end do
 
         end subroutine cbl_get_up
 
@@ -153,7 +180,7 @@ module pbl_dry_contrgradient
         umod = sqrt(bl%u(kmax)*bl%u(kmax) + bl%v(kmax)*bl%v(kmax))
         ustr = sqrt(bl%surf%cm2u /bl%rho(kmax) * umod)
         phys_beta = 0.5 * fluid%g /(bl%theta_v(kmax) + bl%theta_v(kmax-1))
-        cbl%wu(kmax) = 2.46 * ((ustr)**3  &
+        cbl%wu(kmax) = 2.46 * ((ustr)**3   &
                     - 0.4 * phys_beta * bl%surf%hs  &
                         * (grid%z_cell(k))/bl%rho(kmax))**(1.0/3.0) &
                             * sqrt(1.0 - grid%z_cell(k) / bl%hpbla_diag)
@@ -204,7 +231,7 @@ module pbl_dry_contrgradient
             !apply upwind
             rhs_up(:) = 0
             do k = kmax-1, 1, -1
-                if (abs(cbl%wu(k)) > 1.0e-8) then
+                if (abs(cbl%wu(k)) > 1.0e-7) then
                     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))
@@ -217,7 +244,7 @@ module pbl_dry_contrgradient
             rhs_up = 0
             do k = kmax-1, 1, -1
                 if (abs(cbl%wu(k)) > 1.0e-7) then
-                    rhs_up(k) =  0.05 * (cbl%wu(k)) *  &
+                    rhs_up(k) =  0.1 * (cbl%wu(k)) *  &
                             (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
@@ -225,11 +252,11 @@ module pbl_dry_contrgradient
                 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)) *  &
+                    rhs_up(k) =  0.1 * (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
@@ -240,7 +267,7 @@ 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.00 * (cbl%wu(k)) *  &
+                    rhs_up(k) =  0.1 * (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
@@ -251,4 +278,211 @@ module pbl_dry_contrgradient
 
             end subroutine cbl_apply_cntrg
 
+    subroutine new_cntrg(cbl, bl, fluid, grid, dt)
+        use scm_state_data, only : stateBLDataType
+        use pbl_turb_data, only : turbBLDataType
+        use phys_fluid, only: fluidParamsDataType
+        use pbl_grid, only : pblgridDataType
+        implicit none
+        type(pblContrGradDataType), intent(inout) :: cbl
+        type(stateBLDataType), intent(inout):: bl
+        type(fluidParamsDataType), intent(in) :: fluid
+        type(pblgridDataType), intent(in) :: grid
+        real, intent(in) :: dt
+
+        ! local variables
+        real, save, allocatable ::  rhs_up(:), rhs(:)
+        real, save, allocatable :: bouf(:)
+        real wu, phys_beta, umod, rhs_ws, a, b
+        real   dz, dzentrp, dzentrm, temp,dth
+        real sigmaw, sigmath, sigmaqc
+        real wstr, ustr, thstar, thvstar, qstar
+        real z1, hpb, hpbmf
+        real entrint, entexp, entexpu
+        real wp, entw, wn2
+
+
+        real a0w0
+        integer kmax, kpbl, k
+        kmax = grid%kmax
+        kmax = grid%kmax
+        if (.not.(allocated(rhs_up))) then
+            allocate(rhs_up(grid%kmax), source=0.0)
+        end if
+        if (.not.(allocated(rhs))) then
+            allocate(rhs(grid%kmax), source=0.0)
+        end if
+        if (.not.(allocated(bouf))) then
+            allocate(bouf(grid%kmax), source=0.0)
+        end if
+        kpbl = bl%kpbl
+        cbl%theta_up(:) = 0.0
+        cbl%thetav_up(:) = 0
+        cbl%qv_up(:) = 0.0
+        cbl%ua_up(:) = 0.0
+        cbl%va_up(:) = 0.0
+        cbl%entr(:) = 0.0
+        cbl%wu(:) = 0.0
+        rhs_up(:) = 0.0
+        rhs(:) = 0.0
+        ! Check if the flux is convective
+        if (bl%surf%hs > 0.0) then
+            ! get surface properties
+            z1  = grid%z_cell(kmax)
+            hpb = bl%hpbla_diag
+            phys_beta = 0.5 * fluid%g /(bl%theta_v(kmax) + bl%theta_v(kmax-1))
+            wstr = (phys_beta * (bl%surf%hs/bl%rho(kmax)/fluid%cp) &
+                    * bl%hpbla_diag)**(0.33333)
+            wstr = max(wstrmin, wstr)
+
+            ! Initial condition for vertical velocity
+            umod = sqrt(bl%u(kmax)*bl%u(kmax) + bl%v(kmax)*bl%v(kmax))
+            ustr = sqrt(bl%surf%cm2u /bl%rho(kmax) * umod)
+            !cbl%wu(kmax) = 1.0 * ((ustr)**3.0   & !2.46
+            !        + 0.4 * wstr**3.0)**(1.0/3.0) &
+            !        * sqrt(1.0 - grid%z_cell(kmax) / bl%hpbla_diag)
+
+            sigmaw = 1.3 * (ustr**3.0 + 0.6 * wstr**3.0 * z1/hpb) * &
+                    sqrt(1.0 - z1 / hpb)
+
+            thstar = (bl%surf%hs/(bl%rho(kmax))) / wstr
+            thvstar =  (bl%surf%hs/(bl%rho(kmax))) / wstr
+            qstar = (bl%surf%es/bl%rho(kmax)) / wstr
+            cbl%wu(kmax) = 0.5 * sigmaw
+
+
+            ! initial perturbation
+            cbl%thetav_up0 = (bl%surf%hs/(bl%rho(kmax))) / sigmaw
+            cbl%theta_up0 =  (bl%surf%hs/(bl%rho(kmax))) / sigmaw
+            cbl%qv_up0 =  (bl%surf%es/bl%rho(kmax)) / sigmaw
+            cbl%ua_up0 =  0.5 * (bl%u(kmax) + bl%u(kmax-1))
+            cbl%va_up0  =  0.5 *  (bl%v(kmax) + bl%v(kmax-1))
+            cbl%thetav_up(kmax) =  alph0 *cbl%thetav_up0 + bl%theta_v(kmax)
+            cbl%theta_up(kmax) = alph0 * cbl%theta_up0  + bl%theta(kmax)
+            cbl%qv_up(kmax) = alph0 * cbl%qv_up0  + bl%qv(kmax)
+            cbl%ua_up(kmax) = cbl%ua_up0  + bl%u(kmax)
+            cbl%va_up(kmax) = cbl%va_up0  + bl%v(kmax)
+
+
+
+            !Initial entrainment
+            !cbl%entr(kmax) = 1.0 / (cbl%wu(kmax) * tau_dry_scale)  &
+                            !+ c_entr / grid%z_cell(kmax)
+            dz = grid%z_cell(kmax) - grid%z_edge(kmax)
+            cbl%entr(kmax) = c_entr * ( 1.0/(grid%z_cell(kmax) + dz)  &
+                    + 1.0/(bl%hpbla_diag - grid%z_cell(kmax) + dz))
+
+            bl%kpbld_mf = bl%kpbl
+            ! get non-entraiened parcel height
+            !do k = kmax - 1, 2, -1
+            !    dth = cbl%thetav_up(kmax) - bl%theta_v(k)
+            !    if (bl%kpbld_mf == 0 .and. dth < 0.0 ) then
+            !        bl%kpbld_mf = k
+            !    end if
+            !end do
+            hpbmf = 0.5 * (grid%z_cell(bl%kpbld_mf) + &
+                    bl%hpbla_diag)
+            write(*,*) 'hbp, hmf', grid%z_cell(bl%kpbld_mf), bl%hpbla_diag, bl%surf%hs
+            ! get proper entrainment
+            do k = grid%kmax-1,1,-1
+                dz = (grid%z_cell(k) - grid%z_cell(k+1))
+                if (grid%z_cell(k) <= hpbmf) then
+                cbl%entr(k) = c_entr * ( 1.0/(grid%z_cell(k) + dz)  &
+                        + 1.0/(hpbmf - grid%z_cell(k) + dz))
+                else
+                    cbl%entr(k) = 0.0
+                end if
+            end do
+
+            !get updrafts
+            do k = grid%kmax-1,bl%kpbld_mf,-1
+                dz = (grid%z_cell(k) - grid%z_cell(k+1))
+                entexp = exp( - cbl%entr(k) * dz)
+                entexpu = exp( - cbl%entr(k) * dz / 3.0 )
+                cbl%thetav_up(k) = bl%theta_v(k) * (1.0 - entexp) + &
+                        cbl%thetav_up(k+1) * entexp
+                cbl%theta_up(k) = bl%theta(k) * (1.0 - entexp) + &
+                        cbl%theta_up(k+1) * entexp
+                cbl%qv_up(k) = bl%qv(k) * (1.0 - entexp) + &
+                        cbl%qv_up(k+1) * entexp
+                cbl%ua_up(k) = bl%u(k) * (1.0 - entexpu) + &
+                        cbl%ua_up(k+1) * entexpu
+                cbl%va_up(k) = bl%u(k) * (1.0 - entexpu) + &
+                        cbl%va_up(k+1) * entexpu
+                bouf(k) = fluid%g *  &
+                     (0.5 * (cbl%thetav_up(k) + cbl%thetav_up(k+1)) / &
+                        bl%theta_v(k) - 1.0)
+            enddo
+
+            !get vertical speed
+            do k = grid%kmax-1,bl%kpbld_mf,-1
+                dz = (grid%z_cell(k) - grid%z_cell(k+1))
+                wp = b1 * cbl%entr(k)
+                if (wp > 0.0) then
+                    entw = exp( -2.0 * wp * dz )
+                    wn2 = cbl%wu(k+1) * cbl%wu(k+1) * entw + &
+                            b2* bouf(k)*(1.0-entw)/entw
+                else
+                    wn2 = cbl%wu(k+1) * cbl%wu(k+1) + &
+                            2.0 * b2 * bouf(k) * dz
+                end if
+
+                if ( wn2 > 0.0) then
+                    cbl%wu(k) = sqrt(wn2)
+                else
+                    cbl%wu(k) = 0.0
+                end if
+            end do
+
+
+            ! apply upwind
+            rhs_up(:) = 0.0
+            do k = kmax-1, 1, -1
+                if (cbl%wu(k) > 0.0) then
+                    rhs_up(k) =  a0 * (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.0
+                end if
+            end do
+            bl%theta = bl%theta + dt * rhs_up
+
+            rhs_up(:) = 0.0
+            do k = kmax-1, 1, -1
+                if (cbl%wu(k) > 0.0) then
+                    rhs_up(k) =  a0 * (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))
+                else
+                    rhs_up(k)=0.0
+                end if
+            end do
+            bl%qv(:) = bl%qv(:) + dt * rhs_up(:)
+
+            rhs_up(:) = 0.0
+            do k = kmax-1, 1, -1
+                if (cbl%wu(k) > 0.0) then
+                    rhs_up(k) =  a0 * (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.0
+                end if
+            end do
+            bl%u(:) = bl%u(:) + dt * rhs_up(:)
+            rhs_up(:) = 0.0
+            do k = kmax-1, 1, -1
+                if (cbl%wu(k) > 0.0) then
+                    rhs_up(k) =  a0 * (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.0
+                end if
+            end do
+            bl%v(:) = bl%v(:) + dt * rhs_up(:)
+        end if
+    end subroutine new_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 0ed5830daaca73240fcc9ebda3d3525d5120be9f..c6d14be3a4d82f582ffe3171d85c851173c254c5 100644
--- a/src/pbl_fo_turb.f90
+++ b/src/pbl_fo_turb.f90
@@ -22,6 +22,12 @@ module pbl_fo_turb
     public get_fo_coeffs
     public get_fo_coeffs_louis
     public get_turb_length
+    type hpbl_options
+        integer:: hpbl_old = 1
+        integer:: hpbl_fixed = 2
+        integer:: hpbl_rib = 3
+        integer:: hpbl_heffner = 4
+    end type hpbl_options
 !    public :: default_stab_functions
 !    public :: esau_stab_functions
 !    private:: get_ph_efb
@@ -70,14 +76,14 @@ module pbl_fo_turb
         type(turbBLDataType), intent(inout):: turb
 
         real fnriuv, fnritq,rho
-        real ocean_flag
+        integer ocean_flag
 
         integer k, kmax
 
-        ocean_flag = real(bl%land_type)
+        ocean_flag = bl%land_type
         do k = grid%kmax-1, 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%l_turb_h(k), grid%z_cell(k),ocean_flag)
 
             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))
@@ -105,14 +111,15 @@ module pbl_fo_turb
         type(turbBLDataType), intent(inout):: turb
 
         real fnriuv, fnritq,rho
-        real ocean_flag
+        integer ocean_flag
 
         integer k, kmax
 
-        ocean_flag = real(bl%land_type)
+        ocean_flag = 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))
+            call esau_stab_functions( fnriuv,fnritq,turb%rig(k), &
+                    turb%l_turb_m(k),turb%l_turb_h(k),   grid%z_cell(k), &
+                    grid%z_cell(k) - grid%z_cell(k+1), ocean_flag)
 
             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))
@@ -147,7 +154,12 @@ module pbl_fo_turb
         integer k, kmax
         kmax = grid%kmax
         z_surf = grid%z_edge(kmax)
-        if (hbl_option == 1) then
+        if (hbl_option == 0) then
+            call get_hpbl(hpbla_diag, kpbl, bl%theta_v, grid%z_cell, &
+                    z_surf, grid%kmax , bl%cor, turb%ustr)
+            bl%hpbla_diag = hpbla_diag
+            bl%kpbl = kpbl
+        elseif (hbl_option == 1) then
             call get_hpbl(hpbla_diag, kpbl, bl%theta_v, grid%z_cell, &
                     z_surf, grid%kmax , bl%cor, turb%ustr)
             bl%hpbla_diag = hpbla_diag
@@ -157,6 +169,11 @@ module pbl_fo_turb
                     bl%u, bl%v, grid%kmax , grid%z_edge(kmax), du2min)
             bl%hpbla_diag = hpbla_diag
             bl%kpbl = kpbl
+        elseif (hbl_option == 3) then
+            call diag_pblh_rib(hpbla_diag,kpbl, bl%theta_v,grid%z_cell, &
+                    bl%u, bl%v, grid%kmax , grid%z_edge(kmax), du2min)
+            bl%hpbla_diag = hpbla_diag
+            bl%kpbl = kpbl
         else
             write(*,*) 'wrong hbl diagnostics option'
             return
diff --git a/src/pbl_solver.f90 b/src/pbl_solver.f90
index ae92ec25911560046fd821f5a4605ce4f06021f2..d14e36f502406d85a495e2f31ec3311a8fa15c41 100644
--- a/src/pbl_solver.f90
+++ b/src/pbl_solver.f90
@@ -146,7 +146,7 @@ module pbl_solver
 
         ff(kmax) = bl%theta(kmax)  &
                     + (dt/grid%dze(kmax)) * bl%surf%hs / (bl%rho(kmax))
-        write(*,*) 'ff'!,  ff!, aa,bb,cc
+        !write(*,*) 'ff'!,  ff!, aa,bb,cc
 
         call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%theta)
 
@@ -237,7 +237,6 @@ module pbl_solver
         end do
 
         ! bottom boundary
-        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))
diff --git a/src/pbl_turb_common.f90 b/src/pbl_turb_common.f90
index 7cdb9f153e6d56c6438c0dc7aa4f0e0ee8ded9cd..331dda942cfdaa8ca2f617cc6b2e91fc7f50f49a 100644
--- a/src/pbl_turb_common.f90
+++ b/src/pbl_turb_common.f90
@@ -96,7 +96,7 @@ contains
         call get_rig(turb, bl, fluid, grid, du2min, ricr)
 
         call  get_turb_length(turb, bl, fluid, grid, hbl_option)
-        write(*,*)'here5'
+        !write(*,*)'here5'
         if (turb%cl_type == 1) then ! FOM closures
 
             call get_fo_coeffs(turb, bl, fluid,  grid)
diff --git a/src/scm_state_data.f90 b/src/scm_state_data.f90
index b372893aae3c3ca64454cc039106d2df38bf7f08..807be29cdb59ef8c68eb15d63bb768af31b88a7a 100644
--- a/src/scm_state_data.f90
+++ b/src/scm_state_data.f90
@@ -29,6 +29,7 @@ module scm_state_data
         real :: hpbla_diag !< diagnostic boundary layer height [m]
         integer :: ktvdm !< top layer index for diffusion
         integer :: kpbl !< boundary layer height index
+        integer :: kpbld_mf !< boundary layer height index for dry mf
         integer :: land_type !< ocean/land flag
         type(surfBLDataType) :: surf
 
diff --git a/src/tests/cbl_exp.f90 b/src/tests/cbl_exp.f90
index 824cc965013874576db26aee6d443b052e645058..9d83fed06bdd0c7bb13090eca43b37cefc245af5 100644
--- a/src/tests/cbl_exp.f90
+++ b/src/tests/cbl_exp.f90
@@ -59,7 +59,7 @@ program cbl_exp
 
     !io variables
     real, dimension (5):: series_val
-    type (io_struct) :: series_f, scan_f
+    type (io_struct) :: series_f, scan_f, scan_cg_f
 #ifdef USE_NETCDF
     type(variable_metadata) :: meta_z, meta_t
     type(variable_metadata) :: meta_temperature, meta_uwind, meta_vwind
@@ -88,6 +88,10 @@ program cbl_exp
 
     scan_f%fname='time_scan_cbl.dat'
     call set_file(scan_f, scan_f%fname)
+
+    scan_cg_f%fname='time_scan_cgp_cbl.dat'
+    call set_file(scan_cg_f, scan_cg_f%fname)
+
 #ifdef USE_NETCDF
     call set_meta_nc(meta_z, meta_t, meta_temperature, meta_uwind, meta_vwind)
 #endif
@@ -289,7 +293,7 @@ program cbl_exp
     write(*,*)'nt=0', nt
     call to_file_1d_2var('temperature1.dat', grid%z_cell, bl%theta, bl_old%kmax)
     do while (time_current < time_end)
-        call put_tscan( time_current/3600.0, grid%z_cell,bl_old, bl_old%kmax, scan_f)
+        !call put_tscan( time_current/3600.0, grid%z_cell,bl_old, bl_old%kmax, scan_f)
 
 
         tsurf = bl_old%theta(kmax)
@@ -299,7 +303,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
+        !meteo_cell%z0_t = z0/10
 
 
         call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
@@ -307,7 +311,7 @@ program cbl_exp
         bl%surf%cm2u = bl%rho(kmax)* sfx_cell%Cm**2 * meteo_cell%U
         taux =  bl%rho(kmax)* sfx_cell%Cm**2 * meteo_cell%U *  bl_old%u(kmax)
         tauy =  bl%rho(kmax)* sfx_cell%Cm**2 * meteo_cell%U *  bl_old%v(kmax)
-        hflx =  - bl%rho(kmax) &
+        hflx =  bl%rho(kmax) &
                 * tw_flux0
         bl%surf%hs =  hflx
         bl%surf%es = sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U *meteo_cell%dQ
@@ -316,14 +320,17 @@ program cbl_exp
 
         !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)
+        call  get_coeffs_general(turb, bl, fluid_params, 2,  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
         !    write(*,*) 'kdiffs', k, bl%vdctq(k), bl%vdcuv(k)
         !end do
+        call new_cntrg(cbl, bl, fluid_params, grid, dt)
         call solve_diffusion(bl, bl_old, turb, fluid_params, grid, dt)
         !call cbl_apply_cntrg(cbl, bl, fluid_params, grid, dt)
+        call  get_theta_v(bl%theta_v, bl%theta, bl%qv, fluid_params, grid%kmax)
+
         !write(*,*) 'dttheta,', bl%theta- bl_old%theta
         call add_coriolis(bl, bl_old, ug, vg, f_cor, dt, grid)
         !call to_file_1d_2var('temperature.dat', grid%z_cell, bl_old%U - bl%U, bl_old%kmax)
@@ -352,6 +359,7 @@ program cbl_exp
             call write_series(series_val, 5, series_f)
 
             call put_tscan( time_current/3600.0, grid%z_cell,bl_old, bl_old%kmax, scan_f)
+            call put_tscan_cg( time_current/3600.0, grid%z_cell,cbl, bl_old%kmax, scan_cg_f)
 #ifdef USE_NETCDF
             do k = 1, kl
                 ttemp(k,nt) = thold(k)
@@ -459,6 +467,36 @@ subroutine put_tscan( time, z, bl, nl, f)
 
 end subroutine put_tscan
 
+subroutine put_tscan_cg( time, z, cbl, nl, f)
+    use parkinds, only : rf=>kind_rf, im => kind_im
+    use scm_io_default
+    use pbl_dry_contrgradient, only: pblContrGradDataType
+    implicit none
+    type(pblContrGradDataType), intent(in):: cbl
+    real(kind=rf), intent(in), dimension(nl):: z
+    real(kind=rf),intent(in) :: time
+    type (io_struct),intent(in) :: f
+    integer(kind=im), intent(in) :: nl
+
+
+    !local
+    real(kind=rf), dimension(5,nl)::stamp
+    integer k
+
+    ! copy to stamp
+    do k=1,nl
+        stamp(1,k) = time
+        stamp(2,k) = z(k)
+        stamp(3,k) = cbl%wu(k)
+        stamp(4,k) = cbl%theta_up(k)
+        stamp(5,k) = cbl%entr(k)
+    end do
+
+    ! call to write timestamp
+    call write_timescan(stamp,nl, 5, f)
+
+end subroutine put_tscan_cg
+
 subroutine get_surface_from_config(model, surface, z0,zh )
     use iso_c_binding, only: C_NULL_CHAR
     use config_parser
@@ -608,4 +646,5 @@ subroutine set_meta_nc(z_meta, t_meta, theta_meta, uwind_meta, vwind_meta)
                             1.0, &
                             0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
 end subroutine set_meta_nc
-#endif
\ No newline at end of file
+#endif
+
diff --git a/src/tests/gabls1.f90 b/src/tests/gabls1.f90
index 5b22557fb3cd184788065108c3a3537716598b4f..fcb818b8f5c6b73a06ccd1798ec238d351579621 100644
--- a/src/tests/gabls1.f90
+++ b/src/tests/gabls1.f90
@@ -312,7 +312,7 @@ program gabls1
 
         !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)
+        call  get_coeffs_general(turb, bl, fluid_params, 0,  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/gabls2.f90 b/src/tests/gabls2.f90
index 4ad2edf92e10ea7d93ccab8f9d13e034d82f723e..044a88d2eb5b20374ea5d21c5dad0b90f5efede6 100644
--- a/src/tests/gabls2.f90
+++ b/src/tests/gabls2.f90
@@ -18,7 +18,7 @@ program gabls2
     use diag_pbl
     use pbl_grid
     use sfx_data, only: meteoDataType, sfxDataType
-    use sfx_most, only: get_surface_fluxes_most => get_surface_fluxes, &
+    use sfx_esm, only: get_surface_fluxes_most => get_surface_fluxes, &
             numericsType_most => numericsType
     use config_parser
     use config_utils
@@ -254,7 +254,7 @@ program gabls2
         meteo_cell%dQ  = 0.0
         meteo_cell%h = grid%z_cell(kmax)
         meteo_cell%z0_m = z0
-        meteo_cell%z0_t= z0 / 10.0
+        !meteo_cell%z0_t= z0 / 10.0
 
 
         call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
@@ -264,22 +264,24 @@ program gabls2
         tauy =  bl%rho(kmax)* sfx_cell%Cm**2 * meteo_cell%U *  bl_old%v(kmax)
         hflx =  - bl%rho(kmax) &
                 * sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U * meteo_cell%dT
-        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
         !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)
+        call  get_coeffs_general(turb, bl, fluid_params, 2,  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
         !    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)
+        call solve_diffusion(bl, bl_old, turb, fluid_params, grid, dt)
+        call new_cntrg(cbl, bl, fluid_params, grid, dt)
         !write(*,*) 'dttheta,', bl%theta - bl_old%theta
+        call scm_data_copy(bl_old,bl)
         call add_coriolis(bl, bl_old, ug, vg, f_cor, dt, grid)
         !write(*,*) 'aftercoriolis,'
         call    get_wsub(w_sub, time_current, grid)