diff --git a/config-examples/config-gabls1-ex.txt b/config-examples/config-gabls1-ex.txt
index 5cf74b0bbfd195cd76d8b6f1e4ecd18379b80efd..3b20be4fe4cbfc8b49977d133a5674d2faa92c01 100644
--- a/config-examples/config-gabls1-ex.txt
+++ b/config-examples/config-gabls1-ex.txt
@@ -4,6 +4,7 @@ Vg = 0.0;
 hbl_0 = 100.0;
 Tgrad = 0.01;
 cooling_rate = -0.25 * 3600.0;
+fcoriolis = 1.39 * 0.0001;
 
 time{
     begin = 0.0;
@@ -23,7 +24,7 @@ fluid {
 grid {
     type = "uniform";
 
-    nz = 8;
+    nz = 16;
     h_bot = 0.0;
     h_top = 400.0;
 }
diff --git a/src/diag_pbl.f90 b/src/diag_pbl.f90
index da216f1f39d8bb3baeec2e00b8864d913f483401..72df7af554a5de0de7cbcf22fbe6e510e04170b1 100644
--- a/src/diag_pbl.f90
+++ b/src/diag_pbl.f90
@@ -35,7 +35,7 @@ contains
         implicit none
         real, intent(in):: cor, ustar, z_surf
         integer, intent(in):: kl
-        real, intent(in), dimension(KL):: theta_v, z
+        real, intent(in), dimension(*):: theta_v, z
         real, intent(out):: hpbla
         integer, intent(out):: kpbl
 
@@ -44,10 +44,10 @@ contains
         integer:: k, kpbld, kpblc
 
 
-
+        write(*,*) 'here_hpbl'
         dz_low = z(kl) - z_surf
         hdyn = AMIN1(dz_low , 0.5E0 * ustar/cor)
-
+        write(*,*) 'hdyn', hdyn
         kpblc = kl
         kpbld = kl
 
@@ -59,6 +59,7 @@ contains
         enddo
         kpbl = MIN0(kpblc, kpbld, KL-2)
         hpbla = z(kpbl) - z_surf
+        write(*,*) 'hpbla1, kpbl' , hpbla, kpbl
     end subroutine get_hpbl
 
     subroutine diag_pblh_inmcm(z,thetav,lat,zs,ustar,kl,hpbl)
diff --git a/src/pbl_fo_turb.f90 b/src/pbl_fo_turb.f90
index cf617e14f5d90bf4edf58bc11c2966b97631eb01..51762c9262b2f0e7dec5b82bdcdfab64323438b9 100644
--- a/src/pbl_fo_turb.f90
+++ b/src/pbl_fo_turb.f90
@@ -74,19 +74,26 @@ module pbl_fo_turb
         integer k, kmax
 
         ocean_flag = real(bl%land_type)
-        do k = kmax, 1, -1
+        write(*,*) 'here'
+        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%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)
+            !write(*,*) 'here_coeffs', k, turb%l_turb_m(k),  fnriuv, fnritq
         end do
         ! for compliance with old version
+        write(*,*) 'here'
         bl%vdcuv(grid%kmax) = 0.0
         bl%vdcuv(grid%kmax) = 0.0
-        do k = kmax-1,1, -1
+        write(*,*) 'here3.25'
+        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)
+            write(*,*)  bl%vdctq(k), bl%vdcuv(k), turb%km(k)
         end do
+        write(*,*) 'here3.5'
 
     end subroutine get_fo_coeffs_louis
 
@@ -102,25 +109,29 @@ module pbl_fo_turb
         type(pblgridDataType), intent(in)::grid
         type(turbBLDataType), intent(inout):: turb
         integer hbl_option
-        real l_neutral, lm, lh
+        real l_neutral, lm, lh, hpbla_diag, z_surf
+        integer kpbl
 
         integer k, kmax
-
+        kmax = grid%kmax
+        z_surf = grid%z_edge(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)
+            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
         else
             write(*,*) 'wrong hbl diagnostics option'
             return
         end if
 
-        do k = kmax,1, -1
+        do k = kmax, 1, -1
             l_neutral = fluid%kappa * grid%z_cell(k)
-            if ( k <= bl%kpbl) then
+            if ( k < bl%kpbl+1) then
                 lm = 0.0
                 lh = 0.0
             else
-                if (turb%rig(k) >= 0.0) then
+                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
@@ -128,8 +139,8 @@ module pbl_fo_turb
                     lh = 1.0
                 endif
             endif
-            turb%l_turb_m = l_neutral * lm
-            turb%l_turb_h = l_neutral * lh
+            turb%l_turb_m(k) = l_neutral * lm
+            turb%l_turb_h (k) = l_neutral * lh
         end do
 
     end subroutine get_turb_length
diff --git a/src/pbl_grid.f90 b/src/pbl_grid.f90
index cb1fe277857ada38feb288e995357f31dc61d1f1..43bb92417fae6b7674b51a2d7de7f8b5c9fbbca4 100644
--- a/src/pbl_grid.f90
+++ b/src/pbl_grid.f90
@@ -98,9 +98,10 @@ module pbl_grid
             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))
+                grid%z_cell(k+1) = 0.5 * (grid%z_edge(k+1) + grid%z_edge(k))
             end do
             grid%z_edge(0) = h_top
+            grid%z_cell(1) = 0.5 * (grid%z_edge(1) + grid%z_edge(0))
             grid%dzc(1:nk) = dz
             grid%dze(1:nk) = dz
         end subroutine set_pbl_grid_uniform
diff --git a/src/pbl_solver.f90 b/src/pbl_solver.f90
index ce9530e9d99e166bf1a043ef63ba98b9de508187..b5f00fde80a36d693fe95618de5471989fe06d1b 100644
--- a/src/pbl_solver.f90
+++ b/src/pbl_solver.f90
@@ -40,70 +40,138 @@ module pbl_solver
         integer, intent(in)              :: ktvd, kl
         real, intent(in),  dimension(kl) :: aa, bb, cc, ff
         real, intent(in),  dimension(kl) :: prgna, prgnz
-        real, intent(out), dimension(kl) :: y
+        real, intent(inout), dimension(kl) :: y
 
         integer :: k
         real :: prgnb(kl)
-
+        write(*,*) 'here_diff3.1', bb(ktvd)
         prgnb(ktvd) = ff(ktvd) / bb(ktvd)
+        write(*,*) 'here_diff3.1'
         do k = ktvd+1, kl
+            write(*,*) k, ktvd, prgnz(k),  ff(k)
             prgnb(k) = prgnz(k) * (ff(k) + aa(k) * prgnb(k-1))
         end do
+        write(*,*) 'here_diff3'
         y(kl) = prgnb(kl)
         do k = kl-1, ktvd, -1
             y(k) = prgna(k) * y(k+1) + prgnb(k)
         end do
     end subroutine solve_tridiag
 
-    subroutine fill_tridiag(aa, bb, cc, rho, kdiff, kbltop, grid, dt)
+    subroutine fill_tridiag(aa, bb, cc, rho, kdiff, kbltop, cm2u, grid, dt)
         use pbl_grid, only: pblgridDataType
         implicit none
         real, dimension(*), intent(in):: rho, kdiff
-        real, intent(in):: dt
+        real, intent(in):: dt, cm2u
         integer, intent(in):: kbltop
         type(pblgridDataType),intent(in):: grid
 
-        real, dimension(*), intent(out):: aa, bb, cc
+        real, dimension(*), intent(inout):: aa, bb, cc
 
         real:: dtdz
         integer:: k
 
         !nulify before top boundary
-        aa(1:kbltop) = 0.0
-        bb(1:kbltop) = 0.0
-        cc(1:kbltop) = 0.0
+        aa(1:kbltop-1) = 0.0
+        bb(1:kbltop-1) = 0.0
+        cc(1:kbltop-1) = 0.0
         !top boundary condition: flux = 0
+        write(*,*) 'here_diff1.25', kbltop
         k = kbltop
         dtdz = dt / (grid%dzc(k))
         aa(k) = 0
         cc(k) = (kdiff(k)/rho(k)) * dtdz / grid%dze(k)
-
+        write(*,*) 'here_diff1.5', kdiff(k), kbltop
+        write(*,*) 'KTVDM', k, aa(k), bb(k), cc(k), rho(k)
         do k = kbltop + 1, grid%kmax -1
             dtdz = dt / (grid%dzc(k))
             aa(k) = (kdiff(k - 1)/rho(k)) * dtdz / grid%dze(k-1)
             cc(k) = (kdiff(k)/rho(k)) * dtdz / grid%dze(k)
             bb(k) = 1.0 + aa(k) + cc(k)
+            write(*,*) 'fill', k, aa(k), bb(k), cc(k), kdiff(k)
         end do
-
+        write(*,*) 'here_diff1.75'
         !bottom boundary
         k = grid%kmax
         dtdz = dt / (grid%dzc(k))
         aa(k) = (kdiff(k-1)/rho(k)) * dtdz / grid%dze(k-1)
-        bb(k) = 1.0 + aa(k)
+        bb(k) = 1.0 + aa(k) + dtdz * cm2u / rho(k)
         cc(k) = 0.0
     end subroutine fill_tridiag
 
-    subroutine solve_diffusion(bl, bl_old, turb, fluid, grid)
+    subroutine solve_diffusion(bl, bl_old, turb, 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(stateBLDataType), intent(out):: bl
+        type(stateBLDataType), intent(inout):: bl
         type(stateBLDataType), intent(in):: bl_old
         type(turbBLDataType), intent(in):: turb
         type(fluidParamsDataType), intent(in) :: fluid
         type(pblgridDataType), intent(in) :: grid
+        real, intent(in):: dt
 
+        real, allocatable:: aa(:), bb(:), cc(:), ff(:)
+        real, allocatable:: prgna(:), prgnz(:)
+        integer k, integer, ktop, kmax
+        kmax = grid%kmax
+        write(*,*) 'here_diff'
+        if (.not.(allocated(aa))) then
+            allocate(aa(grid%kmax), source=0.0)
+        end if
+        if (.not.(allocated(bb))) then
+            allocate(bb(grid%kmax), source=0.0)
+        end if
+        if (.not.(allocated(cc))) then
+            allocate(cc(grid%kmax), source=0.0)
+        end if
+        if (.not.(allocated(ff))) then
+            allocate(ff(grid%kmax), source=0.0)
+        end if
+        if (.not.(allocated(prgna))) then
+            allocate(prgna(grid%kmax), source=0.0)
+        end if
+        if (.not.(allocated(prgnz))) then
+            allocate(prgnz(grid%kmax), source=0.0)
+        end if
+        write(*,*) 'here_diff1', bl%kpbl
+        ktop = bl%kpbl
+        do k = bl%kpbl-1,1,-1
+            if(bl%vdcuv(k) > 0.e0) then
+                ktop = k
+            end if
+        enddo
+        ! fill for temperature and specific humidity
+        call fill_tridiag(aa, bb, cc, bl%rho, bl%vdctq, ktop, 0.0,  grid, dt)
+        write(*,*) 'here_diff2'
+        call factorize_tridiag(ktop, kmax, aa, bb, cc, prgna, prgnz)
+        do k = ktop,kmax-1
+            ff(k) = bl%theta(k)
+            write(*,*) '2', k, aa(k), bb(k), cc(k), ff(k), grid%dze(k)
+        end do
+        ff(kmax) = bl%theta(kmax) + (dt/kmax) * bl%surf%hs /bl%rho(kmax)
+        write(*,*) 'here_diff3'
+        call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%theta)
+        write(*,*) 'here_diff4'
+        do k = ktop,kmax-1
+            ff(k) = bl%qv(k)
+        end do
+        ff(kmax) = bl%qv(kmax) + (dt/kmax) * bl%surf%es /bl%rho(kmax)
+        call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%qv)
+        write(*,*) 'here_diff4'
+        !velocity
+        call fill_tridiag(aa, bb, cc, bl%rho, bl%vdcuv, ktop, bl%surf%cm2u, grid, dt)
+        call factorize_tridiag(ktop, kmax, aa, bb, cc, prgna, prgnz)
+        do k = ktop,kmax-1
+            ff(k) = bl%u(k)
+        end do
+        ff(kmax) = bl%u(kmax)
+        call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%u)
+        do k = ktop,kmax-1
+            ff(k) = bl%v(k)
+        end do
+        ff(kmax) = bl%v(kmax)
+        call solve_tridiag(ktop, kmax, aa, bb, cc, ff, prgna, prgnz, bl%v)
     end subroutine solve_diffusion
 end module pbl_solver
\ No newline at end of file
diff --git a/src/pbl_turb_common.f90 b/src/pbl_turb_common.f90
index 948abbccd7ee0808301df7cb2e2c85b1cfc51e26..437d4135e4f49020f88673bd2590017943960684 100644
--- a/src/pbl_turb_common.f90
+++ b/src/pbl_turb_common.f90
@@ -43,7 +43,6 @@ contains
         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)) &
@@ -66,7 +65,6 @@ contains
     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)
@@ -95,14 +93,14 @@ contains
         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_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'
+            write(*,*) 'here3'
             call get_fo_coeffs(turb, bl, fluid,  grid)
+            write(*,*) 'here4'
         end if
         if (turb%cl_type == 2) then ! KL closures
             !call get_kl_coeffs(turb, bl, fluid,  grid)
diff --git a/src/scm_state_data.f90 b/src/scm_state_data.f90
index bf65975d3f6e191ff8832058ba80b0325240416e..878dc80442aa1693a5bc1da156244176083be684 100644
--- a/src/scm_state_data.f90
+++ b/src/scm_state_data.f90
@@ -3,6 +3,11 @@
 module scm_state_data
     !< @brief abl state def.
     implicit none
+    type surfBLDataType
+        real :: hs !< surface sensible heat flux W/m2
+        real :: es !< surface specific humidity  flux kg/m2 s
+        real :: cm2u !< Cd * |U| for implicit momentum flux in tridiag solver
+    end type surfBLDataType
     type stateBLDataType
         real, allocatable:: u(:) !< u-componet velocity, [m/s]
         real, allocatable:: v(:) !< v-componet velocity, [m/s]
@@ -25,8 +30,10 @@ module scm_state_data
         integer :: ktvdm !< top layer index for diffusion
         integer :: kpbl !< boundary layer height index
         integer :: land_type !< ocean/land flag
+        type(surfBLDataType) :: surf
 
     end type stateBLDataType
+
     public scm_data_allocate
     public scm_data_deallocate
     public scm_data_copy
@@ -44,6 +51,7 @@ module scm_state_data
             allocate(bl_data%qv(kmax))
             allocate(bl_data%lwp(kmax))
             allocate(bl_data%cloud_frac(kmax))
+            allocate(bl_data%rho(kmax))
             allocate(bl_data%s_e(kmax))
             allocate(bl_data%km(kmax))
             allocate(bl_data%kh(kmax))
@@ -62,6 +70,7 @@ module scm_state_data
             deallocate(bl_data%qv)
             deallocate(bl_data%lwp)
             deallocate(bl_data%cloud_frac)
+            deallocate(bl_data%rho)
             deallocate(bl_data%s_e)
             deallocate(bl_data%km)
             deallocate(bl_data%kh)
diff --git a/src/state_utils.f90 b/src/state_utils.f90
index 8c19a71baa623e664d151ee44297185cacdac1b9..a3404fac4070de5eeec7cfeff92b8340f5eb2e12 100644
--- a/src/state_utils.f90
+++ b/src/state_utils.f90
@@ -184,10 +184,31 @@ module state_utils
         real, intent(in):: p, th, qv
 
         real t
-        t = (p / fluid%p0)**(fluid%R_d / fluid%cp)
-        get_density_theta = (p / (fluid%R_d * t)) * (1.0 + qv) / (1 + qv * fluid%eps_v)
+        t =  th * (p / fluid%p0)**(fluid%R_d / fluid%cp)
+        get_density_theta = (p * 100.0 / (fluid%R_d * t)) * (1.0 + qv) / (1 + qv * fluid%eps_v)
     end function get_density_theta
 
+    subroutine get_density_theta_vector( bl, fluid, grid)
+        use phys_fluid, only: fluidParamsDataType
+        use scm_state_data, only: stateBLDataType
+        use pbl_grid, only: pblgridDataType
+        implicit none
+        type(fluidParamsDataType),intent(in):: fluid
+        type(pblgridDataType),intent(in):: grid
+        type(stateBLDataType),intent(inout)::bl
+
+        real t, p, qv
+        integer k
+
+        do k=grid%kmax, 1, -1
+            p = bl%p0 * grid%sigma(k)
+            t = bl%theta(k) * (p / fluid%p0)**(fluid%R_d / fluid%cp)
+            qv = bl%qv(k)
+            bl%rho(k) = (p * 100.0 / (fluid%R_d * t)) * (1.0 + qv) / (1 + qv * fluid%R_v/fluid%R_d)
+        end do
+
+    end subroutine get_density_theta_vector
+
     subroutine get_theta_v( theta_v, theta, qv,  fluid, kmax)
         use phys_fluid, only: fluidParamsDataType
         implicit none
diff --git a/src/tests/gabls1.f90 b/src/tests/gabls1.f90
index d4b852ac4c5aa7fc2b25ea291aabda70e0ebfa87..f5c2bd1b9b8e8d72d000a3c3cdf3fb7967650b1e 100644
--- a/src/tests/gabls1.f90
+++ b/src/tests/gabls1.f90
@@ -13,7 +13,7 @@ program gabls1
     use pbl_turb_common
     use pbl_solver
     use state_utils, only : geo, theta2ta, ta2theta, get_sigma_from_z, &
-            get_sigma_from_z_theta, get_theta_v
+            get_sigma_from_z_theta, get_theta_v, get_density_theta_vector
     use diag_pbl
     use pbl_grid
     use sfx_data, only: meteoDataType, sfxDataType
@@ -198,6 +198,16 @@ program gabls1
             stop
         end if
     end if
+    call c_config_is_varname("fcoriolis"//C_NULL_CHAR, status)
+    if ((status /= 0)) then
+        !< mandatory in user dataset
+        call c_config_get_float("fcoriolis"//C_NULL_CHAR, f_cor, status)
+        if (status == 0) then
+            write(*, *) ' FAILURE! > could not set:  fcoriolis '
+            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
@@ -216,13 +226,14 @@ program gabls1
     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)
+    !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
+    bl%land_type = 0.0
+    bl%p0=fluid_params%p0
+
 
-    write(*,*) 'F coriolis: =  ', f_cor
-    !fs(1, 1) = 0.0
     !setup surface
     z0 = 0.01
     zh = z0 / 10.0
@@ -247,7 +258,7 @@ program gabls1
     !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 get_theta_v(bl%theta_v, bl%theta, bl%qv, fluid_params, grid%kmax)
     call scm_data_copy(bl,bl_old)
 
     !Initialize turbulent closure
@@ -280,15 +291,19 @@ program gabls1
 
 
         call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
-        cu = 1000.0 * sfx_cell%Cm**2
+        bl%surf%cm2u = sfx_cell%Cm**2 * meteo_cell%U
         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
+        bl%surf%hs = bl%rho(grid%kmax) * hflx
+        bl%surf%es = sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U *meteo_cell%dQ
         turb%ustr = sfx_cell%Cm * meteo_cell%U
 
+        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 solv_diffusion(bl, bl_old, turb, fluid_params, grid)
+        write(*,*) 'bl%kpbl', bl%kpbl
+        call solve_diffusion(bl, bl_old, turb, fluid_params, grid, dt)
         !call add_coriolis(bl, ug, vg, f_cor, dt, grid)
         call scm_data_copy(bl,bl_old)